Search code examples
pythonnumpynumericalexponentialstability

How to order this computation for numerical stability?

I'm trying to compute a vector, whose sum is 1 and whose elements are defined as such:

v[i] = exp(tmp[i])/exp(tmp).sum()

The problem is that the value in the exponential may be large (between -10^2 and 10^2), making the exponential to evaluate to inf or 0.

I tried some variations, like substracting the biggest element or the mean of tmp to the numerator and the denominator, but it's still not enough.

Basically, I'm need of a transformation that reduces the mean value and the dispersion in tmp or of a clever ordering for this computation.

I'm using numpy arrays as containers and exp is numpy.exp.

Solution 
  • >>> tmp = np.array([-10**10, 10**10])
>>> tmp_max = tmp.max()
>>> log_D = log(sum(exp(tmp - tmp_max))) + tmp_max
>>> log_v = tmp - log_D
>>> v = np.exp(log_v)
>>> v
array([ 0.,  1.])

    Or use scipy.misc.logsumexp, which uses the exact same algorithm.