Update
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_] =
basis.LinearSolve[
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_] =
basis.LinearSolve[
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]
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
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"]
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
Other methods for the choice of hyper parameters exist, see e.g.
this page
Update:
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_] :=
Module[{},
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"];
Show[pl1,pl2]
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 !
GaussianFilterwith furtherInterpolation/ListInterpolationwithMethod -> "Spline"be helpful? $\endgroup$GaussianFilter(or any other smoothing algorithm; perhaps a wavelet transform) + Spline interpolation would work as well. $\endgroup$