r函数和rcpp函数计算最近正定矩阵的结果不同

在R中,我们可以使用 Matrix::nearPD() 计算最近的正定矩阵。

我写了一个RCPP版本, nearPD_c 我自己如下(C++文件),

  // [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

using namespace arma;
using namespace Rcpp;

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
vec rep_each(const vec& x, const int each) {
  std::size_t n=x.n_elem;
  std::size_t n_out=n*each;
  vec res(n_out);
  auto begin = res.begin();
  for (std::size_t i = 0, ind = 0; i < n; ind += each, ++i) {
    auto start = begin + ind;
    auto end = start + each;
    std::fill(start, end, x[i]);
  }
  return res;
}

mat mat_vec_same_len(mat mt1, vec v1){
  //do not check the input...
  int t=0;
  for(int i=0;i<mt1.n_cols;i++){
    for(int j=0;j<mt1.n_rows;j++){
      mt1(j,i)=mt1(j,i)*v1(t);
      t++;
    }
  }
  return(mt1);
}

// [[Rcpp::export]]
vec pmax_c(double a, vec b){
  vec c(b.n_elem);
  for(int i=0;i<b.n_elem;i++){
    c(i)=std::max(a,b(i));
  }
  return c;
}

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
mat nearPD_c(mat x, 
             bool corr = false, bool keepDiag = false
               ,bool do2eigen = true  // if TRUE do a sfsmisc::posdefify() eigen step
               ,bool doSym = false // symmetrize after tcrossprod()
               , bool doDykstra = true // do use Dykstra's correction
               ,bool only_values = false // if TRUE simply return lambda[j].
               , double eig_tol   = 1e-6 // defines relative positiveness of eigenvalues compared to largest
               , double conv_tol  = 1e-7 // convergence tolerance for algorithm
               ,double posd_tol  = 1e-8 // tolerance for enforcing positive definiteness
               , int maxit    = 100 // maximum number of iterations allowed
               , bool trace = false // set to TRUE (or 1 ..) to trace iterations
){

  int n = x.n_cols;
  vec diagX0;
  if(keepDiag) {
    diagX0 = x.diag();
  }
  mat D_S;
  if(doDykstra) {
    //D_S should be like x, but filled with '0' -- following also works for 'Matrix':
    D_S = x;
    D_S.zeros(); //set all element
  }

  mat X = x;
  int iter = 0 ;
  bool converged = false; 
  double conv = R_PosInf;


  mat Y;
  mat R;
  mat B;
  while (iter < maxit && !converged) {
    Y = X;
    if(doDykstra){
      R = Y - D_S;
    }

vec d;
mat Q;
if(doDykstra){
  B=R;
}else{
  B=Y;
}

eig_sym(d, Q, B);

// create mask from relative positive eigenvalues
uvec p= (d>eig_tol*d[1]);
if(sum(p)==0){
  //stop("Matrix seems negative semi-definite")
  break;
}

// use p mask to only compute 'positive' part
uvec p_indexes(sum(p));

int p_i_i=0;
for(int i=0;i<p.n_elem;i++){
  if(p(i)){
    p_indexes(p_i_i)=i;
    p_i_i++;
  }
}


Q=Q.cols(p_indexes);

X=mat_vec_same_len(Q,rep_each(d.elem(p_indexes),Q.n_rows))*Q.t();

// update Dykstra's correction D_S = \Delta S_k           
if(doDykstra){
  D_S = X - R;
}

// project onto symmetric and possibly 'given diag' matrices:
if(doSym){
  X = (X + X.t())/2;
}

if(corr){

  X.diag().ones(); //set diagnols as ones
} 
else if(keepDiag){
  X.diag() = diagX0;
} 

conv = norm(Y-X,"inf")/norm(Y,"inf");

iter = iter + 1;
if (trace){
  // cat(sprintf("iter %3d : #{p}=%d, ||Y-X|| / ||Y||= %11g\n",
  // iter, sum(p), conv))
  Rcpp::Rcout << "iter " << iter <<" : #{p}= "<< sum(p) << std::endl;
}

converged = (conv <= conv_tol);       

// force symmetry is *NEVER* needed, we have symmetric X here!
//X <- (X + t(X))/2
if(do2eigen || only_values) {
  // begin from posdefify(sfsmisc)

  eig_sym(d, Q, X);

  double Eps = posd_tol * std::abs(d[1]);
  // if (d[n] < Eps) { //should be n-1?
  if (d(n-1) < Eps) {
    uvec d_comp = d < Eps;
    for(int i=0;i<sum(d_comp);i++){
      if(d_comp(i)){
        d(i)=Eps;
      }
    }

    // d[d < Eps] = Eps; //how to assign values likes this?
    if(!only_values) {

      vec o_diag = X.diag();
      X = Q * (d *Q.t());
      vec D = sqrt(pmax_c(Eps, o_diag)/X.diag());
      x=D * X * rep_each(D,  n);
    }
  }
  if(only_values) return(d);

  // unneeded(?!): X <- (X + t(X))/2
  if(corr) {
    X.diag().ones(); //set diag as ones
  }
  else if(keepDiag){
    X.diag()= diagX0;
  } 

} //end from posdefify(sfsmisc)

}

  if(!converged){ //not converged

    Rcpp::Rcout << "did not converge! " <<std::endl;

  }

  return X;
  // return List::create(_["mat"] = X,_["eigenvalues"]=d,
  //                     
  //                     _["corr"] = corr, _["normF"] = norm(x-X, "fro"), _["iterations"] = iter,
  //                       _["rel.tol"] = conv, _["converged"] = converged);


}

