如何使用Scipy/Python将稀疏线性方程组简化为等价[A]'[x]=[B]'形式

实现思路说明 #

你需要的转换本质是对原增广矩阵做行等价变换，最终得到对角形式的系数矩阵A'，保证变换前后方程组解空间完全一致。由于原矩阵是高稀疏度的整数方阵，我们分规模给出对应的Python实现方案：

小规模矩阵实现（万阶以下） #

直接用Sympy的稀疏矩阵RREF（简化行阶梯形）方法完成变换，符号运算可保证无精度损失：

import numpy as np
from sympy import SparseMatrix

# 构造示例稀疏矩阵A和向量B，可替换为你自己的输入
n = 1000
A = np.zeros((n, n), dtype=np.int8)
# 模拟每行6个非零元素
for i in range(n):
    select_cols = np.random.choice(n, 6, replace=False)
    A[i, select_cols] = np.random.randint(-128, 127, size=6, dtype=np.int8)
B = np.random.randint(-128, 127, size=(n, 1), dtype=np.int8)

# 构造增广矩阵并计算RREF
aug_mat = SparseMatrix(np.hstack([A, B]))
rref_result, pivot_cols = aug_mat.rref()

# 提取变换后的A'和B'
A_prime = rref_result[:, :-1]
B_prime = rref_result[:, -1]

# 如需A'严格每行每列仅一个非零，对列做置换即可
perm_cols = list(pivot_cols) + [col for col in range(n) if col not in pivot_cols]
A_prime_permuted = A_prime[:, perm_cols]

该方案精度高但运算速度较慢，仅适合中小规模矩阵。

大规模稀疏矩阵实现（十万阶及以上） #

用Scipy的稀疏LU分解配合消元实现，全程用稀疏矩阵运算，内存占用低、速度快：

import numpy as np
from scipy.sparse import csr_matrix, hstack
from scipy.sparse.linalg import splu
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

# 构造示例稀疏矩阵A和向量B，可替换为你自己的输入
n = 100000
rows, cols, data = [], [], []
for i in range(n):
    select_cols = np.random.choice(n, 6, replace=False)
    select_data = np.random.randint(-128, 127, size=6, dtype=np.int8)
    rows.extend([i]*6)
    cols.extend(select_cols)
    data.extend(select_data)
A = csr_matrix((data, (rows, cols)), shape=(n, n), dtype=np.int8)
B = csr_matrix(np.random.randint(-128, 127, size=(n, 1), dtype=np.int8))

# LU分解得到上三角矩阵和置换信息
lu_decomp = splu(A.tocsc(), permc_spec="COLAMD")
U_mat = lu_decomp.U.tocsr()
row_perm = lu_decomp.perm_r
col_perm = lu_decomp.perm_c

# 对上三角矩阵做反向消元，得到对角形式A'，同步更新B'
for i in reversed(range(n)):
    pivot_val = U_mat[i, i]
    if pivot_val == 0:
        continue
    # 消去当前列所有上方非零元素
    cur_row_cols = U_mat.indices[U_mat.indptr[i]:U_mat.indptr[i+1]]
    for j in cur_row_cols[cur_row_cols < i]:
        factor = U_mat[j, i] // pivot_val
        U_mat[j, :] -= factor * U_mat[i, :]
        B[j] -= factor * B[row_perm[i]]

# 应用列置换得到最终的A'
A_prime = U_mat[:, col_perm]
B_prime = B

注意事项 #

• 若原矩阵为奇异矩阵，消元后会出现全零行：若对应B'位置不为零则原方程组无解，全零行对应变量为自由变量，完全匹配原方程组的解空间规则
• 若要求A'的非零元素为1，只需每行除以主元值，同步修改B'对应位置即可
• 全程使用CSR/CSC格式稀疏矩阵，不要转为稠密数组，避免内存溢出

内容的提问来源于stack exchange，提问作者J.Doe