数值线性代数（一）：高斯消元、LU 分解与「不要显式求逆」

向量范数：从 L2 到 Lp#

你已经掌握的欧氏长度就是 L2 范数：

其中 $\langle x, y \rangle = \sum_i x_i y_i$ 是内积，余弦相似度 $\cos\theta = \langle x, y \rangle / (\|x\|_2 \|y\|_2)$ 用的就是它。一般地，$L_p$ 范数定义为

常用的是三个特例：

• $p = 1$：L1 范数 $\|x\|_1 = \sum_i |x_i|$，曼哈顿距离；
• $p = 2$：L2 范数，欧氏距离，由内积诱导；
• $p \to \infty$：L∞ 范数 $\|x\|_\infty = \max_i |x_i|$，最大分量。

一个函数要被称为”范数”，必须满足三条公理：非负（$\|x\| \ge 0$，等号当且仅当 $x = 0$）、齐次（$\|cx\| = |c|\,\|x\|$）、三角不等式（$\|x + y\| \le \|x\| + \|y\|$）。直觉上，范数就是”从原点量到 $x$ 的距离”。选择哪种范数取决于场景：稀疏向量常用 L1（压缩感知里是凸松弛），L2 在数值线性代数里最常用，原因是它是酉不变（unitary invariant）的：对任意正交矩阵 $Q$（旋转/反射），$\|Qx\|_2 = \|x\|_2$——旋转不改变长度，这使 L2 范数与矩阵分析配合得最干净。你之后会看到，正交矩阵不放大误差这件事，在 RoPE 那里还会再出现一次。

矩阵范数：单位球被拉伸的最长半径#

向量范数只能量”解”的长度，但误差分析里的放大倍数 $\kappa$ 属于矩阵。定义矩阵的诱导范数（也叫算子范数）为

直观解释：把所有单位长度的向量 $x$ 扔进 $A$，看输出最长能有多长，这个”最大拉伸倍数”就是 $\|A\|$。其中由 L2 范数诱导出的 $\|A\|_2$ 叫谱范数，它有一个非常漂亮的身份：

即 $A$ 的最大奇异值。为什么？因为 $\|Ax\|_2^2 = x^T A^T A x$，而 $A^T A$ 是对称半正定矩阵，其特征值就是奇异值的平方 $\sigma_i^2$；在特征向量方向上 $x^T A^T A x$ 取到最大，开根号就是 $\sigma_{max}$。

另一种常用矩阵范数是 Frobenius 范数，把矩阵当长向量：

注意它和谱范数的区别：Frobenius 范数等于奇异值平方和的平方根，谱范数等于最大奇异值——前者是”整体大小”，后者是”最强拉伸方向”。直觉上，$\|A\|_F$ 像矩阵的”总能量”，$\|A\|_2$ 像矩阵的”最大杠杆”。

铺垫：矩阵 = 线性变换的几何#

范数的几何解读依赖一个视角，这也是你数学档案里缺的一块：矩阵就是线性变换。$A$ 把一个向量”送到”另一个向量，几何效果只有五种基本形态的组合：拉伸、压缩、旋转、反射、剪切。下图把同一张带方向箭头的方块在 5 种矩阵作用下的样子画了出来（图源：Greg Gundersen 的博客文章 The Singular Value Decomposition (SVD) 的 Figure 2）：

从左到右五种变换分别对应：对角元 $> 1$ 的对角矩阵（拉伸）、对角元 $< 1$ 的对角矩阵（压缩）、正交矩阵（旋转）、带 $-1$ 特征值的正交矩阵（反射）、非对称的三角矩阵（剪切）。后面讲条件数、特征值、RoPE 旋转矩阵时，这套几何语言会反复用到。奇异值分解 $A = U\Sigma V^T$ 的几何含义就是：先旋转（$V^T$）、再沿坐标轴拉伸（$\Sigma$）、最后再旋转（$U$）——任何线性变换都只是”转一下、拉一下、再转一下”。拉伸的倍数正是奇异值 $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_n$。

几何直觉：圆变成椭圆#

条件数的几何意义一张图就够。下图（图源同上，Greg Gundersen 的 SVD 文章 Figure 5）左边是单位圆，右边是这个圆被某个矩阵 $A$ 映射后的像：

