Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices
• From: Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>
• Date: Tue, 23 Feb 2010 21:29:10 -0500

```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
```
```#include "eigen/Eigen/Cholesky"
#include<iostream>

using namespace Eigen;
using namespace std;

extern "C" {

int main() {
LDLT<MatrixXd> A_factored;
VectorXd D(5);
D << 1.0, 0.915092047236159,                  0.0, 0.0, 0.0;
MatrixXd L(5,5);
L << 1.0,                 0.0,                0.0, 0.0, 0.0,

0.2913896922745232,  1.0,                0.0, 0.0, 0.0,

0.3506763450086252,  0.9789801109802591, 1.0, 0.0, 0.0,

-0.0551313984496821,  1.0437742813828543, 0.0, 1.0, 0.0,

-0.2269670280948324, -1.0180827163703894, 0.0, 0.0, 1.0;

MatrixXd A = L * D.asDiagonal() * L.transpose();
A_factored = A.ldlt();

cout << "A = LDL' has a unit diagonal (apart from numerical noise)" << endl << A.diagonal() << endl;

cout << "This is a numerical disaster (compare the second cells)." << endl;
cout << "Here is the true D" << endl << D << endl << "Here is the calculated D" <<
endl << A_factored.vectorD() << endl;

for(int i=0; i < 5; i++) A(i,i) = 1.0;
A_factored = A.ldlt();

cout << "This is good (because the diagonal of A is forced to have 1.0 in all its cells)." << endl;
cout << "Here is the true D" << endl << D << endl << "Here is the calculated D" <<
endl << A_factored.vectorD() << endl;

return 0;
}
}

```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/