《算法导论》第四章矩阵相乘

1.采用递归分治的方法进行矩阵相乘，代码如下：

def matrix_divide(mat):
    """将n*n分解为4个n/2*n/2的四个矩阵"""
    n = len(mat) / 2
    a11 = [[0 for i in range(n)]for j in range(n)]
    a12 = [[0 for i in range(n)] for j in range(n)]
    a21 = [[0 for i in range(n)] for j in range(n)]
    a22 = [[0 for i in range(n)] for j in range(n)]
    for i in range(n):
        for j in range(n):
            a11[i][j] = mat[i][j]
            a12[i][j] = mat[i][n+j]
            a21[i][j] = mat[n+i][j]
            a22[i][j] = mat[n+i][n+j]

    return a11, a12, a21, a22


def matrix_add(a, b):
    """将两个矩阵相加"""
    n = len(a)
    c = [[0 for i in range(n)]for j in range(n)]
    for i in range(n):
        for j in range(n):
            c[i][j] = a[i][j] + b[i][j]

    return c

def matrix_merge(a11, a12, a21, a22):
    """将子矩阵合并成一大矩阵"""
    n = len(a11)
    a = [[0 for i in range(2*n)]for j in range(2*n)]
    for i in range(2*n):
        for j in range(2 * n):
            if i <= n-1 and j <= n-1:
                a[i][j] = a11[i][j]
            elif i <= n-1 and j > n-1:
                a[i][j] = a12[i][j-n]
            elif i > n-1 and j <= n-1:
                a[i][j] = a21[i-n][j]
            else:
                a[i][j] = a22[i-n][j-n]
    return a
def matrix_multiply(a, b):
    """两个矩阵递归相乘"""
    n = len(a)
    c = [[0 for i in range(n)] for j in range(n)]
    if n == 1:
        c[0][0] = a[0][0] * b[0][0]
    else:
        a11, a12, a21, a22 = matrix_divide(a)
        b11, b12, b21, b22 = matrix_divide(b)
        c11, c12, c21, c22 = matrix_divide(c)
        c11 = matrix_add(matrix_multiply(a11, b11), matrix_multiply(a12, b21))
        c12 = matrix_add(matrix_multiply(a11, b12), matrix_multiply(a12, b22))
        c21 = matrix_add(matrix_multiply(a21, b11), matrix_multiply(a22, b21))
        c22 = matrix_add(matrix_multiply(a21, b12), matrix_multiply(a22, b22))
        c = matrix_merge(c11, c12, c21, c22)

    return c

if __name__ == "__main__":
    a = [[1, 1, 1, 1], [1, 1, 1, 1], [2, 2, 2, 2], [2, 2, 2, 2]]
    b = a 
    print matrix_multiply(a,b)

2.Strassen算法实现

def matrix_divide(mat):
    """将n*n分解为4个n/2*n/2的四个矩阵"""
    n = len(mat) / 2
    a11 = [[0 for i in range(n)]for j in range(n)]
    a12 = [[0 for i in range(n)] for j in range(n)]
    a21 = [[0 for i in range(n)] for j in range(n)]
    a22 = [[0 for i in range(n)] for j in range(n)]
    for i in range(n):
        for j in range(n):
            a11[i][j] = mat[i][j]
            a12[i][j] = mat[i][n+j]
            a21[i][j] = mat[n+i][j]
            a22[i][j] = mat[n+i][n+j]

    return a11, a12, a21, a22

def matrix_merge(a11, a12, a21, a22):
    """将子矩阵合并成一大矩阵"""
    n = len(a11)
    a = [[0 for i in range(2*n)]for j in range(2*n)]
    for i in range(2*n):
        for j in range(2 * n):
            if i <= n-1 and j <= n-1:
                a[i][j] = a11[i][j]
            elif i <= n-1 and j > n-1:
                a[i][j] = a12[i][j-n]
            elif i > n-1 and j <= n-1:
                a[i][j] = a21[i-n][j]
            else:
                a[i][j] = a22[i-n][j-n]

    return a


def matrix_add(a, b):
    """将两个矩阵相加"""
    n = len(a)
    c = [[0 for i in range(n)]for j in range(n)]
    for i in range(n):
        for j in range(n):
            c[i][j] = a[i][j] + b[i][j]

    return c

def matrix_minus(a, b):
    """将两个矩阵相加"""
    n = len(a)
    c = [[0 for i in range(n)]for j in range(n)]
    for i in range(n):
        for j in range(n):
            c[i][j] = a[i][j] - b[i][j]

    return c

