本征 C++；就地矩阵乘法

Question

使用Eigen C++ Matrix库，如何高效将 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 中有没有办法做到这一点高效并从其优化中获利？

Answer 1

告诉 Eigen 您希望产品发生的最明确方式"in place"可能是：

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

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

在我的 Eigen 副本 (3.2.7) 上，上面的代码直接在 Transpose 表达式上调用了 MatrixBase::applyThisOnTheRight，这又在内部不幸地缩减为 B=A*B；与 applyOnTheLeft 的情况相同，所以在这种情况下你运气不好。

如果你真的需要避免 any nxm 临时，你可以手动向量方式执行乘积，比如：

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

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

Answer 2

您的 B.noalias() = A * B; 示例与 "store the result in A" 不匹配。 that 你只需要A *= B;。如果你 do 打算覆盖 B，那么你说 .noalias()

Answer 3

是的，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 C++; In-place matrix multiplication