Implementation of smoothing splines function

Question

I have some problems in writing a module for spline smoothing. Actually, I have been 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 try to place a knot in every point in my data, numerical errors arise, 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 (let's 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?

Answer 1

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}, 
   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]]]

Answer 2

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 code 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 one compares 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.

Answer 3

Here is a (simplified) implementation of Reinsch's smoothing spline, which is effectively equivalent to csaps() in MATLAB's Curve Fitting Toolbox. Fancier methods have come along since then (e.g. Wahba's cross-validation splines), but this old workhorse has still proved serviceable:

SmoothingSplineFunction[dat_?MatrixQ, p : (_?NumericQ | Automatic) : Automatic] :=
 Module[{n = Length[dat], pv = p, cc, dc, del, h, qg, qm, rh, tm, uv, xa, ya},
        {xa, ya} = Transpose[dat]; h = Differences[xa]; rh = 1/h; 
        del = Differences[ya] rh;
        qm = SparseArray[{Band[{1, 1}] -> Most[rh], 
                          Band[{1, 2}] -> -ListCorrelate[{1, 1}, rh], 
                          Band[{1, 3}] -> Rest[rh]}, {n - 2, n}];
        tm = SparseArray[{Band[{2, 1}] -> Most[Rest[h]],
                          Band[{1, 1}] -> ListCorrelate[{2, 2}, h],
                          Band[{1, 2}] -> Drop[h, -2]}, {n - 2, n - 2}];
        qg = qm.Transpose[qm];
        If[pv === Automatic, pv = 1/(1 + Tr[tm]/(6 Tr[qg]))];
        uv = LinearSolve[6 (1 - pv) qg + pv tm, Differences[del]];
        dc =
        ya - 6 (1 - pv) Differences[ArrayPad[Differences[ArrayPad[uv, 1]]/h, 1]];
        Interpolation[Transpose[{List /@ xa, dc, Append[Differences[dc]/h -
                                 h ListCorrelate[{2, 1}, ArrayPad[pv uv, 1]],
                                 pv Last[uv] Last[h] -
                                 (Subtract @@ Take[dc, -2])/Last[h]]}],
                      InterpolationOrder -> 3, Method -> "Hermite"]]

(It might be a bit confusing that the output is an InterpolatingFunction[] even though no interpolation is being done; it's just that I wanted to use the built-in facility for evaluating piecewise Hermite cubics.)

Let's use it on Andy's example:

BlockRandom[SeedRandom[249304]; 
            data = Table[{i, RiemannSiegelZ[i] + Sin[i] +
                          RandomVariate[NormalDistribution[0, .2]]},
                         {i, 0, 25, 0.05}]]
smth = SmoothingSplineFunction[data, 9/10];
Plot[smth[x], {x, 0, 25}, PlotStyle -> Directive[Thick, Red], 
     Prolog -> {Blue, AbsolutePointSize[5], Point[data]}]

There is an ad hoc element to choosing the "smoothing parameter" in the second argument, however. There have been a number of proposals on how to pick the best smoothing parameter, but I haven't gotten around to evaluating them, as I've mostly used smoothing splines only for looking at the approximate trend.

Answer 4

Anton Antonov has implemented smoothing splines in his Quantile regression with B-splines package (direct link to the M-file). This post (duplicated in this thread) and this WTC2014 talk explain how can it be used. See also this post of mine for an example of use.