为什么 MATLAB 的矩阵乘法速度如此之快？

Question

我正在使用 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 相同，但 2048 与 1024 有一个巨大的跳跃。

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

Answer 1

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

历史：

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

我不是历史专家，但显然在那时，每个人都只是用简单的循环重写了他的 FORTRAN 版本。随后出现了一些标准化，并确定了大多数线性代数问题需要解决的“内核”（基本例程）。这些基本操作随后在称为基本线性代数子程序 (BLAS) 的规范中标准化。然后，工程师可以在他们的代码中调用这些经过良好测试的标准 BLAS 例程，从而使他们的工作更加轻松。

BLAS：

BLAS从level 1（定义标量向量和向量向量运算的第一个版本）到level 2（向量矩阵运算）再到level 3（矩阵矩阵运算），并提供了越来越多的“内核”因此，越来越多的基本线性代数运算标准化。 Netlib's website 上仍然提供原始 FORTRAN 77 实现。

实现更好的性能：

所以多年来（特别是在 BLAS 1 级和 2 级版本之间：80 年代初），随着向量操作和缓存层次结构的出现，硬件发生了变化。这些演变使得大幅提高 BLAS 子程序的性能成为可能。随后，不同的供应商推出了他们的 BLAS 例程实现，这些例程越来越高效。

我不知道所有的历史实现（那时我还不是孩子），但有两个最著名的实现是在 2000 年代初问世的：英特尔 MKL 和 GotoBLAS。您的 Matlab 使用英特尔 MKL，这是一个非常好的优化 BLAS，这解释了您所看到的出色性能。

矩阵乘法的技术细节：

那么为什么 Matlab（MKL）在 dgemm（双精度通用矩阵-矩阵乘法）时如此之快？简单来说：因为它使用矢量化和良好的数据缓存。更复杂的术语：请参阅 Jonathan Moore 提供的 article。

基本上，当您在提供的 C++ 代码中执行乘法运算时，您根本不适合缓存。由于我怀疑您创建了一个指向行数组的指针数组，因此您在内部循环中对“matice2”的第 k 列的访问：matice2[m][k] 非常慢。事实上，当您访问matice2[0][k] 时，您必须获取矩阵的数组 0 的第 k 个元素。然后在下一次迭代中，您必须访问matice2[1][k]，它是另一个数组（数组 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 算法是错误的，由于上述硬件考虑，这两种算法在实践中都无法实现。

Answer 2

这是我在装有 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

在具有 16 个物理内核和 Tesla V100 的 WIN64 机器上使用 R2018b 更新：

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

（注意：在某些时候（我忘了具体时间）gpuArray 从 MAGMA 切换到 cuBLAS - MAGMA 仍用于某些 gpuArray 操作）

Answer 3

This is why。 MATLAB 不会像在 C++ 代码中那样通过循环遍历每个元素来执行简单的矩阵乘法。

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

Answer 4

Matlab 前段时间合并了 LAPACK，所以我假设他们的矩阵乘法至少使用了这么快的东西。 LAPACK 源代码和文档很容易获得。

您还可以查看 Goto 和 Van De Geijn 的论文“Anatomy of High-Performance Matrix 乘法”http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

Answer 5

答案是 LAPACK 和 BLAS 库使 MATLAB 在矩阵运算方面的速度非常快，而不是 MATLAB 人员的任何专有代码。

在您的 C++ 代码中使用 LAPACK 和/或 BLAS 库进行矩阵运算，您应该获得与 MATLAB 相似的性能。这些库应该可以在任何现代系统上免费使用，并且在学术界已经开发了数十年。请注意，有多种实现方式，包括一些封闭源代码，例如Intel MKL。

关于 BLAS 如何获得高性能的讨论is available here.

---

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

Answer 6

在进行矩阵乘法时，您使用的是简单的乘法方法，需要花费O(n^3) 的时间。

存在采用O(n^2.4) 的矩阵乘法算法。这意味着在n=2000，您的算法需要的计算量大约是最佳算法的 100 倍。
您真的应该查看 wikipedia 页面以了解矩阵乘法，以获取有关实现它的有效方法的更多信息。

Answer 7

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

