第 XIII–XXI 讲 · LU 与 Cholesky：解方程的主力军团

解 $\mathbf{Ax}=\mathbf{b}$ 的主力是 LU 分解 $A=LU$。这一组讲把它讲透：高斯消元是 LU 的"过程视角"，LU 是消元的"结果视角"；Cholesky 是对称正定阵的专属快车道。

高斯消元 = LU + 排列

消元的每一步：用主元行把下方各行的首元素清零。把消元乘子存进 $L$（单位下三角）、消元结果存进 $U$（上三角）：

$P$ 是行交换（部分主元：每步选本列绝对值最大的做主元）。为什么要选主元？不选主元时小主元做除数会把误差放大 $1/|l_{ik}|$ 倍（增长因子），甚至中途除零。部分主元在实践中几乎总稳（尽管存在人为的反例矩阵），代价仅 $O(n^2)$ 的比较。求解变成两步三角回代，各 $\frac{1}{2}n^2$ flop；分解本身 $\frac{2}{3}n^3$ flop——一次分解，多次右端是 LU 的商业模式。

Cholesky：对称正定的捷径

对称正定（SPD）阵不需要主元选择，且 $U$ 由 $L$ 转置决定：

代价只有高斯消元的一半（$\frac{1}{3}n^3$），且过程数值稳定到"不需要主元"的程度。SPD 检验本身也是用途：分解中根号下出现负数 = 矩阵非正定 = 你的协方差阵/刚度阵有毛病。

增长因子与稳定性

高斯消元（带部分主元）的向后稳定性由增长因子 $\rho = \max|u_{ij}|/\max|a_{ij}|$ 控制。理论上 $\rho$ 可达 $2^{n-1}$，实践中极罕见——这是数值线性代数里"理论反例与实践安全"的经典错位。稳健主义者的保险：随机扰一下输入或做增长监控。

老师的提醒

• 同一个 $\mathbf{A}$、许多个 $\mathbf{b}$（时序仿真、参数扫描）：先 LU 分解一次，再反复三角求解，省一个数量级。
• 看到"矩阵接近奇异"警告：检查是否本构参数设错（如泊松比 0.5 卡死）、单位混排、约束不足——条件数问题九成是建模错误。

自检

• $PA=LU$ 里三个因子各存什么？为什么不写成 $P A P^T = LU$？
• Cholesky 快在哪两点？
• 增长因子失控的实践风险有多大？

更深一层：LU 的商业本质是"分解一次、反复变现"

$PA = LU$ 的 $\frac23n^3$ flop 是一笔固定投资，而三角回代每次只要 $2n^2$——投资回收点在第二次求解右端。这个商业模式决定了它的统治区：时域仿真每步解同矩阵方程（上万个右端）、参数扫描、灵敏度分析——LU 是不二之选。反过来，矩阵每次都变（非线性牛顿每步新雅可比）时，LU 的投资无法摊销，此时要么接受重复分解（牛顿法标配），要么换 Krylov 类"不用分解"的方法。选分解法还是迭代法，第一问是"矩阵会变吗"——这一个问题就分流了 80% 的场景。

部分主元的另一个身份也值得说：它等价于解方程时自动做行缩放的贪心策略——每步用当前列最大的元素当除数，把乘子压在 $[-1, 1]$。这个贪心没有全局最优性保证（完全主元才是），但 $O(n^2)$ 的额外比较换来实践中的极高安全性——工程算法设计的"性价比美学"典范。

推导细读：Cholesky 为什么敢不要主元

对称正定阵免主元的证明值得走一遍：分解过程中任意前主子式（如 $2\times2$ 块）都是正定的（SPD 的性质），而正定块的主对角元 $\alpha_{jj} - \sum l_{jk}^2 > 0$ 自动成立——根号下永远是正数，不需要任何交换来救场。这个"结构性安全"还可以反向利用：分解过程中根号下出现负数 ⟹ 矩阵非正定。于是 Cholesky 成了SPD 的验钞机：协方差阵合法性、刚度阵约束完整性、核函数合法性，全都可以"Cholesky 一下"来检验。数值分析里很少有一个算法同时是最好的求解器与最好的检验器——Cholesky 是那个例外。

在真实世界里：求解器栈的现代分工

今天的工程计算里，LU/Cholesky 的分工长这样：稠密小中型（$n \le 5\times10^4$）——LAPACK 的 getrf/potrf（多线程、分块、SIMD 全开）；大型稀疏对称——CHOLMOD/MUMPS（排序 + 符号分解 + 数值分解三段式）；迭代法的预条件——不完全 Cholesky（ICT）作为 CG 的前置。还有一类"隐藏的 LU"：隐式 ODE 求解器（BDF）每步的雅可比分解、MPC 控制器的 QP 内点法——工业世界里每秒都有亿级 LU 在运行，只是藏在各求解器内核里。

常见误区与辨析

• "LU 结果异常 = 算法问题"：先查条件数与增长因子；病态问题 LU 照样向后稳定，错的源头常在建模。
• "并行 LU 线性加速"：三角求解部分依赖链长、并行度差——LU 的并行效率天然低于矩阵乘；这是"分块-管道"复杂调度的存在理由。
• "矩阵对称就一定 Cholesky"：还得正定；对称不定硬上 Cholesky 会在根号下遇负数——那正是诊断信息，别绕过去。

专题深入：Cholesky 的"验钞机"身份

Cholesky 分解的深层身份是对称正定矩阵的验钞机：分解过程中根号下出现负数 ⟹ 矩阵非正定——协方差阵合法性、刚度阵约束完整性、核函数合法性，全部可以"Cholesky 一下"来检验。SPD 的结构性安全（任意前主子式正定 ⟹ 根号下永远正数）让 Cholesky 免主元——这个"不需要主元选择"的特权来自正定的代数结构。数值分析里很少有一个算法同时是最好的求解器与最好的检验器——Cholesky 是那个例外。

数学细节展开：$PA = LU$ 的"行交换 = 排列矩阵"

部分主元的代数表述 $PA = LU$ 值得从"排列矩阵"角度再看一遍：$P$ 是行交换的乘积（正交置换矩阵），$L$ 是单位下三角（消元乘子的汇编）、$U$ 是上三角（消元结果）。一次分解、多次右端是 LU 的商业模式：分解 $\frac23n^3$ 一次性投资、回代 $2n^2$ 每次变现——时域仿真与参数扫描的算力结构由这个商业模式决定。迭代精化（残差用双精度算、LU 解修正方程）是"先直接法打底、再精化"的标准工业流程——混合精度迭代精化是超算 Linpack 基准的接力棒技术。

工程案例细节：带状与结构矩阵的"速度账"

特殊结构矩阵的速度账：三对角 Thomas 算法 $8n$ flop（vs 满阵 $O(n^3)$，$n = 10^5$ 时差 10 个数量级）；带状（有限元/差分的稀疏产物——带宽内工作）；Toeplitz（Levinson 递推 $O(n^2)$——信号处理的移位不变性产物）；Vandermonde（插值的特殊结构——快速算法与病态性并存）。"结构换速度"是矩阵计算的第一红利——识别矩阵的结构类型比优化通用算法更重要。

常见误区与辨析（补充）

• "永远不要显式求逆"：$A^{-1}\mathbf{b}$ 的误差是解方程做法的数倍（$A^{-1}$ 的条件作用 + 额外 $O(n^3)$ 浪费）——解方程永远走分解 + 回代。
• "矩阵接近奇异"警告时先查建模：本构参数设错、单位混排、约束不足——条件数问题九成是建模错误。