单位圆上的所有向量 $\|x\| = 1$ 被 $A$ 送到一个椭圆上，椭圆的两个半轴长度恰好是 $\sigma_{max}$ 和 $\sigma_{min}$，方向是左奇异向量 $u_1, u_2$。理由很直接：$x$ 沿右奇异向量 $v_1$ 时，$Ax = \sigma_1 u_1$，长度就是 $\sigma_1$。于是：

• 良态矩阵（$\kappa \approx 1$）：椭圆接近圆，所有方向的拉伸差不多，误差无论从哪个方向进来都不被放大；
• 病态矩阵（$\kappa \gg 1$）：椭圆被压成”雪茄”甚至”针”，$\sigma_{min}$ 方向几乎被压扁。误差一旦落进这个方向，解就会飞出去——Hilbert 矩阵就是这么回事。

配合条件数不等式：$\delta b$ 恰好落在 $u_{min}$（被压扁的方向）时，解的扰动 $\delta x = A^{-1}\delta b$ 落在 $v_{min}$ 方向，被放大了 $1/\sigma_{min}$ 倍。

条件数的几个性质#

• $\kappa(A) \ge 1$：因为 $\|A\|\cdot\|A^{-1}\| \ge \|AA^{-1}\| = \|I\| = 1$。
• 缩放不变：$\kappa(cA) = \kappa(A)$。所以”行列式小”不是病态的判据——把 $A$ 整体乘 1000，行列式变 $10^6$ 倍，条件数纹丝不动。判断病态必须看比率，不是看绝对大小。
• 正交矩阵 $\kappa = 1$：正交矩阵满足 $Q^T Q = I$，所有奇异值都是 1，圆映成圆。旋转/反射是最”温和”的矩阵，这也是为什么 RoPE 用旋转矩阵编码位置（第九节展开）。
• 有效精度：双精度下 $\epsilon \approx 2.2 \times 10^{-16}$，条件数 $\kappa = 10^k$ 意味着解要丢掉约 $k$ 位十进制有效数字。经验法则：$\kappa \cdot \epsilon \gtrsim 1$ 时，解基本是噪声。Hilbert 矩阵 $n=31$ 时 $\kappa \approx 10^{40}$，$\kappa \epsilon \approx 10^{24}$，所以图里那根红线横七竖八——不是求解器的错，问题本身就没有可解性。

顺带把上一篇欠的账补上：$\kappa$ 不仅衡量 $Ax = b$，对任意计算问题都能定义”问题条件数”。函数 $f$ 在输入 $x$ 处的条件数就是 $|f'(x)| \cdot |x|/|f(x)|$，衡量输入相对扰动引起的输出相对扰动。上一篇讲的”相近数相减”（catastrophic cancellation）本质上是减法 $f(a,b) = a - b$ 在 $a \approx b$ 处条件数趋于无穷；在线 softmax 的重缩放则是把一个病态写法重排成良态写法。本篇聚焦矩阵版本。

消元过程与计算量#

高斯消元你课程里学过，这里从”数值稳定性”的角度重新看一遍。目标是把增广矩阵 $[A \mid b]$ 通过行变换变成上三角 $[U \mid \tilde{b}]$，再回代求解。以 $3 \times 3$ 为例：

第一步：用第 1 行消去第 2、3 行的第 1 列。乘数（multiplier）分别是 $-4/2 = -2$ 和 $-2/2 = -1$，第 2 行减去 $-2$ 倍第 1 行（等价于加上 2 倍），第 3 行减去 $-1$ 倍第 1 行，得到

第二步：用第 2 行消第 3 行第 2 列，乘数 $8/(-4) = -2$，得到上三角 $U$。回代两步即可求出 $x$。

第 $k$ 步消元要处理剩下的 $n-k$ 行、每行约 $n-k+1$ 个元素，所以总乘加次数约为

每次乘加是 2 次浮点运算（乘+加），因此高斯消元的总计算量是 $\frac{2}{3} n^3$ flops；回代只有 $\frac{1}{2} n^2$，可以忽略。作为对比，矩阵乘法是 $2n^3$ flops——消元大约只有 GEMM 的三分之一，这是个”便宜”的算法。

