矩阵乘法算法详解

矩阵乘法算法详解

引用

CSDN

1.

https://blog.csdn.net/axecute/article/details/145141345

矩阵乘法在许多数值算法中都是一个核心操作，因此人们投入了大量工作来提高矩阵乘法算法的效率。矩阵乘法在计算问题中的应用涉及许多领域，包括科学计算和模式识别，以及看似无关的问题，例如计算图中的路径。已经设计了许多不同的算法，用于在不同类型的硬件上对矩阵进行乘法，包括并行和分布式系统，其中计算工作分布在多个处理器（可能通过网络）上。

迭代算法

矩阵乘法的定义是，如果 n × m 矩阵 A 和 m × p 矩阵 B 的 C = AB，则 C 是一个 n × p 矩阵，其中包含条目

由此，可以构建一个简单的算法，该算法遍历从 1 到 n 的索引 i 和从 1 到 p 的索引 j，使用嵌套循环计算上述内容：

• 输入：矩阵 A 和 B
• 设 C 为适当大小的新矩阵
• 对于 i，从 1 到 n：
• 对于 j 从 1 到 p：
• 设 sum = 0
• 对于 k，从 1 到 m：
• 设置总和 ← sum + Aik × Bkj
• 集合 Cij ← 和
• 返回 C

该算法需要时间 Θ（nmp）（采用渐近表示法）。为了进行算法分析，一种常见的简化是假设输入都是大小为 n × n 的方阵，在这种情况下，运行时间为 Θ（n3），即维度大小的立方。

缓存行为

迭代矩阵乘法中的三个循环可以任意交换彼此，而不会影响正确性或渐近运行时间。但是，由于算法的内存访问模式和缓存使用，该顺序可能会对实际性能产生相当大的影响；哪种顺序最好还取决于矩阵是按行优先顺序、列优先顺序还是两者混合存储。

特别是，在完全关联缓存由每行缓存行的 M 字节和 b 字节组成的理想化情况下（即⁠M/b⁠cache 行），上述算法对于以行优先顺序存储的 A 和 B 来说是次优的。 什么时候n > ⁠M/b⁠，内部循环的每次迭代（同时扫描 A 的一行和一列 B）在访问 B 的元素时都会产生缓存未命中。这意味着在最坏的情况下，该算法会产生 Θ（n3） 缓存未命中。截至 2010 年，与处理器相比，内存的速度是这样的，缓存未命中，而不是实际计算，主导了大型矩阵的运行时间。

在行优先布局中，A 和 B 迭代算法的最佳变体是平铺版本，其中矩阵被隐式划分为大小为 √M x √M 的方形瓦片：

• 输入：矩阵 A 和 B
• 设 C 为适当大小的新矩阵
• 选择图块大小 T = Θ（√M）
• 对于 I，从 1 到 n，步长为 T：
• 对于 J，从 1 到 p，步长为 T：
• 对于 K，从 1 到 m，以 T 为步长：
• 乘以 AI：I+T、K：K+T 和 BK：K+T， J：J+T 转换为 CI：I+T， J：J+T，即：
• 对于 i 从 I 到 min（I + T， n）：
• 对于 j 从 J 到 min（J + T， p）：
• 设 sum = 0
• 对于 k 从 K 到 min（K + T， m）：
• 设置总和 ← sum + Aik × Bkj
• 集合 Cij ← Cij + 总和
• 返回 C

在理想化缓存模型中，此算法仅产生Θ(⁠注3/b √M⁠)缓存未命中；在现代机器上，除数 b √M 相当于几个数量级，因此实际计算主导运行时间，而不是缓存未命中。

分而治之算法

迭代算法的替代方案是矩阵乘法的分而治之算法。这依赖于块分区

它适用于维度为 2 的幂的所有方阵，即形状为 2n × 2n 表示某些 n。矩阵积现在是

它由子矩阵对的八次乘法组成，然后是一个加法步骤。分而治之算法使用标量乘法 c11 = a11b11 作为其基本情况，递归计算较小的乘法。

该算法作为 n 函数的复杂度由递归]

考虑对大小为 n/2 和 Θ（n2） 的矩阵的 8 次递归调用，以对生成的四对矩阵进行元素求和。分而治之递归的主定理的应用表明，这种递归具有解 Θ（n3），与迭代算法相同。]

非二方矩阵

该算法的一个变体适用于任意形状的矩阵，并且在实践中速度更快，它将矩阵一分为二，而不是四个子矩阵，如下所示。现在，拆分矩阵意味着将其分成大小相等的两个部分，或者在奇数维度的情况下尽可能接近相等大小。

• 输入：大小为 n × m 的矩阵 A，大小为 m × p 的矩阵 B。
• 基本情况：如果 max（n， m， p） 低于某个阈值，则使用迭代算法的展开版本。
• 递归情况：
• 如果 max（n， m， p） = n，则水平拆分 A：
• 否则，如果 max（n， m， p） = p，则垂直拆分 B：
• 否则，max（n， m， p） = m。垂直拆分 A，水平拆分 B：

