浮点运算的基本公理

如果您已经做了足够长时间的软件工程师，您很可能见过这个浮点背信弃义的例子：

我们理解这是因为浮点数只有64位精度存储，不能代表整个实数行。此外，当我们使用这些浮点数执行运算时，其表示形式中固有的错误可能会累积并成倍增加。这个故事的寓意是，永远不要用浮点数来代表金钱。

至少，这是金融应用的寓意所在。在广场，我们用了一大笔钱，我们过得很好。然而，对于大多数有充分理由使用浮点的应用程序来说，这并不是故事的结束。如果浮点数像我曾经认为的那样不可预测、不可靠，我们就不能用数字方法求解微分方程或线性系统，也不能登陆月球。相反，有一门关于浮点误差的科学，它是一般数字误差科学的一部分，它试图驯服和理解误差在我们的计算过程中发生了什么。一般说来，数字误差是可以相当精确地量化的东西，正如我们将看到的。在此过程中，我们将了解浮点算术的基本公理，至少，它听起来比“每个开发人员都应该知道浮点数的10件事”要酷得多。

从根本上说，浮点数的误差是由“舍入”问题引起的。就像在基数10中一样，我们不能表示数字\(1/3\)而不将其舍入到某个位置：

在以2为基数的情况下，我们不能在不舍入的情况下表示许多数字。当然，有些数字我们可以准确地表示出来。例如，然后是1号。或范围\((-2^{53}，2^{53})\)中的任意整数。同样值得注意的是，可以准确地表示许多分数：

但是，像\(1/10\)这样的数字，就像基数10中的\(1/3\)一样，必须截断以适合尾数的24位或53位。当我们在控制台或源代码中输入0.1时，实际存储的值与0.1略有不同。根据这个出色的IEEE-754浮点转换器，实际存储的浮点数(对于64位浮点)是。

所以我们计算的初始输入是有缺陷的！我们不是在计算0.1*3，我们实际上是在计算。

这个误差有多大？我们可以通过计算有效数字\(00055\点0.100\)之间的零来得出一个概念。在这种情况下，有16个。因此，简单地输入一个不能用浮点数精确表示的数字，我们招致了大约\(10^{-16}\)的相对误差。

事实上，在所有情况下，我们都可以预计会招致大约\(10^{-16}\)的相对误差。这个量级称为计算机ε，通常写为\(\epsilon_{\text{MACHINE}}\)。它来自于两个相继浮点数之间的相对差异。与数字行不同的是，对于每个可表示的浮点数\(x\)，都有一个下一个浮点数，它大约是\(x+\epsilon_{\text{Machine}}x\)。因此，对于任意实数\(x_0\)，它落在两个浮点值\(x\)和\(x+\εx\)之间(为简明起见，去掉下标)。当我们以浮点形式表示\(x_0\)时，我们将获得这两个值中的一个。用\(\text{fl}(X_0)\)表示\(x_0\)的浮点表示。绝对误差是。

\BEGIN{公式*}e_{\text{abs}}=|\text{fl}(X_0)-x_0|\leq\max(|x_0-x|，|x_0-(x+\epsilon x)|)\leq|\epsilon x|。\end{公式*}。

凉爽的!。所以我们已经看到，在表示浮点数时，我们能做的最差的相对值大约是\(10^{-16}\)。这在几乎所有的实际目的中都是非常好的。因为记住，我们说的是相对误差。这意味着我们能够非常准确地表示非常小的值。以下是最接近\(10^{-20}\)的浮点数：

如果你好奇，我邀请你数一下9，一共有16个。即使在处理极小的数字时，我们也保持相同的相对精度。

现在你可能在想，等等！当我们表示浮点数时，我们可以获得很好的相对精度，这是很好的，但是当我们要对它们做些什么的时候呢？这里我们有两个浮点数，它们都不精确，我们要把它们相乘！谁知道会发生什么事呢？

这种担心是有根据的，因为浮点运算的算法必须以有限精度实现。如果我们被要求用笔和纸把两个大数相乘，我们大多数人将使用的算法是我们在学校学到的算法，它涉及到大量的加法、进位、减法等运算。如果所有这些中间步骤都使用某种二进制表示法，那么在我们进行的过程中，中间产品可能会失去精度！这看起来像是灾难的秘诀。幸运的是，浮点实现有一个我们可以要求的特性，IEEE-754和大多数其他流行的浮点标准都满足这个特性，这将把我们从完全的无政府状态中拯救出来。这就是Trefeten和Bau，数值线性代数，称为浮点算术的基本公理：

