计算机算术

目录

• Floating-point number systems
• Rounding error analysis
  • Rounding a single number
  • Binary floating-point operations
  • The propagation of rounding errors
• Accuracy and stability
  • Avoiding catastrophic cancellation
  • Backward stability and numerical stability
  • Condition numbers: scalar functions
  • Condition numbers: vector functions
  • Condition numbers: algorithms
  • Overall error of a computer solution
• Problem

 

Floating-point number systems

Definition 4.1. 位置数字系统的基 base or radix of a positional numeral system 是表示数字的唯一记号的数量.

Definition 4.3. 位 bit 是计算机基本单位，只能有 $0$ 和 $1$ 两种值.

Definition 4.4. 字节 byte 是计算机信息的单位，通常包含 $8$ 位；它是计算机储存的最小单位.

 

Definition 4.7 (Floating point numbers). 浮点数 floating point number (FPN) 具有形式

$$x = \pm m\times \beta^e \tag{4.1}$$

其中 $\beta$ 是基， $e\in[L,U]$ ，有效数字 significand or mantissa $m$ 具有形式

$$m = \left(d_0+\dfrac{d_1}{\beta}+\cdots+\dfrac{d_{p-1}}{\beta^{p-1}}\right) \tag{4.2}$$

其中整数 $d_i$ 满足 $\forall i\in[0,p-1],\ d_i\in[0,\beta-1]$ ； $d_0$ 和 $d_{p-1}$ 称为最大有效数字 most significant digit 和最小有效数字 least significant digit ， $m$ 的数字串 string of digits 是字符串 $d_0.d_1d_2\cdots d_{p-1}$ ，其中 $.d_1d_2\cdots d_{p-1}$ 称为 $m$ 的小数部分 fraction .

 

Algorithm 4.8. 十进制整数转为二进制的方法

• 每次除以 $2$ ，然后记录余数，直到结果为 $0$ ，从后向前排列余数

$$97 = ((((((0\times 2+1)\times 2 + 1)\times 2 + 0) \times 2 + 0)\times 2 + 0)\times 2 + 0)\times 2 + 1 = (1100001)_2$$

十进制小数转为二进制的方法

• 每次乘 $2$ ，检查整数部分是否不大于 $1$ ，如果是就记录 $1$ ，否则记录 $0$ ，直到结果为 $0$ ，从前向后排列

$$\dfrac{2}{3} = \left(\left(\dfrac{1}{3}+1\right)\times\dfrac{1}{2^2}+1\right) \times \dfrac{1}{2} = (0.1010\cdots)_2$$

写成括号嵌套的形式更为直观.

 

Definition 4.11 (FPN systems). 浮点数系统 floating point number system $\mathcal{F}$ 是有理函数 $\mathbb{Q}$ 的真子集，由 $(\beta,p,L,U)$ 描述.

 

Definition 4.12. 一个 FPN 是正规的 normalized 如果它的有效数字 $1\le m<\beta$ .

 

Definition 4.13. 非正规数 subnormal 是 FPNs 中 $e=L,\ m\in(0,1)$ 的数，正规 FPN 系统可以通过包含非正规数扩张 extended .

 

Definition 4.14 (IEEE standard 754-2019). 当前的 IEEE 754 标准中的单精度 single precision 和双精度 double precision 浮点数有三种二进制格式和两种十进制格式；这里主要考虑 $32$ 位比特： $1$ 位符号位 $+$ $8$ 位指数 $+$ $23$ 位有效数字，精度为 $24$ ，因为我们可以在正规二进制浮点系统中取 $d_0=1$ ，然后不用存储 $d_0$ .

 

Definition 4.16. 正规 FPN 系统 $\mathcal{F}$ 的机器精度 machine precision 是 $1.0$ 到 $\mathcal{F}$ 中下一个较大的 FPN 的距离

$$\epsilon_M = \beta^{1-p} \tag{4.3}$$

按照格式整数加小数部分一共有 $p$ 位有效数字，最后一位为 $\beta^{1-p}$ .

 

Definition 4.17. 下溢极限 UFL 和上溢极限 OFL 为

$$UFL(\mathcal{F}) = \min|\mathcal{F}\setminus\{0\}| = \beta^L\\ OFL(\mathcal{F}) = \max|\mathcal{F}| = \beta^U(\beta-\beta^{1-p}) \tag{4.4}$$

下溢极限是最小的非零浮点数，上溢极限则是最大的浮点数.

 

Corollary 4.19 (Cardinality of $\mathcal{F}$ ). 对于二进制正规 FPN 系统，有

$$\#\mathcal{F} = 2^p(U-L+1) + 1 \tag{4.5}$$

$Proof.$ 注意在二进制储存时不用存 $d_0$ ，而是在这里存放符号，因此有 $2\times 2^{p-1}$ 种有效数字，同时有 $U-L+1$ 种指数，再加上 $0$ 就是所有的正规浮点数.

 

Definition 4.20. 正规 FPN 系统的范围 range 是 $\mathbb{R}$ 的子集

$$\mathcal{R}(\mathcal{F}) = \{x:x\in\mathbb{R},\ UFL(\mathcal{F})\le |x|\le OFL(\mathcal{F})\} \tag{4.6}$$

它是两段区间，分别在正半轴和负半轴；注意 $\mathcal{F}$ 不能取得范围内的所有实数.

 

