M10 · 第 1 章 · 第一编 误差与数值线性代数

浮点数、条件数与误差传播

从有限精度舍入模型出发,辨认吸收与消去等误差机制,用条件数区分问题敏感性,用前向和后向误差区分结果偏差与输入扰动,并据此判断算法是否数值稳定。

报告页面错误
预备知识极限与连续性矩阵

本章目标

  1. 用基数、有效位数和舍入方式描述浮点系统,并区分单位舍入误差与相邻数间距。
  2. 使用单次运算舍入模型和累计误差界解释吸收、溢出、下溢与消去。
  3. 计算标量问题的绝对与相对条件数,区分输入敏感性和算法造成的误差。
  4. 分别报告前向误差、后向误差和残差,并用条件数连接这些量。
  5. 通过等价改写、尺度调整和独立复算改善数值稳定性。
页面阅读位置0% · 仅保存在此浏览器
章节未开始
本册完成进度0/6 章 · 0%
本页目录

实数公式进入计算机以后

数学公式通常在实数域中书写,而计算机只能保存有限多个状态。一个十进制数即使输入时看起来很短,转成二进制后也可能是无限展开;一次乘法或开方的精确结果也常常不在可表示集合中。于是计算过程不是“精确实数运算的较快版本”,而是一串运算后紧跟舍入的离散过程。

误差分析要回答四个不同问题:输入数据本身是否精确,问题对输入扰动是否敏感,所用算法是否放大了本可避免的舍入误差,最终数值离精确答案有多远。把这些问题混成一句“有浮点误差”无法指导计算。可靠的报告应指出误差的尺度、来源和传播路径。

浮点表示与 IEEE 风格舍入模型

规格化浮点数

给定基数 β2\beta\ge2、有效位数 pp 和指数范围,非零规格化浮点数写成

x=±(d0.d1dp1)ββe,d00.x=\pm(d_0.d_1\cdots d_{p-1})_\beta\,\beta^e, \qquad d_0\ne0.

舍入到最近可表示数记为 fl(x)\operatorname{fl}(x)。若中间结果既不溢出也不落入渐进下溢区,则

fl(x)=x(1+δ),δu,\operatorname{fl}(x)=x(1+\delta), \qquad |\delta|\le u,

其中 uu 称为单位舍入误差。对“舍入到最近值”的系统, u=12β1pu=\tfrac12\beta^{1-p}

在二进制 IEEE 风格系统中,有限数还包括带符号零、靠近零的次正规数,并与正负无穷和 NaN 一起规定异常结果。次正规数让绝对间距逐渐缩小,避免在最小规格化数处突然跳到零;但此时统一的相对误差界可能失效。溢出得到无穷或触发环境规定的异常,NaN 会沿后续运算传播。因此,舍入模型的使用前提必须写明“结果处于正常有限范围”。

常见的两个量不应混淆。uu 是一次正确舍入的最大相对误差;从 11 到下一个更大的浮点数的距离常记为 εmach\varepsilon_{\mathrm{mach}}。二进制舍入到最近值时,

u=2p,εmach=21p=2u.u=2^{-p}, \qquad \varepsilon_{\mathrm{mach}}=2^{1-p}=2u.

binary64 有 p=53p=53,故 u=2531.11×1016u=2^{-53}\approx1.11\times10^{-16},而 11 的下一浮点邻居与它相差 2522.22×10162^{-52}\approx2.22\times10^{-16}。某些软件把后一个量也叫 machine epsilon,使用前应先核对定义。

对基本运算 {+,,×,/}\circ\in\{+,-,\times,/\},在精确结果正常且除数非零时,可写成

fl(xy)=(xy)(1+δ),δu.\operatorname{fl}(x\circ y)=(x\circ y)(1+\delta), \qquad |\delta|\le u.

这是对“每一次已舍入运算”的局部描述,不表示一长串运算只有一次误差。若 nn 个小舍入因子相乘且 nu<1nu<1,常用