缓存行为

递归矩阵乘法的缓存未命中率与平铺迭代版本的缓存未命中率相同，但与该算法不同的是，递归算法是缓存无关的：不需要调整参数即可获得最佳缓存性能，并且在多编程环境中表现良好，由于其他进程占用缓存空间，缓存大小实际上是动态的。（简单的迭代算法也是缓存无关的，但如果矩阵布局不适应算法，在实践中会慢得多。

在具有 M 行理想缓存（每行大小为 b 字节）的机器上，此算法产生的缓存未命中数受 的限制

亚三次算法

更多信息：矩阵乘法的计算复杂度

存在比直接算法提供更好的运行时间的算法。首先被发现的是 Strassen 算法，由 Volker Strassen 于 1969 年设计，通常被称为“快速矩阵乘法”。它基于一种将两个 2 × 2 矩阵相乘的方法，只需要 7 次乘法（而不是通常的 8 次），但代价是几个额外的加法和减法运算。递归应用此算法会得到一个乘法成本为

Strassen 算法更复杂，与朴素算法相比，数值稳定性降低，但在 n > 100 左右的情况下，它会更快并出现在多个库中，例如 BLAS。它对于精确域（如有限域）上的大型矩阵非常有用，其中数值稳定性不是问题。

由于 Strassen 算法实际上用于实际数值软件和计算机代数系统，因此改进隐藏在大 O 表示法中的常数有其优点。下表比较了基于 2x2 块矩阵的递归乘法的改进版本的关键方面，通过 7 个块矩阵乘法。 照常N给出矩阵的维度，M指定内存大小。

类似 Strassen 的递归 2x2 块矩阵乘法的进展

众所周知，具有 2x2 块矩阵步长的类似 Strassen 的算法至少需要 7 次块矩阵乘法。1976 年，Probert表明，这种算法至少需要 15 次加法（包括减法），但是，一个隐藏的假设是块和 2x2 块矩阵以相同的基表示。Karstadt 和 Schwartz 以不同的基数进行计算，并将 3 次加法换取更便宜的基变换。他们还证明，使用不同的基数每步不能低于 12 次加法。在随后的工作中，Beniamini 等人将这种碱基变化技巧应用于比 2x2 块矩阵更一般的分解，并改进了其运行时间的前导常数。

在理论计算机科学中，Strassen 算法在渐近复杂度方面的改进程度是一个悬而未决的问题。矩阵乘法指数，通常表示为w是任何n*n字段上的矩阵可以使用field 操作。当前最佳边界w是w <2.371552，由 Williams、Xu、Xu 和 周 编写。该算法与该研究领域的所有其他最新算法一样，是 Don Coppersmith 和 Shmuel Winograd 于 1990 年给出的 Coppersmith-Winograd 算法的推广。这些算法的概念思想类似于 Strassen 的算法：设计了一种方法，用于将两个 k × k 矩阵乘以少于 k3 次的乘法，并且该技术是递归应用的。但是，大 O 表示法隐藏的常数系数如此之大，以至于这些算法仅对太大而无法在当今计算机上处理的矩阵有价值。Victor Pan 提出了所谓的可行次三次矩阵乘法算法，其指数略高于 2.77，但返回的隐藏常数系数要小得多。

Freivalds 算法是一种简单的蒙特卡洛算法，给定矩阵 A、B 和 C，如果 AB = C，则在 Θ（n2） 时间内进行验证。

AlphaTensor

2022 年，DeepMind 推出了 AlphaTensor，这是一种神经网络，它使用单人游戏类比发明了数千种矩阵乘法算法，包括一些以前由人类发现的算法和一些尚未发现的算法。操作仅限于非交换接地场（正态算术）和有限域（Mod 2 算术）。发现的最佳“实用”（矩阵乘法张量的显式低秩分解）算法在 O（n2.778） 中运行。找到此类张量（及更高）的低秩分解是 NP 困难的；即使在交换场中，即使对于 3x3 矩阵，最优乘法也是未知的。在 4x4 矩阵上，AlphaTensor 出乎意料地发现了一个具有 47 个乘法步骤的解决方案，比 1969 年 Strassen 算法所需的 49 个步骤有所改进，尽管仅限于 mod 2 算术。同样，AlphaTensor 用 96 个步骤而不是 Strassen 的 98 个步骤求解 5x5 矩阵。基于存在这种改进的惊人发现，其他研究人员很快就找到了类似的独立 4x4 算法，并分别将 Deepmind 的 96 步 5x5 算法在 mod 2 算术中降低到 95 步，在普通算术中降低到 97一些算法是全新的：例如，（4， 5， 5） 在正常和 mod 2 算术中都从基线 80 改进到 76 步。