Definition 4.22. 两个正规 FPNs $a,b$ 在 $\mathcal{F}$ 中相邻 adjacent 当且仅当

$$\forall c\in\mathcal{F}\setminus\{a,b\},\quad |a-b|<|a-c|+|b-c| \tag{4.7}$$

也就是两个数之间没有其它浮点数.

 

Lemma 4.23. 令两个正规 FPNs $a,b$ 在 $\mathcal{F}$ 中相邻，满足 $|a|<|b|,\ ab>0$ ，则

$$\beta^{-1}\epsilon_M|a| < |a-b|\le\epsilon_M|a| \tag{4.8}$$

$Proof.$ 考虑 $a = m\times\beta^e > 0$ ，则有 $|a-b| = \epsilon_M\beta^e$ ，由 $1.0\le m<\beta$ 有

$$\beta^{-1}\epsilon_M < \dfrac{|a-b|}{|a|} \le \epsilon_M$$

即证.

 

Rounding error analysis

Rounding a single number

Definition 4.24 (Rounding). 舍入 Rounding 是映射 $\mathrm{fl}:\mathbb{R}\to\mathcal{F}\cup\{+\infty,-\infty,NaN\}$ ，默认舍入是舍入到最近 round to nearest 的浮点数；在距离相等 tie 的情况下， $\mathrm{fl}(x)$ 选择舍入到最后一位为偶数 round to even .

 

Definition 4.25. 舍入数 $\mathrm{fl}(x)$ 上溢 overflows 若 $|x|>OFL(\mathcal{F})$ ，此时 $\mathrm{fl}(x) = NaN$ ；下溢 underflows 若 $0<|x|<UFL(\mathcal{F})$ ，此时 $\mathrm{fl}(x) = 0$ ；扩张的 FPN 系统的下溢称为渐进下溢 gradual underflow .

 

Definition 4.26. $\mathcal{F}$ 的单位舍入 unit roundoff 是

$$\epsilon_u = \dfrac{1}{2}\epsilon_M = \dfrac{1}{2}\beta^{1-p} \tag{4.9}$$

Theorem 4.27 (Range of round-off errors). 对 $x\in\mathcal{R}(\mathcal{F})$ 有

$$\mathrm{fl}(x) = x(1+\delta) \tag{4.10},\quad |\delta| < \epsilon_u$$

$Proof.$ 将浮点范围的实数舍入到 $\mathcal{F}$ 中，则我们考虑相邻的浮点数 $x_L,x_R\in\mathcal{F},\ x_L\le x\le x_R$ ，相等的情况显然；否则

$$|\mathrm{fl}(x)-x| \le \dfrac{1}{2}|x_R-x_L| \le \epsilon_u\min(|x_L|,|x_R|) < \epsilon_u|x|$$

其中第一步由舍入到最近的特性得到，第二步由 Lemma 4.23 得到.

 

Theorem 4.28. 对 $x\in\mathcal{R}(\mathcal{F})$ 有

$$\mathrm{fl}(x) = \dfrac{x}{1+\delta},\quad |\delta| \le \epsilon_u \tag{4.11}$$

$Proof.$ 类似的有

$$|\mathrm{fl}(x)-x| \le \dfrac{1}{2}|x_R-x_L| \le \epsilon_u\min(|x_L|,|x_R|) \le \epsilon_u|\mathrm{fl}(x)|$$

注意最后一步可以取等式.

 

Binary floating-point operations

Definition 4.30 (Addition/subtraction of two FPNs). 对 $a,b\in\mathcal{F},\ a = M_a\times\beta^{e_a},\ b = M_b\times \beta^{e_b},\ M_a = \pm m_a,\ M_b = \pm m_b$ ，不妨设 $|a|\ge|b|$ ，在精度至少为 $2p$ 的寄存器中计算 $c=\mathrm{fl}(a+b)\in\mathcal{F}$ .

• 比较指数
  • 若 $e_a-e_b > p+1$ ，相差超过范围，直接返回 $c=a$
  • 否则 $e_c\gets e_a,\ M_b\gets M_b/\beta^{e_a-e_b}$ ，这里对 $M_b$ 右移即可
• 加法 $M_c\gets M_a+M_b$ 并在寄存器中进行舍入，精度为 $2p$
• 正规化
  • 若 $M_c = 0$ ，返回 $0$
  • 否则移位使 $M_c\in[1,\beta)$
• 检查范围
  • 上溢 $NaN$
  • 下溢 $0$
• 舍入到精度 $p$
• 返回 $c\gets M_c\times \beta^e_c$

需要注意，由于精度为 $p$ ，舍入在第 $p$ 位小数进行；对于加法，最多错开 $p$ 位相加，然后舍入；对于减法，最多错开 $p+1$ 位相减，特殊情况为

$$1.000 \times 10^{0}-9.000 \times 10^{-5} = 1.000 - 0.00009 = 0.9999100\to 9.999 \times 10^{-1},\quad p=4$$

上面错开 $5$ 位，这是因为减法可能使得 $d_0=0$ ，然后左移一位，因而舍入在原先的第 $p+1$ 位小数进行.

 

Lemma 4.34. 对 $a,b\in\mathcal{F},\ a+b\in\mathcal{R}(\mathcal{F})$ 有

$$\mathrm{fl}(a+b) = (a+b)(1+\delta),\quad |\delta|<\epsilon_u \tag{4.12}$$

