矩阵积的 Cholesky 的特征和 C++11 类型推理失败

Eigen and C++11 type inference fails for Cholesky of matrix product

本文关键字：类型 C++11 推理 失败 特征 Cholesky      更新时间：2023-10-16

我正在尝试使用特征和 C++11"auto" 类型对矩阵乘积进行乔莱斯基分解及其转置。当我尝试做时，问题就来了 

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是否会导致问题。

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

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

让我总结一下正在发生的事情以及为什么它是错误的。首先，让我们用它们所采用的类型实例化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);

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

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，那么您仍然必须自己做。

我不是Eigen的专家，但像这样的库经常从操作中返回代理对象，然后使用隐式转换或构造函数来强制实际工作。（表达式模板就是一个极端的例子。在许多情况下，这避免了对完整数据矩阵的不必要复制。不幸的是，auto很乐意只创建一个代理类型的对象，通常库的用户永远不会显式声明。由于您最终需要计算数字，因此从铸造到MatrixXd本身不会对性能产生影响。（Scott Meyers在"Effective Modern C++"中给出了将auto与vector<bool>一起使用的相关示例，其中operator[](size_t i)返回一个代理。

不要

将auto与特征表达式一起使用。我之前遇到了更"戏剧性"的问题，请参阅

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

并被其中一位特征创建者（Gael）建议不要在特征表达式中使用auto。

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

相关文章：