I know this topic has come up a few times, but I am still a little confused about how permutations are applied to the SimplicialL(D)LT decompositions, and I cannot find a clear reference in the documentation.

Referring to the code and output that I pasted below, it appears that P'AP = LDL. That of course implies that A = PLDLP'. And we see that A and PLDLPinv are the same. So far, so good.

Now, in the documentation for the twistedBy() function, we should expect that LDL = A.twistedBy(P). However, we see that LDL=A.twistedBy(Pinv) instead.

Also, my understanding is that Cholmod has PAP' = LDL, and Eigen is the other way around. Is this correct?

In any case, I'm not quite sure what I'm doing wrong here. Can you help? I'm also open to suggestions on how to make some of the sparse matrix code more efficient.

Thanks.

Here's the code:

Code: Select all

#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/Sparse>
#include <iostream>

int main() {

  using std::cout;
  using std::endl;
  using Eigen::SparseMatrix;
  using Eigen::VectorXi;
  using Eigen::VectorXd;
  using Eigen::MatrixXd;
  using Eigen::Dynamic;
  using Eigen::Upper;
  using Eigen::Lower;
  using Eigen::SimplicialLLT;
  using Eigen::SimplicialLDLT;
  using Eigen::PermutationMatrix;

  typedef SimplicialLDLT<SparseMatrix<double> > SparseLDLT;
  Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");

  int k = 5;

  typedef SparseMatrix<double> SparseMatrixXd;
  SparseMatrixXd A;
  A.resize(k,k);
  A.reserve(VectorXi::Constant(k,2));
  for (int j=0; j<k; j++) {
    A.insert(j,j)=1.;
  }
  A.insert(3,0) = .2;
  A.insert(4,3) = -.1;
  A.insert(0,3) = A.coeff(3,0);
  A.insert(3,4) = A.coeff(4,3);
  A.makeCompressed();

  SparseLDLT LDLT(A);
  SparseMatrixXd L = LDLT.matrixL();
  VectorXd D = LDLT.vectorD();
  PermutationMatrix<Dynamic> perm = LDLT.permutationP();
  PermutationMatrix<Dynamic> permInv = LDLT.permutationPinv();
  MatrixXd P = perm.toDenseMatrix().cast<double>();
  MatrixXd Pinv = permInv.toDenseMatrix().cast<double>();
  SparseMatrixXd ApermInv;
  A.twistedBy(permInv).evalTo(ApermInv);
  SparseMatrixXd Aperm;
  A.twistedBy(perm).evalTo(Aperm);

  SparseMatrixXd LDL = L * D.asDiagonal() * L.transpose();
  MatrixXd PinvLDLP = Pinv * LDL * P;
  MatrixXd PLDLPinv = P * LDL * Pinv;
  MatrixXd PinvAP = Pinv * A * P;

  cout << "A = \n" << A.toDense().format(CleanFmt) << endl << endl;
  cout << "PLDLPinv = \n" << PLDLPinv.format(CleanFmt) << endl << endl;

  cout << "LDL = \n" << LDL.toDense().format(CleanFmt) << endl << endl;
  cout << "A.perm = \n" << Aperm.toDense().format(CleanFmt) << endl << endl;
  cout << "A.permInv = \n" << ApermInv.toDense().format(CleanFmt) << endl << endl;
 
  return(0);
}


and here's the output:

Code: Select all

A =
[   1,    0,    0,  0.2,    0]
[   0,    1,    0,    0,    0]
[   0,    0,    1,    0,    0]
[ 0.2,    0,    0,    1, -0.1]
[   0,    0,    0, -0.1,    1]

PLDLPinv =
[   1,    0,    0,  0.2,    0]
[   0,    1,    0,    0,    0]
[   0,    0,    1,    0,    0]
[ 0.2,    0,    0,    1, -0.1]
[   0,    0,    0, -0.1,    1]

LDL =
[   1,    0,    0,    0,    0]
[   0,    1,    0,    0,    0]
[   0,    0,    1,    0, -0.1]
[   0,    0,    0,    1,  0.2]
[   0,    0, -0.1,  0.2,    1]

A.perm =
[   1,  0.2,    0, -0.1,    0]
[ 0.2,    1,    0,    0,    0]
[   0,    0,    1,    0,    0]
[-0.1,    0,    0,    1,    0]
[   0,    0,    0,    0,    1]

A.permInv =
[   1,    0,    0,    0,    0]
[   0,    1,    0,    0,    0]
[   0,    0,    1,    0, -0.1]
[   0,    0,    0,    1,  0.2]
[   0,    0, -0.1,  0.2,    1]



alright, first of all there is a mistake in the documentation of twistedBy:

A.twistedBy(P) = P A P'


Now regarding Simplicial*, I think it's probably a good idea to be consistant with other factorizations, and thus factorize P A P^-1. Just in time before we release 3.1.