[R] Beyond double-precision?

[Stavros Macrakis]

On Sat, May 9, 2009 at 12:17 PM, Berwin A Turlach
<berwin at maths.uwa.edu.au> wrote:
> log(H) = log(n) - log( 1/x_1 + 1/x_2 + ... + 1/x_n)
>...But we need to calculate the logarithm of a sum from the logarithms of the individual terms.

> ...The way to calculate log(x+y) from lx=log(x) and ly=log(y) ...
>  max(lx,ly) + log1p(exp(-abs(lx-ly)))

Agreed completely so far. But instead of calculating the logsum
pairwise, you can do it all in one go, which is both more efficient
and more accurate.

Here are some timing and accuracy measurements of the one-shot logsum
compared to the loop and the Reduce versions. (Full code at the bottom
of this email.) The vector version is much faster and much more
accurate in general.  There must be cases where the log1p method
increases accuracy, but I couldn't find them.

                -s

Large examples to test accuracy and speed

Test case: runif(1e+06)
  function. time    error
1    logsum 0.22 9.31e-16
2  logsum_s 0.15 9.31e-16
3  logsum_r 9.75 3.10e-13

Test case: rexp(1e+06)
  function.  time     error
1    logsum  0.21 -1.40e-15
2  logsum_s  0.15 -1.40e-15
3  logsum_r 10.13 -1.38e-14

Test case: abs(rnorm(1e+06))
  function.  time     error
1    logsum  0.24 -4.38e-16
2  logsum_s  0.14 -4.38e-16
3  logsum_r 10.01 -8.74e-14

Test case: rep(1, 1e+05)
  function. time    error
1    logsum 0.01 1.46e-16
2  logsum_s 0.01 1.46e-16
3  logsum_r 0.96 6.24e-14

Test case: rep(10^-(1:10), each = 10000)
  function. time     error
1    logsum 0.02  6.14e-16
2  logsum_s 0.01  6.14e-16
3  logsum_r 0.95 -6.96e-12

More accurate even for small cases

Test case: 1:100
  function. time     error
1    logsum    0 -3.60e-16
2  logsum_s    0 -3.60e-16
3  logsum_r    0  3.24e-15

Test case: abs(rnorm(100))
  function. time     error
1    logsum    0 -3.48e-16
2  logsum_s    0 -3.48e-16
3  logsum_r    0 -2.09e-15


##########################
# Fast, accurate sum in log space
#
logsum <- function(l) {
     maxi <- which.max(l)
     maxl <- l[maxi]
     maxl + log1p(sum(exp(l[-maxi]-maxl))) }
##########################

##########################
# Simpler, perhaps less accurate sum in log space
#
logsum_s <- function(l) {
     maxl <- max(l)
     maxl + log(sum(exp(l-maxl))) }
##########################



# Pairwise reduction
logsum_r <- function(x) Reduce( function(lx, ly) max(lx, ly) +
log1p(exp(-abs(lx-ly))), x )


function_names <- c("logsum","logsum_s","logsum_r")

logsum_test <- function(l) {
  cat("\nTest case:",deparse(substitute(l)),"\n")
  realsum <- sum(l)
  logl <- log(l)
  results <- times <- list()
  lapply( function_names, function(f) times[[f]] <<- system.time(
results[[f]] <<- getFunction(f)(logl))[1])
  data.frame(`function`=function_names,
             time=as.numeric(times),
             error=(exp(as.numeric(results))-realsum)/realsum )
}

set.seed(1)

cat("\n\nLarge examples to test accuracy and speed\n\n")
logsum_test(runif(1000000))
logsum_test(rexp(1000000))
logsum_test(abs(rnorm(1000000)))
logsum_test(rep(1,100000))
logsum_test(rep(10^-(1:10),each=10000))

cat("\n\nMore accurate even for small cases\n\n")
logsum_test(1:100)
logsum_test(abs(rnorm(100)))