特征 C++；就地矩阵乘法答案

【问题标题】：Eigen C++; In-place matrix multiplication特征 C++；就地矩阵乘法
【发布时间】：2018-01-17 08:14:11
【问题描述】：

使用 Eigen C++ 矩阵库，如何高效将 n x n 矩阵 A 与 n x m 矩阵 B 相乘并将结果存储在 A 中？也就是说，我该如何避免生成一个临时 n x m 矩阵并存储结果直接在B？

对于 m 非常大的应用程序（例如 100000）比 n （例如 3），这绝对是有道理的，因为它避免超大数组的应用。

以下代码，我无法工作：

     B.noalias() = A * B;

我认为，内部必须发生的事情是以下。 B 的每一列都应该分开处理。正在考虑的列column_i 必须复制到备份专栏column_tmp。那么，

     B(row_i, column_i) = A.row(row_i) * column_tmp; // dot-product

对于所有column_i = 0 到 m。 Eigen 中有没有办法做到这一点高效并从其优化中获利？

【问题讨论】：

【解决方案1】：

告诉 Eigen 您希望产品“就地”发生的最明确方式可能是：

B.transpose() *= A.transpose();
// or B.applyOnTheLeft(A);

也就是说，不能保证这不会在任何暂时发生，您必须为此相信内部的 Eigen 成本逻辑（Eigen 设计师可能更了解 :) ... 或检查在适当的分析表明这是一个真正的问题，而不仅仅是过早的优化之后，你可以通过调试器自行解决）。

在我的 Eigen 副本 (3.2.7) 上，上面直接在 Transpose 表达式上调用 MatrixBase::applyThisOnTheRight，而内部又可悲地减少为 B=A*B； applyOnTheLeft 的情况与此相同，所以在这种情况下你很不走运。

如果您确实需要暂时避免任何 nxm，您可以手动执行矢量乘积，例如：

for(auto i=0;i<B.cols();++i)
    B.col(i) = A * B.col(i);

假设B.rows()B.cols()，这将消耗更少的额外内存，但是你可能会错过一些重要的优化；确实，我想在这里有一个临时的仍然可以提供最好的权衡。

【讨论】：

在 '*=' 赋值中调用 '.transpose()' 的目的是什么。
@Frank-ReneSchäfer 使其等同于B = A*B（取两边的转置）
对不起，我不明白。 B *= A和B.transpose() *= A.transpose()这两个操作有什么区别？
@Frank-ReneSchäfer B*=A 是B=B*A 而我们想要B=A*B，即B^T=(A*B)^T，即B^T=B^T*A^T，即B^T*=A^T

【解决方案2】：

是的，使用 3x3 乘以 3xHuge 您的逐列评估确实有意义，但在一般情况下并非如此。例如，如果 n=m=1000，那么逐列求值将比当前的 Eigen 逻辑慢几个数量级。

如果你写：

B.noalias() = A * B;

Eigen 将遵循 col-by-col 评估（因为 A 在编译时很小），但结果会错误，因为 B 做了别名，基本上它会生成：

for j = 0..m-1
  B.col(j).noalias() = A * B.col(j);

为了优雅地解决这个问题，我们需要一种方法来说明只有不同的列不别名...建议：

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

确实是让 Eigen 在编译时知道这种信息的一个选项，尽管必须转置双方有点麻烦。并且正确的评估逻辑仍然需要在 Eigen 方面实现。目前此信息未被利用。

【讨论】：

另一种更明确的方式（当前 Eigen 无法识别为有效表达式）可能是 B.tranpose().rowwise() *= A.transpose()，这将表明调用者的意图......（以更糟糕的优化为代价）一般情况下）
我希望在未来的 Eigen 版本中看到的另一件事是通过表达式提供“优化提示”，例如 B.transpose() *= smallish(A).transpose(); ...
@MassimilianoJanes 应该可以将 applyThisOnTheLeft 和 applyThisOnTheRight 专门用于编译时的小型矩阵，以便逐列或逐行评估（它需要有人不过，是时候实现这一点了）

【解决方案3】：

您的 B.noalias() = A * B; 示例与“将结果存储在 A”中不匹配。那个你只需要A *= B;。如果您确实打算覆盖B，那么您在撒谎 .noalias()

【讨论】：