Sum over very small exponentials: Underflow

Question

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.

Answer 1

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

[1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm

[2] https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format

[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/

Answer 2

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.