Dear all,

I am starting to use Eigen library in the context of statistical modelling. I need to solve a linear system using only the lower triangular matrix from the Cholesky decomposition of a SPD matrix. My goal is to compute C^{-1/2} V^{-1/2} where V^{-1/2} and C are known matrices.
For general matrix class MatrixXd I can do using the code below. But, now I would like to do the same using the SparseMatrix class. I tried to adapt my code, but it did not work. Any help will be welcome.
## Using MatrixXd

Code: Select all

const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));
const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));
MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV);


## Using SparseMatrix

Code: Select all

const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));
const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver;
SparseMatrix<double>  LT = solver.compute(Omega).matrixL().solve(invSqrtV); // It does not work.


Thank you in advance.

For sparse matrice, the input matrix has to be reordered to reduce fill-in and make the Cholesky decomposition tractable. Two options:

1 - your matrix Omega is already well structured (e.g., non-zeros are nearby the diagonal) in which case you can use NaturalOrdering<int> as the third template parameter of SimplicialLLT to by-pass this reordering.

2 - you need to take into account the permutation returned by solver.permutationP() or solver.permutationPinv(); Recall that P is such that: P A P^-1 = L L^T.