def strassen_multiply(a, b):
    n = len(a)
    c = [[0 for i in range(n)]for j in range(n)]

    if n == 1:
        c[0][0] = a[0][0] * b[0][0]
    else:
        a11, a12, a21, a22 = matrix_divide(a)
        b11, b12, b21, b22 = matrix_divide(b)
        c11, c12, c21, c22 = matrix_divide(c)

        s1 = matrix_minus(b12, b22)
        s2 = matrix_add(a11, a12)
        s3 = matrix_add(a21, a22)
        s4 = matrix_minus(b21, b11)
        s5 = matrix_add(a11, a22)
        s6 = matrix_add(b11, b22)
        s7 = matrix_minus(a12, a22)
        s8 = matrix_add(b21, b22)
        s9 = matrix_minus(a11, a21)
        s10 = matrix_add(b11, b12)

        p1 = strassen_multiply(a11, s1)
        p2 = strassen_multiply(s2, b22)
        p3 = strassen_multiply(s3, b11)
        p4 = strassen_multiply(a22, s4)
        p5 = strassen_multiply(s5, s6)
        p6 = strassen_multiply(s7, s8)
        p7 = strassen_multiply(s9, s10)

        c11 = matrix_add(matrix_minus(matrix_add(p5, p4), p2), p6)
        c12 = matrix_add(p1, p2)
        c21 = matrix_add(p3, p4)
        c22 = matrix_minus(matrix_minus(matrix_add(p5, p1), p3), p7)
        c = matrix_merge(c11, c12, c21, c22)
    return c

if __name__ == "__main__":
    a = [[1, 1, 1, 1], [1, 1, 1, 1], [2, 2, 2, 2], [2, 2, 2, 2]]
    b = a
    print(strassen_multiply(a, b))

3.前两个都需要矩阵是2的n次方，且必须是方阵，如果要求任意矩阵，可以考虑先用零填充到2的n次方，最后再删除到原大小。

# -*- coding:utf-8 -*-
#_author_ = 'xuqn'
from math import *

def matrix_divide(mat):
    """将n*n分解为4个n/2*n/2的四个矩阵"""
    n = len(mat) / 2
    a11 = [[0 for i in range(n)]for j in range(n)]
    a12 = [[0 for i in range(n)] for j in range(n)]
    a21 = [[0 for i in range(n)] for j in range(n)]
    a22 = [[0 for i in range(n)] for j in range(n)]
    for i in range(n):
        for j in range(n):
            a11[i][j] = mat[i][j]
            a12[i][j] = mat[i][n+j]
            a21[i][j] = mat[n+i][j]
            a22[i][j] = mat[n+i][n+j]

    return a11, a12, a21, a22


def matrix_add(a, b):
    """将两个矩阵相加"""
    n = len(a)
    c = [[0 for i in range(n)]for j in range(n)]
    for i in range(n):
        for j in range(n):
            c[i][j] = a[i][j] + b[i][j]

    return c

def matrix_merge(a11, a12, a21, a22):
    """将子矩阵合并成一大矩阵"""
    n = len(a11)
    a = [[0 for i in range(2*n)]for j in range(2*n)]
    for i in range(2*n):
        for j in range(2 * n):
            if i <= n-1 and j <= n-1:
                a[i][j] = a11[i][j]
            elif i <= n-1 and j > n-1:
                a[i][j] = a12[i][j-n]
            elif i > n-1 and j <= n-1:
                a[i][j] = a21[i-n][j]
            else:
                a[i][j] = a22[i-n][j-n]
    return a

def matrix_expand(mat):
    """将a扩展成2的n次方，0填充"""
    n = max(len(mat),len(mat[0]))
    m = ceil(log(n, 2))
    n_expand = int(2 ** m)
    c = [[0 for i in range(n_expand)] for j in range(n_expand)]
    for i in range(n_expand):
        for j in range(n_expand):
            if i < len(mat) and j < len(mat[0]):
                c[i][j] = mat[i][j]
            else:
                c[i][j] = 0
    return c

def matrix_shrink(mat, high, weight):
    """将a中之前填充的0去掉，输入需要降到的维度"""
    c = [[0 for i in range(weight)] for j in range(high)]
    for i in range(high):
        for j in range(weight):
            c[i][j] = mat[i][j]
    return c


def matrix_multiply(a, b):
    """两个矩阵递归相乘"""
    n = len(a)
    c = [[0 for i in range(n)] for j in range(n)]
    if n == 1:
        c[0][0] = a[0][0] * b[0][0]
    else:
        a11, a12, a21, a22 = matrix_divide(a)
        b11, b12, b21, b22 = matrix_divide(b)
        c11, c12, c21, c22 = matrix_divide(c)
        c11 = matrix_add(matrix_multiply(a11, b11), matrix_multiply(a12, b21))
        c12 = matrix_add(matrix_multiply(a11, b12), matrix_multiply(a12, b22))
        c21 = matrix_add(matrix_multiply(a21, b11), matrix_multiply(a22, b21))
        c22 = matrix_add(matrix_multiply(a21, b12), matrix_multiply(a22, b22))
        c = matrix_merge(c11, c12, c21, c22)

    return c

if __name__ == "__main__":
    a = [[1, 1, 1, 1], [1, 1, 1, 1], [2, 2, 2, 2]]
    b = [[1, 1], [2, 2], [2, 1], [1, 1]]
    print(matrix_shrink(matrix_multiply(matrix_expand(a), matrix_expand(b)), len(a), len(b[0])))

编程时看了网友的程序，在此表示感谢：

https://blog.csdn.net/weixin_42206504/article/details/81042798