CODEWITHSUNDEEP
rfloating-pointprecisionfloating-accuracy
Ben
Ben

Reputation: 1556

R: Computing maximum floating-point error for log(exp(...))

I am working on some programming problems where I have to transition probabilities between standard space and log-space. For this purpose, I am trying to find out the maximum absolute error for the floating-point error in R for the computation log(exp(...)) where the input is a log-probability (i.e., a non-positive number).

At the moment I have computed the answer by using a grid search (see code and plot below), but I'm not sure if the value I've computed is correct. (I checked some other ranges but the range shown in the plot seems to get the maximum absolute error.)

#Set function for computing floating-point error of log(exp(...))
fp.error <- function(x) { abs(log(exp(x)) - x) }

#Compute and plot floating-point error over a grid of non-negative values
xx <- -(0:20000/10000)
ff <- fp.error(xx)
plot(xx, ff, col = '#0000FF10',
     main = 'Error in computation of log(exp(...))', 
     xlab = 'x', ylab = 'Floating-Point Error')

#Compute maximum floating-point error
fp.error.max <- max(ff)
fp.error.max
[1] 1.110223e-16

enter image description here

From this analysis, my estimated maximum absolute error is half the size of .Machine$double.eps (which is 2.220446e-16). I am not sure if there is a theoretical reason for this, or if I am getting the wrong answer.

Question: Is there any way to determine if this is really the maximum floating-point error for this computation? Is there any theoretical way to compute the maximum, or is this kind of grid-search method sufficient?

Upvotes: 5

Views: 370

Answers (3)

Sam Mason
Sam Mason

Reputation: 16213

In general, I'd suggest using a randomised approach to generate a lot more xs, e.g.:

x <- runif(10000000, 0, 2)

Your regularly spaced values might happen to stumble on a pattern that "just works".

Also it depends on whether you care about absolute error or relative error. The absolute error should be close to .Machine$double.xmax while relative error increases as x approaches zero. E.g. log(exp(1e-16)) gets truncated to zero.

Upvotes: 2

GKi
GKi

Reputation: 39717

The error of log(exp(x)) depends on the value of x. If you use a float also x has an precision which depending on its value. The precission could be calculated with nextafter from C:

library(Rcpp)
cppFunction("double getPrec(double x) {
  return nextafter(x, std::numeric_limits<double>::infinity()) - x;}")

getPrec(2)
#[1] 4.440892e-16

getPrec(exp(2))
#[1] 8.881784e-16

or not using Rcpp:

getPrecR <- function(x) {
  y <- log2(pmax(.Machine$double.xmin, abs(x)))
  ifelse(x < 0 & floor(y) == y, 2^(y-1), 2^floor(y)) * .Machine$double.eps
}

Have also a look at: What is the correct/standard way to check if difference is smaller than machine precision?.

Upvotes: 2

ThomasIsCoding
ThomasIsCoding

Reputation: 102680

I think you got the right answer. Here I refined the step as small as sqrt(.Machine$double.eps), you will see

> x <- seq(0, 2, by = sqrt(.Machine$double.eps))

> max(abs(log(exp(x)) - x))
[1] 1.110725e-16

However, once your x is extremely large, you will have Inf error, e.g.,

> (x <- .Machine$double.xmax)
[1] 1.797693e+308

> max(abs(log(exp(x)) - x))
[1] Inf

Upvotes: 2

Related Questions