How to calculate sums in log-space without underflow?

In short, use the following expression:

fmax(la, lb) + log1p(exp(-fabs(lb - la)))

where I am using la and lb as the variables storing log(a) and log(b) respectively, and function names are from C's math.h. Most other languages will also have these functions, but they might be named differently, e.g. abs and max instead of fabs and fmax. (Note: These are conventions I will use throughout this answer.)

Speaking of functions, Mark Dickinson has a good point: you might want to check whether you have access to one that will do this for you directly. For example, if you are using Python and NumPy, you have access to logaddexp, and SciPy has logsumexp.

If you'd like more detail on where the above came from, how to add more than two numbers, and how to subtract, keep reading.

In more detail

There isn't a rule as simple as the ones for multiplication and division, but there is a mathematical identity that will help:
log(a + b) = log(a) + log(1 + b/a)

We can play with that identity a little to get the following:

log(a + b) = log(a) + log(1 + b/a)
           = log(a) + log(1 + exp(log(b) - log(a)))
           = la + log(1 + exp(lb - la))

There's still a problem here. If la is much greater than lb, you'll have 1 + 0.000...000something inside the log. There's not enough digits in the floating-point mantissa to store the something, so you'll just get log(1), losing lb entirely.

Luckily, most programming languages have a function in their standard library for solving just this problem, log1p, which calculates the logarithm of one plus its argument. That is, log1p(x) returns log(1 + x) but in a way that's accurate for very small x.

So, now we have:

log(a + b) = la + log1p(exp(lb - la))

We are almost there. There's one other thing to consider. In general, you want la to be greater than lb. It won't always matter, but sometimes this will gain you extra precision.* If the difference between lb and la is really big, this will save you from overflow in exp(lb - la). In the most extreme case, the calculation works when lb is negative infinity (i.e. b is 0), but not when la is.

Sometimes, you'll know which one is greater, and you can just use that one as la. But when it could be either one, you can use maximum and absolute value to work around it:

log(a + b) = fmax(la, lb) + log1p(exp(-fabs(lb - la)))

Sum of a collection

If you need to take the sum of more than two numbers, we can derive an expanded version of the identity above:

log(a[0] + a[1] + a[2] + ...)
    = log(a[0] * (a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...))
    = log(a[0]) + log(a[0]/a[0] + a[1]/a[0] + a[2]/a[0] + ...)
    = la[0] + log(1 + exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

We will want to use similar tricks as when we were adding two numbers. That way, we get the most accurate answer and avoid over- and underflow as much as possible. First, log1p:

log(a[0] + a[1] + a[2] + ...)
    = la[0] + log1p(exp(la[1]-la[0]) + exp(la[2]-la[0]) + ...)

The other consideration is which operand to pull out in front of the log1p. I used la[0] for demonstration so far, but you want to use whichever one is greatest. This is for all the same reasons we wanted la > lb when adding two numbers. For example, if la[1] were greatest, you would do the calculation like this:

log(a[0] + a[1] + a[2] + ...)
    = la[1] + log1p(exp(la[0]-la[1]) + exp(la[2]-la[1]) + ...)

Putting this into proper code, it would look something like this (this is C, but it should translate well to other languages):

double log_sum(double la[], int num_elements)
{
    // Assume index_of_max() finds the maximum element
    // in the array and returns its index
    int idx = index_of_max(la, num_elements);

    double sum_exp = 0;
    for (int i = 0; i < num_elements; i++) {
        if (i == idx) {
            continue;
        }
        sum_exp += exp(la[i] - la[idx]);
    }

    return la[idx] + log1p(sum_exp);
}

Calculating difference in log space

This wasn't part of the question, but since it may be useful nonetheless: subtraction in log space can be performed similarly. The basic formula is like this:

log(a - b) = la + log(1 - exp(lb - la))

Note that this is still assuming la is greater than lb, but for subtraction, it's even more important. If la is less than lb, you are taking the logarithm of a negative number!

Similar to addition, this has an accuracy issue that can be fixed by using specialized functions, but it turns out there are two ways. One uses the same log1p function as above, but the other uses expm1, where expm1(x) returns exp(x) - 1. Here are both ways:

log(a - b) = la + log1p(-exp(lb - la))
log(a - b) = la + log(-expm1(lb - la))

Which one you should use depends on the value of -(lb - la). The first is more accurate when -(lb - la) is greater than roughly 0.693 (i.e. log(2)), and the second is more accurate when it's less. For more detail on why that is and where log(2) came from, see this note from the R project evaluating the two methods.

The end result would look like this:

(lb - la < -0.693) ? la + log1p(-exp(lb - la)) : la + log(-expm1(lb - la))

or in function form:

double log_diff(double la, double lb)
{
    if (lb - la < -0.693) {
        return la + log1p(-exp(lb - la));
    } else {
        return la + log(-expm1(lb - la));
    }
}

---

* There is a bit of a sweet spot for this. When the difference between la and lb is small, the answer will be accurate either way. When the difference is too big, the result will always equal the greater of the two, since floating point numbers don't have enough precision. But when the difference is just right, you'll get better accuracy when la is greater.