Does Mathematica have 2D smoothing spline interpolation built in? I requires an interpolation method with smooth first derivatives and cubic bivariate splines fulfill this nicely. In python I would use RectBivariateSpline or SmoothBivariateSpline.

A quick search only revealed this answer, which I guess could be adapted to 2D with some effort.

Here is some test data:

F[x_, y_] := 
 Sin[.5 y] Cos[.9 x]/Sec[0.1 x y] - 
  0.01 (x^2 + y^2) RiemannSiegelZ[1.5 Sqrt[x^2 + y^2]]
data = N[Flatten[
    Table[{x, y, F[x, y]}, {x, -RANGEX, RANGEX, 1}, {y, -RANGEY, 
      RANGEY, 1}], 1]];
(*add some noise*)
data[[All, 3]] = 
  data[[All, 3]] + 
   RandomVariate[NormalDistribution[0, 0.1], Length[data]];

PlotPointsAndSurface[points_, surface_, label_] := Module[{},
     PlotStyle -> {Directive[PointSize[0.01], Red], 
       Directive[PointSize[0.01], Green]}, PlotLabel -> label, 
     ImageSize -> Medium],
    Plot3D[surface, {x, -RANGEX, RANGEX}, {y, -RANGEY, RANGEY}, 
     PlotStyle -> Directive[Purple, Opacity[0.2]]]
PlotPointsAndSurface[data, F[x, y], "orig and noise"]

enter image description here

  • 1
    $\begingroup$ Can GaussianFilter with further Interpolation/ListInterpolation with Method -> "Spline" be helpful? $\endgroup$
    – ybeltukov
    Commented Oct 15, 2014 at 18:57
  • $\begingroup$ @ybeltukov Yes, GaussianFilter (or any other smoothing algorithm; perhaps a wavelet transform) + Spline interpolation would work as well. $\endgroup$
    – Ajasja
    Commented Oct 15, 2014 at 20:56
  • $\begingroup$ @RahulNarain Indeed, but the surface also goes through every point, which if there is noise is not desirable. $\endgroup$
    – Ajasja
    Commented Oct 16, 2014 at 11:11

3 Answers 3



Since version 12, this functionality in integrated in Mathematica via the Option FitRegularization

Following on @Ajasja's answer in the spirit of this answer one can in fact provide controlled smoothing va an explicit Tichonov like penalty as follows:

ff = Function[{x, y}, basis // Evaluate]; 
a = ff @@ # & /@ (Most /@ data);

so that

fit[x_, y_] = 
    Transpose[a]. a + 0 IdentityMatrix[Length[basis]], 
    Transpose[a].( Last /@ data )];
 pl0 = PlotPointsAndSurface[data, fit[x, y], "fit and data"];

reproduces exactly @Ajasja's fit, whereas e.g.

 fit[x_, y_] = 
    Transpose[a]. a + 10^1 IdentityMatrix[Length[basis]], 
    Transpose[a].( Last /@ data )];
 pl1 = PlotPointsAndSurface[data, fit[x, y], "fit and data"];

would correspond to a smoother solution.

  Show[pl1, pl0]

Mathematica graphics

Note the hyper parameter (here 10^1), which fixes the sought level of smoothness imposed onto the solution, by effectively correlating the coefficients of the basis expansion.

The main advantage is that one need not focus too much on the exact properties of the chosen basis.

For instance we could also use BSplineBasis

knots = Range[-RANGEX - 2, RANGEX + 2];
basis = Flatten@ Table[BSplineBasis[{3, knots}, i, x]
   BSplineBasis[{3, knots}, j,  y], {i, 0, 2 RANGEX}
    , {j, 0, 2 RANGEX}];

Then, as previously

ff = Function[{x, y}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);

Then one could use a penalty function based on second derivatives:

s = SparseArray[{{i_, i_} -> -1, {i_, j_} /; i - j == 1 -> 2,
     {i_, j_} /; i - j == 2 -> -1}, {17, 15}] // Transpose;
s1 = ArrayFlatten[TensorProduct[s, s]]; 
pen = Transpose[s1].s1; pen//ArrayPlot

Mathematica graphics