$Proof.$ 倒数第二步的舍入误差占主要地位，因为第二步舍入精度为 $2p$ ，只有在错开 $p+1$ 位时才有误差；然后由 Theorem 4.27 即证.

 

Definition 4.35 (Multiplication of two FPNs). 对 $a,b\in\mathcal{F},\ a = M_a\times\beta^{e_a},\ b = M_b\times \beta^{e_b},\ M_a = \pm m_a,\ M_b = \pm m_b$ ，在精度至少为 $2p$ 的寄存器中计算 $c = \mathrm{fl}(ab)\in\mathcal{F}$ .

• 指数求和 $e_c\gets e_a+e_b$
• 乘法 $M_c\gets M_aM_b$ 并在寄存器中进行舍入
• 正规化：移位使 $M_c\in[1,\beta)$
• 检查范围
  • 上溢 $NaN$
  • 下溢 $0$
• 舍入到精度 $p$
• 返回 $c\gets M_c\times \beta^e_c$

两个 $p$ 位精度的浮点数相乘得到 $2p$ 位小数.

 

Lemma 4.37. 对 $a,b\in\mathcal{F},\ |ab|\in\mathcal{R}(\mathcal{F})$ 有

$$\mathrm{fl}(ab) = (ab)(1+\delta),\quad |\delta|<\epsilon_u \tag{4.13}$$

$Proof.$ 只有倒数第二步的舍入误差，由 Theorem 4.27 即证.

 

Definition 4.38 (Division of two FPNs). 对 $a,b\in\mathcal{F},\ a = M_a\times\beta^{e_a},\ b = M_b\times \beta^{e_b},\ M_a = \pm m_a,\ M_b = \pm m_b$ ，在精度至少为 $2p+1$ 的寄存器中计算 $c = \mathrm{fl}(a/b)\in\mathcal{F}$ .

• 如果 $m_b=0$ ，直接返回 $NaN$ ；指数作差 $e_c\gets e_a-e_b$
• 除法 $M_c\gets M_a/M_b$ 并在寄存器中进行舍入
• 正规化：移位使 $M_c\in[1,\beta)$
• 检查范围
  • 上溢 $NaN$
  • 下溢 $0$
• 舍入到精度 $p$
• 返回 $c\gets M_c\times \beta^e_c$

 

Lemma 4.39. 对 $a,b\in\mathcal{F},\ a/b\in\mathcal{R}(\mathcal{F})$ 有

$$\mathrm{fl}\left(\dfrac{a}{b}\right) = \dfrac{a}{b}(1+\delta),\quad |\delta|<\epsilon_u \tag{4.14}$$

$Proof.$ 当 $|M_a|=|M_b|$ 时，没有误差；考虑 $|M_a|>|M_b|$ ，则由 $|M_a|,|M_b|\in[1,\beta)$ ，尽可能让商更小，有

$$\beta > \beta-\epsilon_M > \beta-2\epsilon_M,\quad \left|\dfrac{M_a}{M_b}\right| \ge \dfrac{\beta-\epsilon_M}{\beta-2\epsilon_M} > 1+\beta^{-1}\epsilon_M$$

这意味着第三步正规化不需要进行；设有 $p+k$ 精度，则舍入单位

$$\dfrac{1}{2}\beta^{1-p-k} = \dfrac{1}{2}\beta^{1-p}\beta^{1-p}\beta^{p-1-k} = \epsilon_u\epsilon_M\beta^{p-1-k}$$

令 $k=p+1$ ，则舍入单位为 $\epsilon_u\epsilon_M\beta^{-2}$ ；记第二步舍入结果为 $M_{c_1}$ ，倒数第二步舍入结果为 $M_{c_2}$ ，从而

$$\begin{aligned} M_{c_2} &= M_{c_1} + \delta_2,\quad |\delta_2| < \epsilon_u\\ &= \dfrac{M_a}{M_b} + \delta_1 + \delta_2,\quad |\delta_1|<\epsilon_u\epsilon_M\beta^{-2}\\ &= \dfrac{M_a}{M_b}(1+\delta)\\ |\delta| &= \left|\dfrac{\delta_1+\delta_2}{M_a/M_b}\right| < \dfrac{\epsilon_u(1+\epsilon_M\beta^{-2})}{1+\epsilon_M\beta^{-1}} < \epsilon_u \end{aligned}$$

考虑 $|M_a|<|M_b|$ ，则类似的有

$$\left|\dfrac{M_a}{M_b}\right| \le \dfrac{\beta-2\epsilon_M}{\beta-\epsilon_M} = 1 - \dfrac{\epsilon_M}{\beta-\epsilon_M} < 1 - \beta^{-1}\epsilon_M$$

因此第三步正规化一定会进行， $M_{c_1}$ 需要左移一位，因此

$$M_{c_1} = \dfrac{M_a}{M_b} + \delta_1,\quad |\delta_1|<\epsilon_u\epsilon_M\beta^{-2}\\ M_{c_2} = \beta M_{c_1} + \delta_2 = \beta\dfrac{M_a}{M_b}\left(1+\dfrac{\beta\delta_1+\delta_2}{\beta M_a/M_b}\right),\quad |\delta_2|<\epsilon_u$$

我们考虑分母上的部分，进行缩放则有

$$\beta\left|\dfrac{M_a}{M_b}\right| \ge \dfrac{\beta}{\beta-\epsilon_M} = 1 + \dfrac{\epsilon_M}{\beta-\epsilon_M} > 1 + \beta^{-1}\epsilon_M$$

