25
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ This is a thing. Google for 'logsumexp'. $\endgroup$
    – user1026
    Feb 22, 2012 at 16:35

4 Answers 4

15
$\begingroup$

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$.

$\endgroup$
6
  • $\begingroup$ Thanks for the clear notation -- but I believe this is essentially what I have proposed(?) If I need to avoid underflow errors when some $a_i$ are small, I gather I need the Kahan summation approach proposed by @gareth? $\endgroup$
    – cboettig
    Feb 2, 2012 at 0:53
  • $\begingroup$ Ah, I now see what you were getting at. You actually don't need to worry about underflow, as adding exceptionally tiny results to your solution shouldn't change it. If there was an exceptionally large number of them, then you should sum the small values first. $\endgroup$ Feb 2, 2012 at 0:57
  • $\begingroup$ To the downvoter: would you mind letting me know what is wrong with my answer? $\endgroup$ Feb 13, 2012 at 0:06
  • $\begingroup$ what if you have many very small terms? It could happen that $e^{a_i - K} \approx 0$ for these. If there are many terms like this, you would have a large error. $\endgroup$
    – a06e
    Aug 4, 2016 at 17:43
  • $\begingroup$ scicomp.stackexchange.com/q/24624/988 $\endgroup$
    – a06e
    Aug 4, 2016 at 19:39
10
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ Kahan summation is but one of a family of "compensated summation" methods. If for some reason Kahan doesn't work quite right, there are a number of other methods for adding up terms of varying magnitudes and opposite signs properly. $\endgroup$
    – J. M.
    Feb 2, 2012 at 1:00
  • $\begingroup$ @J.M. could you provide me with the names of those other methods, I would be quite intrested to read into them. Thanks. $\endgroup$ Feb 2, 2012 at 1:17
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.