R: Accurately computes the logarithm of the sum of exponentials
logSumExp
R Documentation
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)).
This function, which avoid numerical underflow, is often used when
computing the logarithm of the sum of small numbers (|x| << 1)
such as probabilities.
Usage
logSumExp(lx, idxs=NULL, na.rm=FALSE, ...)
Arguments
lx
A numericvector.
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, any missing values are ignored, otherwise not.
...
Not used.
Details
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.
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))