Dear all,

I am trying to update Cholesky decomposition using Eigen. In the problem Ax = b where A is symmetric positive definite (so we have A = L*L^T), I want to update the lower triangular matrix L in the simplest case when the structure of A (and L) is not changed. For example I change some values a_ij = a_ij + delta_ij and a_ji = a_ji + delta_ij if a_ij != 0. New matrix A is still SPD, so it is very similar with rank-one update problem.

My question is: does Eigen support update Cholesky decomposition? Assume that I can recompute matrix L using my code, is there any way to input it directly into a SimplicialLLT solver instead of calling the solver's compute()? If it is not supported yet in Eigen, could you please suggest me alternative C++ libraries which support update Cholesky decomposition like MATLAB? Thanks.

For dense factorization, LLT and LDLT supports rank-k updates. However, this is not the case yet for sparse matrices (SimplicialLLT). You can by-pass compute() as follows:

SimplicialLLT<SparseMatrix<double>,Lower> llt(A);
llt.matrixL().nestedExpression().const_cast_derived().coeffRef(i,j) = ...

However, beaware that by default, SimplicialLLT applies a symmetric permutation to reduce fill-in. You can get it with llt.permutationP(). If A is already well ordered with respect to fill-in, you can also disable permutation with:

SimplicialLLT<SparseMatrix<double>,Lower,NaturalOrdering<int> >

btw, Cholmod supports rank-one updates of the form A += v * v';

ggael wrote:For dense factorization, LLT and LDLT supports rank-k updates. However, this is not the case yet for sparse matrices (SimplicialLLT). You can by-pass compute() as follows:

SimplicialLLT<SparseMatrix<double>,Lower> llt(A);
llt.matrixL().nestedExpression().const_cast_derived().coeffRef(i,j) = ...

However, beaware that by default, SimplicialLLT applies a symmetric permutation to reduce fill-in. You can get it with llt.permutationP(). If A is already well ordered with respect to fill-in, you can also disable permutation with:

SimplicialLLT<SparseMatrix<double>,Lower,NaturalOrdering<int> >

Hi ggael,

Could you tell me which ordering algorithm is used for the permutation in SimpliciallLLT?

The default is AMD (Approximate Minimum Degree), as explained in the doc: http://eigen.tuxfamily.org/dox-devel/cl ... alLLT.html