Numeric Hessian using gradient function R

Question

I need to calculate Hessian of my function numerically using my gradient function (programmed by formula, not numeric). Packages like numDeriv or rootSolve calculate hessian using numerical gradient that do not satisfy my needs. The thing I need performed internally (no separate method I can call) in optim package, however the only method that handles my optimization task well implemented in nlopt package and passing its optimal value to optim in order to get hessian is too expensive for my program.

So I need some function that calculates Hessian using not numeric gradient (see for example these formulas https://neos-guide.org/content/difference-approximations). I can't make such a function as I don't understand how to choose a parameter h (increment), that my function is very sensitive to. Could I find such a function in R or retrieve it somehow from optim package? Or may someone at least explain how to choose error minimizing value for h so then I will post this function by myself?

Answer 1

If I understand you correctly, you should use numDeriv::jacobian(), which takes a vector-valued function and computes a matrix (the derivative of each element with respect to each input), and apply it to your analytical gradient function. jacobian() does use numerical approximation (Richardson extrapolation, to be precise), but I don't see any other way that you could get from a black-box gradient function to a Hessian?

You do need to specify (or use the default value) of a numerical 'delta' function (1e-4 by default). On the other hand, the internal code that optim() uses to compute the Hessian also uses finite differences: see here and here ...

Below I define a function, its gradient function, and its Hessian; this code shows that jacobian(grad(x)) is the same as the Hessian (for a particular test case).

library(numDeriv)
x1 <- c(1,1,1)

Test that I didn't screw up the gradient and Hessian functions:

all.equal(grad(f,x1),g(x1)) ## TRUE
all.equal(h(x1),hessian(f,x1)) ## TRUE

Numerical equivalence of my Hessian and the Jacobian of the gradient:

all.equal(h(x1),jacobian(g,x1)) ## TRUE

Test functions:

f <- function(x) {
  sin(x[1])*exp(x[2])*x[3]^2
}

g <- function(x) {
  c(cos(x[1])*exp(x[2])*x[3]^2,
    sin(x[1])*exp(x[2])*x[3]^2,
    sin(x[1])*exp(x[2])*2*x[3])
}  

h <- function(x) {
    m <- matrix(NA,3,3)
    m[lower.tri(m,diag=TRUE)] <-
        c(-sin(x[1])*exp(x[2])*x[3]^2,
          cos(x[1])*exp(x[2])*x[3]^2,
          cos(x[1])*exp(x[2])*2*x[3],
    # col 2
           sin(x[1])*exp(x[2])*x[3]^2,
           sin(x[1])*exp(x[3])*2*x[3],
    # col 3
           sin(x[1])*exp(x[2])*2)
    m <- Matrix::forceSymmetric(m,"L")
    m <- unname(as.matrix(m))
    return(m)
}