# logSumExp: Accurately computes the logarithm of the sum of exponentials In matrixStats: Methods that apply to rows and columns of a matrix

## 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

 `1` ```logSumExp(lx, na.rm=FALSE, ...) ```

## Arguments

 `lx` A `numeric` `vector`. Typically `lx` are log(x) values. `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. Analoguously 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.

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', https://inst.eecs.berkeley.edu/~ee127a/book/login/def_lse_fcn.html)
[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].
 ``` 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``` ```## 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)) ```