对于矩阵乘积的Cholesky,Eigen和C++ 11类型推断失败

Question

我试图使用Eigen和C++ 11"auto"类型对其转置矩阵的乘积进行cholesky分解.我试图做的时候出现问题

auto c = a * b
auto cTc = c.tranpose() * c;
auto chol = cTc.llt();

我正在使用XCode 6.1,Eigen 3.2.2.我得到的类型错误在这里.

这个最小的例子显示了我的机器上的问题.更改的类型c从auto到MatrixXd看到它的工作.

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


int main(int argc, const char * argv[]) {
    MatrixXd a = MatrixXd::Random(100, 3);
    MatrixXd b = MatrixXd::Random(3, 100);
    auto c = a * b;
    auto cTc = c.transpose() * c;
    auto chol = cTc.llt();
    return 0;
}

有没有办法让这个工作仍然使用自动？

作为一个侧面问题,是否有一个表现理由不断言矩阵是MatrixXd在每个阶段？使用auto可以让Eigen将答案保持为它所想象的任何奇怪的模板表达式.我不确定是否输入MatrixXd会导致问题.

Answer 1

问题是第一次乘法返回a Eigen::GeneralProduct而不是a MatrixXd并且auto正在拾取返回类型.您可以隐式地创建一个MatrixXd从Eigen::GeneralProduct,所以当你明确声明它正常工作的类型.见http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

编辑:我不是铸造的特征产品或性能特征方面的专家.我刚从错误消息中推测出答案并从在线文档中确认.分析始终是检查代码不同部分性能的最佳选择.

Answer 2

让我总结一下发生了什么以及为什么是错误的。首先，让我们auto用它们所采用的类型来实例化关键字：

typedef GeneralProduct<MatrixXd,MatrixXd> Prod;
Prod c = a * b;
GeneralProduct<Transpose<Prod>,Prod> cTc = c.transpose() * c;

请注意，Eigen 是一个表达式模板库。这里，GeneralProduct<>是代表产品的抽象类型。不执行任何计算。因此，如果您复制cTc到MatrixXdas：

MatrixXd d = cTc;

这相当于：

MatrixXd d = c.transpose() * c;

那么产品a*b将会进行两次！因此，无论如何，最好a*b在显式临时内进行评估，对于以下情况也是如此c^T*c：

MatrixXd c = a * b;
MatrixXd cTc = c.transpose() * c;

最后一行：

auto chol = cTc.llt();

也是相当错误的。如果 cTc 是抽象产品类型，那么它会尝试实例化在抽象产品类型上工作的 Cholesky 分解，这是不可能的。现在，如果cTc是 a MatrixXd，那么您的代码应该可以工作，但这仍然不是首选方式，因为该方法llt()而是实现单行表达式，例如：

VectorXd b = ...;
VectorXd x = cTc.llt().solve(b);

如果您想要一个命名的 Cholesky 对象，请使用它的构造函数：

LLT<MatrixXd> chol(cTc);

甚至：

LLT chol(c.transpose() * c);

除非您必须c.transpose() * c在其他计算中使用，否则这是等效的。

最后，根据a和的大小b，最好计算cTc为：

MatrixXd cTc = b.transpose() * (a.transpose() * a) * b;

在未来（即 Eigen 3.3），Eigen 将能够看到：

auto c = a * b;
MatrixXd cTc = c.transpose() * c;

作为四个矩阵的乘积m0.transpose() * m1.transpose() * m2 * m3并将括号放在正确的位置。但是，它无法知道m0==m3和，因此如果首选方法是临时m1==m2评估，那么您仍然需要自己进行评估。a*b