# Best parameter for smoothing splines

I am trying to implement a module to determine the best value for the smoothing parameter of a smoothing spline interpolation of a given dataset. To calculate the cross-validation score of a given interpolation I have the following code:

CrossValidateScore[dat_?MatrixQ, p : (_?NumericQ | Automatic) : Automatic] :=
Module[{n = Length[dat], pv = p, i, cvs, xa, ya}, {xa, ya} = Transpose[dat];
Off[InterpolatingFunction::dmval];
cvs = 0;
For[i = 1, i <= n, i++,
cvs +=
Power[ya[[i]] - SmoothingSplineFunction[Delete[dat, {i}], p][xa[[i]]], 2]
];
cvs/n
]


where SmoothingSplineFunction is J.M. excellent implementation of Reinsch's smoothing spline provided here. Basically, what the code does is set aside the data (xi,yi) in turn, spline interpolate the remaining data, and sums fits to the left out values to get a cross–validated error sum of squares.

Now what I do is run through a certain range of parameters with the module

SplineCrossValidation[dat_?MatrixQ,pmax_,pmin_,psteps_] :=
Module[{n = Length[dat], pv, i, cvspl, stp},
cvspl = {};
stp = (pmax - pmin)/psteps;
For[i = 1, i <= psteps + 1, i++,
AppendTo[
cvspl, {pmin + (i - 1)*stp,
CrossValidateScore[dat, pmin + (i - 1)*stp]}]
];
cvspl
]


and then select the parameter p that minimizes the score, for example:

cvs = SplineCrossValidation[data, 1, 0.5, 1000];
bestsp = cvs[[Flatten[Position[Part[cvs , All, 2], Min[Part[cvs , All, 2]]]][[1]], 1]]


The problem is that this procedure becomes very slow pretty soon, especially when the size of the data set grows.

Any suggestion on how to improve performance please? Or has any of you implemented an efficient cross validation procedure?

A possible answer was actually suggested here for a 2D case. Let's generate some data first:

SeedRandom[2491304];
dat =  Table[{i,RiemannSiegelZ[20 i] + Sin[25 i] + RandomReal[NormalDistribution[0, .2]]},
{i, -0.5, 0.5, .001}];


Next, one can run the following code

(* construct knot interval *)
data = Sort[dat];
ll = Part[First@data, 1];
ul = Part[Last@data, 1];
intl = Abs[ul - ll];
If[ll == 0, ll = -intl 0.2, If[ll <  0 , ll = ll 1.2, ll = ll 0.8]];
If[ul == 0, ul = intl 0.2, If[ul <  0, ul = ul 0.8, ul = ul 1.2]];

(* construct knots and spline basis *)
nk = 50;
knots = Chop[Range[ll, ul, intl/nk]];
basis = Flatten@
Table[BSplineBasis[{3, knots}, i, x] , {i, 0, Length[knots] - 5}];

(* evaluate smoothing matrix *)
ff = Function[{x}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);

(* construct penalty matrix *)
s = SparseArray[{{i_, i_} -> -1, {i_, j_} /; i - j == 1 ->
2, {i_, j_} /; i - j == 2 -> -1}, {Length[basis], Length[a]}] //
Transpose;
pen = Transpose[s].s;

(* evalute best smoothing parameter *)
spdat = Table[at = a.Inverse[Transpose[a].a + 10^i pen].Transpose[a];
{i, ((IdentityMatrix[Length[a]] - at).(Last /@ data) // #.# &)/
Tr[IdentityMatrix[Length[a]] - at]^2}, {i, -3, 3, 0.01}];
ispdat = Interpolation[spdat, Method -> "Spline"];
pv = FindMinimum[{ispdat[x], -3 <= x <= 3}, x][[2, 1, 2]];


In this particular case, we should get a pv value of 0.602096 (25 seconds on my system, which is OK for my purposes). Then, finally:

ocss[x_] = basis.LinearSolve[Transpose[a].a + 10^pv pen, Transpose[a].(Last /@ data)];
Plot[ocss[x], {x, Part[First@data, 1], Part[Last@data, 1]}, PlotStyle -> Directive[Thick, Red],
Frame -> True, AspectRatio -> 0.75,  Prolog -> {Blue, AbsolutePointSize[5], Point[data]}]


giving

• Ideally I would like to put all the code above in a Module accepting as arguments the data points the number of nodes nk (which is hard coded above) and (optionally) a smoothing parameter p. However, if I try to do that the code hangs up. This looks like a problem in the scope of the x variable in basis, ff and a. Even though must be somehow elementary, I could not figure out how to solve it, and I would appreciate help on this point.
• It is not certain that the ispdat function will have a minimum in the hard-coded interval [-3,3], but I would not know how one could handle the minimum search in a general fashion.
• Have you seen this? mathematica.stackexchange.com/q/218764/1089 Apr 24 '20 at 11:35
• @chris Nope, I completely miss it. Thank you for the pointer, digging into it right now. Apr 24 '20 at 13:57
• it has a significant feature/possible drawback : the splines basis knots are set by the data itself. But conversely it would be useful to update it to find optimal regularisation by itself. Apr 24 '20 at 13:59

This corresponds to a quick and dirty encapsulation.

GCV[data_] :=
Module[{i, j, ll, ul, intl, nk, knots, basis, ff, a, s, pen, spdat,
ispdat, pv},
ll = Part[First@data, 1];
ul = Part[Last@data, 1];
intl = Abs[ul - ll];
If[ll == 0, ll = -intl 0.2, If[ll < 0, ll = ll 1.2, ll = ll 0.8]];
If[ul == 0, ul = intl 0.2, If[ul < 0, ul = ul 0.8, ul = ul 1.2]];
(*construct knots and spline basis*)
nk = 50;
knots = Chop[Range[ll, ul, intl/nk]];
(*evaluate interpolating matrix*)
ff = (Evaluate@
Flatten@Table[
BSplineBasis[{3, knots}, i, #], {i, 0,
Length[knots] - 5}]) &;
a = ff @@ # & /@ (Most /@ data);
(*construct penalty matrix*)
s = SparseArray[{{i_, i_} -> -1, {i_, j_} /; i - j == 1 ->
2, {i_, j_} /; i - j == 2 -> -1}, {Length[knots] - 4,
Length[a]}] // Transpose;
pen = Transpose[s].s;
(*evalute best smoothing parameter*)
spdat = Table[at = a.Inverse[Transpose[a].a + 10^i pen].Transpose[a];
{i, ((IdentityMatrix[Length[a]] - at).(Last /@ data) // #.# &)/
Tr[IdentityMatrix[Length[a]] - at]^2}, {i, -3, 3, 0.01}];
ispdat = Interpolation[spdat, Method -> "Spline"];
pv = FindMinimum[{ispdat[x], -3 <= x <= 3}, x][[2, 1, 2]];
pv
]