Subtracting very small probabilities - How to compute?

Question

This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $\ell_1 \geqslant \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-diffference of these probabilities, which we denote by:

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

Questions: How can you effectively compute this log-difference? Can this be done in the base R? If not, what is the simplest way to do it with package extensions?

Answer 1

To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:

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

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

$$\begin{equation} \begin{aligned} \ell_- &= \ln \big( \exp(\ell_1) - \exp(\ell_2) \big) \\[6pt] &= \ln \big( \exp(\ell_1) (1 - \exp(-(\ell_1 - \ell_2))) \big) \\[6pt] &= \ell_1 + \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 0 = -\infty$. Using the Maclaurin series expansion for $\ln(1-x)$ we obtain the formula:

$$\begin{equation} \begin{aligned} \ell_- &= \ell_1 - \sum_{k=1}^\infty \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 possible to compute this log-difference accurately in base R using the log1p function. 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-difference:

logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }

Implementation with VGAM package: Machler (2012) analyses accuracy issues in evaluating the function $\ln(1-\exp(-|x|))$, and suggests that use of the base R functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp function in the VGAM package. This gives you the an alternative function for the log-difference:

logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }

Answer 2

The following workaround is often very useful for these sorts of problems:

1. subtract the smaller of l1 and l2 from each of l1 and l2 (effectively, we have multiplied the probabilities by some constant z = exp(min(l1,l2)))
2. Now compute the sum using standard functions (you can use log1p and expm1 if you want).
3. Afterwards, add the quantity you subtracted in step 1.

A simple example:

l1 <- -2000 ## exp(-2000) is computational zero
l2 <- -2002
z <- min(l1,l2)
l1 <- l1 - z
l2 <- l2 - z
## now we are guaranteed that one of them is zero
y <- log(exp(l1) - exp(l2))
y + z

returns the correct value -2000.145, which is equal to -2000 + log(1-exp(-2)).