1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/logSumExp.R
\name{logSumExp}
\alias{logSumExp}
\title{Accurately computes the logarithm of the sum of exponentials}
\usage{
logSumExp(lx, idxs = NULL, na.rm = FALSE, ...)
}
\arguments{
\item{lx}{A \code{\link[base]{numeric}} \code{\link[base]{vector}}.
Typically \code{lx} are \eqn{log(x)} values.}
\item{idxs}{A \code{\link[base]{vector}} indicating subset of elements to
operate over. If \code{\link[base]{NULL}}, no subsetting is done.}
\item{na.rm}{If \code{\link[base:logical]{TRUE}}, missing values are
excluded.}
\item{...}{Not used.}
}
\value{
Returns a \code{\link[base]{numeric}} scalar.
}
\description{
Accurately computes the logarithm of the sum of exponentials, that is,
\eqn{log(sum(exp(lx)))}. If \eqn{lx = log(x)}, then this is equivalently to
calculating \eqn{log(sum(x))}.
}
\details{
This function, which avoid numerical underflow, is often used when computing
the logarithm of the sum of small numbers (\eqn{|x| << 1}) such as
probabilities.
This is function is more accurate than \code{log(sum(exp(lx)))} when the
values of \eqn{x = exp(lx)} are \eqn{|x| << 1}. The implementation of this
function is based on the observation that \deqn{ log(a + b) = [ la = log(a),
lb = log(b) ] = log( exp(la) + exp(lb) ) = la + log ( 1 + exp(lb - la) ) }
Assuming \eqn{la > lb}, then \eqn{|lb - la| < |lb|}, and it is less likely
that the computation of \eqn{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 \code{lx} and subtracts it from the remaining
values in \code{lx}.
}
\section{Benchmarking}{
This method is optimized for correctness, that
avoiding underflowing. It is implemented in native code that is optimized
for speed and memory.
}
\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))
}
\references{
[1] R Core Team, \emph{Writing R Extensions}, v3.0.0, April 2013. \cr
[2] Laurent El Ghaoui, \emph{Hyper-Textbook: Optimization Models
and Applications}, University of California at Berkeley, August 2012.
(Chapter 'Log-Sum-Exp (LSE) Function and Properties') \cr
[3] R-help thread \emph{logsumexp function in R}, 2011-02-17.
\url{https://stat.ethz.ch/pipermail/r-help/2011-February/269205.html}\cr
}
\seealso{
To compute this function on rows or columns of a matrix, see
\code{\link{rowLogSumExps}}().
For adding \emph{two} double values in native code, R provides the C
function \code{logspace_add()} [1]. For properties of the
log-sum-exponential function, see [2].
}
\author{
Henrik Bengtsson
}
|
