为什么MATLAB在矩阵乘法方面如此之快？

Question

我用CUDA、C++、C语言和java做了一些基准，并用MATLAB进行验证和矩阵生成。但是当我用MATLAB进行乘法运算时，2048x2048甚至更大的矩阵几乎都是立即乘法的。

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

< >只有CUDA是竞争性的，但我认为至少C++会有点接近，而不是{{CD2}}较慢。

所以我的问题是-MATLAB是怎么做到这么快的？

C++代码：

float temp = 0;
timer.start();
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

编辑： 我也不知道该怎么看待C结果。该算法与C++和java完全一样，但从^ {< CD4> }有一个巨大的跳转^ {CD3}}？

编辑2： 更新了MATLAB和4096x4096结果

Answer 1

这种问题反复出现，应该比“MATLAB使用高度优化的库”或“MATLAB使用MKL”一次堆栈溢出更清楚地回答。

历史记录：

矩阵乘法（连同矩阵向量、向量乘法和许多矩阵分解）是线性代数中最重要的问题。工程师们从早期就开始用计算机解决这些问题。

我不是历史方面的专家，但显然在那时，每个人都只是用简单的循环重写他的FORTRAN版本。随后出现了一些标准化，并确定了大多数线性代数问题需要解决的“核心”（基本例程）。这些基本操作随后在一个称为：基本线性代数子程序（BLAS）的规范中标准化。工程师可以在他们的代码中调用这些标准的、经过良好测试的BLAS例程，使他们的工作更加容易。

布拉斯：

BLAS从1级（定义标量向量和向量向量运算的第一个版本）发展到2级（向量矩阵运算）到3级（矩阵矩阵运算），并提供了越来越多的“核”，使越来越多的基本线性代数运算标准化。最初的FORTRAN 77实现仍然可以在Netlib's website上使用。

提高性能：

因此，多年来（特别是在BLAS 1级和2级版本之间：80年代早期），随着矢量操作和缓存层次结构的出现，硬件发生了变化。这些演进使得大幅度提高BLAS子例程的性能成为可能。随后，不同的供应商推出了他们的BLAS例程的实现，这些例程的效率越来越高。

我不知道所有的历史实现（那时我还没有出生或是个孩子），但最著名的两个实现出现在21世纪初：英特尔MKL和GotoBLAS。您的Matlab使用Intel MKL，这是一个非常好的、优化的BLAS，这解释了您看到的出色性能。

矩阵乘法的技术细节：

那么，为什么Matlab（MKL）在dgemm（双精度通用矩阵乘法）上速度如此之快呢？简单来说：因为它使用矢量化和良好的数据缓存。用更复杂的术语来说：请参阅Jonathan Moore提供的article。

基本上，当你在你提供的C++代码中执行乘法时，你根本就不喜欢缓存。因为我怀疑您创建了一个指向行数组的指针数组，所以您在内部循环中对“matice2”：matice2[m][k]的第k列的访问非常慢。实际上，当您访问matice2[0][k]时，必须获得矩阵数组0的第k个元素。然后在下一次迭代中，必须访问matice2[1][k]，这是另一个数组（数组1）的第k个元素。然后在下一次迭代中访问另一个数组，依此类推。。。由于整个矩阵matice2不能放在最高的缓存中（它的8*1024*1024字节大），程序必须从主内存中获取所需的元素，从而浪费大量时间。

如果您只是转置矩阵，以便访问位于连续的内存地址中，那么您的代码将运行得更快，因为现在编译器可以同时将整行加载到缓存中。请尝试此修改版本：

timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
    for (int q = 0; q < rozmer; q++)
    {
        tempmat[p][q] = matice2[q][p];
    }
}
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * tempmat[k][m];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

所以您可以看到缓存局部性是如何大大提高代码的性能的。现在，真正的dgemm实现将这一点扩展到了一个非常广泛的层次：它们在由TLB（Translation lookaside buffer，long story short:what can effective be cached）大小定义的矩阵块上执行乘法，以便它们将它能够处理的数据量准确地流到处理器。另一个方面是矢量化，它们使用处理器的矢量化指令来获得最佳的指令UpPoT，你不能从你的跨平台C++代码中真正做到。

最后，人们声称这是因为斯特拉森或铜匠-温诺格拉德算法是错误的，这两种算法在实践中都无法实现，因为硬件考虑上述。

Answer 2

This is why。Matlab不会像你在C++代码中那样对每个元素进行循环，从而执行一个天真的矩阵乘法。

当然，我假设您只是使用C=A*B，而不是自己编写乘法函数。

Answer 3

下面是我在一台装有特斯拉C2070的机器上使用MATLAB R2011a+Parallel Computing Toolbox的结果：

>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.

MATLAB使用高度优化的库进行矩阵乘法，这就是为什么简单的MATLAB矩阵乘法如此快速的原因。gpuArray版本使用MAGMA。

在具有Tesla K20c和新的timeit和gputimeit功能的计算机上使用R2014a进行更新：

>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
    0.0324
>> gputimeit(@()gA*gA)
ans =
    0.0022

在具有16个物理核心和Tesla V100的WIN64计算机上使用R2018b进行更新：

>> timeit(@()A*A)
ans =
    0.0229
>> gputimeit(@()gA*gA)
ans =
   4.8019e-04