1

I need to compute the np.sum(np.exp(L)) where L is a long list of numbers which can be as high as 3000. np.float64 will overflow around exp(700). However, np.logaddexp is immune to overflow but only works for two elements.

I wrote the following recursive function to compute it for a list.

import numpy as np
from smartprint import smartprint as sprint 

def logsum(l):
    if len(l) == 1:
        return l[0]
    return logsum([np.logaddexp(l[0], l[1])] + l[2:])


# odd length
l = [-10, 1, 2, 45, 100]

assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))

# even length
l = [-10, 1, 2, 43]

assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))

Output:

logsum(l) : 100.0
logsum(l) : 43.0

Additionally, logsum([1,2,3,4,50000]) outputs 50000 as expected whereas the np.exp would overflow. Now, the problem is that my list L is huge and has upto 10 million elements and I need to compute the logsum at least 1000 times. So, I wonder if there is a way to vectorize using the np.logaddexp so that I can use it for a longer list instead of just two elements.

lifezbeautiful
  • 1,160
  • 8
  • 13
  • See: https://stackoverflow.com/questions/59850770/fastest-python-log-sum-exp-in-a-reduceat – slothrop Jun 01 '23 at 07:51
  • 1
    Thanks @slothrop, yes the scipy version does the job for me, The question you referenced is aiming for something faster, which is definitely a niche requirement. Mine was more generic. The scipy solution works perfectly fine for my range of data values. and as expected is order of magnitudes faster than computing it in python explicitly. – lifezbeautiful Jun 01 '23 at 08:07

1 Answers1

1

Based on user: slothrop's comment, the scipy function for logsumexp works well and is already vectorised as expected. My recursive solution is no longer needed. Posting here for completeness.

from scipy.special import logsumexp
logsumexp(l)
lifezbeautiful
  • 1,160
  • 8
  • 13