Smooth delta functions to compute Jacobian and Hessian from samples

If I could afford more evaluations of the function, the best will probably be to use a Smolyak grid to construct an interpolant and then use automatic differentation?