Best parameter for smoothing splines

Question

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?

Answer 1

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

Some comments:

• 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.

Answer 2

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
  ]