Hi all,
In the course of my research, I need to compute the Cholesky decomposition of a covariance matrix (to draw joint multivariate normal samples).
I came across a strange phenomenon that I cannot quite explain: when computing the decomposition for certain matrix (symmetric definite positive and small), the decomposition is not accurate.
Here is a sample of the code that results in strange output (I tried to simplify it to the point where the error occurs):

Code: Select all

#include <Eigen/Eigenvalues>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include <time.h>
#include <string>
#include <iostream>

using Eigen::MatrixXd;
using Eigen::VectorXd;
using Eigen::SparseMatrix;
typedef Eigen::Triplet<double> T;

int main(int argc, char *argv[])
{
    // DENSE

    MatrixXd sigma_dense = MatrixXd::Zero(8,8);

    sigma_dense << 1,0.940222,0.940222,0,0,0,0,0,
                   0.940222,1,0,0.940222,0,0,0,0,
                   0.940222,0,1,0.940222,0,0,0,0,
                   0,0.940222,0.940222,1,0,0,0,0,
                   0,0,0,0,1,0.940222,0.940222,0,
                   0,0,0,0,0.940222,1,0,0.940222,
                   0,0,0,0,0.940222,0,1,0.940222,
                   0,0,0,0,0,0.940222,0.940222,1;

    std::cout<<"sigma = "<<std::endl;
    std::cout<<sigma_dense<<std::endl;
    std::cout<<"Determinant = "<<sigma_dense.determinant()<<std::endl;

    // Cholesky of the Covariance matrix

    Eigen::LLT<MatrixXd> llt(sigma_dense);
    MatrixXd L = llt.matrixL();
    MatrixXd U = llt.matrixU();

    std::cout<<"L*Lt = "<<std::endl;
    std::cout<< L*U<<std::endl;

    MatrixXd a = sigma_dense - L*U;
    std::cout<<"Error = "<<sqrt((a*a.transpose()).trace())<<std::endl;

    // SPARSE

    SparseMatrix<double> sigma(8,8);

    std::vector<T> tripletList;

    tripletList.push_back(T(0,0,1));
    tripletList.push_back(T(1,1,1));
    tripletList.push_back(T(2,2,1));
    tripletList.push_back(T(3,3,1));
    tripletList.push_back(T(4,4,1));
    tripletList.push_back(T(5,5,1));
    tripletList.push_back(T(6,6,1));
    tripletList.push_back(T(7,7,1));

    tripletList.push_back(T(0,1,0.940222));
    tripletList.push_back(T(2,3,0.940222));
    tripletList.push_back(T(4,5,0.940222));
    tripletList.push_back(T(6,7,0.940222));

    tripletList.push_back(T(1,0,0.940222));
    tripletList.push_back(T(3,2,0.940222));
    tripletList.push_back(T(5,4,0.940222));
    tripletList.push_back(T(7,6,0.940222));

    tripletList.push_back(T(0,2,0.940222));
    tripletList.push_back(T(1,3,0.940222));
    tripletList.push_back(T(4,6,0.940222));
    tripletList.push_back(T(5,7,0.940222));

    tripletList.push_back(T(2,0,0.940222));
    tripletList.push_back(T(3,1,0.940222));
    tripletList.push_back(T(6,4,0.940222));
    tripletList.push_back(T(7,5,0.940222));

    sigma.setFromTriplets(tripletList.begin(), tripletList.end());

    std::cout<<"sigma = "<<std::endl;
    std::cout<<sigma<<std::endl;

    // Cholesky of the Covariance matrix

    Eigen::SimplicialLLT< SparseMatrix<double> > llt_sparse(sigma);
    SparseMatrix<double> L_sparse = llt_sparse.matrixL();
    SparseMatrix<double> U_sparse = llt_sparse.matrixU();

    std::cout<<"L*Lt = "<<std::endl;
    std::cout<< L_sparse*U_sparse<<std::endl;

    MatrixXd a_sparse = sigma - L_sparse*U_sparse;
    std::cout<<"Error = "<<(a_sparse*a_sparse.transpose()).trace()<<std::endl;

    return 0;
}


This gives me the following output:

