C++ Cholesky factorization

Question

I need to rewrite some MatLab code using C++.

Inside the Matlab code, we are calling the function chol to calculate an upper triangular matrix.

For the C++ part, I'm starting to look at Eigen. However, I struggle to get an equivalent of Matlab's chol function.

I tried to use the LDLT class of Eigen:

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

int main() {

  MatrixXd matA(2, 2);
  matA << 1, 2, 3, 4;

  MatrixXd matB(4, 4);
  matB << matA, matA/10, matA/10, matA;
  matB = matB*matB.transpose();

  Eigen::LDLT<MatrixXd> tmp(matB);
  MatrixXd U = tmp.matrixU();
  cout << U << endl;
}

but the result is different compared to the Matlab code:

matB = [  1   2 0.1 0.2
          3   4 0.3 0.4
        0.1 0.2   1   2
        0.3 0.4   3   4];
matB = matB*matB';
D = chol(matB);

Answer 1

Using your code example and the Matlab documentation, I get the same result when using LLT instead of LDLT (online):

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using std::cout;

int main()
{
  MatrixXd matA(3,3);
  matA << 1, 0, 1, 0, 2, 0, 1, 0, 3;
  cout << matA << "\n\n";
  Eigen::LDLT<MatrixXd> tmp(matA);
  cout << ((tmp.info() == Success) ? "succeeded" : "failed") << "\n\n";
  MatrixXd U = tmp.matrixL();
  cout << U << "\n\n";
  // Using LLT instead
  cout << MatrixXd(matA.llt().matrixL()) << "\n\n";
  cout << MatrixXd(matA.llt().matrixU()) << "\n\n";
}

Outputs:

1 0 1
0 2 0
1 0 3

succeeded

1 0 0
0 1 0
0.333333 0 1

1 0 0
0 1.41421 0
1 0 1.41421

1 0 1
0 1.41421 0
0 0 1.41421