如何实现Scipy Linear Operator与CSC稀疏矩阵的矩阵乘积？

Scipy并没有内置直接计算LinearOperator和CSC稀疏矩阵乘积并返回稀疏矩阵的方法，但可以利用CSC矩阵的列结构手动实现——因为CSC矩阵的每一列本质是稀疏向量，我们可以对每一列调用LinearOperator的matvec方法，再将结果组合为CSC矩阵。

实现代码 #

import numpy as np
from scipy.sparse import csc_matrix, hstack
from scipy.sparse.linalg import LinearOperator

def csc_linop_mult(A: csc_matrix, lin_op: LinearOperator) -> csc_matrix:
    # 维度合法性校验
    if A.shape[1] != lin_op.shape[0]:
        raise ValueError("维度不匹配：A的列数必须等于LinearOperator的行数")
    
    # 遍历A的每一列，计算与LinearOperator的乘积
    sparse_cols = []
    for col_idx in range(A.shape[1]):
        # 提取当前列并转为一维数组
        current_col = A[:, col_idx].toarray().ravel()
        # 调用LinearOperator的matvec方法计算乘积
        product_col = lin_op.matvec(current_col)
        # 转为稀疏列矩阵后收集
        sparse_cols.append(csc_matrix(product_col.reshape(-1, 1)))
    
    # 堆叠所有稀疏列得到最终CSC矩阵
    return hstack(sparse_cols, format='csc')

测试你的示例 #

# 定义CSC矩阵A
A = csc_matrix([[1, 0], [0, 1]])

# 基于稠密矩阵创建LinearOperator（模拟你的H²矩阵场景）
B = np.array([[10, 0], [0, 10]] )
lin_op_B = LinearOperator(B.shape, matvec=B.dot)

# 计算乘积
result = csc_linop_mult(A, lin_op_B)
print(result.todense())
# 输出：
# [[10  0]
#  [ 0 10]]
print(type(result))  # 确认结果为CSC格式：<class 'scipy.sparse._csc.csc_matrix'>

优化说明 #

如果你的CSC矩阵列数极多，循环单列的方式可能效率偏低，可以利用CSC矩阵的indptr和indices属性批量处理非零元素，减少循环次数。但对于大多数H²矩阵的使用场景，上述简单实现的性能已经足够匹配h2tools的高效matvec操作。

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