C++ Cholesky 分解

Question

我需要用 C++ 重写一些 MatLab 代码。

在 Matlab 代码中，我们调用函数 chol 来计算上三角矩阵。

对于 C++ 部分，我开始研究 Eigen。但是，我很难获得与 Matlab 的 chol 函数等效的函数。

我尝试使用 Eigen 的 LDLT class：

#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;
}

但结果与 Matlab 代码不同：

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

根据 Eigen::LDTL 的文档，第二个模板参数 _UpLo 默认为 Lower 但您省略了该参数并想要计算上三角矩阵。

因此您的 class 实例化应该类似于此（不知道此处正确的 Eigen-define 是否是 Upper）：

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

编辑 @chtz 是对的-使用 Upper 不会给您预期的结果，因为 LDLT class 用于通过旋转进行稳健的 cholesky 分解。因此，除了@Avi 的正确答案之外，您还可以使用正确的 class 进行标准的 cholesky 分解：

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

Answer 2

使用您的代码示例和 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";
}

输出：

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

C++ Cholesky 分解

C++ Cholesky factorization

c++

matlab

eigen