如何计算SciPy稀疏CSR矩阵的幂次并解决np.matmul运算报错问题

报错原因 #

numpy.matmul仅适配Numpy稠密数组，无法识别Scipy稀疏矩阵的结构，会将输入的稀疏矩阵判定为0维对象，触发维度不匹配的ValueError。

解决方案 #

方案1：直接使用稀疏矩阵原生幂运算（推荐） #

Scipy的稀疏矩阵对象原生支持**运算符完成幂运算，内部已经做了快速幂和稀疏计算优化，最适合10001×10001这类大尺寸稀疏矩阵场景，代码示例如下：

import scipy.sparse as sparse
import numpy as np

# 测试矩阵构造
A = [[1,1],[1,0]]
AS = sparse.csr_matrix(A)

# 直接计算10次幂
AS_10 = AS ** 10

# 小矩阵可以转稠密数组验证结果，大矩阵请勿执行该操作，避免内存溢出
print(AS_10.todense())

运行后输出结果符合预期，AS^2为[[2 1],[1 1]]，AS^10为[[89 55],[55 34]]。

方案2：修改自定义乘法函数 #

如果需要自行实现幂运算逻辑，将np.matmul替换为稀疏矩阵自带的.dot()方法即可完成矩阵乘法：

def matrixMul(AS, n):
    if n <= 1:
        return AS
    else:
        # 调用稀疏矩阵dot方法完成乘法
        return matrixMul(AS, n-1).dot(AS)

注意：该递归实现是O(n)时间复杂度，幂次较高时性能远低于官方原生的O(log n)快速幂实现，仅适合测试场景使用。

大稀疏矩阵使用注意事项 #

• 非必要不要调用todense()/toarray()方法将稀疏矩阵转为稠密矩阵，10001×10001的稠密矩阵会占用近800MB内存（单精度浮点），远高于稀疏矩阵的内存占用。
• 官方原生幂运算会自动保留矩阵的稀疏结构，仅计算非零元素相关运算，计算效率和内存占用都最优。

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