1 朴素算法

这个算法就是矩阵乘法的定义：

• 很容易看出这个算法复杂度是$\Theta (n^3)$。

2 递归算法

分治法首先是从分割问题开始的，得到数学上的递归关系后，然后使用递归的方式实现。

由上面的数学性质，可以使用递归实现：

• $T(n)=8T(n/2)+\Theta (n^2)$，应用主定理可知，$T(n)=\Theta (n^3)$，没有什么进步，不过这里分割思路是没有问题的，只不过还有一步技巧性较强的变换在里面。

3 Strassen’s algorithm for matrix multiplication

1969年，斯特拉森（V.Strassen）利用分治策略并加上数学处理设计出了一种时间复杂度是$O(n^{2.81})$（准确地说是$O(n^{log7})$）的矩阵相乘算法，宣称在时间复杂度数量级上有所突破。此结果一发布，立即震动了整个数学界。

整个算法的核心是在分块的基础上做了一步变换，这个变换甚至高中生都可以看懂，但是只有牛人想的到。

得到上面的问题分割方式后，仿照递归算法容易实现。

• $T(n)=7T(n/2)+\Theta (n^2)$，应用主定理可知，$T(n)=\Theta (n^{log7})$。

4 实现

为便于说明问题这里仅考虑矩阵阶数为$2^n$的情况。

# -* coding: utf-8 -*-import numpy as np
from timeit import default_timer as timerdef STRASSEN_MATRIX(A, B):n = A.shape[0]C = np.array([0] * (n * n)).reshape((n, n))d = int(n / 2)if n == 1:C[0, 0] = A[0, 0] * B[0, 0]else:S1 = B[0:d, d:n] - B[d:n, d:n]S2 = A[0:d, 0:d] + A[0:d, d:n]S3 = A[d:n, 0:d] + A[d:n, d:n]S4 = B[d:n, 0:d] - B[0:d, 0:d]S5 = A[0:d, 0:d] + A[d:n, d:n]S6 = B[0:d, 0:d] + B[d:n, d:n]S7 = A[0:d, d:n] - A[d:n, d:n]S8 = B[d:n, 0:d] + B[d:n, d:n]S9 = A[0:d, 0:d] - A[d:n, 0:d]S10 = B[0:d, 0:d] + B[0:d, d:n]P1 = STRASSEN_MATRIX(A[0:d, 0:d], S1)P2 = STRASSEN_MATRIX(S2, B[d:n, d:n])P3 = STRASSEN_MATRIX(S3, B[0:d, 0:d])P4 = STRASSEN_MATRIX(A[d:n, d:n], S4)P5 = STRASSEN_MATRIX(S5, S6)P6 = STRASSEN_MATRIX(S7, S8)P7 = STRASSEN_MATRIX(S9, S10)C[0:d, 0:d] = P5 + P4 - P2 + P6C[0:d, d:n] = P1 + P2C[d:n, 0:d] = P3 + P4C[d:n, d:n] = P5 + P1 - P3 - P7return Cdef run_time():tic = timer()# 待测试的代码toc = timer()print(toc - tic)  # 输出的时间，秒为单位if __name__ == '__main__':A = np.array(list(range(0, 16))).astype(int).reshape((4, 4))B = np.array(list(range(2, 18))).astype(int).reshape((4, 4))print('strassen:\n', STRASSEN_MATRIX(A, B))print('dot:\n', np.dot(A, B))

运行结果：

strassen:[[ 68  74  80  86][196 218 240 262][324 362 400 438][452 506 560 614]]
dot:[[ 68  74  80  86][196 218 240 262][324 362 400 438][452 506 560 614]]