Smoothing splines (clamped spline)

Alan · July 17th, 2019, 7:15 pm

There are a zillion resources on splines, but I can't seem to find this.
I need an algorithm for a smoothing (cubic) spline where the endpoint derivatives are specified (the so-called 'clamped spline').

With no smoothing (so just ordinary interpolation), there are many treatments.
With so-called 'natural' endpoints (second derivative is zero) + smoothing, there are many treatments.

Trying to not reinvent the wheel.
Does anyone know a link or a cite?

Cuchulainn · July 17th, 2019, 9:32 pm

Alan.
My 2018 C++ book has a chapter on this + code. I remember a number of boundary conditons.
I don't have the book here. 1200 pages is a bit too much across the sea

BTW what do you mean by 'smoothing', precisely?

enum CubicSplineBC {SecondDeriv, FirstDeriv};

I have just sent you the code

The best references are Bulirch and Stoer and the magical Ahlberg, Nielsen and Walsh (in my bibliography). I have implemented them in my C++ code.

Alan · July 18th, 2019, 12:47 am

Hi Daniel,

Thank you for the code. I haven't opened it yet. Let me know if I will find the solution to the following problem, which will answer your question.

I want to find the piecewise cubic polynomial [$]g(t)[$] that solves the following problem:

(*) minimize [$] \left( \sum_{i=1}^n \{ y_i - g(t_i) \}^2 + \lambda \int_a^b \{ g''(u) \}^2 \, du \right)[$],

subject to
(A) given values for [$]g'(a)[$] and [$]g'(b)[$], where [$] a < t_1 < t_2 < \cdots < t_n < b[$],
(B) given data [$]\{t_i,y_i\}[$] and
(C) given [$]\lambda > 0[$], the 'smoothing' parameter.

Without (A), this is the standard smoothing spline problem, solved by Reinsch (1967) (fulltext here) and with many subsequent treatments.

BTW, I've got your 2018 book. Correct me if I'm wrong, but it looks like the cubic splines discussed there are just *interpolating* splines; i.e., the solution to (*) when [$]\lambda=0[$]?
.

Cuchulainn · July 18th, 2019, 10:23 am

That's clear, Alan. Thanks

So, you want SSSP + (A) together? In that case I reckon that the algorithm can be modified as an optimisation problem. But instead of solving a tridiagonal matrix system we must minimise a nonlinear function for the coefficients of the cubic spline. Looks benign enough especially since the integral of the 2nd derivative can be computed exactly.
The 'easiest' would be a derivative-free method, e.g. DE.

Does this make sense?

Cuchulainn · July 18th, 2019, 1:19 pm

Here is Python code (MM can also do it for sure). The bigger the smooth the better the fit.

import numpy as np
import matplotlib.pyplot as plt

import csaps

np.random.seed(1234)

x = np.linspace(-5., 5., 25)
y = np.exp(-(x/2.5)**2) + (np.random.rand(25) - 0.2) * 0.3

sp = csaps.UnivariateCubicSmoothingSpline(x, y, smooth=0.85)

xs = np.linspace(x[0], x[-1], 150)
ys = sp(xs)

plt.plot(x, y, 'o', xs, ys, '-')
plt.show()

FaridMoussaoui · July 18th, 2019, 1:36 pm

Alan, check "Numerical Analysis" by Burden, Faires & Reynods (1981). There is something about a clamped cubic spline algorithm.

rmodafferi · July 18th, 2019, 3:21 pm

Hi Alan, you can find the theory in Green & Silverman (1994), and Turlach (2005) the theory and a very readable implementation in Matlab in Martinez & Martinez (2016)

Alan · July 18th, 2019, 5:30 pm

Thanks, guys but it looks like you are all citing either (i) cubic spline interpolation (the interpolant goes through the data points) or (ii) standard smoothing splines (Green & Silverman, for example, which have a so-called "natural" boundary condition [$]g''(a) = g''(b) = 0[$]).

