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?