Adding smoothing splines

[Kevin-Mattheus-Moerman]

@jipolanco Do you think it would be in scope for this project to add smoothing splines?

Here are some other/related implementations:
https://github.com/nignatiadis/SmoothingSplines.jl
https://github.com/robertdj/SmoothSpline.jl
https://github.com/JuliaMath/Dierckx.jl (Based on a Fortran wrapper)
https://github.com/francescoalemanno/KissSmoothing.jl

MATLAB implementation:
https://www.mathworks.com/help/curvefit/csaps.html

R Implementation:
https://www.rdocumentation.org/packages/DescTools/versions/0.99.49/topics/SmoothSpline

[jipolanco]

Hi Kevin, yes this is definitely in the scope of BSplineKit.jl and I've been wanting to add this for a while.

I'm not yet very familiar with how to implement smoothing splines based on a smoothing parameter, but I guess it shouldn't be too difficult. Also, I'm not yet sure about how they generalise to non-cubic splines (most implementations seem to be for cubic splines only). But that's ok, having cubic smoothing splines would be better than nothing.

As you probably know, an alternative approach consists in simply reducing the number of degrees of freedom (= number of knots) of the resulting spline compared to the input data. It would be quite easy to provide an interface for this. The problem is that the result is quite dependent on the choice of knots. I think this is sometimes called regression splines.

[Kevin-Mattheus-Moerman]

@jipolanco Great. Something like the R or MATLAB implementation there would be very useful! They other projects I linked to might offer insight into how to implement, and having it just for cubic splines for the moment seems fine. I am not an expert on the maths involved but happy to help create/add test cases etc.

[jipolanco]

In fact I got a working implementation in #103. It seems to give exactly the same results as the SmoothingSplines.jl package (even though I didn't look at their code in detail), which is a good sign. I just need to add some tests and then it should be merged.

[jipolanco]

The implementation is now included in the latest release, v0.18.0. It's quite simple to use, you can find some examples in the docs.

Let me know if you have any thoughts!

[Kevin-Mattheus-Moerman]

@jipolanco This is great thanks! I'll test it shortly!

What do you make of MATLAB's implementation with the $p$ parameter varying from 0 to 1. It feels like a handy way to blend between smooth and plain cubic, or is it mathematically mad?

[jipolanco]

@Kevin-Mattheus-Moerman In fact MATLAB's formulation is equivalent to the one we use as long as their $\lambda(t)$ is constant (by default it's constant and equal to 1). In that case, our $\lambda$ is equivalent to their $(1 - p) \lambda / p$.

What can be interesting about MATLAB's formulation is that their roughness parameter $\lambda(t)$ can vary over the curve. I don't know how useful this is, since the weights $w_j$ can already be used for more or less the same purpose, i.e. to choose regions where the curve needs to be particularly close to the data or particularly smooth.

A minor issue is that having an arbitrary $\lambda(t)$ function means that the roughness integral can no longer be computed exactly, but I think it would be OK to approximate $\lambda$ to a constant within each interval $[x_j, x_{j + 1})$ (which, by looking at their docs, is what MATLAB basically does).

In practice it wouldn't be too difficult to accept a variable $\lambda(t)$ as well as a $p$ parameter in BSplineKit's implementation, so I would be OK with supporting that if it's useful or convenient.

[Kevin-Mattheus-Moerman]

I have never used the varying $\lambda(t)$ but I find the use of $p$ quite useful. I think I understood what you said but I am not sure if the formulation is entirely equivalent under constant $\lambda$ i.e. the $p$ in front of the first term is not accounted for right? In the MATLAB formulation the $p$ helps to "lerp" between rough and smooth. I think the difference is clear also when we consider the case of $p=0$, this cannot be achieved with your formulation right?

[jipolanco]

Sorry I wasn't very clear. What I meant is that, if one takes MATLAB's expression (let's call it X) and then divides the whole thing by $p$, one ends up with a $\lambda (1 - p) / p$ coefficient in front of the integral (assuming $\lambda$ is constant), and the first term is the same as in BSplineKit's formulation (which is copied from other implementations). Note that minimising X or X/p is the same, since p is a constant from the point of view of the minimisation problem (we're optimising the shape of the spline for a fixed p), and therefore both formulations are equivalent.

I think the difference is clear also when we consider the case of $p = 0$, this cannot be achieved with your formulation right?

Well, if I'm not wrong the problem is actually ill-posed for $p = 0$. Looking again at MATLAB's expression, setting $p = 0$ means that the data points $(x_j, y_j)$ are simply not taken into account at all, and the problem has infinite solutions (any straight line). What is well defined is the limit $p \to 0$, in which case the spline should indeed tend to the "least-squares straight-line fit" quoting from the MATLAB docs. I guess MATLAB has a special implementation for the case $p = 0$ which directly solves the least-squares problem. In the current BSplineKit's implementation, setting a very large value of $\lambda$ should allow to get something close to that. In the future one could also imagine a special implementation for the case λ = Inf equivalent to MATLAB's behaviour.

[Kevin-Mattheus-Moerman]

Hmmm I get what you meant now. Your implementation seems to make more sense then. No strong feelings on either implementations. I'll test this now shortly. Actually, on that, I'm going to test it on closed loops. Do you have an example for these smoothing splines whereby the input is treated as periodic?

[jipolanco]

Well, the Periodic boundary condition is not yet implemented. It needs a bit more work. It should be easy to get it working, but not sure it will be very efficient at first (basically, I can't rely on BandedMatrices.jl due to out-of-bands terms).

[jipolanco]

Periodic conditions are now supported on v0.18.1. You can also find some examples here. Let me know if you have any remarks!

[Kevin-Mattheus-Moerman]

Fantastic! The force is strong with this one.

[Kevin-Mattheus-Moerman]

@jipolanco so you can see what it can be used for: smoothing contour lines for medical image segmentation. Works like a charm.

segment_01.mp4