# Implementation of smoothing splines function

I have some problems by writing a module for spline smoothing. Actually i am trying for about two weeks. My listing is here:

SplSmooth[data_, knots_, lambda_, degree_] :=
Module[{M, Knots, NKnots, NBasis, X, Dsq, a},
M = Length@data;
Knots = Flatten@{Table[1, {i, 1, degree}], knots,Table[M, {i,1,degree}]};
NKnots = Length@Knots;
NBasis = NKnots - degree - 1;
X = Table[
Evaluate@BSplineBasis[{degree, Knots}, n, t] // N, {t, 1, M},
{n, 0, NBasis - 1}];
Dsq = Differences[X, 2];
a=Inverse[Transpose[X].X + lambda*Transpose[Dsq].Dsq // N].Transpose[X].data // N;
Return[X.a]
];


When I am trying to place a knot in every point in my data, numerical errors are arising such as:

Inverse::luc: Result for Inverse of badly conditioned matrix {{1.251,-0.1255,-0.251,0.0836667,0.0418333,0.,0.,0.,0.,0.,<<72>>},<<9>>,<<72>>} may contain significant numerical errors. >>

obviously the corresponding result is wrong (I can see it from the plot). It seems that the matrix to be inverted, is ill conditioned:

a = Inverse[Transpose[X].X + lambda*Transpose[Dsq].Dsq // N].Transpose[X].data // N;


but now comes the other problem. I use equidistant knots (lets say with 7 points distance) to overcome this problem. But then sometimes the algorithm works with:

Knots = Flatten@{Table[1, {i, 1, degree}], knots,Table[M, {i,1,degree}]};


and some other times works with

Knots = Flatten@{Table[1, {i, 0, degree}], knots,Table[M, {i,0,degree}]};


Now, I think that there is some kind of problem in BSplineBasis function.

Q: Can you spot the problem please? Or has anyone of you implemented a simillar function in the past with BSplineBasis function?

-
Always avoid calling Inverse when you can use LinearSolve. Solving a linear system is much faster and stabler than calculating an inverse. –  ssch Sep 29 '13 at 16:47
have you seen this? mathematica.stackexchange.com/… –  belisarius Sep 29 '13 at 17:00
Thank you very much ssch. The algorithm is more stable now. But I am still in doubt for the duplication of the first and the last point as knots. Anyway, you helped me a lot. –  jojosthegreat Sep 29 '13 at 17:09
Actually the problem persists: LinearSolve::luc: "Result for ... of badly conditioned matrix may contain significant numerical errors. –  jojosthegreat Sep 29 '13 at 17:21
Dear Belissarius thank you for your response. All these posts state that we must duplicate the external knots d times, where d is the degree of the spline. I know that fact, but why my second problem persists? Anyway, this correction seems to cure many instabilities, but not all: a = LinearSolve[ Transpose[X].X + lambda*Transpose[Dsq]. Dsq, Transpose[X].data, Method -> "Krylov"] // N; –  jojosthegreat Sep 29 '13 at 17:51

I'm not sure if this addesses all of the issues you are having but here is an implementation I put together some time ago that allows us to use LinearModelFit and BSplineBasis to do spline regression.

The benefit of this approach is that all of the properties of FittedModel are immediately available to us. This allows for checking for fit, residual diagnostics etc.

SplineModel[data_, deg_, knots_] :=
Block[{basis, allKnots, n, kmin, kmax},
n = Length[knots] + 2;
kmin = 0;
kmax = Ceiling[Max[data[[All, 1]]]] + 1;

basis =
Array[\[FormalX]^# &, deg + 1, 0]~Join~
Table[BSplineBasis[{deg, knots}, i, \[FormalX]],
{i, 0, Length[knots] - deg - 2}];

LinearModelFit[data, basis, \[FormalX]]
];


Lets generate some interesting data...

SeedRandom[249304]; data =
Table[{i,
RiemannSiegelZ[i] + Sin[i] +
RandomReal[NormalDistribution[0, .2]]}, {i, 0, 25, .05}];


And now we pick some knots and smooth the data using cubic splines.

knots = Range[0, 25, 1];
mod = SplineModel[data, 3, knots];

Show[ListPlot[data], Plot[mod[x], {x, 0, 25}, PlotStyle -> Directive[Red, Thick]]]


-
Thank you Andy. It seems that my only problem is the ill conditioning. I tried your approach and everything works fine, except from the absence of the smoothing parameter. Now, I think that in my approach, the term lambda*Transpose[Dsq].Dsq , witch is a faster alternative for roughness penalty, causes all the trouble. –  jojosthegreat Sep 30 '13 at 20:53
Dear Andy, I've noticed that you are using as basis the terms 1,x,...x^degree and also the BSplineBasis function. Namely, you are using the B-splines and some functions of the truncated power basis. Is that approach correct? I am asking you because I've never seen this before. –  jojosthegreat Oct 1 '13 at 7:21
What I have found is that adding these has the effect of minimizing edge effects. Feel free to remove these extra basis functions but also keep in mind that there is nothing inherently wrong with adding any basis functions we wish. –  Andy Ross Oct 1 '13 at 13:33

I have long been looking for a good implementation of Cubic Spline Smoothing with adjustable roughness penalty parameter for Mathematica. Your module gave me enough hints to understand how to make this work in Mathematica so I basically made a Cubic spline smoothing from your code with minor adjustments (about knots, a little bit about performance)

CubicSplSmooth[data_, lambda_] :=
Module[{M, Knots, X, Dsq, a},
M = Length @ data;
Knots = Flatten@{ 1, 1, 1, Range @ M, M, M, M};
X = Table[ Evaluate @ N @ BSplineBasis[{3, Knots}, n, t],
{t, 1, M}, {n, 0, M + 1}];
Dsq = Differences[X, 2];
a = LinearSolve[ Transpose[X].X + lambda*Transpose[Dsq].Dsq,
Transpose[X].data, Method -> "Multifrontal"];
Return[X.a]
];


This is restricted to Cubic degree but can be generalized to arbitrary degree as in your example. Manipulate is a nice way to get a feeling for the performance by moving the slider around:

Manipulate[
smoothdata = CubicSplSmooth[data, 10^lambda];
Show[ ListPlot[ data, PlotRange -> {-5, 3}],
ListLinePlot[ smoothdata, Mesh -> All, PlotStyle -> Red]],
{{lambda, 0}, -5, 5}]


The smoother behaves very naturally, yielding the original data for low (close to 0) values of lambda and a linear fit to data for extremely high ones.

If comparing this to the performance Labview achieves with the Cubic Spline Fit VI, it is still slower especially for large datasets. But the source of that is not accessible I think. Anyway it works well up to now but I think that performance can be surely improved.

-
Welcome to mathematica.SE and thanks for sharing :) –  ssch Oct 27 '13 at 13:56
The speed is the last thing that can bother us. We can use the Compile function, so the code will be transformed in the background into C code. And if Mathematica can reach gcc, then the result is really fast. –  jojosthegreat Oct 29 '13 at 11:48
Also I've found that LeastSquares is a faster and more reliable command for matrix inversion problems, so it is preferred than LinearSolve for arbitrary data. –  jojosthegreat Oct 29 '13 at 11:51
I think that the problem is the B-Spline functions themselves. They have the advantage to transform data in such way that all the math from least squares fitting can be applied without problem. But they still have problems with the stability and they aren't well documented. I am planning to use the "old school" approach from the book: Pollock, Handbook of Time Series Analysis, Signal Processing, and Dynamics. –  jojosthegreat Oct 29 '13 at 12:57
@Tobi In the line: X = Table[ Evaluate @ N @ BSplineBasis[{3, Knots}, n, t],{t, 1, M}, {n, 0, M + 1}]; n & t have no definition. Was there a reason for these or is this a typo? –  R Hall Jan 5 at 0:00