# Sum over very small exponentials: Underflow

I am trying to compute (in C) a sum like

$S = \sum_i \exp( - a_i )$,

where $10^{4} < a_i < 10^{5}$ are approximately normal distributed. So even if I do the Log-Sum-Exp trick

$S = \exp(\log[\sum_i \exp( - a_i + K )] - K)$,

with $K = \min_i(a_i)$, I have exponents of order $-10^5$. That means that almost every term vanishes during summation, since only the far right side of the normal-distribution contributes (only about $3$ points of $10^6$!). Is there any clever way to handle a sum like this? Thanks in advance for any suggestions.

• Are you worried about losing precision, or about inefficiently adding numbers that contribute negligibly to the final result? After all, if your sum is dominated by only a few terms, you only need to sum those terms and the rest can be safely ignored without affecting precision. – Kirill Dec 9 '15 at 22:41

You may use Kahan summation algorithm [1] The idea is to reschedule the sum operations in such a way precision loss is limited. The code is very simple (reproduced from [1] below). If this does not suffice, you may use multiprecision representations, such as quad doubles [2]. They are supported by several languages / compilers (including GNU c). Finally, if this does not suffice, you may use multi-precision arithmetics based on expansions [3]. I have written an implementation with a C++ class that can be used like a number [4]. Now if the underflow occurs within a single term (i.e., if the exp() underflows), then you may need to use floating point representations with larger number of bits. For instance, you can use MPFR [5]. It has implementations of all classical operators and functions with arbitrary precision (very good library, implemented by friends of mine :-).

function KahanSum(input)
var sum = 0.0
var y,t                      // Temporary values.
var c = 0.0                  // A running compensation for lost low-order bits.
for i = 1 to input.length do
y = input[i] - c         // So far, so good: c is zero.
t = sum + y              // Alas, sum is big, y small, so low-order digits of y are lost.
c = (t - sum) - y        // (t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
sum = t                  // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
next i                       // Next time around, the lost low part will be added to y in a fresh attempt.
return sum


[3] J. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete and Computational Geometry (Impact Factor: 0.69). 07/1996; 18(3). DOI: 10.1007/PL00009321

[4] https://gforge.inria.fr/frs/?group_id=5833, MultiPrecision_psm

[5] MPFR: http://www.mpfr.org/

• How is this going to help when computing the summands already underflows? At that point you are summing only zeros. – Federico Poloni Dec 4 '15 at 18:33
• Yes, you are right, if the underflow occurs right when computing the exp() and before the +, then more digits are required to compute the exp(). Can be done with MPFR (I am editing the answer). – BrunoLevy Dec 5 '15 at 14:15

You are adding positive terms, so you don't need to worry about the imprecision of the final result as if you had negative terms: (10^100+5)-(10^100+3)=2 but (10^100+5)+(10^100+3)=2*10^100 (as long you don't need 100 digits of precision).

In your case I'd sort the exponents in decreasing order (in the sense of magnitude). Then I'd keep adding just in case there was a carry from the smallest terms.

• Thanks for the hint with sorting! I can save quite a bit of time with that. – fastfforward Dec 3 '15 at 8:59
• Anyway, as $\exp(-10^{6}) \approx 3.56 \times 10^{-43430}$ and $\exp(-10^{5}) \approx 1.14 \times 10^{-4343}$ you'll need $10^{39000}$ samples of the smaller values to influence the sum of the 3 greater numbers, $10^{6}$ will make no influence. – N74 Dec 3 '15 at 13:46