SimplicialLLT returns wrong cholesky factors

Question

I want to use Eigen to compute the cholesky decomposition of a sparse matrix. However, the result is incorrect and I am not able to find a reason for it. How do I get the correct answer?

And are there special routines implemented in Eigen which exploit the structure of sparse matrices in order to boost performance (e.g. for banded matrices as in the example below or triangular matrices)?

#include <iostream>
#include <Eigen/Sparse>
#include <Eigen/Dense>

int main() 
{

    // create sparse Matrix
    int n = 5; 
    std::vector<Eigen::Triplet<double> > ijv; 
    for(int i = 0; i < n; i++)
    {
        ijv.push_back(Eigen::Triplet<double>(i,i,1));
        if(i < n-1)
        {
            ijv.push_back(Eigen::Triplet<double>(i+1,i,-0.9));
        }
    }
    Eigen::SparseMatrix<double> X(n,n);
    X.setFromTriplets(ijv.begin(), ijv.end());
    Eigen::SparseMatrix<double> XX = X * X.transpose();

    // Cholesky decomposition 
    Eigen::SimplicialLLT <Eigen::SparseMatrix<double> > cholesky;
    cholesky.analyzePattern(XX);
    cholesky.factorize(XX);

    std::cout << Eigen::MatrixXd(XX) << std::endl;
    std::cout << Eigen::MatrixXd(cholesky.matrixL()) << std::endl;

}

The matrices look as follows:

Input XX:

   1 -0.9    0    0    0
-0.9 1.81 -0.9    0    0
   0 -0.9 1.81 -0.9    0
   0    0 -0.9 1.81 -0.9
   0    0    0 -0.9 1.81

Output (cholesky.matrixL()):

  1.34536         0         0         0         0
-0.668965   1.16726         0         0         0
        0 -0.771039    1.1025         0         0
        0         0         0         1         0
        0         0 -0.816329      -0.9  0.577587

What it should look like (X):

   1    0    0    0   0
-0.9    1    0    0   0
   0 -0.9    1    0   0
   0    0 -0.9    1   0
   0    0    0 -0.9   1

Answer 1

Don't forget that SimplicialLLT does not factorize A = L * L^T but P * A * P^T = L * T^T with P a permutation matrix. If you need P to be the identity then use NaturalOrdering:

Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::NaturalOrdering<int> > cholesky;