我想计算 log(a + b)，已知log(a)和log(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本身，并避免下溢？
简而言之，使用以下表达式：
fmax(la, lb) + log1p(exp(-fabs(lb - la)))
我使用la和lb作为变量来存储log(a)和log(b)，函数名称来自于C语言的math.h库。大多数其他编程语言也会有这些函数，但它们可能命名不同，例如用abs和max代替fabs和fmax。（注意：这是我在本回答中将使用的约定）
谈到函数，Mark Dickinson提出了一个很好的观点：您可能希望检查您是否可以直接访问可执行此操作的函数。例如，如果您正在使用Python和NumPy，则可以访问logaddexp，而SciPy则有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，它计算其参数的对数。也就是说，log1p(x)返回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代码，但其他语言也应该能很好地转换）：
把这个写成正确的代码，它看起来会像这样（这是 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);
}
尽管这并不是问题的一部分，但由于可能仍然有用：对数空间中的减法可以类似地执行。基本公式如下：
log(a - b) = la + log(1 - exp(lb - la))
注意，这仍假设la大于lb，但对于减法更为重要。如果la小于lb，则会取负数的对数！
与加法类似，这也存在精度问题，可以通过使用特定的函数来解决，但有两种方法。一种方法使用与上文相同的log1p函数，另一种使用expm1，其中expm1（x）返回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)的来历，请参阅 R 项目中这份说明文档比较这两种方法。
最终结果如下：
(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更大时，您将获得更好的精度。
不失一般性，假设 b<a。
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 = exp(log(b) - log(a))