c++ - 对于矩阵乘法，Eigen + MKL 比 Matlab 慢

Question

我在 C++ 程序中进行了大量矩阵乘法运算，我使用与英特尔的 MKL (2018.3.222) 链接的 Eigen (3.3.5)。我使用 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 秒。在带有 g++ 6.4 的 Ubuntu 16.04 上使用 -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 头文件中发现了一个错误:

Core/products/GeneralMatrixMatrixTriangular_BLAS.h

在第 55 行。有一个错位的括号。我改变了这个:

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

为此:

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

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

Answer 1

真正的答案，因为你已经弄清楚了问题，但有一些评论:

1. 问题 Core/products/GeneralMatrixMatrixTriangular_BLAS.h 已经在 devel 分支中修复，但事实证明它从未被移植到 3.3 分支。

   <

2. 现在的问题是fixed在 3.3 分支中。该修复将成为 3.3.6 的一部分。

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