Hello,

I am new to here and to Eigen. I am trying to decompose large possibly semidefinitie covariance matricies with LDLT from eigen in order to use to draw multivariate normal samples for monte-carlo. For a lot of my smaller matrices less than 1,000 rows I am able to take the components and multiply them together to recover the original matrix. But for some large examples I seem to run into numerical issues as I get a lot of really small numbers spuriously included on the reconstructed matrix (which should have a lot of identical 0's) as well as the numbers that should be non-zero in some cases being very wrong.
Is there anything I can do to increase the precision of the calculation?? This code should return (P^T * L * sqrt(D))

mat is a const Eigen::MatrixXd &

Eigen::MatrixXd res = mat.ldlt().vectorD().asDiagonal();
_sqrt_of_diag_mat(res);
Eigen::MatrixXd L = mat.ldlt().matrixL();
res = L * res;
res = mat.ldlt().transpositionsP().transpose() * res;
return res;

Make sure you trie with the latest version of Eigen, and if it's still unstable could you share matrix allowing to reproduce the issue? thanks.

Hi ggael,
Thank you very much for your reply. I am using the latest version 3.2.4 downloaded yesterday. Here is the original matrix which leads to issues. (its actually only 627,627)

http://wikisend.com/download/266802/output.xlsx

thanks

For some reason excel keeps messing up the format of my matrix. Here is the raw matrix file.

http://wikisend.com/download/268404/output

Sorry for the multiple posts. I am still having problems with the stability of my LDLT decomposition on my positive semi-definite matrix. I believe a good way to view the matrix leading to the issue is below.

http://w1.wikisend.com/node-fs/download ... a1b/output

Thanks again

Your matrices are not symmetric.

btw, instead of performing the LDLT factorization 3 times, you should create a LDLT object: LDLT<MatrixXd> ldlt(mat);