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.