代入得到

$$|\delta| = \left|\dfrac{\beta\delta_1+\delta_2}{\beta M_a/M_b}\right| < \dfrac{\epsilon_u\epsilon_M\beta^{-1}+\epsilon_u}{1+\beta^{-1}\epsilon_M} = \epsilon_u$$

即证.

 

Theorem 4.40 (Model of machine arithmetic). 记 $\mathcal{F}$ 为精度 $p$ 的正规 FPN 系统，则对任意运算 $\odot = +,-,\times,/$ ，我们有

$$\forall a,b\in\mathcal{F},\ a\odot b\in\mathcal{R}(\mathcal{F}),\quad \mathrm{fl}(a\odot b) = (a\odot b)(1+\delta),\quad |\delta|<\epsilon_u \tag{4.15}$$

当且仅当在 $2p+1$ 精度的寄存器中计算.

 

The propagation of rounding errors

Theorem 4.41. 若 $\forall i=0,1,\cdots,n,\ a_i\in\mathcal{F},\ a_i>0$ ，则

$$\mathrm{fl}\left(\sum_{i=0}^na_i\right) = (1+\delta_n)\sum_{i=0}^na_i \tag{4.16}$$

其中 $|\delta_n|<(1+\epsilon_u)^n-1\approx n\epsilon_u$ .

 

Theorem 4.44. 对给定的 $\mu\in\mathbb{R}^+$ 及正整数 $n\le \lfloor \frac{\ln2}{\mu}\rfloor$ ，设 $|\delta_i|\le \mu,\ i=1,2,\cdots,n$ ，则

$$1-n\mu\le \prod_{i=1}^n(1+\delta_i) \le 1+n\mu+(n\mu)^2 \tag{4.17}$$

或者等价地对 $I_n=[-\frac{1}{1+n\mu},1]$ 有

$$\exists\theta\in I_n,\quad \mathrm{s.t.}\quad \prod_{i=1}^n(1+\delta_i) = 1+\theta(n\mu+n^2\mu^2) \tag{4.18}$$

$Proof.$ 由 $|\delta_i|\le \mu,\ i=1,2,\cdots,n$ 可得

$$(1-\mu)^n \le \prod_{i=1}^n(1+\delta_i) \le (1+\mu)^n \le e^{n\mu}$$

我们将 $f(\mu) = (1-\mu)^n$ 在 $\mu=0$ 泰勒展开的拉格朗日余项有

$$(1-\mu)^n = 1-n\mu+\frac{1}{2}n(n-1)\theta^2 \ge 1-n\mu$$

从而左边不等式得证，而

$$e^{x} = 1+x+\dfrac{x^2}{2}+\dfrac{x^3}{3!} +\cdots \le 1+x+\dfrac{x^2}{2}e^x$$

令 $x=n\mu\le \ln2$ ，则有

$$e^{n\mu} \le 1+n\mu+(n\mu)^2$$

右边不等式得证；最后，由于 $\prod_{i=1}^n(1+\delta_i)$ 在连续函数 $f(\tau) = 1+\tau(1+n\mu)n\mu,\ \tau\in I_n$ 的值域中，由中值定理即证.

 

Accuracy and stability

Avoiding catastrophic cancellation

Definition 4.45. 令 $\hat{x}$ 为 $x\in\mathbb{R}$ 的逼近，则逼近精度可由绝对误差 absolute error

$$E_{abs}(\hat{x}) = |\hat{x}-x| \tag{4.18}$$

或者相对误差 relative error

$$E_{rel}(\hat{x}) = \dfrac{|\hat{x}-x|}{|x|} \tag{4.19}$$

进行度量.

 

Definition 4.46. 对于通过 $\hat{y} = \hat{f}(x)$ 逼近 $y=f(x)$ ，向前误差 forward error 是逼近 $y$ 的相对误差，向后误差 backward error 是满足 $f(\hat{x})=\hat{f}(x)$ 的用 $\hat{x}$ 逼近 $x$ 的相对误差的最小值.

Note. 关于向后误差，我们原本想要的是 $y=f(x)$ ，但必须通过 $\hat{y} = \hat{f}(x)$ 来近似，也就是说我们计算的是 $\hat{f}(x)$ ，有 $f(\hat{x}) = \hat{f}(x)$ ，即我们实际计算的是 $f$ 在 $\hat{x}$ 处的值。总的来说，可以看作我们实际上是在用 $f$ 在 $\hat{x}$ 处的值逼近 $f$ 在 $x$ 处的值，向后误差考虑的就是用于逼近的 $\hat{x}$ 关于实际的 $x$ 的相对误差。

事实上，我们不需要考虑结果的误差，只要 $\hat{x}$ 离 $x$ 足够近，也可以认为这种逼近是有效的。因为 $f(\hat{x}) = \hat{f}(x)$ ，所以我们保证了用 $\hat{f}(x)$ 逼近 $f(\hat{x})$ 是完全准确的，当 $\hat{x}$ 离 $x$ 足够近，$f(\hat{x})$ 离 $f(x)$ 也足够近， $\hat{f}(x)$ 可以看做是对 $f(x)$ 的良好逼近.

 

Definition 4.47 (Accuracy). 计算 $y=f(x)$ 的算法 $\hat{y} = \hat{f}(x)$ 是精确的 accurate 如果它的向前误差对任意的 $x$ 都很小