sigma =
1 0.940222 0.940222 0 0 0 0 0
0.940222 1 0 0.940222 0 0 0 0
0.940222 0 1 0.940222 0 0 0 0
0 0.940222 0.940222 1 0 0 0 0
0 0 0 0 1 0.940222 0.940222 0
0 0 0 0 0.940222 1 0 0.940222
0 0 0 0 0.940222 0 1 0.940222
0 0 0 0 0 0.940222 0.940222 1

Determinant = 6.43165

L*Lt =
1 0.940222 0.940222 0 0 0 0 0
0.940222 1 0 0.940222 0 0 0 0
0.940222 0 8.62198 -6.22613 0 0 0 0
0 0.940222 -6.22613 9.506 0 0 0 0
0 0 0 0 1 0.940222 0.940222 0
0 0 0 0 0.940222 1.88402 0.884017 0.940222
0 0 0 0 0.940222 0.884017 1.88402 0.940222
0 0 0 0 0 0.940222 0.940222 2.76803

Error = 15.4729

sigma =
1 0.940222 0.940222 0 0 0 0 0
0.940222 1 0 0.940222 0 0 0 0
0.940222 0 1 0.940222 0 0 0 0
0 0.940222 0.940222 1 0 0 0 0
0 0 0 0 1 0.940222 0.940222 0
0 0 0 0 0.940222 1 0 0.940222
0 0 0 0 0.940222 0 1 0.940222
0 0 0 0 0 0.940222 0.940222 1

Segmentation fault (core dumped)

=====> The decomposition using the standard Cholesky gives me an l2 error of 15 and the sparse decomposition do not seem to work
=====> In my experiment, if the matrix is not too sparse, then the sparse solvers works fine

Any suggestions would be greatly appreciated,

Thanks,

Ok I am going to answer my own questions with another question,
The matrix that I am using here:
1 0.940222 0.940222 0 0 0 0 0
0.940222 1 0 0.940222 0 0 0 0
0.940222 0 1 0.940222 0 0 0 0
0 0.940222 0.940222 1 0 0 0 0
0 0 0 0 1 0.940222 0.940222 0
0 0 0 0 0.940222 1 0 0.940222
0 0 0 0 0.940222 0 1 0.940222
0 0 0 0 0 0.940222 0.940222 1
is a thresholded covariance matrix. I was probably too harsh on the threshold (in order to replicate the error on small size) and this matrix is actually not positive with 2 clearly negative eigenvalues.
However, my problem still occurs with larger matrices (10611x10611). Specifically the reconstruction error of the Cholesky decomposition is very large (The l2 error is roughly 10^7, which means a mean error of 10^7 / 10611^2 = 0.2).
Could it be the case that if the eigenvalues get closer to 0, the Cholesky decomposition ends up performing very poorly ?

How do you know that your matrix 10 000*10 000 is positive symmetric and definite ?
Checking the symmetry is quite easy. But checking if such a big matrix is positive and definite is not that easy.
So how can you be sure ?

The matrix is the covariance matrix of a point cloud using a Gaussian kernel.
Specifically, for 2 points x and y, the matrix C is defined as C(x,y) = exp( - d(x,y)^2 / w^2) with d(.,.) the euclidean norm and w the scale parameter.
Theoretically, it should be positive definite. In practice, when the number of points grows large, a lot of eigenvalues are extremely small.
My guess is that the Cholesky decomposition do not perform so well in that case. An alternative being to use the eigen solver to compute C = UDU^(-1), threshold the extremelly small eigenvalues and use L = U sqrt(D). But this takes much more time.

A very basic, and probably not so useful idea. Did you try to multiply all the values of your matrix by a large coefficient so that your very small eiganvalues become larger ?
I guess it might not work since the ratio between eigenvalues will still be large, but it does not take much to try it.

So I tried this idea and it does not improve the reconstruction error.

In my case, the matrix 10000x10000 which should theoretically by positive definite numerically has very small negative eigenvalues: the largest eigenvalue is around 50 and the smallest is around -10^(-17).

I tried the different algorithms available and my conclusion is that the Cholesky decomposition (LLt, LDLt and the incomplete Cholesky decomposition from the unsupported eigen) do not perform well. A solution is to compute the actual svd UDU^(-1) and threshold the negative eigenvalue to zero.

The main drawback is that I lose a lot in computational efficiency by doing so.

For nearly rank deficient problems, PartialPivLU or HouseholderQR will perform much better than Cholesky-based methods while being much fatser than SVD. If the matrix is sparse, try the SparseLU and SparseQR variants.