Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

note that the LDLT decomp does pivoting, so it is normal that you don't get the same diagonal matrix. If you reconstruct the matrix from the decomp you will see that in this case the decomp is correct.


In the devel branch I've just added a reconstructedMatrix() to LDLT for that precise purpose, so you can do:

cout << (A - A_factored.reconstructedMatrix()) << endl;


As a side note,


A.diagonal().setOnes();

works as expected ;) (simpler than a home made for loop)

gael.

On Wed, Feb 24, 2010 at 3:29 AM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote:

Hi,

I have found another robustness problem with ldlt() in the current
development branch. Some combination of rank-deficiency, constant
diagonal, pivoting, and double precision can make it impossible to
roundtrip a matrix through the LDLt decomposition. I am attaching
another test case that yields the output below.

Thank you,
Ben

A = LDL' has a unit diagonal (apart from numerical noise)
1
1
1
1
1
This is a numerical disaster (compare the second cells).
Here is the true D
1
0.915092
0
0
0
Here is the calculated D
1
0.0785762
2.01553e-16
7.97973e-17
1.39591e-15
This is good (because the diagonal of A is forced to have 1.0 in all its cells).
Here is the true D
1
0.915092
0
0
0
Here is the calculated D
1
0.915092
-6.06069e-16
-8.0296e-16
-8.4893e-16