Python instability due to exponentials

Question

Setup: I have the following function in python, where x can get very large:

import numpy as np

def function(x, pi):
  d = len(pi)
  output = 0
  for r in range(d):
    output += pi[r] * np.exp(-x)
  return output

Input description: x can be very large causing np.exp(-x) to evaluate to zero which results in the entire function being zero, and pi is just a vector of probabilities (e.g., [0.5, 0.5]).

Question: Is there a more stable way to implement this function such that it wouldn't lead to the output being zero? Thanks.

Edit: I have decided to give more details since it was asked in the comments. The entire function is

def entire_function(x_array, pi, r):
  d = len(pi)
  numerator = np.exp(-x_array[r])
  denominator = 0
  for r_prime in range(d):
    denominator += pi[r_prime] * np.exp(-x_array[r_prime])
  return numerator / denominator

Even trying to use np.log doesn't really help. For example:

a = np.array([np.exp(-900), np.exp(-800)])
print(np.log(a[0]+a[1]))

This gives me -Inf. The summation in the denominator is the nasty part that is giving me trouble since it is preventing me from accessing the exponents (to make the computation more numerically stable). I guess this issue is similar to the logsumexp examples in machine learning with the extra pi[r] factors in front.

The function can be simplified to just `np.sum(pi) * np.exp(-x)`, but that doesn't help if x is so big that the result is too small for floating point to handle. The smallest storeable float is on the order of 1e-308: https://stackoverflow.com/questions/1835787/what-is-the-range-of-values-a-float-can-have-in-python. If you need to represent values smaller than that, you could think about storing the logarithm of the value rather than the value itself. — slothrop, May 01 '23 at 09:39
nothing can stop `exp(-x)` being too small for a Python float if `x` is large enough and positive, but you've got a few options: avoid having to evaluate `exp(-x)` by rearranging whatever overall math expression you are computing, transform the data you're working with somehow so that it lives in a different range, increase precision by using something like [mpmath](https://mpmath.org/)... — Matthew Towers, May 01 '23 at 09:48
x is around 700 to 900. That is around how large x needs to be. — Resu, May 01 '23 at 13:59
Right, so exp(-700) is about 1e-304 and exp(-900) is about 1e-391. You're likely dealing with numbers too small for floating point to represent, so it needs an alternative strategy. — slothrop, May 01 '23 at 14:07
I wonder if the result is going to be normalized afterwards (e.g. something like a softmax function). If so, consider handling the normalization in the same place as the exp. Divide through by the exp(maximum value of x) -- that will make the leading term 1, and everything else smaller. If the difference between x values is very large, some of the terms will underflow to 0.0 -- no problem there, since you're guaranteed at least one term is nonzero. — Robert Dodier, May 01 '23 at 19:40
Its the opposite, this function is the normalising constant. So i am going to use this function to divide one of the terms in the for loop. I.e., pi[r]*np.exp(-x) for a chosen r is the numerator and this function is the denominator. Hence, I couldn't use what you suggested unfortunately... — Resu, May 01 '23 at 23:55
**Note:** I have added the entire function to the question for greater clarity. It is under the final section titled "Edit:". — Resu, May 02 '23 at 01:05
And, just for the record, this is not a Python problem. These are the inherent limits of IEEE-754 floating point numbers. — Tim Roberts, May 02 '23 at 01:44
Apart from the *range* of exponents that floats can handle, *precision* would also be an issue here. For example, `exp(-900) + exp(-800)` is equal to `[exp(-800) * (exp(-100) + 1)]`, which to about 43 significant figures is just equal to `exp(-800)`. — slothrop, May 02 '23 at 08:16

score 1 · Accepted Answer · answered Jun 04 '23 at 10:54

Note that in general we have:

pe^x = e^log(p) e^x = e^{log(p) + x}

Using this we can apply the log-sum-exp trick you linked

import numpy as np

xs = np.array([700, 900])
ps = np.array([0.6, 0.4])

def original(xs, ps, r):
    ex = np.exp(-xs) 
    return ex[r] / (ps*ex).sum()

def log_sum_exp(x):
    c = x.max()
    return c + np.log(np.sum(np.exp(x - c)))

def adjusted(xs, ps, r):
    return np.exp(-xs[r] - log_sum_exp(-xs + np.log(ps)))

Which we can check with

def check(xs, ps, r):
    # calculate to 1,000 decinal places to check result against
    from decimal import Decimal, getcontext
    getcontext().prec = 1000
    ex = [Decimal.exp(-Decimal(float(xi))) for xi in xs]
    return ex[r] / sum(Decimal(float(pi))*ei for pi,ei in zip(ps, ex))

print([adjusted(xs, ps, i) for i in range(2)])       # [1.6666666666666516, 2.306494211227875e-87]
print([float(check(xs, ps, i)) for i in range(2)])   # [1.6666666666666667, 2.306494211227896e-87]

Python instability due to exponentials

1 Answers1