为什么MATLAB在矩阵乘法中如此快?

206

我正在使用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();

14
可能是由于你使用的算法不同所导致的问题。 - Robert J.
26
确保Matlab没有缓存你的结果,它是一个棘手的东西。首先确保实际进行了计算,然后进行比较。 - rubenvb
29
LAPACK和向量化。http://www.mathworks.com/company/newsletters/news_notes/clevescorner/winter2000.cleve.htmlLAPACK是用于数值线性代数计算的软件库,它可以高效地求解线性方程组、特征值问题和奇异值分解等。向量化可以使这些计算更快速地执行,它指的是将代码重写为使用向量操作,以便充分利用现代处理器的并行能力。在Matlab中,向量化是一种常见的优化技巧,可以显著提高代码的执行速度。 - James
11
我认为这篇文章确实很有趣,但我想看到更适当的基准测试。例如,我知道Matlab R2011a会自动使用多线程,并且矩阵乘法是使用Intel的MKL / BLAS库实现的。因此,如果使用mkl调用来进行矩阵乘法,我猜c ++会更快。问题在于Matlab的开销是多少。我知道这取决于矩阵乘法的其他细节,但目前上述数字相当无意义。 - Lucas
1
你可以使用运行时间为O(n^2.81)的Strassen算法来进行大型矩阵乘法,这比原生乘法快10倍左右,原生乘法的运行时间为O(n^3)。此外,SSE/AVX可以帮助你实现8-20倍的代码执行速度提升。总体而言,你可以拥有一个比Matlab更快的C语言实现。 - DU Jiaen
显示剩余14条评论
12个回答

211

这种问题经常出现,需要在Stack Overflow上更清楚地回答,而不是像“MATLAB使用高度优化的库”或“MATLAB使用MKL”那样回答。

历史:

矩阵乘法(以及矩阵向量、向量向量乘法和许多矩阵分解)是线性代数中最重要的问题之一。自早期以来,工程师们一直在使用计算机解决这些问题。

我不是历史专家,但显然当时每个人都只是用简单的循环重写他的FORTRAN版本。随后出现了一些标准化,其中包括识别大多数线性代数问题所需的“核心”(基本例程)。这些基本操作然后在规范中被标准化,称为:基本线性代数子程序(BLAS)。工程师可以在他们的代码中调用这些标准的、经过良好测试的BLAS例程,使他们的工作变得更加容易。

BLAS:

BLAS从级别1(第一个定义标量向量和向量向量操作的版本)演变到级别2(向量矩阵操作)到级别3(矩阵矩阵操作),并提供越来越多的“核心”,使更多基本的线性代数操作被标准化。最初的FORTRAN 77实现仍然可以在Netlib网站上获取。

向更好的性能:

因此,多年来(尤其是在BLAS级别1和级别2发布之间:80年代早期),硬件发生了变化,随着向量运算和缓存层次结构的出现。这些演变使得BLAS子程序的性能大幅提高。然后出现了不同的供应商,他们的BLAS例程实现越来越高效。

我不知道所有的历史实现(那时我还没有出生或还是孩子),但最著名的两个实现于2000年初推出:英特尔MKL和GotoBLAS。你的Matlab使用了英特尔MKL,这是一个非常好的、优化的BLAS,这解释了你看到的很好的性能。

矩阵乘法的技术细节:

那么为什么Matlab(MKL)在dgemm(双精度通用矩阵-矩阵乘法)上如此快?简单地说:因为它使用了向量化和良好的数据缓存。在更复杂的术语中:请参见Jonathan Moore提供的文章

基本上,当您在提供的C++代码中执行乘法时,您并不是很友好于缓存。由于我怀疑您创建了一个指向行数组的指针数组,您在内部循环中对“matice2”的第k列进行访问:matice2[m][k]非常缓慢。实际上,当您访问matice2[0][k]时,您必须获取您的矩阵的数组0的

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 算法的人是错误的,这两个算法都不能在实践中实现,因为上面提到的硬件考虑。


3
我刚刚观看了Scott Meyers的一段视频,内容讲述了缓存大小和将数据适配到缓存行大小的重要性,以及在多线程解决方案中可能出现的问题,即使源代码中没有共享数据,但在硬件/核心线程级别上仍可能存在共享数据的情况。视频链接:https://youtu.be/WDIkqP4JbkE - WillC
链接到文章已经失效了。您能否更新一下呢? - Our
我们的链接现在对我来说有效... - reverse_engineer

97

下面是我使用MATLAB R2011a和Parallel Computing Toolbox在一台装有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使用高度优化的矩阵乘法库,这就是为什么纯MATLAB矩阵乘法如此快的原因。gpuArray版本使用MAGMA

R2014a更新在一台搭载Tesla K20c的机器上,使用新的timeitgputimeit函数:

>> 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

(注: 在某个时刻(我忘记具体是什么时候)gpuArray从MAGMA转换到了cuBLAS —— 但仍然用于一些 gpuArray 操作中)

使用 R2022a 更新 WIN64 机器,拥有 32 个物理核心和 A100 GPU:

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