$$\exist c>0,\quad \mathrm{s.t.}\quad \forall x\in \mathrm{dom}(f),\ E_{rel}(\hat{f}(x)) = \left|\dfrac{\hat{f}(x)-f(x)}{f(x)}\right| \le c\epsilon_u \tag{4.20}$$

 
Note. 糟糕的消去：我们知道乘法和除法是精确的，但是加法和减法不一定是精确的；事实上，通过计算

$$\mathrm{fl}(\mathrm{fl}(x)\odot \mathrm{fl}(y)) = (1+C(x,y))(x\odot y)$$

对于乘法和除法 $C(x,y) = C$ ；而对于加法和减法 $C(x,y)\to \infty$ 当 $x+y\to 0$ .

 

Theorem 4.49 (Loss of most significant digits). 设 $x,y\in\mathcal{F},\ x>y>0$ ，且

$$\beta^{-t}\le 1-\dfrac{y}{x} \le \beta^{-s} \tag{4.21}$$

则在计算减法 $x-y$ 时最大有效数字的位数至多丢失 $t$ ，至少丢失 $s$ .

$Proof.$ 我们写成浮点格式 $x=m_x\times \beta^n,\ y=m_y\times \beta^m,\ 1\le m_x,m_y< \beta$ ，根据减法过程，首先对较小的 $y$ 进行移位

$$y = (m_y\times \beta^{m-n})\times\beta^n \Rightarrow x-y = (m_x-m_y\times \beta^{m-n})\times\beta^n$$

把相减后的有效数字记为 $m_{x-y}$ ，则m

$$m_{x-y} = m_x-m_y\times \beta^{m-n} = m_x\left(1-\dfrac{m_y\times \beta^m}{m_x\times \beta^n}\right) = m_x\left(1-\dfrac{y}{x}\right)$$

根据题中的不等式以及 $m_x$ 的条件有

$$\beta^{-t}\le m_{x-y} < \beta^{1-s}$$

因此进行正规化时，要左移至多 $t$ 位，至少 $s$ 位.

 

Note. 当 $x,y$ 是两个非常接近的数时，进行减法操作就会出现上面的问题，因此要尽量避免这种情况出现.

 

Backward stability and numerical stability

Definition 4.52 (Backward stability). 计算 $y=f(x)$ 的算法 $\hat{y} = \hat{f}(x)$ 是向后稳定的 backward stable 若

$$\exist c>0,\quad \mathrm{s.t.}\quad \forall x\in \mathrm{dom}(f),\ \exists\hat{x}\in\mathrm{dom}(f),\quad \mathrm{s.t.}\quad \hat{f}(x)=f(\hat{x})\ \Rightarrow\ E_{rel}(\hat{x}) \le c\epsilon_u \tag{4.20}$$

Definition 4.53. 计算 $y=f(x_1,x_2)$ 的算法 $\hat{y} = \hat{f}(x_1,x_2)$ 是向后稳定的 backward stable 若

$$\exist c_1,c_2>0,\quad \mathrm{s.t.}\quad \forall (x_1,x_2)\in \mathrm{dom}(f),\ \exists(\hat{x}_1,\hat{x}_2)\in\mathrm{dom}(f),\quad \mathrm{s.t.}\\ \hat{f}(x_1,x_2)=f(\hat{x}_1,\hat{x}_2)\ \Rightarrow\ E_{rel}(\hat{x}_1) \le c_1\epsilon_u,\ E_{rel}(\hat{x}_2) \le c_2\epsilon_u \tag{4.21}$$

也就是说对任何 $x$ 的向后误差都足够小.

 

Lemma 4.54. 对 $f(x_1,x_2) = x_1-x_2,\ x_1,x_2\in\mathcal{R}(\mathcal{F})$ ，则算法 $\hat{f}(x_1,x_2) = \mathrm{fl}(\mathrm{fl}(x_1)- \mathrm{fl}(x_2))$ 是向后稳定的.

$Proof.$ 直接计算

$$\hat{f}(x_1,x_2) = (x_1(1+\delta_1)-x_2(1+\delta_2))(1+\delta_3) = x_1(1+\delta_1+\delta_3+\delta_1\delta_3)-x_2(1+\delta_2+\delta_3+\delta_2\delta_3)$$

因此分别有

$$\hat{x}_1 = (1+\delta_1+\delta_3+\delta_1\delta_3)x_1,\quad \hat{x}_2 = (1+\delta_2+\delta_3+\delta_2\delta_3)x_2$$

从而向后稳定.

 

Definition 4.56. 计算 $y=f(x)$ 的算法 $\hat{y} = \hat{f}(x)$ 是稳定的 stable 或数值稳定 numerically stable 当且仅当

$$\exist c,c_f>0,\quad \mathrm{s.t.}\quad \forall x\in \mathrm{dom}(f),\ \exists\hat{x}\in\mathrm{dom}(f),\quad \mathrm{s.t.}\quad \left|\dfrac{\hat{f}(x)-f(\hat{x})}{f(\hat{x})}\right| \le c_f\epsilon_u,\ E_{rel}(\hat{x}) \le c\epsilon_u \tag{4.22}$$

注意这里第一个条件是关于 $f(\hat{x})$ 的相对误差，第二个条件是关于 $\hat{x}$ 的相对误差.

 

Lemma 4.57. 如果一个算法向后稳定，则它是数值稳定的.

