Python中求和的对数的数值函数

Question

给定 log(a) 和 log(b)，我想计算 log(a+b)（要以一种数值稳定的方式进行）。

我写了一个小函数来实现这个：

def log_add(logA,logB):
    if logA == log(0):
        return logB
    if logA<logB:
        return log_add(logB,logA)
    return log( 1 + math.exp(logB-logA) ) + logA

我写了一个程序，其中这部分代码是 绝对 最耗时的。显然，我可以尝试优化它（比如消除递归调用）。

你知道有没有标准的 math 或 numpy 函数可以用来从 log(a) 和 log(b) 计算 log(a+b) 吗？

如果没有，你知道有没有简单的方法可以为这个函数做一个 C++ 的接口吗？这个函数并不复杂（它使用浮点数），正如我所说，它占用了我大部分的运行时间。

提前谢谢你，数值方法高手！

Answer 1

def log_add(logA, logB): 
    return math.log(math.exp(logA) + math.exp(logB))

是不是太慢了？或者不准确呢？

Answer 2

注意：到目前为止，最好的答案就是直接使用 numpy.logaddexp(logA, logB)。

为什么要和 log(0) 比较呢？这个值等于 -numpy.inf，在这种情况下，你会得到 log(1 + math.exp(-inf-logB)) + logB，这最后会简化成 logB。这个调用总是会发出警告信息，而且速度非常慢。

我想出了一个一行代码的写法。不过你需要真正测量一下，看看这个方法是否真的更快。它只用了一个“复杂”的计算函数，而不是你用的两个，而且没有递归，if 语句仍然存在，但在 fabs/maximum 中被隐藏（可能还进行了优化）。

def log_add(logA,logB):
    return numpy.logaddexp(0,-numpy.fabs(logB-logA)) + numpy.maximum(logA,logB)

编辑：

我用 timeit() 快速测试了一下，结果如下：

1. 你原来的版本大约花了 120 秒
2. 我的版本大约花了 30 秒
3. 我从你的版本中去掉了与 log(0) 的比较，时间降到了 20 秒
4. 我编辑了我的代码，保留了 logaddexp，同时也用了你的递归 if，时间降到了 18 秒。

更新后的代码，你也可以用一个内联更新的公式替换递归调用，但在我的时间测试中，这几乎没有什么区别：

def log_add2(logA, logB):
    if logA < logB:
        return log_add2(logB, logA)
    return numpy.logaddexp(0,logB-logA)+logA

编辑 2：

正如 pv 在评论中提到的，你其实可以直接使用 numpy.logaddexp(logA, logB)，这实际上是在计算 log(exp(logA) + exp(logB))，当然这等于 log(A + B)。我在同一台机器上测试，时间进一步降到了大约 10 秒。所以我们把时间缩短到了大约 1/12，效果不错 ;)。