Hi,

I'm using the Cholesky module (Eigen 2.0) for computing matrix square roots, and it sometimes happens that for apparently well-defined matrices, LLT will produce inaccurate results, with a few values diverging.

It's a matter of fact that these matrices are never entirely symmetric, which is due to numerical issues, but for most of them LLT computes useful results nonetheless (though the documentation says it has to be symmetric). At least I can't figure out a difference between the stable and non-stable ones.

Example: the matrix

Code: Select all

999.2 -0.8884 -3.105e-24 -1.067e-15 -5.577e-16 0.1595 0 0 0 0 0
-0.8884 998.4 -7.252e-24 -2.491e-15 -8.866e-16 -0.1081 0 0 0 0 0
1.263e-25 -2.746e-24 1e-07 -2.111e-11 3.21e-12 -9.319e-26 0 0 0 0 0
4.338e-17 -9.431e-16 -2.111e-11 0.0007476 0.001103 -3.201e-17 0 0 0 0 0
-7.264e-16 -1.122e-15 3.21e-12 0.001103 0.001626 -7.78e-17 0 0 0 0 0
0.1595 -0.1081 1.202e-26 4.128e-18 -8.33e-17 9.034e-05 0 0 0 0 0
0 0 0 0 0 0 1e-06 0 0 0 0
0 0 0 0 0 0 0 1e-07 0 0 0
0 0 0 0 0 0 0 0 1e-07 0 0
0 0 0 0 0 0 0 0 0 1e-07 0
0 0 0 0 0 0 0 0 0 0 1.884e-09


will produce the result (denoted R) :

Code: Select all

31.61 0 0 0 0 0 0 0 0 0 0
-0.02811 31.6 0 0 0 0 0 0 0 0 0
-9.823e-26 -2.296e-25 0.0003162 0 0 0 0 0 0 0 0
-3.374e-17 -7.887e-17 -6.677e-08 0.02734 0 0 0 0 0 0 0
-1.764e-17 -2.807e-17 1.015e-08 0.04033 0 0 0 0 0 0 0
0.005045 -0.003416 -2.956e-22 -1.174e-15 0 0.007295 0 0 0 0 0
0 0 0 0 31.62 0 0 0 0 0 0
0 0 0 0 0 0 0 0.0003162 0 0 0
0 0 0 0 -9.309e-28 0 0.08944 0 0 0 0
0 0 0 0 -8.231e-28 0 3.72e-28 0 0 0.0003162 0
0 0 0 0 -1.332e-23 0 1.876e-24 0 0.418 -3.687e-47 0


Note the bottom left corner, these values may be even smaller or also very high. Obviously, R * R.transpose() will not produce the original matrix. I originally found this to be a problem of elements in the zero blocks (I always use block matrices), and tried setting them to zero, but that will produce similar errors near the diagonal in the long term. I also tried using the LDLT, but it seems LDLT directly returns NaN values in these cases.

I wonder, is this an error in the Cholesky module, or are my matrices somehow ill-conditioned? A symmetry issue indeed?

Thank you in advance.

Some important fixes just made it into the Cholesky module, including in the 2.0 branch. Can you check out the latest 2.0 branch (instructions on the main page of the wiki) and tell us if it works better for you?

EDIT: oops, you said LLt, the recent fixes were in LDLt. I'll have a look at your stuff.

OK, now that I realized that you're using LLt, not LDLt, the reason for your trouble is clear: your matrix is almost singular, so LLt is going to be very unstable on it. There's no way around that. My best recommendation would be to try LDLt instead, and for that, update to the latest 2.0 branch.

Thanks, I'll try using the LDLt instead.

By "latest branch", do you mean the development branch or the 2.0.12 release?

By "latest 2.0 branch" i mean the 2.0 branch, not the development branch. It's what's going to be released as 2.0.13 eventually.