第 5–7 章 · 线性方程组的直接解法

与 Trefethen 讲故事的写法不同，Golub-Van Loan 把直接法的每个细节都摊开：存储格式、循环次序、分块、主元策略。这三章是"想把求解器写对写快"的人的必经之路。

高斯消元的完整工程账

$PA = LU$ 的三种实现次序（kij/ijk…）对应"消元/分解/右看"的不同循环组织；flop 均为 $\frac{2}{3}n^3$，但缓存命中率天差地别。现代实现一律分块：把矩阵切成 $b\times b$ 块，内层用高缓存效率的矩阵-矩阵乘（BLAS3），外层处理更新。这就是"教科书算法"与"LAPACK 差 50 倍速度"的原因。

主元策略的完整谱系

• 部分主元：行内选最大，$PA = LU$，实践中稳。
• 完全主元：全余子阵选最大，更稳但搜索代价高，仅病态特例使用。
• 对称正定：免主元 Cholesky，$\frac{1}{3}n^3$。
• 对称不定：对角主元不可靠，Bunch-Kaufman $1\times1/2\times2$ 主元，$LDL^T$ 形。
• 稀疏：主元以"保持稀疏性"优先（MARKOWITZ 准则），数值稳定性让位，配合阈值主元折中。

病态与迭代精化

解完不满意？迭代精化：残差 $\mathbf{r} = \mathbf{b}-A\mathbf{x}_c$ 用 LU 再解修正方程 $A\mathbf{d} = \mathbf{r}$，$\mathbf{x} \leftarrow \mathbf{x}_c+\mathbf{d}$，两三轮就能把精度拉回（条件数不太离谱时）到机器精度——前提是残差用双精度计算。这是"先直接法打底、再精化"的标准工业流程。

老师的提醒

• 永远不要显式求 $A^{-1}$ 再乘：$\mathbf{x} = A^{-1}\mathbf{b}$ 的误差是解方程做法的数倍（$A^{-1}$ 的条件作用 + 额外 $O(n^3)$ 浪费）。
• 带状矩阵要用带状存储接口，别当满阵存——三对角与满阵的求解成本差三个数量级。

自检

• 分块 LU 优化的是 flop 还是别的什么？
• 对称不定矩阵为什么不能照搬 Cholesky？
• 迭代精化的两个前提是什么？

更深一层：主元策略的"安全-成本"光谱

主元选择不是开关而是光谱，一端是数值安全、另一端是成本与稀疏性。完全主元最安全（增长因子有界），但每步要扫全余子阵——$O(n^2)$ 次比较的搜索税；部分主元只扫当前列，实践中安全（增长因子爆炸的矩阵在自然界几乎不存在）；稀疏世界再退一步：先保稀疏后保数值（MARKOWITZ 准则选"填充最少"的主元，加数值阈值兜底）——因为稀疏 LU 的成本由填充量主导，主元数值上的小风险换填充的大节约是笔好买卖。对称不定情形还有专属方案（Bunch-Kaufman 的 $1\times1/2\times2$ 混合主元）。主元策略 = 问题的数值难度与结构的工程价值之间的合约条款——读第 3 章时抓住这条光谱，各种变体就串成一条线。

推导细读：迭代精化为什么"几乎免费"

迭代精化的账本：残差 $r = b - Ax_c$（$O(n^2)$，用原矩阵）、解修正方程（已有 LU，$O(n^2)$ 回代）、更新 $x$。每轮 $O(n^2)$——对比重新分解的 $O(n^3)$ 几乎免费。它救回来的精度有多大？单精度分解 + 双精度残差可以逼近双精度直接求解的精度（这是"混合精度迭代精化"（HPL-AI/Exascale 项目）的核心思想，超级计算机 Linpack 基准的接力棒技术）。工程 takeaway：精度不是一次分解决定的"一次购"，而是可以事后用精化迭代"分期补齐"的预算项。

在真实世界里：稀疏直接法的三段式工业流水线

稀疏 LU/Cholesky 的现代实现是三段式：排序（嵌套剖分/AMD 决定填充结构，只做一次）、符号分解（按排序预算出非零位置，纯图计算）、数值分解（按既定结构算数，可并行可复用）——同一个矩阵多次分解（时域仿真）时，一、二段成本全摊销。工业求解器（MUMPS、PARDISO、SuperLU、CHOLMOD）的差异就在这三段的实现深度。选型速查：对称正定 → CHOLMOD 系；一般对称 → MA57 类；非对称 → MUMPS/SuperLU；配合迭代 → ILU 系预条件。

专题深入：主元策略的"完整光谱"再展开

主元策略的完整光谱值得从"数值安全-成本"轴重新排列：完全主元（全余子阵选最大——最稳但搜索税 $O(n^2)$ 每步）→ 部分主元（行内选最大——实践安全、成本低）→ 稀疏主元（MARKOWITZ 准则——保持稀疏性优先，数值稳定性让位，阈值主元折中）→ 对称不定（Bunch-Kaufman $1\times1/2\times2$ 主元——对称结构的保存者）。主元策略 = 问题的数值难度与结构的工程价值之间的合约条款。

数学细节展开：平衡（balancing）的"对角相似"

平衡（balancing）是对角相似缩放 $D^{-1}AD$——先验地把"范数失衡"洗掉。$A = \begin{pmatrix}10^6 & 10^6\\10^{-6} & 2\times10^{-6}\end{pmatrix}$ 这类矩阵的行列尺度相差 $10^{12}$ 倍，特征值计算（尤其非对称的特征向量）对这种失衡极敏感。平衡的工程地位：LAPACK 的 gebal 在特征值求解前自动做这一步——花小钱（$O(n^2)$）办大事（精度提升几个数量级）。

工程案例细节：带状矩阵的"存储红利"

带状矩阵的存储红利值得算清：满阵存 $10^6$ 自由度问题需要 TB 级内存，按带宽或稀疏 CSR 格式只需 MB-GB 级。网格节点编号决定带宽——前处理里的编号优化（RCM 排序）直接省内存、提速。三对角与满阵的求解成本差三个数量级——带状矩阵要用带状存储接口，别当满阵存。

常见误区与辨析（补充）

• "稀疏直接法免费"：三维问题的填充增长让纯直接法吃紧——排序只能缓解不能消除；大三维问题最终都走到迭代法门口。
• "条件数问题是算法问题"：条件数是问题的属性——换算法不换问题；先量 $\kappa$ 再决定对策。