logSumExp {matrixStats}

Accurately computes the logarithm of the sum of exponentials

Description

Accurately computes the logarithm of the sum of exponentials, that is, log(sum(exp(lx))). If lx = log(x), then this is equivalently to calculating log(sum(x)).

Usage

logSumExp(lx, idxs = NULL, na.rm = FALSE, ...)

Arguments

lx

A numeric vector. Typically lx are log(x) values.

idxs

A vector indicating subset of elements to operate over. If NULL, no subsetting is done.

na.rm

If TRUE, missing values are excluded.

...

Not used.

Details

This function, which avoid numerical underflow, is often used when computing the logarithm of the sum of small numbers (|x| << 1) such as probabilities.

This is function is more accurate than log(sum(exp(lx))) when the values of x = exp(lx) are |x| << 1. The implementation of this function is based on the observation that

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

Assuming la > lb, then |lb - la| < |lb|, and it is less likely that the computation of 1 + exp(lb - la) will not underflow/overflow numerically. Because of this, the overall result from this function should be more accurate. Analogously to this, the implementation of this function finds the maximum value of lx and subtracts it from the remaining values in lx.

Value

Returns a numeric scalar.

Benchmarking

This method is optimized for correctness, that avoiding underflowing. It is implemented in native code that is optimized for speed and memory.

Author(s)

Henrik Bengtsson

References

[1] R Core Team, Writing R Extensions, v3.0.0, April 2013.
[2] Laurent El Ghaoui, Hyper-Textbook: Optimization Models and Applications, University of California at Berkeley, August 2012. (Chapter 'Log-Sum-Exp (LSE) Function and Properties')
[3] R-help thread logsumexp function in R, 2011-02-17. https://stat.ethz.ch/pipermail/r-help/2011-February/269205.html

See Also

To compute this function on rows or columns of a matrix, see rowLogSumExps().

For adding two double values in native code, R provides the C function logspace_add() [1]. For properties of the log-sum-exponential function, see [2].

Examples

## EXAMPLE #1
lx <- c(1000.01, 1000.02)
y0 <- log(sum(exp(lx)))
print(y0) ## Inf

y1 <- logSumExp(lx)
print(y1) ## 1000.708


## EXAMPLE #2
lx <- c(-1000.01, -1000.02)
y0 <- log(sum(exp(lx)))
print(y0) ## -Inf

y1 <- logSumExp(lx)
print(y1) ## -999.3218


## EXAMPLE #3
## R-help thread 'Beyond double-precision?' on May 9, 2009.

set.seed(1)
x <- runif(50)

## The logarithm of the harmonic mean
y0 <- log(1 / mean(1 / x))
print(y0)  ## -1.600885

lx <- log(x)
y1 <- log(length(x)) - logSumExp(-lx)
print(y1)  ## [1] -1.600885

# Sanity check
stopifnot(all.equal(y1, y0))