为什么MATLAB在matrix乘法中如此之快？

我正在用CUDA，C ++，C＃和Java做一些基准testing，并使用MATLAB进行validation和matrix生成。 但是，当我乘以MATLAB，2048×2048，甚至更大的matrix几乎立即倍增。

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是有竞争力的，但我认为至lessC ++会稍微接近一点，而不是慢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＃的结果。 这个algorithm和C ++和Java一样，但是从1024开始有2048的巨大跳跃？

编辑2：更新MATLAB和4096×4096结果

Answer 1

下面是使用MATLAB R2011a + 并行计算工具箱在带有Tesla C2070的机器上的结果：

 >> 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使用高度优化的库进行matrix乘法，这就是为什么普通的MATLABmatrix乘法速度如此之快的原因。 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 

Answer 2

这种问题是反复出现的，应该在Stackoverflow上比“Matlab使用高度优化的库”或“Matlab使用MKL”更清楚地回答。

历史：

matrix乘法（与matrix向量，向量向量乘法和许多matrix分解一起）是线性arrays中最重要的问题。 自从早期以来，工程师一直在用电脑解决这些问题。

我不是历史上的专家，但显然当时，每个人都用简单的循环重写了他的Fortran版本。 随后出现了一些标准化问题，为了解决大多数线性代数问题所需要的“内核”（基本例程）。 然后这些基本操作被标准化为称为基本线性代数子程序（BLAS）的规范。 然后工程师可以在他们的代码中调用这些标准，经过充分testing的BLAS例程，使他们的工作更容易。

BLAS：

BLAS从级别1（定义标量vector和vectorvector运算的第一个版本）发展到级别2（vectormatrix运算）到级别3（matrix运算），并提供了越来越多的“内核”和更多的基本线性代数运算。 最初的Fortran 77实现仍然可以在Netlib的网站上find 。

迈向更好的performance：

所以多年来（特别是在BLAS 1级和2级版本之间：80年代初），随着vector操作和caching层次的出现，硬件发生了变化。 这些演变使得可以大大提高BLAS子程序的性能。 不同的供应商随之带来了更高效的BLAS例程的实现。

我不知道所有的历史实现（当时我还没有出生或者是一个小孩），但是在二十一世纪初出现了两个最着名的：英特尔MKL和GotoBLAS。 您的Matlab使用英特尔MKL，这是一个非常好的，优化的BLAS，并解释了您看到的卓越性能。

matrix乘法的技术细节：

那么为什么Matlab（MKL）在dgemm （双精度一般matrix – matrix乘法）上这么快？ 简而言之：因为它使用vector化和良好的数据caching。 更复杂的说法是：参见Jonathan Moore提供的文章 。

基本上，当你在你提供的C ++代码中进行乘法运算时，你并不总是caching友好的。 因为我怀疑你创build了一个指向行数组的指针数组，所以你在内层循环访问到“matice2”的第k列： matice2[m][k]是非常缓慢的。 事实上，当你访问matice2[0][k] ，你必须得到matrix数组0的第k个元素。 然后在下一次迭代中，您必须访问另一个数组（数组1）的第k个元素matice2[1][k] 。 然后在下一次迭代中访问另一个数组，等等…由于整个matrixmatice2不能适应最高的caching（它的大小是8*1024*1024个字节），所以程序必须从main记忆，失去了很多时间。

如果你只是转置了matrix，所以访问将是在连续的内存地址，你的代码已经运行得更快了，因为现在编译器可以同时加载caching中的整个行。 试试这个修改版本：

 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(); 

所以你可以看到caching本地如何提高了你的代码的性能。 现在，真正的dgemm实现将其用于非常广泛的层面：它们在由TLB（翻译旁视缓冲区，长话短说：可以有效地被caching的）的大小定义的matrix块上执行乘法，以便它们stream向处理器正是它可以处理的数据量。 另一个方面是vector化，他们使用处理器的向量化指令来获得最佳的指令吞吐量，而这在跨平台的C ++代码中是无法做到的。

最后，声称这是因为Strassen或Coppersmith-Winogradalgorithm是错误的，由于上面提到的硬件考虑，这两种algorithm在实践中都是不可实现的。

Answer 3

这是为什么 。 MATLAB不会像你在C ++代码中那样循环遍历每一个元素来执行一个简单的matrix乘法。

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

Answer 4

Matlab在一段时间内引入了LAPACK，所以我假设它们的matrix乘法至less使用了这样的速度。 LAPACK源代码和文档是随时可用的。

你也可以看看Goto和Van De Geijn的论文“高性能matrix乘法剖析” http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

Answer 5

在进行matrix乘法时，使用O(n^3)时间的朴素乘法。

存在matrix乘法algorithm，其中O(n^2.4) 。 这意味着在n=2000您的algorithm需要的计算量是最佳algorithm的100倍。
你应该真的检查维基百科页面的matrix乘法的进一步信息的有效方法来实现它。

Answer 6

您需要小心使用C ++进行公平的比较。 你可以发布C ++代码，显示你用于matrix乘法的核心内部循环吗？ 大多数情况下，我关心你的内存布局，以及你是否在做浪费的事情。

我已经编写了与Matlab一样快的C ++matrix乘法，但是需要注意。 （编辑：在Matlab之前使用GPU的。）

你几乎可以保证Matlab在这些“内置”function上浪费很less的周期。 我的问题是，你在哪里浪费周期？ （没有冒犯的意思）

Answer 7

答案是LAPACK和BLAS库使得MATLAB在matrix运算上非常快速，而不是MATLAB中的任何专有代码。

