征和C++ 11的类型推断失败的矩阵产品

的乔莱斯基我试图把其转置矩阵的产品的Cholesky分解，利用特征向量和C++ 11“自动”类型。我试图做的时候出现问题征和C++ 11的类型推断失败的矩阵产品

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

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

这个小例子，显示了我的机器上的问题。将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？使用自动将允许Eigen保留答案，无论它幻想哪种奇怪的模板表达。我不确定是否将它键入MatrixXd会导致问题。

来源

2014-11-24 c0g

让我总结什么是怎么回事，为什么这是错的。首先，让我们来实例化auto关键字他们正在采取的类型：

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

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

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是一个MatrixXd，那么你的代码应该工作，但是这仍然没有方法llt()的首选方法是相当实行一班轮表达，如：

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

如果你想要一个名为乔莱斯基对象，然后而使用它的构造：

LLT<MatrixXd> chol(cTc);

甚至：

LLT CHOL（c.transpose（）* C）;

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

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

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

未来（即本征3.3），本征将能够看到：

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

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

来源

2014-11-25 08:57:23 ggael

谢谢 - 从图书馆的开发人员那里听到真的很棒！我的理由是看Eigen是否可以在他们有用的形状时正确地优化'm0.transpose（）* m1.transpose（）* m2 * m3' - 因此我想把所有东西放在表达式空间中，直到最后一刻。是否由于模板的工作原因，我无法对GeneralProduct进行cholesky分解，难道只有没有人有足够的关注将其添加到Eigen中吗？还是有理由这么做是愚蠢的？ – c0g 2014-11-25 13:00:00

问题是第一次乘法返回Eigen::GeneralProduct而不是MatrixXd而auto正在拾取返回类型。你可以从Eigen::GeneralProduct隐含地创建MatrixXd，所以当你显式声明它的类型时它能正常工作。见http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

编辑：我不是做铸件的本征产品或性能方面的专家。我刚刚从错误信息中猜出了答案，并从在线文档中得到了确认。性能分析始终是您检查代码不同部分性能的最佳选择。

来源

2014-11-24 20:45:55

对于投射到“MatrixXd”有没有性能上的提升？显然，在这个受限制的例子中，它会很小，但在现实生活中呢？你会说什么解决方案在这里 - 使用ProductReturnType来做这件事似乎在执行'c.tranpose（）* c'步骤时会有些生气。 – c0g 2014-11-24 21:02:42

我不是在Eigen专家，但像这样的图书馆经常回运营代理对象，然后用隐式转换或构造迫使实际工作。（表达式模板是这方面的一个极端例子。）这样可以避免在许多情况下不必要地复制整个数据矩阵。不幸的是，auto很高兴创建一个代理类型的对象，通常库的用户永远不会显式声明。由于您需要最终计算出数字，因此从铸造到MatrixXd本身不会造成任何性能下降。（斯科特迈尔斯，在“有效的现代C++”，给出了使用auto与vector<bool>，其中operator[](size_t i)返回代理的相关示例。）

来源

2014-11-24 21:21:59 sfjac

不要使用auto与本征表达式。我碰上这个更“戏剧化”的问题之前，请参阅

eigen auto type deduction in general product

，并通过本征的创作者（盖尔）不使用auto与本征的表现之一建议。

从表达式转换为MatrixXd这样的特定类型应该非常快，除非您想懒惰评估（因为在执行投射时评估结果）。

来源

2014-11-24 23:09:46 vsoftco

我认为这与你的稍有不同。你从'id'返回一个对象，然后由'autoresult'引用 - 然后'Id'的对象消失，'autoresult'引用一些不会退出的东西。 'c.transpose（）'引用'c'，它仍然存在，所以我的'auto cTc'不应该引用任何不存在的东西 - 类似于ggael作为一种解决方案提供的使用'auto id = Id（Foo， 4）'。我想测试Eigen是否可以正确地注意到'c.transpose（）* c = b.transpose（）* a.transpose（）* a * b'简化了数学运算，所以理想情况下它将保留为表达式。 – c0g 2014-11-24 23:35:29

是的，你是对的，但是，如果我没有得到解释，我完全不知道为什么代码表现异常。由于Eigen有很多类型（基本上每个表达式都是不同的类型），加上间接引用，'auto'会让事情变得复杂，因为你并不真正意识到底层发生了什么，以及如何评估这些表达式。 – vsoftco 2014-11-24 23:43:39

征和C++ 11的类型推断失败的矩阵产品

回答

相关问题