Discrete Thin Plat Splines
Work with Stephen Roberts
The thinplate spline for a general domain Ω, as formulated by Wahba [
2] and Duchon [
1], is the function, f, that minimises the functional
J_{α}(f) = 
1
n


n ∑
i=1

(f(x^{(i)})−y^{(i)})^{2} + α  ⌠ ⌡

Ω


∑
ν=2






( D^{ν} f(x) )^{2} dx, 
 (1) 
where ν = (ν
_{1},...,ν
_{d}) is a d
dimensional multiindex, ν = ∑
_{s=1}^{d} ν
_{s},
x is a predictor variable in
R^{d}, and
x^{(i)} and y
^{(i)} are the corresponding ith predictor and response data value (1 ≤ i ≤ n).
The smoothing parameter α controls the tradeoff between smoothness and
fit. In the limit α→ 0 the function f becomes an
interpolant. If α is large, f becomes very smooth but may not reflect
the data very well. The value of α may be automatically calculated using generalised cross validation [
2].
Radial basis functions are often used to represent f as they give an analytical solution to the minimiser of the functional in Equation (
1). The system of equations resulting from the radial basis functions is dense and its size is directly proportional to the number of data points. In applications with a large number of data points the system is very expensive to compute.
We propose a discrete thinplate spline method that uses polynomial basis function with local support defined on a finite element mesh. The resulting system of equations is sparse and its size depends only on the number of grid points in the finite element mesh. A preconditioned conjugate gradient (CG) method is used to solve the system of equations.
The smoothing problem from Equation (
1) can be approximated with finite elements
so that the discrete smoother f is a piecewise multilinear
function. In vector notation f will be of the form
The idea is to minimise J
_{α} over all f of this
form. The smoothing term is not defined for
piecewise multilinear functions, but the nonconforming
finite element principle can be used to introduce piecewise multilinear functions
u
= (
b^{T}g_{1},...,
b^{T}g_{d})
to represent the gradient of f.
The functions f and
u satisfy the relationship
 ⌠ ⌡

Ω

∇f(x) ·∇v(x) dx =  ⌠ ⌡

Ω

u(x) ·∇v(x) dx, 
 (2) 
for all piecewise multilinear function v.
This is equivalent to the relationship
L c = 
d ∑
s=1

G_{s} g_{s}, 
 (3) 
where L is a discrete approximation to the negative Laplace
operator and (G
_{1} ,..., G
_{d}) is a discrete approximation to the transpose of the gradient operator.
We now consider the minimiser of the functional
 


1
n


n ∑
i=1

(b(x^{(i)})^{T}c−y^{(i)})^{2} + α  ⌠ ⌡

Ω


d ∑
s=1

∇(b^{T} g_{s})·∇(b^{T} g_{s}) dx 
 
 


1
n


n ∑
i=1

(b(x^{(i)})^{T}c−y^{(i)})^{2} + α 
d ∑
s=1

g_{s}^{T} L g_{s} . 
  (4) 

Our smoothing problem consists of minimising this functional over
all vectors
c,
g_{1},...,
g_{d}, defined on the domain Ω
_{h}, subject to
the Constraint (
3).
In the 3D case the discrete minimisation problem (
4)
is equivalent to finding the minimum of
J_{α}(c,g_{1}, g_{2}, g_{3}) = c^{T} A c − 2d^{T}c + y^{T}y/n + α(g_{1}^{T} L g_{1} + g_{2}^{T} L g_{2} + g_{3}^{T} L g_{3}), 
 (5) 
subject to
L c − G_{1} g_{1} − G_{2} g_{2} − G_{3} g_{3} = 0. 
 (6) 
The matrices L, G
_{1}, G
_{2} and G
_{3} are
independent of the data points but the matrix
A = 
1
n


n ∑
i=1

b(x^{(i)})b(x^{(i)})^{T}, 

and vector
d = 
1
n


n ∑
i=1

b(x^{(i)})y^{(i)}, 

must be assembled by sweeping through the data points. The matrix
A is symmetric indefinite and sparse.
By using Lagrange multipliers Minimisation Problem (
5) may be rewritten as the solution of the linear system:









= 



, 
 (7) 
where
w is a Lagrange multiplier associated with Constraint
(
6). We have assumed zero Dirichlet boundary conditions.
One way to solve
Equation (
7) is to eliminate all of the variables except
g_{1},
g_{2} and
g_{3}, which gives








= 



, 

where Z = L
^{−1}A L
^{−1} and
c = L
^{−1}(G
_{1} g_{1} + G
_{2} g_{2} + G
_{3} g_{3} −
h_{5} ).
We now briefly present some examples of data fitting in 3D.
Isosurface plot of the sphere on a grid with 189 and 68705 nodes.
The first data set contained 10
^{6} points generated randomly on a sphere with centre (0.5, 0.5, 0.5) and radius 1/3. The data points were assigned a value of 1. The smoothing parameter α was set to 10
^{−3}.
Figure
1 shows an Isosurface plot of the discrete thin plate spline with 189 nodes and 68705 nodes respectively.
Isosurface representation of the original data set for the semisphere example and approximation on a grid with 68705 nodes.
The next data set is designed to represent the case where there is missing data. Specifically, we took the data from the example described above and removed all of those points with xcoordinate less than 0.5. See Figure
2. The data set contained approximately 1/2 ×10
^{6} points. Note that the spline attempted to fill in the missing data. The value for α was set to 10
^{−3}.
Isosurface representation of the original data set for the semisphere example with an additional layer of data points placed in the interior and the approximation on a grid with 68705 nodes.
In the final example an inner layer of data points were added to force the discrete thin plate spline approximation to more closely follow the shape of the semisphere. Along with the data points from the previous example, additional data points sitting on the semisphere with centre (0.5, 0.5, 0.5) and radius = 1/6 were added and given a value of 0. See Figure
3. There are approximately 2×(1/2 ×10
^{6}) data points. The smoothing parameter, α, was reduced to 10
^{−7}.
References
 [1]

J. Duchon, Splines minimizing rotationinvariant., in Lecture Notes
in Math, vol. 571, SpringerVerlag, 1977, pp. 85100.
 [2]

G. Wahba, Spline models for observational data, vol. 59 of CBMSNSF
Regional Conference Series in Applied Mathematics, Society for Industrial and
Applied Mathematics (SIAM), Philadelphia, PA, 1990.
File translated from
T_{E}X
by
T_{T}H,
version 3.67.
On 12 Oct 2005, 23:18.