I am asking about what is apparently a non-standard case: smoothing splines with a first derivative condition at endpoints. It's problem (*) above with conditions (A) + (B) + (C).

I will guess this non-standard case has been solved, but if I can't find a reference in a week or so, I'll try it from scratch.

FaridMoussaoui · July 18th, 2019, 8:40 pm

Well, as your are looking for [$]N[$] locally cubic spline, each spline is defined by [$]4[$] parameters. (e.g [$]4 N[$] unknowns).
The cubics need to match at the knots points ([$] 2 N[$] equations). Make also the first derivative & second derivatives continuous at the interior knots points ([$] N -1[$] equations each). That [$] 4 N - 2[$] equations.
The "natural" boundary conditions or your set first derivatives provide two more equations. That's [$]4 N [$] equations.

Cuchulainn · July 19th, 2019, 3:00 pm

My feeling is that I understand the requirement * + A + B + C.
The solver compute the coefficients of the cubic spline. Natural BCs are not very difficult.

So put Farid's 4N unknowns in a (some) NL solver, yes? I don't think it is all that difficult but you probably need to do preprocessing of the input.

FaridMoussaoui · July 19th, 2019, 4:30 pm

For Alan specific problem with the smoothing parameter, the problem was already solved in
Carl de Boor (1978), A Practical Guide to Splines (Chapter XIV). While it is solved with "natural" BC
(allowing to have a tridiagonal matrix and adding 2 equations to the linear system), solving for your specified equations adds also 2 equations and now a band matrix.

Just let me know if you want to receive a copy of the book at your optioncity email.

Cuchulainn · July 19th, 2019, 4:39 pm

Well, as your are looking for [$]N[$] locally cubic spline, each spline is defined by [$]4[$] parameters. (e.g [$]4 N[$] unknowns).
The cubics need to match at the knots points ([$] 2 N[$] equations). Make also the first derivative & second derivatives continuous at the interior knots points ([$] N -1[$] equations each). That [$] 4 N - 2[$] equations.
The "natural" boundary conditions or your set first derivatives provide two more equations. That's [$]4 N [$] equations.

I get (N = number of intervals)

4N unknowns

N+1 interpolation points
3(N-1) cnty conditions
2 BCs

So, the same.

Cuchulainn · July 19th, 2019, 4:41 pm

Well, as your are looking for [$]N[$] locally cubic spline, each spline is defined by [$]4[$] parameters. (e.g [$]4 N[$] unknowns).
The cubics need to match at the knots points ([$] 2 N[$] equations). Make also the first derivative & second derivatives continuous at the interior knots points ([$] N -1[$] equations each). That [$] 4 N - 2[$] equations.
The "natural" boundary conditions or your set first derivatives provide two more equations. That's [$]4 N [$] equations.

I get (N = number of intervals)

4N unknowns

N+1 interpolation points
3(N-1) cnty conditions
2 BCs

So, the same.

Alan · July 19th, 2019, 4:46 pm

For Alan specific problem with the smoothing parameter, the problem was already solved in
Carl de Boor (1978), A Practical Guide to Splines (Chapter XIV). While it is solved with "natural" BC
(allowing to have a tridiagonal matrix and adding 2 equations to the linear system), solving for your specified equations adds also 2 equations and now a band matrix.

Just let me know if you want to receive a copy of the book at your optioncity email.

Thanks, Farid. I've got that book and will (re)-read that chapter with your comment in mind. I bought these books (de Boor and Green & Silverman) many years ago for a project, but now need to revisit them for this variant problem.

jherekhealy · July 21st, 2019, 10:13 pm

Yes, the BFIT routine of de Boor will do what you're after. FYI, this is part of the NSWC Fortran library https://raw.githubusercontent.com/jacob ... src/nswc.f

Smoothing splines (clamped spline)

Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)

Re: Smoothing splines (clamped spline)