避免溢出、下溢和精度损失

• 下载源代码 - 7.38 KB

引言

计算数学函数的最直接方法可能不起作用。例如，如何计算 log(x+1)/x？最直接的方法是将 x 加 1，取对数，然后除以 x。对于某些 x 值，这将产生一个准确的结果。而对于其他值，结果可能会偏离 100%。

本文提供了三个函数直接计算可能失败的示例。在此过程中，它给出了精确计算函数的四条通用规则。

本文附带的代码将文章中的原理应用于九个统计学中经常出现的函数的计算。

示例 1：二次方程

假设您想找到 3x2 + 109x + 5 = 0 的根。这很简单：只需应用二次方程。有什么可能出错？让我们看看。

double a = 3.0, b = 1.0e9, c = 5.0;
double d = b*b - 4.0*a*c;
double r1 = (-b - sqrt(d))/(2.0*a);
double r2 = (-b + sqrt(d))/(2.0*a);

这段代码正确地计算了 r1 到所有显示位数，但将 r2 计算为 0，即使其精确到机器精度的正确值为 -4.5e-8。计算 r2 的绝对误差很小，但相对误差为 100%。为什么一个根计算准确而另一个不准确？

数值分析中的一个基本规则是避免减去接近的数。两个数越接近，减法损失的精度就越多。由于 b 远大于 a 和 c，因此变量 sqrt(d) 与 b 非常接近。变量 b 和 sqrt(d) 在小数点后 12 位相同，因此减法损失了大约 12 位精度。我们如何更准确地计算 r2？

规则 1：使用代数技巧避免精度损失。

二次方程有一个等价但不太为人知的形式。将标准形式的分子和分母乘以 ± 符号相反的分子，然后进行简化。这个新形式在分母中将 b 和 sqrt(d) 相加，避免了减去接近的数的问题。修改后的代码如下，完全准确。

double a = 3.0, b = 1.0e9, c = 5.0;
double d = b*b - 4.0*a*c;
double r1 = (-b - sqrt(d))/(2.0*a);
double r2 = -2.0*c/(b + sqrt(d));

仍然存在一些不准确之处，因为在机器精度下 d 等于 b*b，但我们将错误率从 100% 提高到了 10%。为了获得更高的精度，请放弃二次方程，而使用求根算法。

同样，请参阅 比较计算标准差的三种方法，其中包含数值示例，说明了三个代数上等价的公式如何表现出截然不同的行为。在 拟合一条直线到数据 时也会发生类似的现象。

示例 2：阶乘之比

概率计算通常涉及将非常大的数相除以得到一个中等大小的数。最终结果可能在 double 的范围内绰绰有余，但中间结果会发生溢出。例如，假设您需要计算从 200 个元素集合中选择 10 个元素的方法数。这等于 200!/(190! 10!)，约等于 2.2 e16。但是 200! 和 190! 会溢出 double 的范围。

有两种方法可以解决这个问题。两者都遵循以下规则。

规则 2：使用代数技巧避免溢出。

第一个技巧是认识到 200! = 200*199*198*...*191*190!，因此 200!/(190! 10!) = 200*199*198*...*191/10!。这肯定有效，但仅限于阶乘。更通用的技术是使用对数来避免溢出：取要计算的表达式的对数，然后对结果进行指数运算。在此示例中，log( 200!/(190! 10!) ) = log(200!) - log(190!) - log(10!)。如果您有直接计算阶乘对数而不先计算阶乘的代码，您可以使用它来找到所需结果的对数，然后应用 exp 函数。有关此示例的更多信息，请参阅 浮点计算的五个技巧。

示例 3：log(1 + x)

现在让我们来看介绍中提到的计算 log(x+1) 的示例。考虑以下代码

double x = 1e-16;
double y = log(1 + x)/x;

这段代码将 y 设置为 0，尽管正确值等于 1 到机器精度。出了什么问题？

