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

news2024/10/7 10:15:56

文章目录

@[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 np


def strassen_matrix_multiply(a, b):
    n = a.shape[0]

    # 如果输入矩阵的维度小于等于阈值, 使用传统的矩阵乘法
    if n <= 128:
        return np.dot(a, b)

    # 将输入矩阵划分为四个子矩阵
    mid = n // 2

    a11 = 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 + m6
    c12 = m1 + m2
    c21 = m3 + m4
    c22 = m5 + m1 - m3 - m7

    # 组合四个子矩阵形成结果矩阵
    c = np.zeros((n, n))

    c[:mid, :mid] = c11
    c[:mid, mid:] = c12
    c[mid:, :mid] = c21
    c[mid:, mid:] = c22

    return c


a = 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.]]