添加两个对数时防止下溢

Question

我正在使用维基百科 log probability 文章中描述的 log space 方程中的加法，但是在计算非常大的负对数的 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) 计算，并且会因除法而遭受精度损失。对吗？

Answer 1

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

For IEEE-compatible type double, overflow is guaranteed if 709.8 < arg, and underflow is guaranteed if arg < -708.4

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

If a range error occurs due to underflow, the correct result (after rounding) is returned.

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

但是，请注意

1 + exp(n)

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

所以由于 exp 造成的精度损失远小于添加 1.0 + exp(...)

时损失的精度

Answer 2

这里的问题是在没有中间under/overflow的情况下准确计算表达式log(1+exp(x))。幸运的是，Martin Maechler（R 核心开发人员之一）在 section 3 of this vignette.

中详细介绍了如何做到这一点

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

base 的选择不太可能对准确性（或性能）产生任何影响，并且大多数合理的数学库能够为两个函数提供低于 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 and log2），因此使用其中一个并没有真正的好处。

添加两个对数时防止下溢

Prevent underflow when adding two logarithms

c++

floating-point

precision

numerical-methods