$Proof.$ 由向后稳定，有 $\hat{f}(x)=f(\hat{x})$ ，故第一个条件显然；第二个条件显然成立.

 

Condition numbers: scalar functions

Definition 4.59. 函数 $y=f(x)$ 的相对条件数 condition number 是输入的微小变化对输出的相对变化的度量

$$C_f(x) = \left|\dfrac{xf^{\prime}(x)}{f(x)}\right| \tag{4.23}$$

Definition 4.60. 条件数小的问题称为是良态的 well-conditioned 否则称为是病态的 ill-conditioned .

 

Lemma 4.63. 考虑在单根 $r$ 附近求解 $f(x)=0$ ，即有 $f(r)=0,\ f^{\prime}(r)\neq 0$ ，设有微小扰动 $F=f+\epsilon g,\ f,g\in\mathcal{C}^2,\ g(r)\neq 0$ 且满足 $|\epsilon g^{\prime}(r)|\ll |f^{\prime}(r)|$ ，则 $F$ 有根 $r+h$ ，其中

$$h\approx -\epsilon\dfrac{g(r)}{f^{\prime}(r)} \tag{4.24}$$

$Proof.$ 直接泰勒展开

$$0 = F(r+h) = f(r) + hf^{\prime}(r) + \epsilon [g(r)+hg^{\prime}(r)] + O(h^2)$$

从而有

$$h\approx -\epsilon\dfrac{g(r)}{f^{\prime}(r)+\epsilon g^{\prime}(r)} \approx -\epsilon\dfrac{g(r)}{f^{\prime}(r)}$$

即证.

 

Example 4.64 (Wilkinson). 定义

$$f(x) = \prod_{k=1}^p(x-k),\quad g(x) = x^p$$

考虑根 $x=p$ 受到扰动 $f+\epsilon g$ 后的影响.

由上面的引理，我们有

$$h\approx -\epsilon\dfrac{g(r)}{f^{\prime}(r)} = -\epsilon \dfrac{p^p}{(p-1)!}$$

当 $p$ 很大时，扰动越来越大，因此求解高阶多项式的根几乎不可能.

 

Condition numbers: vector functions

Definition 4.65. 向量函数 $\mathbf{f}:\mathbb{R}^m\to\mathbb{R}^n$ 的条件数为

$$\mathrm{cond}_{\mathbf{f}}(\mathbf{x}) = \dfrac{\|\mathbf{x}\|\|\nabla\mathbf{f}\|}{\|\mathbf{f}(\mathbf{x})\|} \tag{4.25}$$

其中 $\|\cdot\|$ 表示欧几里得范数.

 

Definition 4.68. 向量函数 $\mathbf{f}:\mathbb{R}^m\to\mathbb{R}^n$ 的分量条件数 component wise condition number 为

$$\mathrm{cond}_{\mathbf{f}}(\mathbf{x}) = \|A(\mathbf{x})\| \tag{4.26}$$

其中矩阵 $A(\mathbf{x}) = [a_{ij}(\mathbf{x})]$ ，每一个部分为

$$a_{ij}(\mathbf{x}) = \left|\dfrac{x_j\frac{\partial f_i}{\partial x_j}}{f_i(\mathbf{x})}\right| \tag{4.27}$$

也就是说 $a_{ij}$ 是 $f_i(\mathbf{x})$ 关于 $x_j$ 的条件数.

 

Definition 4.70. 希尔伯特矩阵 Hilbert matrix $H_n\in\mathbb{R}^{n\times n}$ 是

$$h_{i,j} = \dfrac{1}{i+j-1} \tag{4.28}$$

 
Definition 4.72. 范德蒙德矩阵 Vandermonde matrix $V_n\in\mathbb{R}^{n\times n}$ 是

$$v_{i,j} = t_j^{i-1} \tag{4.29}$$

其中 $t_1,\cdots,t_n$ 是参数.

 

Condition numbers: algorithms

Definition 4.74. 向量函数 $\mathbf{f}:\mathbb{R}^m\to\mathbb{R}^n$ 通过算法 $\mathbf{f}_A:\mathcal{F}^m\to\mathcal{F}^n$ 逼近，设

$$\forall \mathbf{x}\in\mathcal{F}^m,\ \exist\mathbf{x}_A\in\mathbb{R}^m,\quad \mathrm{s.t.}\quad \mathbf{f}_A(\mathbf{x}) = \mathbf{f}(\mathbf{x}_A) \tag{4.30}$$

则算法 $\mathbf{f}_A$ 的条件数定义为

$$\mathrm{cond_A(\mathbf{x})} = \dfrac{1}{\epsilon_u}\inf_{\{\mathbf{x}_A\}}\dfrac{\|\mathbf{x}_A-\mathbf{x}\|}{\|\mathbf{x}\|} \tag{4.31}$$

也就是向后误差的下确界.

 

Theorem 4.76. 设光滑函数 $f:\mathbb{R}\to\mathbb{R}$ 由算法 $A:\mathcal{F}\to\mathcal{F},\ f_A(x) = f(x)(1+\delta(x)),\ |\delta(x)|\le\varphi(x)\epsilon_u$ 逼近，若 $\mathrm{cond}_f(x)$ 有界非零，则有

$$\forall x\in\mathcal{F},\quad \mathrm{cond}_A(x)\le \dfrac{\varphi(x)}{\mathrm{cond}_f(x)} \tag{4.32}$$

