How to add large exponential terms reliably without overflow errors?

Question

A very common problem in Markov Chain Monte Carlo involves computing probabilities that are sum of large exponential terms,

$ e^{a_1} + e^{a_2} + ... $

where the components of $a$ can range from very small to very large. My approach has been to factor out the largest exponential term $K := \max_{i}(a_{i})$ so that:

$$a' =K + log\left( e^{a_1 - K} + e^{a_2 - K } + ... \right)$$ $$e^{a'} \equiv e^{a_1} + e^{a_2} + ...$$

This approach is reasonable if all elements of $a$ are large, but not such a good idea if they aren't. Of course, the smaller elements aren't contributing to the floating-point sum anyway, but I'm not sure how to reliably deal with them. In R code, my approach looks like:

if ( max(abs(a)) > max(a) )
  K <-  min(a)
else
  K <- max(a)
ans <- log(sum(exp(a-K))) + K

It seems a common enough problem that there should be a standard solution, but I'm not sure what it is. Thanks for any suggestions.

Answer 1

There is a straightforward solution with only two passes through the data:

First compute $$K := \max_i\; a_i,$$

which tells you that, if there are $n$ terms, then $$\sum_i e^{a_i} \le n e^K.$$

Since you presumably don't have $n$ anywhere near as large as even $10^{20}$, you should have no worry about overflowing in the computation of $$\tau := \sum_i e^{a_i-K} \le n$$ in double precision.

Thus, compute $\tau$ and then your solution is $e^K \tau$.

Answer 2

To keep precision while you add doubles together you need to use Kahan Summation, this is the software equivalent to having a carry register.

This is fine for most values, but if you are getting overflow then you are hitting the limits of IEEE 754 double-precision which would be about $e^{709.783}$. At this point you need a new representation. You can detect an overflow at addition time by doubleMax - sumSoFar < valueToAdd and also detect exponents to large to evaluate by exponent > 709.783. At this point you can modify the interpretation of a double by shifting the exponent and keeping track of this shift.

This for the most part is the similar to your exponent offseting approach, but this version is kept it in base 2 and does not require an initial search to find the largest exponent. Hence $\mathit{value}\times2^\mathit{shift}$.

#!/usr/bin/env python
from math import exp, log, ceil

doubleMAX = (1.0 + (1.0 - (2 ** -52))) * (2 ** (2 ** 10 - 1))

def KahanSumExp(expvalues):
  expvalues.sort() # gives precision improvement in certain cases 
  shift = 0 
  esum = 0.0 
  carry = 0.0 
  for exponent in expvalues:
    if exponent - shift * log(2) > 709.783:
      n = ceil((exponent - shift * log(2) - 709.783)/log(2))
      shift += n
      carry /= 2*n
      esum /= 2*n
    elif exponent - shift * log(2) < -708.396:
      n = floor((exponent - shift * log(2) - -708.396)/log(2))
      shift += n
      carry *= 2*n
      esum *= 2*n
    exponent -= shift * log(2)
    value = exp(exponent) - carry 
    if doubleMAX - esum < value:
      shift += 1
      esum /= 2
      value /= 2
    tmp = esum + value 
    carry = (tmp - esum) - value 
    esum = tmp
  return esum, shift

values = [10, 37, 34, 0.1, 0.0004, 34, 37.1, 37.2, 36.9, 709, 710, 711]
value, shift = KahanSumExp(values)
print "{0} x 2^{1}".format(value, shift)

Answer 3

Your approach is solid.

You don't need to know $K$ exactly, just well enough to avoid overflow. So you may be able to estimate $K$ analytically before you do any MCMC sampling.

Answer 4

There is an R package which supplies a fast and efficient implementation of the "log-sum-exp trick"

http://www.inside-r.org/packages/cran/matrixStats/docs/logSumExp

The logSumExp function accepts a numeric vector lX and outputs log(sum(exp(lX))) while avoiding underflow and overflow problems using the method you described.