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

Question

有人知道MATLAB用于矩阵乘法的算法及其时间复杂度吗？

- pikachu

---

这个问题在2009年在Matlab Central上得到了回答，链接在这里（特别是查看Tim Davis的第二个回复）。不确定自那以后有没有什么变化... - Colin T Bowers

Answer 1

为了完整起见--如此线程所提到的，Matlab使用BLAS（基本线性代数子程序）中的DGEMM（双精度一般矩阵乘法）例程。

请注意，BLAS没有一个单一的实现--它是针对特定处理器架构进行调整的。因此，如果不找出正在使用哪个版本的BLAS，则无法确定在您的计算机上使用哪个算法。

BLAS的规范指定了每个子程序的输入和输出，并为每个子程序的输出提供可接受的误差边界。实现可以自由地使用任何算法，只要遵循规范即可。

BLAS的参考实现在DGEMM中使用块矩阵乘法算法，其时间复杂度为O(n^3)，用于计算两个n x n矩阵的乘积。我认为可以合理地假设大多数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

这是因为通常整个矩阵无法适应本地内存。如果数据不断地在本地内存中移进移出，算法将会变慢。块矩阵算法将操作分解成小块，使得每个块都足够小以适应本地内存，从而减少了内存移进移出的次数。

存在渐近更快的矩阵乘法算法，例如Strassen算法或Coppersmith-Winograd算法，其速率略快于O(n^3)。然而，它们通常不考虑缓存并忽略局部性 - 这意味着数据需要不断地在内存中移动，因此对于大多数现代架构来说，整体算法实际上比优化的块矩阵乘法算法更慢。

维基百科指出，对于矩阵大小大于几千的单核CPU，Strassen算法可能会提供加速，但加速可能只有约10％左右，并且BLAS的开发人员可能不认为这种情况值得（话说，这篇1996年的this paper声称，在n大约为200以上时，与DGEMM相比，速度提高了约10％-尽管我不知道它是否过时）。另一方面，Coppersmith-Winograd算法“仅在矩阵非常大以至于现代硬件无法处理时才提供优势”。

因此，答案是Matlab使用一个朴素但高效且缓存感知的算法来获得其快速的矩阵乘法。

---

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

在以下每个视频中，我们将可视化两个8x8矩阵A和B的乘法，以创建乘积C = A x B。黄色高亮显示了算法每一步中正在处理的每个矩阵A、B和C中的元素。您可以看到块矩阵乘法仅逐个处理矩阵的小块，并多次重复使用每个块，以使数据在本地内存中移入和移出的次数最小化。

• 简单矩阵乘法
• 块矩阵乘法，块大小=2
• 块矩阵乘法，块大小=4

- Chris Taylor

---

好的回答 +1。我修改了你的措辞“比O(n^3)更少的操作”，因为严格来说，两个例程都可以是O(n^3)，但一个例程可以比另一个例程有更少的操作。 - Colin T Bowers