为何Python实现的Strassen算法比标准矩阵乘法慢？

我需要在不依赖NumPy的Python环境中实现高效矩阵乘法，因此实现了Strassen算法，但测试显示它比循环实现的标准矩阵乘法慢3.5倍。测试仅针对N为2的幂的N×N方阵，使用from time import perf_counter()计时，仅统计函数执行耗时，排除矩阵创建和结果输出时间，矩阵元素由randint(1, 9)生成。测试数据如下：

希望优化现有Strassen算法实现，而非完全重写。

现有实现代码 #

Strassen算法实现 #

from time import perf_counter

t_start = perf_counter()

def submatrices(n, matrix): # 将矩阵划分为子块
    A = [[j for j in matrix[i][:int(n / 2)]] for i in range(int(n / 2))]
    B = [[j for j in matrix[i][int(n / 2):]] for i in range(int(n / 2))]
    C = [[j for j in matrix[i][:int(n / 2)]] for i in range(int(n / 2), n)]
    D = [[j for j in matrix[i][int(n / 2):]] for i in range(int(n / 2), n)]
    return [A, B, C, D]


def addition(n, matrix1, matrix2):  # 矩阵加法
    res = [[matrix1[i][j] + matrix2[i][j] for j in range(n)] for i in range(n)]
    return res


def subtraction(n, matrix1, matrix2):   # 矩阵减法
    res = [[matrix1[i][j] - matrix2[i][j] for j in range(n)] for i in range(n)]
    return res


def strassen(n, matrix1, matrix2):

    if n == 2:  # 算法最后一步采用标准矩阵乘法
        xy = [[0] * n for i in range(n)]
        for i in range(n):
            for j in range(n):
                for x in range(n):
                    xy[i][j] += matrix1[i][x] * matrix2[x][j]
    else:
        A, B, C, D = submatrices(n, matrix1)    # 划分原始矩阵1
        E, F, G, H = submatrices(n, matrix2)    # 划分原始矩阵2

        n = int(n / 2)  # 此时矩阵尺寸减半

        p1 = strassen(n, A, subtraction(n, F, H))
        p2 = strassen(n, addition(n, A, B), H)
        p3 = strassen(n, addition(n, C, D), E)
        p4 = strassen(n, D, subtraction(n, G, E))
        p5 = strassen(n, addition(n, A, D), addition(n, E, H))
        p6 = strassen(n, subtraction(n, B, D), addition(n, G, H))
        p7 = strassen(n, subtraction(n, A, C), addition(n, E, F))

        xy1 = addition(n, addition(n, p5, p6), subtraction(n, p4, p2))  # 构建结果矩阵子块
        xy2 = addition(n, p1, p2)
        xy3 = addition(n, p3, p4)
        xy4 = subtraction(n, addition(n, p1, p5), addition(n, p3, p7))

        xy = [xy1[i] + xy2[i] for i in range(n)] + [xy3[i] + xy4[i] for i in range(n)]  # 拼接子块得到结果矩阵
    return xy

print(f'Time: {perf_counter() - t_start} sec')

# 打印结果
# for raw in strassen(n, matrix1, matrix2):
#     print(*raw)

标准矩阵乘法实现 #

from time import perf_counter

t_start = perf_counter()


def multiply(n, matrix1, matrix2):
    res = [[0]*n for i in range(n)]
    for i in range(n):
        for j in range(n):
            for x in range(n):
                res[i][j] += matrix1[i][x] * matrix2[x][j]

    return res

print(f'Time: {perf_counter() - t_start} sec')

# 打印结果
# for raw in multiply(n, matrix1, matrix2):
#     print(*raw)

优化方案 #

1. 消除子矩阵复制开销 #

当前submatrices函数会创建全新的子矩阵列表，涉及大量数据复制，这是主要性能瓶颈之一。改为通过索引范围描述子矩阵，而非复制数据：

def get_submatrix_bounds(n):
    half = n // 2
    # 返回四个子矩阵的(行起始, 行结束, 列起始, 列结束)
    return (
        (0, half, 0, half),    # A
        (0, half, half, n),    # B
        (half, n, 0, half),    # C
        (half, n, half, n)     # D
    )

# 调整加法函数，支持按索引范围操作原矩阵
def add_submatrices(mat1, mat2, bounds1, bounds2, result_bounds, result_mat):
    r1_start, r1_end, c1_start, c1_end = bounds1
    r2_start, r2_end, c2_start, c2_end = bounds2
    res_r_start, res_r_end, res_c_start, res_c_end = result_bounds
    
    for i in range(res_r_end - res_r_start):
        for j in range(res_c_end - res_c_start):
            result_mat[res_r_start + i][res_c_start + j] = mat1[r1_start + i][c1_start + j] + mat2[r2_start + i][c2_start + j]

# 减法函数同理
def sub_submatrices(mat1, mat2, bounds1, bounds2, result_bounds, result_mat):
    r1_start, r1_end, c1_start, c1_end = bounds1
    r2_start, r2_end, c2_start, c2_end = bounds2
    res_r_start, res_r_end, res_c_start, res_c_end = result_bounds
    
    for i in range(res_r_end - res_r_start):
        for j in range(res_c_end - res_c_start):
            result_mat[res_r_start + i][res_c_start + j] = mat1[r1_start + i][c1_start + j] - mat2[r2_start + i][c2_start + j]

递归过程中不再复制子矩阵，直接操作原矩阵的指定区域，大幅减少内存操作开销。

2. 调整递归终止阈值 #

当前递归终止于n=2，但Python函数调用有固定开销，当矩阵尺寸较小时，标准乘法的实际速度更快。建议将终止阈值调整为n ≤ 64（可根据实际测试调整），此时直接切换为标准乘法：

def strassen(n, matrix1, matrix2):
    # 调整终止阈值
    if n <= 64:
        return multiply(n, matrix1, matrix2)
    # ... 剩余递归逻辑

3. 优化矩阵拼接与临时对象创建 #

原代码中通过列表相加拼接结果矩阵，且多次创建临时矩阵。改为预先分配结果矩阵，直接填充数据：

# 替换原拼接逻辑
half_n = n // 2
xy = [[0]*n for _ in range(n)]
# 填充xy1和xy2区域
for i in range(half_n):
    xy[i][:half_n] = xy1[i]
    xy[i][half_n:] = xy2[i]
# 填充xy3和xy4区域
for i in range(half_n):
    xy[half_n + i][:half_n] = xy3[i]
    xy[half_n + i][half_n:] = xy4[i]

4. 减少全局变量查找开销 #

在递归函数内部，将全局函数（如addition、subtraction）赋值为局部变量，Python局部变量查找速度远快于全局：

def strassen(n, matrix1, matrix2):
    # 局部变量缓存全局函数
    add = addition
    sub = subtraction
    # ... 后续递归中使用add和sub替代原函数名

5. 避免重复计算n/2 #

在递归分支中，多次计算int(n/2)，改为提前计算一次：

else:
    half_n = n // 2
    A, B, C, D = submatrices(n, matrix1)
    E, F, G, H = submatrices(n, matrix2)
    # 后续用half_n替代int(n/2)

内容的提问来源于stack exchange，提问作者kdduha