# Adding Regularisation to Built-in Interpolation Function

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

ExperimentalNumericalFunction::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:

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.

• It was drawn to my attention 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.

• 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. Apr 4, 2020 at 13:02
• 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. Apr 4, 2020 at 13:04
• @HenrikSchumacher in these troubled times I suggest youtube.com/watch?v=lqnxP_2zyDc Apr 5, 2020 at 12:13
• 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. Apr 7, 2020 at 6:49
• 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?
– JimB
Apr 8, 2020 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 *)

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.

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 @@
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[BoxFormArrangeSummaryBox[
"RegularisedInterpolatingFunction", "",
Graphics3D[Cone[],
ImageSize ->
20], {BoxFormMakeSummaryItem[{"Dimensions: ", a["dimension"]},
StandardForm],
BoxFormMakeSummaryItem[{"Domain: ", MatrixForm@a["domain"]},
StandardForm]}, {BoxFormMakeSummaryItem[{"Regularisation \
method: ", MatrixForm@a["method"]}, StandardForm],
BoxFormMakeSummaryItem[{"Interpolation order: ",
a["order"] // Shallow}, StandardForm],
BoxFormMakeSummaryItem[{"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).

• 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? Apr 7, 2020 at 20:32
• @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 Apr 8, 2020 at 7:12

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!