浮点数，误差估计与分析

众说周知，由于内存的有限，计算机只能通过近似的方法表示一般的实数。其中最为广泛使用的方法就是浮点数。本文介绍浮点数的表示方法，重点放在对其误差的分析和理解上。本文的诞生来自于个人科研需要的总结，势必无法面面俱到。

定义

所谓浮点数表示，就是一种编码(近似)实数的方式，它通过尾数和幂数来表示：

这里的 \(m\) 和 \(e\) 都是有限的（有固定的存储位数)，也就导致一些数字只能使用近似的方式存储。我们只关注基数为2的情况，即假设 \(\beta = 2\)。这里列出 IEEE 754 规定的一些浮点数类型以及表示尾数：

相对误差与 ULP

参考 [1] 1.2节。

由于误差天然的存在于浮点数运算之中，找到一种衡量误差的方式就尤为重要。一种显然的方式是相对误差，也就是 \(\frac{x - \hat{x}}{x}\).

而另一种衡量误差的方式源自于其表示内部。考虑 \(1.d_1d_2\dots d_m \times 2^e\), 这一表示对应的实数实际上可能是 \(1.d_1d_2\dots d_m 1 \times 2^e\), 而由于浮点数表示的性质，这一被忽略了。不难发现，这会导致 1 个单位的误差，也就是 unit in the last place (ULP), 在上述例子中也就是 \(2^{-m + e}\).

一个有趣的事实是如果 ULP 误差会 bound 住相对误差, 那么有相对误差 \(\frac{x - \hat{x}}{x} \le k \epsilon, \epsilon=2^{-m}\), 证明如下：

考虑正数 \(x\) 与其浮点数表示 \(1.d_1d_2\dots d_m \times 2^e\), ULP 为 \(0.00\dots01\), 那么有 \(\begin{aligned} |x - \hat{x}| \le k \epsilon &= k 2^{-m + e} \\ \frac{x-\hat{x}}{x} & \le \frac{k 2^{-m + e}}{x} \\ &\le \frac{2^{-m}}{1.d_1d_2\dots d_m} \le 2^{-m} \end{aligned}\) 这表明了 ULP 误差实际上限制的是相对误差而不是绝对误差。

我们列出常见浮点数类型的 1ULP 对应的误差 ( \(\epsilon\) ) [3].

误差分析

本节我们讨论在浮点数误差框架下，基本运算，复合运算以及求和与矩阵乘法的误差估计。

基本运算

我们实际使用的各种库对例如加法，乘法和基本的数学运算都有基于 ULP 的相对误差界限规范，例如 GNU C 库中要求这些运算的相对误差小于 2 ULP [4]，CUDA中也有类似规定 [5].

复合运算与条件数

在获得了基本运算的相对误差估计后，一个自然的问题就是考虑一系列运算下，最终结果的误差是如何被这些基本误差bound的。例如 \(f(x) = \frac{1 - \cos x}{x^2}\) 就是以下一组基本运算的结果。

y = x * x
z = cos(x)
k = 1 - z
t = k / y

假设基本运算可微，有如下结果。 \(\begin{aligned} \operatorname{Err}_{r e l}(f(x), f(x+\Delta x)) & =\frac{f(x+\Delta x)-f(x)}{f(x)} \\ & =\frac{f(x+\Delta x)-f(x)}{\Delta x} \cdot \frac{\Delta x}{f(x)} \\ & =\left(f^{\prime}(x)+\frac{f^{\prime \prime}(x+\theta \Delta x)}{2!} \Delta x\right) \cdot \frac{\Delta x}{f(x)}, \theta \in(0,1) \\ & =\frac{\Delta x}{x} \cdot\frac{x f^{\prime}(x)}{f(x)}+O\left((\Delta x)^2\right) \\ & =\operatorname{Err}_{r e l}(x, x+\Delta x) \cdot\frac{x f^{\prime}(x)}{f(x)}+O\left((\Delta x)^2\right) \end{aligned}\)

我们称 \(\Gamma(x) = \frac{x f'(x)}{f(x)}\) 为条件数, 也就是放大误差的倍数。直觉上来看，就是自变量的相对误差被放大了条件数倍。显然，利用这一结论，我们可以给出整个运算的误差估计。

求和

参考 [2] 4节 Summation

给定一列数字 \(a[0], a[1], ..., a[n-1]\), 求和是极其常见的操作，对于它的误差估计也很重要。最直观的求和方法如下：

s = 0
for i in range(len(a)):
    s += a[i]

也就是每次都使用求和的结果和数组里的下一个元素相加，称为 递归求和。事实上我们也可以将这些数字两两组合然后相加，还有其他等等方法。

无论如何相加，我们都可以用下面的算法来表示，也就是每次在所有数字（包括中间求和的结果）中选两个数字相加：

S = {a[0], ..., a[n-1] }.
while len(S) > 1
	x, y = 选择 S 中两个元素移除
	z = x + y
	S.add(z)
end

我们的误差分析也就基于上面的一般性算法，从而会偏高。

考虑第 \(i\) 次加法 \(T_i = x_{i_1} + y_{i_1}\), 满足： \(\hat{T_i} = \frac{x_{i_1} + y_{i_1}}{1 + \delta_i}, |\delta_i| \le u\) 这里的 \(u\) 就是加法的误差界？显然一次加法运算的局部误差为 \(\delta_i \hat{T}_i\), 全局误差就是局部误差的总和 \(E_n = S_n - \hat{S_n} = \sum_{i=1}^{n-1}\delta_i\hat{T_i}\) 可想而知，可能的最小误差界是 \(|E_n| \le u\sum_{i=1}^{n-1}|\hat{T_i}|\). 容易发现 \(|\hat{T}_i| \le \sum_{j=1}^n |x_j| + O(u)\), 所以我们可以得到误差 \(|E_n| \le (n-1) u \sum_{i=1}^n|x_i| + O(u^2)\)

而相对误差也就小于等于 \((n-1)u\).

这一分析实际上还给了我们一个启示：想要最小化求和的误差，就要最小化这些中间结果。然而，这一问题一般来说是 NP-hard 的。对于递归求和的情况，如果所有元素都非负的，显然最好的方法就是从最小到大加。

矩阵乘法

参考 [2] 3 节

最后我们来讨论矩阵乘法的误差。在此之前，我们需要对于向量内积的分析。

援引 [2].3.1 的结论 \(|x^T y - \widehat{x^T y} | \le \gamma_n\sum_i^n|x_iy_i| = \gamma_n|x|^T|y|, \gamma_n = \frac{nu}{1-nu}\) 如果 \(nu \le 0.01\), 可以加强到 \(1.01n u |x|^T|y|\). 这里的 \(u\) 如上，也就是 \(\delta_i\) 的最大值。

基于此，对于矩阵乘法 \(C = A B\) 我们有 \(|C - \hat{C}| \le \gamma_n |A| |B|.\)

奇怪的东西

参考文献（按重要程度排序）

[1] David Goldberg. 1991. What every computer scientist should know about floating-point arithmetic. ACM Comput. Surv. 23, 1 (March 1991), 5–48. https://doi.org/10.1145/103162.103163

[2] Accuracy and Stability of Numerical Algorithms, Second Edition (2002)

[3] https://en.wikipedia.org/wiki/Machine_epsilon

[4] https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html

[5] https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#standard-functions