4. 稀疏Cholesky分解与排序算法

对对称正定矩阵 A 做 Cholesky 分解，目标是得到下三角矩阵 L，使 A = L·Lᵀ。稀疏情形下，真正的难点不在分解公式本身，而在于填充（fill-in）：原本为零的位置在消元过程中变成非零元。填充会直接抬高 L 的存储量和浮点运算量，因此排序算法的核心任务就是通过重排未知量顺序，尽量抑制填充。

先回顾一个关键事实：Cholesky 分解的结果与排序强相关。若对 A 施加对称置换 P·A·Pᵀ，分解得到的 L 只是原 L 的行列重排，数学上等价，但非零元个数可能相差数倍。这就是“同一个矩阵、不同排序、不同代价”的根源。工程上常用的两类排序是最小度排序（AMD，Approximate Minimum Degree）和嵌套剖分（ND，Nested Dissection）。

填充从哪来：消去树与符号分解

在数值分解之前，可以先做一次符号分解，只依据矩阵的稀疏结构预测 L 的非零模式，不进行浮点运算。符号分解的产物是消去树：每个节点对应一个未知量，父节点是它在消元时“传染”到的第一个更高编号节点。消去树刻画了填充的传播路径——某个节点的所有后代在消元后都会在该节点处产生填充。因此，排序本质上是在重塑消去树的形状：树越“矮胖”、子树越均衡，填充通常越少。

AMD 的思路是贪心的：每一步选择当前图中度最小的节点消去，并近似估计消去后产生的团（clique）规模，避免精确计算度的高开销。它对不规则稀疏矩阵、有限元网格都表现稳健，是很多稀疏直接求解器的默认排序。ND 则采用分治：递归地寻找小顶点分隔集，把图切成若干互不连通的子图，先分别对子图排序，最后处理分隔集。对来自二维、三维网格的矩阵，ND 的渐近填充更优，尤其适合结构规整的大规模问题。

实战：用 AMD 排序优化 Cholesky 分解并比较填充

下面用 Python 生态中常用的 scipy.sparse 与 sksparse.cholmod（或 scipy.sparse.linalg.splu 的对称变体）演示完整流程。若环境没有 sksparse，可用 scipy.sparse.csgraph.reverse_cuthill_mckee 作为对照排序，但 AMD 的填充抑制通常更好。示例矩阵采用二维五点差分拉普拉斯矩阵，规模取 n=40×40=1600，这是有限元与网格计算中最典型的对称正定稀疏结构。

import numpy as np
import scipy.sparse as sp
from sksparse.cholmod import analyze, cholesky

# 构造二维五点拉普拉斯矩阵（对称正定，加对角占优保证正定）
n = 40
I = sp.identity(n, format='csr')
T = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csr')
A = sp.kron(I, T) + sp.kron(T, I)          # 1600 x 1600
A = A + sp.identity(A.shape[0]) * 0.1       # 轻微对角加强

# 1) 不排序（自然顺序）
f_nat = analyze(A, mode='simplicial')
L_nat = f_nat.cholesky(A)
nnz_nat = L_nat.L().nnz

# 2) AMD 排序
f_amd = analyze(A, mode='simplicial', ordering_method='amd')
L_amd = f_amd.cholesky(A)
nnz_amd = L_amd.L().nnz

print("自然顺序 L 非零元:", nnz_nat)
print("AMD 排序 L 非零元:", nnz_amd)
print("填充减少比例: %.1f%%" % (100*(nnz_nat-nnz_amd)/nnz_nat))

这段代码的关键在于 analyze 与 cholesky 的分离：analyze 只做符号分析，确定排序与消去树，cholesky 才做数值分解。这样可以在同一结构上反复更换排序、比较填充，而不必重复数值计算。运行后应观察三个量：L 的非零元个数、分解耗时、以及回代求解的残差 ‖Ax−b‖/‖b‖。前两者反映填充收益，第三者确认排序没有破坏正确性——理论上残差应在 1e-12 量级（双精度）。

AMD 与 ND 的选用判断

两者并非互斥，很多求解器会先做一次 ND 得到粗粒度分隔，再在子图内部用 AMD 细化。实际选型可参考矩阵来源：二维/三维网格、规则有限元问题优先 ND；不规则图、电路网络、一般稀疏对称系统优先 AMD。若不确定，最稳妥的做法是两种排序都跑一次符号分析，直接比较 L 的非零元个数，选小者。符号分析本身很便宜，这个“试跑比较”在工程上完全可接受。

最后提醒一个易错点：排序只影响 L 的稀疏模式，不影响 A 的数值。若分解失败（出现负主元或数值不稳定），问题在矩阵本身是否正定或条件数，而非排序。排序是“省资源”的手段，不是“修数值”的手段，二者要分清。