How to interpolate {{xi,yi},... ,{xn,yn}} with minimum bending energy?

Question

Thin Plate Splines are explained here, and a Mathematica implementation of that for real samples { { x1, y2, z1 }, { x2, y2, z2 }, ... , { xn, yn, zn } } is given here. As I understand it, the algorithm gives the function $f(x, y)$ that interpolates through the samples { { xi, yi, zi }, ... } and minimizes the bending energy $\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }(D[f[x,y],\{x,1\}]^2+D[f[x,y],\{x,2\}]^2)dx$ $dy$. Now suppose I have the following:

data = {{1.0, 4.2}, {2.0, 9.3}, {4.0, 7.1}, {5.0, 4.3}, {8.0, 
5.7}, {9.0, 10.2}, {10.0, 10.6}, {12.0, 7.2}};

How would I find the function $f(x)$ that interpolates the {xi,yi} data samples, and minimizes the bending energy $\int _{-\infty }^{\infty } (f''[x])^2 dx.$ ?

****** Edit ******

The Mathematica implementation I refer to above is the answer by J.M.

I give an improvement on that implementation below:

ThinPlateInterpolation::usage = "ThinPlateInterpolation[{{x1,y1,z1},{x2,y2,z2}..}] returns a CompiledFunction that gives a smooth interpolation between the points provided. The function returned is a function of two real numbers.";

ThinPlateInterpolation[points_?(And[Last[Dimensions@#] === 3, 3 <= First[Dimensions@#],MatrixQ[#, Element[#, Reals] &]] &)]:=
Module[{xs, ys, zs, xys, len, func, m, b, xx,yy},
{xs, ys, zs} = N@Transpose[points];
xys = Transpose[{xs, ys}];
len = Length[xs];
func = Compile[{{x0, _Real}, {y0, _Real},{xsys, _Real,2},{n, _Integer}},
Join[
 With[{abs = Norm[#]}, 
   If[abs == 0.0, 0.0, abs^2* Log[abs]]] & /@ (Table[{x0, y0},n] - xsys),
 {1.0, x0, y0}]
];
m = func[#1, #2, xys, len] & @@@ xys;
m = m~Join~{ConstantArray[1.0, len]
~Join~{0., 0., 0.}, 
xs~Join~{0., 0., 0.}, ys~Join~{0., 0., 0.}};
b=zs~Join~{0., 0., 0.};
Compile @@ {{{xx, _Real}, {yy, _Real}}, 
  Quiet[func[xx, yy, xys,len].
  LinearSolve[m, b], {LinearSolve::luc,CompiledFunction::cfsa}]}
 ]

Minimizing $E[f]$ requires $f$ to be harmonic, right? so it is basically a piece-wise-quadratic interpolating polynomial, no? — AccidentalFourierTransform, Jan 26 '20 at 01:13
@AccidentalFourierTransform, I don't think so. The implementation that does this for samples {{x1,y1,z1}, ..., {xn,yn,zn}} more or less finds a linear combination of logarithms! See the implementation of that which I added. Doesn't make sense to me, but that is what it does. — Ted Ersek, Jan 26 '20 at 01:59
Well, in two dimensions the logarithm is the fundamental solution to the Laplace equation, while in one dimension the solution is a parabola. So I am not surprised the implementation combines logarithms, and I still think that in your case the functions should be parabolas. But I might be wrong here... — AccidentalFourierTransform, Jan 26 '20 at 02:34
@AccidentalFourierTransform, Parabola's have zero bending energy, but the second derivative will be undefined where the function changes from one parabola to another. The goal is to find the interpolating function with second derivative defined at ALL real numbers that minimizes Integrate[ D[f[x],{x,2}]^2, {x,-Infinity,Infinity}]. — Ted Ersek, Jan 26 '20 at 20:55

score 2 · Answer 1 · answered Jan 28 '20 at 11:20

In the very simple case you have given, the expression for the required thin plate spline is very easily derived, and the elaborate machinery in my other answer is not necessary. Using formulae 28, 29, and 31 in Eberly's note, we have the following:

pts = {{1., 4.2}, {2., 9.3}, {4., 7.1}, {5., 4.3}, {8., 5.7}, {9., 10.2},
       {10., 10.6}, {12., 7.2}};

{xa, ya} = Transpose[pts];
lsf = LinearSolve[DistanceMatrix[xa, xa]^3/12];
nmat = PadLeft[{xa}, {2, Automatic}, 1];

bc = LinearSolve[nmat.lsf[Transpose[nmat]], nmat.lsf[ya]];
ac = lsf[ya - Transpose[nmat].bc];

tps[x_] = bc.{1, x} + ac.Abs[x - pts[[All, 1]]]^3/12;

Plot[tps[x], {x, 1, 12},
     Epilog -> {Directive[AbsolutePointSize[6], ColorData[97, 4]], Point[pts]}]

1D thin plate spline

How to interpolate {{xi,yi},... ,{xn,yn}} with minimum bending energy?

1 Answers1