4. 稀疏矩阵与大规模线性系统

在科学计算中，很多矩阵的绝大多数元素都是零。例如有限差分网格上的拉普拉斯算子、图结构的邻接矩阵、推荐系统的用户—物品交互矩阵，非零元往往只占全部元素的千分之一甚至更少。若仍用 NumPy 的稠密二维数组存储，一个 10⁶×10⁶ 的矩阵即使每元素 8 字节也要 8 TB 内存，完全不可行。scipy.sparse 正是为此设计的稀疏矩阵工具集：它只保存非零元及其位置，并提供专门的矩阵—向量乘法和线性求解器，使问题规模从“内存装不下”变为“可以计算”。

三种核心存储格式：COO、CSR、CSC

稀疏矩阵的存储格式决定了它的构建效率、运算效率和可修改性。SciPy 中最常用的三种格式各有明确的适用场景。

COO（坐标格式）用三个等长数组描述矩阵：row、col 记录每个非零元的行、列下标，data 记录数值。它最直观，允许重复坐标（后续求和合并），构建时无需预知每行非零元个数，因此非常适合“边遍历边收集三元组”的场景，例如从有限元单元组装整体刚度矩阵。缺点是随机访问和算术运算效率低，通常作为中间格式，构建完成后转换为 CSR 或 CSC。

CSR（压缩稀疏行）按行存储：data 存放非零元值，indices 存放对应列号，indptr 长度为行数加一，indptr[i] 到 indptr[i+1] 之间的切片就是第 i 行的非零元。CSR 的行切片极快，矩阵—向量乘法 A @ x 按行累加，缓存友好，是大多数运算和求解器的首选格式。CSC（压缩稀疏列）是 CSR 的转置结构，按列压缩，适合按列访问的场景，例如稀疏 LU 分解中的列主元选择。二者可通过 .tocsr()、.tocsc() 互相转换，转换本身有 O(nnz) 开销，应避免在循环中反复转换。

稀疏矩阵的构建与运算

构建稀疏矩阵有两条主线。第一条是“三元组收集”：先用 Python 列表累积 rows、cols、vals，最后一次性调用 coo_matrix((vals, (rows, cols)), shape=(n, n))。这种方式在有限元组装中非常自然——遍历每个单元，把局部刚度矩阵的贡献按全局自由度编号写入三元组，重复坐标会在转 CSR 时自动求和。第二条是“按结构生成”：对于规则网格，可直接用 scipy.sparse.diags 构造对角线，例如一维拉普拉斯算子 diags([-1, 2, -1], [-1, 0, 1], shape=(n, n))，无需循环，效率极高。

稀疏矩阵支持与 NumPy 类似的算术：A + B、A * 2、A @ x、A.multiply(B)（逐元素乘）。需要特别注意 * 在稀疏矩阵上是矩阵乘法还是逐元素乘，取决于对象类型：scipy.sparse 的矩阵对象中 * 表示矩阵乘法，而 @ 始终表示矩阵乘法，逐元素乘应显式使用 .multiply()。此外，A.toarray() 会把稀疏矩阵转为稠密数组，仅在小规模调试时使用；A.nnz 返回非零元个数，是判断稀疏度、估算内存的第一手指标。

稀疏线性求解器：直接法与迭代法

求解 A x = b 时，SciPy 提供两类求解器。直接法以 spsolve 为代表，内部对矩阵做稀疏 LU 或 Cholesky 分解。它的优点是稳健、无需调参，对中等规模（例如十万未知量以内、二维结构问题）通常一次调用即可得到高精度解；缺点是分解过程会产生“填充”（fill-in），即原本为零的位置在消元中变为非零，内存和计算量可能急剧增长。对三维问题或非结构化网格，填充往往使直接法不可承受。

迭代法不分解矩阵，而是从初值出发逐步逼近解，每次迭代主要成本是一次稀疏矩阵—向量乘法。常用接口包括：cg（共轭梯度法，仅适用于对称正定矩阵）、gmres（广义最小残差法，适用于一般非对称矩阵）、bicgstab、minres 等。调用形式为 cg(A, b, rtol=1e-8, maxiter=1000)，返回 (x, info)：info=0 表示在容差内收敛，info>0 表示达到最大迭代次数仍未收敛，info<0 表示 breakdown。判断是否成功，不能只看函数是否返回，必须检查 info 并计算残差 norm(b - A @ x) / norm(b)。

预条件子：让迭代法真正可用

原始迭代法对条件数大的矩阵收敛极慢。预条件子的思想是把原方程等价变换为 M⁻¹A x = M⁻¹b，其中 M 近似 A 但易于求逆，使变换后矩阵的条件数显著降低。SciPy 在 scipy.sparse.linalg 中提供多种预条件子：jacobi（对角缩放，最便宜）、ilu（不完全 LU 分解，对非对称问题效果好）、spilu（更灵活的 ILU 实现）、LinearOperator 自定义预条件。使用方式是把预条件对象作为 M= 参数传入：

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg, spilu, LinearOperator

n = 200
A = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csr')
b = np.ones(n)

# 无预条件
x0, info0 = cg(A, b, rtol=1e-8, maxiter=500)

# ILU 预条件
ilu = spilu(A.tocsc(), drop_tol=1e-4, fill_factor=10)
M = LinearOperator(A.shape, ilu.solve)
x1, info1 = cg(A, b, rtol=1e-8, maxiter=500, M=M)

print(info0, info1)
print(np.linalg.norm(b - A @ x1) / np.linalg.norm(b))

这段示例中，drop_tol 控制 ILU 分解时丢弃小元素的阈值，fill_factor 限制填充增长的上限，二者共同决定预条件子的质量与构建成本。实际调参时，先观察无预条件时的迭代次数与 info，再加入 Jacobi 或 ILU 对比；若 info 仍为正，可适当放宽 drop_tol 或增大 fill_factor，同时监控预条件子构建时间是否已成为瓶颈。对于对称正定问题，优先尝试 cg 加不完全 Cholesky 类预条件；对于非对称问题，gmres 配 ILU 是常见组合。

总体而言，稀疏计算的关键不在于记住 API，而在于建立“格式—运算—求解器—预条件”的匹配意识：COO 负责灵活组装，CSR 负责高效运算，直接法负责稳健小中规模，迭代法加预条件负责大规模问题。理解这条链路，才能在面对百万级未知量时做出正确的工程选择。