Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Numeric function for log of sum in Python

Given log(a) and log(b), I want to compute log(a+b) (in a numerically stable way).

I wrote a little function for this:

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

I wrote a program where this is by far the most time-consuming piece of code. Obviously I could try to optimize it (eliminate the recursive call, for instance).

Do you know of a standard math or numpy function for computing log(a+b) from log(a) and log(b)?

If not, do you know of a simple way to make a single C++ hook for this function? It's not a complicated function (it uses floats), and as I said, it's taking up the majority of my runtime.

Thanks in advance, numerical methods ninja!

like image 292
user Avatar asked Sep 20 '11 06:09

user


People also ask

Is there log function in Python?

The math. log() method returns the natural logarithm of a number, or the logarithm of number to base.

How do you sum log values?

The rule is that you keep the base and add the exponents. Well, remember that logarithms are exponents, and when you multiply, you're going to add the logarithms. The log of a product is the sum of the logs.

How do you do log base 2 in Python?

log2(a) : This function is used to compute the logarithm base 2 of a. Displays more accurate result than log(a,2). Syntax : math. log2(a) Parameters : a : The numeric value Return Value : Returns logarithm base 2 of a Exceptions : Raises ValueError if a negative no. is passed as argument.


1 Answers

Note: Best answer until now is to simply use numpy.logaddexp(logA,logB).

Why exactly do you compare with log(0)? This is equal to -numpy.inf, in this case you come to log(1 + math.exp(-inf-logB) ) + logB Which reduces itself to logB. This call always will give an warning message which is extremely slow.

I could come up with this one-liner. However you'll need to really measure to see if this is actually faster. It does only use one 'complex' calculation function instead of the two that you use, and no recursion is happening, the if is still there but hidden (and maybe optimized) in fabs/maximum.

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

edit:

I did a quick timeit() with following results :

  1. Your original version took about 120s
  2. My version took about 30s
  3. I removed the compare with log(0) from your version and it came down to 20s
  4. I edited my code to keep the logaddexp but also worked with your recursive if and it went down to 18s.

Updated code, you could also switch the recursive call with an inline updated formula but this made little difference in my timing tests:

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

Edit 2:

As pv noted in comments, you could actually just do numpy.logaddexp(logA, logB) which comes down to calculating log(exp(logA)+exp(logB)) which is of course equal to log(A+B). I timed it (on the same machine as above) and it went further down to about 10s. So we've come down to about 1/12, not bad ;).

like image 61
KillianDS Avatar answered Nov 06 '22 01:11

KillianDS