防止在添加两个对数时下溢

Prevent underflow when adding two logarithms

本文关键字：两个 下溢 添加      更新时间：2023-10-16

我正在使用维基百科对数概率文章中描述的对数空间方程中的加法，但是在计算非常大的负对数的exp时，我遇到了下溢。结果，我的程序崩溃了。

示例输入为a = -2和b = -1033.4391885529124。

我的代码直接从维基百科文章中实现，如下所示：

double log_sum(double a, double b)
{
double min_ab = std::min(a, b);
a = std::max(a, b);
b = min_ab;
if (isinf(a) && isinf(b)) {
return -std::numeric_limits<double>::infinity();
} else if (isinf(a)) {
return b;
} else if (isinf(b)) {
return a;
} else {
return a + log2(1 + exp2(b - a));
}
}

我想出了以下想法，但无法决定哪个是最好的：

  • 评估前检查超出范围的输入。
  • 禁用(以某种方式)异常，并在评估后冲洗或箝位输出
  • 实现不会引发异常并自动刷新或固定结果的自定义日志和 exp 函数。
  • 还有其他方式吗？

此外，我很想知道选择对数基数对计算有什么影响。我之所以选择二基数，是因为我相信其他对数基数会根据log_n(x) = log_2(x) / log_2(n)计算，并且会因除法而遭受精度损失。这是对的吗？

根据http://en.cppreference.com/w/cpp/numeric/math/exp：

对于与IEEE兼容的双精度型，如果参数<709.8，则保证溢出，如果参数<-708.4，则保证下溢

所以你无法阻止下溢。然而：

如果由于下溢而发生范围错误，则返回正确的结果(舍入后)。

所以不应该有任何程序崩溃 - "只是"精度的损失。

但是，请注意

1 + exp(n)

会更快地失去精度，即已经在 n = -53 处。这是因为1.0之后的下一个可表示数字是1.0 + 2^-52。

因此，由于exp造成的精度损失远远小于添加1.0 + exp(...)时损失的精度

这里的问题是在没有中间下溢/溢出的情况下准确计算表达式log(1+exp(x))。幸运的是，Martin Maechler(R核心开发人员之一)在本小插曲的第3节中详细介绍了如何做到这一点。

他使用自然基函数：应该可以通过适当缩放函数将其转换为 base-2，但它在一个部分中使用log1p函数，而且我不知道有任何数学库提供 base-2 变体。

基数的选择不太可能对准确性(或性能)产生任何影响，大多数合理的数学库能够为这两个函数提供sub 1-ulp保证(即您将拥有最接近确切答案的两个浮点值之一)。一种非常常见的方法是将浮点数分解为其以 2 为底的指数k和有效数1+f，这样1/sqrt(2) < 1+f < sqrt(2)，然后使用多项式近似来计算log(1+f)：由于一些数学怪癖(基本上，泰勒级数的第 2 项可以精确表示的事实)，事实证明，在自然底而不是以 2 为底的情况下这样做更准确， 因此，典型的实现如下所示：

log(x) = k*log2 + p(f)
log2(x) = k + p(f)*invlog2

(例如，参见 openlibm 中的 log 和 log2)，因此使用一个没有真正的好处。

相关文章：