Accurately computes the logarithm of the sum of exponentials
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)).
logSumExp(lx, idxs = NULL, na.rm = FALSE, ...)
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.
Returns a numeric scalar.
This method is optimized for correctness, that avoiding underflowing. It is implemented in native code that is optimized for speed and memory.
Henrik Bengtsson
[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
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].
## 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))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.