致命弱点：除以小主元#

消元的关键操作是除以主元（对角线元素）。如果主元很小，会发生什么？用一个经典的教科书例子（几乎所有数值线性代数教材都会讲，比如 Trefethen & Bau 的《Numerical Linear Algebra》）：

精确解是 $x_1 \approx 0.9999$，$x_2 \approx 1.0001$。现在假设我们只有三位十进制浮点精度（$10^{-4}$ 量级机器精度，相当于”极寒版”的 FP16 环境，用来放大问题）。

不换行直接消：乘数 $m = 1 / 0.0001 = 10000$。第二行减去 $m$ 倍第一行：$(1 - 10000 \times 1) x_2 = 2 - 10000 \times 1$，即 $-9999 x_2 = -9998$，三位精度舍入为 $-1.00 \times 10^4 \cdot x_2 = -1.00 \times 10^4$，解得 $x_2 = 1.00$。回代：$0.0001 x_1 = 1 - 1.00 = 0$，于是 $x_1 = 0$。

结果：$x_1 = 0$，真解是 $x_1 \approx 1$——完全错了。 问题出在哪？第一步 $2 - 10000 \times 1 = -9998$ 这一步：$10000$ 和 $1$ 相差四个数量级，大数加小数把 $x_1$ 的信息直接舍入掉了（正是上一篇讲的”大数加小数丢精度”），于是 $x_1$ 的信息在消元第一步就彻底消失。

先换行再消（把第一行和第二行交换）：主元变成 $1$，乘数 $m = 0.0001$。第二行变为 $(1 - 0.0001 \times 1)x_2 = 1 - 0.0001 \times 2$，即 $0.9999 x_2 = 0.9998$，三位精度下 $x_2 = 1.00$，回代 $x_1 = 2 - 1.00 = 1.00$。

结果正确。 同样的机器、同样的精度，只是行交换了一下，结果从全错变全对。这就是部分主元（partial pivoting）的全部动机。

部分主元：把最大元换到对角线#

部分主元策略：消元到第 $k$ 步时，在 $A[k:n, k]$（第 $k$ 列主元位置往下的所有元素）中挑绝对值最大的，把它所在的行和当前第 $k$ 行交换，再以它为除数。这样保证所有乘数满足 $|m| \le 1$。

为什么”选最大的”？因为消元中真正危险的是除以很小的数：除数越小，舍入误差被放大的倍数越大。把主元选为列中绝对值最大者，就把每一步的放大倍数压到最小。而且行交换不改变方程组的解（只是重排方程顺序），换行是零成本、零风险的操作。

更严格的稳定性保证来自 Wilkinson 的经典分析：带部分主元的高斯消元是向后稳定的（第八节详细说），其误差界中的”增长因子”（growth factor）$\rho = \max_{i,j,k} |u_{ij}| / \max_{i,j} |a_{ij}|$ 在部分主元下通常很小（实践中 $\rho$ 很少超过 10），而不换主元时 $\rho$ 可以呈 $2^{n-1}$ 指数爆炸——理论最坏情况由 Wilkinson 构造的特殊矩阵达到，工程上无主元消元是绝对禁区。

还有更强的完全主元（complete pivoting）：每一步在右下角整个子矩阵里挑最大元。它理论上更稳，但每次要扫 $O(n^2)$ 个元素找最大值，总代价 $O(n^3)$ 的额外比较，几乎没人用——部分主元”够好且免费”，完全主元”略好但昂贵”，这是工程上经典的取舍。

从消元到分解#

设 $E_1, E_2, \dots, E_{n-1}$ 依次是各轮消元的行变换矩阵，消元完成后有

其中 $U$ 是上三角。把所有的 $E$ 移到右边：

关键事实：初等行变换矩阵的逆不需要算，直接读出来。$E_k$ 的逆就是把第 $k$ 轮乘数的符号反过来，于是 $(E_{n-1}\cdots E_1)^{-1}$ 恰好是把各轮乘数（变号后）填到对角线下方的下三角矩阵 $L$，且对角线全为 1：

