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

Question

我可以访问许多矩阵库，但是对于这个项目，我使用 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，但我想知道是否存在使用最少周期的高效现有实现。

Answer 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% 正确。