Adding very small probabilities—How to compute?

Question

In some problems, probabilities are so small that they are best represented in computational facilities as log-probabilities. Computational problems can arise when you try to add these small probabilities together, since some computational facilities (e.g., base R) can't distinguish very small probabilities from zero. In these cases it is necessary to solve the problem either by finding a workaround that avoids the small numbers, or by using a computational facility that can deal with very small numbers.

The archetypal problem of this kind is as follows. Suppose you have log-probabilities $\ell_1$ and $\ell_2$, where the corresponding probabilities $\exp(\ell_1)$ and $\exp(\ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R). We want to find the log-sum of these probabilities, which we denote by:

$$\ell_+ \equiv \ln \big( \exp(\ell_1) + \exp(\ell_2) \big)$$

Assume that we are—at least initially—working in an environment where $\exp(\ell_1)$ and $\exp(\ell_2)$ cannot be computed, since they are so small that they are indistinguishable from zero.

Questions: How can you effectively compute this log-sum? Can this be done in the base R? If not, what is the simplest way to do it with package extensions? (See this related question for extension to the vectorised case).

Answer 1

Some preliminary mathematics: As a preliminary matter it is worth noting that this function is the LogSumExp (LSE) function, which is often encountered when performing computation on values that are represented in the log-scale. To see how to deal with sums of this kind, we first note a useful mathematical result concerning sums of exponentials:

$$\begin{equation} \begin{aligned} \exp(\ell_1) + \exp(\ell_2) &= \exp(\max(\ell_1,\ell_2)) + \exp(\min(\ell_1,\ell_2)) \\[6pt] &= \exp(\max(\ell_1,\ell_2)) (1 + \exp(\min(\ell_1,\ell_2)-\max(\ell_1,\ell_2)) \\[6pt] &= \exp(\max(\ell_1,\ell_2)) (1 + \exp(-|\ell_1 - \ell_2|)). \\[6pt] \end{aligned} \end{equation}$$

This result converts the sum to a product, which allows us to present the log-sum as:

$$\begin{equation} \begin{aligned} \ell_+ &= \ln \big( \exp(\ell_1) + \exp(\ell_2) \big) \\[6pt] &= \ln \big( \exp(\max(\ell_1,\ell_2)) (1 + \exp(-|\ell_1 - \ell_2|)) \big) \\[6pt] &= \max(\ell_1, \ell_2) + \ln (1 + \exp(-|\ell_1 - \ell_2|)). \\[6pt] \end{aligned} \end{equation}$$

In the case where $\ell_1 = \ell_2$ we obtain the expression $\ell_+ = \ell_1 + \ln 2 = \ell_2 + \ln 2$, so the log-sum is easily computed. In the case where $\ell_1 \neq \ell_2$ this expression still reduces the problem to a simpler case, where we need to find the log-sum of one and $\exp(-|\ell_1 - \ell_2|)$.$^\dagger$

Now, using the Maclaurin series expansion for $\ln(1+x)$ we obtain the formula:

$$\begin{equation} \begin{aligned} \ell_+ &= \max(\ell_1, \ell_2) + \sum_{k=1}^\infty (-1)^{k+1} \frac{\exp(-k|\ell_1 - \ell_2|)}{k} \quad \quad \quad \text{for } \ell_1 \neq \ell_2. \\[6pt] \end{aligned} \end{equation}$$

Since $\exp(-|\ell_1 - \ell_2|) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $|\ell_1 - \ell_2|$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.

---

Implementation in base R: It is actually possible to compute this log-sum accurately in base R through creative use of the log1p function. (Alternatively, one can even use the log1pexp function in extension packages, which is essentially just a shifted version of logsumexp; here we will use only the base functions.) This is a primitive function in the base package that computes the value of $\ln(1+x)$ for an argument $x$ (with accurate computation even for $x \ll 1$). This primitive function can be used to give a simple function for the log-sum:

logsum <- function(l1, l2) { max(l1, l2) + log1p(exp(-abs(l1-l2))) }

Implementation of this function succeeds in finding the log-sum of probabilities that are too small for the base package to deal with directly. Moreover, it is able to calculate the log-sum to a high level of accuracy:

l1 <- -3006
l2 <- -3012

logsum(l1, l2)
[1] -3005.998

sprint("%.50f", logsum(l1, l2))
[1] "-3005.99752431486240311642177402973175048828125000000000"

As can be seen, this method gives a computation with 41 decimal places for the log-sum. It only uses the functions in the base package, and does not involve any changes to the default calculation settings. It gives as high a level of accuracy as you are likely to need in most cases.

It is also worth noting that there are many packages in R that extend the computational facilities of the base program, and can be used to deal with sums of very small numbers. It is possible to find the log-sum of small probabilities using packages such as gmp or Brobdingnag, but this requires some investment in learning their particular syntax.

---

$^\dagger$ From this result we can also see that if $|\ell_1 - \ell_2|$ is itself large (i.e., if one of the probabilities is very small compared to the other) then the exponential term in this equation will be near zero, and we will then have $\ell_+ \approx \max(\ell_1, \ell_2)$ to a very high degree of accuracy.