slow level set curvature via numDeriv with mgcv::gam thin plate spline

[aavogt]

I have two models. I would like to know the curvature of the level set:

ms <- d %>% group_by(f) %>% group_map( ~ gam(elapsed ~ s(n, delta, width), data= .) )

curvature <- function(p) {
        f <- function(q) { newdata <- as.list(q) %>% setNames(c("n", "delta", "width"))
                           predict(ms[[1]], newdata) - predict(ms[[2]], newdata)
        }
        H <- numDeriv::hessian(f, p)
        G <- numDeriv::grad(f, p)
        G %*% (det(H) * solve(H)) %*% G / (sum(G^2))**(length(p)/2+1)
}
tic()
curvature(c(1029,180,2))
toc()

Calculating the above curvature takes 1.5s, while predicting a single value from one gam takes 0.01s. I'm also worried about accuracy since the final and many intermediate values are small (1e-14).
Is there a model or package I could use with or instead of mgcv::gam to get these derivatives faster?