Eigen vs Matlab：并行矩阵乘法答案

【问题标题】：Eigen vs Matlab: parallelized Matrix-MultiplicationEigen vs Matlab：并行矩阵乘法
【发布时间】：2015-02-02 19:12:32
【问题描述】：

我想比较 Matlab 在矩阵乘法中的速度与 Eigen 3 在 Intel(R) Core(TM) i7-4770 CPU 上的速度@ 3.40GHz。包含Eigen的代码：

#include <iostream>
#include "Eigen/Dense"
#include <chrono>
#include <omp.h>


using namespace std;
using namespace Eigen;

const int dim=100;

int main()
{
    std::chrono::time_point<std::chrono::system_clock> start, end;

    int n;
    n = Eigen::nbThreads();
    cout<<n<<"\n";

    Matrix<double, Dynamic, Dynamic> m1(dim,dim);
    Matrix<double, Dynamic, Dynamic> m2(dim,dim);
    Matrix<double, Dynamic, Dynamic> m_res(dim,dim);

    start = std::chrono::system_clock::now();

    for (int i = 0 ; i <100000; ++i) {
        m1.setRandom(dim,dim);
        m2.setRandom(dim,dim);
        m_res=m1*m2;

    }

    end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;

    std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";

    return 0;
}

它使用g++ -O3 -std=c++11 -fopenmp 编译并使用OMP_NUM_THREADS=8 ./prog 执行。在我正在使用的 Matlab 中

function mat_test(N,dim)
%
% N:    how many tests
% dim:  dimension of the matrices

tic
parfor i=1:N
     A = rand(dim);
     B = rand(dim);
     C = A*B;
end
toc

结果是：Matlab 为 9 秒，Eigen 为 36 秒。在 Eigen 案例中我做错了什么？我可以排除矩阵的动态分配。此外，只使用了 3 个线程而不是 8 个。

编辑：

也许我说得不够清楚：任务是将dim=100 的100000 倍双值矩阵相乘每次随机填充，而不仅仅是一次。使用 Eigen 尽可能快地完成它。如果 Eigen 无法应对 Matlab，您会建议什么选择？

【问题讨论】：

您实际上想比较什么？除了矩阵乘法之外，您还在测量随机数生成的时间。 Eigen 可能使用不同且速度慢得多的随机数生成器。在开始计时之前，您也没有在 Matlab 中预先分配您的数组。我建议对每个部分进行分析或计时。您是否首先在不使用 OMP/parfor 的情况下比较了两者？最后，为什么你认为你做错了什么，只是因为 Eigen 比较慢？
比较 Eigen 中固定大小的矩阵。如果matlab使用惰性计算，则在使用C之前不必进行A*B，并且永远不会使用C，因此它可以消除......好吧，乘法。您可以通过 auto res = std::async( std::launch::deferred, [&]{return m1*m2;} ); 在 C++ 中进行模拟。那么，在共享指针中使用不可变矩阵并使用惰性求值（以及惰性矩阵）在它们上实现*？简而言之，你必须做一些事情来进行有意义的比较。
@Yakk：使用固定大小的矩阵不会改变结果。
@ggael：在循环外调用 setRandom 与 Matlab 版本不同。因此，这里都是关于随机生成器的，对吧？
@Yakk：实际上不在这个例子中。在循环中添加vec(i)=trace(A*B) 和下面的display(sum(vec)) 后，它归结为~10s。所以它必须生成矩阵A 和B 并计算它们的乘积。使用同一台计算机。

标签： c++ matlab c++11 eigen

【解决方案1】：

以下是合理使用 Eigen 的代码的更好版本。总结一下：

将setRandom() 移到基准测试循环之外。 setRandom()调用系统的rand()函数，比较慢。
使用.noalias() 来避免创建临时对象（仅当右侧是产品时才有意义）
将 OMP_NUM_THREADS 设置为真实的内核数，而不是超线程数。（在你的情况下是 4 个）
您的 CPU 支持仅由 Eigen 的 devel 分支支持的 AVX 和 FMA（将变为 3.3），因此请使用 devel 分支并使用 -mavx 和 -mfma 编译器选项启用它们（大约 x3.5 速度仅与上证所相比）

代码：

#include <iostream>
#include "Eigen/Dense"
#include <chrono>

using namespace std;
using namespace Eigen;

const int dim=100;

int main()
{
    std::chrono::time_point<std::chrono::system_clock> start, end;

    int n;
    n = Eigen::nbThreads();
    cout << n << "\n";

    Matrix<double, Dynamic, Dynamic> m1(dim,dim);
    Matrix<double, Dynamic, Dynamic> m2(dim,dim);
    Matrix<double, Dynamic, Dynamic> m_res(dim,dim);

    start = std::chrono::system_clock::now();

    m1.setRandom();
    m2.setRandom();

    for (int i = 0 ; i <100000; ++i) {
      m_res.noalias() = m1 * m2;
    }

    end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;

    std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";

    return 0;
}

【讨论】：

完全定时 setRandom 意味着代码不仅仅是定时乘法。
现在真的有办法以类似的速度在 Matlab 代码中做同样的事情吗？我看到这不是矩阵乘法，但是在 C++ 代码中只执行一次随机化不是在 Matlab 中完成的。换句话说：没有可以与matlabs代码竞争的随机生成器吗？
请澄清您的问题。您是否正在寻找非常快速的随机数生成器或快速矩阵乘法？如果您正在寻找快速随机数，那么还要说明您的要求，因为通常需要在生成器的质量和速度之间进行权衡。
@Goz, setRandom() 只被调用一次，因此与 10000 个矩阵产品相比，它的成本完全可以忽略不计。

【解决方案2】：

除了在循环外移动随机化（在 Eigen 和 Matlab 中），正如 ggael 建议的那样，在 Matlab 中将 parfor 替换为 for，因为在 Eigen 代码中，您按顺序处理矩阵。

我不确定 Matlab 是如何并行化其代码的：也许多个线程在同一对矩阵上工作，并在完成后切换到下一个；也许每个线程都处理自己的一对矩阵。有人可能会争辩说，后者可能会更快，因为更好地使用了特定于核心的缓存。

【讨论】：