matlab - MATLAB 中的矩阵乘法时间复杂度

Question

有谁知道 MATLAB 使用哪种算法进行矩阵乘法及其时间复杂度是多少？

Answer 1

为了完整性——如 this thread 中所述, Matlab 使用 DGEMM (双一般矩阵乘法)例程来自 BLAS (基本线性代数子程序)。

请注意，BLAS 的实现不是单一的——它针对特定的处理器架构进行了调整。因此，如果不查明正在使用哪个版本的 BLAS，您就无法绝对确定您的机器上正在使用哪种算法。

BLAS 的规范指定了每个子程序的输入和输出，并为每个子程序的输出提供了可接受的错误范围。实现可以自由使用他们喜欢的任何算法，只要它们遵循规范即可。

BLAS 的引用实现使用 block matrix multiplication algorithm在 DGEMM 中，将两个 n x n 矩阵相乘的时间复杂度为 O(n^3)。我认为可以合理地假设大多数 BLAS 实现将或多或少地遵循引用实现。

请注意，它不使用朴素的矩阵乘法算法

for i = 1:N
    for j = 1:N
        for k = 1:N
            c(i,j) = c(i,j) + a(i,k) * b(k,j);
        end
    end
end

这是因为，通常情况下，整个矩阵不会适合 local memory .如果数据不断地移入和移出本地内存，算法就会变慢。 block 矩阵算法将操作分解成小块，这样每个 block 都足够小以适合本地内存，从而减少移入和移出内存的次数。

存在渐进更快的矩阵乘法算法，例如 Strassen algorithm或 Coppersmith-Winograd algorithm它的速度比 O(n^3) 稍快。但是，它们通常不感知缓存并忽略局部性 - 这意味着数据需要不断地在内存中分流，因此对于大多数现代架构而言，整体算法实际上比优化的 block 矩阵乘法算法慢。

维基百科指出，Strassen 算法可以在单核 CPU 上为大于几千的矩阵提供加速，但是加速可能在 10% 左右，BLAS 的开发人员可能认为不值得对于这种罕见的情况(也就是说，1996 年的 this paper 声称对于 n 大约 200 以上的 DGEMM 速度提高了大约 10%——尽管我不知道那是多么的过时)。另一方面，Coppersmith-Winograd 算法“只为现代硬件无法处理的大矩阵提供优势”。

所以答案是 Matlab 使用一种朴素但高效且缓存感知的算法来获得其超快的矩阵乘法。

---

我通过创建一些视频来更新这个答案，这些视频展示了与朴素算法相比， block 矩阵乘法算法的局部性。

在以下每个视频中，我们都在可视化两个 8x8 矩阵 A 和 B 的乘积以创建乘积 C = <强>A x B。黄色突出显示表示在算法的每个步骤中正在处理每个矩阵 A、B 和 C 中的哪个元素。您可以看到 block 矩阵乘法如何一次仅对矩阵的小块起作用，并多次重复使用这些 block 中的每一个，从而最大限度地减少数据必须移入和移出本地内存的次数.

• Simple matrix multiplication
• Block matrix multiplication, block size = 2
• Block matrix multiplication, block size = 4