06 好算法不给准答案

▷ 猜一下差多少

要算 e⁻⁵⁰。手边没有 exp，就用泰勒级数，教科书第一页那个：

e^x = 1 + x + x²/2! + x³/3! + x⁴/4! + …

代 x = −50，一路加到项小得可以忽略为止（169 项）

这个级数在全实轴上收敛，收敛得还很快，每一项都用 double 算，误差不超过半个 ulp。

问：加完 169 项，得到什么？（真值是 1.93 × 10⁻²²）

A 1.9287498479639178e-22。十几位全对，级数收敛得这么快，没理由不对 B 大约 1.9287e-22，对前四五位 C 量级还在 10⁻²² 附近，但一位有效数字都不对 D 得到 2041.8329628976246

顺便想想：为什么会这样？级数每一项都没算错啊。

两种「误差」，问的是两件事

假设你想算 y = f(x)，实际算出来的是 ŷ。有两种问法：

叫什么	问什么	怎么写
前向误差	我的答案离真答案多远？ ——这是你真正关心的	\|ŷ − y\| / \|y\|
后向误差	我的答案，精确地是哪个输入的答案？那个输入离我的真输入多远？ ——这是你能验证的	min \|x̃ − x\| / \|x\|，使得 f(x̃) = ŷ

后向误差这个提法第一次见到都会觉得别扭：我要的是答案，你却跟我讨论输入？

它之所以有用，是因为它把一件事拆干净了：

◆ 主线

前向误差 ≲ 条件数 × 后向误差

这条式子把「答案错了」这件事分成了责任明确的两块：

后向误差是算法的成绩单。它衡量「这个算法有多老实」。好算法的后向误差在 eps 量级（10⁻¹⁶），跟问题多难无关。
条件数是问题的成绩单，上一章讲过，你改不了。

一个算法如果它的后向误差总在 eps 量级，就叫后向稳定（backward stable）。这是数值分析里对「好算法」的正式定义——注意它完全不提答案准不准。

为什么要这么定义？因为这是唯一能做到的承诺。问题病态的时候，谁都给不出准答案（上一章）。但「我给出的答案是一个和你输入几乎一样的问题的精确答案」这件事，是可以做到、可以证明、也可以检验的。

还有一层更实际的理由：你的输入本来就有误差。测量数据、上一步的计算结果、存进 double 时的抹零——输入早就不是「真输入」了。既然如此，一个算法只要把误差控制在「和输入自带的不确定度同一量级」，它就已经做到了它能做的一切。再往下就是问题的事了。

把这条公式验一遍

拿一个真会病态的例子：希尔伯特矩阵。它的第 i 行第 j 列是 1/(i+j−1)，长这样：

        ⎡ 1    1/2  1/3  1/4 ⎤
H₄ =    ⎢ 1/2  1/3  1/4  1/5 ⎥
        ⎢ 1/3  1/4  1/5  1/6 ⎥
        ⎣ 1/4  1/5  1/6  1/7 ⎦

看起来无比温和，实际上是数值分析界最有名的病态矩阵。

做法：取真解 x = (1, 1, …, 1)，算出 b = Hx，然后假装不知道 x，用 numpy.linalg.solve（内部是 LAPACK 的带部分主元高斯消元）把它解回来。三个量都能量出来：

n	κ(H)	后向误差 β	前向误差	κ × β
6	1.495 × 10⁷	4.53 × 10⁻¹⁷	4.41 × 10⁻¹⁰	6.78 × 10⁻¹⁰	✓ 前向 < κβ
8	1.526 × 10¹⁰	1.63 × 10⁻¹⁶	4.84 × 10⁻⁷	2.49 × 10⁻⁶	✓
10	1.602 × 10¹³	7.58 × 10⁻¹⁷	4.20 × 10⁻⁴	1.21 × 10⁻³	✓
12	1.611 × 10¹⁶	9.78 × 10⁻¹⁷	5.42 × 10⁻¹	1.58	✓ 答案已经 54% 是错的
14	3.002 × 10¹⁷	2.25 × 10⁻¹⁷	5.65	6.76	✓ 答案比真值还大好几倍

看这张表最该看的是后向误差那一列：从 n=6 到 n=14，它一直待在 10⁻¹⁷ 到 10⁻¹⁶，纹丝不动。矩阵变得再病态，LAPACK 的成绩单一直是满分。

而前向误差从 4×10⁻¹⁰ 一路爬到 5.65（也就是答案比真值差了五倍多）。这份恶化百分之百来自 κ，跟算法没有任何关系。

