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