3.5 高斯消元、初等矩阵与 LU 分解

本节位置：模块 3 · 矩阵与线性映射 ｜ 前置：3.4 ｜ 预计：40 分钟

一句话：高斯消元用三种行变换把矩阵化成三角，LU 分解则把这三步「录下来」存成一个下三角矩阵，供以后反复使用。

0. 场景引入

结构设计院，一栋 60 米高的钢框架要做抗震验算。简化成层模型后自由度是 480 个（80 个节点 × 6 个自由度），刚度矩阵 $K\in {\mathbb{R}}^{480\times 480}$，对称正定。

规范要求做包络设计：恒载与活载并成一个组合，加上 8 个方向的风、3 条地震波，一共 12 个荷载工况。每个工况对应一个右端项 ${\mathbf{f}}_k$，都要解一次

$K$ 只跟结构本身有关，12 个工况共用它。

决定：是对每个 ${\mathbf{f}}_k$ 从头消元 12 次，还是把 $K$ 分解一次（$K=LU$）、后面只做 12 次三角求解？

代价：$n=480$ 时，一次消元约 $\frac{2}{3}{n}^3=7.4\times 10^7$ FLOPs，12 次就是 $8.8\times 10^8$；分解一次加 12 次回代是 $7.4\times 10^7+12\cdot 2\cdot 480^2\approx 8.0\times 10^7$，省 11 倍。这个模型再加密到 $n=1.2\times 10^6$（精细壳单元），一次分解本身就是 $10^{18}$ FLOPs 量级——重复 12 次直接不可行。

比慢更可怕的是漏算：12 个工况少算一个，某根柱的配筋就不足，地震时是倒塌。

问题：消元的过程为什么能「存下来」？主元为 0、或者只是很小时，为什么必须换行？

读完本节，你应该能回答三件事：

1. 三种行变换为什么都不改变解集，却可能改变行列式；
2. $A=LU$ 里 $L$ 存的到底是什么，为什么它恰好是各步的乘数；
3. 为什么实际实现都是 $PA=LU$，不选主元会坏在哪。

1. 直觉

1.1 消元与存档

高斯消元（Gaussian elimination）就是「用一行去减另一行的倍数」，制造零，让矩阵变成阶梯形：主元逐行右移，主元下方全为零。每一步都是一次初等行变换（elementary row operation），而每一种行变换都能写成一个**初等矩阵（elementary matrix）**左乘。消元的过程就是一串初等矩阵连乘。

「高斯消元」这个名字来自 Gauss 1809：他在《Theoria motus》里系统地用逐次消去法解观测方程。

LU 分解把结论记下来：$A=LU$，$L$ 记录「每一步减了谁的第几倍」，$U$ 是消元后的上三角。解方程时先解 $L\mathbf{y}=\mathbf{b}$（前代），再解 $U\mathbf{x}=\mathbf{y}$（回代），都只要 $O(n^2)$。

💡 一句话直觉：消元是「做」，LU 是「把做法存档」；存档一次，换右端项时不用重做。

类比边界：这个存档只对同一个 $A$、不同的 $\mathbf{b}$ 有用。若 $A$ 本身在变（如训练中每步都更新权重），LU 的复用价值就消失了。

1.2 现实对应：拆机器时录下的顺序

消元像拆一台机器：你先拧下 4 号螺栓（消掉第 1 列的下面几项），再拧 7 号，再拧 12 号。每一步都让结构「松一点」，最后机器散成一堆零件（上三角）。

LU 是把拆解顺序录成视频：$L$ 记录「每一步拧的是哪颗螺栓、拧了几圈」，$U$ 是拆完后的零件堆。下次拆同型号的机器，照视频做就行，不用重新摸索。

在结构分析里，这台「机器」就是 $K$：它编码了「哪个节点和哪个节点通过一根梁相连、这根梁有多硬」。$L$ 记录的是「用哪个节点去消哪个节点的约束」，$U$ 是消完之后彼此独立的自由度。

这个类比在哪里失效：视频只对同型号机器有效。结构一改（换一根梁的截面、加一层楼），$K$ 就变了，视频作废，必须重新分解。这也是为什么优化设计里「改一个参数就重算一次」代价很高——每次都要重新录一遍视频。

1.3 量级直觉

