Adding Regularisation to Built-in Interpolation function

Question

tl;tr; How to improve performance of regularised Interpolation?

Context

Since Version 12, Mathematica now incorporates a range of (underrated IMHO) regularisation methods to Fit and FindFit.

The option FitRegularization is particularly useful to regularise a fit, i.e. allow the BSpline basis not to go through the data points but provide a smoother fit. More generally these two functions also provide different NormFunction, which specify what norm[residual] should be minimised when doing the fit, which is also useful and general.

Hence I find it would be great that these two options be added to Interpolation transparently so that when interpolating, the option of not going exactly through the points can be controlled, e.g. if the purpose is to later differentiate the interpolation function.

(fairly successful) Attempt

Hence I have wrapped a new function, RegularisedInterpolation, inspired by this (elegant) post by @J.M. which does the interpolation explicitly using Fit and therefore inherits these functionalities.

Clear[RegularisedInterpolation];
RegularisedInterpolation[dat_, opt : OptionsPattern[]] := 
 Module[{fspl, fb, sd, scpts, sk, var, dim, x, y, z , t},
  fspl = Interpolation[dat, Method -> "Spline", 
    FilterRules[{opt}, Options[Interpolation]]];
  If[Length[FilterRules[{opt}, Options[Fit]]] == 0, Return[fspl]];
  fb = First[Cases[fspl, _BSplineFunction, \[Infinity]]];
  {sd, scpts, sk} = fb /@ {"Degree", "ControlPoints", "Knots"};
  dim = Last@Dimensions[dat];
  var = Which[dim == 2, {x},
    dim == 3, {x, y}, dim == 4, {x, y, z}, dim == 5, {x, y, z, t}];
  Compile[var // Evaluate, 
   Fit[dat, 
      Flatten[Outer[Times, 
         Sequence @@ MapThread[
           Table[BSplineBasis[{#1, #2}, k - 1, #3], {k, #4}] &, {sd, 
            sk, var // Evaluate, Dimensions[scpts]}]]] // Evaluate, 
      var // Evaluate,
      FilterRules[{opt}, Options[Fit]]
      ] // Chop // Evaluate]
  ]

It works as expected (in arbitrary dimension) as I will illustrate, up to 3 glitchs for which I am seeking help.

Validation

---

Let's first look at one dimensional data

 dat1 = Table[{i, Sin[i^2]}, {i, 0, 2 Pi, Pi/64}] // N;
 dat1 = dat1 /. {x_, y_} :> {x, y (1 + RandomVariate[NormalDistribution[0, 0.05]])};

First note that by default it does a standard interpolation.

g1 = RegularisedInterpolation[dat1];

Then I can add options from Fit and from Interpolation

gr = RegularisedInterpolation[dat1, 
FitRegularization -> {"Curvature", 10^-0.5}, InterpolationOrder -> 3]

If I plot the 2 Interpolations (pink and yellow) they look fairly similar

Show[Plot[{Sin[x^2], gr[x], g1[x]}, {x, 0, 2 Pi}], ListPlot[dat1], PlotRange -> All]

but when I differentiate it the regularised version (in pink) remains closer to the un-noised data (in red), when compared to the (yellow) interpolation.

Plot[{2 x Cos[x^2], D[gr[x], x], D[g1[x], x]} // Evaluate, {x, 0,  2 Pi}]

This is why I believe regularised interpolation is useful (even though it is a bit of a contradiction in terms).

---

Let's now look at two or three dimensional data

dat2 = Flatten[
Table[{i, j, Sin[i*j]}, {i, 0, 5, 0.5}, {j, 0, 5, 0.5}], 1] // N;
g2 = RegularisedInterpolation[dat2]
g2r = RegularisedInterpolation[dat2, 
  FitRegularization -> {"Curvature", 10^-0.5}]
pl2 = Plot3D[{D[g2[x, y], x],D[g2r[x, y], x] } // Evaluate, {x, 0, 5}, {y, 0, 5}, 
   PlotStyle -> Opacity[0.4], PlotPoints -> 20, PlotRange -> All];

Note the difference between the two curves: one is slightly smoother than the other.

Similarly in 3D

dat3 = Flatten[
 Table[{i, j, k, Sin[i*j*k]}, {i, 0, 2, 0.5}, {j, 0, 2, 0.5}, {k, 0, 2, 0.5}], 2];
g3 = RegularisedInterpolation[dat3, 
 FitRegularization -> {"Curvature", 10^-7.5}];

The regularised fit does not go exactly through the 3D data (as it should)

dat3 - Flatten[
 Table[{x, y, z, g3[x, y, z]}, {x, 0, 2, 0.5}, {y, 0, 2, 0.5}, {z,
    0, 2, 0.5}], 2] // Transpose // Last // ListPlot

While I am fairly happy with the result, I still have 3 issues

---

Questions

I have in fact three questions:

1. The code, while working complains about this

Experimental`NumericalFunction::dimsl: {y} given in {x,y} 
 should be a list of dimensions for a particular argument.

I have no idea what the issue is. Does anyone?

2. The performance is not very good in 2 and 3 dimensions when evaluating the fit (not when doing the inversion).

pl2 = Plot3D[D[g2[x, y], x] // Evaluate, {x, 0, 5}, {y, 0, 5}, 
    PlotStyle -> Opacity[0.4], PlotPoints -> 20, 
    PlotRange -> All]; // Timing

(* {0.179184,Null} *)

pl2r = Plot3D[D[g2r[x, y], x] // Evaluate, {x, 0, 5}, {y, 0, 5}, 
    PlotStyle -> Opacity[0.4], PlotPoints -> 20, 
    PlotRange -> All]; // Timing

(* {2.12889,Null} *)

It may be because the plotting routine refine on the boundaries of the splines? Would you know how to get back to Interpolation's native performance?

3. I would ideally prefer that the result of the fit show the same kind of wrapper as Interpolation (specifying the boundary of the interpolation)

I.e. it should return something like this:

Do you have any idea how to 'hide' information in this manner? Through associations??

Thank you for your help. I am hoping that many people will find this generic wrapper useful eventually.

---

Comments/ improvements

• It was drawn to me by @Henkik that the formulation is semantically incorrect since interpolation is meant to refer to 'going through the points'. The present function could/should be called BSplineRegularisedFit. On the other hand it has the same structure as Interpolation, and could seamlessly be integrated to the built-in function by simply increasing the number of options.

• @Kuba pointed out that this and that question addresses the wrapping issue. He also pointed out that

  Needs@"GeneralUtilities`" PrintDefinitions@InterpolatingFunction

yields the wrapper for InterpolatingFunction.

• This post shows how it could be generalised to un evenly sampled data as well. Here the main advantage is to inherit all the nice properties of Interpolation.

• It might be useful to be able to impose extra knots at the boundary of the domain?

Complement

To illustrate the benefit of regularisation let's add a glitch to the data and see how different penalty operate while changing the method and its amplitude.

Table[Plot[
   Table[RegularisedInterpolation[dat1, 
       FitRegularization -> {meth, 10^i}][x], {i, -2, 2, 0.5}] // 
    Evaluate,
   {x, 0, 2 Pi}], {meth, {"Curvature", "Variation", "TotalVariation", 
    "LASSO"}}] // Partition[#, 2] &

---

To illustrate the issue with the end point let us consider the successive derivative of the sine wave

dat1 = Table[{i, Sin[3 i]}, {i, 0, 2 Pi, Pi/8}] // N;
Clear[g]; 
g[x_] = RegularisedInterpolation[dat1, 
   FitRegularization -> {"Curvature", 10^-9.5}, 
   InterpolationOrder -> 12][x];
Table[Plot[D[{g[x], Sin[3 x]}, {x, i}] // Evaluate, {x, 0, 2 Pi}]
  , {i, 0, 5}] // Partition[#, 2] &

The match at the edge becomes worse with order of the derivative.

Answer 1

Regarding item 2/

A possible performance improvement (which still not scale very well, so please feel free to provide better answers!). is achieved by re-ordering the PieceWise function behind the spline.

Let us start with some 1D data and find the Regularised Interpolation

dat1 = Table[{i, Sin[3 i]}, {i, 0, 2 Pi, Pi/32}] // N;
Clear[g]; g[x_] = RegularisedInterpolation[dat1, 
  FitRegularization -> {"Variation", 10^-3.5}, InterpolationOrder -> 3][x];

Then If I reorder the PieceWise function

tt = PiecewiseExpand /@ g[x] // Simplify;
tt1 = Sort[Select[tt[[1]], FreeQ[#[[2]], Equal] &], #1[[2, 1]] < #2[[2, 1]] &] /. 
   Less -> LessEqual;
g2 = Compile[x, Piecewise[tt1] // Evaluate]

then the performance when plotting is quite improved:

Plot[g2[x], {x, 0, 2 Pi}]; // Timing

(* 0.024671` *)

over

Plot[g[x], {x, 0, 2 Pi}]; // Timing

(* 0.296903 *)

---

In 2D it works as well (but the Simplify takes a little while)

dat2 = Flatten[
    Table[{i, j, Sin[i*j]}, {i, 0, 5, 0.5}, {j, 0, 5, 0.5}], 1] // N;
Clear[g]; g[x_, y_] = RegularisedInterpolation[dat2, 
    FitRegularization -> {"Curvature", 10^-2.5}][x, y];
tt = PiecewiseExpand /@ g[x, y] // Simplify;
tt1 = Sort[Select[tt[[1]],FreeQ[#[[2]], 
     Equal] &], #1[[2, 1]] < #2[[2, 1]] &] /. Less -> LessEqual;
g2 = Compile[{x, y}, Piecewise[tt1] // Evaluate]

Then

  Plot3D[g2[x, y], {x, 0, 5}, {y, 0, 5}, 
  PlotStyle -> Opacity[0.4], PlotPoints -> 20, PlotRange -> All]

(* 0.268701 *)

is 16 times faster than

Plot3D[g[x, y], {x, 0, 5}, {y, 0, 5}, PlotStyle -> Opacity[0.4],
     PlotPoints -> 20, PlotRange -> All]; // Timing

(* 8.43557 *)

So one can add a option PiecewiseSimplify to RegularisedInterpolation, see the full code below.

---

Regarding item 3, following closely this post the trick is to define a new head, RegularisedInterpolatingFunction and return an Association which contains more than just the compiled code (Domain, Regularisation method and parameter etc..). One adds the rule that

  RegularisedInterpolatingFunction[a_Association][b__] := a["code"][b]

i.e. that the Object applied to some data or symbol b applies the Piecewise code a["code"]. Then we can use BoxFormArrangeSummaryBox BoxForm MakeSummaryItem to wrap a Summary box around the RegularisedInterpolationFunction object.

 gr = RegularisedInterpolation[dat1, FitRegularization -> {"Curvature", 10^-0.5}]

which can be expanded as

---

Finally I have added RegularisedListInterpolation which as ListInterpolation takes tables as arguments as in

  dat=Table[Sin[x^2+y]//N,{x,1,5,0.5},{y,1,5,0.5}];
  gr = RegularisedListInterpolation[dat, 
     FitRegularization -> {"Curvature", 10^-0.5}];

It should work in dimensions 1 to 4.

---

The full code now reads

ClearAll[RegularisedInterpolation, RegularisedListInterpolation,
RegularisedInterpolatingFunction];

RegularisedInterpolation::usage="Works like Interpolation but also inherits
options from Fit including FitRegularization.
 Example: \n dat1=Table[{x,Sin[x ]},{x,0,2,0.2}];
gr1=RegularisedInterpolation[dat1,FitRegularization->{\"Curvature\",0.01}];
dat2=Table[{x,y,Sin[x y]},{x,0,2,0.2},{y,0,2,0.2}]//Flatten[#,1]&;
gr2=RegularisedInterpolation[dat2
,FitRegularization->{\"Curvature\",0.01}];";

 RegularisedListInterpolation::usage="Works
like ListInterpolation but also inherits options from Fit including
FitRegularization.
 Example: \n dat1=Table[Sin[x ],{x,0,2,0.2}];
       gr1=RegularisedListInterpolation[dat1,FitRegularization->{\"Curvature\",0.01}];
dat2=Table[Sin[x y],{x,0,2,0.2},{y,0,2,0.2}];
gr2=RegularisedListInterpolation[dat2 ,FitRegularization->{\"Curvature\",0.01}]\
;";


Options[RegularisedInterpolation] = 
  Union[{PiecewiseSimplify -> False}, Options[Interpolation],
   Options[Fit]];
RegularisedInterpolation[dat_, opt : OptionsPattern[]] := 
 Module[{fspl, fb, sd, scpts, sk, var, dim, x, y, z, t, tt},
  fspl = Interpolation[dat, Method -> "Spline", 
    FilterRules[{opt}, Options[Interpolation]]];
  If[Length[FilterRules[{opt},Options[Fit]]]==0, Return[fspl]];
  fb = First[Cases[fspl, _BSplineFunction, \[Infinity]]];
  {sd, scpts, sk} = fb /@ {"Degree", "ControlPoints", "Knots"};
  dim = Last@Dimensions[If[Depth[dat] > 3, 
      Flatten /@ dat, dat]] - 1;
  var = Which[dim == 1, {x},
    dim == 2, {x, y}, dim == 3, {x, y, z}, dim == 4, {x, y, z, t}];
  tt = Compile[var // Evaluate, 
    Fit[If[Depth[dat] > 3, Flatten /@ dat, dat], 
       Flatten[Outer[Times, 
          Sequence @@ 
           MapThread[
            Table[BSplineBasis[{#1, #2}, k - 1, #3], {k, #4}] &, {sd, 
             sk, var // Evaluate, Dimensions[scpts]}]]] // Evaluate, 
       var // Evaluate,
       FilterRules[{opt}, Options[Fit]]
       ] // Chop // Evaluate];
  RegularisedInterpolatingFunction[<|
    "dimension" -> dim,
    "method" -> OptionValue[ FitRegularization],
    "order" -> OptionValue[InterpolationOrder],
    "domain" -> Most@( MinMax /@ Transpose[dat]), 
    "code" -> If[OptionValue[PiecewiseSimplify],
       tt = PiecewiseExpand /@ tt @@ var // Simplify;
       Compile[var // Evaluate, tt// Evaluate, 
        CompilationTarget -> "C"],
       tt] // Evaluate|>]
  ]

RegularisedListInterpolation[dat_, opt : OptionsPattern[]] :=RegularisedInterpolation[ 
   Flatten[MapIndexed[Flatten[{#2, #1}] &, dat, {TensorRank@dat}], 
    TensorRank[dat] - 1], opt];


RegularisedInterpolatingFunction /: 
  Format[b : RegularisedInterpolatingFunction[a_Association]] := 
  RawBoxes[BoxForm`ArrangeSummaryBox[
    "RegularisedInterpolatingFunction", "", 
    Graphics3D[Cone[], 
     ImageSize -> 
      20], {BoxForm`MakeSummaryItem[{"Dimensions: ", a["dimension"]}, 
      StandardForm], 
     BoxForm`MakeSummaryItem[{"Domain: ", MatrixForm@a["domain"]}, 
      StandardForm]}, {BoxForm`MakeSummaryItem[{"Regularisation \
method: ", MatrixForm@a["method"]}, StandardForm],
     BoxForm`MakeSummaryItem[{"Interpolation order: ", 
       a["order"] // Shallow}, StandardForm],
     BoxForm`MakeSummaryItem[{"Compiled Code: ", a["code"]}, 
      StandardForm]}, StandardForm]];    



RegularisedInterpolatingFunction[a_Association][b__] := 
 a[b] /; (b == "dimension" || b == "order" || b == "domain" || 
    b == "method")
RegularisedInterpolatingFunction[a_Association][b__] := a["code"][b]

With this code we have 2 new functions:

 ?RegularisedInterpolation

  ?RegularisedListInterpolation

Concluding remark

The present implementation still remains quite sub optimal, mainly because the Fit does not scale well, see this question. It could be enhanced by adding resampling (of the knots wrt data points) and Automated choice of regularisation (via GCV).

Answer 2

The function CurveToBSplineFunction Does something closely related:

dat1 = Table[{i, Sin[i^2]}, {i, 0, 2 Pi, Pi/64}] // N;
dat1 = dat1 /. {x_, y_} :> {x, 
     y (1 + RandomVariate[NormalDistribution[0, 0.05]])};
pl1 = ListPlot[dat1, PlotStyle -> Red];
g = ResourceFunction["CurveToBSplineFunction"][dat1, 1];
Show[pl1, ParametricPlot[g[u], {u, 0, 1}]]

And by changing the control parameter:

g = ResourceFunction["CurveToBSplineFunction"][dat1, 20];
Show[pl1, ParametricPlot[g[u], {u, 0, 1}]]

Note that CurveToBSplineFunction has a few useful options:

Options[ResourceFunction["CurveToBSplineFunction"]]

(* {CurveClosed->False,FillGaps->True,AssignWeights->False,DuplicateNodes->False} *)

All credits goes to Alexey Popkov who wrote this function!