为什么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 ++会稍微靠近，而不是慢60倍。
所以我的问题是 - 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相同，但是从1024开始有2048的巨大跳跃？
编辑2：更新了MATLAB和4096x4096结果
这是我在使用Tesla 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的机器上使用R2014a进行更新，以及新的timeit和gputimeit函数：
>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
0.0324
>> gputimeit(@()gA*gA)
ans =
0.0022
这种问题反复出现，应该在Stackoverflow上比“Matlab使用高度优化的库”或“Matlab使用MKL”更清楚地回答。
历史：
矩阵乘法（与矩阵向量，向量向量乘法和许多矩阵分解一起）是线性阵列中最重要的问题。 自早期以来，工程师一直在用计算机解决这些问题。
我不是历史专家，但显然当时，每个人都用简单的循环重写了他的Fortran版本。 随后出现了一些标准化问题，为了解决大多数线性代数问题需要识别的“内核”（基本例程）。 然后将这些基本操作标准化为称为基本线性代数子程序（BLAS）的规范。 然后工程师可以在代码中调用这些标准的，经过充分测试的BLAS例程，使他们的工作变得更容易。
BLAS：
BLAS从级别1（定义标量矢量和矢量矢量运算的第一个版本）发展到级别2（矢量矩阵运算）到级别3（矩阵运算），并提供越来越多的“内核”和更多的基本线性代数运算。 最初的Fortran 77实现仍然可以在Netlib的网站上找到。
迈向更好的表现：
所以多年来（特别是BLAS 1级和2级发布之间：80年代初），随着向量操作和缓存层次的出现，硬件发生了变化。 这些演变使得可以大幅提高BLAS子程序的性能。 不同的供应商随之带来了更高效的BLAS例程的实现。
我不知道所有的历史实现（当时我还没有出生或者是一个小孩），但是最出名的两个最早出现在21世纪初：英特尔MKL和GotoBLAS。 你的Matlab使用英特尔MKL，这是一个非常优秀的，经过优化的BLAS，并解释了你所看到的卓越性能。
矩阵乘法的技术细节：
那么为什么Matlab（MKL）在dgemm （双精度通用矩阵 - 矩阵乘法）上如此之快？ 简而言之：因为它使用矢量化和良好的数据缓存。 更复杂的说法是：请参阅Jonathan Moore提供的文章。
基本上，当您在您提供的C ++代码中执行乘法运算时，您完全不是缓存友好的。 因为我怀疑你创建了一个指向行数组的指针数组，所以你在内层循环中访问“matice2”的第k列： matice2[m][k]非常缓慢。 事实上，当你访问matice2[0][k] ，你必须得到矩阵的数组0的第k个元素。 然后在下一次迭代中，您必须访问另一个数组（数组1）的第k个元素matice2[1][k] 。 然后在下一次迭代中访问另一个数组，等等...由于整个矩阵matice2不能适应最高的缓存（它的大小为8*1024*1024个字节），程序必须从main记忆，失去了很多时间。
如果你只是转置了矩阵，所以访问会在连续的内存地址中，你的代码已经运行得更快了，因为编译器现在可以同时加载缓存中的整个行。 试试这个修改后的版本：
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（翻译后备缓冲区，长话短说：可以有效缓存的内容）大小定义的矩阵块上执行乘法，以便它们流向处理器正是它可以处理的数据量。 另一方面是矢量化，他们使用处理器的矢量化指令来获得最佳的指令吞吐量，而这是跨平台C ++代码无法实现的。
最后，声称这是因为Strassen或Coppersmith-Winograd算法是错误的，由于上面提到的硬件考虑，这两种算法在实践中都是不可实现的。
这就是为什么。 MATLAB不会像你在C ++代码中那样遍历每一个元素来执行一个简单的矩阵乘法。
当然，我假设你只是使用C=A*B而不是自己写一个乘法函数。