Eigen 中的高效矩阵转置矩阵乘法答案

【问题标题】：Efficient matrix transpose matrix multiplication in EigenEigen 中的高效矩阵转置矩阵乘法
【发布时间】：2017-03-10 07:07:08
【问题描述】：

我可以访问许多矩阵库，但对于这个项目，我使用的是 Eigen，因为它的编译时间定义和包含 SVD。

现在，我正在做以下操作：

Eigen::Matrix<double,M,N> A;     // populated in the code

Eigen::Matrix<double,N,N> B = A.transpose() * A;

据我所知，这会复制 A 并形成转置，然后再乘以 A。此操作在相对较小的矩阵（M=20-30,N=3）上执行，但每秒数百万次，这意味着它必须尽可能快。

我读到使用以下方法更快：

B.noalias() = A.transpose() * A;

我可以编写自己的子例程，接受 A 作为输入并填充 B，但我想知道是否存在使用最少周期的有效的现有实现。

【问题讨论】：

【解决方案1】：

首先，由于 Eigen 依赖于模板表达式，A.transpose() 不会计算为临时值。

其次，在：

Matrix<double,N,N> B = A.transpose() * A;

Eigen 知道B 不能出现在表达式的右侧（因为这里编译器调用了 B 的构造函数），因此根本不会创建临时变量。这相当于：

Matrix<double,N,N> B;             // declare first
B.noalias() = A.transpose() * A;  // eval later

最后，对于这么小的矩阵，我不希望使用 B.selfadjointView().rankUpdate(A) 会有所帮助（正如 kennytm 评论中所建议的那样）。

另一方面，当 N=3 时，可能值得尝试惰性实现：

B = A.transpose().lazyProduct(A)

只是为了确定。 Eigen 具有内置的启发式方法来选择最佳的产品实现，但由于启发式方法必须简单快速地进行评估，因此它可能不是 100% 正确的。

【讨论】：

谢谢。懒惰的项目技巧是非常如何。现在，我最终做了一些完全不同的事情，因为我在提出这个后发现 Eigen 在 GPU 上的 cuda 中不起作用。不过我喜欢图书馆。另外，根本不构建 A 是最有效的，这就是我所做的。