并行和分布式算法

共享内存并行性

前面草拟的分而治之算法可以通过两种方式对共享内存多处理器进行并行化。这些是基于这样一个事实，即

可以彼此独立地执行，四个求和也可以（尽管算法需要在进行求和之前“连接”乘法）。利用问题的完全并行性，人们获得了一种可以用 fork-join 样式伪代码表示的算法：

过程乘法 （C， A， B）：

• 基本情况：如果 n = 1，则将 c11 ← a11 × b11（或乘以一个小块矩阵）。
• 否则，为形状为 n × n 的新矩阵 T 分配空间，然后：
• 将 A 分割成 A11、A12、A21、A22。
• 将 B 划分为 B11、B12、B21、B22。
• 将 C 划分为 C11、C12、C21、C22。
• 将 T 划分为 T11、T12、T21、T22。
• 并行执行：
• 叉乘（C11， A11， B11）。
• 叉乘（C12， A11， B12）。
• 叉乘（C21， A21， B11）。
• 叉乘（C22， A21， B12）。
• 叉乘（T11， A12， B21）。
• 叉乘（T12， A12， B22）。
• 叉乘（T21， A22， B21）。
• 叉乘（T22， A22， B22）。
• 加入（等待并行分叉完成）。
• add（C， T） 的
• 解除分配 T。

过程 add（C， T） 将 T 添加到 C 中，按元素：

• 基本情况： 如果 n = 1，则设置 c11 ← c11 + t11 （或做一个短循环，也许展开）。
• 否则：
• 将 C 划分为 C11、C12、C21、C22。
• 将 T 划分为 T11、T12、T21、T22。
• 并行：
• 分叉 add（C11， T11）。
• 分叉 add（C12， T12）。
• 叉子 add（C21， T21）。
• 分叉 add（C22， T22）。
• 加入。

这里，fork 是一个关键字，表示计算可以与函数调用的其余部分并行运行，而 join 则等待所有先前“fork”的计算完成。partition 仅通过指针操作来实现其目标。

该算法的临界路径长度为 Θ（log2 n） 步，这意味着在具有无限数量处理器的理想机器上需要花费大量时间；因此，在任何真实计算机上，它的最大可能加速比为 Θ（n3/log2 n）。由于将数据移入和移出临时矩阵 T 的固有通信成本，该算法并不实用，但更实用的变体可以在不使用临时矩阵的情况下实现 Θ（n2） 加速。

通信避免和分布式算法

在具有分层内存的现代架构上，加载和存储输入矩阵元素的成本往往主导算术成本。在单台机器上，这是在 RAM 和缓存之间传输的数据量，而在分布式内存多节点机器上，它是节点之间传输的数据量；在任何一种情况下，它都称为通信带宽。使用三个嵌套循环的朴素算法使用 Ω（n3） 通信带宽。

Cannon 算法，也称为 2D 算法，是一种避免通信的算法，它将每个输入矩阵划分为一个块矩阵，该矩阵的元素是大小为 √M/3 x √M/3 的子矩阵，其中 M 是快速内存的大小。然后，在块矩阵上使用朴素算法，完全在快速内存中计算子矩阵的乘积。这将通信带宽降低到 O（n3/√M），这是渐近最优的（对于执行 Ω（n3） 计算的算法）。

在分布式设置中， 个处理器以 x 2D 网格排列，可以将结果的一个子矩阵分配给每个处理器，并且可以计算乘积每个处理器传输 字，假设每个节点存储最少的 元素，这是渐近最优的。在分布式设置中，p proces.sors 由 √p 2D 网格排列在 √p 中，结果的一个子矩阵可以分配给 each 处理器，并且乘积可以在每个处理器传输 O（n2/√p） 字的情况下计算，假设每个节点存储最小 O（n2/p）元素。 这可以通过 3D 算法p1/3来改进，该算法将处理器排列在 3D 立方体网格中，将两个输入子矩阵的每个乘积分配给单个处理器。结果子矩阵然后通过对每一行执行缩减来生成。<b.2031 这个算法rithm 传输 O（n2/p2/3）字每个处理器，这是渐近最优的。<b.19> 但是，这需要将每个输入矩阵元素 p 复制 1/3 次，因此需要

网格算法

网格上的乘法有多种算法。对于使用 2D Cannon 算法在标准二维网格上乘以两个 n×n，可以在 3n-2 步中完成乘法，尽管对于重复计算，这可以减少到这个数字的一半。标准数组效率低下，因为来自两个矩阵的数据不会同时到达，必须用零填充。

在两层交叉布线网格上，结果甚至更快，其中只需要 2个 n-1 步骤。重复计算的性能进一步提高，从而达到 100% 的效率。交叉布线网格阵列可以被视为非平面（即多层）处理结构的特例。