Hi,

I just have a simple question. In the documentation of Eigen 3.1, it says that, for dense matrices:
LDLT requires a SPsD matrix (http://eigen.tuxfamily.org/dox/TopicLin ... tions.html),
while, for the sparse version,
SimplicialLDLT requires a SPD matrix (http://eigen.tuxfamily.org/dox/TutorialSparse.html).

Does the sparse version of LDLT requires strict SPD? So, if I have a sparse SPsD matrix, the only built-in option is BiCGSTAB right?

Thanks!

hm, I think it should work too. You can still test solver.info()==Success after the factorization, and fall back to BiCGSTAB if needed.

It seems like SimplicialLDLT is not working for symmetric positive semi definite input. Here's a small test program:

Code: Select all


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

int main(int argc,char * argv[])
{
  using namespace Eigen;
  using namespace std;

  MatrixXd NAf(5,5);
  NAf<<17,4,-20,0,1,4,42,-10,-20,0,-20,-10,42,4,0,0,-20,4,17,1,1,0,0,1,0;
  LDLT<MatrixXd > ldlt_solver;
  ldlt_solver.compute(NAf);
  switch(ldlt_solver.info())
  {
    case Eigen::Success:
      break;
    case Eigen::NumericalIssue:
      cerr<<"Error: LDLT::info(): Numerical issue."<<endl;
    default:
      return 1;
  }

  SparseMatrix<double> NA(5,5);
  vector<Triplet<double> > ijv;
  for(int i = 0;i<5;i++)
  {
    for(int j = 0;j<5;j++)
    {
      // Commenting out this zero check will cause SimplicialLDLT to succeed
      // (but be dense)
      if(NAf(i,j) != 0)
      {
        ijv.push_back(Triplet<double>(i,j,NAf(i,j)));
      }
    }
  }
  NA.setFromTriplets(ijv.begin(),ijv.end());
  SimplicialLDLT<SparseMatrix<double> > simplicial_ldlt_solver;
  simplicial_ldlt_solver.compute(NA);
  switch(simplicial_ldlt_solver.info())
  {
    case Eigen::Success:
      cerr<<"Error: SimplicialLDLT::info(): success."<<endl;
      break;
    case Eigen::NumericalIssue:
      cerr<<"Error: SimplicialLDLT::info(): Numerical issue."<<endl;
    default:
      return 1;
  }
  return 0;
}


The hard-coded matrix is positive semi-definite (tested in matlab, too). The dense version of LDLT succeeds but the sparse SimplicialLDLT fails. Interesting if I create a "dense" sparse matrix (insert 0 elements) then SimplicialLDLT succeeds.

Hope this helps track down this issue. It would be great to have an LDLT that works on positive semi-definite sparse matrices.

There is a major difference between the dense and sparse LDL^T factorization: in the dense version we perform numerical pivoting while the sparse version performs fill-in reducing pivoting. It seems that for this particular example the fill-in pivoting does a bad job on the numerical side that is why disabling fill-in pivoting make it work, but that's just luck. For another matrix it might be the opposite. As a workaround, the SparseLU solver does numerical column pivoting.

BTW, note that you can convert NAF to NA with NA = NAF.sparseView();

EDIT: after fill-in ordering the matrix looks like:

Code: Select all

  0   1   0   0   1
  1  17   4 -20   0
  0   4  42 -10 -20
  0 -20 -10  42   4
  1   0 -20   4  17


which is pretty bad! There must exist a way to make the fill-in ordering numerically more robust.

For the record, a few months back, I've improved this aspect by making AMD reordering aware of structural zeros: https://bitbucket.org/eigen/eigen/commits/ce076b6d66ce

The example given above by alecjacobson does work now.

There is also a respective bug entry: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=751

Hi everyone,

Gaël, does it mean that SimplicialLDLT should now (version 3.2.5) work on any symmetric positive semi-definite matrix?
I just tried and got an error...
I took a symmetric positive definite matrix (for which the factorization is working) and I cleared a row and a col to make it semi-definite. Is it a (bad) special case?

Thanks.

No, it still does not work on problems with null eigenvalues, but it now properly handles problems of the form:

Code: Select all


| A  C^T |
| C   0  |


with A SPD, and C of rank row(C). The example provided by Alec is of this form, but such problems are usually not positive-semidefinite as they have both positive and negative eigenvalues, but they are invertible.