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

Question

有谁知道MATLAB用于矩阵乘法的算法是什么,它的时间复杂度是多少？

Answer 1

为了完整性 - 如本主题所述,Matlab使用BLAS(基本线性代数子程序)中的DGEMM(双通用矩阵乘法)例程.

请注意,BLAS没有一个单独的实现 - 它针对特定的处理器体系结构进行了调整.因此,如果没有找到正在使用的BLAS版本,您无法确定您的计算机上使用了哪种算法.

BLAS规范规定了每个子程序的输入和输出,并为每个子程序的输出提供了可接受的误差范围.实现可以随意使用他们喜欢的任何算法,只要它们遵循规范.

BLAS的参考实现使用块矩阵乘法算法,DGEMM其具有时间复杂度O(n ^ 3),用于乘以两个n × 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)略快的速率.但是,它们通常不会识别缓存并忽略局部性 - 这意味着数据不断需要在内存中分流,因此对于大多数现代架构而言,整体算法实际上比优化的块矩阵乘法算法慢.

维基百科指出,Strassen算法可以在单个核心CPU上为大于几千的矩阵提供加速,但是加速可能在10%左右左右,而且BLAS的开发人员可能认为不值得这种罕见的情况下(说,本文从1996年声称的10%左右,在速度提升DGEMM的ň以上约200 -尽管我不知道怎么搞的是过时的).另一方面,Coppersmith-Winograd算法"只为矩阵提供了一个优势,使得它们无法被现代硬件处理".

所以答案是Matlab使用一种天真但有效的缓存感知算法来获得其快速的矩阵乘法.

我通过创建一些演示块矩阵乘法算法的局部性的视频来更新这个答案,与天真算法相比较.

在下面每一个视频,我们可视化两个8×8矩阵的乘法甲和乙创建产品Ç = 甲 X 乙.黄色突出显示指示在算法的每个步骤中处理每个矩阵A,B和C中的哪个元素.您可以看到块矩阵乘法一次只对矩阵的小块有效,并且多次重复使用这些块中的每一个,以便最小化数据必须移入和移出本地存储器的次数.

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