如何在没有下溢的情况下计算 log-space 中的总和？

Question

我正在尝试计算给定 log(a) 和 log(b) 的 log(a + b)。问题是，log(a) 和 log(b) 太负了，以至于当我尝试自己计算 a 和 b 时，它们下溢，我得到 log(0)，这是未定义。

对于 log(a * b) 和 log(a / b)，这不是问题，因为 log(a * b) = log(a) + log(b) 和 log(a / b) = log(a) - log(b)。是否有类似的方法来计算 log(a + b) 而不需要 a 和 b 本身，避免下溢？

Answer 1

简而言之，使用以下表达式：

fmax(la, lb) + log1p(exp(-fabs(lb - la)))

这里我使用la和lb作为变量，分别存储log(a)和log(b)，函数名来自C的math.h。大多数其他语言也将具有这些功能，但它们的名称可能不同，例如abs 和 max 而不是 fabs 和 fmax。 （注意：这些是我将在整个回答中使用的约定。）

说到函数，: you might want to check whether you have access to one that will do this for you directly. For example, if you are using Python and NumPy, you have access to logaddexp, and SciPy has logsumexp。

如果您想了解以上内容的来源、如何将两个以上的数字相加以及如何相减的更多详细信息，请继续阅读。

更详细

没有像乘法和除法那样简单的规则，但有一个数学恒等式可以提供帮助：
log(a + b) = log(a) + log(1 + b/a)

我们可以稍微玩一下这个身份以获得以下内容：

log(a + b) = log(a) + log(1 + b/a)
           = log(a) + log(1 + exp(log(b) - log(a)))
           = la + log(1 + exp(lb - la))

这里还有问题。如果 la 比 lb 大得多，您将在 log 中得到 1 + 0.000...000something。浮点数尾数中没有足够的数字来存储 something，所以你只会得到 log(1)，完全失去 lb。

幸运的是，大多数编程语言在它们的标准库中都有一个函数来解决这个问题，log1p，它计算 1 加上它的自变量的对数。也就是说，log1p(x) returns log(1 + x) 但在某种程度上对于非常小的 x.

是准确的

所以，现在我们有：

log(a + b) = la + log1p(exp(lb - la))

我们快到了。还有一件事要考虑。通常，您希望 la 大于 lb。它并不总是很重要，但有时这会让你获得额外的精度。*如果 lb 和 la 之间的差异真的很大，这将使你免于 exp(lb - la) 中的溢出。在最极端的情况下，当 lb 为负无穷大（即 b 为 0）时计算有效，但当 la 为

时无效。

有时候，你会知道哪个更大，你可以把那个当作la。但是当它可能是其中之一时，您可以使用最大值和绝对值来解决它：

log(a + b) = fmax(la, lb) + log1p(exp(-fabs(lb - la)))

一个集合的总和

如果需要取两个以上数的和，我们可以推导出上面恒等式的扩展版本：

log(a[0] + a[1] + a[2] + ...)
    = log(a[0] * (a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...))
    = log(a[0]) + log(a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...)
    = la[0] + log(1 + exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

我们将要使用与将两个数字相加时类似的技巧。这样，我们就能得到最准确的答案，并尽可能避免溢出和下溢。首先，log1p:

log(a[0] + a[1] + a[2] + ...)
    = la[0] + log1p(exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

另一个考虑因素是在log1p 前面拉出哪个操作数。到目前为止，我使用 la[0] 进行演示，但您想使用最大的那个。这与我们在添加两个数字时想要 la > lb 的所有相同原因。例如，如果 la[1] 最大，您可以这样计算：

log(a[0] + a[1] + a[2] + ...)
    = la[1] + log1p(exp(la[0]-la[1]) + exp(la[2]-la[1]) + ...)

将其放入适当的代码中，它看起来像这样（这是 C，但它应该很好地翻译成其他语言）：

double log_sum(double la[], int num_elements)
{
    // Assume index_of_max() finds the maximum element
    // in the array and returns its index
    int idx = index_of_max(la, num_elements);

    double sum_exp = 0;
    for (int i = 0; i < num_elements; i++) {
        if (i == idx) {
            continue;
        }
        sum_exp += exp(la[i] - la[idx]);
    }

    return la[idx] + log1p(sum_exp);
}

正在计算对数差异 space

这不是问题的一部分，但由于它可能仍然有用：log space 中的减法可以类似地执行。基本公式是这样的：

log(a - b) = la + log(1 - exp(lb - la))

请注意，这仍然假设 la 大于 lb，但对于减法，它更重要。如果la小于lb，你就是取负数的对数！

与加法类似，这有一个准确性问题，可以通过使用专门的函数来解决，但事实证明有两种方法。一个使用与上面相同的 log1p 函数，但另一个使用 expm1，其中 expm1(x) returns exp(x) - 1。以下是两种方式：

log(a - b) = la + log1p(-exp(lb - la))
log(a - b) = la + log(-expm1(lb - la))

您应该使用哪一个取决于 -(lb - la) 的值。当 -(lb - la) 大于大约 0.693（即 log(2)）时，第一个更准确，而第二个则更准确。有关为什么会这样以及 log(2) 从何而来的更多详细信息，请参阅此 note from the R project 评估这两种方法。

最终结果如下所示：

(lb - la < -0.693) ? la + log1p(-exp(lb - la)) : la + log(-expm1(lb - la))

或函数形式：

double log_diff(double la, double lb)
{
    if (lb - la < -0.693) {
        return la + log1p(-exp(lb - la));
    } else {
        return la + log(-expm1(lb - la));
    }
}

* 这有一点甜蜜的地方。当 la 和 lb 之间的差异很小时，无论哪种方式，答案都是准确的。当差异太大时，结果将始终等于两者中的较大者，因为浮点数没有足够的精度。但是当差异恰到好处时，la 越大，您将获得更好的准确性。

Answer 2

不失一般性地假设 b

log(a+b) = log(a) + log(1 + b/a)
         = log(a) + b/a - 1/2(b/a)^2 + 1/3(b/a)^3 etc.

涉及b/a的项可以根据已知量计算，因为

b/a = exp(log(b) - log(a))

如果 b/a 的计算下溢，则 log(a+b) = log(a) 到机器精度。

如何在没有下溢的情况下计算 log-space 中的总和？

How to calculate sums in log-space without underflow?

floating-point

language-agnostic

precision

logarithm

underflow

更详细

一个集合的总和

正在计算对数差异 space