# Thin-Plate-Spline understanding and solution

This is a migrated question from: Thin-Plate-Spline understanding and solution. If I need to delete one of the questions let me know. I was suggested to post it here as well.

As I understand it a Thin-Plate-Spline in 2D is an interpolant function that minimizes $$\left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} f(x_i,y_i) = z_i & i=1,\ldots,n \end{array} \right.$$

However according to the link above the minimizer of such problem is equivalent to solve the biharmonic equation

$$\Delta^2f = 0$$

Which I don't really understand. If the problem was unconstrained I agree the two are equivalent (by using the Euler Lagrange equations).

However in case of constrained interpolation and assuming R.K.H.S. I think I need to minimize the following actually

$$U(f,\vec\lambda) = \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy - \sum_{i=1}^n \lambda_i\left(z_i - L_i f\right)$$

where $$L_i$$ is the evaluation functional, which exists due to the R.K.H.S. hypothesis. However I don't know how to minimize such functional (and I assume it wouldn't just yield the biharmonic equation at the end).

Can anyone either give me a reference or if it is easy enough just show the derivation of the minimizer?

Update: Based on my reading Calculus of Variations - Gelfand & Fomin I've tried to formulate correctly my problem.

I am still assuming RKHS as setting (so I can use evaluation functional), although I see (including in one of the answer given) commonly weak solutions and distribution theory is used. Suppose I want to solve:

$$\left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} f(x_i,y_i) = z_i & i=1,\ldots,n \end{array} \right.$$

Observe that for $$i = 1, \ldots, n$$ we have can rewrite each constraint as

$$L_i f = z_i \text{for i = 1,\ldots,n}$$

where $$L_i$$ is the evaluation functional, which exists by my hypothesis of RKHS. Using Reisz representation theorem each constraint can be rewritten in theform

$$L_i f = \int_{\mathbb{R^2}} \omega_i(x,y) f(x,y)dxdy = z_i \text{for i = 1,\ldots,n}.$$ I can rewrite the constrained optimization problem as

$$\left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} \int_{\mathbb{R}^2} \omega_i(x,y)f(x,y)dxdy = z_i & i=1,\ldots,n \end{array} \right.$$

I am looking at this last formulation as a special case of the isoperimetric problem (see Calculus of Variations - Gelfand & Fomin) therefore using the lagrange multipliers therefore this constrained problem should be equivalent to solve

$$\int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2 + \sum_{i=1}^n \lambda_i\omega_i(x,y)f(x,y) \right)dxdy$$

If now I did my calculations using the Euler Equations I get the following PDE to solve

$$\Delta^2 f + \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i = 0 \hspace{10mm} (1)$$

which make sense to me now since in the distribution theory way those $$\omega_i$$ would be Dirac deltas probably.

Update:

I gave it more thought so I'll add an update (more for reference for people like me who are still reading on the subject). Equation $$(1)$$ can be rewritten as

$$\Delta^2 f = - \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i$$

And I observe that such equation is defined on $$\mathbb{R}^d$$ (so not a boundary value problem). I've been reading through fundamental solutions and Green function, from Evans - Partial Differential Equation. In my setting I believe the use of fundamental solution is more appropriate than the Green function.

To find the fundamental solution we can solve the equivalent following system

$$\left\{ \begin{array}{l} \Delta g = -\frac{1}{2}\sum_{i=1}^n\lambda_i \omega_i \\ \Delta f = g \end{array} \right.$$

So essentially we need to solve twice the laplace equation, non homogeneous. In my reference (Evans) how to solve the laplace equation is explained and I think a very similar observation can be done here (i.e. the fundamental solution is a function) $$\Phi(x) = \Phi(\left| x \right|)$$ (so radial) and solution of the non homogenous equation can be obtained as

$$f(x) = - \frac{1}{2} \int_{\mathbb{R}^d} \Phi(x - y) \left( \sum_{i=1}^n \lambda_i \omega_i(y) \right) dy = - \frac{1}{2} \sum_{i=1}^n \lambda_i\int_{\mathbb{R}^d} \Phi(x - y) \omega_i(y) dy. \hspace{10mm} (2)$$

Remember that each $$\omega_i$$ is a representation of the evaluation functional $$L_i$$ therefore I can write $$\int_{\mathbb{R}^d} \Phi(x - y)\omega_i(y) dy = \Phi(x - x_i)$$. This leads to write the solution $$f$$ as

$$f(x) = -\frac{1}{2} \sum_{i=1}^n \lambda_i \Phi(x - x_i).$$

Observe that for $$j = 1 ... n$$ we can write

$$f(x_j) = -\frac{1}{2} \sum_{i=1}^n \lambda_i \Phi(x_j - x_i) = \sum_{i=1}^n \gamma \Phi(x_j - x_i)$$

with $$\gamma_i = - \frac{\lambda_i}{2}$$ and this leads to a square linear system of the form

$$\begin{pmatrix} \Phi(0) & \Phi(x_2 - x_1) & \dots & \Phi(x_n - x_1) \\ \Phi(x_1 - x_2) & \Phi(0) & \dots & \Phi(x_n - x_2) \\ \vdots & \vdots & \ddots & \vdots \\ \Phi(x_1 - x_n) & \Phi(x_2 - x_n) & \dots & \Phi(0) \end{pmatrix} \begin{pmatrix} \gamma_1 \\ \gamma_2 \\ \vdots \\ \gamma_n \end{pmatrix} = \begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix}$$

At this point I am not very sure how to show the system is solvable. But I know it can be solved, and I think this is overall the theory of the TPS.

• Part of this is presumably related to the property of the biharmonic equation in the plane where isolated points have positive capacity, i.e. you can prescribe a boundary condition at points, unlike the Laplace equation (this property can be read off the fundamental solution: it's continuous at the origin). So you can minimize over the class of functions with prescribed values at the given set of points, and you'll get a minimizer in this class (and it will satisfy the EL equation on the complement of the set of points, analogously to what you derive with Lagrange multipliers formally). May 22, 2023 at 4:22
• Actually minimizing would perhaps require more care, though, or at least some consideration of how your functions behave at infinity. May 22, 2023 at 4:26
• What if I assume maybe $f \in W(\mathbb{R}^2)^{3,2} \cap C^3(\mathbb{R}^2)$? (so first 3 order derivatives in $L^2$ but also continuous?) Is it too restrictive? May 22, 2023 at 4:37
• I don't know, and maybe someone can answer more conclusively (it's clear you can minimize over some homogeneous space like $\dot{H}^2$, but you have to check whether you get unique solutions that way). I'd suggest thinking through the 1D Laplacian analogy: minimize $\int_{\mathbb{R}} (u')^2$ over $u$ with some finite number of points prescribed. You should get a piecewise linear function which has slope $0$ at the rightmost and leftmost (infinite) intervals as the only minimizer, so the energy form "forced" this nontrivial boundary condition at infinity. May 22, 2023 at 6:50
• I'll wait for a more conclusive answer. May 22, 2023 at 7:36