引言
计算机中的浮点算术不是实数算术：实数集是连续完备的，而浮点数是有限的离散集合，每次运算都可能引入舍入。数值分析的根本问题因此是双重的——除了数学模型离散化带来的截断误差，还要面对有限精度运算带来的舍入误差。一个数学上等价的公式，在浮点算术中可能一个稳定、一个完全失效。本文系统梳理误差的来源与度量：浮点表示、前向/后向误差、条件数与数值稳定性的关系，并用三个经典数值实验展示理论如何指导安全算法的设计。
理论基础
浮点表示与舍入模型
IEEE 754 双精度数由 1 位符号、11 位指数与 52 位尾数构成，机器精度（单位舍入）为 $u = 2^{-53} \approx 1.11\times10^{-16}$（相对舍入误差上界；常说的 $\varepsilon_{\text{mach}} = 2u$）。标准舍入模型断言：对算术运算 $\mathrm{op} \in {+, -, \times, /}$，
$$\mathrm{fl}(x \mathrm{op} y) = (x \mathrm{op} y)(1 + \delta), \qquad |\delta| \le u$$
一切舍入误差分析都建立在这个模型上。需警惕的是：该模型假设 $x, y$ 本身已是浮点数，且下溢/溢出区域不成立。
截断误差与舍入误差的竞争
截断误差来自用有限步逼近无限过程（差商代替导数、级数截断、迭代中止），随步长 $h$ 按幂律下降；舍入误差随 $h$ 减小反而累积。用中心差商 $f'(x) \approx [f(x+h) - f(x-h)]/(2h)$ 求导时，截断误差 $\mathcal{O}(h^2)$ 而舍入误差 $\mathcal{O}(u/h)$，总误差在 $h \approx u^{1/3} \approx 10^{-5}$ 处取极小——「步长越小越精确」在浮点世界中并不成立。
条件数：问题的固有属性
条件数度量问题本身对数据扰动的敏感程度。对标量可微函数，相对条件数为
$$\kappa_f(x) = \left|\frac{x f'(x)}{f(x)}\right|$$
对线性方程组 $Ax = b$，扰动理论给出
$$\frac{|\delta x|}{|x|} \le \kappa(A), \frac{|\delta b|}{|b|}, \qquad \kappa(A) = |A|,|A^{-1}|$$
$\kappa(A)$ 大意味着右端的微小相对扰动可被放大 $\kappa(A)$ 倍——这是问题的内在性质，任何算法都无法规避。大致的实用判读是：$\kappa(A) \approx 10^k$ 时，双精度解约有 $16 – k$ 位有效数字。
数值稳定性：算法的属性
算法 $f$ 的近似实现 $\tilde{f}$ 称为后向稳定的，若计算解恰为邻近输入的精确解：$\tilde{f}(x) = f(\tilde{x})$ 且 $|\tilde{x} - x| = \mathcal{O}(u)|x|$。前向误差、后向误差与条件数由基本不等式联系：
$$\text{前向误差} less\sim \text{条件数} \times \text{后向误差}$$
于是「好问题 + 稳定算法 = 可靠结果」：列主元 LU、Householder QR 都是后向稳定算法，而代数余子式求逆、Gram-Schmidt 经典形式则不稳定。灾难性抵消（两个相近数相减损失有效位）是最常见的稳定性陷阱：若 $x \approx y$，则 $x – y$ 的相对误差被放大 $|x|/|x-y|$ 倍。
数值方法与实现
下面三个实验分别展示：病态矩阵（Hilbert 矩阵 $H_{ij} = 1/(i+j-1)$）条件数的爆炸、二次方程求根公式的灾难性抵消及其修复、Kahan 补偿求和对累积舍入的抑制。
import numpy as np
# 实验1: Hilbert 矩阵条件数随阶数指数增长
for n in (5, 8, 11, 14):
i = np.arange(1, n + 1)
H = 1.0 / (i[:, None] + i[None, :] - 1.0)
print(f"n={n:2d} cond2(H)={np.linalg.cond(H):.3e}")
# 实验2: x^2 - 2b x + 1 = 0 的小根, b >> 1 时直接公式灾难性抵消
b = 1e8
x_naive = b - np.sqrt(b**2 - 1.0) # 两个相近大数相减
x_fixed = 1.0 / (b + np.sqrt(b**2 - 1.0)) # 分子有理化改写
print(f"小根 直接公式: {x_naive:.12e} 改进公式: {x_fixed:.12e}")
# 实验3: float32 顺序求和 vs Kahan 补偿求和
def kahan32(x):
s = np.float32(0.0); c = np.float32(0.0)
for v in x:
y = v - c
t = s + y
c = (t - s) - y
s = t
return float(s)
xs = np.full(1_000_000, np.float32(0.1), dtype=np.float32)
s_naive = np.float32(0.0)
for v in xs:
s_naive += v
print(f"真值 100000.0 朴素: {float(s_naive):.4f} Kahan: {kahan32(xs):.4f}")实验 2 中 $b = 10^8$ 时 $\sqrt{b^2 - 1}$ 与 $b$ 前 16 位几乎相同，直接公式只剩约 8 位有效数字，有理化改写 $\frac{1}{b + \sqrt{b^2-1}}$ 则完全避开抵消。实验 3 中朴素 float32 累加在 $N$ 较大时相对误差可达数个百分点，而 Kahan 算法以每步存一份「补偿项」的代价把误差控制在 $\mathcal{O}(u)$ 级。
工程案例分析
舍入与溢出的后果远不止论文里的数值实验。1991 年海湾战争中，Patriot 导弹防御系统用 24 位定点数存储 0.1 秒时钟刻度，二进制截断每拍引入约 $9.5\times10^{-8}$ 秒误差，连续运行约 100 小时后时钟偏差累积至约 0.34 秒——足以让跟踪雷达错过以高速飞来的目标，造成 Dhahran 军营被 Scud 导弹击中的事故。1996 年 Ariane 5 首飞爆炸，直接原因是惯性制导系统把 64 位浮点水平速度量转换为 16 位有符号整数时溢出，而该段代码沿用了 Ariane 4 的参数范围、未对新火箭的更大速度做保护。两个案例的共同教训是：数值问题必须在系统级设计中被显式建模与验证，而不能依赖「过去一直没问题」。
| 运算 | 相对误差传播 | 风险点 | ||||||
|---|---|---|---|---|---|---|---|---|
| $x \pm y$ | $sim frac{ | x | + | y | }{ | x pm y | }(varepsilon_x + varepsilon_y)$ | 相近数相减（灾难性抵消） |
| $x \cdot y$、$x / y$ | $\sim \varepsilon_x + \varepsilon_y$ | 溢出/下溢 | ||||||
| $\sqrt{x}$ | $\sim \frac{1}{2}\varepsilon_x$（良态） | 负数输入 | ||||||
| 长序列求和 | $\sim N u$（朴素顺序） | 应改用 Kahan 或成对求和 | ||||||
| $Ax = b$ 求解 | 放大 $\kappa(A)$ 倍 | 病态矩阵需正则化 |
结论与展望
误差理论把「结果可信吗」这一模糊问题分解为可计算的量：问题条件数给出精度上限，算法后向稳定性给出能否逼近该上限的保证，二者相乘即前向误差的估计。值得进一步关注的方向包括：区间算术与仿射算术（自动给出严格误差界）、随机舍入与混合精度算法（在低精度硬件上恢复统计意义下的精度），以及自动微分中的数值稳定性分析——它已成为深度学习系统可信性的基础议题。
发表回复