j=1n(1+δj)=1+θn,θnγn:=nu1nu\prod_{j=1}^{n}(1+\delta_j)=1+\theta_n, \qquad |\theta_n|\le\gamma_n:=\frac{nu}{1-nu}

把累计效应压缩成一个界。它说明短计算通常只有 O(nu)O(nu) 的相对扰动,但也提醒我们:运算次数、数据尺度和相减结构同样重要。

例 1:四位十进制系统中的吸收

β=10\beta=10p=4p=4,采用舍入到最近值,则 u=5×104u=5\times10^{-4}。按给定次序计算

(10,000+1)10,000.(10,000+1)-10,000.

第一步的精确结果是 10,001=1.0001×10410,001=1.0001\times10^4。四位系统只能保留 1.000×1041.000\times10^4,所以

fl(10,000+1)=10,000.\operatorname{fl}(10,000+1)=10,000.

该步相对误差为 1/10,0019.999×105<u1/10,001\approx9.999\times10^{-5}<u,完全符合局部舍入模型。第二步却得到零,而原表达式等于一。小数 11 被大数的间距“吸收”,随后相减使丢失的信息无法恢复。局部误差很小并不保证最终相对误差很小。

消去误差:相减暴露已有误差

两个接近数相减时,前导有效位相互抵消,结果本身很小。减法的最终舍入未必特别大;危险在于两个操作数此前各自携带的绝对误差,在小差值的尺度上会成为巨大的相对误差。这叫消去性相减。它与精确恒等式中的正常约分不同,关键是被减数已经过近似计算。

例 2:等价改写避开平方根相减

在只能保留八位十进制有效数字的环境中,计算

f(x)=1+x1,x=108.f(x)=\sqrt{1+x}-1, \qquad x=10^{-8}.

精确的平方根约为 1.0000000051.000000005,舍入到八位后成为 1.00000001.0000000,直接相减给出 00。先作有理化则有

f(x)=x1+x+1.f(x)=\frac{x}{\sqrt{1+x}+1}.

即使分母按八位舍入为 2.00000002.0000000,也得到 5.0000000×1095.0000000\times10^{-9}。更精细的值为

1082.0000000054.9999999875×109.\frac{10^{-8}}{2.000000005} \approx4.9999999875\times10^{-9}.

两个公式在实数中完全等价,但第二个计算图没有把两个接近的一阶量相减,因此保留了小结果的有效数字。

常见修复方式不是盲目提高精度,而是先改写公式:用有理化替换相近根式之差,用 log1p 型计算替换 log(1+x)\log(1+x) 在小 xx 处的直接求值,用 expm1 型计算替换 ex1e^x-1,或把二次方程的一个根稳定算出后由根的乘积得到另一个根。提高精度能推迟信息丢失,结构改写才消除诱因。

条件数描述问题本身的敏感性

设精确问题为 y=f(x)y=f(x)。绝对条件数衡量输入绝对扰动造成的输出绝对变化;当 xxf(x)f(x) 非零时,相对条件数衡量相对变化。对可微标量函数,局部条件数为

κabs(x)=f(x),κrel(x)=xf(x)f(x).\kappa_{\mathrm{abs}}(x)=|f'(x)|, \qquad \kappa_{\mathrm{rel}}(x) =\left|\frac{x f'(x)}{f(x)}\right|.

于是小扰动满足一阶近似

f(x+Δx)f(x)f(x)κrel(x)Δxx.\frac{|f(x+\Delta x)-f(x)|}{|f(x)|} \approx \kappa_{\mathrm{rel}}(x)\frac{|\Delta x|}{|x|}.

条件数属于问题和所选输入、输出尺度,不属于某段程序。κ\kappa 大表示问题病态:即使算法完全稳定,输入末位的不确定性也可能变成输出前几位的变化。κ\kappa 适中而结果仍差,才应重点检查算法和实现。

例 3:靠近极点时输入扰动被放大

