3. 直接法基础

直接法的核心思想是：通过有限步的代数消元，把原方程组 Ax = b 化为易于回代求解的三角形式，且在不计舍入误差的前提下得到精确解。对稠密矩阵而言，高斯消元与 LU 分解本质上是同一件事的两种表述——前者作用于增广矩阵，后者把消元过程中产生的乘数保存下来，形成 A = LU。理解稠密情形是理解稀疏情形的起点，因为稀疏 LU 的全部困难，都来自"消元会制造新的非零元"这一事实。

稠密 LU 分解回顾

设 A 为 n×n 非奇异矩阵。带部分主元的高斯消元在第 k 步做两件事：先在列 k 的第 k 行至第 n 行中选出绝对值最大的元素作为主元并交换行，然后用主元行消去下方各行的第 k 列。若把第 k 步使用的乘数 l_ik = a_ik / a_kk 记录到矩阵的严格下三角部分，消元结束时就得到

P A = L U
L：单位下三角，元素为各步乘数
U：上三角，含主元与消元后的右上块
P：行置换矩阵，记录主元交换

求解时依次做前代 L y = P b 与回代 U x = y。稠密 LU 的计算量约为 2n³/3 次浮点运算，存储量为 n²。关键在于：即使 A 中大量元素为零，只要按常规消元顺序执行，L 和 U 的"结构零"位置也可能被填入非零值，这就是填充。

稀疏 LU 的填充（fill-in）现象

考虑一个 5×5 的稀疏矩阵，其非零结构为三对角加一条远端耦合。消去第 1 列时，第 1 行在第 3 列有非零元，而第 2 行在第 1 列有非零元，于是第 2 行第 3 列原本为零的位置被更新为非零——这就是一次填充。填充一旦产生，后续消元又会在新非零元之间继续制造填充，形成连锁反应。

填充的代价是双重的：存储上，L、U 的非零元数量可能远超 A；计算上，消元只对非零元操作，填充越多，浮点运算越多。对二维泊松矩阵这类典型稀疏结构，若按自然行序消元，填充量可达 O(n² log n) 量级；而采用嵌套剖分等更优排序，可降到接近 O(n^{3/2})。因此，稀疏直接法的工程核心不是"会不会消元"，而是"如何排序以减少填充"。

符号分解与数值分解

稀疏 LU 通常拆成两个阶段。符号分解只处理非零模式：给定矩阵结构、选定的排序（如 AMD、COLAMD、嵌套剖分）与主元策略，计算出 L、U 的稀疏结构、各列的非零元个数，并预先分配存储。数值分解则在已知结构上执行实际的浮点消元，只访问符号阶段确定的位置。

这种分离带来两个好处：其一，结构分析可以复用——当矩阵模式不变而数值改变时（如非线性迭代、时步推进中的雅可比矩阵），只需重做数值分解；其二，内存可以一次性精确分配，避免动态扩容。工程上常用 Symbolic LU 与 Numeric LU 两个接口来体现这一划分。需要注意的是，若采用阈值部分主元（threshold pivoting），主元选择依赖数值，符号阶段只能给出保守上界，实际结构可能略小。

实战：对二维泊松矩阵做 LU 分解并可视化填充

取单位正方形上的五点差分泊松方程，网格 m×m，按行优先编号，得到 n = m² 阶矩阵。以 m = 5 为例，n = 25，矩阵为块三对角：对角块是三对角矩阵，非对角块为负单位阵。下面用 Python 的 SciPy 演示符号与数值分解，并统计填充量。

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu

m = 5
n = m * m
# 构造二维五点差分泊松矩阵
T = sp.diags([-1, 4, -1], [-1, 0, 1], shape=(m, m))
I = sp.identity(m)
A = sp.kron(I, T) + sp.kron(sp.diags([-1, 1], [-1, 1], shape=(m, m)), I)
A = A.tocsc()

# 数值 LU（内部先做符号分析）
lu = splu(A, permc_spec="NATURAL")   # 自然序，便于观察填充
L, U = lu.L, lu.U

print("A 非零元:", A.nnz)
print("L 非零元:", L.nnz, " U 非零元:", U.nnz)
print("填充总量:", L.nnz + U.nnz - A.nnz - n)

在自然序下，m = 5 时 A 有 105 个非零元，而 L 与 U 合计约 200 余个非零元，填充量接近一倍。若把 permc_spec 改为 "MMD_AT_PLUS_A"（近似最小度排序），填充量会明显下降。判断填充是否严重，可比较 (L.nnz + U.nnz) / A.nnz 这一比值：接近 1 说明几乎无填充，比值越大说明排序越需要改进。

可视化填充最直接的方式是把 A、L+U 的非零模式画成散点图，用不同颜色区分原有非零元与填充元。填充元的位置恰好是 L+U 模式中不属于 A 的位置，可用集合差 set(zip(*LU.nonzero())) - set(zip(*A.nonzero())) 提取。对 m = 20（n = 400）重复上述实验，自然序的填充会沿带宽方向形成明显的"带状膨胀"，而嵌套剖分排序则呈现分块递归的填充结构——这正是排序算法要优化的对象。