$Proof.$ 设有 $x_A$ 满足 $f(x_A)=f_A(x)$ ，将 $x_A$ 写作 $x(1+\epsilon_A)$ ，则

$$f(x)(1+\delta) = f(x_A) = f(x(1+\epsilon_A)) = f(x) + x\epsilon_Af^{\prime}(x) + O(\epsilon_A^2)\ \Rightarrow\ \epsilon_A = \dfrac{f(x)}{xf^{\prime}(x)}\delta$$

计算相对误差

$$\left|\dfrac{x_A-x}{x}\right| = |\epsilon_A| = \left|\dfrac{f(x)}{xf^{\prime}(x)}\right||\delta(x)|\ \Rightarrow\ \dfrac{1}{\epsilon_u}\left|\dfrac{x_A-x}{x}\right| = \left|\dfrac{\delta(x)}{\epsilon_u\mathrm{cond}_f(x)}\right| \le \dfrac{\varphi(x)}{\mathrm{cond}_f(x)}$$

对 $x_A$ 取下确界， $x$ 取上确界即证.

 

Overall error of a computer solution

Theorem 4.79. 考虑正规 FPN 求解数学问题

$$\mathbf{f}:\mathbb{R}^m\to\mathbb{R}^n,\quad \mathbf{y} = \mathbf{f}(\mathbf{x}) \tag{4.33}$$

将计算机的输入输出表示为

$$\mathbf{x}^*\approx \mathbf{x},\quad \mathbf{y}^*_A = \mathbf{f}_A(\mathbf{x}^*) \tag{4.34}$$

其中 $\mathbf{f}_A$ 是逼近 $\mathbf{f}$ 的算法，则逼近 $\mathbf{y}$ 与 $\mathbf{y}^*$ 的相对误差有界

$$E_{rel}(\mathbf{y}^*_A) \lessapprox E_{rel}(\mathbf{x}^*_A)\mathrm{cond}_{\mathbf{f}}(\mathbf{x}) + \epsilon_u\mathrm{cond}_{\mathbf{f}}(\mathbf{x}^*)\mathrm{cond}_{A}(\mathbf{x}^*) \tag{4.35}$$

$Proof.$ 证明省略.

 

Problem

$\mathrm{III}.$ 令 $x=\beta^e,\ e\in\mathbb{Z},\ L<e<U$ 为正规 FPN ，而 $x_L,x_R$ 为两个相邻的正规 FPN 满足 $x_L<x<x_R$ ，则 $x_R-x=\beta(x-x_L)$ .

$Proof.$ 由 $x = 1 \times \beta^e$ ，则有 $x_R = (1 + 1 \times \beta^{1-p}) \times \beta^e$ ，从而 $x_R - x = \beta^{1-p} \times \beta^e = \beta^{e+1-p}$ ，而

$$x_L = (\beta - 1) \times (1 + \beta^{-1} + \beta^{-2} + \cdots + \beta^{1-p}) \times \beta^{e-1}$$

从而有

$$\begin{aligned} x - x_L &= (\beta - (\beta - 1) \times (1 + \beta^{-1} + \beta^{-2} + \cdots + \beta^{1-p})) \times \beta^{e-1}\\ &= \left(\beta - (\beta - 1) \times \dfrac{1-\beta^{-p}}{1-\beta^{-1}}\right) \times \beta^{e-1}\\ &= (\beta - (1-\beta^{-p}) \times \beta) \times \beta^{e-1}\\ &= \beta^{e-p} \end{aligned}$$

故 $x_R - x = \beta(x - x_L)$ .

 

$\mathrm{IV}.$ 在 IEEE 754 单精度下进行舍入.

$$\begin{aligned} \dfrac{3}{5} = (1.00110011\cdots)_2 \times 2^{-1}\\ x_L = (1.0011\cdots001)_2 \times 2^{-1}\\ x_R = (1.0011\cdots010)_2 \times 2^{-1} \end{aligned}$$

且有

$$\begin{aligned} x - x_L &= \dfrac{3}{5} \times 2^{-24}\\ x_R - x_L &= 2^{-24}\\ x_R - x &= (x_R - x_L) - (x - x_L) = \dfrac{2}{5} \times 2^{-24} \end{aligned}$$

故 $\mathrm{fl}(x) = x_R = (1.0011\cdots010)_2 \times 2^{-1}$ ，

$$\left|\dfrac{x_R-x}{x}\right| = \dfrac{2}{3} \times 2^{-24}$$

注意舍入精度为 $24$ ，即保留到小数点后 $23$ 位.

 

$\mathrm{V}.$ 如果 IEEE 754 单精度不再舍入到最近，而是直接截断，求单位舍入.

机器精度 $\epsilon_M = 2^{-23}$ ，由于直接截断，则舍入的最大相对误差

$$\epsilon_u = \max\left|\dfrac{\mathrm{fl}(x)-x}{x}\right| < \max\left|\dfrac{2^{-23}}{x}\right| = 2^{-23}$$

从而 $\epsilon_u = 2^{-23}$ .

 

$\mathrm{VI}.$ 在计算 $1-\cos \frac{1}{4}$ 的精度损失.

根据定理 Theorem 4.49 我们计算

$$1-\dfrac{y}{x} = 1 - \dfrac{\cos \frac{1}{4}}{1} = 1-\cos \frac{1}{4}$$

有范围

