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);

That's because "`matB << matA, matA/10, matA/10, matA;`" does not do what you think it does. The only thing that happens here is `matB << MatA`, and everything else gets completely ignored. If you do not know how to do something in C++, just trying random syntax until it compiles is unlikely to produce the expected results. P.S. Same thing with "`matA << 1, 2, 3, 4;`", looks to me like [you could use a good C++ book](https://stackoverflow.com/questions/388242/the-definitive-c-book-guide-and-list). — Sam Varshavchik, Feb 13 '20 at 13:25
@SamVarshavchik no need to be offensive. Please provide more information. If I execute `cout << matB` it looks perfectly fine. — user7431005, Feb 13 '20 at 13:29
Now try `cout << matA, matB` and see if you note any difference. That's my entire point. — Sam Varshavchik, Feb 13 '20 at 13:29
@SamVarshavchik According to the very basic examples of Eigen it is perfectly fine syntax. https://eigen.tuxfamily.org/dox/group__TutorialAdvancedInitialization.html — user7431005, Feb 13 '20 at 13:30
But according to people who know C++, [the comma operator should not be overloaded](https://stackoverflow.com/questions/5602112/when-to-overload-the-comma-operator), so it's doing something that has so many pitfalls, that it's hard to avoid them. — Sam Varshavchik, Feb 13 '20 at 13:34
@SamVarshavchik do you want to say that Eigen is doing things wrong and you should avoid Eigen? However, this does not help me with my problem at all. You can think of the initialization done differently if you wish. — user7431005, Feb 13 '20 at 13:36
@user7431005 Eigen is perfectly fine in overloading `operator,`. The answer Sam is quoting actually says (plus the [comment](https://stackoverflow.com/questions/5602112/when-to-overload-the-comma-operator#comment93931744_5602136)) that if you don't know what you're doing, then you shouldn't overload `operator,`. — Avi Ginsburg, Feb 13 '20 at 13:40
To those who are focusing on the operator overloading, Eigen did it fine. In fact, [here](https://godbolt.org/z/eQj2kX) is proof that the matrix OP thinks was computed matches the one later in the question. — Avi Ginsburg, Feb 13 '20 at 14:18

score 4 · Accepted Answer · answered Feb 13 '20 at 14:07

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

Odysseus · Answer 2 · 2020-02-13T17:33:31.780

2

According the documentation of Eigen::LDTL the 2nd template parameter _UpLo defaults to Lower but you omitted that parameter and want to calculate the upper triangular matrix.

So your class instantiation should look similar to this (don't know if the correct Eigen-define here is Upper):

Eigen::LDLT<MatrixXd, Upper> tmp(matB);

Edit @chtz is right - using Upperwont give you the result you expect because LDLT class is for robust cholesky decomposition with pivoting. So in in Addition to the correct answer of @Avi you can also use the right class for standard cholesky decomposition:

Eigen::LLT<MatrixXd> tmp(matA);
cout <<  MatrixXd(tmp.matrixU()) << "\n\n";

edited Feb 13 '20 at 17:33

answered Feb 13 '20 at 14:02

Odysseus

1,213
4
12

3

The `_UpLo` parameter describes which half of the input matrix is considering. You can output either factor using `.matrixL()` or `.matrixU()` (these are just different views to the same data though). The correct answer is to use `LLT` instead of `LDLT` as @Avi pointed out. – chtz Feb 13 '20 at 16:26

C++ Cholesky factorization

2 Answers2