这为什么很重要? - Mad Physicist
为什么这很重要呢?我试图针对 MATLAB 在各种情况下使用的库提供一些见解,以解释为什么 MATLAB 的性能很好 - 即因为它使用高度优化的数值库。 - Edric
5
哇,感谢您这些年来的更新! - dpdp

42

这就是为什么,MATLAB不会像你在C++代码中那样通过循环遍历每个元素进行简单矩阵乘法。

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


20

7
MATLAB使用英特尔MKL库,该库提供了对BLAS/LAPACK例程的优化实现:https://dev59.com/13LYa4cB1Zd3GeqPTA0y#16723946 - Amro

13

答案是 LAPACKBLAS 库使 MATLAB 在矩阵运算方面非常快,而不是来自 MATLAB 团队的任何专有代码。

在你的 C++ 代码中使用 LAPACK 和/或 BLAS 库进行矩阵运算,你应该能够获得类似于 MATLAB 的性能。这些库应该在任何现代系统上都可以免费使用,其中一些部分是在学术界经过数十年的发展而成。请注意,有多个实现版本,包括一些闭源版本,例如 Intel MKL

关于 BLAS 如何获得高性能的讨论 可在此处获得。


顺便提一下,从C语言直接调用LAPACK库是我个人的痛苦经历(但是非常值得)。你需要非常精确地阅读文档。


9

在执行矩阵乘法时,您使用了朴素乘法方法,其时间复杂度为O(n^3)

存在一种矩阵乘法算法,其时间复杂度为O(n^2.4)。这意味着当n=2000时,您的算法需要比最佳算法多计算约100倍。
您应该确实查看维基百科页面以获取有关高效实现的更多信息。


而MATLAB可能使用了这种算法,因为10241024矩阵乘法的时间小于20482048矩阵乘法时间的8倍!做得好MATLAB团队。 - Renaud
6
尽管具有理论优势,但我相当怀疑他们使用“高效”乘法算法。即使是Strassen算法也存在实现困难,而您可能已经了解到的Coppersmith-Winograd算法就目前而言根本不切实际。此外,还有相关的SO线程:https://dev59.com/GXTYa4cB1Zd3GeqPwpkq - Ernir
那个算法仅适用于非常大的矩阵。 - user1911226
@Renaud. 这是相对恒定开销的定义。 - Mad Physicist

7
根据您的Matlab版本,我相信它可能已经在使用您的GPU了。另外,Matlab会跟踪矩阵的许多属性;例如是否为对角线矩阵、共轭矩阵等,并根据此专门优化其算法。也许它是基于您传递给它的零矩阵进行优化的,或者类似的情况?也许它正在缓存重复的函数调用,从而破坏了您的时间记录?也许它优化掉了重复未使用的矩阵乘积?
为了防止发生这种情况,请使用随机数矩阵,并确保通过将结果打印到屏幕或磁盘等方式强制执行。

4
作为一个频繁使用机器学习的用户,我可以告诉你,他们还没有使用GPGPU。Matlab的新版本确实使用了SSE1/2指令集(终于),但是我已经做了测试。一个执行逐元素相乘操作的MexFunction比A.*B命令运行快两倍。因此,问题提出者几乎肯定在某些地方出错了。 - KitsuneYMG
6
Matlab与Parallel Computing Toolbox可以使用CUDA GPU,但需要显式地将数据发送到GPU中。 - Edric
我使用 M1 = single(rand(1024,1024)*255); M2 = single(rand(1024,1024)*255); 和 M3 = M1 * M2; ... 然后将其写入浮点数的二进制文件,所有操作都非常快速完成。 - Wolf

4
MATLAB使用由英特尔提供的高度优化的LAPACK实现,称为Intel Math Kernel Library(Intel MKL),特别是dgemm function函数。该库利用处理器功能,包括SIMD指令和多核处理器。他们没有记录使用哪种具体算法。如果您从C ++调用Intel MKL,则应看到类似的性能。
我不确定MATLAB用于GPU乘法的库是什么,但可能类似于nVidia CUBLAS

1
你说得没错,但是你看过这个答案吗?然而,IPP并不是MKL,而MKL在线性代数性能方面比IPP要好得多。此外,IPP在最近的版本中已经弃用了它们的矩阵数学模块。 - chappjc
抱歉,我的意思是 MKL 而不是 IPP。 - gregswiss
你说得对,另一个答案已经涵盖了它。它太冗长了,我错过了它。 - gregswiss

4
“为什么Matlab在处理xxx方面比其他程序更快?”的一般答案是Matlab具有许多内置的、经过优化的函数。
使用的其他程序通常没有这些功能,因此人们会应用自己的创造性解决方案,这些解决方案比专业优化代码慢得令人惊讶。
这可以有两种解释:
1)普遍/理论上:Matlab并没有显著地更快,只是你的基准测试做错了。
2)现实中:对于这些内容,Matlab在实践中更快,因为像C++这样的语言很容易被用于低效的方式。

8
他正在比较MATLAB的速度和他在两分钟内编写的函数的速度。我可以在10分钟内编写出更快的函数,或者在两个小时内编写出更快的函数。MATLAB工程师花费了超过两个小时来提高矩阵乘法的速度。 - gnasher729