built so that s.( Range[17]*0 + 1) and s.Range[17] are both null (i.e. there is no penalty to have a constant or linear function of x and y.

Then, as previously

fit3[x_, y_] = 
  basis.LinearSolve[Transpose[a]. a + 10^1 pen, 
    Transpose[a].( Last /@ data )];
pl1 = PlotPointsAndSurface[data, fit3[x, y], "fit and data"]

Mathematica graphics

The main advantage of this second approach is that it is the penalty which sets smoothing, not the sampling of the basis function. Even if the conditioning of Transpose[a]. a is poor, the inverse will be well conditioned thanks to the regularisation terms pen.

Note that for the sake of efficiency and memory one could fill the a matrix using sparse matrices following this answer.

With[{xOrder = Ordering[Join[data[[All, 1]], knots]],
      yOrder = Ordering[Join[data[[All, 2]], knots]]},
 With[{xPar = xOrder[[# + 1 ;; #2 - 1]] & @@@ Partition[Ordering[xOrder, -Length[knots]], 2, 1],
       yPar = yOrder[[# + 1 ;; #2 - 1]] & @@@ Partition[Ordering[yOrder, -Length[knots]], 2, 1]},
  nonzero = Join @@ Outer[Intersection, Union @@@ Partition[xPar, 4, 1], 
                                        Union @@@ Partition[yPar, 4, 1], 1];]]

colIndex = Range[Length[basis]];

a2 = SparseArray[Join @@ MapThread[Thread[Thread[{#2, #3}] -> 
    Function[{x, y}, #] @@@ data[[#2, {1, 2}]]] &, {basis, nonzero,colIndex}]];  a == a2

(* True *)

The choice of optimal level of smoothing can be done via generalised cross validation, i.e. by choosing the penalty weight to correspond to the minimum of $$ \hat \lambda = {\rm min}_\lambda\left\{ \frac{||( \mathbf{1}- \tilde{\mathbf{a}}) \cdot {\mathbf{y}} ||^2}{ \left[{\rm trace}( \mathbf{1}- \tilde{\mathbf{a}}) \right]^2} \right\} \,. $$ having defined $$ \tilde{\mathbf{a}}(\lambda) =\mathbf{a} \cdot ({\mathbf{a}^{\rm T}} \cdot \mathbf{a} + \lambda\, \mathbf{s}^{\rm T}\cdot \mathbf{s})^{-1} \cdot {\mathbf{a}^{\rm T}} $$

Table[at = a.Inverse[Transpose[a]. a + 10^i pen].Transpose[a];
  {i, ((IdentityMatrix[289] - at).(Last /@ data) // #.# &)/
    Tr[IdentityMatrix[289] - at]^2}, {i, -3, 3, 1/2}] // ListLinePlot

Mathematica graphics

Other methods for the choice of hyper parameters exist, see e.g. this page


Note that if smoothing is not an issue, then in version 10 and above mathematica can deal directly with the data as demonstrated here

PlotPointsAndSurface2[points_, surface_, label_] := 
   Show[ListPlot3D[points, PlotLabel -> label, ImageSize -> Medium,
     PlotStyle -> Directive[Orange, Opacity[0.5]]], 
    Plot3D[surface, {x, -RX, RX}, {y, -RY, RY}, 
     PlotStyle -> Directive[Purple, Opacity[0.1]]]]];

 pl2 = PlotPointsAndSurface2[data, fit3[x, y], "fit and data"];

Mathematica graphics

as can be seen the regularised and un regularised surfaces look quite similar.

It would be great if mathematica allowed for adding a penalty to the built in function behind ListPlot3D, ListContourPlot or ListInterpolate !

  • $\begingroup$ Is there a way to evaluate a sufficient number of knots without reference to the range of the dataset? For example, imagine a dataset within the interval [-0.5,0.5] in both x and y and highly varying in this interval $\endgroup$ Commented Apr 22, 2020 at 19:16
  • $\begingroup$ 2@DanieleBinosi of course you can put the knots where you want. You need just to decide how you couple the coefficients for the penalty. Should it reflect this uneven sampling or not is up to you. $\endgroup$
    – chris
    Commented Apr 22, 2020 at 19:44
  • $\begingroup$ I am new into this spline issue, so I am a bit lost. Allow me just 1D for the moment. In your original example the data extends in [-8,8] and you have chosen knots = Range[-8- 2, 8+ 2]. Now if I have data say in [0,0.02] how do I $\endgroup$ Commented Apr 22, 2020 at 20:25
  • $\begingroup$ it depends :-) knots = Range[1,n]/n*0.025 -0.001 or something like that? you could have a look at this mathematica.stackexchange.com/a/34880/1089 and the link in the answer above $\endgroup$
    – chris
    Commented Apr 23, 2020 at 8:00
  • $\begingroup$ Thank you @chris. I am now trying to put everything in a module (still 1D). But I have problems with the following block: ff = Function[{x}, basis // Evaluate]; a = ff @@ # & /@ (Most /@ data). The matrix a does not get evaluated. It seems to be related to the scope of the x variable in basis = Flatten@ Table[BSplineBasis[{3, knots}, i, x] , {i, 0, 2RANGEX}] and ff. I guess a solution could be to fill the matrix a using sparse matrices but I cannot reduce your 2D example to the 1D case $\endgroup$ Commented Apr 23, 2020 at 8:36

Here is a Fourier Basis approach:


FourierBasis2D[{numx_, numy_}, {λx_, λy_}, x_, y_] := 
  N[With[{ωn = 2 π/λx, ωm = 
      2 π/λy},
      Table[ {Cos[ n ωn x] Cos[m  ωm y], 
        Cos[ n ωn x] Sin[  m  ωm y], 
        Sin[ n  ωn x] Cos[m ωm  y], 
        Sin[ n ωn x] Sin[ m  ωm y]}, {n, numx}, {m, 
FourierBasis2D[num_, λ_, x_, y_] := 
  FourierBasis2D[{num, num}, {λ, λ}, x, y];

Clear[basis, fit];
basis = FourierBasis2D[5, 20, x, y];
fit[x_, y_] = Fit[data, basis, {x, y}];
PlotPointsAndSurface[data, fit[x, y], "fit and data"]

Fit and data

Plot derivatives of original and interpolated function:

Clear[DFx, DFy, DfitX, DfitY]
DFx[x_, y_] = Simplify@D[F[x, y], x];
DFy[x_, y_] = Simplify@D[F[x, y], y];
DfitX[x_, y_] = Simplify@D[fit[x, y], x];
DfitY[x_, y_] = Simplify@D[fit[x, y], y];

Plot3D[{DFx[x, y], DfitX[x, y]}, {x, -RANGEX, RANGEX}, {y, -RANGEY, 
  RANGEY}, PlotStyle -> {Directive[Black, Opacity[0.5]], 
   Directive[Red, Opacity[0.5]]}]
Plot3D[{DFy[x, y], DfitY[x, y]}, {x, -RANGEX, RANGEX}, {y, -RANGEY, 
  RANGEY}, PlotStyle -> {Directive[Black, Opacity[0.5]], 
   Directive[Red, Opacity[0.5]]}]

It's not perfect, but it works...

  • 1
    $\begingroup$ Thanks for the formatting @Mr.W $\endgroup$
    – Ajasja
    Commented Oct 16, 2014 at 8:37

You can use ListPlot3D and its InterpolationOrder option. That has built in smoothing.

ListPlot3D[points, InterpolationOrder->3]

This will give you a cubic interpolation. Close to what you are looking for and built into Mma.

EDIT: Well, if the points are evenly spaced which they appear to be. But if not then no as this will not produce a scatter plot. Just realised the difference between ListPlot3D and ListPointPlot3D. I'm still a bit new here.


  • $\begingroup$ Hi! Perhaps the question or the example code was a bit misleading. I need to obtain a differentiable interpolated and smoothed function from the data. In the test case it should be as close as possible to the original. $\endgroup$
    – Ajasja
    Commented Oct 15, 2014 at 17:59
  • $\begingroup$ A duplicate question with code to solve this problem is mathematica.stackexchange.com/questions/11765/…. That gives the smoothest possible interpolation. The only problem is that the algorithm quickly gets slow as the number of points interpolated grows. $\endgroup$
    – Ted Ersek
    Commented Jun 6, 2019 at 0:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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