3
$\begingroup$

Why does the finite-differences gradient/hessian approximation nlme::fdHess return NaNs when seeking the gradient/hessian at 0?

nlme::fdHess(0,function(x)x^2)  
> $mean  [1] 0
> 
>  $gradient  [1] NaN
> 
>  $Hessian
>       [,1]  [1,]  NaN

nlme::fdHess(1,function(x)(x-1)^2)  
> $mean  [1] 0
> 
>  $gradient  [1] -1.110224e-16
> 
>  $Hessian
>       [,1]  [1,]    2

I traced through the source and can indeed see a 0/0, but couldn't make sense of the logic (the man page says "This function uses a second-order response surface design known as a Koschal design", so it apparently isn't using the usual central differences type of approximation to the derivative). If this is just a bug, and not me fudging something, is there another well-tested finite-differences grad/hessian approximation implemented in R?

I don't know if this post belongs here or Stackoverflow or operations research or elsewhere but chose it here since the nlme package is mainly for statistics (even though the question has nothing to do with the main modeling part of nlme).

$\endgroup$

1 Answer 1

3
$\begingroup$

Looking at the code it seems as though the basic problem here is that the function, by default, chooses the finite difference step size relative to the absolute value of the parameters: in the following code fragment, the parameters have default values .relStep = .Machine$double.eps^(1/3), minAbsPar = 0.

This code effectively computes a $\Delta x$ factor incr, and a numerator frac that will rescale zeroth, first, and second differences:

incr <- pmax(abs(pars), minAbsPar) * .relStep
frac <- c(1, incr, incr^2)

So, when the max(abs(pars)) is 0 and minAbsPar is also 0, the function uses $\Delta x = 0$ ... which will lead to $(f(x+\Delta x) - f(x))/\Delta x = 0/0$ = NaN.

You could choose an alternative version of minAbsPar that makes sense for your problem (for reasonably scaled functions, setting minAbsPar = 1 would probably be safe).

The optimHess() function in base R uses the ndeps component of the control list (default value is 0.001) to set the step. Users can use the parscale component to change the scale.

numDeriv::hessian is another alternative (it uses Richardson extrapolation rather than simple finite differences).

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.