本页目录
数值 I · 误差与数值稳定性
数学恒等式给出同一个实数，机器的运算路径却可能给出不同近似。本页学习把输入误差、截断误差、舍入误差分开，再判断：应当改写算法、提高精度，还是重新审视问题本身？
学习层：同一个公式，哪一步把信息丢了？
1. 从一个能亲手检查的矛盾开始
在 binary64 的常用“舍入到最近、正中取偶”模式下，\(10^{16}\) 附近相邻数差 2。因此先算 \(10^{16}+1\)，会舍入回 \(10^{16}\)；再减去 \(10^{16}\)，得到 0。实数的精确和却是 1。
这不是减法那一步必然算错：小项在前面的加法中已经丢了。把顺序改为 \((10^{16}-10^{16})+1\)，机器就能得到 1。说明“问题的答案”和“求答案的路径”需要分别检查。
2. 先预测，再看三个不同的账
选择根式、余弦、差分或求和。先回答当前结果是否为 0（求和时是否精确），再判断算法是否具有普遍保证；揭示后才能看结果和误差图。换参数或改答案都会重新隐藏结果。
- 输入账：本实验的 \(x,h\) 取固定 binary64 网格；参考也针对这个二进制数，而不是尚未舍入的十进制数。
- 算法账：实际双精度输出，与独立的 80 位参考比较。稳定改写仍可能有非零误差。
- 截断账：差分公式本身也只是导数的近似，即使每步实数运算完全准确，也不等于导数。
3. 两个恒等改写
对本实验 \(x>0\)，
根式左边可能先把 \(1+x\) 或它的平方根舍入成 1，接着相减只剩 0；右边直接算一个接近 \(1/2\) 的数。余弦左边也可能把 \(\cos x\) 舍入成 1；右边保留了小角度的信息。
这里的“更稳定”指常见小量情境下误差更受控，不是保证每个输入都比直接式精确，更不是保证零误差。当 \(x\ne0\) 时根式两式相等；右边还能连续延拓到 \(x=0\)，左边原式在那里没有定义。
4. 差分：把步长缩小也有代价
用
近似 \(f'(1)=e\)。第二式是这个特殊函数的差分改写，不是把参考值 \(e\) 直接放进“更好算法”一列。两种计算都与同一个导数参考比较。
实数 Taylor 展开给出
直接式通常还受约 \(O(u/h)\) 的舍入放大影响，粗略包络为 \(Ah^2+Bu/h\)。它解释了“太大、太小都可能不好”，但实际采样有锯齿和偶然抵消，不应画成保证单调的光滑 U 形。特殊的 sinh 改写避开了这条相消路径；不能据此假设任意黑箱函数都有这种改写。
5. 求和：一个成功例子和一个失败例子
比较朴素、Kahan、Neumaier 三种求和。对 \(10^{16}\) 后跟 20000 个 1 再抵消的大数，Kahan 能保存朴素方法丢掉的小项。但是：
精确和为 4。Neumaier 在这个例子成功，并不保证对任何有限序列都精确。也试试奇数个 1：先聚合小项再加大数，仍可能在最后一次舍入时丢失或多出 1。
无 JavaScript 时的静态读法：
| 情境 | 直接路径 | 改写或补偿 | 检查重点 |
|---|---|---|---|
| 根式，\(x\approx10^{-16}\) | 先形成 \(1+x\)，可能变成 1 | 有理化结果接近 \(1/2\) | 参考值不能由被评估算法自己定义 |
| 余弦，\(x\approx10^{-8}\) | \(1-\cos x\) 可能为 0 | \(2\sin^2(x/2)\approx5\times10^{-17}\) | 相消发生前的求值舍入 |
| 导数，\(h\approx0.1\) | 中央差分有截断误差 | sinh 改写仍有同一截断误差 | 改算法不能消除公式近似 |
| \(10^{16},1,-10^{16},3\) | 结果 3 | Kahan 3，Neumaier 4 | 一个成功样本不是普遍定理 |
图中“精确 0”单独成行，不把 0 假装成 \(10^{-18}\)。曲线只连接所列采样；同坐标的点可能重合。参考以高、低两个浮点数保存来估算残差；这能看见低于一个 binary64 间距的误差，但不是严格区间证明。浏览器的超越函数实现允许末位差异。
6. 用结论时，先带上边界
条件数依赖你如何描述输入和度量扰动；后向稳定需要一个明确的“邻近输入”定义。比较近似数常用绝对与相对容差组合，但计数和经过证明的精确整数可以直接比较。提高输出精度不能恢复输入阶段已经抹掉的信息。
7. 迁移题
题 1：小心“减法丢位”。设输入是两个已知的精确 binary64 数 \(x=1+2^{-52}\)、\(y=1\)。求机器差值；再说明为什么这不否定“相近的测量值相减可能很不可靠”。
展开推导
差为可精确表示的 \(2^{-52}\)，此减法精确。若原始测量值却各带绝对误差至多 \(\eta\)，差的输入误差可达 \(2\eta\)，相对差值便可能很大。相近浮点数在 Sterbenz 条件下能精确相减，仍无法修复此前的输入或函数求值误差。
题 2：步长怎么选？对正数 \(A,B,u\)，最小化误差包络 \(E(h)=Ah^2+Bu/h\)，\(h>0\)。这能否精确预测每个浏览器里最好的采样步长？
展开推导
\(E'(h)=2Ah-Bu/h^2\)，故 \(h_*=(Bu/(2A))^{1/3}\)。且 \(E''(h)=2A+2Bu/h^3>0\)，这是唯一全局最小点。它预测量级，不能给出实际离散误差的精确极小位置：\(A,B\) 与函数和实现有关，真实误差不是两个正项的精确和，还会有偶然抵消。
题 3：稳定 softmax。输入为有限实数 \(x_1,\ldots,x_n\)，证明同减 \(m=\max x_i\) 不改变 softmax；说明它保护哪一步，以及已经量化的输入能否恢复。
展开推导
实数中分子分母同乘 \(e^{-m}\)，故 \(e^{x_i-m}/\sum_j e^{x_j-m}=e^{x_i}/\sum_j e^{x_j}\)。平移后最大指数为 1，其余在 \((0,1]\)，避免大正指数溢出；机器中很小的项仍可能下溢为 0，因此不能声称每个概率都严格正。已量化成相同值的两个输入不会因平移恢复差别。这个推导假定输入有限；含无穷或 NaN 的程序接口需要单独规定行为。
1. 误差与计算问题
四类来源是模型误差、观测误差、截断误差、舍入误差。给定精确目标 \(y\) 与近似 \(\widehat y\)，绝对误差为 \(|\widehat y-y|\)；当 \(y\ne0\)，相对误差为 \(|\widehat y-y|/|y|\)。目标为 0 时应使用绝对误差或明确的外部尺度，不能偷偷用“最小正数”替代分母。
有效数字必须结合首位和采用的定义判断；“相对误差约为 \(10^{-p}\)”可以描述精度量级，不自动等于任意定义下恰有 \(p\) 位有效数字。
2. binary64：间距、舍入和范围
正常数写作 \((-1)^s(1.b_1\cdots b_{52})_2\,2^e\)，有 53 位二进制有效精度。在区间 \([2^e,2^{e+1})\)，相邻数间距是 \(2^{e-52}\)。因此 1 的左邻点是 \(1-2^{-53}\)，右邻点是 \(1+2^{-52}\)。
- 机器 epsilon 通常指 \(1\) 到右邻点的间距：\(\varepsilon_{\rm mach}=2^{-52}\)。
- 舍入到最近时常用单位舍入 \(u=2^{-53}\)。在无溢出且精确结果落在正常数范围的标准条件下， [ \operatorname{fl}(a\circ b)=(a\circ b)(1+\delta),\qquad|\delta|\le u. ]
- 非正常数采用固定间距 \(2^{-1074}\)，能逐渐接近 0，却不能维持同一个相对误差界；最小正正常数是 \(2^{-1022}\)。
- 超越函数不应无条件套用“四则运算正确舍入”结论；语言规范和具体数学库需要分别核对。
“正中取偶”的偶，是保留尾数最低位为 0，不是看十进制末位。有限精度还破坏结合律；容差比较本身也通常不具传递性，不能直接作为排序等价关系。
3. 相消：误差可能早于减法发生
设输入近似值为 \(\widehat x=x+\Delta x,\widehat y=y+\Delta y\)。即使减法精确，
分母很小，就会放大先前的误差。不能只数差值剩几个非零数字，就断言减法自身丢掉了多少准确位。
Sterbenz 引理给出一个精确性窗口：同一二进制格式的非负浮点数 \(x,y\)，若 \(y/2\le x\le2y\)，且采用支持渐进下溢的相应运算，则 \(x-y\) 可精确计算。它与灾难性相消不冲突：精确减去两个已经错误的近似值，仍可得到很差的真实差值近似。
稳定改写的三个例子：
| 直接表达式与范围 | 改写 | 仍须检查 |
|---|---|---|
| \(\sqrt{x+1}-\sqrt x\)，\(x\ge0\) | \(1/(\sqrt{x+1}+\sqrt x)\) | 中间结果溢出及输入舍入 |
| \(1-\cos x\)，小角度 | \(2\sin^2(x/2)\) | 下溢和三角函数求值误差 |
| 二次式 \(at^2+bt+c\)，\(a\ne0,\Delta=b^2-4ac\ge0\) | \(q=-\tfrac12(b+\operatorname{sgn}_+(b)\sqrt\Delta)\)，根为 \(q/a,c/q\) | \(q=0\) 特例、判别式和乘积溢出；\(\operatorname{sgn}_+(0)=1\) |
二次式中 \(q=0\) 对应 \(b=c=0\)（实数条件下），两个根均为 0，不能做 \(c/q\)。大尺度系数还需缩放；稳定选符号本身不保证判别式计算安全。小量对数与指数可使用专门的 log1p、expm1 算法，避免先形成接近 1 的量再相消。
4. 条件数来自问题与扰动模型
对可微标量函数 \(f\)，在 \(x\ne0,f(x)\ne0\) 的点，Taylor 展开给出
因此局部相对条件数是 \(\kappa_f(x)=|xf'(x)/f(x)|\)；绝对条件数是 \(|f'(x)|\)。零输入或零输出应改用适合的绝对/混合尺度，而非机械代入。
对 \(g(x,y)=x-y\ne0\)，若 \(|\Delta x|\le\epsilon|x|,|\Delta y|\le\epsilon|y|\)，分量相对条件数为
同样，求和 \(s=\sum x_i\ne0\) 的分量相对条件数为 \(\sum|x_i|/|s|\)。这解释了正负大项抵消的敏感度。即使某个整数样本恰好精确，附近问题仍可能病态。
两个容易混淆的对数问题：
[ f(x)=\log x:\quad\kappa_f=\frac1{|\log x|}\quad(x>0,x\ne1); ] [ g(y)=\log(1+y):\quad \kappa_g=\left|\frac{y}{(1+y)\log(1+y)}\right|\longrightarrow1. ]
接近 1 的 \(x\) 与接近 0 的原始小量 \(y\) 使用了不同的相对扰动尺度。log1p 改善的是给定 \(y\) 的求值路径；若手上只有已经舍入为 1 的 \(x\)，它不能恢复未知的 \(y\)。
5. 前向、后向与稳定性
前向误差直接比较 \(\widehat y\) 与 \(f(x)\)。在指定输入空间和范数下，后向误差寻找最小扰动，使 \(\widehat y=f(x+\Delta x)\)；这种解释有时不存在，不能默认任何输出都对应附近输入。
后向稳定算法使该扰动按问题规模受 \(O(u)\) 控制。在可微、扰动足够小、条件数适用的范围内，前向相对误差一阶上界约为“条件数 × 后向相对误差”。若还存在输入误差，要把输入误差也算进总扰动；这不是无条件等式，也不是对任意大扰动的保证。
一个简单证明模板是将 \(\widehat y=f(x+\Delta x)\) 代入上一节 Taylor 展开：稳定性控制 \(\Delta x\)，条件数控制这段输入变化如何传到输出。病态问题可能让一个后向稳定算法仍产生较大前向误差。
6. 差分的误差估计为何有三分之一次方
若 \(f\) 在邻域内三次连续可微，中央差分的截断误差至多 \(\max|f'''|h^2/6\)。若两端函数求值分别带至多 \(\eta\) 的绝对误差，它们对差分的贡献至多 \(\eta/h\)；实际形成 \(x\pm h\) 的误差和最后除法还需纳入估计。
于是 \(Ah^2+Bu/h\) 是典型尺度模型，最优量级 \(h\sim u^{1/3}\) 来自平衡两项。函数值、坐标尺度、导数量级改变后，常数就变了；不能把某个固定步长当成通用建议。
7. 求和算法到底补偿了什么？
朴素求和每次舍入。标准模型下，长度 \(n\) 的递归求和（无溢出等）满足
证明要点：把每一步写成乘 \(1+\delta_j\)，展开后每个输入带若干舍入因子；用 \(|\prod(1+\delta_j)-1|\le\gamma_k\) 控制累积偏差。除以 \(|s|\) 后，求和条件数自然出现。
Kahan 的每步操作为 \(y=x_i-c,\ t=s+y,\ c=(t-s)-y,\ s=t\)。补偿量 \(c\) 试图保存本次未进入主和的低位。Neumaier 变体根据 \(|s|\) 与 \(|x_i|\) 谁更大，分别积累 \((s-t)+x_i\) 或 \((x_i-t)+s\)，最后返回 \(s+c\)。这些补偿量自身也是浮点数；不能任意重排括号，否则算法意义会改变。
成对求和缩短累积链，补偿求和多保存低位，更多精度增加表示能力；它们解决的路径不同。检查精度时应给出格式、顺序、样本和参考，不能凭“逆序通常更好”宣称固定多出几位。对于整数、精确有理数、高精度参考或区间界，分别说明依据。
8. 从基础走向更深入的数值计算
接下来可把这里的条件数、后向误差和补偿思想带到 LU 选主元、迭代改进、病态最小二乘与微分方程。遇到新算法时，先找明确的输入模型和输出误差目标，再看误差分析或可复现对照；一次曲线漂亮的演示还不构成稳定性证明。
进一步阅读
- Goldberg：What Every Computer Scientist Should Know About Floating-Point Arithmetic：浮点表示、相消与舍入分析的原始论文重印。
- Higham：Accuracy and Stability of Numerical Algorithms：第 2–4 章系统讨论舍入模型和求和。
- ECMAScript 数值与 Math 函数规范：核对浏览器中的 Number 语义与超越函数实现边界。
下一页：线性方程组的数值解——把“残差小”与“解准确”分开。