$$0.015625 = 2^{-6} < 1 - \cos\frac{1}{4} < 2^{-5} = 0.03125$$

故精度损失 $6$ 位.

 

$\mathrm{VII}.$ 两种避免 $1-\cos x$ 的精度损失的方法.

• 计算泰勒展开

$$\begin{aligned} 1 - \cos x &= 1 - \left(1 - \dfrac{x^2}{2} + \dfrac{x^4}{4!} - \dfrac{x^6}{6!} \cdots \right)\\ &= \dfrac{x^2}{2} - \dfrac{x^4}{4!} + \dfrac{x^6}{6!} - \cdots \end{aligned}$$

• 计算

$$1 - \cos x = \dfrac{1 - \cos^{2}x}{1 + \cos x} = \dfrac{\sin^{2}x}{1 + \cos x}$$

用多个乘除法替代.

 

$\mathrm{IX}.$ 算法 $A$ 估计 $f(x)=1-e^{-x}$ ，其中指数函数的相对误差为机器精度，对 $[0,1]$ 估计 $\mathrm{cond}_A(x)$ .

利用估计

$$\forall x\in\mathcal{F},\quad \mathrm{cond}_A(x)\le \dfrac{\varphi(x)}{\mathrm{cond}_f(x)}$$

不妨假设估计为 $f_A(x)$ ，在计算时忽略 $\delta_i\delta_i$ ，则有

$$\begin{aligned} f_A(x) &= (1-e^{-x}(1+\delta_1))(1+\delta_2)\\ &= 1 + \delta_2 - e^{-x}(1+\delta_1+\delta_2 + \delta_1\delta_2)\\ &\approx f(x)\left(1+\dfrac{\delta_2-e^{-x}(\delta_1+\delta_2)}{1-e^{-x}}\right)\\ &= f(x)(1+\delta(x)) \end{aligned}$$

考虑 $|\delta(x)|<\epsilon_u\varphi(x)$ ，从而有

$$\varphi(x) = \dfrac{1-2e^{-x}}{1-e^{-x}},\quad \mathrm{cond}_f(x) = \left|\dfrac{xf^{\prime}(x)}{f(x)}\right| = \dfrac{x}{e^x-1}$$

其中 $\mathrm{cond}_f(x)$ 有界非零，代入估计式

$$\mathrm{cond}_A(x) \le \dfrac{e^x-2}{x}$$

从而在 $x=0$ 附近算法条件数趋于无穷.

 

$\mathrm{X}.$ 寻找多项式的根

$$q(x) = \sum_{i=0}^na_ix^i,\quad a_n=1,\ a_0\neq 0,\ a_i\in\mathbb{R}$$

可以考虑向量函数 $f:\mathbb{R}^n\to\mathbb{C}$

$$r=f(a_0,a_1,\cdots,a_{n-1})$$

推导 $f$ 基于 $1$ 范数的分量条件数.

令 $\mathbf{a} = (a_0,a_1,\cdots,a_{n-1})^{T}$ ，根据 $q(x)$ 得到的函数 $f$ 为

$$f(\mathbf{a}) = \sum_{i=0}^{n}a_ix^i$$

条件矩阵 $A(\mathbf{a})$ 是 $1\times n$ 矩阵

$$A(\mathbf{a}) = \dfrac{1}{|f(\mathbf{a})|}\left(|a_0|, |a_1x|, \cdots, \left|a_{n-1}x^{n-1}\right|\right)$$

从而有

$$\begin{aligned} \mathrm{cond}_f(\mathbf{a}) &= \dfrac{1}{|f(\mathbf{a})|}\left\|\left(|a_0|, |a_1x|, \cdots, \left|a_{n-1}x^{n-1}\right|\right)\right\|_1\\ &= \dfrac{\max\left(|a_0|, |a_1x|, \cdots, \left|a_{n-1}x^{n-1}\right|\right)}{\left|\sum_{i=0}^{n}a_ix^i\right|} \end{aligned}$$

对于 Wilkinson 情况，

$$f(x) = \prod_{k=1}^p(x-k) = x^p + \sum_{i=0}^{p-1}b_ix^i,\quad g(x) = x^p$$

则求 $f+\epsilon g$ 的根等价于

$$f+\epsilon g= x^p + \sum_{i=0}^{p-1}b_ix^i + \epsilon x^{p} = 0\ \Leftrightarrow\ x^p + \sum_{i=0}^{p-1}c_ix^i = 0,\ c_i = \dfrac{b_i}{1+\epsilon}$$

也就是说，求解 $f+\epsilon g$ 的根实际上相当于对 $f$ 的系数做出相对扰动

$$\dfrac{\left|c_i-b_i\right|}{|b_i|} = \dfrac{\epsilon}{1+\epsilon}$$

得到的首一多项式求根，应用之前首一多项式系数的条件数公式，它在 $x=p$ 附近条件数趋于无穷，因此系数受到的扰动很大，它的根受到的扰动也很大.

 

$\mathrm{XII}.$ 在 IEEE 754 单精度 FPN 中，对区间 $[128,129]$ 应用二分法，则绝对精度能否 $<10^{-6}$ ？

要满足 $|x-x^*| < 10^{-6}$ ，但有整数部分 $128 = 2^7$ ，而

$$2^{-20} < 1 \times 10^{-6} < 2^{-19}$$

从而整数部分 $8$ 位，小数部分 $19$ 位，正规化后需要至少 $27$ 位，因此不能达到要求.