Numeric gradient in custom op

Hello, I am quite new to the community, so please excuse if my question is too stupid.
I want to solve an inverse problem. More specifically, I want to infer a set of parameters from a number of simulated values using variational inference.
My current program works as follows:
I choose points in parameter space using sparse grids and evaluate my simulation in these points.
Afterwards, I calculate a B-Spline interpolant in these points.
The resulting B-Spline is now the surrogate I want to use with variational inference.
One big advantage of B-Splines are the easy-to-compute gradients, that can later on be used by a NUTS sampler for example.
So I define a custom theano op class for my B-Spline surrogate, where the perform method returns the value of my B-Spline in a given, numeric coordinate.
Unfortunately I can’t do the same thing with the grad(), since it needs to be an analytic function.
Is there an easy way to define the grad() - function in the same way, where I can simply return the numeric gradient in a given coordinate?

Help would be very much appreciated, thank you in advance.
Cheers,

VWagner

Isnt B-Spline has a gradient implementation: https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-derv.html