Overflow in numpy.exp()

Question

I have to calculate the exponential of the following array for my project:

w  = [-1.52820754859, -0.000234000845064, -0.00527938881237, 5797.19232191, -6.64682108484,
       18924.7087966, -69.308158911, 1.1158892974, 1.04454511882, 116.795573742]

But I've been getting overflow due to the number 18924.7087966.

The goal is to avoid using extra packages such as bigfloat (except "numpy") and get a close result (which has a small relative error).

1.So far I've tried using higher precision (i.e. float128):

def getlogZ_robust(w):

    Z = sum(np.exp(np.dot(x,w).astype(np.float128)) for x in iter_all_observations())
    return np.log(Z)

But I still get "inf" which is what I want to avoid.

I've tried clipping it using nump.clip():

def getlogZ_robust(w):

    Z = sum(np.exp(np.clip(np.dot(x,w).astype(np.float128),-11000, 11000)) for x in iter_all_observations())
    return np.log(Z)

But the relative error is too big.

Can you help me solving this problem, if it is possible?

Can you use SciPy? If so, check out logsumexp. If not, you can grab the source code from github.com/scipy/scipy/blob/master/scipy/special/_logsumexp.py, or search a bit here on stackoverflow. There are probably many implementations of the same method scattered about stackoverflow answers. — Warren Weckesser
– Warren Weckesser, Commented Nov 10, 2019 at 16:52
I can't use SciPy, but I've checked out logsumexp and found out that numpy also has a similar function called 'logaddexp' which computes logaddexp(x1, x2) == log(exp(x1) + exp(x2)) without explicitly computing the intermediate exp() values. This way it avoids the overflow. So thank you after all. — erenerogullari
– erenerogullari, Commented Nov 10, 2019 at 17:59

Gianluca · Accepted Answer · 2019-11-10 18:16:42Z

Only significantly extended or arbitrary precision packages will be able to handle the huge differences in numbers. The exponential of the largest and most negative numbers in w differ by 8000 (!) orders of magnitude. float (i.e. double precision) has 'only' 15 digits of precision (meaning 1+1e-16 is numerically equal to 1), such that adding the small numbers to the huge exponential of the largest number has no effect. As a matter of fact, exp(18924.7087966) is so huge, that it dominates the sum. Below is a script performing the sum with extended precision in mpmath: the ratio of the sum of exponentials and exp(18924.7087966) is basically 1.

w  = [-1.52820754859, -0.000234000845064, -0.00527938881237, 5797.19232191, -6.64682108484,
       18924.7087966, -69.308158911, 1.1158892974, 1.04454511882, 116.795573742]

u = min(w)
v = max(w)

import mpmath
#using plenty of precision
mpmath.mp.dps = 32768
print('%.5e' % mpmath.log10(mpmath.exp(v)/mpmath.exp(u)))
#exp(w) differs by 8000 orders of magnitude for largest and smallest number

s = sum([mpmath.exp(mpmath.mpf(x)) for x in w])

print('%.5e' % (mpmath.exp(v)/s))
#largest exp(w) dominates such that ratio over the sums of exp(w) and exp(max(w)) is approx. 1

Gianluca · Accepted Answer · 2019-11-10 18:33:50Z

4

If the issues of loosing digits in the final results due to hugely different orders of magnitudes of added terms in not a concern, one could also mathematically transform the log of sums over exponentials the following way avoiding exp of large numbers:

log(sum(exp(w)))
= log(sum(exp(w-wmax)*exp(wmax)))
= wmax + log(sum(exp(w-wmax)))

In python:

import numpy as np
v = np.array(w)
m = np.max(v)
print(m + np.log(np.sum(np.exp(v-m))))

Note that np.log(np.sum(np.exp(v-m))) is numerically zero as the exponential of the largest number completely dominates the sum here.

answered Nov 10, 2019 at 18:33

Gianluca

5464 silver badges7 bronze badges

Comments

erenerogullari · Accepted Answer · 2019-11-10 18:05:09Z

1

Numpy has a function called logaddexp which computes

logaddexp(x1, x2) == log(exp(x1) + exp(x2))

without explicitly computing the intermediate exp() values. This way it avoids the overflow. So here is the solution:

def getlogZ_robust(w):

    Z = 0
    for x in iter_all_observations():
        Z = np.logaddexp(Z, np.dot(x,w))
    return Z

answered Nov 10, 2019 at 18:05

erenerogullari

272 silver badges8 bronze badges

Collectives™ on Stack Overflow

Overflow in numpy.exp()

3 Answers 3

Comments

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

Comments

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related