02 三角求解、Gaussian 消元与 LU 分解

1. 三角系统是直接法的基本内核

对非奇异下三角矩阵 $L$，前代公式为

$$x_i=\frac{1}{{\ell }_{ii}}(b_i-\sum _{j=1}^{i-1}{\ell }_{ij}{x}_j),i=1,\dots ,n.$$

运算量约为 $n^2$。标准浮点实现具有分量后向稳定性：计算结果 $\hat{x}$ 满足

$$(L+\mathrm{\Delta }L)\hat{x}=b,|\mathrm{\Delta }L|\le {\gamma }_n|L|.$$

上三角系统用回代，逻辑完全对称。

前代有两种常见实现方式：逐行做内积，或解出 $x_j$ 后立刻更新后续右端项（saxpy 形式）。二者在精确算术中等价，但访存模式不同；后者容易和分块更新统一，前者则直接对应递推公式。

2. Gaussian 消元就是逐步构造 LU

消元矩阵 $L_k^{-1}$ 把第 $k$ 列主元下方元素消为零：

$$L_{n-1}^{-1}\cdots L_1^{-1}A=U.$$

于是

$$A=LU,L=L_1{L}_2\cdots L_{n-1}.$$

稠密 $n\times n$ 矩阵分解约需 $\frac{2}{3}{n}^3$ 次浮点运算；分解后，每个新的右端项只需两次三角求解，即 $O(n^2)$。

3. 为什么要选主元

若当前主元为零，消元无法继续；若主元很小，乘子会很大，中间元素可能被显著放大。

部分选主元在当前列选择绝对值最大的候选，并交换行，得到

$$PA=LU.$$

完全选主元同时搜索剩余子矩阵并交换行列，代价更高。实践中部分选主元通常兼顾成本与可靠性，是通用稠密求解器的默认方案。

若每一步都在当前列选最大元，则乘子满足 $|{\ell }_{ik}|\le 1$。这控制了 $L$，却不能阻止 Schur 补中的 $U$ 元素增长。rook 选主元在行、列之间交替搜索，是部分选主元与完全选主元之间的折中。

若矩阵严格对角占优，或对称正定并采用 Cholesky，则无需普通的行选主元也有稳定性保证。选主元策略应与矩阵结构匹配。

4. 增长因子与后向稳定性

定义消元过程的增长因子

$$\rho =\frac{\underset{i,j,k}{max}|a_{ij}^{(k)}|}{\underset{i,j}{max}|a_{ij}|}.$$

部分选主元 LU 的后向误差受 $\rho$ 控制：若 $\rho$ 温和，则计算得到的 $\hat{L},\hat{U}$ 可解释为邻近矩阵 $A+\mathrm{\Delta }A$ 的精确分解。理论最坏情形可能很大，但常见数据上通常远好于该上界。

5. 迭代改进

已有近似解 $x_k$ 时，计算高质量残差

$$r_k=b-A{x}_k,$$

再用已有 LU 因子解

$$A{d}_k=r_k,x_{k+1}=x_k+d_k.$$

迭代改进用较低成本纠正解。如果残差用更高精度计算，甚至可以从低精度分解获得更高精度解；能否成功仍取决于条件数和因子分解质量。

6. 结构与实现

• 多个右端项：分解一次，重复前代/回代。
• 带状矩阵：若消元不产生额外填充，可把复杂度从 $O(n^3)$ 降到与带宽相关。
• 不要显式计算 $A^{-1}$ 再乘 $b$；直接解方程更快且更稳。
• 高性能实现应尽量用矩阵-矩阵更新，把计算交给 Level-3 BLAS。

7. 分块 LU 与 Schur 补

将矩阵分块为

$$A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right).$$

在 $A_{11}=L_{11}{U}_{11}$ 后，计算

$$U_{12}=L_{11}^{-1}{A}_{12},L_{21}=A_{21}{U}_{11}^{-1},$$

再更新 Schur 补

$$S=A_{22}-L_{21}{U}_{12}.$$

最后递归分解 $S$。绝大部分运算集中在矩阵乘法 $L_{21}{U}_{12}$，这正是分块算法速度高的原因。

8. BLAS 层次与算术强度

层次	典型操作	运算量	数据复用
Level 1	向量加法、内积	$O(n)$	低
Level 2	矩阵-向量乘法	$O(n^2)$	中
Level 3	矩阵-矩阵乘法	$O(n^3)$	高

现代处理器常受内存带宽限制。分块 LU 并没有改变 $O(n^3)$ 运算量，却把大量 Level-2 更新改成 Level-3 更新，使数据留在缓存中重复使用。数值算法的实际性能不能只看浮点运算次数，还要看通信量。

9. 后向误差与增长因子的联系

计算因子通常满足

$$PA+\mathrm{\Delta }A=\hat{L}\hat{U},|\mathrm{\Delta }A|\lesssim {\gamma }_n|\hat{L}||\hat{U}|.$$

所以稳定性的关键不是单看 $\parallel L\parallel$ 或 $\parallel U\parallel$，而是看它们的分量乘积相对于 $|A|$ 有多大。增长因子是便于理解的代理指标。最终解还要经过两次三角求解，其后向误差可继续吸收到 $\mathrm{\Delta }A$ 中。

10. 迭代改进何时有效

设低精度分解的单位舍入误差为 $u_f$，工作精度为 $u$。粗略地说，若

$$\kappa (A)u_f<1,$$

校正方程仍能提供有效方向；若残差用更高精度计算，混合精度迭代改进有机会达到工作精度。实作时应监控缩放后的分量后向误差，而不是只看校正量变小。

11. 自检

• [ ] 能写出前代/回代并估算复杂度。
• [ ] 能从消元矩阵推出 $PA=LU$。
• [ ] 能解释小主元、增长因子与稳定性的关系。
• [ ] 知道迭代改进复用的是哪一组因子。
• [ ] 能从分块消元写出 Schur 补更新。
• [ ] 能解释为什么 Level-3 BLAS 更接近机器峰值。

下一章：Cholesky 分解与结构利用。