这意味着对于任意两个浮点数，例如\(x\)和\(y\)，涉及它们的任何操作都将给出一个浮点结果，该结果在真实结果的因子\(1+\epsilon_{\text{machine}}\)内。如果我们使用特殊的符号来表示浮点版本，比如说，\(+\)，这是最容易理解的。让我们写\(\Oplus\)表示浮点加法，\(+\)表示精确加法。那么基本公理告诉我们，

请记住，\(x\)和\(y\)完全可以表示为浮点数，但当然它们在数学上也是实数。因此，(x+y)是一个实数(数字线上的点)，它可能不是精确表示的推理点。基本公理告诉我们，假设我们的计算机可以访问数学的无限精度对象\(x+y\)，浮点运算\(\Oplus\)不会比简单地尝试将数字\(x+y\)表示为浮点值\(x+y\)做得更差。

换句话说，我们通过浮点运算不会产生额外的错误，而不是使用神奇的计算机对我们的浮点值进行精确的运算，然后再转换回浮点。用数学术语来说，

\BEGIN{公式*}\FRAC{\LEFT|(\text{fl}(A)\Oplus\text{fl}(B))-(a+b)\right|}{a+b}\Approx\frac{\Left|\text{fl}(\text{fl}(A)+\text{fl}(B))-(a+b)\right|}{a+b}\Approx\epsilon_{\text{机器}}。\。

回想一下，当我们试图用浮点表示\(A)时，\(\text{fl}(A)\)是我们实际表示的数学数字。因此，第一个分数给出了对浮点表示\(a\)和\(b\)执行\(Oplus)运算的相对误差，第二个分数给出了对浮点表示\(a\)和\(b\)执行精确算术运算\(+\)，然后将结果转换为浮点运算的相对误差。

浮点算术的基本公理使我们能够分析即使复杂的算术运算可能引起的数值误差。为了了解这是如何工作的，让我们考虑一下计算2D向量\([x，y]\)的长度的问题。从数学上讲，它的形式是。

计算有几个阶段，在每个阶段我们都会产生一点数值误差：

在每一步，我们得到的结果将等于真实结果，乘以某个误差因子\((1+\epsilon)\)，其中每个\(\epsilon\)非常小，在\(\epsilon_\text{机器}\)的另一端。每个操作可能有不同的错误\(\epsilon\)，但我们将对所有操作使用相同的符号。我们不太关心\(\\ε\)的精确值，只关心它非常小。

我们将把这些\(\ε\)通过计算，看看它们是如何影响最终结果的。为了方便起见，我们可以使用特殊的算术规则来操作\(\epsilon\)：

\(\Epsilon^2=0\)。如果\(epsilon\约10^{-16}\)，则\(\epsilon^2\约10^{-32}\)如此之小以至于我们决定完全忽略它。

\(\sqrt{1+\epsilon}=1+\frac{\epsilon}{2}-\frac{\epsilon^2}{8}+\cdots=1+\frac{\epsilon}{2}\)。这是通过在点1附近对sqrt进行泰勒展开而得到的。另一种观点是，当\(x\)非常接近1时，\(\sqrt{x}\近似x\)，且\(\sqrt{x}\)在\(x=1\)处的图的斜率是\(1/2)。

有了这些规则，我们就可以进行误差计算了。红色表达式是每个计算的新误差因子，在下一个表达式中显示为蓝色：

