模2二进制矩阵乘法算法性能优化及相关技术咨询

已明确问题答案 #

• 在操作总数不变的前提下，用与/逻辑异或替代加法/乘法后取模2的操作，确实可节省部分时间，但该时间是否“显著”尚存争议。
• 若想了解numpy底层BLAS矩阵乘法的具体实现思路，无需解读古老FORTRAN代码，可参考提出当前BLAS矩阵乘法算法的原始论文。

初始问题原文 #

以下是两种模2二进制矩阵乘法（即计算“点积”）的算法。第一种（“默认”方法）直接使用numpy矩阵乘法，再取模2；第二种（“替代”方法）尝试用异或操作替代加法以提升速度。

import timeit
import numpy as np
import matplotlib.pyplot as plt

def mat_mult_1(A,B):
    return A@B%2

def mat_mult_2(A,B):
    return np.logical_xor.reduce(A[:,:,None]&B[None,:,:],axis = 1)

与预期相反，对于较大的二进制矩阵，替代方法的耗时约为默认方法的4倍。这是为何？是否有方法可加速替代算法？

以下是测试上述两种方法的脚本：

n_vals = np.arange(5,205,5)
times = []

for n in n_vals:
    s_1 = f"mat_mult_1(np.random.randint(2,size = ({n},{n}))\
                        ,np.random.randint(2,size = ({n},{n})))"
    s_2 = f"mat_mult_2(np.random.randint(2,size = ({n},{n})),\
                        np.random.randint(2,size = ({n},{n})))"

    times.append((timeit.timeit(s_1, globals = globals(), number = 100),
              timeit.timeit(s_2, globals = globals(), number = 100)))

测试结果的图表如下：

后续更新 #

1. 渐近复杂度测试
   测试了最大1000x1000的矩阵，明确了算法的渐近复杂度：默认算法复杂度为O(n².⁷)，替代方法符合预期的O(n³)（实际观测斜率为2.703和3.133）。
   

2. 与“教科书式”实现对比
   对比了替代算法与以下“教科书式”矩阵乘法后取模的实现：

def mat_mult_3(A,B):
    return np.sum(A[:,:,None]*B[None,:,:],axis = 1)%2

结果该实现性能优于基于与/异或的方法！

3. 布尔矩阵运算优化测试
   针对评论修改了mat_mult_2，使用布尔矩阵运算：

def mat_mult_2(A,B):
    return np.logical_xor.reduce(A.astype(bool)[:,:,None]  
            & B.astype(bool)[None,:,:],axis = 1).astype(int)

尽管存在类型转换开销，但性能未显著改变，修改后mat_mult_2性能略优于mat_mult_3，符合预期。

4. 默认方法复杂度复测
   针对评论再次测试@方法的渐近复杂度，因设备限制无法处理远大于2000的矩阵，观测到的斜率为2.93。
   

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