《算法导论》——矩阵乘法Strassen算法

注：本文为《算法导论》中分治相关内容的笔记。对此感兴趣的读者还望支持原作者。

矩阵乘法

接触过线性代数的读者，对于矩阵乘法想必一定不陌生。若$A=(a_{ij})$和$B=(b_{ij})$是$n\ast n$的方阵，则对$i,j,\dots ,n$，定义乘积$C=A\cdot B$中的元素$c_{ij}$为：

$$c_{ij}=\sum _{k=1}^n{a}_{ik}{b}_{kj}$$

因此，我们可以根据矩阵乘法的定义给出矩阵乘法的伪代码。它接收$n\ast n$的矩阵$A$和$B$，返回它们的乘积——$n\ast n$的矩阵$C$，并且假设每个矩阵都有一个属性$rows$，表示矩阵的行数。

不难看出，由于三重for循环都恰好执行$n$步，而第7行每次执行都花费常量时间。因此，SQUARE-MATRIX-MULTIPLY的时间复杂度为$\theta (n^3)$，即矩阵乘法的朴素实现需要花费$\theta (n^3)$时间。你可能因此认为任何矩阵乘法都要花费$\Omega (n^3)$时间，因为矩阵乘法的自然定义就需要进行这么多次的标量乘法。而在学术界，也的确在很长一段时间内，很少人敢设想一个算法能渐近快于平凡算法SQUARE-MATRIX-MULTIPLY，直至Strassen大神的出现。

算法流程

Strassen算法采用分治法解决矩阵乘积问题，并通过排列组合的技巧使得分治法产生的递归树不那么“茂盛”以减少矩阵乘法的次数。Strassen算法并不直观，它包含4个步骤：

1. 将输入矩阵$A、B$和输出矩阵$C$通过以下方式分解为$\frac{n}{2}\ast \frac{n}{2}$的子矩阵；
   $$A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right],B=\left[\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right],C=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]$$

2. 创建10个$\frac{n}{2}\ast \frac{n}{2}$的矩阵$S_1,S_2,\dots ,S_{10}$，每个矩阵保存步骤1中创建的两个子矩阵的和或差，时间复杂度为$\Theta (n^2)$；

3. 用步骤1中创建的子矩阵和步骤2中创建的10个矩阵，递归地计算7个矩阵积$P_1,P_2,\dots ,P_7$。每个矩阵$P_i$都是$\frac{n}{2}\ast \frac{n}{2}$的；

4. 通过$P_i$矩阵的不同组合进行加减计算，计算出矩阵$C$的子矩阵$C_{11},C_{12},C_{21},C_{22}$，时间复杂度为$\Theta (n^2)$。

是不是感觉很抽象？一顿猛如虎的操作，就能完成矩阵乘积计算了？没错，就是这么秀。接下来，为了帮助大家掌握这种操作，就再看看Strassen算法的细节。在步骤2中，创建如下10个矩阵：

$$S_1=B_{12}-B_{22}$$

$$S_2=A_{11}+A_{12}$$

$$S_3=A_{21}+A_{22}$$

$$S_4=B_{21}-B_{11}$$

$$S_5=A_{11}+A_{22}$$

$$S_6=B_{11}+B_{22}$$

$$S_7=A_{12}-A_{22}$$

$$S_8=B_{21}+B_{22}$$

$$S_9=A_{11}-A_{21}$$

$$S_{10}=B_{11}+B_{22}$$

由于必须进行10次$\frac{n}{2}\ast \frac{n}{2}$的加减法，因此，该步骤花费$\Theta (n^2)$。

在步骤三中，递归地计算7次$\frac{n}{2}\ast \frac{n}{2}$矩阵的乘法，如下所示：

$$P_1=A_{11}\cdot S_1=A_{11}\cdot B_{12}-A_{11}\cdot B_{22}$$

$$P_2=S_2\cdot B_{22}=A_{11}\cdot B_{22}+A_{12}\cdot B_{22}$$

$$P_3=S_3\cdot B_{11}=A_{21}\cdot B_{11}+A_{22}\cdot B_{11}$$

$$P_4=A_{22}\cdot S_4=A_{22}\cdot B_{21}-A_{22}\cdot B_{11}$$

$$P_5=S_5\cdot S_6=A_{11}\cdot B_{11}+A_{11}\cdot B_{22}+A_{22}\cdot B_{11}+A_{22}\cdot B_{22}$$

$$P_6=S_7\cdot S_8=A_{12}\cdot B_{21}+A_{12}\cdot B_{22}-A_{22}\cdot B_{21}-A_{22}\cdot B_{22}$$

$$P_7=S_9\cdot S_{1}0=A_{11}\cdot B_{11}+A_{11}\cdot B_{12}-A_{21}\cdot B_{11}-A_{21}\cdot B_{12}$$

步骤4对步骤3创建的$P_i$矩阵进行加减法运算，计算出$C$的4个$\frac{n}{2}\ast \frac{n}{2}$的子矩阵。

$$C_{11}=P_5+P_4-P_2+P_6=A_{11}\cdot B_{11}+A_{12}\cdot B_{21}$$

$$C_{12}=P_1+P_2=A_{11}\cdot B_{12}+A_{12}\cdot B_{22}$$

$$C_{21}=P_3+P_4=A_{21}\cdot B_{11}+A_{22}\cdot B_{21}$$

$$C_{22}=P_5+P_1-P_3-P_7=A_{22}\cdot B_{22}+A_{21}\cdot B_{12}$$

如此，我们便获得矩阵$A$和$B$的乘积——矩阵$C$。

算法分析

之前说过，Strassen算法的时间复杂度是优于朴素计算的，可是，它到底是多少呢？我们不妨再回到Strassen算法的流程。当$n>1$时，步骤1、2和4共花费$\theta (n^2)$时间，步骤3要求7次$\frac{n}{2}\ast \frac{n}{2}$矩阵的乘法。因此，我们得到如下描述Strassen算法运行时间$T(n)$的递归式：

$$T(n)=\{\begin{array}{rlr}& \theta (1)& 若n=1\\ & 7T(n/2)+\theta (n^2)& 若n>1\end{array}$$

求解上式可得，$T(n)=\theta (n^{\lg 7})$。

算法实现

废话千句，不如代码两行，接下来直接上Strassen算法的实现。（注意，如果$n$不是2的幂，可以采取对原矩阵填充0的方式，使$n$扩展到2的幂）。

算法总结

Strassen算法发表于1969年，它的发表引起了很大的轰动。在此之前，很少人敢设想一个算法能渐近快于平凡算法SQUARE-MATRIX-MULTIPLY。矩阵乘法的上界自此被改进了。到目前为止，$n\ast n$矩阵相乘的渐近复杂性最优的算法是Coppersmith和Winograd提出的，运行时间是$O(n^{2.376})$。