消去误差与稳定重写

形式陈述 ​

设精确量为 $x,y$，而参与计算的数已经带有扰动

$$\hat{x}=x(1+{\delta }_x),\hat{y}=y(1+{\delta }_y),|{\delta }_x|,|{\delta }_y|\lesssim u.$$

当 $x\approx y$ 时，真差 $d=x-y$ 远小于两个操作数。即使最后一次减法符合浮点算术标准误差模型公理库浮点算术标准误差模型Standard floating-point arithmetic model以每次基本运算的小相对扰动和 gamma 记号组织多步浮点误差分析。，输入扰动在差中的相对影响仍约为

$$\frac{|x{\delta }_x-y{\delta }_y|}{|x-y|}\le u\frac{|x|+|y|}{|x-y|}.$$

分母越小，先前舍入、测量或函数求值误差越容易吞没结果的有效数字。所谓灾难性消去，指的正是这种相对于小结果的误差放大；减法指令本身可能仍然正确舍入。

相减并非一概危险。Sterbenz 引理说明，在通常采用渐进下溢的浮点系统中，若两个同号浮点数满足 $y/2\le x\le 2y$，且结果可表示，则 $x-y$ 可以精确计算。若两个操作数本来就是准确数据，这类消去甚至可能是良性的；危险来自它把操作数中已经存在的误差暴露在更小的结果尺度上。

稳定重写不改变实数中的数学值，而是改变有限精度的中间量。常用做法包括提取公共因子、用共轭有理化避免相近根式直接相减，以及先缩放再计算，防止中间溢出或下溢。接近特殊点时，应优先使用为小增量设计的函数，例如以 log1p(x) 计算 $\log (1+x)$、以 expm1(x) 计算 $e^x-1$；它们避免先把 $1+x$ 或 $e^x$ 舍入到 $1$。

有些公式必须按数据分支。分支的目的不是猜测答案，而是选择不会让两个大而接近的量相减的代数形式。判断一项重写是否更稳定，需要结合输入范围、条件性和允许的异常值，而不能只看表达式更长或使用了更多位数。

直觉 ​

消去像用两次大额读数求一笔很小的差额。若两次读数都只精确到个位，拿“100000001 减 100000000”得到的个位差并没有继承七位可靠数字；前面的共同高位消失后，原先藏在末位的误差成了结果的主要部分。

稳定公式的工作方式是绕开这次信息对撞。它利用根与系数的关系、因式分解或局部专用函数，让小量直接由同尺度的量产生。这样做不能改善问题本身的条件性，却可以避免算法额外浪费输入仍然携带的精度。

例子与边界 ​

考虑方程

$$x^2-{10}^{8}x+1=0.$$

直接用二次公式计算小根，binary64 运算会先形成

$$x_{-}=\frac{{10}^8-\sqrt{{10}^{16}-4}}{2},$$

其中两个约为 ${10}^8$ 的数几乎相等。一次 binary64 实验得到

$${\hat{x}}_{-}^{\mathrm{direct}}=7.450580596923828\times {10}^{-9},$$

而高精度参考约为

$$x_{-}=1.0000000000000001\times {10}^{-8}.$$

直接公式的相对误差约为 $0.25494$。这不是因为最终除以 $2$ 困难，而是相减前两个操作数已各自舍入，差值只剩很少可靠位。

两根乘积为 $1$。先用加号稳定地计算大根

$$x_{+}=\frac{{10}^8+\sqrt{{10}^{16}-4}}{2},$$

再令 $x_{-}=1/x_{+}$，实验得到 ${10}^{-8}$，相对于高精度参考的误差约为 ${10}^{-16}$。对一般方程 $a{x}^2+bx+c=0$，常令

$$q=-\frac{1}{2}(b+\mathrm{sign}(b)\sqrt{b^2-4ac}),$$

再计算 $q/a$ 与 $c/q$；判别式、$q=0$ 和系数缩放仍须分别处理。

根式差也有同样结构：

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

当 $x$ 很大时，左式可能先把两个根式舍入成同一个数，右式却直接形成正确尺度的小量。不过，有理化并非无条件更好；若分母可能溢出、输入含无穷或表达式处在其他区域，还要配合缩放和分支。

提高工作精度能推迟精度耗尽，却不能代替结构重写。只要两个操作数的相对误差最终仍被 $|x-y|$ 的小尺度放大，同一问题会在更苛刻的数据上再次出现。

推论与应用 ​

导数公理库导数Derivative函数增量比在步长趋零时的极限。的有限差商会相减相近函数值，几何中的方向判定公理库方向判定Orientation test用二维或高维行列式符号判断点组转向或仿射定向的基本谓词。会从若干大乘积求一个很小的行列式符号，二者都需要把“公式在实数中正确”与“计算路径保留了多少信息”分开。前者常通过步长权衡或更换求导方式处理，后者在接近退化时可能需要自适应精度或精确谓词。

求根、特殊函数、方差计算和几何距离也经常存在多种代数等价形式。选择公式时，应先定位可能发生的消去、溢出与下溢，再用数值稳定性公理库数值稳定性Numerical stability · Backward stability以允许的小输入扰动刻画算法的有限精度行为，并与问题条件性及其他稳定性概念分开。所要求的扰动模型检验重写，而不是把“少用减法”当成机械规则。

参考资料

• David Goldberg, “What Every Computer Scientist Should Know About Floating-Point Arithmetic,” ACM Computing Surveys 23(1), 1991.
• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Chs. 1–2.
• IEEE, IEEE Standard for Floating-Point Arithmetic, IEEE Std 754-2019.