\BEGIN{ALIGN*}\text{fl}(X)&；={\color{red}{(1+\epsilon)}}x\text{fl}(Y)&；={\color{red}{(1+\epsilon)}}y\text{fl}(X)\oTimes\text{fl}(X)&；=(1+\epsilon)\Left[{\color{Blue}{(1+\epsilon)}}^2 x^2\右]=(1+\epsilon)^3 x^2\\&；=(1+3\epsilon+3\epsilon^2+\epsilon^3)x^2\\&；={\color{red}{(1+3\epsilon)}}x^2\text{fsilon。=(1+\epsilon)\Left[{\color{Blue}{(1+\epsilon)}}^2 y^2\右]=(1+\epsilon)^3 y^2\\&；=(1+3\epsilon+3\epsilon^2+\epsilon^3)y^2\&；={\color{red}{(1+3\epsilon)}}y^2\[\text{fl}(X)\oTimes\text{fl}(X)]\oplus[\text{fl}(Y)\oTimes\text{fl}(Y)]&；=(1+\epsilon)\Left[{\color{Blue}{(1+3\epsilon)}}x^2+{\color{Blue}{(1+3\epsilon)}}y^2\right]\\&；=(1+\epsilon+3\epsilon+3\epsilon^2)[x^2+y^2]\&；={\color{red}{(1+4\epsilon)}}(x^2+y^2)\sqrt[\text{fl}]{\text{fl}(X)\oTimes\text{fl}(X)\oplus\text{fl}(Y)\oTimes\text{fl}(Y)}&；=(1+\epsilon)\sqrt{{\color{Blue}{(1+。=(1+\epsilon)\sqrt{(1+4\epsilon)}\sqrt{x^2+y^2}\&；=(1+\epsilon)(1+2\epsilon)\sqrt{x^2+y^2}\&；={\color{red}{(1+3\epsilon)}\sqrt{x^2+y^2}\end{align*}。

哟！将我们的误差因素带入计算的所有步骤后，我们发现总的数值误差至多是机器精度的3倍。请记住，这是上限。我们真的不能说实际误差是多少，只能说它不大于\(3\epsilon_{\text{机器}}\)。

我们上面所做的误差分析非常随意地操纵了神奇的误差值\(\epsilon\)。我们不必太小心，因为使用这些规则，正值(如\(x^2)和\(y^2\))的乘法和加法是安全的。然而，当我们考虑减法时，我们必须更加小心。假设\(X\)和\(Y\)是非常大的数，但是它们非常接近：\(X-Y\)是小的。让我们试着用浮点数来计算\(X-Y\)，看看结果如何。现在我们可能会想写，

\BEGIN{ALIGN*}\text{fl}(X)\ominus\text{fl}(Y)&；=(1+\epsilon)[(1+\epsilon)(X)-(1+\epsilon)(Y)]\&；=(1+\epsilon)^2(X-Y)\&；=(1+2\epsilon)(X-Y)。\end{Align*}。

然而，我们不能从减法中提取\(1+\ε\)的因子！我们必须记住，我们只是在假装\(\epsilon\)在这些操作中是相同的。事实上，更准确地说，我们应该写。

\BEGIN{ALIGN*}\text{fl}(X)\ominus\text{fl}(Y)&；=(1+\epsilon_\ominus)[(1+\epsilon_X)(X)-(1+\epsilon_Y)(Y)]\end{Align*}

为了区分我们在不同步骤得到的误差值。现在，考虑如果\(\epsilon_Y=0\)，但\(\epsilon_X\NEQ 0\)会发生什么。嗯，那么，

\BEGIN{ALIGN*}(1+\epsilon_X)(X)-(1+\epsilon_Y)(Y)=X+\epsilon_X-Y=(X-Y)+\epsilon_X X.\end{Align*}。

上面我们假设(X-Y)很小，但是(X)很大。现在我们已经证明，减去\(X\)和\(Y\)的浮点表示所产生的误差是\(\epsilon_\text{Machine}X\)量级的，这可能比\(\epsilon_\text{Machine}\)大得多！下面是一段Python代码片段来说明：

我们应该得到0.125，但是最终得到了0.1，误差为25%！这种效应被称为“灾难性抵消”，数值分析人员非常谨慎地设计他们的算法来避免它。

浮点算术的基本公理为我们分析数值方法提供了一个起点，但这绝不是全部。最好的数值算法有一种称为“向后稳定”的神奇特性，它本质上确保它们准确地解决几乎是正确的问题--它们所解决的问题完全在你想要解决的问题的机器精度之内。这些算法的精确度远远超过它们应有的精确度。一旦你开始使用数值线性代数的算法来求解随时间演变的微分方程，就有了一套完整的科学来确保每个时间步产生的舍入误差不会以指数级膨胀。有了这些工具，我们就可以预测和控制它，而不是被数字误差吓倒。虽然看到你的计算给出几乎完全正确的答案是一种很好的感觉，但计算你得到的误差，并验证它也和预期的一样大，还有一些更令人满意的事情。

我所知道的关于这些东西的一切，我都是从劳埃德·特雷费森和大卫·鲍尔写的“数值线性代数”这本优秀的书中学到的，还有华盛顿大学无与伦比的安妮·格林鲍姆(Anne Greenbaum)。