Bsplines are fully non-local

[antoine-levitt]

This is easily seen by interpolating a vector with only zeros except for a single one: the basis function is exponentially decaying. If I understand correctly this is not standard Bspline behavior, which are compactly supported. This looks related to the "prefiltering" mentioned in the code (which solves a tridiagonal system and is therefore nonlocal). Is it possible to turn it off?

[tomasaschan]

Does this happen for all degrees of B-splines? Could you post some example code?

Although a tridiagonal solver in general is nonlocal, I think (IIRC) the solution should, if the math and implementation are correct, yield the coefficients for the compactly supported B-splines - if they don't we have a bug, a mistake in the mathematical derivations, or both...

[antoine-levitt]

Degrees 2 and up. For what it's worth Dierckx does the same thing.

using Interpolations
using Dierckx
N=20
x = linspace(0,1,N)
y = zeros(N)
y[div(end,2)] = 1
f_interp = Spline1D(x,y,k=2)
f_interp2 = scale(interpolate(y, BSpline(Quadratic(Free())), OnGrid()),x)

x = linspace(0,1,1000)
semilogy(x,abs([f_interp(xx) for xx in x]))
semilogy(x,abs([f_interp2[xx] for xx in x]))

edit: so it seems the splines have an extension (up to 1e-16) of 20 grid points.

[tomasaschan]

Hm. The non-zero coefficients you're seeing are all on the same order of magnitude as eps()? I think that's unavoidable; we don't have infinite precision, after all.

I was worried you were seeing perturbations that can not be explained by floating-point rounding errors.

[antoine-levitt]

No they're not, see the code. They decay exponentially from 1 to 1e-16 in about 20 grid points (eg ten grid points from the perturbation the error is 1e-8)

[tomasaschan]

I wasn't able to get the code to run; semilogy was not defined (and I've been out of the loop for a while so I'm not sure where it should come from). Plotting with plot(x, log(abs([f_interp2[xx] for xx in x]))) I see your point, though:

I'm not sure if this is an inherent limitation of higher-order B-splines, or if there's a problem in our derivations and/or code, though. Needs more investiagtion (for which I sadly don't have enough time on my hands - but the derivations are documented here if you want to dig into it yourself).

[antoine-levitt]

Sorry about the semilogy, it's a PyPlot thing.

I'm definitely not an expert in this, but from what I see the filtering (and hence the non-locality) is to make the spline interpolating, correct? So it might not be easy to have compactly supported interpolating splines (which is the sort of interpolations I'm used to in e.g. finite elements)

[dkarrasch]

Perhaps, this is because B-splines are not meant for interpolation but for fitting, after all? http://mathworld.wolfram.com/B-Spline.html

I suspect a classic spline interpolation (no B-splines!) as in Dierckx.jl (note that by default s = 0.0 in Spline1D, which corresponds to "an interpolating spline") is what is done here, which would explain the nonlocality. If that's the case, one should think about avoiding the "B-spline" terminology?

[ssfrr]

In the audio world when people refer to "cubic" interpolation we usually mean a cubic hermite spline where the derivatives at the endpoints of the region being interpolated are estimated just using the adjacent points (I think just averaging the 1st-order discrete derivatives). So if you're interpolating in the range [5,6] the support is just the range [4,7]. This gives an interpolation (the curve passes through the given points) that's continuous and once-differentiable. There's also a version where instead of computing the derivatives you construct a cubic curve that passes through all 4 points surrounding the interpolation location, but it gives you a less-smooth curve (non-differentiable).

I'm guessing maybe Interpolations.jl is doing something more clever to estimate the derivatives which is where the wider support comes in?

Here's a plot of the situation @antoine-levitt mentioned. Normally I'd expect a cubic spline interpolation to be zero outside of the range [8,12], like the green cubic line here.

The Interpolations version is better in the sense that it is closer to the ideal bandlimited continuous function (the sinc function), but I wonder if the more compactly-supported cubic interpolation would be useful to some folks and is likely much faster (at least for the initial construction).

using Interpolations
using Plots

# code adapted from http://paulbourke.net/miscellaneous/interpolation/
function cubic(A, i)
    y0 = i < 2 ? A[1] : A[floor(Int, i)-1]
    y1 = i < 1 ? A[1] : A[floor(Int, i)]
    y2 = i > length(A) ? A[end] : A[ceil(Int, i)]
    y3 = i > length(A)-1 ? A[end] : A[ceil(Int, i)+1]

    i -= floor(i)
    a0 = -0.5*y0 + 1.5*y1 - 1.5*y2 + 0.5*y3
    a1 = y0 - 2.5*y1 + 2*y2 - 0.5*y3
    a2 = -0.5*y0 + 0.5*y2
    a3 = y1

    a0*i^3 + a1*i^2 + a2*i + a3
end

A = zeros(20)
A[10] = 1
itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())

x = linspace(1, 20, 500)
plot(x, itp[x], label="itp")
plot!(x, x->sinc(x-10), label="sinc")
plot!(x, x->cubic(A, x, true), label="cubic")
scatter!(A, label="A")

[tomasaschan]

Hermite splines have been discussed before and would be a welcome addition to the package (see eg #6) but someone has to sit down and sort out the maths, and so far no one has.

(I have an idea for how to make it easier to add new behaviors like that and get access to things like scaling, extrapolation and ND-from-1D, but I can't find time to sit down and work with it. I suspect it's one of those things I'll think about now and then but maybe never actually do... :/)