Hi,
I wasn't able to find clear documentation on LDLT. I was hoping you could help me please. In this code I'm trying to find the squre root of a positive semidefinite matrix:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
   MatrixXd A(2, 2);
   A << 2, -1, -1, 3;
   cout << "Here is the matrix A:\n" << A << endl;
   
   LDLT<MatrixXd> ldltOfA(A);
   VectorXd d = ldltOfA.vectorD();
   d = d.array().sqrt();
   MatrixXd L;
   // L = ldltOfA.matrixL()*d.asDiagonal(); // fails to compile
   // L = ldltOfA.matrixL().matrix()*d.asDiagonal(); // fails to compile
   MatrixXd Ll = ldltOfA.matrixL();
   L = ldltOfA.transpositionsP().inverse()*Ll*d.asDiagonal();
   
   cout << "Here is the matrix L:\n" << L << endl;
   
   // MatrixXd P =  ldltOfA.transpositionsP(); // fails to compile
}


These lines fails to compile:

Code: Select all

   // L = ldltOfA.matrixL()*d.asDiagonal(); // fails to compile
   // L = ldltOfA.matrixL().matrix()*d.asDiagonal(); // fails to compile


The operator * is undefined on that object. I had to assign it to a temporary matrix first. Is that expected bahaviour?

Also, I wondered why LDLT provides the method transpositionsP() of type TranspositionType rather than a PermutationMatrix type. I assume that I can still use it as a regular a matrix in this line

Code: Select all

   L = ldltOfA.transpositionsP().inverse()*Ll*d.asDiagonal();


Is there a way to get an actual matric object from it? Like this...

Code: Select all

MatrixXd P =  ldltOfA.transpositionsP(); // fails to compile


Finally, what's the difference between the methods transpose() and inverse() for type TranspositionType?

Is there a build in to find the sqrt of a PSD matrix??

Thanks!
Raman

For teh sqrt matrix operator look at this function:

http://eigen.tuxfamily.org/dox/classEig ... 8b2d2b0775

Regarding your issue, triangular times diagonal is not supported. You need to copy the L factor to a full matrix first:

Code: Select all

MatrixXd L = ldlt.matrixL();

You can copy a Transpositions object to a PermutationMatrix through ctor or operator=:

Code: Select all

MatrixXd P = PermutationMatrix<Dynamic,Dynamic>(ldlt.transpositionsP());

but why would you sore a permutation matrix into a dense matrix?

Hi ggael,

Thanks a lot. I didn't need a dense representation of P. I only wanted to print (cout) the permutation matrix our some representation of the permutation. So is there a way to do that. I guess I could multiply it to an identity matrix and print the result.

Is there a difference between inverse and transpose for the TranspositionType?

thanks!
Raman

The solution I showed you is faster than multiplying by the identity.

transpose and inverse are obviously the same for a permutation matrix. having both permits to write generic code.

And, again don't miss the operatorSqrt function that does what you're looking for.