12.3 条件数与数值稳定性：让误差不被算法放大

有限输入也可能让模型损失变成 NaN：直接计算 softmax 时，指数先上溢为无穷大，再形成 inf/inf。相近数相减造成的有效位丢失则是另一类问题，往往先得到不准确的有限结果。应沿计算路径检查，不能只看纸面公式。

数值分析要分清两件事：问题本身是否敏感，以及算法是否额外放大了敏感性。

本课目标

• 区分前向误差、后向误差与问题条件；
• 识别灾难性消减、上溢、下溢和累计误差；
• 实现稳定的 softmax、log-sum-exp 与求和；
• 为数值代码设计边界测试和误差检查。

1. 条件描述问题，稳定性描述算法

若输入 $x$ 有小扰动 $\mathrm{\Delta }x$，输出 $f(x)$ 的变化可能很小，也可能很大。条件数概括问题对输入扰动的敏感程度。

对可逆矩阵的线性系统 $A\mathbf{x}=\mathbf{b}$，相应诱导矩阵范数下的条件数是：

$$\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert .$$

条件数很大时，输入或舍入的微小相对误差可能在解中被显著放大。条件数为无穷对应奇异矩阵。条件数依赖所选范数，也不单独给出每次计算的实际误差。

这里用后向稳定性来衡量算法：它问：计算结果是否接近某个轻微扰动输入的精确答案。一个病态问题无法靠稳定算法变得不敏感，但不稳定算法会让原本条件良好的问题也产生糟糕结果。

2. 前向误差与后向误差

设精确答案为 $y=f(x)$，计算答案为 $\hat{y}$：

• 前向误差比较 $\hat{y}$ 与 $y$；
• 后向误差寻找一个小扰动 $\mathrm{\Delta }x$，使 $\hat{y}=f(x+\mathrm{\Delta }x)$。

后向稳定意味着算法算出的结果，可以解释为“对几乎相同的输入做了精确计算”。若问题条件良好，小后向误差通常导向小前向误差；若问题病态，即使后向稳定，前向误差也可能很大。

这解释了为什么只打印最终小数位不够。还要检查残差、输入尺度和条件数。例如线性系统可检查：

$$\mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}.$$

小残差是重要信号，但病态系统中“小残差”和“解接近真实参数”仍不是同一件事。

下面用 NumPy 2.x 构造对角系统。右端向量相对扰动只有约 1e-8，第二个解分量却从 1 变成 2，相对于原方程的残差也只有 1e-8：

import numpy as np

A = np.diag([1.0, 1e-8])
b = np.array([1.0, 1e-8])
perturbed_b = b + np.array([0.0, 1e-8])
x = np.linalg.solve(A, b)
changed_x = np.linalg.solve(A, perturbed_b)
relative_input = np.linalg.norm(perturbed_b - b) / np.linalg.norm(b)
residual = np.linalg.norm(b - A @ changed_x)
print(f"{np.linalg.cond(A):.0e} {relative_input:.0e} {residual:.0e}")
print(x, changed_x)
# 1e+08 1e-08 1e-08
# [1. 1.] [1. 2.]

这里第二个方程自身的相对扰动很大；按整个向量的范数衡量才显得很小。范数式与逐分量误差模型回答不同问题，单位缩放也会影响这种判断。不能只看一个裸残差就宣布参数准确。

3. 灾难性消减

两个非常接近的浮点数相减时，前面的有效位被抵消，输入中已有的舍入误差相对于结果会被放大。

典型例子是小 $x$ 时计算：

$$\sqrt{1+x}-1.$$

把分子有理化：

$$\sqrt{1+x}-1=\frac{x}{\sqrt{1+x}+1}.$$

数学上等价，数值上后者避免了两个接近 1 的数相减。

import math

def unstable(x: float) -> float:
    return math.sqrt(1.0 + x) - 1.0

def stable(x: float) -> float:
    return x / (math.sqrt(1.0 + x) + 1.0)

for x in [1e-8, 1e-12, 1e-16]:
    print(x, unstable(x), stable(x))

“避免大数减小数”不是准确规则；真正危险的是两个接近的数相减，尤其结果远小于操作数时。

两种实数公式都要求 $x\ge -1$。当 x=1e-16 时，binary64 的 1+x 已舍入为 1，直接公式返回零；改写仍保留了分子中的 x，给出约 5e-17。相减本身有时可以精确，灾难性消减放大的是操作数此前已有的误差。

4. 稳定的 softmax 与 log-sum-exp

直接计算

$$\mathrm{softmax}(x_i)=\frac{e^{x_i}}{\sum _j{e}^{x_j}}$$

可能因大正数上溢。利用平移不变性，令 $m=\underset{i}{max}{x}_i$：

$$\mathrm{softmax}(x_i)=\frac{e^{x_i-m}}{\sum _j{e}^{x_j-m}}.$$

import math

def softmax(values: list[float]) -> list[float]:
    if not values:
        raise ValueError("values must not be empty")
    if not all(math.isfinite(value) for value in values):
        raise ValueError("only finite inputs are accepted")

    maximum = max(values)
    shifted = [math.exp(value - maximum) for value in values]
    denominator = math.fsum(shifted)
    return [value / denominator for value in shifted]