验证一个 $2 \times 2$ 例子：$A = \begin{pmatrix} 2 & 1 \\ 8 & 7 \end{pmatrix}$，乘数 $m_{21} = 8/2 = 4$，消元得 $U = \begin{pmatrix} 2 & 1 \\ 0 & 3 \end{pmatrix}$，于是

为什么这个平凡的重述如此重要？因为分解一次，可以反复使用：

• 求解 $Ax = b$ 变成两步：先前代解 $Ly = b$（$O(n^2)$，下三角从顶往下算），再回代解 $Ux = y$（$O(n^2)$）；
• 换一个新的右端项 $b'$，只需要重做两次三角求解 $O(n^2)$，分解的 $2n^3/3$ 不用重算。

带部分主元时，行交换记成一个置换矩阵 $P$，分解写作 $PA = LU$，求解流程变为 $Ly = Pb$，$Ux = y$。这就是 LAPACK 例程 getrf（分解，输出 $L, U, P$）和 getrs（用分解结果求解）的由来——永远两步走，永远不显式求逆。

为什么不显式求逆：成本与稳定性的双重理由#

把”求逆再相乘”和”LU 分解”对比，这是本系列最重要的工程结论，也是 John D. Cook 的经典博客文章 Don’t invert that matrix 的标题——不要对矩阵求逆。

• 成本：显式求 $A^{-1}$ 等价于解 $n$ 个右端项 $A X = I$，成本 $O(n^3)$（约 $2n^3$ 到 $4n^3$ flops，取决于实现），而 LU 分解只要 $2n^3/3$，再解任意多个右端项每个只需 $O(n^2)$。对一个右端项来说，求逆贵 3 倍以上。
• 稳定性：求逆本身是”解 $n$ 个方程组”，每一步都会引入舍入误差；求完逆再乘 $b$，误差又被乘一遍。用 John Cook 的比喻：为了回答”这家餐厅有什么菜”，你把整本菜单背下来再逐页翻——直接把菜单翻到那一页（LU 求解）又快又稳。
• 结构：$A$ 如果稀疏、带状、对称正定，$L$ 和 $U$ 会继承部分结构（比如 Cholesky 的 $L$ 保持带状），而 $A^{-1}$ 几乎总是稠密的——求逆会破坏所有可以利用的结构。

数值上还有一个更微妙的点：$A^{-1}$ 的条件数等于 $A$ 的条件数（$\kappa(A^{-1}) = \kappa(A)$，因为奇异值互为倒数），所以”先求逆再乘”并不会让问题变得更好；而直接 LU 求解只需承担一次”向后误差 × 条件数”。结论是明确的：工程代码里出现 inv(A) @ b 或 A.inverse() 的地方，几乎都可以替换为分解求解，又稳又快。

与 GEMM 的联系：分块 LU#

GPU 上怎么把 LU 分解跑得快？答案是把消元分块化（blocked / panel factorization）：每一轮对外层的 $k$ 列（panel）做一次小规模消元，然后用矩阵乘更新剩余的大块子矩阵 $A_{22} \leftarrow A_{22} - L_{21} U_{12}$。这一步是标准的 GEMM——正好复用你在这两篇《GPU GEMM 优化完全拆解（一）》《GPU GEMM 优化完全拆解（二）》里学到的分块、共享内存、寄存器分块技术。LAPACK 的 getrf 带 DGETRF 的块化实现，cuSOLVER 的 cusolverDnDgetrf 是它的 GPU 版本——GPU 上 LU 的加速路径和 GEMM 是同一套打法。这也是为什么 ”$\frac{2}{3}n^3$ 的消元” 在硬件上能以接近 GEMM 的效率运行：大部分工作都被重构成了 GEMM。

1. RoPE 为什么数值友好：正交矩阵 κ = 1#

RoPE（旋转位置编码，论文 arXiv:2104.09864）把位置信息编码为对 $q$、$k$ 向量施加的旋转：2 维子空间上用 $\begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}$ 旋转。这类旋转矩阵是正交矩阵：$Q^T Q = I$，所有奇异值为 1，$\kappa = 1$。