双精度数精确到大约 15 位小数，因此在机器精度下 1 + x 等于 1。1 的对数是零，所以 y 被设置为零。但对于小的 x 值，log(1 + x) 近似于 x，因此 log(1+x)/x 近似于 1。这意味着上面计算 log(1 + x)/x 的代码返回的结果相对误差为 100%。如果 x 不小到 1 + x 在机器中等于 1，我们仍然可能遇到问题。如果 x 中等大小，那么在计算 1 + x 时 x 的位没有完全丢失，但其中一些丢失了。x 越接近 0，丢失的位越多。

这个问题与二次方程的问题类似，所以我们可能会尝试一个类似的解决方案。不幸的是，代数技巧无济于事。所以我们尝试我们的第二个规则。

规则 3：使用解析近似值避免精度损失。

数值分析师最喜欢的近似形式是幂级数。log(1 + x) 的幂级数是 x + x2/2! + x3/3! + ... 对于小的 x 值，简单地返回 x 作为 log(1 + x) 是一个改进。这对于最小的 x 值效果很好，但对于一些不太小的值，这仍然不够准确，直接计算 log(1 + x) 也是如此。本文附带的示例代码对小的 x 值使用有理逼近，对较大的值使用直接计算。

示例 4：逆 Logit 函数

接下来我们来看计算 f(x) = ex/(1 + ex)。(统计学家称之为“逆 Logit”函数，因为它与他们称之为“Logit”函数的反函数。) 最直接的方法是计算 exp(x)/(1 + exp(x))。让我们看看它可能在哪里出错。

double x = 1000;
double t = exp(x);
double y = t/(1.0 + t);

打印 y 的结果是 -1.#IND。这是因为 t 的计算发生了溢出，产生了一个大于 double 范围的数。但我们可以轻松地算出 y 的值。如果 t 大到我们无法存储，那么 1 + t 就基本上等于 t，比率就非常接近 1。这表明我们需要找出 f(x) 在机器精度下等于 1 的 x 值，然后对于这些 x 值直接返回 1，以避免溢出的可能性。这是我们的第三条规则。

规则 4：不要计算您能准确预测的结果。

头文件 float.h 有一个常量 DBL_EPSILON，它是我们可以加到 1 而不会得到 1 的最小 double 精度。一些代数运算表明，如果 x 大于 -log(DBL_EPSILON)，那么 f(x) 在机器精度下将等于 1。所以这里有一段计算大 x 值的 f(x) 的代码片段。

const double x_max = -log(DBL_EPSILON);
if (x > x_max) return 1.0;

Using the Code

本文提供的代码计算了统计学中出现的七个函数。每个函数都避免了在参数为大负数、大正数或接近零时可能发生的溢出、下溢或精度损失问题。

• LogOnePlusX 计算 log(1 + x)，如示例 3 所示。
• ExpMinusOne 计算 ex-1。
• Logit 计算 log(x/(1-x))。
• LogitInverse 计算 ex/(1 + ex)，如示例 4 所述。
• LogLogitInverse 计算 log(ex/(1 + ex))。
• LogitInverseDifference 计算 LogitInverse(x) - LogitInverse(y)。
• LogOnePlusExpX 计算 log(1 + exp(x))。
• ComplementaryLogLog 计算 log( -log(1 - x) )。
• ComplementaryLogLogInverse 计算 1.0 - exp(-exp(x))。

关注点

这里提出的解决方案起初看起来不必要，甚至错误。如果这类代码没有得到良好的注释，有人会过来“简化”它，而且是错误的。他们会为删除了所有不必要的混乱而感到自豪。如果他们不测试极端值，他们的新代码似乎就能正常工作。错误的答案和 NaN 只会在之后出现。

历史

• 2008年4月16日：初稿
• 2008 年 4 月 18 日：修订文章以澄清几点。
• 2008 年 10 月 9 日：添加了进一步示例的链接。
• 2008 年 10 月 29 日：添加了示例。