Smoothed Spline Interpolation compatible with ForwardDiff?

Specific Domains Numerics
,
1

Hi,

I am wondering whether someone has insights on the following:
I have a complicated function inside which I perform smoothed (1-D) Cubic Spline interpolation using an irregular grid. Right now I am using the Dierckx.jl package for this.
In particular, the function calls something equivalent to

spl = Spline1D(x,y,k=3,s=1e-4)
y_eval = evaluate(spl, x_eval)

where x is an irregularly spaced grid and s the smoothing parameter.
Now, I wish to obtain a derivative of said function w.r.t a small number of input parameter. My usual go-to package for this would be ForwardDiff, but it is unfortunately not compatible with Dierckx.
So, I am wondering whether there are other packages with equivalent functionalities that would be. I know that e.g. Interpolations.jl or BasicInterpolators.jl are compatible with ForwardDiff, but they don’t seem to have the smoothing functionality.

I have some freedom in choosing the grid x so if there is a package that would require an regular grid, that may be helpful as well. Of course, as fallback I could resort to finite differences for getting the derivatives, but I would prefer ForwardDiff if possible.

Any help is very appreciated!

2

For one of my projects I had to do something similar.
I needed the first and second derivative of a (3D) cubic spline and was using forwarddiff.

I wrote my own implementation of the spline, and the tricky bit was letting derivatives work when selecting which bezier curve.

Essentially, if my spline is formed of 3 bezier curves, and I ask evaluate(spl, 1.2), it has to figure out “call evaluate(curves[2], 0.2)”. The problem was that doing modulo on a dual number didn’t work which made finding out the index hard.

My solution was as follows:

function get_curve_and_t(t::T1, s::CubicSpline{D, T2}) where {D, T1, T2}
    i = Int(ceil(t))
    t_curve = t-i+1
    curve = get_curve(s, i)
    return curve, i, t_curve 
end

Here, t is what you have called x. Comparing with my example, t=1.2, i==2 and t_curve==0.2

ceil when used on a dual number throws away the partials, which in this case is what we want. When we do t_curve - i , we get a dual number with the correct partials.

1 Like
3

You might have a look at ApproxFun : Home · ApproxFun.jl

This was also just posted in another thread: Other Interpolation Packages · Interpolations.jl

4

BSplineApprox should do it. It should be ForwardDiff compatible out of the box, but it also has an analytical derivative setup since differentiating the spline can be less numerically stable:

There are chain rules setup already exists https://github.com/SciML/DataInterpolations.jl/blob/master/ext/DataInterpolationsChainRulesCoreExt.jl, but it would be nice to setup Dual overloads so it uses the directly defined derivatives automatically.

1 Like
5

Thank you all for your help! Given what I want to do, DataInterpolations.jl seems like a good place to start and substantially simpler than an approach requiring an own implementation of the spline.

6

BSplineKit.jl has a keyword to obtain derivatives.

7

One relatively easy way to get a derivative for an arbitrary 1-dimensional cubic spline:

using BasicInterpolators
using ForwardDiff

# irregularly spaced knots
x = [0, 0.2, 1.3, 4]
# values to interpolate
y = [0.5, -2.6, -2, 0.2]

# initialize the spline
spline = CubicSplineInterpolator(x, y)

# evaluate the derivative at x = 1
ForwardDiff.derivative(spline, 1))

# evaluate the derivative at x = 1 and x = 3
map(x -> ForwardDiff.derivative(spline, x), [1, 3])