另一件事; Matlab 跟踪矩阵的许多属性；无论是对角线、赫尔墨斯等，并以此为基础专门化其算法。也许它的专业化是基于你传递给它的零矩阵，或者类似的东西？也许它正在缓存重复的函数调用，这会打乱你的时间？也许它优化了重复未使用的矩阵产品？

为防止此类事情发生，请使用随机数矩阵，并确保通过将结果打印到屏幕或磁盘等来强制执行。

Answer 8

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

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

这可以用两种方式解释：

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

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

Answer 9

MATLAB 使用来自 Intel 的高度优化的 LAPACK 实现，称为 Intel Math Kernel Library (Intel MKL) - 特别是 dgemm function。速度 这个库利用了处理器特性，包括 SIMD 指令和多核处理器。他们没有记录他们使用的具体算法。如果您要从 C++ 调用英特尔 MKL，您应该会看到类似的性能。

我不确定 MATLAB 使用什么库来进行 GPU 乘法，但可能类似于 nVidia CUBLAS。

Answer 10

形成鲜明对比的原因不仅在于 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^3 表示大小为 n 的方阵）

我要求您对其进行计时，以便结果与您之前的结果（在同一台机器上）相当。通过比较，您将准确显示内存访问时间的重要性！

Answer 11

在 C++ 中速度很慢，因为您没有使用多线程。本质上，如果 A = BC，它们都是矩阵，则 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]，多线程是内置的，线程数是可调的。

Answer 12

因为 MATLAB 最初是为数值线性代数（矩阵运算）开发的编程语言，它具有专门为矩阵乘法开发的库。 now MATLAB 也可以为此额外使用 GPUs (Graphics processing unit)。

如果我们看看你的计算结果：

             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 在 1970 年代后期开始开发 MATLAB。他设计它是为了让他的学生可以访问 LINPACK（一个用于执行数值线性代数的软件库）和 EISPACK（一个用于线性数值计算的软件库）代数），而他们不必学习 Fortran。它很快传播到其他大学，并在应用数学界找到了强大的受众。工程师杰克·利特尔 (Jack Little) 在 1983 年莫勒访问斯坦福大学时接触到了它。认识到它的商业潜力，他与莫勒和史蒂夫·班格特一起加入。他们用 C 语言重写了 MATLAB，并于 1984 年创立了 MathWorks 以继续其发展。这些重写的库被称为 JACKPAC。 2000 年，MATLAB 被改写为使用一组更新的矩阵运算库 LAPACK（一个用于数值线性代数的标准软件库）。

Source

什么是 CUDA C

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

CUDA platform 旨在与 C、C++ 和 Fortran 等编程语言配合使用。与之前的 API（如 Direct3D 和 OpenGL）相比，这种可访问性使并行编程专家更容易使用 GPU 资源，后者需要高级图形编程技能.此外，CUDA 还支持 OpenACC 和 OpenCL 等编程框架。

CUDA 处理流程示例：

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

比较 CPU 和 GPU 执行速度

我们运行了一个基准测试，测量了在英特尔至强处理器 X5650 上，然后使用 NVIDIA Tesla C2050 GPU 执行 50 个时间步长，网格大小为 64、128、512、1024 和 2048 所需的时间.

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

Source

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

在市场对实时、高清 3D 图形的永不满足的需求的推动下，可编程图形处理器单元或 GPU 已发展为具有巨大计算能力和非常高内存带宽的高度并行、多线程、多核处理器，如 @ 所示987654354@和Figure 2。

图 1. CPU 和 GPU 的每秒浮点运算次数

图 2。 CPU 和 GPU 的内存带宽

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

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

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

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

Source

进阶阅读

• GPUs (Graphics processing unit)
• MATLAB
• CUDA C Programming Guide
• Using GPUs in MATLAB

• Basic Linear Algebra Subprograms (BLAS)

• Anatomy of High-Performance Matrix Multiplication，来自 Kazushige Goto 和 Robert A. Van De Geijn

---

一些有趣的事实

我已经编写了与 Matlab 一样快的 C++ 矩阵乘法，但它需要一些小心。 （在 Matlab 使用 GPU 之前）。

来自this answer的Сitation。