In R, we can use Matrix::nearPD() to calculate nearest positive definite matrix.
I have written a Rcpp-version, nearPD_c, myself as follows (c++ file),
// [[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);
}
However, although nearPD and nearPD_c give similar results, they are not identical. For example (in 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
There are some differences in 7th or 8th decimal, which make nearPD(mt0) positive define while nearPD_c(mt0) not.
> 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
I sense that there might be some rounding issue in Rcpp. But I couldn't identify it. Any insights of what goes wrong?
There is at least one logic error in your post-processing. In R we have:
e <- eigen(X, symmetric = TRUE)
d <- e$values
Eps <- posd.tol * abs(d[1])
if (d[n] < Eps) {
d[d < Eps] <- Eps
[...]
While you have:
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;
}
}
According to the Armadillo docs, eigen values are in ascending order, while they are in decreasing order in R. So R builds Eps based on the largest eigen value, while you use the second(!) smallest. Then R compares with the smallest eigen value, while you compare with the largest. Something like this should give the same results as R (untested):
eig_sym(d, Q, X);
double Eps = posd_tol * std::abs(d[n-1]);
if (d(0) < Eps) {
uvec d_comp = d < Eps;
for(int i=0;i<sum(d_comp);i++){
if(d_comp(i)){
d(i)=Eps;
}
}
BTW, you only need // [[Rcpp::export]] for functions that you want to call from R.