最小二乘把条件数平方了
线性回归是这个星球上被执行次数最多的数值算法。它的教科书解法叫「正规方程」，形式优美、推导简单、三行就能实现。它也有一个几乎从不被提起的性质：它把问题的条件数平方了。这意味着你手上 16 位的精度，一进门就只剩 8 位——而这件事在数据有共线性时会立刻要命。
一个 3×2 的最小二乘问题（这是数值分析里的标准测试例，叫 Läuchli 矩阵）：
⎡ 1 1 ⎤ ⎡ 1 ⎤
A = ⎢ ε 0 ⎥ b = ⎢ 0 ⎥ ε = 10⁻⁸
⎣ 0 ε ⎦ ⎣ 0 ⎦
真解（精确算出来的）：x = [0.5, 0.5]
κ(A) = 1.414 × 10⁸ —— 有点病态，但离 double 的极限 10¹⁶ 还差八个数量级
用教科书的正规方程解：x = solve(AᵀA, Aᵀb)。
问：会得到什么？
第二问：改用 QR 分解解同一个问题呢？
把 AᵀA 写出来
AᵀA = ⎡ 1 + ε² 1 ⎤ ε = 10⁻⁸，所以 ε² = 10⁻¹⁶
⎣ 1 1 + ε² ⎦
而 1 + 10⁻¹⁶ 在 double 里是多少？
ulp(1) = 2.22 × 10⁻¹⁶
10⁻¹⁶ 比半格（1.11 × 10⁻¹⁶）还小
⇒ ⟨1 + 10⁻¹⁶ = 1⟩
所以 double 里的 AᵀA 是：
⎡ 1 1 ⎤
⎣ 1 1 ⎦ ← 两行完全相同，奇异矩阵
numpy 抛 LinAlgError: Singular matrix。答案是 D。
ε 出现在 A 里的时候好好的（10⁻⁸ 离 double 的极限远得很）。一进 AᵀA 就被平方成 10⁻¹⁶，正好卡在机器 epsilon 上，整个消失。
第二问：QR 分解给出 [0.5, 0.5]，和精确解逐位相同。numpy.linalg.lstsq（内部走 SVD）也给出 [0.5, 0.5]。
κ(AᵀA) = κ(A)²，这是个定理
把 ε 扫一遍，这条关系会以最直白的方式显示出来：
| ε | κ(A) | κ(AᵀA) | 正规方程解 | QR 解 |
|---|---|---|---|---|
| 10⁻² | 1.414 × 10² | 2.000 × 10⁴ | [0.4999750012, 0.4999750013] | [0.4999750012, 0.4999750012] |
| 10⁻⁴ | 1.414 × 10⁴ | 2.000 × 10⁸ | [0.4999999978, 0.4999999972] | [0.4999999975, 0.4999999975] |
| 10⁻⁶ | 1.414 × 10⁶ | 2.000 × 10¹² | [0.5000000000, 0.5000000000] | [0.5000000000, 0.5000000000] |
| 10⁻⁸ | 1.414 × 10⁸ | 5.962 × 10¹⁶ | 奇异矩阵，解不出来 | [0.5000000000, 0.5000000000] |
κ(AᵀA) 那一列，每一行都是左边那一列的平方：(1.414×10²)² = 2×10⁴，(1.414×10⁴)² = 2×10⁸，(1.414×10⁶)² = 2×10¹²。
这不是巧合，是一行代数：A 的奇异值是 σᵢ，那么 AᵀA 的特征值就是 σᵢ²，于是最大比最小的比值被平方。
正规方程要求 κ(A)² < 1/eps，也就是 κ(A) < 10⁸。QR 只要求 κ(A) < 10¹⁶。
换一种说法，也是更该记住的说法：用正规方程，等于把你的 double 当 float 用。float32 的 eps 是 1.19×10⁻⁷，正好是 double 开根号的量级。
这条不等式是这一章的全部内容。它有一个很实用的推论：如果你的特征之间有强相关（共线性），κ(A) 会轻松突破 10⁸——而那正是正规方程失效的地方。共线性和病态是同一件事的两个名字。
那为什么正规方程还到处都是
因为它有真实的优点，不是纯粹的历史包袱：
| 正规方程 | QR | SVD | |
|---|---|---|---|
| 运算量（m×n，m≫n） | mn² + n³/3 | 2mn² − 2n³/3 | 约 2mn² + 11n³ |
| 大致快慢 | 最快（约快一倍） | 中 | 最慢 |
| 能承受的 κ(A) | 10⁸ | 10¹⁶ | 10¹⁶ |
| 能处理秩亏 | 不能 | 要带列主元才行 | 能，而且能告诉你秩是多少 |
| 能流式／增量更新 | 能（只要累加 AᵀA 和 Aᵀb） | 要用 Givens 旋转 | 难 |
最后一行是它在工程里活下来的真正原因：AᵀA 只有 n×n，和样本数无关。一亿行数据、二十个特征，你只需要累加出一个 20×20 的矩阵，数据扫一遍就能丢掉。这在分布式和流式场景里是不可替代的（每个分片各自累加，最后把矩阵加起来就行）。
所以正确的态度不是「永远别用」，而是：用之前先看 κ(A)，或者至少先做标准化。
三件几乎免费的补救
- 中心化 + 标准化。每一列减均值、除标准差。这一步能把「因为量纲不同而虚高的条件数」全部消掉——一列是「年龄（20–80）」、一列是「收入（10⁴–10⁶）」，不标准化的话 κ 里有四个数量级纯粹来自单位。统计软件默认做这件事，理由就是这个。
- 正则化（岭回归）。
(AᵀA + λI)x = Aᵀb。加上λI把最小的特征值从σₘᵢₙ²抬到σₘᵢₙ² + λ，条件数应声下降。岭回归的数值意义就是买条件数——它的统计意义（缩减方差）和数值意义（降 κ）是同一件事的两面，这不是巧合。 - 直接调
lstsq。numpy、scipy、R 的lm()、statsmodels 内部走的都是 QR 或 SVD，不是正规方程。你只要不自己写solve(A.T@A, A.T@b)，就已经在正确的路上了。
共线性（collinearity）、病态（ill-conditioning）、秩亏（rank deficiency）——统计课和数值课用不同的词说同一件事：
- 统计说「两个特征高度相关，系数不稳定，VIF 很大」。
- 数值说「设计矩阵的两列几乎平行，最小奇异值接近 0，κ 很大」。
- 秩亏是极限情况：完全平行，σₘᵢₙ = 0，κ = ∞。
方差膨胀因子 VIF 和条件数是可以互相换算的（都在量同一件事）。知道这一点之后，统计里那些「处理共线性」的手段（删特征、PCA、岭回归）和数值里那些「降条件数」的手段（换基、正则化、SVD 截断）就会显示出它们本来就是同一批工具。
Läuchli 矩阵只要三行就能构造，值得亲手看一次「1 + ε² 变成 1」：
import numpy as np
print('%-8s %-12s %-12s %-24s %s' % ('ε', 'κ(A)', 'κ(AᵀA)', '正规方程解', 'QR 解'))
for e in (1e-2, 1e-4, 1e-6, 1e-8):
A = np.array([[1.0, 1.0], [e, 0.0], [0.0, e]])
b = np.array([1.0, 0.0, 0.0])
AtA = A.T @ A
try:
xn = np.linalg.solve(AtA, A.T @ b)
s = '[%.10f %.10f]' % (xn[0], xn[1])
except np.linalg.LinAlgError:
s = '奇异矩阵，解不出来'
q, r = np.linalg.qr(A)
xq = np.linalg.solve(r, q.T @ b)
print('%-8.0e %-12.3e %-12.3e %-24s [%.10f %.10f]'
% (e, np.linalg.cond(A), np.linalg.cond(AtA), s, xq[0], xq[1]))
# 看清楚 AᵀA 变成了什么
A = np.array([[1.0, 1.0], [1e-8, 0.0], [0.0, 1e-8]])
print('AᵀA =', (A.T @ A).tolist()) # [[1.0, 1.0], [1.0, 1.0]]
print('1 + 1e-16 ==', repr(1 + 1e-16)) # 1.0
print('真解恒为 [0.5, 0.5]')
会打出：
ε κ(A) κ(AᵀA) 正规方程解 QR 解 1e-02 1.414e+02 2.000e+04 [0.4999750012 0.4999750013] [0.4999750012 0.4999750012] 1e-04 1.414e+04 2.000e+08 [0.4999999978 0.4999999972] [0.4999999975 0.4999999975] 1e-06 1.414e+06 2.000e+12 [0.5000000000 0.5000000000] [0.5000000000 0.5000000000] 1e-08 1.414e+08 5.962e+16 奇异矩阵，解不出来 [0.5000000000 0.5000000000]
注意 ε=10⁻⁸ 那一行的 κ(AᵀA) 是 5.962×10¹⁶，而不是理论值 2×10¹⁶——因为这时候连「计算 κ(AᵀA)」这件事本身都不可靠了。条件数很大时，条件数的估计值本身也不准。
python3 lauchli.py
在线：Google Colab。
线性回归。R 的 lm()、Python 的 statsmodels、scikit-learn 的 LinearRegression——全部走 QR 或 SVD（sklearn 用的是 scipy.linalg.lstsq，也就是 SVD）。你如果自己写 np.linalg.inv(X.T@X) @ X.T@y，就同时踩了第 15 章和第 16 章两个坑。这行代码在教程和面试题里出现的频率高得惊人。
三维重建与 SLAM。Bundle adjustment 解的是一个巨大的非线性最小二乘，每次迭代都要解一个正规方程（那个 JᵀJ 就是 AᵀA）。这里不得不用正规方程，因为 JᵀJ 的稀疏结构（舒尔补）是唯一能让问题算得动的东西。代价是必须靠 Levenberg–Marquardt 的阻尼项（也就是那个 λ）来把条件数压住——LM 算法的阻尼项在数值上就是岭回归。
卡尔曼滤波（又是它）。协方差矩阵 P 本质上是一个 AᵀA，所以它天生条件数平方。平方根滤波之所以叫这个名字，就是因为它维护的是 P 的 Cholesky 因子——相当于一直待在 QR 那一侧，从不形成 AᵀA。
深度学习里的二阶方法。K-FAC、Shampoo 这些方法要估计并求逆 Fisher 信息矩阵（也是个 AᵀA）。它们全都带阻尼项，理由完全相同。
「κ(A) = 10⁸ 不算大，double 有 16 位，还剩 8 位，够用。」
「还剩 8 位」这个算法只在你直接解 A 时成立。走正规方程，实际面对的条件数是 κ(A)² = 10¹⁶——剩 0 位。
这条直觉错在把「问题的条件数」和「我实际求解的那个问题的条件数」当成了同一个。你选的解法会创造出一个新的问题，而新问题有自己的条件数。这是这本书的一条暗线：第 8 章「先收后支再相减」把条件数从 1 变成 8.4×10¹⁰，第 15 章的求逆丢掉后向稳定，这一章的正规方程把 κ 平方——三次都是「解法制造了新的病态」。
可操作的判据：算完 κ(A) 之后，问一句「我的解法会不会动这个数」。会动，就用动完之后那个数去判断精度够不够。
正确答案是 D：numpy 抛 LinAlgError: Singular matrix，因为 AᵀA 在 double 里变成了 [[1,1],[1,1]]。QR 给出 [0.5, 0.5]，精确。
A 做 QR，这句话完全正确（QR 给出的答案是精确的）。它错在没有把解法本身算进去。κ(A) 描述的是问题，κ(AᵀA) 描述的是你实际在解的那个东西。
B 「掉几位而已」——在 ε = 10⁻⁴ 那一行，B 是对的（正规方程给出 [0.4999999978, 0.4999999972]，掉了几位但可用）。ε 只再小四个数量级，就从「掉几位」跳到「解不出来」——因为平方让退化速度也翻了倍。
C 「一个离谱但合法的数」——这是最危险的那种结果，而这道题上你运气好，没碰到。AᵀA 恰好退化成两行完全相同，numpy 才检测得出来。如果 ε 是 10⁻⁷.⁵ 这类值，AᵀA 会是「几乎奇异但不完全奇异」——numpy 不会报错，会安静地返回一个错得离谱的数。C 描述的是更常见的情形，只是这个特定例子恰好跳过了它。
x = np.linalg.inv(A.T @ A) @ A.T @ b —— 第 15 章 + 第 16 章，两个坑一起踩
x = np.linalg.lstsq(A, b, rcond=None)[0] —— 走 SVD，能承受 κ(A) 到 10¹⁶，还能处理秩亏
或者 x = scipy.linalg.lstsq(A, b)[0]；要快就 lapack_driver='gelsy'（带列主元 QR）
数据没标准化就直接拟合
先中心化 + 标准化。免费，而且往往能把 κ 降好几个数量级——因为大部分条件数来自量纲，不是来自数据结构
发现结果不稳定，于是加大迭代次数／换优化器
先算 np.linalg.cond(A)。如果它是 10¹⁰，那不是优化器的问题，加多少迭代都没用——要么加正则化，要么删特征，要么承认这几个系数就是估不准
非要用正规方程（因为要流式累加、或者要利用稀疏结构）时，至少做两件事：① 累加 AᵀA 时用更高精度或补偿求和（第 8 章）；② 加一个小的 λI。这两条加起来还是比 QR 便宜，而且把最坏情况堵住了。
这一章的一句话
正规方程的推导只有三行，代价却是把条件数平方——它等于悄悄把你的 double 降级成 float，而这份降级从不出现在任何报错里。
卷 IV 到此结束。前十六章的全部内容，都发生在同一个格子里：IEEE 754 双精度。
卷 V 换格子。接下来四章会看到：钱为什么必须离开这个格子（第 17 章）；AI 为什么主动搬到一个只有 7 位尾数的格子里去，还训得更好（第 18 章）；以及一个让所有人都不舒服的事实——就算格子不变，同一份源代码在同一台机器上，只要加一个编译开关，答案就会变。第 19 章会在这台机器上用真的 C 编译器把它跑出来。