Python环境下存储稀疏矩阵LU分解结果的实现方法咨询

scipy的scipy.sparse.linalg.splu返回的SuperLU对象无法直接序列化，是因为它封装了底层C实现的结构体，Python的pickle无法直接处理C层面的数据，但该对象实际上暴露了所有重构LU分解所需的Python侧属性，不需要额外重复造轮子。以下是两种符合Python风格的实现方案：

方案1：直接提取SuperLU属性重构对象（推荐） #

该方案完全复用scipy原生的solve逻辑，性能和原生分解对象完全一致，代码量最小。

存储步骤 #

import pickle
import numpy as np
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import splu, SuperLU

# 示例稀疏矩阵，建议优先用CSC格式获得最佳性能
A = csc_matrix([[1, 2, 0], [0, 3, 4], [5, 0, 6]])
lu_obj = splu(A)

# 提取所有核心可序列化属性
lu_params = {
    "L": lu_obj.L,
    "U": lu_obj.U,
    "perm_c": lu_obj.perm_c,
    "perm_r": lu_obj.perm_r,
    "shape": lu_obj.shape
}

# 存储参数，也可以用numpy.savez、scipy.sparse.save_npz获得更好的跨版本兼容性
with open("lu_decomp.pkl", "wb") as f:
    pickle.dump(lu_params, f)

加载使用步骤 #

# 加载存储的参数
with open("lu_decomp.pkl", "rb") as f:
    lu_params = pickle.load(f)

# 重构SuperLU对象
restored_lu = SuperLU(
    shape=lu_params["shape"],
    L=lu_params["L"],
    U=lu_params["U"],
    perm_c=lu_params["perm_c"],
    perm_r=lu_params["perm_r"]
)

# 和原生对象完全一致调用solve方法
b = np.array([1, 2, 3])
x = restored_lu.solve(b)

方案2：手动管理分解组件（可控性最高） #

如果担心SuperLU对象的跨版本兼容性，可以直接用scipy提供的lu函数获取置换矩阵、L、U矩阵，自行实现求解逻辑，所有组件都是标准的稀疏矩阵/NumPy数组，天然支持序列化。

from scipy.sparse.linalg import lu, spsolve_triangular

# 直接获得分解组件
P, L, U = lu(A)

# 三个矩阵都支持pickle、save_npz等任意序列化方式
# 求解逻辑自行实现
def custom_lu_solve(P, L, U, b):
    permuted_b = P @ b
    y = spsolve_triangular(L, permuted_b, lower=True)
    x = spsolve_triangular(U, y, lower=False)
    return x

# 调用求解
x = custom_lu_solve(P, L, U, b)

注意事项 #

• 稀疏LU分解默认对输入矩阵做列置换和行置换来提升数值稳定性，存储时必须同时保存置换参数，否则求解结果会出错
• 跨scipy大版本迁移时，优先使用方案2，避免SuperLU构造参数变更导致的加载失败
• 如果矩阵是对称正定结构，优先选择Cholesky分解，存储成本更低、求解速度更快

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