消元的代价是 $\frac{2}{3}{n}^3$，三角求解只有 $2n^2$。两者之比是 $n/3$：$n=480$ 时约 160 倍，$n=10^6$ 时约 33 万倍。

这个比值就是「分解一次、复用多少次才划算」的临界点：右端项只有 1 个时，直接消元和 LU 一样贵；有 $k$ 个时，复用大约省 $k$ 倍——直到 $2k{n}^2$ 赶上 $\frac{2}{3}{n}^3$，即 $k\approx n/3$。

2. 定义

2.1 增广矩阵与三种行变换

把 $A\mathbf{x}=\mathbf{b}$ 写成增广矩阵（augmented matrix） $[A\mid \mathbf{b}]$。允许三种行变换：

1. 交换两行：$R_i\leftrightarrow R_j$；
2. 数乘一行：$R_i\leftarrow c{R}_i$，$c\ne 0$；
3. 倍加：$R_i\leftarrow R_i+c{R}_j$，$i\ne j$。

三种变换都不改变方程组的解集。

2.2 初等矩阵

对单位阵做一次行变换得到的矩阵叫初等矩阵。左乘初等矩阵 $E$ 等于对 $A$ 做对应的行变换。例如:

$E_{21}=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ 0& 0& 1\end{array}\right]$

左乘 $A$，就是把 $A$ 的第 2 行减去第 1 行的 2 倍。 初等矩阵都可逆，逆矩阵就是把倍数变号：

$E_{21}^{-1}$ 的第 $(2,1)$ 项是 $+2$

2.3 行阶梯形与简化行阶梯形

• 行阶梯形（row echelon form, REF）：非零行在上；每行首个非零元（主元 pivot）所在列严格右于上一行；主元下方全为 0。
• 简化行阶梯形（reduced row echelon form, RREF）：在 REF 基础上，主元全为 1，且主元所在列的其它元素全为 0。

RREF 唯一，消元路径不唯一。

2.4 LU 分解

若方阵 $A$ 可写成

其中 $L$ 是单位下三角（对角线为 1，下三角为消元乘数），$U$ 是上三角，则称其为 LU 分解（LU decomposition）。若消元过程需要交换行，则写成

$P$ 是置换矩阵（permutation matrix），每行每列恰有一个 1。

3. 为什么这样定义

为什么不直接求 $A^{-1}$？因为求逆是 $O(n^3)$，而且显式逆的数值误差通常比直接消元更大。LU 把一次 $O(n^3)$ 的消元与每次 $O(n^2)$ 的三角求解分开：$k$ 个右端项的总代价是 $\frac{2}{3}{n}^3+2k{n}^2$，而逐个求逆再乘是 $2n^3+2k{n}^2$。

为什么需要 $P$？如果主元位置是 0（或非常小），消元会除以 0 或放大误差。部分主元法（partial pivoting） 每次选当前列绝对值最大的行交换上来，这就是 $PA=LU$ 的来历。把消元写成三角分解、并分析它如何放大舍入误差，最早的系统工作是 Turing 1948。

为什么 $L$ 存的是乘数而不是矩阵乘积？因为消元的每一步都是「单位阵 + 一个乘数」，它们的逆只是把乘数变号；把逆按顺序乘起来，恰好得到一个单位下三角，对角线下方就是各步乘数。这就是「存档」的全部内容。

3.1 反事实：如果定义成「不许换行的 $A=LU$」

$A=LU$ 写起来干净，于是有人规定：LU 分解不允许出现 $P$。取一个不算极端的例子：

精确解：$\det =10^{-5}-1=-0.99999$，所以

不换行，用 4 位有效数字算：

• 乘数 $l_{21}=1/10^{-5}=10^5$；
• $u_{22}=1-l_{21}\cdot 1=1-10^5$。这个 $1$ 在 4 位有效数字里被「吃掉」，得 $-1.000\times 10^5$；
• 前代：$y_2=2-10^5\cdot 1\approx -1.000\times 10^5$；
• 回代：$x_2=y_2/u_{22}=1.000$，$x_1=(1-1\cdot x_2)/10^{-5}=(1-1.000)/10^{-5}=0$。

$x_1$ 得 $0$，真值 $1.00001$。相对误差 100%——不是「精度差一点」，是整个答案没了。

换行之后（部分主元法，交换两行）：$l_{21}=10^{-5}$，$u_{22}=1-10^{-5}\approx 1$，回代得 $x=(1.000,1.000)$，相对误差 $10^{-5}$，正好是 4 位有效数字的极限。