考虑 f(x)=1/(1x)f(x)=1/(1-x)。由

f(x)=1(1x)2,κrel(x)=x1x.f'(x)=\frac1{(1-x)^2}, \qquad \kappa_{\mathrm{rel}}(x)=\left|\frac{x}{1-x}\right|.

x=0.99x=0.99 处,条件数为 9999,精确输出是 100100。若输入变为 x^=0.9901\widehat x=0.9901,输入相对改变量为

0.00010.991.0101×104,\frac{0.0001}{0.99}\approx1.0101\times10^{-4},

f(x^)=1/0.0099101.0101f(\widehat x)=1/0.0099\approx101.0101,输出相对变化约 1.0101×1021.0101\times10^{-2}。放大倍数约为 100100,与局部条件数给出的预测一致。这里更换等价程序不能消除极点附近的敏感性;应提高输入精度、远离奇异参数,或诚实报告不确定度。

求和也有数据相关的条件数。若 s=ixi0s=\sum_i x_i\ne0,对各项独立相对扰动可用

κsum=ixiixi\kappa_{\mathrm{sum}}=\frac{\sum_i|x_i|}{|\sum_i x_i|}

衡量敏感性。大量正负项抵消时分母很小,求和问题本身就可能病态;改变求和次序、从小到大累加或补偿求和可以减少算法误差,却不能恢复输入中原本不存在的精度。

前向误差、后向误差与稳定性

前向误差与后向误差

对问题 y=f(x)y=f(x) 和算法输出 y^\widehat y,前向误差是

y^f(x),\|\widehat y-f(x)\|,

它直接比较答案。后向误差寻找最小输入扰动 Δx\Delta x,使

y^=f(x+Δx),\widehat y=f(x+\Delta x),

并度量 Δx\|\Delta x\|。若一个算法对每个允许输入都能把计算结果解释为邻近问题的精确解,且相对后向误差与 uu 同阶,就称它后向稳定。

后向稳定不等于前向误差一定小。局部上可概括为

相对前向误差κrel(x)×相对后向误差.\text{相对前向误差} \lesssim \kappa_{\mathrm{rel}}(x) \times \text{相对后向误差}.

因此,“稳定算法 + 良态问题”才导向高相对精度。稳定算法遇到病态问题时可能正确地显示数据不够精确;不稳定算法遇到良态问题时则会浪费本可保留的数字。

例 4:同一结果的前向与后向读法

f(x)=x2f(x)=x^2x=2x=2,算法返回 y^=4.0004\widehat y=4.0004。精确值为 44,相对前向误差是

4.000444=104.\frac{|4.0004-4|}{4}=10^{-4}.

把结果视为某个邻近输入的精确平方,取

x^=4.0004=21.00012.0000999975.\widehat x=\sqrt{4.0004} =2\sqrt{1.0001} \approx2.0000999975.

相对后向误差约为

x^224.9999×105.\frac{|\widehat x-2|}{2}\approx4.9999\times10^{-5}.

f(x)=x2f(x)=x^2 的相对条件数恒为 22,所以“条件数乘后向误差”约为 9.9998×1059.9998\times10^{-5},正好解释前向误差的量级。三种量各有角色:后向误差评价算法,条件数评价问题,前向误差评价最终答案。

在线性方程、最小二乘或微分方程离散中,残差常是可直接计算的后向信息,但残差不是误差本身。必须结合算子条件数、数据尺度和离散误差,才能把小残差转成“解接近”的结论。

一套可复用的误差检查顺序

面对数值表达式,可依次检查:先写清输入的绝对或相对不确定度;再找奇点、重根、接近抵消和极端尺度以判断条件数;然后画出实际计算顺序,标出每次舍入、下溢与溢出风险;最后用更高精度、等价公式或独立算法复算。只报告许多小数位不是精度证明,最后稳定的位数应由误差界或交叉核验支持。

练习

练习 1:区分单位舍入误差与相邻间距