这就是那条公式的实际用法：算法交出了它能交的一切，剩下的账要记在问题头上。知道这一点之后，你就不会再去「优化」一个已经后向稳定的求解器了——那是白费力气。

那什么样的算法是后向不稳定的

回到开头那个泰勒级数。e⁻⁵⁰ 这个问题一点也不病态：f(x) = eˣ 的条件数是 |x|，在 x=−50 处是 50。输入错一点，输出最多错 50 倍——非常温和。

可是算出来是这样：

x	泰勒直接展开	改写成 1/eˣ	相对误差
1	0.36787944117144245	0.36787944117144233	3.0 × 10⁻¹⁶
5	0.006737946999084638	0.006737946999085467	1.2 × 10⁻¹³
10	4.539992967040021e-5	4.539992976248485e-5	2.0 × 10⁻⁹
20	6.147561828914626e-9	2.061153622438558e-9	1.98（错 198%）
30	6.1030424788918156e-6	9.357622968840174e-14	6.5 × 10⁷
50	2041.8329628976246	1.9287498479639178e-22	1.1 × 10²⁵

e⁻⁵⁰ 算出来是 2041.83。一个必定落在 0 和 1 之间的量，算出了两千多。

原因是第 3 章那件事的极端版本。这个级数是交错的（正负交替），而中间项大得吓人：

x = −50 时，最大的那一项是 50⁵⁰/50! ≈ 2.9 × 10²⁰

于是要做的事情是：
  用一堆 10²⁰ 量级的数，正负相加，
  最后得到 10⁻²² 的结果。

中间那些 10²⁰ 的数，每一个自带 10⁴ 量级的绝对误差（ulp(1e20) ≈ 16384）。
它们加起来的时候，正确的部分互相抵消了，
⟨误差的部分没有抵消⟩，于是剩下的全是误差。

要的答案比噪声小 42 个数量级。它没有任何机会。

用后向误差的语言说：你算出来的 2041.83，不是任何一个「邻近的 x」的 eˣ 值。不存在一个接近 −50 的数 x̃ 使得 e^x̃ = 2041.83（那需要 x̃ ≈ 7.62）。这个算法的后向误差是灾难级的，所以它不是后向稳定的——尽管它每一步都只抹了半个 ulp。

✎ 术语正名

稳定（stable）这个词在不同领域里意思完全不同，容易串台：

数值稳定（这一章）：算法不放大误差。反义词是不稳定。
控制系统稳定（《过冲》那本书）：系统受扰动后回到平衡点，不发散。
排序稳定：相等元素保持原有相对顺序。
API 稳定：接口不变。

还有一对更细的区分，看论文时会碰到：后向稳定（backward stable，后向误差 ~ eps）比数值稳定（numerically stable，前向误差不比 κ·eps 差太多）要求更强。前者蕴含后者。这本书只在后向稳定的意义上用「稳定」。

⌨ 自己跑一遍

这一章最值得亲手跑的是那个泰勒级数——因为「每一步都没错，结果错了 25 个数量级」这件事必须自己看一遍才信：

import math
def taylor_exp(x, n=400):
    term, s = 1.0, 1.0
    biggest = 0.0
    for k in range(1, n):
        term *= x / k
        biggest = max(biggest, abs(term))
        s += term
        if abs(term) < 1e-20 * abs(s): break
    return s, k, biggest

for x in (-1, -5, -10, -20, -30, -50):
    s, k, big = taylor_exp(float(x))
    print(f'x={x:>4}  {k:>3} 项  泰勒={s:>24.17g}  真值={math.exp(x):>24.17g}'
          f'  最大项={big:.3e}')

会打出：

x=  -1   23 项  泰勒=     0.36787944117144245  真值=     0.36787944117144233  最大项=1.000e+00
x= -10   65 项  泰勒=  4.5399929670400212e-05  真值=  4.5399929762484854e-05  最大项=2.756e+03
x= -30  127 项  泰勒=  6.1030424788918156e-06  真值=  9.3576229688401749e-14  最大项=7.716e+11
x= -50  169 项  泰勒=      2041.8329628976246  真值=  1.9287498479639178e-22  最大项=2.926e+20

最后一列是关键：中间项最大到 2.9×10²⁰，而答案是 1.9×10⁻²²。这两个数中间隔了 42 个数量级——那就是被吞掉的信息量。

python3 taylor.py

在线：python.org/shell。JS 版逐位相同。

▸ 在现实里

