【分治算法】Strassen矩阵乘法Python实现

文章目录

@[toc]
问题描述
基础算法
时间复杂性

`Strassen`算法
时间复杂性

问题时间复杂性
`Python`实现

因上努力

个人主页：丷从心.

系列专栏：Python基础

学习指南：Python学习指南

果上随缘

问题描述

设 $A$ 和 $B$ 是两个 $\times n$ 矩阵， $A$ 和 $B$ 的乘积矩阵 $C$ 中元素 $c_{ij} = \displaystyle\sum\limits_{k = 1}^{n}{a_{ik} b_{kj}}$
每计算 $C$ 的一个元素 $c_{ij}$ ，需要做 $n$ 次乘法和 $n - 1$ 次加法，求出矩阵 $C$ 的 $n^{2}$ 个元素所需的时间为 $O(n^{3})$

基础算法

假设 $n$ 是 $2$ 的幂，将矩阵 $A$ 、 $B$ 和 $C$ 中每个矩阵都分块成 $4$ 个大小相等的子矩阵，每个子矩阵都是 $n/2 \times n / 2$ 的方阵

$\begin{vmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{vmatrix} = \begin{vmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{vmatrix} \begin{vmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{vmatrix}$

$C_{11} = A_{11} B_{11} + A_{12} B_{21} \\ C_{12} = A_{11} B_{12} + A_{12} B_{22} \\ C_{21} = A_{21} B_{11} + A_{22} B_{21} \\ C_{22} = A_{21} B_{12} + A_{22} B_{22}$

时间复杂性

计算 $2$ 个 $n$ 阶方阵的乘积转化为计算 $8$ 个 $n /2$ 阶方阵的乘积和 $4$ 个 $n /2$ 阶方阵的加法， $2$ 个 $n/2 \times n / 2$ 矩阵的加法显然可以在 $O(n^{2})$ 时间内完成

$\begin{cases} O(1) & n = 2 \\ 8 T(n / 2) + O(n^{2}) & n > 2 \end{cases}$

$T(n) = O(n^{3})$

`Strassen`算法

Strassen算法只用了 $7$ 次乘法运算，但增加了加减法的运算次数

$M_{1} = A_{11} (B_{12} - B_{22}) \\ M_{2} = (A_{11} + A_{12}) B_{22} \\ M_{3} = (A_{21} + A_{22}) B_{11} \\ M_{4} = A_{22} (B_{21} - B_{11}) \\ M_{5} = (A_{11} + A_{22})(B_{11} + B_{22}) \\ M_{6} = (A_{12} - A_{22})(B_{21} + B_{22}) \\ M_{7} = (A_{11} - A_{21})(B_{11} + B_{12})$

$C_{11} = M_{5} + M_{4} - M_{2} + M_{6} \\ C_{12} = M_{1} + M_{2} \\ C_{21} = M_{3} + M_{4} \\ C_{22} = M_{5} + M_{1} - M_{3} - M_{7}$

时间复杂性

Strassen算法用了 $7$ 次对于 $n /2$ 阶矩阵乘积的递归调用和 $18$ 次 $n /2$ 阶矩阵的加减运算

$\begin{cases} O(1) & n = 2 \\ 7 T(n / 2) + O(n^{2}) & n > 2 \end{cases}$

$O(n^{\log{7}}) \approx O(n^{2.81})$

问题时间复杂性

$Ho p cro f t$ 和 $Kerr$ 已经证明计算 $2$ 个 $\times 2$ 矩阵的乘积， $7$ 次乘法是必要的
目前最好的计算时间上界是 $O(n^{2.376})$ ，所知的矩阵乘法的最好下界仍是它的平凡下界 $\Omega(n^{2})$

`Python`实现

import numpy as npdef strassen_matrix_multiply(a, b):n = a.shape[0]# 如果输入矩阵的维度小于等于阈值, 使用传统的矩阵乘法if n <= 128:return np.dot(a, b)# 将输入矩阵划分为四个子矩阵mid = n // 2a11 = a[:mid, :mid]a12 = a[:mid, mid:]a21 = a[mid:, :mid]a22 = a[mid:, mid:]b11 = b[:mid, :mid]b12 = b[:mid, mid:]b21 = b[mid:, :mid]b22 = b[mid:, mid:]# 递归地计算七个矩阵乘法m1 = strassen_matrix_multiply(a11, b12 - b22)m2 = strassen_matrix_multiply(a11 + a12, b22)m3 = strassen_matrix_multiply(a21 + a22, b11)m4 = strassen_matrix_multiply(a22, b21 - b11)m5 = strassen_matrix_multiply(a11 + a22, b11 + b22)m6 = strassen_matrix_multiply(a12 - a22, b21 + b22)m7 = strassen_matrix_multiply(a11 - a21, b11 + b12)# 计算结果矩阵的四个子矩阵c11 = m5 + m4 - m2 + m6c12 = m1 + m2c21 = m3 + m4c22 = m5 + m1 - m3 - m7# 组合四个子矩阵形成结果矩阵c = np.zeros((n, n))c[:mid, :mid] = c11c[:mid, mid:] = c12c[mid:, :mid] = c21c[mid:, mid:] = c22return ca = np.random.randint(0, 10, (256, 256))
b = np.random.randint(0, 10, (256, 256))res = strassen_matrix_multiply(a, b)print('矩阵 a:')
print(a)print('\n矩阵 b:')
print(b)print('\n乘积矩阵 c:')
print(res)

矩阵 a:
[[2 6 1 ... 9 7 7][7 1 8 ... 1 0 9][5 0 8 ... 6 2 5]...[6 7 2 ... 6 0 1][1 4 7 ... 0 5 2][3 3 3 ... 9 6 8]]矩阵 b:
[[3 5 1 ... 9 8 9][9 6 0 ... 8 1 5][3 9 6 ... 9 0 5]...[5 7 7 ... 3 8 8][9 5 3 ... 1 8 1][3 3 9 ... 4 7 0]]乘积矩阵 c:
[[5090. 5517. 4863. ... 4977. 4769. 4939.][5148. 5909. 5747. ... 5603. 5070. 5260.][4376. 5175. 4717. ... 4968. 4668. 4526.]...[4708. 5294. 4991. ... 4945. 4681. 5202.][4641. 5307. 4955. ... 5087. 5157. 4795.][4722. 5173. 5144. ... 5050. 4845. 4829.]]