结论：$A=LU$ 这个定义在两个地方不够用——主元为 0 时直接除零，主元很小但非 0 时精度悄悄崩掉。所以工程上的定义是 $PA=LU$：$P$ 不是装饰，它是「先把最大的主元换上来」这个动作的代数记录。部分主元法保证所有乘数满足 $|l_{ij}|\le 1$，一旦算出来超过 1，就说明该换行了。

4. 性质与推导

4.1 行变换不改变解集

设 $E$ 是初等矩阵，则 $A\mathbf{x}=\mathbf{b}⟺EA\mathbf{x}=E\mathbf{b}$（两边左乘 $E$）。反过来，因为 $E$ 可逆，$EA\mathbf{x}=E\mathbf{b}\Rightarrow A\mathbf{x}=\mathbf{b}$。所以解集完全不变。

4.2 从消元到 LU

对 $A$ 做消元，相当于左乘一串初等矩阵：

移项得 $A=E_{21}^{-1}{E}_{31}^{-1}{E}_{32}^{-1}U$。记 $L=E_{21}^{-1}{E}_{31}^{-1}{E}_{32}^{-1}$。每个 $E_{ij}^{-1}$ 是单位下三角，单位下三角对乘法封闭，所以 $L$ 也是单位下三角，且其第 $(i,j)$ 项就是消去 $A_{ij}$ 时用的乘数。于是 $A=LU$。

4.3 复杂度

• 消元：第 $k$ 步要更新约 $(n-k{)}^2$ 个元素，总计约 ${\sum }_k(n-k{)}^2\approx n^3/3$ 次乘加（每次乘加 2 FLOPs），即约 $\frac{2}{3}{n}^3$ FLOPs。
• 前代 + 回代：每个右端项约 $n^2$ 次乘加，即 $2n^2$ FLOPs。
• 显式求逆：约 $2n^3$ FLOPs。

所以 $n=4096$ 时，一次 LU 约 $4.6\times 10^{10}$ FLOPs，而每个右端项只要约 $3.4\times 10^7$ FLOPs。这套 $\frac{2}{3}{n}^3+2k{n}^2$ 的计数与 $PA=LU$ 的标准表述见 Golub & Van Loan 第 3 章。

4.4 行列式与逆

由 $PA=LU$，$\det (P)\det (A)=\det (U)$。$\det (P)=\pm 1$（交换次数的奇偶），$\det (U)$ 是对角线乘积。所以 $\det (A)=\pm {\prod }_i{u}_{ii}$。求逆则可以逐列解 $A\mathbf{x}={\mathbf{e}}_j$，复用同一个 LU。

5. 例子

5.1 手算：$3\times 3$ 消元与 LU

取

消元：$R_2\leftarrow R_2-2R_1$ 得 $(0,1,1)$；$R_3\leftarrow R_3-3R_1$ 得 $(0,1,0)$；再 $R_3\leftarrow R_3-R_2$ 得 $(0,0,-1)$。于是

验证 $LU$ 的第三行：$3(1,2,3)+1(0,1,1)+1(0,0,-1)=(3,7,9)$，正是 $A$ 的第三行。

前代 $L\mathbf{y}=\mathbf{b}$：$y_1=1$；$y_2=3-2=1$；$y_3=4-3\cdot 1-1\cdot 1=0$。

回代 $U\mathbf{x}=\mathbf{y}$：$-x_3=0\Rightarrow x_3=0$；$x_2+x_3=1\Rightarrow x_2=1$；$x_1+2x_2+3x_3=1\Rightarrow x_1=-1$。

所以 $\mathbf{x}=(-1,1,0{)}^{\top }$。验证：$A\mathbf{x}=(-1+2,-2+5,-3+7{)}^{\top }=(1,3,4{)}^{\top }$。

行列式：$\det (A)=1\cdot 1\cdot (-1)=-1$，与直接展开一致。

5.2 手算：需要换行的 $PA=LU$

取 $A=\left[\begin{array}{ccc}0& 2& 1\\ 1& 1& 1\\ 2& 1& 0\end{array}\right]$，首个主元为 0。交换第 1、2 行，$P=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$，$PA=\left[\begin{array}{ccc}1& 1& 1\\ 0& 2& 1\\ 2& 1& 0\end{array}\right]$。消元：$R_3\leftarrow R_3-2R_1$ 得 $(0,-1,-2)$；$R_3\leftarrow R_3+\frac{1}{2}{R}_2$ 得 $(0,0,-\frac{3}{2})$。因此

