C++ Cholesky 分解答案

【问题标题】：C++ Cholesky factorizationC++ Cholesky 分解
【发布时间】：2020-05-29 05:18:32
【问题描述】：

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

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

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

我尝试使用 Eigen 的 LDLT 类：

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

【问题讨论】：

那是因为“matB << matA, matA/10, matA/10, matA;”并没有像你想象的那样做。这里唯一发生的是matB << MatA，其他所有内容都被完全忽略了。如果您不知道如何在 C++ 中做某事，那么在编译之前尝试随机语法不太可能产生预期的结果。附言与“matA << 1, 2, 3, 4;”相同，在我看来就像you could use a good C++ book。
@SamVarshavchik 根据 Eigen 的非常基本的示例，它是非常好的语法。 eigen.tuxfamily.org/dox/…
@SamVarshavchik 你是想说 Eigen 做错事了，你应该避免 Eigen 吗？但是，这根本不能帮助我解决我的问题。如果您愿意，您可以考虑以不同的方式进行初始化。
@user7431005 Eigen 在重载operator, 时非常好。 Sam 引用的答案实际上是说（加上comment）如果你不知道自己在做什么，那么你不应该超载operator,。
对于那些关注运算符重载的人来说，Eigen 做得很好。事实上，here 证明了 OP 认为计算的矩阵与问题后面的矩阵相匹配。

标签： c++ matlab eigen

【解决方案1】：

使用您的代码示例和Matlab documentation，我在使用 LLT 而不是 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

成功

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

【讨论】：

【解决方案2】：

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

所以你的类实例化应该和这个类似（不知道这里正确的 Eigen-define 是不是Upper）：

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

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

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

【讨论】：

_UpLo 参数描述了输入矩阵的哪一半正在考虑。您可以使用.matrixL() 或.matrixU() 输出任一因子（尽管这些只是对相同数据的不同视图）。正如@Avi 指出的那样，正确的答案是使用LLT 而不是LDLT。