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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

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;
    }
}