矩阵乘法的Strassen算法

最新推荐文章于 2022-09-01 20:14:50 发布

原创最新推荐文章于 2022-09-01 20:14:50 发布 · 2.6k 阅读

4 ·

CC 4.0 BY-SA版权

文章标签：

#矩阵算法

算法同时被 2 个专栏收录

25 篇文章

订阅专栏

算法导论

17 篇文章

订阅专栏

本文介绍了矩阵乘法的基本算法及其实现，并提出了一种基于分治的优化算法。进一步介绍了Strassen算法，通过减少递归调用次数降低时间复杂度。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

若 $A=(a_{ij})$ 和 $B=(b_{ij})$ 是 $n*n$ 方阵,则对于i,j=1,2,…,n；定义乘积C = A * B中的元素 $c_{ij}$ 为：

c i j = \sum k = 1 n a i k b k j

$c_{ij}= \sum_{k=1}^{n} a_{ik}b_{kj}$
我们需要计算

n2 $n^2$ 个矩阵元素，每个元素是n个值的和。下面过程接收n*n矩阵A和B，返回它们的乘积 n* n矩阵C.假设每个矩阵都有一个属性rows，给出矩阵的行数。
SQUARE-MATRIX-MULTIPLY(A, B)

n = A.rows
let C be a new n*n matrix
for i = 1 to n
    for j = 1 to n
        c(ij) = 0
        for k = 1 to n
            c(ij) = c(ij)+a(ik)*b(kj)
return C

C语言代码实现如下：

void Mul(int **matrixA,int **matrixB, int **matrixC, int n)
{
    for (int i=0; i < n; ++i){
        for (int j=0; j < n; ++j){
            matrixC[i][j] = 0;
            for (int k=0; k < n; ++k)
                matrixC[i][j] += matrixA[i][k] + matrixB[k][j];
        }
    }
}

由于三重for循环的每一重都恰好执行n步，而第7步每次执行都花费常量时间，因此过程SQUARE-MATRIX-MULTIPLY(A, B)花费 $\varTheta(n^3)$

试着用一个简单的分治算法优化刚刚的矩阵算法， C = A * B,假定三个矩阵均为 n * n 矩阵，其中 n 为 2 的幂。我们做这个假定因为在每个分解步骤中， n * n 矩阵都被划分为 4 个 n/2 * n/2 的子矩阵，如果假定 n 是 2 的幂，则只要 $n \geq 2$ 即可保证子矩阵规模 n/2 为整数。
假定将 A、B 和 C 均分解为 4 个 n/2 * n/2 的子矩阵:

A = $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}$ , B = $\begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}$ , C = $\begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix}$ （公式4.9）

因此可以将公式C = A * B 改写为

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

等价于如下 4 个公式：

$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}$

直接的递归分治算法：

SQUARE-MATRIX-MULTIPLY-RECURSIVE(A, B)

n = A.rows
let C be a new n * n matrix
if n == 1
    c11 = a11 * b11
else partion A, B, and C as in equations (4, 9)
    C11 = SQUARE-MATRIX-MULTIPLY-RECURSIVE(A11, B11)
    + SQUARE-MATRIX-MULTIPLY-RECURSIVE(A12, B21)
    C12 = SQUARE-MATRIX-MULTIPLY-RECURSIVE(A11, B12)
    + SQUARE-MATRIX-MULTIPLY-RECURSIVE(A12, B22)
    C21 = SQUARE-MATRIX-MULTIPLY-RECURSIVE(A21, B11)
    + SQUARE-MATRIX-MULTIPLY-RECURSIVE(A22, B21)
    C22 = SQUARE-MATRIX-MULTIPLY-RECURSIVE(A21, B12)
    + SQUARE-MATRIX-MULTIPLY-RECURSIVE(A22, B22)
return C

由于每次递归调用完成两个 n/2 * n/2 矩阵的乘法，因此花费时间为T(n/2), 8 次递归调用总时间为8T(n/2). 我们还需要计算第6-9行的4次矩阵加法。每个矩阵包含 $n^2/4$ 个元素，因此，每次矩阵加法花费 $\varTheta(n^2)$ 时间其他的时间开销为 $\varTheta(1)$ 故矩阵加法时间之和：
$T(n) = \varTheta(1)+ 8T(n/2) + \varTheta(n^2) = 8T(n/2) + \varTheta(n^2)$

T (n) = {Θ (1) 8 T (n / 2) + Θ (n 2) = = n = 1 n > 1

$T(n) =\left\{ \begin{aligned} \varTheta(1) &=& n=1\\ 8T(n/2) + \varTheta(n^2) & = & n>1 \end{aligned} \right.$
利用主方法求解上式得：

T(n)=Θ(n3) $T(n) = \varTheta(n^3)$
(

T(n)=aT(n/b)+f(n)这里的情况符合Θ(nlogba)>f(n),Θ(nlog28)=Θ(n3)>Θ(n2) $T(n) = aT(n/b) + f(n) 这里的情况符合 \varTheta(n^{log_ba} ) > f(n) , \varTheta(n^{log_28} ) = \varTheta(n^3 ) > \varTheta(n^2)$ )
可见用了一般的分治算法并没有使算法的 T(n) 减少

所以才有了天才般的Strassen方法，Strassen算法的核心思想是令递归树稍微不那么茂盛一点儿，即只进行7次而不是8次递n/2 * n/2的矩阵乘法。它包含4个步骤：
1. 按公式（4.9）将输入矩阵A、B和输出矩阵C分解为n/2 * n/2 的子矩阵。采用下标计算方法，此步骤花费 $\varTheta(1)$ 时间，与SQUARE-MATRIX-MULTIPLY-RECURSIVE相同
2. 创建10个n/2 * n/2的矩阵 $S_1, S_2, ..., S_{10}$ 每个矩阵保存步骤 1 中创建的两个矩阵的和或差。花费时间为 $\varTheta(n^2)$
3. 用步骤1中创建的子矩阵和步骤2中创建的10个矩阵，递归的计算7个矩阵积 $P_1,P_2,..P_7$ 。每个矩阵 $P_i$ 都是n/2 * n/2 的。
4. 通过 $P_i$ 矩阵的不同组合进行加减运算，计算出结果矩阵C的子矩阵 $C_{11},C_{12}, C_{21}, C_{22}$ .花费时间 $\varTheta(n^2)$ .