Reply to topic

SimplicialLDLT for symmetric positive semi definite

n2k
Registered Member
Posts
41
Karma
0
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!
User avatar ggael
Moderator
Posts
3447
Karma
19
OS
hm, I think it should work too. You can still test solver.info()==Success after the factorization, and fall back to BiCGSTAB if needed.
User avatar alecjacobson
Registered Member
Posts
26
Karma
0
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.
User avatar ggael
Moderator
Posts
3447
Karma
19
OS
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.
User avatar ggael
Moderator
Posts
3447
Karma
19
OS
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
matthieun
Registered Member
Posts
5
Karma
0
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.
User avatar ggael
Moderator
Posts
3447
Karma
19
OS
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.

 
Reply to topic

Bookmarks



Who is online

Registered users: Bing [Bot], Felix Ernst, Google [Bot], markhm, mrperfecttt, shyamns, Sogou [Bot], Yahoo [Bot], zaydplays