一个十进制浮点系统有 p=5p=5 位有效数字并舍入到最近值。求 uu11 的下一较大浮点数到 11 的距离,并给出舍入 x=3.1415926x=3.1415926 的结果和相对误差上界。

查看提示
舍入到最近值时,单位舍入误差是 1 附近向上间距的一半。
查看解答

u=121015=5×105,εmach=1015=104.u=\tfrac12 10^{1-5}=5\times10^{-5}, \qquad \varepsilon_{\mathrm{mach}}=10^{1-5}=10^{-4}.

五位舍入给 fl(x)=3.1416\operatorname{fl}(x)=3.1416。实际相对误差为

3.14163.14159263.14159262.36×106<u.\frac{|3.1416-3.1415926|}{3.1415926} \approx2.36\times10^{-6}<u.
练习 2:稳定计算很小的根式差

9+h3\sqrt{9+h}-3 改写为避免消去的形式,并在 h=106h=10^{-6} 时计算到十位有效数字。

查看提示
乘以共轭式,把小量留在分子。
查看解答

乘以共轭式得到

9+h3=h9+h+3.\sqrt{9+h}-3 =\frac{h}{\sqrt{9+h}+3}.

代入 h=106h=10^{-6},分母约为 6.00000016676.0000001667,故

9.00000131.666666620×107.\sqrt{9.000001}-3 \approx1.666666620\times10^{-7}.

直接计算需要先保存两个约为 33 的数再相减,改写后则直接把小参数按正常尺度传递到结果。

练习 3:幂函数的相对条件数

f(x)=xnf(x)=x^nx0x\ne0,求相对条件数。若 n=20n=20 且输入相对误差不超过 3×1063\times10^{-6},估计输出的一阶相对误差上界。

查看提示
计算 xf'(x)/f(x),并注意 x 非零。
查看解答

因为 f(x)=nxn1f'(x)=n x^{n-1},所以

κrel(x)=xnxn1xn=n.\kappa_{\mathrm{rel}}(x) =\left|\frac{x n x^{n-1}}{x^n}\right|=|n|.

n=20n=20 时,一阶相对误差约不超过

20×3×106=6×105.20\times3\times10^{-6}=6\times10^{-5}.

这只描述输入扰动的传播;若计算 x20x^{20} 的算法另有舍入误差,还应再加入算法误差界。

练习 4:给极限公式选择稳定计算图

函数 g(x)=(1cosx)/x2g(x)=(1-\cos x)/x^2 在小 xx 处直接计算容易失去有效数字。给出等价稳定形式,求 x0x\to0 的极限,并说明为什么改写有效。

查看提示
使用半角公式 1cosx=2sin2(x/2)1-\cos x=2 \sin^{2}(x/2)
查看解答

由半角公式,

g(x)=2sin2(x/2)x2=12(sin(x/2)x/2)2.g(x)=\frac{2\sin^2(x/2)}{x^2} =\frac12\left(\frac{\sin(x/2)}{x/2}\right)^2.

因此

limx0g(x)=12.\lim_{x\to0}g(x)=\frac12.

直接形式先把两个接近 11 的数相减;改写后计算的是接近 11 的比值再平方,没有小差值承载全部有效数字。实际实现还可在极小 xx 时用专门的 sinc 近似。

前后联系与资料

课程 · 2012

Introduction to Numerical Analysis

Laurent Demanet

用于核对 M10 各类算法的误差阶、稳定性条件、停止准则和可复算例题。

打开官方来源

MIT 18.330 Introduction to Numerical Analysis 提供浮点误差、条件数和数值线性代数的课程脉络,可用于复核本章的舍入模型、稳定性术语与后续算法分析。

有限精度不是计算结束后附加的一条免责声明,而是算法设计的一部分。先区分问题是否病态,再判断算法是否稳定,最后用前向误差说明答案质量,才能知道应当改写公式、提高输入精度,还是接受问题本身无法提供更多可靠数字。