10. 稀疏矩阵与定常迭代

大规模离散问题中的矩阵通常很稀疏。此时真正昂贵的不是一个标量运算，而是存储、访存以及不必要的填充。本章先说明如何保存稀疏矩阵，再讨论最基本的迭代解法。

10.1 CSC 存储与稀疏矩阵乘法

压缩列存储（CSC）使用三个数组：

• values：按列保存非零元；
• row_indices：对应的行号；
• col_ptr：每一列在前两个数组中的起止位置。

计算 $y=Ax$ 时，对第 $j$ 列的非零元做

$$y_i+=a_{ij}{x}_j.$$

若矩阵有 $\mathrm{nnz}(A)$ 个非零元，乘法的运算量和存储量都是 $O(\mathrm{nnz}(A))$。计算 $A^{T}x$ 时仍可逐列遍历，只是每列累加一个内积。

稀疏算法的第一原则：围绕非零元组织计算，不要先把矩阵转成稠密格式。

坐标格式（COO）用三元组 $(i,j,a_{ij})$ 交换数据，易于增量装配；CSC/CSR 则把列指针或行指针压缩为长度 $n+1$ 的数组，更适合计算。有限元装配常先生成 COO，再排序、合并重复位置并转为压缩格式。索引从 0 还是 1 开始必须在接口处明确。

选择格式要服从操作：CSR 适合逐行计算 $Ax$，CSC 适合逐列更新、访问 $A^T$ 或许多稀疏直接法。格式转换本身有成本，迭代中不要来回转换。

10.2 矩阵分裂

将线性方程组写成

$$A=M-N,$$

其中 $M$ 容易求解。于是

$$x^{(k+1)}=M^{-1}N{x}^{(k)}+M^{-1}b=B{x}^{(k)}+c.$$

误差满足 $e^{(k+1)}=B{e}^{(k)}$，因此对任意初值收敛的充要条件是

$$\rho (B)<1.$$

这个判据精确但未必容易直接验证，实际常用对角占优、正定性或范数估计给出充分条件。

10.3 Jacobi、Gauss–Seidel 与 SOR

按约定写

$$A=D-L-U,$$

其中 $D$ 是对角部分，$L,U$ 分别是严格下、上三角部分。

Jacobi

$$D{x}^{(k+1)}=(L+U)x^{(k)}+b,B_J=D^{-1}(L+U).$$

各分量只使用上一轮结果，容易并行。

Gauss–Seidel

$$(D-L)x^{(k+1)}=U{x}^{(k)}+b,B_{GS}=(D-L{)}^{-1}U.$$

新得到的分量会立即参与后续计算，串行依赖更强，但通常比 Jacobi 快。

SOR

引入松弛参数 $\omega$：

$$(D-\omega L)x^{(k+1)}=[(1-\omega )D+\omega U]x^{(k)}+\omega b.$$

$\omega =1$ 即 Gauss–Seidel；$0<\omega <1$ 为欠松弛，$1<\omega <2$ 为超松弛。对称正定情形中，$0<\omega <2$ 是 SOR 收敛的基本范围，但最佳参数仍依赖问题结构。

10.4 何时可以保证收敛

• 若 $A$ 严格行对角占优，则 Jacobi 与 Gauss–Seidel 均收敛。
• 若 $A$ 对称正定，则 Gauss–Seidel 收敛。
• 即便残差单调下降，误差也未必按同样速度下降；两者由 $A$ 的条件数联系。

严格行对角占优下，可用最大模分量证明。若 $Bx=\lambda x$ 且 $|x_i|=\parallel x{\parallel }_{\mathrm{\infty }}$，把第 $i$ 行移项并利用

$$|a_{ii}|>\sum _{j\ne i}|a_{ij}|$$

即可推出 $|\lambda |<1$。这类证明揭示了收敛与迭代矩阵谱半径的联系，而不是仅凭“使用新值应该更快”的直觉。

10.5 二维 Poisson 方程

在单位正方形上用五点差分离散

$$-\mathrm{\Delta }u=f$$

得到

$$\frac{1}{h^2}(4u_{ij}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1})=f_{ij}.$$

所得矩阵稀疏、对称、正定。随着网格加密，条件数约按 $h^{-2}$ 增长，所以简单定常迭代会明显变慢。这解释了为什么它们更常被用作多重网格的光滑器或 Krylov 方法的预处理步骤。

在规则网格上，Jacobi 对高频误差衰减快、对低频误差衰减慢。红黑排序把棋盘相邻点染成不同颜色：同色点之间不直接耦合，可以并行更新；完成红、黑两次更新相当于一种重排后的 Gauss-Seidel。

10.6 停机标准与实现检查

常用相对残差

$$\frac{\parallel b-A{x}^{(k)}\parallel }{\parallel b\parallel }$$

作为停机指标，并同时设置最大迭代次数。实现时应记录残差历史；若出现停滞或增长，检查矩阵分裂、索引和松弛参数，而不是只增加迭代次数。

10.7 SOR 的参数与谱

SOR 的迭代矩阵为

$$B_{\omega }=(D-\omega L{)}^{-1}[(1-\omega )D+\omega U].$$

最佳 $\omega$ 依赖谱信息。对规则 Poisson 问题可以从 Jacobi 谱半径估计最佳参数；一般问题则很难事先知道。$\omega >1$ 并不自动加速，过度松弛可能导致振荡甚至发散。

10.8 定常迭代就是预条件 Richardson

从分裂 $A=M-N$ 得

$$x_{k+1}=x_k+M^{-1}(b-A{x}_k).$$

所以每一步都在用 $M^{-1}$ 近似 $A^{-1}$ 作用于残差。这个视角把 Jacobi、Gauss-Seidel、不完全分解和现代预条件统一起来。单独使用时它是定常迭代；把 $M^{-1}$ 嵌入 Krylov 方法时，它改变搜索空间中的谱几何。

10.9 多重网格思想

细网格上的低频误差在粗网格上看起来更高频，因此可以通过三步循环消除：

1. 在细网格做少量平滑，压低高频误差；
2. 把残差限制到粗网格，近似求解误差方程；
3. 插值校正回细网格，再做后平滑。

若各层工作量按几何级数下降，一个 V-cycle 可达到近 $O(n)$ 成本。多重网格既可独立求解，也可作为 CG/GMRES 的强预处理器。

10.10 稀疏实现的验证

• 用小型稠密矩阵核对 COO 到 CSC/CSR 的转换和转置乘法。
• 检查重复坐标是否正确求和、空行/空列的指针是否一致。
• 区分递推残差与显式残差，周期性重算 $b-A{x}_k$。
• 报告每步 $O(\mathrm{nnz}(A))$ 成本之外，还要报告迭代次数和预处理成本。

自检

1. 为什么 $\rho (B)<1$ 与对任意初值收敛等价？
2. CSC 格式下如何分别实现 $Ax$ 与 $A^{T}x$？
3. 为什么细网格上的 Poisson 问题会让 Jacobi 变慢？
4. 如何从矩阵分裂看出定常迭代本质上是预条件残差校正？
5. 多重网格为什么能同时处理高频与低频误差？