验证第三行：$2(1,1,1)-\frac{1}{2}(0,2,1)+1(0,0,-\frac{3}{2})=(2,1,0)$，等于 $PA$ 的第三行。$\det (PA)=1\cdot 2\cdot (-\frac{3}{2})=-3$，而 $\det (P)=-1$，所以 $\det (A)=3$。

5.3 场景闭环：12 个工况、480 个自由度

输入：$K\in {\mathbb{R}}^{480\times 480}$ 对称正定，fp64 每元素 8 字节，整块占 $480^2\times 8\approx 1.8$ MB；12 个右端项各 $480\times 8=3.8$ KB。按节点编号优化后带宽约 60，非零元约 $480\times 60\approx 2.9\times 10^4$ 个——是稀疏的。

计算：因为 $K$ 对称正定，实际用 Cholesky（$K=L{L}^{\top }$）而不是 LU，代价约 $\frac{1}{3}{n}^3=3.7\times 10^7$ FLOPs，正好是 LU 的一半。12 次三角求解是 $12\times 2n^2=5.5\times 10^6$ FLOPs。合计约 $4.2\times 10^7$。

输出：12 组位移 ${\mathbf{u}}_k$，进而算杆件内力、画应力云图，再按规范做包络——每个构件取 12 个工况里最不利的组合，得到配筋量。

成本：$4.2\times 10^7$ FLOPs 在单核 5 GFLOP/s 下约 8 ms。但稀疏 Cholesky 的瓶颈不是 FLOPs，而是填充（fill-in）：分解过程中原本为零的位置会变成非零，非零元从 $2.9\times 10^4$ 涨到 $10^5$ 量级，内存和访存都涨 3 倍以上。这就是工程上要花力气做节点编号重排（AMD、RCM）的原因——把非零元尽量挤到对角线附近，减少填充。

失败条件与兜底：

• 该换行时没换：$K$ 对称正定时 Cholesky 天然稳定，可以不选主元；但一旦模型里有接触、索、支座沉降，$K$ 变成不定的，Cholesky 会开根号出负数，直接崩。兜底：检测到非正主元就退回带主元的 LU。
• 病态：薄壳、或者钢与橡胶垫并存导致刚度差 $10^6$ 倍时，$\kappa (K)$ 能到 $10^{10}$，位移解的有效位数掉到 6 位以下。兜底：用双精度、并对自由度做尺度化（每个自由度除以其典型量级）再分解。
• 模型误差远大于数值误差：比起 $\kappa$，更大的误差来源是「梁的边界条件假设错了」。数值上解得很准，物理上是错的。兜底：做收敛性检查（网格加密一倍，看应力是否稳定），并与实测或规范简化公式对照。

6. 常见误区

1. 误区：行变换会改变方程组的解。 纠正：三种行变换都是可逆的，解集不变；但数乘一行会把行列式乘 $c$，交换两行会变号——行列式变了，解没变。

2. 误区：LU 分解总是不用换行。 纠正：主元为 0 时必须换行；主元很小（即使非 0）也应换行以控制误差，所以实际实现都是 $PA=LU$。

3. 误区：$L$ 的对角线也存乘数。 纠正：$L$ 是单位下三角，对角线全是 1；乘数在对角线下方。

7. 概念网络

• 前置：3.4 逆矩阵与可逆性 —— 可逆 ⇔ 主元都不为 0。
• 前置：3.3 矩阵乘法 —— 初等矩阵左乘就是行变换。
• 后续：3.7 四个基本子空间 —— 主元列给出列空间的基。
• 后续：4.2 解的存在性与通解 —— RREF 读出特解与零空间。
• 后续：4.6 条件数与病态 —— 主元选择与误差。
• 关联：10.3 BLAS 与 GEMM —— 为什么直接法在大矩阵上让位于迭代法。

8. 通往大模型：求解器藏在哪一层

出现位置：Transformer 的前向没有显式消元，但训练与分析的许多环节会解线性系统：二阶优化（Newton/K-FAC）、高斯过程回归、线性探针、概率校准、以及某些归一化流。