使用C ++代码中的LAPACK和/或BLAS库进行matrix运算，您应该获得与MATLAB相似的性能。 这些图书馆应该可以在任何现代系统上自由使用，并且在学术界已经发展了几十年。 请注意，有多种实现方式，包括一些封闭的源代码，如“ 英特尔MKL” 。

为什么BLAS / LAPACK这么快？ （一）有效的algorithm和（二）微调，利用CPU架构。 例如。 事实certificatematrix乘法可以用O（n ^ 2.807）代替O（n ^ 3）algorithm来完成，并且这被合并到几个BLAS实现中。 巧妙的操作分组最大限度地减less了将数字移入或移出CPU寄存器等的需求。

---

顺便说一句，根据我的经验，从c直接调用LAPACK库是一个严重的痛苦（但值得）。 你需要非常精确地阅读文档。

Answer 8

根据你的Matlab版本，我相信它可能已经在使用你的GPU了。

另一件事; Matlab跟踪matrix的许多属性; 其对angular线，hermetian等等，并且基于它的algorithm专门化。 也许它是基于你传递的零matrix的专业，或类似的东西？ 也许是caching重复的函数调用，这会弄乱你的时间？ 也许它优化了重复使用的matrix产品？

为防止发生这种情况，请使用随机数字matrix，并确保您通过将结果打印到屏幕或磁盘或其他部件来强制执行。

Answer 9

使用双打和一个固定的数组，而不是三个单独的线索我的C＃代码几乎相同的结果作为C + + / Java（与您的代码：1024 – 有点快一点，2048年 – 约140和4096 – 约22分钟）

                 1024x1024 2048x2048 4096x4096
                 --------- --------- ---------
你的C ++（ms）6137.10 64369.29 551390.93
我的C＃（ms）9730.00 90875.00 1062156.00

这里是我的代码：

  const int rozmer = 1024; double[][] matice1 = new double[rozmer * 3][]; Random rnd = new Random(); public Form1() { InitializeComponent(); System.Threading.Thread thr = new System.Threading.Thread(new System.Threading.ThreadStart(() => { string res = ""; Stopwatch timer = new Stopwatch(); timer.Start(); double temp = 0; int r2 = rozmer * 2; for (int i = 0; i < rozmer*3; i++) { if (matice1[i] == null) { matice1[i] = new double[rozmer]; { for (int e = 0; e < rozmer; e++) { matice1[i][e] = rnd.NextDouble(); } } } } timer.Stop(); res += timer.ElapsedMilliseconds.ToString(); int j = 0; int k = 0; int m = 0; timer.Reset(); timer.Start(); for (j = 0; j < rozmer; j++) { for (k = 0; k < rozmer; k++) { temp = 0; for (m = 0; m < rozmer; m++) { temp = temp + matice1[j][m] * matice1[m + rozmer][k]; } matice1[j + r2][k] = temp; } } timer.Stop(); this.Invoke((Action)delegate { this.Text = res + " : " + timer.ElapsedMilliseconds.ToString(); }); })); thr.Start(); } 

Answer 10

你是否检查过所有的实现对algorithm使用了multithreading优化？ 他们是否使用相同的乘法algorithm？

我真的怀疑这一点。

Matlab本身并不是很快，你可能使用的是慢速实现。

有效matrix乘法的algorithm

Answer 11

“为什么matlab比其他程序更快做xxx”的一般答案是matlab有很多内置的优化function。

其他使用的程序通常不具备这些function，因此人们使用自己的创意解决scheme，这比专业优化的代码慢得多。

这可以用两种方式来解释：

1）通用/理论方式：Matlab没有明显更快，你只是做了基准testing错误

2）现实的方法：对于这个东西，Matlab在实践中速度更快，因为c ++语言太简单易用了。

Answer 12

这种鲜明的对比不仅是由于Matlab的惊人的优化（正如许多其他答案已经讨论过的），而且也是由于您将matrix制定为对象。

你似乎是matrix列表的列表？ 列表列表包含指向包含matrix元素的列表的指针。 包含的列表的位置是任意指定的。 当你循环你的第一个索引（行号？）时，内存访问的时间非常重要。 相比之下，为什么不尝试使用下面的方法实现matrix作为单个列表/vector？

 #include <vector> struct matrix { matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {} int n_row; int n_col; std::vector<double> M; double &operator()(int i, int j); }; 

和

 double &matrix::operator()(int i, int j) { return M[n_col * i + j]; } 

应该使用相同的乘法algorithm，以使翻转次数相同。 （对于n的matrix，n ^ 3）

我要求你把它计算出来，以便结果可以和你以前（在同一台机器上）的结果相比较。 通过比较，您将显示准确的内存访问时间！

Answer 13

MATLAB使用来自英特尔（被称为英特尔MKL）的LAPACK（特别是dgemm函数 ）的高度优化的实现。 速度该库利用了包括SIMD指令和多核处理器在内的处理器function。 他们不logging他们使用的具体algorithm。 如果您要从C ++调用“英特尔MKL”，则应该看到类似的性能。

我不确定什么库MATLAB用于GPU乘法，但可能像nVidia CUBLAS 。

Answer 14

C ++的速度很慢，因为你没有使用multithreading。 本质上，如果A = BC，它们都是matrix，那么A的第一行可以独立于第二行计算，如果A，B和C都是n乘nmatrix，则可以通过乘以一个因素n ^ 2，如

a_ {i，j} = sum_ {k} b_ {i，k} c_ {k，j}

如果使用Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ]，则内置multithreading，并且线程数可以调整。