print([round(p, 6) for p in softmax([1000.0, 1001.0, 1002.0])])
# [0.090031, 0.244728, 0.665241]

对应的 log-sum-exp：

$$\log \sum _i{e}^{x_i}=m+\log \sum _i{e}^{x_i-m}.$$

减去最大值让最大指数参数为 0、对应指数函数值为 1，避免正向上溢。非常小的项仍可能下溢为零，但它们对总和的相对贡献本来就可能低于当前精度；是否可接受要由误差要求决定。

继续使用上面的 Python 3.9+ 环境，可把 log-sum-exp 也写成完整函数。这里统一拒绝空输入与非有限值，不替其他库定义特殊值语义：

def logsumexp(values: list[float]) -> float:
    if not values or not all(math.isfinite(x) for x in values):
        raise ValueError("require a nonempty list of finite values")
    maximum = max(values)
    return maximum + math.log(math.fsum(math.exp(x - maximum) for x in values))

print(f"{logsumexp([1000.0, 1001.0, 1002.0]):.6f}")  # 1002.407606
try:
    math.log(math.fsum(math.exp(x) for x in [1000.0, 1001.0, 1002.0]))
except OverflowError:
    print("direct exponentiation overflowed")

对非空有限输入，至少一项为 1，所以 softmax 分母不会因所有项下溢而变成零。极小概率仍可能变零；如果随后需要对数概率，应保留平移后的计算 (x_i-m) - log(sum_j exp(x_j-m))，避免 log(0)，也避免先把小的对数修正加到巨大偏移量上，再相减时丢失。

5. 累计求和也需要算法

大量不同数量级的数相加时，小量可能不断被舍掉。改进方法包括：

• 从绝对值较小的数开始求和；
• 成对求和，降低误差随项数增长的速度；
• Kahan 等补偿求和，保留部分丢失低位；
• 使用语言或库提供的高质量实现，如 Python math.fsum。

这些方法提高精度，不保证任意输入上的精确结果。并行系统还需在速度、确定性和可复现性之间权衡。

Python 3.12 已改进内置浮点 sum，因此旧版本的 sum(values)==0 不能作为跨版本预期。下面用显式逐项加法展示丢失，另打印当前版本的内置结果。math.fsum 使用多个部分和改善精度，仍不是任意精度算术；见 Python 3.12 变更 和 fsum 文档。

import math
import sys

values = [1e16, 1.0, -1e16]
sequential = 0.0
for value in values:
    sequential += value
print(sequential)         # 0.0
print(math.fsum(values))  # 1.0
print(sys.version_info[:2], sum(values))  # Observed with CPython 3.12: (3, 12) 1.0

6. 选择公式时先看边界

数值稳定的实现经常使用专门函数：

• 小 $x$ 时用 log1p(x) 计算 $\log (1+x)$；
• 小 $x$ 时用 expm1(x) 计算 $e^x-1$；
• 用 hypot(x, y) 计算 $\sqrt{x^2+y^2}$，降低中间上溢或下溢风险；
• 解线性方程用分解和 solve，不显式求逆；
• 概率损失尽量从 logits 直接计算，不先舍入成概率。

“数学等价”不保证“浮点执行等价”。算法设计要把中间值范围、消减、分支连续性和舍入次数纳入考虑。

7. 数值代码的验证清单

1. 使用解析答案或高精度结果建立小规模基准；
2. 测试零、极小量、极大量、相近值和符号变化；
3. 明确是否接受 NaN、无穷大和次正规数；
4. 检查绝对误差、相对误差、残差或守恒量；
5. 改变运算顺序、批大小和并行度，观察敏感性；
6. 对迭代算法记录收敛历史，不只看最后一次输出；
7. 记录数据类型、舍入规则、库版本和硬件假设。

误差阈值应与问题尺度相连。只断言“小数点后六位相等”，可能对零附近太宽松，对大数又太严格。

常见误区

• 公式正确，实现就正确：中间运算可能溢出、下溢或消减。
• 高精度能修好任何问题：它推迟误差暴露，但不修复病态模型或错误算法。
• 残差小就说明参数准确：病态问题中二者可能相差很远。
• 稳定算法会给出精确答案：稳定性控制误差放大，不消除输入和舍入误差。

练习

1. 比较 sqrt(1+x)-1 两种公式在不同 $x$ 下的相对误差。
2. 实现稳定的 logsumexp，并与直接公式比较大输入结果。
3. 构造一个条件数较大的 $2\times 2$ 线性系统，轻微扰动 $b$ 后比较解。
4. 为批量聚合服务设计“快速非确定求和”和“可复现求和”两种模式。

小结

条件数描述问题对扰动的先天敏感性，稳定性描述算法是否额外放大误差。识别消减、控制中间值范围、选择稳定函数并验证极端输入，是把数学公式变成可靠软件的必要步骤。

• LAPACK: Error Bounds for Linear Equation Solving