Smolyak/sparse grids aren’t a panacea. They don’t work well if the Hessian has lots of non-negligible off-diagonal elements.
In high dimensions, there is no free lunch for interpolation — there has to be something special about your function. That’s why I said that it may be easier just to figure out how to differentiate G directly.