为什么MATLAB在矩阵乘法中速度这么快？

Question

为什么MATLAB在矩阵乘法中速度这么快？

我正在使用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+会有点接近，而不是60x慢点。

所以我的问题是-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完全相同，但是有一个巨大的飞跃2048从…1024?

编辑2：最新的MATLAB和4096x4096结果

Answer 1

这是我使用MATLABR2011a+的结果并行计算工具箱在一台带有特斯拉C 2070的机器上：>> A = rand(1024); gA = gpuArray(A);% warm up by executing the operations a couple of times, and then:>> tic, C = A * A;  tocElapsed time is 0.075396 seconds.>> tic, gC = gA * gA; tocElapsed time is 0.008621 seconds.MATLAB在矩阵乘法中使用了高度优化的库，这就是为什么普通MATLAB矩阵乘法速度如此之快的原因。这个gpuArray版本使用岩浆.使用R2014a更新在一台带有特斯拉K20c的机器上timeit和gputimeit职能：>> A = rand(1024); gA = gpuArray(A);>> timeit(@()A*A)ans =     0.0324>> gputimeit(@()gA*gA)ans =     0.0022使用R2018b更新在一台WIN 64机上，它有16个物理内核和一个Tesla V 100：>> timeit(@()A*A)ans =     0.0229>> gputimeit(@()gA*gA)ans =    4.8019e-04

Answer 2

这类问题是反复出现的，应该比“Matlab使用高度优化的库”或“Matlab使用MKL”一次更清楚地回答堆栈溢出。历史：矩阵乘法(与矩阵向量、向量乘法和许多矩阵分解一起)是线性代数中最重要的问题。工程师们从早期就开始用计算机解决这些问题。我不是历史专家，但很显然，当时每个人都用简单的循环重写了他的Fortran版本。随后出现了一些标准化，识别了大多数线性代数问题需要解决的“核”(基本例程)。然后，这些基本操作在一个名为：基本线性代数子程序(BLAS)的规范中标准化。然后，工程师可以在他们的代码中调用这些标准的、经过良好测试的blas例程，从而使他们的工作更加容易。BLAS：BLAS从第一级(定义标量向量和向量操作的第一个版本)发展到第二级(向量矩阵运算)到第三级(矩阵运算)，并提供了越来越多的“核”，使越来越多的基本线性代数运算标准化。最初的Fortran 77实现仍然可以在Netlib网站.争取更好的业绩：因此，多年来(特别是在BLAS第1级和第2级发布之间：80年代初)，随着向量操作和缓存层次结构的出现，硬件发生了变化。这些进化使BLAS子例程的性能大大提高成为可能。然后，不同的供应商来实现越来越有效率的BLAS例程。我不知道所有的历史实现(我不是天生的，也不是那个时候的孩子)，但最著名的两个实现出现在21世纪初：Intel MKL和GotoBLAS。您的Matlab使用Intel MKL，这是一个非常好的，优化的BLAS，这解释了您看到的伟大性能。矩阵乘法的技术细节：那么，为什么Matlab(MKL)在dgemm(双精度一般矩阵-矩阵乘法)？简单地说：因为它使用了矢量化和良好的数据缓存。更复杂的术语：参见文章由乔纳森·摩尔提供。基本上，当您在所提供的C+代码中执行乘法时，您对缓存一点也不友好。由于我怀疑您创建了一个指向行数组的指针数组，所以您在内部循环中访问“matice 2”的第k列：matice2[m][k]都很慢。实际上，当你访问matice2[0][k]，您必须得到矩阵数组0的k元素。然后在下一次迭代中，您必须访问matice2[1][k]，它是另一个数组(数组1)的第k个元素.然后在下一次迭代中访问另一个数组，依此类推.因为整个矩阵matice2不能放在最高的缓存中(它是8*1024*1024)，程序必须从主内存中获取所需的元素，从而损失大量时间。如果您只是转换了矩阵，以便访问将位于连续的内存地址中，那么您的代码将运行得更快，因为现在编译器可以同时在缓存中加载整个行。只需尝试这个修改后的版本：timer.start();float&nbsp;temp&nbsp;=&nbsp;0;//transpose&nbsp;matice2for&nbsp;(int&nbsp;p&nbsp;=&nbsp;0;&nbsp;p&nbsp;<&nbsp;rozmer;&nbsp;p++){ &nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(int&nbsp;q&nbsp;=&nbsp;0;&nbsp;q&nbsp;<&nbsp;rozmer;&nbsp;q++) &nbsp;&nbsp;&nbsp;&nbsp;{ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;tempmat[p][q]&nbsp;=&nbsp;matice2[q][p]; &nbsp;&nbsp;&nbsp;&nbsp;}}for(int&nbsp;j&nbsp;=&nbsp;0;&nbsp;j&nbsp;<&nbsp;rozmer;&nbsp;j++){ &nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(int&nbsp;k&nbsp;=&nbsp;0;&nbsp;k&nbsp;<&nbsp;rozmer;&nbsp;k++) &nbsp;&nbsp;&nbsp;&nbsp;{ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;temp&nbsp;=&nbsp;0; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(int&nbsp;m&nbsp;=&nbsp;0;&nbsp;m&nbsp;<&nbsp;rozmer;&nbsp;m++) &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;{ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;temp&nbsp;=&nbsp;temp&nbsp;+&nbsp;matice1[j][m]&nbsp;*&nbsp;tempmat[k][m]; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;} &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;matice3[j][k]&nbsp;=&nbsp;temp; &nbsp;&nbsp;&nbsp;&nbsp;}}timer.stop();因此，您可以看到缓存局部性如何极大地提高了代码的性能。现在是真实的dgemm实现将其利用到了一个非常广泛的层次：它们在由TLB(TransferingLookAbout缓冲器，长话短说：可以有效缓存的内容)定义的矩阵块上执行乘法，从而将处理的数据量准确地流到处理器。另一方面是矢量化，它们使用处理器的向量化指令来优化指令吞吐量，这在跨平台C+代码中是无法做到的。最后，人们声称这是因为Strassen‘s或CoppersSmith-Winograd算法是错误的，这两种算法在实践中都是不可实现的，因为上面提到的硬件考虑因素。

Answer 3

这就是为什么..MATLAB不会像在C+代码中那样循环每个元素来执行简单的矩阵乘法。当然，我假设你刚刚用了C=A*B而不是自己写乘法函数。