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: stackoverflow.com/questions/1835787/…. 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
– slothrop, Commented May 1, 2023 at 9: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... — Matthew Towers
– Matthew Towers, Commented May 1, 2023 at 9:48
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
– slothrop, Commented May 1, 2023 at 14:07
And, just for the record, this is not a Python problem. These are the inherent limits of IEEE-754 floating point numbers. — Tim Roberts
– Tim Roberts, Commented May 2, 2023 at 1:44

hbwales · Accepted Answer · 2023-06-04 10:54:08Z

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]

Collectives™ on Stack Overflow

Python instability due to exponentials

1 Answer 1

Comments

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Linked

Related