数值库的文档就是用这套语言写的。LAPACK 每个例程的文档里都有一段误差界，写的都是后向误差（形如「计算出的解是 (A+E)x̂ = b 的精确解，其中 ‖E‖ ≤ c(n)·eps·‖A‖」）。看懂这一句，你就知道这个库承诺了什么、没承诺什么。它从不承诺答案准，它承诺自己老实。

机器学习里的同构说法。「我们的模型在训练集上完美拟合」——这是前向误差为零。而真正该问的是后向问题：这个模型是「哪一份数据」的正确答案？如果换掉一个样本模型就大变，那这个学习问题的条件数很大，训练算法再好也救不了。这正是「模型稳定性」「影响函数」这一整支研究在做的事。

工程判断上的用法。下次看到「结果不对」，先分两步问：① 把我的结果代回原式，残差大不大？（后向误差）② 这个问题本身敏感吗？（条件数）。残差小 + 问题敏感 ⇒ 你的代码没问题，是问题太难，去改问法；残差大 ⇒ 是代码的锅，去查代码。这两步能省掉大量瞎找。

✗ 这个直觉是错的

「好算法就是给出准确答案的算法。答案不准，算法就有问题。」

在病态问题上，准确答案是不存在的——上一章讲过，输入里就没有那么多信息。一个坚持要给出「准确答案」的算法，只会是在自欺欺人。

正确的标准是后向稳定：我给出的答案，是一个和你的输入几乎无法区分的问题的精确答案。希尔伯特那张表里，n=14 时答案比真值大五倍多，但 LAPACK 完全没有失职——它的后向误差是 2.25×10⁻¹⁷，比 eps 还小。

反过来也成立，而且更常被忽略：答案看起来对，不等于算法稳定。泰勒级数在 x=−1 处给出 16 位全对的答案，看起来完美；它在 x=−50 处给出 2041.83。一个算法在温和输入上的表现，说明不了它在极端输入上的表现。这就是为什么要看误差界，而不是看几个测试用例。

◇ 对账

正确答案是 D：泰勒级数算出 2041.8329628976246，而真值是 1.9287498479639178 × 10⁻²²。相对误差 1.1 × 10²⁵。

A 「级数收敛得快，每项误差不超过半个 ulp」——两句都对，结论错，而且错在一个很典型的地方：「每一步的相对误差都小」推不出「结果的相对误差小」。因为「每一步」的相对误差是相对那一步的中间结果（最大到 2.9×10²⁰），而最终答案是 1.9×10⁻²²。分母换了 42 个数量级。这是第 3 章那个错误推理的极端版。 B 「对前四五位」——这个答案背后的想法是「误差会累积，但也就累积一点」。在同号相加时这个想法基本成立（第 4 章）。这里的级数是交错的，正确的部分互相抵消了，而误差没有——抵消掉的越多，剩下的噪声占比越大。交错级数是相消的重灾区。 C 「量级还在，但一位不对」——这在 x=−30 附近是对的（算出 6.1×10⁻⁶，真值 9.4×10⁻¹⁴，还剩个「很小的数」的样子）。到 x=−50 就连量级都保不住了，因为剩下的纯噪声大小由中间项决定（10²⁰ × eps ≈ 10⁴ 量级的绝对误差），和真答案完全脱钩。C 的直觉是对的，只是低估了退化的速度。

◈ 换一种写法

把一个后向不稳定的算法改成稳定的，通常只需要绕开相消：

taylor_exp(-50)——交错级数，中间项 2.9×10²⁰，答案 1.9×10⁻²²。 1.0 / taylor_exp(50)——同一个级数，x 取正号，所有项同号，一次相消都没有。答案 16 位全对。

一个负号的事。这也正是所有标准库里 exp 的实现思路：先把指数规约到一个小区间（x = k·ln2 + r，|r| ≤ ln2/2），在小区间上用多项式逼近，最后乘上 2ᵏ——永远不在大数上做交错求和。

这一栏的通则，第 3 章已经说过一次，这里再说一遍因为它太重要：让式子里不要出现「大数抵消成小数」的步骤。做不到就换表示，还做不到就承认自己算不了，别硬算。

这一章的一句话

好算法的承诺不是「答案准」，而是「我的答案精确地属于一个和你的输入几乎一样的问题」；把这份承诺乘上问题的条件数，就得到你实际能拿到的准确度上界。

下一章把这条公式推到一个让人不舒服的结论上。你会看到一个方程组：解出来的答案百分之五十四是错的，而把它代回原方程算残差，残差是 4.44 × 10⁻¹⁶——小到不能再小。更狠的是，这一章会把精度提到 60 位十进制再算一遍，答案照样错。残差小，从来不是答案对的证据。