2
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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

enter image description here

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.
| improve this answer | |
$\endgroup$
  • $\begingroup$ Have you seen this? mathematica.stackexchange.com/q/218764/1089 $\endgroup$ – chris Apr 24 at 11:35
  • $\begingroup$ @chris Nope, I completely miss it. Thank you for the pointer, digging into it right now. $\endgroup$ – Daniele Binosi Apr 24 at 13:57
  • $\begingroup$ 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. $\endgroup$ – chris Apr 24 at 13:59
0
$\begingroup$

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
  ]
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.