为什么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倍。我也不知道如何看待c#的结果。该算法与c++和Java一样，但从1024到2048有一个巨大的跳跃。

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();
60948 次浏览
小开

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

另一件事;Matlab可以跟踪矩阵的许多性质;无论是对角线的，还是赫尔密斯的，等等，并在此基础上专门设计算法。也许它的专门化是基于你传递给它的0矩阵，或者类似的东西?也许它正在缓存重复的函数调用，这会打乱您的计时?也许它优化了重复未使用的矩阵积?

为了防止这样的事情发生，使用一个随机数矩阵，并确保通过将结果打印到屏幕或磁盘或其他地方来强制执行。

小开
最佳答案

以下是我在一台特斯拉C2070的机器上使用MATLAB R2011a + 并行计算工具箱的结果:

>> 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版本使用岩浆。

在Tesla K20c的机器上使用R2014a进行更新，以及新的timeit和gputimeit函数:

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

使用R2018b更新在WIN64机器上，有16个物理核和特斯拉V100:

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

(注意:在某些时候(我忘记确切的时间)gpuArray从MAGMA切换到cuBLAS -岩浆仍然用于一些gpuArray操作)

使用R2022a进行更新在有32个物理核和A100 GPU的WIN64机器上:

>> timeit(@()A*A)
ans =
0.0076
>> gputimeit(@()gA*gA)
ans =
2.5344e-04
小开

Matlab在一段时间前集成了LAPACK，所以我假设他们的矩阵乘法至少用了这么快的速度。LAPACK源代码和文档是现成的。

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

小开

这就是为什么。MATLAB不像在c++代码中那样，通过遍历每一个元素来执行简单的矩阵乘法。

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

小开

对于“为什么matlab在做xxx时比其他程序快”的一般答案是，matlab有很多内建的优化函数。

使用的其他程序通常没有这些功能，因此人们应用自己的创造性解决方案，这比专业优化的代码慢得多。

这有两种解释:

1)常见的/理论的方法:Matlab并没有明显更快，你只是做错了基准测试

2)现实的方法:对于这些东西，Matlab在实践中更快，因为像c++这样的语言太容易以无效的方式使用。

小开

当做矩阵乘法时，你使用朴素乘法方法，它的时间为O(n^3)。

存在矩阵乘法算法，该算法取O(n^2.4)。这意味着在n=2000，你的算法需要的计算量是最佳算法的100倍。
你真的应该去维基百科的矩阵乘法页面上查看，以获得更多关于有效实现矩阵乘法的信息

小开

这种鲜明的对比不仅是由于Matlab的惊人优化(正如许多其他答案已经讨论过的那样)，而且是由于您将矩阵作为一个对象来表述的方式。

看起来你把矩阵变成了列表的列表?列表的列表包含指向列表的指针，然后包含您的矩阵元素。所包含列表的位置是任意分配的。当循环遍历第一个索引(行号?)时，内存访问时间非常重要。相比之下，为什么不尝试实现矩阵作为一个单一的列表/向量使用下面的方法?

#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];
}

应该使用相同的乘法算法，以使触发器的数量相同。(大小为n的方阵为n^3)

我要求您对它进行计时，以便结果与前面(在同一台机器上)的结果相比较。通过比较，您将准确地显示内存访问时间有多么重要!

小开

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

历史:

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

我不是历史方面的专家，但显然在那个时候，每个人都只是用简单的循环重写了他的FORTRAN版本。随后出现了一些标准化，确定了大多数线性代数问题需要解决的“内核”(基本例程)。这些基本操作后来在一个叫做:基本线性代数子程序(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++代码中执行乘法时，您根本不是缓存友好型的。由于我怀疑您创建了一个指向行数组的指针数组，因此您在内部循环中访问“matie2”: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，长话短说:可以有效缓存的东西)大小定义的矩阵块执行乘法运算，这样它们就可以精确地向处理器传输它可以处理的数据量。另一个方面是向量化，它们使用处理器的向量化指令来优化指令吞吐量，这在跨平台的c++代码中是无法实现的。

最后，有人声称这是因为Strassen's或Coppersmith-Winograd算法是错误的，这两个算法在实践中都是不可实现的，因为上面提到的硬件方面的考虑。

小开

MATLAB使用高度优化的LAPACK实现，从英特尔称为英特尔数学内核库(英特尔MKL) -特别是dgemm函数。这个库充分利用了处理器的特性，包括SIMD指令和多核处理器。他们没有记录他们使用的具体算法。如果从c++调用Intel MKL，应该会看到类似的性能。

我不确定MATLAB使用什么库来进行GPU乘法，但可能是类似英伟达CUBLAS的东西。

小开

它在c++中很慢，因为你没有使用多线程。本质上，如果A = B C，其中它们都是矩阵，则A的第一行可以独立于第二行计算，等等。如果A、B和C都是n × n矩阵，您可以将乘法运算速度提高一个因子n^2，如

A_ {i,j} = sum_{k} b_{i,k} c_{k,j}

