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+b)log(a)log(b)

如果没有,您是否知道为此函数创建单个C++挂钩的简单方法？它不是一个复杂的函数(它使用浮点数),正如我所说,它占据了我运行时的大部分时间.

在此先感谢,数值方法忍者!

Answer 1

注意:到目前为止,最好的答案就是简单地使用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和它下降到18s.

更新了代码,您也可以使用内联更新公式切换递归调用,但这对我的时序测试没有什么影响:

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,不差;).