Eigen + MKL矩阵乘法比Matlab慢

Question

我正在一个C++程序中进行许多矩阵乘法，并且我使用Eigen (3.3.5)与Intel的MKL (2018.3.222)链接。我使用MKL的顺序版本，禁用OpenMP。问题是它比Matlab慢。

一些示例代码：

#define NDEBUG
#define EIGEN_USE_MKL_ALL

#include <iostream>
#include <chrono>
#include <Core>

using namespace Eigen;
using namespace std;

int main(){
    MatrixXd jac = 100*MatrixXd::Random(10*1228, 2850);
    MatrixXd res = MatrixXd::Zero(2850, 2850);

    for (int i=0; i<10; i++){
        auto begin = chrono::high_resolution_clock::now();
        res.noalias() = jac.transpose()*jac;
        auto end = chrono::high_resolution_clock::now();

        cout<<"time: "<<chrono::duration_cast<chrono::milliseconds>(end-begin).count() <<endl;
    }

    return 0;
}

平均报告时间约为8秒。使用Ubuntu 16.04上的g++ 6.4编译时，采用-O3和无调试符号。

Matlab代码：

m=100*(-1+2*rand(10*1228, 2850));
res = zeros(2850, 2850);
tic; res=m'*m; toc

它报告大约4秒，比原来快了两倍。我在同一系统上使用的是maxNumCompThreads(1)的Matlab R2017a版本。Matlab使用MKL 11.3。

如果没有MKL，只使用Eigen，则需要大约18秒。有什么方法可以将C++运行时间降到与Matlab相同的值？

谢谢。

后续编辑： 正如@Qubit所建议的那样，Matlab认识到我正在尝试将矩阵乘以它的转置，并进行了一些“隐藏”的优化。当我在Matlab中相乘两个不同的矩阵时，时间增加到了8秒。

所以，现在问题变成了：如何告诉Eigen这个矩阵乘积是“特殊的”并且可以进一步优化？

后续编辑2： 我尝试像这样做：

MatrixXd jac = 100*MatrixXd::Random(10*1228, 2850);
MatrixXd res = MatrixXd::Zero(2850, 2850);

auto begin = chrono::high_resolution_clock::now();
res.selfadjointView<Lower>().rankUpdate(jac.transpose(), 1);
res.triangularView<Upper>() = res.transpose();
auto end = chrono::high_resolution_clock::now();

MatrixXd oldSchool = jac.transpose()*jac;
if (oldSchool.isApprox(res)){
    cout<<"same result!"<<endl;
}
cout<<"time: "<<chrono::duration_cast<chrono::milliseconds>(end-begin).count() <<endl;

但现在它需要 9.4 秒（这是没有 MKL 的 Eigen 进行经典乘积所需时间的一半）。禁用 MKL 不会对此计时产生时间影响，因此我认为“rankUpdate”方法不使用 MKL ?!?

最后编辑： 我在 eigen 头文件中发现了一个 bug：

Core/products/GeneralMatrixMatrixTriangular_BLAS.h

在第55行有一个括号放错了位置。我进行了更改：

if ( lhs==rhs && ((UpLo&(Lower|Upper)==UpLo)) ) { \

转换为：

if ( lhs==rhs && ((UpLo&(Lower|Upper))==UpLo) ) { \

现在，我的C++版本和Matlab具有相同的执行速度（在我的系统上约为4秒）。

- Costin Florian Ciusdel

---

MATrix LABoratory是一款专门针对矩阵进行高度优化的软件。在处理矩阵方面，您不会找到比MATLAB更快的工具。 - Ander Biguri

我认为首先有几件事情需要检查。Matlab是否真的只使用了一个线程？如果你自己使用LAPACK会发生什么（确保用单个调用分配所有矩阵，这样它们就都是连续的）？同时，用一个一般的矩阵乘法运行相同的测试，因为你的矩阵乘以自身的转置是比较特殊的情况，Matlab可能已经意识到了，而C++可能没有。 - Qubit

我真的不明白这怎么会是 https://dev59.com/aG025IYBdhLWcg3wW0oc 的重复。那个问题比较了简单的3层循环实现和 MKL，而这里 C++ 和 Matlab 的方法都使用了相同的库。 - Qubit

@Wolfie 最高票答案直接指向了MKL作为原因，这也是OP在他的C++实现中使用的。此外，新问题似乎足够具体。他不想重写Matlab，只是希望在这个特殊例子中改进他的C++/Eigen实现，如果可能的话。类似于“你能说服Eigen仅计算矩阵乘法的上三角部分吗？”，因为结果矩阵将是对称的。但我确实同意这可以作为一个新问题提出。 - Qubit

Costin：请将您的解决方案发布为答案（当然，除了向Eigen项目报告错误之外）。 - Cris Luengo

Answer 1

既然您已经找出问题了，这里只是一些评论：

1. Core/products/GeneralMatrixMatrixTriangular_BLAS.h的问题已在devel分支中得到解决，但事实证明它从未被向3.3分支回溯。

2. 该问题现在在3.3分支中已得到解决。修复将成为3.3.6的一部分。

3. 在单线程模式下，内置Eigen和MKL之间的加速系数x2是没有意义的。请确保通过编译时使用-march=native来启用CPU支持的所有功能，除了-O3 -DNDEBUG之外。在我的Haswell 2.6GHz上，我得到3.4秒与3秒之间的时间。

- ggael