形状与运算：二阶方法要解 $H\mathbf{p}=-\mathbf{g}$，$H\in {\mathbb{R}}^{D\times D}$ 是 Hessian 近似。$D=4096$ 时一次 LU 约 $\frac{2}{3}{D}^3\approx 4.6\times 10^{10}$ FLOPs，显存要存 $D^2=16.7$M 个浮点数；这就是全量二阶方法在 LLM 上不可行的原因。实际做法是用对角/低秩近似（如 K-FAC 的分块），把 $D\times D$ 拆成小块分别做 Cholesky/LU。

为什么需要：LU/Cholesky 把「解方程」与「分解」分离。若同一个 $H$ 要解多个右端项（多个 batch 的梯度），分解一次即可；若 $H$ 每步都变，就只能用共轭梯度等迭代法（第 10.4 节）。可检验结论：用 NumPy 对随机 $4096\times 4096$ 矩阵计时 lu_factor 与 solve，LU 时间约与 $n^3$ 成正比，每个额外右端项只与 $n^2$ 成正比；$n$ 从 1024 增到 2048，前者约 8 倍，后者约 4 倍。

如果这一步算错，用户会看到什么现象：主元没选好，应力云图会出现棋盘格状跳变或负应力——工程师第一反应是「结构有问题」，实际是舍入噪声被放大。LU 复用搞错（$K$ 改了却还用旧的 $L,U$），位移会系统性偏大或偏小，配筋量跟着偏，这类错误最危险：它不报错，只安静地给出一个「看起来合理」的答案。训练侧的同构现象是二阶优化里用了一个过期的曲率近似，loss 不炸，但收敛停滞，调参的人会误以为是学习率问题。

9. 练习

1. 概念：为什么行变换不改变解集，却可能改变行列式？各举一种变换说明。
2. 计算：对 $A=\left[\begin{array}{ccc}1& 1& 0\\ 2& 1& 1\\ 0& 1& 3\end{array}\right]$ 做 LU 分解（不需换行），并解 $A\mathbf{x}=(2,5,4{)}^{\top }$。
3. 计算：接上题，用 $U$ 的对角线求 $\det (A)$，并与直接展开对比。
4. 计算：$A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]$，写出 $P$、$L$、$U$，验证 $PA=LU$。
5. 开放：如果 $A$ 每 10 步才更新一次，而每步有 1000 个右端项，LU 复用能省多少？请用 $n^3$ 与 $n^2$ 的量级估算。

10. 一页速记

• 三种行变换：交换、数乘、倍加；都不改变解集，但可能改变行列式。
• 物理意义 / 现实对应：消元是「逐个消掉约束」，LU 是把消约束的顺序存档；结构分析里 $K$ 是刚度，$L$ 是约束的消除顺序，$U$ 是解耦后的自由度。
• 初等矩阵左乘 = 行变换；它们都可逆。
• REF/RREF：主元阶梯；RREF 唯一。
• $A=LU$（或 $PA=LU$）：$L$ 单位下三角存乘数，$U$ 上三角。
• 消元 $O(n^3)$，三角求解每个右端项 $O(n^2)$；$\det (A)=\pm \prod u_{ii}$。

11. 参考文献与致谢

这一页的知识不是凭空来的。下面列出本页直接借用的原始工作与权威整理；正文中的人名、定理名已链接到对应条目。链接优先指向 DOI、arXiv 或可长期访问的原文存档。

• 高斯（Carl Friedrich Gauss, 1809）：Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium（1809 原著，Cambridge 2011 影印）。系统使用逐次消去法解观测方程，本页 §2.1 的三种行变换是它的现代形式。Cambridge University Press
• 图灵（Alan M. Turing, 1948）：Rounding-off Errors in Matrix Processes。把消元写成三角分解（$A=LU$，必要时换行）并分析舍入误差，本页「LU + 主元选择」的起点。Quart. J. Mech. Appl. Math.
• Golub & Van Loan（2013）：Matrix Computations（第 4 版）。第 3 章给出 $PA=LU$、部分主元法与 $\frac{2}{3}{n}^3$ FLOPs 的标准表述，本页 §4.3 的计数以它为准。Johns Hopkins University Press

下一篇 → 3.6 转置、对称与分块矩阵

源文件：03-矩阵与线性映射/05-高斯消元与LU.md