# How to add large exponential terms reliably without overflow errors?

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.

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

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

• 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? Feb 2, 2012 at 0:53
• 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. Feb 2, 2012 at 0:57
• To the downvoter: would you mind letting me know what is wrong with my answer? Feb 13, 2012 at 0:06
• 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.
– a06e
Aug 4, 2016 at 17:43
• scicomp.stackexchange.com/q/24624/988
– a06e
Aug 4, 2016 at 19:39

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)

• 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. Feb 2, 2012 at 1:00
• @J.M. could you provide me with the names of those other methods, I would be quite intrested to read into them. Thanks. Feb 2, 2012 at 1:17

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.