然而,尽管 nearPD 和 近距离接触 给出相似的结果,它们是不相同的。例如(在r中):

> mt0=matrix(c(0.5416,  -0.0668 , -0.1538,  -0.2435,
+              -0.0668 ,  0.9836 , -0.0135 , -0.0195,
+              -0.1538 , -0.0135  , 0.0226 ,  0.0334,
+              -0.2435,  -0.0195 ,  0.0334  , 0.0487),4,byrow = T)
> nearPD(mt0)$mat
4 x 4 Matrix of class "dpoMatrix"
            [,1]        [,2]        [,3]        [,4]
[1,]  0.55417390 -0.06540967 -0.14059121 -0.22075966
[2,] -0.06540967  0.98375373 -0.01203943 -0.01698557
[3,] -0.14059121 -0.01203943  0.03650733  0.05726836
[4,] -0.22075966 -0.01698557  0.05726836  0.08983952
> nearPD_c(mt0)
            [,1]        [,2]        [,3]        [,4]
[1,]  0.55417390 -0.06540967 -0.14059123 -0.22075967
[2,] -0.06540967  0.98375373 -0.01203944 -0.01698557
[3,] -0.14059123 -0.01203944  0.03650733  0.05726837
[4,] -0.22075967 -0.01698557  0.05726837  0.08983952

小数点后7或8有一些不同,这使得 nearPD(mt0) 肯定定义while nearPD_c(mt0) 不是。

> chol(nearPD(mt0)$mat)
4 x 4 Matrix of class "Cholesky"
     [,1]          [,2]          [,3]          [,4]         
[1,]  7.444286e-01 -8.786561e-02 -1.888579e-01 -2.965491e-01
[2,]             .  9.879440e-01 -2.898297e-02 -4.356729e-02
[3,]             .             .  1.029821e-04  1.014128e-05
[4,]             .             .             .  1.071201e-04
> chol(nearPD_c(mt0))
Error in chol.default(nearPD_c(mt0)) : 
  the leading minor of order 3 is not positive definite

我感觉在RCPP中可能存在一些舍入问题。但我无法识别。有什么问题吗?