3
由于MATLAB最初是为数值线性代数(矩阵操作)而开发的编程语言,具有专门开发的用于矩阵乘法的库。现在,MATLAB还可以额外使用GPU(图形处理器)
如果我们看一下您的计算结果:
             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简史 克里夫·莫勒是新墨西哥大学计算机科学系主席,他在1970年代末开始开发MATLAB。他设计它是为了让他的学生能够使用LINPACK(用于执行数值线性代数的软件库)和EISPACK(用于数值计算线性代数的软件库),而不必学习Fortran。它很快传播到其他大学,并在应用数学界找到了强大的听众。1983年,工程师杰克·利特在莫勒访问斯坦福大学期间接触到了它。他意识到它的商业潜力,与莫勒和史蒂夫·班格特合作。他们用C重写了MATLAB,并于1984年创立了MathWorks以继续其发展。这些重写的库被称为JACKPAC。2000年,MATLAB被重写以使用一组更新的矩阵操作库,LAPACK(是用于数值线性代数的标准软件库)。 来源

CUDA C是什么?

CUDA C还使用了专门用于矩阵乘法的库,例如OpenGL(开放图形库)。它还使用GPU和Direct3D(在MS Windows上)。
CUDA平台旨在与诸如C、C ++和Fortran等编程语言一起使用。这种可访问性使得并行编程专家更容易使用GPU资源,与之前的API(例如Direct3D和OpenGL)相比,这需要先进的图形编程技能。此外,CUDA支持编程框架,例如OpenACC和OpenCL。
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,该算法的计算时间从 CPU 上的一分钟以上降至 GPU 上的不到 10 秒,降低了 7.5 倍。对数刻度图显示,CPU 实际上在小网格大小时更快。然而,随着技术的发展和成熟,GPU 解决方案越来越能够处理较小的问题,这是我们预计会继续的趋势。

来源

来自 CUDA C 编程指南介绍:

由于市场对实时高清 3D 图形的需求不可满足,可编程图形处理器或 GPU 已经演变成高度并行、多线程、众核心的处理器,具有巨大的计算能力和非常高的内存带宽,如 图1图2 所示。 图1。CPU 和 GPU 的浮点运算次数每秒。 enter image description here 图2。CPU 和 GPU 的内存带宽。 enter image description here GPU 的浮点性能之所以比 CPU 更强,是因为 GPU 专门进行计算密集、高度并行的计算,这正是图形渲染所需的,因此设计中更多的晶体管用于数据处理而不是数据缓存和流量控制,如 图3 所示。 图3。GPU 更多晶体管用于数据处理。 enter image description here 更具体地说,GPU 特别适合解决可以表示为数据并行计算的问题,即同一程序在许多数据元素上并行执行,并具有高算术强度,即算术运算与内存操作的比值。由于每个数据元素都执行相同的程序,因此对复杂流控制的要求较低,并且由于它在许多数据元素上执行并具有高算术强度,因此可以通过计算来隐藏内存访问延迟,而不是使用大型数据缓存。
数据并行处理将数据元素映射到并行处理线程。许多处理大数据集的应用程序可以使用数据并行编程模型加速计算。在 3D 渲染中,大量像素和顶点被映射到并行线程。类似地,图像和媒体处理应用程序,如渲染图像后处理、视频编码和解码、图像缩放、立体视觉和模式识别等,可以将图像块和像素映射到并行处理线程。实际上,许多图像渲染和处理领域之外的算法都通过数据并行处理进行加速,从一般的信号处理或物理模拟到计算金融或计算生物学。 来源

高级阅读


一些有趣的事实:
引用:this answer">这个答案中提到:“我写了一个 C++ 矩阵乘法,速度和 Matlab 一样快,但需要一些小心。(在 Matlab 开始使用 GPU 之前)。”

2
那句话并不是“事实”,而是空洞的吹嘘。自从他发布那篇文章以来,已经有好几个人请求他提供代码,但是却没有看到任何代码。 - Cris Luengo
1
你关于GPU计算速度的描述完全没有回答问题。我们都知道,128个小核心可以比2个大核心更快地完成相同、单调的工作。“现在MATLAB也可以额外使用GPU(图形处理器)。”是的,但不是默认情况下。正常的矩阵乘法仍然使用BLAS。 - Cris Luengo
@CrisLuengo,好的,这不是事实!也许你对他的“吹嘘”是正确的 - 我们不知道这一点,也不知道他为什么不回答。对于第二条评论:关于GPU上的计算描述回答了问题,因为在线性代数中的矩阵乘法中使用浮点运算。也许并不是所有人都能理解,但我认为他们必须先理解这些基础知识,否则,他们在阅读有关矩阵的文章之前必须先学习这些基础知识。如果其他人给我写信询问这方面的问题,那么我会加上这些细节。谢谢! - Bharata
@CrisLuengo,我写了单词“additionally”。它的意思是:它可以被使用。这也意味着普通矩阵乘法仍然使用软件库。你认为我需要改变我的帖子以使其更易理解吗?感谢您的评论! - Bharata

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接