[eigen] Solving positive semidefinite systems

[Edward Lam <edward@xxxxxxxxxx>]

Hi,

First of all, I'd like to thank you all for the wonderful library!

In trying to convert some particular old code to use Eigen, I've run into a case 
where Eigen's Cholesky solvers don't have a away to ignore A(k,k) if it's 0 (or 
sufficiently close to 0). Instead, they either silently give me bad results or 
give up altogether (depending on if it's the sparse or dense variants).

Golub and Van Loan's Matrix Computations, chapter 4.2.8 (The semdidefinite case) 
gives a very easy modification to the regular Cholesky decomposition that allows 
it to be used for positive semidefinite matrices. All it does is skip the update 
iteration if A(k,k) is found to be 0. In practice, the implementation that I'm 
used to checks if A(k,k) is less than some epsilon, and if it is, sets the 
current diagonal entry and the rest of the row (or column) to be exactly 0 and 
then skip the rest of the loop. Then during back substitution, it checks whether 
the diagonal value is 0, and if so, then it sets the corresponding component of 
the solution to 0 instead of dividing by it. This trick is similarly used in the 
LU matrix solvers as well with the exception that we just set the 0 diagonal 
entries to 0 during factorization and simply skip the divide by 0's during back 
substitution.

Can such a modification be considered for Eigen's dense/sparse Cholesky/LU 
solvers? I have managed to work around the problem by using other Eigen solvers 
like FullPivLU (for dense matrices) or SparseQR (for sparse matrices). However, 
they're much slower.

Thanks,
-Edward