永远不要求逆矩阵
线性代数课上，方程组的解写作 x = A⁻¹b。于是几乎每个人第一次写代码时都会照着抄一遍——inv(A) @ b。所有的数值线性代数教材都会在某一页写「不要这样做」，但很少解释清楚为什么。这一章给三个理由，每个都有实测数字；顺便讲一个更基础的东西：换两行的位置，能把一个答案从「全错」救成「全对」。
解这个 2×2 方程组：
10⁻²⁰ · x₁ + 1 · x₂ = 1
1 · x₁ + 1 · x₂ = 2
真解：x₁ ≈ 1，x₂ ≈ 1（精确到 20 位以内都是这样）
用最标准的高斯消元：拿第一行去消第二行的首列。不做行交换，就按矩阵原来的样子做——这是每本线性代数教材的第一个算法。
问：算出来的 [x₁, x₂] 是什么？
第二问：numpy.linalg.solve 解这个方程组，会给出什么？
把消元过程写出来
用第一行消第二行：乘数 f = 1 / 10⁻²⁰ = 10²⁰
第二行变成：
a₂₂' = 1 − 10²⁰ × 1 = −10²⁰ + 1
= −100000000000000000000 ← ⟨那个 1 被抹掉了⟩
ulp(1e20) = 16384，1 连一格都不到
b₂' = 2 − 10²⁰ × 1 = −10²⁰ + 2
= −100000000000000000000 ← ⟨那个 2 也被抹掉了⟩
回代：
x₂ = b₂' / a₂₂' = 1 ← 这个恰好对
x₁ = (1 − x₂) / 10⁻²⁰ = 0 / 10⁻²⁰ = 0 ← ✗ 应该是 1
结果 [0, 1]。答案是 C。
注意这里发生了什么：第二个方程 x₁ + x₂ = 2 里的信息被完全抹掉了。消元之后那一行变成了 −10²⁰·x₂ = −10²⁰，也就是 x₂ = 1——这句话里已经没有 x₁ 的任何痕迹。而 x₁ 只能靠第一行回代求出来，那一行的系数是 10⁻²⁰，把误差放大了 10²⁰ 倍。
修法只有一步：换行。
先把绝对值最大的那一行换到上面（这叫「部分主元」）：
1 · x₁ + 1 · x₂ = 2
10⁻²⁰ · x₁ + 1 · x₂ = 1
乘数 f = 10⁻²⁰ / 1 = 10⁻²⁰（很小，安全）
第二行变成：
a₂₂' = 1 − 10⁻²⁰ × 1 = 1 ← 抹掉的是 10⁻²⁰，无所谓
b₂' = 1 − 10⁻²⁰ × 2 = 1
回代：x₂ = 1，x₁ = (2 − 1) / 1 = 1 ✓ [1, 1]
没有换更高精度，没有换算法，只是把两行调了个位置。
选主元的作用是：让所有的乘数都不超过 1。
消元的每一步都是「第 i 行减去 f 倍的主元行」。如果 |f| > 1，这一步就在放大——被减数的信息会被淹没在放大后的主元行里。f = 10²⁰ 意味着放大一百万亿亿倍，那一行原有的信息一点不剩。
把绝对值最大的那一行选作主元，保证 |f| ≤ 1，于是每一步最多「保持」，不会放大。这是把一个可以任意不稳定的算法，变成后向稳定算法的全部代价：一次比较和一次行交换。
第二问的答案：numpy.linalg.solve 给出 [1, 1]。LAPACK 的 dgesv 内部就是带部分主元的 LU 分解——你调库的时候已经在享受这一节了。
那为什么不能求逆
回到标题。x = inv(A) @ b 和 x = solve(A, b) 在数学上完全一样。实测三个理由：
理由一：慢一倍多
# n = 2000 的随机矩阵，本机实测 numpy.linalg.solve(A, b) 0.037 秒 numpy.linalg.inv(A) @ b 0.079 秒 ★ 慢 2.14 倍
算术上也说得通：LU 分解是 n³/3 次运算，求逆是 n³ 次（相当于解 n 个方程组）——三倍的工作量，只为了得到一个你根本不需要的中间结果。
理由二：残差差十一个数量级
这个才是真正的理由。前面十四章一直在说「后向稳定是算法唯一能给的承诺」——而求逆把这个承诺弄丢了：
| 希尔伯特矩阵 | solve 的残差 | inv @ b 的残差 | 相差 |
|---|---|---|---|
| n = 10（κ = 1.6×10¹³） | 2.220 × 10⁻¹⁶ | 4.908 × 10⁻⁵ | 10¹¹ 倍 |
| n = 12（κ = 1.6×10¹⁶） | 4.441 × 10⁻¹⁶ | 6.847 × 10⁻² | 10¹⁴ 倍 |
看第二行：solve 的残差是 4.44 × 10⁻¹⁶（机器精度，后向稳定），inv @ b 的残差是 0.068。
这不是「稍微不准一点」。这意味着：用 inv 求出来的解，连「它是某个邻近问题的精确解」都不成立了。第 6 章那条 前向 ≲ κ × 后向 的公式，在这里给出的上界从 κ×10⁻¹⁶ 变成 κ×10⁻²——凭空多了十四个数量级的不确定。
原因不难理解：inv(A) 的每一列都是解一个方程组得到的，各自带着自己的误差；这些误差在后面那个矩阵—向量乘法里被混在一起，再也没有任何一个「原方程」能把它们约束住。而 solve 从头到尾只面对一个方程组，误差始终被那个方程组管着。
理由三：它会把结构毁掉
# 一个 8×8 的三对角矩阵（只有主对角和上下各一条） 原矩阵的非零元素：22 / 64 ← 稀疏，可以 O(n) 求解 它的逆的非零元素：64 / 64 ← ★ 完全稠密
三对角、带状、稀疏、对称正定——真实世界里的大矩阵几乎都有结构，而结构可以把 O(n³) 降到 O(n) 或 O(n log n)。求逆会毫无例外地把所有结构毁掉：稀疏矩阵的逆一般是稠密的。
一个 10⁶ × 10⁶ 的稀疏矩阵（有限元、图拉普拉斯、推荐系统里都是这个规模），存下来可能只要几百 MB；它的逆有 10¹² 个元素，8 TB。
LU 分解是这一节所有事情的实际执行者，值得说清它是什么：
带部分主元的高斯消元，本质上是把 A 拆成 PA = LU——P 是行交换（选主元），L 是下三角，U 是上三角。拆完之后，解方程只要两次三角回代，各 O(n²)。
关键的实用价值在这里：分解一次，可以解无数个不同的右端项。如果你要对同一个 A 解一百个 b，正确的做法是分解一次（n³/3）+ 一百次回代（100n²），而不是分解一百次，更不是求一次逆。
这也是「求逆」这个念头最常见的来源：「我要解很多次，所以先把逆算出来存着」。正确的对应物不是逆，是存下 LU 分解（scipy.linalg.lu_factor / cho_factor）。同样是「预计算一次」，但保住了稳定性和结构。
按结构选求解器
知道矩阵的结构，能省的不只是时间，还有精度：
| 矩阵是什么 | 该用什么 | 代价 |
|---|---|---|
| 一般方阵 | LU + 部分主元（solve） | n³/3 |
| 对称正定 | Cholesky（cho_factor） | n³/6，而且不用选主元就稳定 |
| 三对角 / 带状 | Thomas 算法 / solve_banded | O(n) |
| 稀疏 | scipy.sparse.linalg.spsolve | 取决于稀疏结构 |
| 超大 / 只有矩阵—向量乘法 | 迭代法（CG、GMRES） | 每步一次乘法 |
| 长方形（最小二乘） | QR 或 SVD，绝不是正规方程 | 下一章 |
「对称正定不用选主元」这一条值得记：正定性本身就保证了乘数不会爆炸，所以 Cholesky 天生后向稳定，还省一半运算。协方差矩阵、图拉普拉斯、核矩阵、二阶优化里的 Hessian（正定时）——全都属于这一档。
三个实验，一个个看：
import numpy as np, time
# ① 不选主元 vs 选主元
eps = 1e-20
f = 1 / eps
a22 = 1 - f * 1; b2 = 2 - f * 1
x2 = b2 / a22; x1 = (1 - x2) / eps
print('不选主元 =', [x1, x2]) # [0.0, 1.0]
print('numpy =', np.linalg.solve([[eps, 1.], [1., 1.]], [1., 2.])) # [1. 1.]
# ② solve vs inv：速度与残差
rng = np.random.default_rng(7); n = 2000
A = rng.standard_normal((n, n)); b = A @ np.ones(n)
t = time.perf_counter(); np.linalg.solve(A, b); print('solve %.3fs' % (time.perf_counter()-t))
t = time.perf_counter(); np.linalg.inv(A) @ b; print('inv@b %.3fs' % (time.perf_counter()-t))
for n in (10, 12):
H = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)])
b = H @ np.ones(n)
x1 = np.linalg.solve(H, b); x2 = np.linalg.inv(H) @ b
print(f'Hilbert {n}: solve 残差 {np.abs(b-H@x1).max():.3e} '
f'inv@b 残差 {np.abs(b-H@x2).max():.3e}')
# ③ 求逆毁结构
T = (np.diag(np.full(8, 2.0)) + np.diag(np.full(7, -1.0), 1)
+ np.diag(np.full(7, -1.0), -1))
print('三对角非零元', int((T != 0).sum()), '/ 64 '
'它的逆非零元', int((np.abs(np.linalg.inv(T)) > 1e-12).sum()), '/ 64')
python3 solve.py # n=2000 的那一步大约几十毫秒
在线：Google Colab。把 n 调到 4000 会看到差距拉得更开（求逆是 O(n³)，常数更大）。
卡尔曼滤波。教科书上的更新式里赫然写着一个矩阵求逆。而所有严肃的实现都不这么做——它们用 Joseph 形式（保证协方差矩阵保持对称正定）或者平方根滤波（直接维护协方差的 Cholesky 因子，永远不形成协方差本身）。阿波罗登月的导航计算机用的就是平方根滤波，理由正是这一章：机载计算机精度低，直接求逆会让协方差矩阵失去正定性，滤波器发散。
高斯过程。预测公式里有 K⁻¹y，而实现里一律是 cho_solve(cho_factor(K + σ²I), y)。那个 + σ²I（jitter）本身也是数值手段——把最小特征值抬起来，降条件数，和第 7 章的正则化是同一招。
图形学的矩阵求逆是个例外。4×4 的变换矩阵求逆是标准操作，因为矩阵小、结构已知（旋转部分正交、平移单独一列），而且逆本身要被复用几百万次。「不要求逆」是一条针对大矩阵和一次性求解的建议，不是绝对禁令。
如果你真的需要逆本身（比如要报告参数的协方差矩阵、算标准误），那当然就得算。但要意识到你是在为「看它」付钱，而不是为「解方程」付钱。
「x = A⁻¹b 是线性代数的定义，照着写总没错。要是慢，那也是常数因子的事。」
数学定义描述的是结果，不是路径——这是第 9 章那句话的又一次应用。A⁻¹b 定义了 x 是什么，它没有说要先把 A⁻¹ 造出来。
而代价不只是常数因子。后向稳定性丢了：希尔伯特 12 上残差从 4.4×10⁻¹⁶ 变成 6.8×10⁻²，这不是慢，这是错。
一个好用的翻译规则：看到 A⁻¹B，读成「解 AX = B」；看到 BA⁻¹，读成「解 XA = B」。数学式子里的 ⁻¹ 几乎总是「解一个方程」的简写，而不是「造一个矩阵」的指令。numpy 里对应的写法是 np.linalg.solve(A, B)；MATLAB 干脆把它做成了运算符 A \ b。
正确答案是 C：[0, 1]，第一个分量整个丢了。numpy.linalg.solve 给出 [1, 1]，因为它内部做了行交换。
1 − 10²⁰ 里那个 1 不是被削弱了，是完全没进入结果（因为 ulp(10²⁰) = 16384）。信息要么在，要么不在——很多时候没有中间状态。
D 「NaN 或 Infinity」——除以 10⁻²⁰ 完全合法（结果是 10²⁰，离 double 上限 1.8×10³⁰⁸ 还远得很）。这一章的失败又是一次「安静的失败」：没有异常、没有警告，返回一个格式完全正确的 [0, 1]。这已经是这本书第五次出现同一种失败形态了——它们从不报错，这就是为什么要靠事前分析。
x = np.linalg.inv(A) @ b
x = np.linalg.solve(A, b) 快两倍，残差差十一个数量级
Ainv = np.linalg.inv(A)，然后循环里 Ainv @ b_i
lu = scipy.linalg.lu_factor(A)，然后循环里 lu_solve(lu, b_i) —— 同样是「预计算一次」，但保住稳定性
solve(A, b)，而 A 是协方差矩阵
cho_solve(cho_factor(A), b) —— 对称正定专用，运算量减半，不用选主元
自己写高斯消元（作业、面试、嵌入式里没有 BLAS 的场合）
一定要加选主元。三行代码：找出这一列绝对值最大的行、交换、记住交换记录。不加这三行，你的实现在某些输入上会给出零位有效数字。
还有一条元层面的建议：看到数学式子里的 ⁻¹，先翻译成「解一个方程」再动手。这个翻译动作能挡掉这一章的绝大部分问题，也能挡掉下一章的全部问题。
这一章的一句话
矩阵求逆是一个数学记号，不是一个计算步骤；把它照字面执行，你会多花两倍时间、丢掉后向稳定性、并且毁掉矩阵的全部结构。
下一章是这一招在数据分析里最常见的形态。用最小二乘拟合时，教科书给的是「正规方程」AᵀA x = Aᵀb——它看起来很干净，而且到处都在用。可是它有一个致命的性质：κ(AᵀA) = κ(A)²。条件数被平方了。这意味着：一个只需要 8 位精度的问题，用正规方程去解就要 16 位；而只要 κ(A) 超过 10⁸，double 就已经不够了。下一章会给出一个 3×2 的矩阵，numpy 在它上面直接拒绝求解。