如果你使用Eigen [http://eigen.tuxfamily.org/dox/GettingStarted.html]，多线程是内置的，线程的数量是可调的。

小开

答案是LAPACK和布拉斯特区库使MATLAB在矩阵运算方面速度惊人，而不是MATLAB的任何专有代码。

在你的c++代码中使用LAPACK和/或布拉斯特区库来进行矩阵运算，你应该得到与MATLAB类似的性能。这些库在任何现代系统上都应该是免费的，其中一部分是学术界在几十年里开发出来的。注意，有多个实现，包括一些封闭源代码，如英特尔MKL。

讨论BLAS如何获得高性能在这里可以找到。

顺便说一句，在我的经验中，直接从c调用LAPACK库是一种严重的痛苦(但值得)。你需要非常精确地阅读文档。

小开

因为MATLAB .最初是为数值线性代数(矩阵操作)开发的编程语言，它有专门为矩阵乘法开发的库。并且now MATLAB也可以使用图形处理器(图形处理单元)来实现这一点。

如果我们看一下计算结果:

             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

然后我们可以看到不仅MATLAB在矩阵乘法方面如此之快:CUDA C(来自NVIDIA的编程语言)有一些比MATLAB更好的结果。CUDA C也有专门为矩阵乘法开发的库，它使用gpu。

MATLAB简史

Cleve Moler是新墨西哥大学(University of New Mexico)计算机科学系的系主任，他在20世纪70年代末开始开发MATLAB。他设计它是为了让他的学生能够访问LINPACK(一个用于执行数值线性代数的软件库)和EISPACK(一个用于线性代数数值计算的软件库)，而不需要学习Fortran。它很快传播到其他大学，并在应用数学界找到了强大的受众。1983年，工程师杰克·利特尔(Jack Little)在莫勒访问斯坦福大学(Stanford University)时接触到了它。认识到它的商业潜力，他与莫勒和史蒂夫·班格特一起加入了公司。他们用C语言重写了MATLAB，并在1984年创建了MathWorks以继续其发展。这些重新编写的库被称为JACKPAC。2000年，MATLAB被重写，以使用一组新的矩阵操作库LAPACK(数值线性代数的标准软件库)。

Source .

什么是CUDA C

CUDA C还使用专门为矩阵乘法开发的库，如OpenGL(开放图形库)。它还使用GPU和Direct3D(在MS Windows上)。

CUDA平台被设计用于C、c++和Fortran等编程语言。这种可访问性使并行编程专家更容易使用GPU资源，相比之下，以前的api如Direct3D和OpenGL需要图形编程的高级技能。同时，CUDA支持诸如OpenACC和OpenCL之类的编程框架。

enter image description here

CUDA处理流程示例:

  1. 将数据从主存复制到GPU内存
  2. CPU启动GPU计算内核
  3. GPU的CUDA内核并行执行内核
  4. 将结果数据从GPU内存复制到主内存

比较CPU和GPU的执行速度

我们运行了一个基准测试，在Intel Xeon处理器X5650上，然后使用NVIDIA Tesla C2050 GPU，测量了在网格大小为64、128、512、1024和2048的情况下执行50个时间步骤所花费的时间。

enter image description here

对于2048的网格大小，该算法显示计算时间减少了7.5倍，从CPU上的超过1分钟到GPU上的不到10秒。对数刻度图显示，对于较小的网格大小，CPU实际上更快。然而，随着技术的发展和成熟，GPU解决方案越来越能够处理更小的问题，我们预计这一趋势将继续下去。

Source .

来自CUDA C编程指南的介绍:

由于市场对实时、高清3D图形的需求无法满足，可编程图形处理器单元或GPU已经发展成为高度并行、多线程、多核处理器，具有巨大的计算能力和非常高的内存带宽，如Figure 1和Figure 2所示。

图1所示。 CPU和GPU每秒浮点运算数

enter image description here

图2。内存CPU和GPU的带宽

enter image description here

CPU和GPU之间的浮点运算能力差异背后的原因是GPU专门用于计算密集型、高度并行的计算——这正是图形渲染的意义所在——因此被设计成更多的晶体管用于数据处理，而不是数据缓存和流控制，如Figure 3所示。

图3。GPU将更多的晶体管用于数据处理

enter image description here

更具体地说，GPU特别适合解决可以表示为数据并行计算的问题——相同的程序并行地在许多数据元素上执行——具有高算术强度——算术运算与内存运算的比率。由于对每个数据元素执行相同的程序，因此对复杂的流控制要求较低，并且由于它在许多数据元素上执行，具有较高的算术强度，因此可以用计算而不是大数据缓存来隐藏内存访问延迟。

数据并行处理将数据元素映射到并行处理线程。许多处理大型数据集的应用程序可以使用数据并行编程模型来加快计算速度。在3D渲染中，大量的像素和顶点被映射到并行线程。类似地，图像和媒体处理应用程序，如渲染图像的后处理、视频编码和解码、图像缩放、立体视觉和模式识别，可以将图像块和像素映射到并行处理线程。事实上，许多图像渲染和处理领域之外的算法都被数据并行处理加速了，从一般的信号处理或物理模拟到计算金融学或计算生物学。

Source .

先进的阅读

一些有趣的面孔

我写过c++矩阵乘法，它和Matlab一样快，但它需要一些小心。(在Matlab使用图形处理器之前)。

Сitation from 这个答案 .。


提交答案