Introduction

If numbers are stored in log space, you cannot add them by simply adding the logs. The correct operation for log(a + b) is the log-sum-exp transformation, which preserves numerical stability without forcing you to work directly with huge or tiny raw values.

This matters in statistics, probabilistic models, hidden Markov models, Bayesian inference, and softmax calculations. In those domains, direct exponentiation can overflow, underflow, or throw away most of the useful precision.

Why Plain Exponentiation Is Risky

Suppose you have x = log(a) and y = log(b). The naive expression is:

• compute exp(x)
• compute exp(y)
• add them
• take log again

That works mathematically, but not always numerically. If x or y is very large, exp may overflow. If they are very negative, exp may underflow to zero and erase meaningful relative differences.

The Log-Sum-Exp Identity

The stable identity is:

log(exp(x1) + exp(x2) + ... + exp(xn)) = m + log(sum(exp(xi - m)))

where m is the maximum log value in the list.

Subtracting m rescales the terms so every exponent is less than or equal to 1. That avoids overflow and keeps the smaller contributions representable.

A Runnable Python Implementation

1import math
2
3
4def logsumexp(values):
5    if not values:
6        raise ValueError("values must not be empty")
7
8    m = max(values)
9
10    if math.isinf(m):
11        return m
12
13    total = sum(math.exp(v - m) for v in values)
14    return m + math.log(total)
15
16
17logs = [-1000.0, -1001.0, -999.5]
18print(logsumexp(logs))

This code is stable even when the original probabilities are far too small to represent directly in standard floating-point form.

Two-Value Version

For only two terms, the formula is especially compact:

log(a + b) in log space becomes:

max(x, y) + log(1 + exp(-abs(x - y)))

where x = log(a) and y = log(b).

That form is useful in tight loops and dynamic programming tables because it avoids constructing temporary lists.

1import math
2
3
4def logaddexp(x, y):
5    m = max(x, y)
6    return m + math.log1p(math.exp(-abs(x - y)))
7
8
9print(logaddexp(-1000.0, -1001.0))

Using log1p improves precision when the second term is tiny relative to the first.

When Libraries Are Better

If you are already using numerical libraries, prefer the built-in implementation. For example, SciPy provides scipy.special.logsumexp, and numerical frameworks such as NumPy, PyTorch, and TensorFlow expose similar functionality.

Library implementations usually handle:

• vectorized inputs
• infinities
• axis reductions
• better low-level performance

The custom implementation is still worth understanding because it explains why the stable formula works.

Important Edge Cases

The most common special case is a list containing only negative infinity values, which corresponds to a sum of zero in the original domain. In that case, the log result should remain negative infinity rather than producing a misleading finite number.

You should also be careful with mixed finite and infinite inputs. A positive infinity term dominates the sum immediately, while negative infinity terms simply contribute nothing. Good library implementations already handle these cases correctly, which is another reason to prefer them for production work.

Common Pitfalls

• Adding log values directly when the goal is to add the original quantities.
• Exponentiating huge negative values and silently underflowing to zero.
• Forgetting to subtract the maximum term before exponentiation.
• Reimplementing the formula without log1p in cases where precision near zero matters.
• Ignoring special cases such as all values being negative infinity.

Summary

• Summing quantities in log space requires log-sum-exp, not plain addition of logs.
• The stable trick is to subtract the maximum log value before exponentiating.
• This avoids overflow and reduces underflow risk.
• A two-value logaddexp form is useful for tight loops and dynamic programming.
• Prefer library implementations in production, but understand the identity behind them.