8

我正在尝试使用 Eigen 和 C++11“自动”类型对矩阵乘积及其转置进行 cholesky 分解。当我尝试做的时候问题就来了

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

我正在使用 XCode 6.1，Eigen 3.2.2。我得到的类型错误是here。

这个最小的例子显示了我机器上的问题。c将from的类型更改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在每个阶段都是 a？使用 auto 将允许 Eigen 将答案保留为它喜欢的任何奇怪的模板表达式。我不确定将其键入为 MatrixXd 是否会导致问题。

4

4 回答 4

6

问题是第一次乘法返回 aEigen::GeneralProduct而不是 aMatrixXd并且auto正在选择返回类型。MatrixXd您可以从 a隐式创建 a ，Eigen::GeneralProduct因此当您显式声明类型时它可以正常工作。见http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

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

于 2014-11-24T20:45:55.527 回答
4

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

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

请注意，Eigen 是一个表达式模板库。这里，GeneralProduct<>是代表产品的抽象类型。不执行任何计算。因此，如果您复制cTc到 a 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临时评估，那么您仍然需要自己进行。

于 2014-11-25T08:57:23.700 回答
2

我不是这方面的专家Eigen，但是像这样的库通常会从操作中返回代理对象，然后使用隐式转换或构造函数来强制执行实际工作。（表达式模板是一个极端的例子。）这避免了在许多情况下不必要地复制完整的数据矩阵。不幸的是，auto很高兴只创建一个代理类型的对象，通常库的用户永远不会明确声明。由于您最终需要计算数字，因此从强制转换到MatrixXd. （Scott Meyers 在“Effective Modern C++”中给出了使用autowith的相关示例vector<bool>，其中operator[](size_t i)返回一个代理。）

于 2014-11-24T21:21:59.540 回答
0

不要auto与 Eigen 表达式一起使用。我之前遇到过更“戏剧性”的问题，请参阅

一般产品的本征自动类型扣除

并被一位 Eigen 创造者 (Gael) 建议不要auto与 Eigen 表达式一起使用。

从表达式到特定类型的转换MatrixXd应该非常快，除非您想要延迟评估（因为在进行转换时会评估结果）。

于 2014-11-24T23:09:46.737 回答