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


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

enter image description here

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.

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


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]

enter image description here

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

enter image description here

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];

enter image description here

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

enter image description here

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


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?

  1. 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?

  1. 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:

enter image description here

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?


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.

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

enter image description here

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;
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] &

enter image description here

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

  • 6
    $\begingroup$ I apologize for having been too harsh. I did not mean to be. I admit that I lost contenance over the frustration. Still, somebody who says "dog" all the time but actually means to say "duck" must not be surprised that their audience is confused. I am pretty sure that the term interpolation is always used in the same way in all discliplines (if correctly used): It denotes methods that reproduce the input data exactly. In contrast, for a fit ones is satisfied with reproducing the input approximately and typically focuses on simplicity or regularity of the model. $\endgroup$ – Henrik Schumacher Apr 4 '20 at 13:02
  • 2
    $\begingroup$ Sometimes, interpolation is used for fitting, but Runge's phenomenon tells us that this is actually not a good idea. So in some sense, interpolation mathods can be considered a subset of fitting methods. But the other inclusion is definitely not true: A fit that does not reproduce the input is by definition not an interpolation. $\endgroup$ – Henrik Schumacher Apr 4 '20 at 13:04
  • 1
    $\begingroup$ @HenrikSchumacher in these troubled times I suggest youtube.com/watch?v=lqnxP_2zyDc $\endgroup$ – chris Apr 5 '20 at 12:13
  • 1
    $\begingroup$ I must say don't understand the focus on linguistic in upvoting comments. The relevant issue should be,' is this function useful or not?' not wether is it legitimate to give the function a name which still contains the word Interpolation, for which there is some logic: adding some options to the interpolation function so that derivatives remain as smooth as they should be. $\endgroup$ – chris Apr 7 '20 at 6:49
  • 1
    $\begingroup$ By "performance" do you mean "goodness-of-fit" or "speed of calculations" ? Or both? If "goodness-of-fit", how are you characterizing goodness-of-fit AND what value of goodness-of-fit is adequate for your needs? $\endgroup$ – JimB Apr 8 '20 at 16:52

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` *)


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]


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

enter image description here

which can be expanded as

enter image description here

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

  gr = RegularisedListInterpolation[dat, 
     FitRegularization -> {"Curvature", 10^-0.5}];

It should work in dimensions 1 to 4.

The full code now reads

ClearAll[RegularisedInterpolation, RegularisedListInterpolation,

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}];
dat2=Table[{x,y,Sin[x y]},{x,0,2,0.2},{y,0,2,0.2}]//Flatten[#,1]&;

like ListInterpolation but also inherits options from Fit including
 Example: \n dat1=Table[Sin[x ],{x,0,2,0.2}];
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],
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], 
          Sequence @@ 
            Table[BSplineBasis[{#1, #2}, k - 1, #3], {k, #4}] &, {sd, 
             sk, var // Evaluate, Dimensions[scpts]}]]] // Evaluate, 
       var // Evaluate,
       FilterRules[{opt}, Options[Fit]]
       ] // Chop // Evaluate];
    "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]] := 
    "RegularisedInterpolatingFunction", "", 
     ImageSize -> 
      20], {BoxForm`MakeSummaryItem[{"Dimensions: ", a["dimension"]}, 
     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:


enter image description here


enter image description here

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

  • 2
    $\begingroup$ I know you're not too interested in linguistics, but maybe RegularisedInterpolatingFunction would be consistent with Interpolation -> InterpolatingFunction? Actually, could the result of RegularisedInterpolation be an InterpolatingFunction? $\endgroup$ – Chris K Apr 7 '20 at 20:32
  • $\begingroup$ @ChrisK you are absolutely right (twice). I wish RegularisedInterpolation Could simply be an InterpolatingFunction so I would not have to write special wrapper. And I will rename RegularisedInterpolationFunction to RegularisedInterpolatingFunction $\endgroup$ – chris Apr 8 '20 at 7:12

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.