Numerical blowup problem in a fractional function in R

Question

Ciao,

I am working with this function in R:

betaFun = function(x){
  if(x == 0){
    return(0.5)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

The function is smooth and well defined for every x (at least from a theoretical point of view) and in 0 the limit approach to 0.5 (you can convince yourself about this by using Hopital theorem).

I have the following problem:

i.e. the fact that, due to the limit, R wrongly compute the values and I get a blowup in 0.

Here I report the numerical issue:

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13)  
sapply(x, betaFun)

[1] 5.000083e-01 5.000442e-01 2.220446e+00 0.000000e+00 0.000000e+00 1.111111e+10

As you can see the evaluation is pretty weird, in particular last one. I thought that I could solve this problem by defining the missing value in 0 (as you can see from the code) but it is not true.

Do you know how can I solve this numerical blow up problem?

I need high precision for this function since I have to invert it around 0. I will do it using nleqslv function from nleqslv library. Of course the inversion will return wrong solutions if the function has numerical problems.

I don't see that with `curve((1+exp(x)*(x-1) )/( x*(exp(x)-1)), -10, 10, n = 1e5)`. Please provide a full reproducible example. — Roland, Feb 08 '19 at 12:14
@Roland be aware that the points you select are only upto 1E-4. The OP computes 1E-13 — kvantour, Feb 17 '19 at 16:13
Could you please make the plot compatible with the formula used? For the domain `[-10,10]` as depicted, the given function does not come that close to the asymptotic lines, the asymptotic lines are at `y=0` and `y=1`, and the singularity is much tighter around `x=0`, where the value is `0.5`, not `0`. — Lutz Lehmann, Feb 18 '19 at 08:45

dmuir · Answer 1 · 2019-02-09T09:29:14.267

I think that you are losing accuracy in the evaluation of exp(x)-1 for x close to 0. In C if I evaluate your function as

double  f2( double x)
{   return (x==0)   ? 0.5
            : (x*exp(x) - expm1(x))/( x*expm1(x));
}

The problem goes away. Here expm1 is a math library function that computes exp(x) - 1, without losing accuracy for small x. I'm afraid I don't know if R has this, but you'd hope it would.

I think, though, that you would be better to test for |x| was sufficiently small, rather than 0.0. The point is that for small enough x both x*exp(x) and expm1(x) will be, as doubles, x, so their difference will be 0. To keep maximum accuracy may need to add a linear term to the 0.5 you return. I've not worked out precisely what 'sufficiently small should be, but it's somewhere around 1e-16 I think.

Yes, there is a `expm1` function in R. – Roland Feb 08 '19 at 13:37 — Roland, Feb 08 '19 at 13:37

Roland · Accepted Answer · 2019-02-08T14:29:16.060

Your problem is that you take the quotient of two numbers with very small absolute values. Such numbers are only represented to floating point precision.

You don't specify why you need these function values for x values close to zero. One easy option would be coercion to high precision numbers:

library(Rmpfr)  
betaFun = function(x){
  x <- mpfr(as.character(x), precBits = 256) 
  #if x is calculated, you should switch to high precision numbers for its calculation
  #this step could be removed then

  #do calculation with high precision, 
  #then coerce to normal precision (assuming that is necessary)
  ifelse(x == 0, 0.5, as((1 + exp(x) * (x - 1)) / (x * (exp(x) - 1)), "numeric"))
}  

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13, 0) 
betaFun(x)
#[1] 0.5000083 0.5000001 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000

Thx, is works fine except for x = 0. It is enough to add a if-statement to manually assign value at x = 0. — clarkmaio, Feb 08 '19 at 14:21

score 0 · Answer 3 · answered Feb 17 '19 at 16:27

As you notice, you are encountering the problem near zero. The roots of both the numerator and denominator are zero. And as the OP mentioned, using L'Hôpitcal, you notice that in that f(x) = 1/2.

From a numerical point of view, things go slightly different. Floating points will always have an error as not every Real number can be represented as a floating point number. For example:

exp(1E-3)  -1 = 0.0010005001667083845973138522822409868               # numeric
exp(1/1000)-1 = 0.001000500166708341668055753993058311563076200580... # true
                                  ^

The problem in evaluating numerically exp(1E-3)-1 already starts at the beginning, i.e. 1E-3

1E-3 = x   = 0.0010000000000000000208166817117216851
exp(x)     = 1.0010005001667083845973138522822409868
exp(x) - 1 = 0.0010005001667083845973138522822409868

1E-3 cannot be represented as a floating point, and is accurate upto 17 digits.
IEEE will give the closest floating point value possible to the true value of x, which already has an error due to (1). Still exp(x) is only accurate upto 17 digits.
By subtracting 1, we get a bunch of zero's in the beginning, and now our result is only accurate upto 14 digits.

So now that we know that we cannot represent everything exactly as a floating point, you should realize that near zero, it becomes a bit awkward and both numerator and denominator become less and less accurate, especially near 1E-13.

numerator_numeric(1E-13) = 1.1102230246251565E-16
numerator_true(1E-13)    = 5.00000000000033333333333...E-27

Generally, what you do near such a point is use a Taylor expansion around zero, and the normal function everywhere else:

betaFun = function(x){
  if(-1E-1 < x && x < 1E-1){
    return(0.5 + x/12. - x^3/720. + x^5/30240.)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

The above expansion is accurate upto 13 digits for x in the small region

Numerical blowup problem in a fractional function in R

3 Answers3