《算法导论》——Strassen算法

矩阵乘法

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

c_{i j} = \sum_{k = 1}^{n} a_{i k} b_{k j}

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

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

算法流程

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

将输入矩阵 $A 、 B$ 和输出矩阵 $C$ 通过以下方式分解为 $\frac{n}{2} * \frac{n}{2}$ 的子矩阵；
$A = [\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}], B = [\begin{matrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{matrix}], C = [\begin{matrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{matrix}]$
创建10个 $\frac{n}{2} * \frac{n}{2}$ 的矩阵 $S_{1}, S_{2}, \dots, S_{10}$ ，每个矩阵保存步骤1中创建的两个子矩阵的和或差，时间复杂度为 $Θ (n^{2})$ ；
用步骤1中创建的子矩阵和步骤2中创建的10个矩阵，递归地计算7个矩阵积 $P_{1}, P_{2}, \dots, P_{7}$ 。每个矩阵 $P_{i}$ 都是 $\frac{n}{2} * \frac{n}{2}$ 的；
通过 $P_{i}$ 矩阵的不同组合进行加减计算，计算出矩阵 $C$ 的子矩阵 $C_{11}, C_{12}, C_{21}, C_{22}$ ，时间复杂度为 $Θ (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} * \frac{n}{2}$ 的加减法，因此，该步骤花费 $Θ (n^{2})$ 。

在步骤三中，递归地计算7次 $\frac{n}{2} * \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} * \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共花费 $θ (n^{2})$ 时间，步骤3要求7次 $\frac{n}{2} * \frac{n}{2}$ 矩阵的乘法。因此，我们得到如下描述Strassen算法运行时间 $T (n)$ 的递归式：

T (n) = {\begin{aligned} θ (1) & 若 n = 1 \\ 7 T (n / 2) + θ (n^{2}) & 若 n > 1 \end{aligned}

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

算法实现

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

算法总结

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

矩阵乘法——Strassen算法