这意味着两件事：其一，旋转保长度（$\|Qx\| = \|x\|$），注意力点积的模长不受位置编码破坏；其二，误差不会被旋转放大——即使 FP16/BF16 下 $q$、$k$ 的表示有舍入误差，旋转也不会把它们放大 $\kappa$ 倍。用前面的话说，RoPE 选择了一种”条件数最优”的线性变换。对比之下，如果位置编码用非正交矩阵（比如某些线性外推方案），数值风险是隐性的。

2. 在线 softmax 的重缩放：把病态重排成良态#

上一篇推导过在线 softmax：为防溢出，减去当前最大值 $m$，用 $\exp(x_i - m)$ 代替 $\exp(x_i)$，再乘回重缩放因子。用本篇的语言：直接算 $\exp$ 求和是病态写法（大数加小数），减去最大值把中间量压到 $[0,1]$，等价于把一个 $\kappa$ 巨大的求和重排成 $\kappa \approx 1$ 的良态求和——同一个数学对象，不同的计算顺序，条件数天差地别。FlashAttention 的数值重缩放（论文 arXiv:2205.14135）是同一思想的流式版本。这是”算法稳定性”最重要的工程模板：不改问题，只改写法。

3. 不显式求逆：MLA 吸收、LoRA 合并、量化#

你的方向上”求逆”出现的方式，全是被刻意回避的：

• MLA 的矩阵吸收（详见《MLA 完全拆解》）：把 $W^{UK}$ 吸收进 $W^Q$、$W^U$ 吸收进 $W^O$，是把小矩阵乘到大权重里，是纯乘法重组，不涉及任何求逆；
• LoRA 合并：$W' = W_0 + BA$，加法；
• 量化/反量化：scale 因子是标量，$x / s$ 等价于乘 $1/s$——1×1 的”对角矩阵求逆”是平凡且 $\kappa = 1$ 的。

如果哪一天你在推理管线里看到显式 inv()，几乎都可以问一句：这里是不是可以用分解或直接重组替代？这是”Don’t invert that matrix”精神在 LLM 部署中的版本。

4. 奇异值谱与低秩压缩：σ_min 决定”可逆性”#

SVD 视角下，权重矩阵的低秩可压缩性来自奇异值谱衰减：大模型权重矩阵的奇异值通常快速衰减，$\sigma_{min}$ 接近零，因此矩阵”数值上不可逆”——但这恰恰是好事，说明它可以被低秩近似（截断 SVD 砍掉小奇异值方向）。MLA 的低秩 KV 压缩、SVD-LLM（arXiv:2403.17414，用 SVD 做 LLM 权重压缩）、XFP 那类码本量化，数学上都是同一个事实的利用：$\sigma_{min}$ 小 → 精确求逆无意义 → 但低秩/伪逆有意义。伪逆 $A^{+} = V\Sigma^{+}U^T$（$\Sigma^{+}$ 对非零奇异值取倒数）是”在秩亏矩阵上最有意义的求逆”——它恰好跳过 $\sigma \approx 0$ 的方向，等价于数值上最稳的逆。这也是为什么数值库提供 gesvd 而不是”安全求逆”。

5. 状态空间模型与线性注意力：特征值、分块消元与谱半径#

最后两个预告性的连接，细节留给后面的文章：

• S4 论文（arXiv:2111.00396）明确讨论了数值问题：状态空间模型的 $A$ 矩阵直接对角化在长序列下数值不稳，原因是特征向量矩阵可能病态（$\kappa$ 大），于是改用正规矩阵（正规矩阵可被正交对角化，特征向量矩阵 $\kappa = 1$）。这是”特征值问题也会病态”在 LLM 系统的真实案例，也是下一篇系列文章（特征值与谱）的引子。Mamba 的 ZOH 离散化要计算 $(\Delta A)^{-1}(e^{\Delta A} - I) \Delta B$——这个”除以 $\Delta A$“的小型矩阵运算在实现中同样用级数展开回避显式求逆。
• RetNet 的分块并行（arXiv:2307.08621）在形式上与分块消元同构：序列被切成块，块内并行、块间递推，需要计算的块矩阵组合与分块高斯消元中的 Schur 补结构一致。分块、消元、递推这三件事在数值线性代数里本来就是一家。