Polyharmonic spline explained
In applied mathematics, polyharmonic splines are used for function approximation and data interpolation. They are very useful for interpolating and fitting scattered data in many dimensions. Special cases include thin plate splines[1] [2] and natural cubic splines in one dimension.[3]
Definition
A polyharmonic spline is a linear combination of polyharmonic radial basis functions (RBFs) denoted by
plus a polynomial term:
where
} (
denotes matrix transpose, meaning
is a column vector) is a real-valued vector of
independent variables,
ci=[ci,1 ci,2 … ci,d]rm{T
} are
vectors of the same size as
(often called centers) that the curve or surface must interpolate,
} are the
weights of the RBFs,
} are the
weights of the polynomial.
The polynomial with the coefficients
improves fitting accuracy for polyharmonic smoothing splines and also improves extrapolation away from the centers
See figure below for comparison of splines with polynomial term and without polynomial term.
The polyharmonic RBFs are of the form:
\begin{align}
\varphi(r)&=\begin{cases}
rk&withk=1,3,5,\ldots,\\
rkln(r)&withk=2,4,6,\ldots
\end{cases}\\[5mm]
r&=|x-ci|=\sqrt{(x-
(x-ci)}.
\end{align}
Other values of the exponent
are not useful (such as
), because a solution of the interpolation problem might not exist. To avoid problems at
(since
), the polyharmonic RBFs with the natural logarithm might be implemented as:
\varphi(r)=\begin{cases}
rk-1ln(rr)&forr<1, (thisworksbecause00isdefined)\\
rkln(r)&forr\ge1.
\end{cases}
or, more simply adding a continuity extension in
\varphi(r)=\begin{cases}
0&forr<\epsilon, (forsomeverysmallvalueof\epsilon,e.g.ifusingfloatingpointnumbersindoubleprecisions,\epsilon=10-200)\\
rkln(r)&forr\ge\epsilon.
\end{cases}
The weights
and
are determined such that the function interpolates
given points
(for
) and fulfills the
orthogonality conditions
All together, these constraints are equivalent to the symmetric linear system of equations
where
Ai,j=\varphi(|ci-cj|),
B=\begin{bmatrix}
1&1& … &1\\
c1&c2& … &cN
\end{bmatrix}rm{T
}, \quad \mathbf = [f_1, f_2, \ldots, f_N]^.
In order for this system of equations to have a unique solution,
must be full rank.
is full rank for very mild conditions on the input data. For example, in two dimensions, three centers forming a non-degenerate triangle ensure that
is full rank, and in three dimensions, four centers forming a non-degenerate tetrahedron ensure that B is full rank. As explained later, the linear transformation resulting from the restriction of the domain of the linear transformation
to the
null space of
is positive definite. This means that if
is full rank, the system of equations always has a unique solution and it can be solved using a linear solver specialised for symmetric matrices. The computed weights allow evaluation of the spline for any
using equation . Many practical details of implementing and using polyharmonic splines are explained in Fasshauer.
[4] In Iske
[5] polyharmonic splines are treated as special cases of other multiresolution methods in scattered data modelling.
Discussion
The main advantage of polyharmonic spline interpolation is that usually very good interpolation results are obtained for scattered data without performing any "tuning", so automatic interpolation is feasible. This is not the case for other radial basis functions. For example, the Gaussian function
needs to be tuned, so that
is selected according to the underlying grid of the independent variables. If this grid is non-uniform, a proper selection of
to achieve a good interpolation result is difficult or impossible.
Main disadvantages are:
- To determine the weights, a dense linear system of equations must be solved. Solving a dense linear system becomes impractical if the dimension
is large, since the memory required is
and the number of operations required is
- Evaluating the computed polyharmonic spline function at
data points requires
operations. In many applications (image processing is an example),
is much larger than
and if both numbers are large, this is not practical.
Fast construction and evaluation methods
One straightforward approach to speeding up model construction and evaluation is to use a subset of
nearest interpolation nodes to build a local model every time we evaluate the spline. As a result, the total time needed for model construction and evaluation at
points changes from
to
.This can yield better timings if
is much less than
. Such an approach is advocated by some software libraries, the most notable being
scipy.interpolate.RBFInterpolator.The main drawback is that it introduces small discontinuities in the spline and requires problem-specific tuning: a proper choice of the neighbors count,
.Recently, methods have been developed to overcome the aforementioned difficulties without sacrificing main advantages of polyharmonic splines.
First, a bunch of methods for fast
evaluation were proposed:
- Beatson et al.[6] present a method to interpolate polyharmonic splines with
being a basis function at one point in 3 dimensions or less
- Cherrie et al. present a method to interpolate polyharmonic splines with
as a basis function at one point in 4 dimensions or less
Second, an accelerated model construction by applying an iterative solver to an ACBF-preconditioned linear system was proposed by Brown et al. This approach reduces running time from
to
,and further to
when combined with accelerated evaluation techniques.
The approaches above are often employed by commercial geospatial data analysis libraries and by some open source implementations (e.g. ALGLIB).Sometimes domain decomposition methods are used to improve asymptotic behavior,reducing memory requirements from
to
, thus making polyharmonic splines suitable for datasets with more than 1.000.000 points.
Reason for the name "polyharmonic"
A polyharmonic equation is a partial differential equation of the form
for any natural number
, where
is the
Laplace operator. For example, the
biharmonic equation is
and the triharmonic equation is
. All the polyharmonic radial basis functions are solutions of a polyharmonic equation (or more accurately, a modified polyharmonic equation with a
Dirac delta function on the right hand side instead of 0). For example, the thin plate radial basis function is a solution of the modified 2-dimensional biharmonic equation.
[7] Applying the 2D Laplace operator (
\Delta=\partialxx+\partialyy
) to the thin plate radial basis function
either by hand or using a
computer algebra system shows that
. Applying the Laplace operator to
(this is
) yields 0. But 0 is not exactly correct. To see this, replace
with
(where
is some small number tending to 0). The Laplace operator applied to
yields
. For
the right hand side of this equation approaches infinity as
approaches 0. For any other
, the right hand side approaches 0 as
approaches 0. This indicates that the right hand side is a Dirac delta function. A computer algebra system will show that
\limh
8h2/(x2+y2+h2)2dxdy=8\pi.
So the thin plate radial basis function is a solution of the equation
\Delta2ftp=8\pi\delta(x,y)
.
Applying the 3D Laplacian (
\Delta=\partialxx+\partialyy+\partialzz
) to the biharmonic RBF
yields
and applying the 3D
operator to the triharmonic RBF
ftri(x,y,z)=(x2+y2+z2)3/2
yields
. Letting
and computing
\Delta(1/\rho)=-3h2/\rho5
again indicates that the right hand side of the PDEs for the biharmonic and triharmonic RBFs are Dirac delta functions. Since
\limh
-3h2/(x2+y2+z2+h2)5/2dxdydz=-4\pi,
the exact PDEs satisfied by the biharmonic and triharmonic RBFs are
\Delta2fbi=-8\pi\delta(x,y,z)
and
\Delta3ftri=-96\pi\delta(x,y,z)
.
Polyharmonic smoothing splines
Polyharmonic splines minimize
where
is some box in
containing a neighborhood of all the centers,
is some positive constant, and
is the vector of all
th order partial derivatives of
For example, in 2D
and
\nabla2f=(fxx fxy fyx fyy)
and in 3D
\nabla2f=(fxx fxy fxz fyx fyy fyz fzx fzy fzz)
. In 2D
making the integral the simplified
thin plate energy functional.
To show that polyharmonic splines minimize equation, the fitting term must be transformed into an integral using the definition of the Dirac delta function:
}\sum_^N (f(\mathbf) - f_i)^2 \delta(\mathbf - \mathbf_i) \,d\mathbf.
So equation can be written as the functional
} F(\mathbf,f, \partial^f, \partial^f, \ldots, \partial^f) \,d\mathbf = \int_ \left[\sum_{i=1}^N (f(\mathbf{x}) - f_i)^2 \delta(\mathbf{x} - \mathbf{c}_i) + \lambda |\nabla^m f|^2 \right]\,d\mathbf.
where
is a
multi-index that ranges over all partial derivatives of order
for
In order to apply the
Euler–Lagrange equation for a single function of multiple variables and higher order derivatives, the quantities
{\partialF\over\partialf}=
(f(x)-fi)\delta(x-xi)
and
{\partialF\over\partial
f)}=2λ\Deltamf
are needed. Inserting these quantities into the E−L equation shows that
A weak solution
of satisfies
for all smooth test functions
that vanish outside of
A weak solution of equation will still minimize while getting rid of the delta function through integration.
[8] Let
be a polyharmonic spline as defined by equation . The following calculations will show that
satisfies . Applying the
operator to equation yields
\Deltamf=
wiCm,d\delta(x-ci)
where
and
So is equivalent to
The only possible solution to for all test functions
is
(which implies interpolation if
). Combining the definition of
in equation with equation results in almost the same linear system as equation except that the matrix
is replaced with
where
is the
identity matrix. For example, for the 3D triharmonic RBFs,
is replaced with
Explanation of additional constraints
In, the bottom half of the system of equations (
}\mathbf = 0) is given without explanation. The explanation first requires deriving a simplified form of
when
is all of
First, require that This ensures that all derivatives of order
and higher of
vanish at infinity. For example, let
and
and
be the triharmonic RBF. Then
\varphizzy=3y(x2+y2)/(x2+y2+z2)3/2
(considering
as a mapping from
to
). For a given center
On a line
for arbitrary point
and unit vector
\varphizzy(x-P)=
| 3(a | | 2t-P2)((a2+b2t-P+(a1+b1t-P | | 2+b | |
|
((a | | 1t-P+(a2+b2t-P+(a3+b3t-P3/2 | | 1+b | |
|
.
Dividing both numerator and denominator of this by
shows that
a quantity independent of the center
So on the given line,
\limt\toinftyfzyy(x)=\limt\toinfty
wi\varphizyy(x-ci)=
wi\right)3b2(b
+
/
+
+
3/2=0.
It is not quite enough to require that because in what follows it is necessary for
to vanish at infinity, where
and
are multi-indices such that
For triharmonic
wiuj\varphi\alpha(x-ci)\varphi\beta(x-dj)
(where
and
are the weights and centers of
) is always a sum of total degree 5 polynomials in
and
divided by the square root of a total degree 8 polynomial. Consider the behavior of these terms on the line
as
approaches infinity. The numerator is a degree 5 polynomial in
Dividing numerator and denominator by
leaves the degree 4 and 5 terms in the numerator and a function of
only in the denominator. A degree 5 term divided by
is a product of five
coordinates and
The
(and
) constraint makes this vanish everywhere on the line. A degree 4 term divided by
is either a product of four
coordinates and an
coordinate or a product of four
coordinates and a single
or
coordinate. The
constraint makes the first type of term vanish everywhere on the line. The additional constraints
will make the second type of term vanish.
Now define the inner product of two functions
defined as a linear combination of polyharmonic RBFs
with
and
as
\langlef,g\rangle=
(\nablamf) ⋅ (\nablamg)dx.
Integration by parts shows that
For example, let
and
Then
Integrating the first term of this by parts once yields
fxxgxxdxdy=
fxgxx
dy-
fxgxxxdxdy=-
fxgxxxdxdy
since
vanishes at infinity. Integrating by parts again results in
So integrating by parts twice for each term of yields
\langlef,g\rangle=
f(gxxxx+2gxxyy+gyyyy)dxdy=
f(\Delta2g)dxdy.
Since shows that
\begin{align}
\langlef,f\rangle&=(-1)m
f(x)
wi(-1)mCm,d
dx
=(-1)mCm,d
wif(ci)\\
&=(-1)mCm,d
wiwj\varphi(ci-cj)=(-1)mCm,dwrm{T
} A \mathbf. \end
So if and
Now the origin of the constraints
}\mathbf = 0 can be explained. Here
is a generalization of the
defined above to possibly include monomials up to degree
In other words,
where
is a column vector of all degree
monomials of the coordinates of
The top half of is equivalent to
So to obtain a smoothing spline, one should minimize the scalar field
defined by
F(w,v)=|Aw+Bv-f|2+λCwrm{T
} A \mathbf.
The equations
=2Ai*(Aw+Bv-f)+2λCAi*w=0 rm{for} i=1,2,\ldots,N
and
}_ (A\mathbf + B\mathbf - \mathbf)=0 \quad\textrm \ i=1,2,\ldots,d+1
(where
denotes row
of
) are equivalent to the two systems of linear equations
and
}(A\mathbf + B\mathbf - \mathbf) = 0. Since
is invertible, the first system is equivalent to
So the first system implies the second system is equivalent to
}\mathbf = 0. Just as in the previous smoothing spline coefficient derivation, the top half of becomes
This derivation of the polyharmonic smoothing spline equation system did not assume the constraints necessary to guarantee that But the constraints necessary to guarantee this, and are a subset of
}w=0 which is true for the critical point
of
So
is true for the
formed from the solution of the polyharmonic smoothing spline equation system. Because the integral is positive for all
the linear transformation resulting from the restriction of the domain of linear transformation
to
such that
must be positive definite. This fact enables transforming the polyharmonic smoothing spline equation system to a symmetric positive definite system of equations that can be solved twice as fast using the Cholesky decomposition.
Examples
The next figure shows the interpolation through four points (marked by "circles") using different types of polyharmonic splines. The "curvature" of the interpolated curves grows with the order of the spline and the extrapolation at the left boundary (x < 0) is reasonable. The figure also includes the radial basis functions φ = exp(−r2) which gives a good interpolation as well. Finally, the figure includes also the non-polyharmonic spline phi = r2 to demonstrate, that this radial basis function is not able to pass through the predefined points (the linear equation has no solution and is solved in a least squares sense).
The next figure shows the same interpolation as in the first figure, with the only exception that the points to be interpolated are scaled by a factor of 100 (and the case phi = r2 is no longer included). Since φ = (scale·r)k = (scalek)·rk, the factor (scalek) can be extracted from matrix A of the linear equation system and therefore the solution is not influenced by the scaling. This is different for the logarithmic form of the spline, although the scaling has not much influence. This analysis is reflected in the figure, where the interpolation shows not much differences. Note, for other radial basis functions, such as φ = exp(−kr2) with k = 1, the interpolation is no longer reasonable and it would be necessary to adapt k.
The next figure shows the same interpolation as in the first figure, with the only exception that the polynomial term of the function is not taken into account (and the case phi = r2 is no longer included). As can be seen from the figure, the extrapolation for x < 0 is no longer as "natural" as in the first figure for some of the basis functions. This indicates, that the polynomial term is useful if extrapolation occurs.
See also
External links
Computer Code
Notes and References
- R.L. Harder and R.N. Desmarais: Interpolation using surface splines. Journal of Aircraft, 1972, Issue 2, pp. 189−191
- J. Duchon: Splines minimizing rotation-invariant semi-norms in Sobolev spaces. Constructive Theory of Functions of Several Variables, W. Schempp and K. Zeller (eds), Springer, Berlin, pp. 85−100
- Book: Wendland . Holger . Scattered Data Approximation . limited . 2005 . Cambridge University Press . 0521843359 . 9.
- G.F. Fasshauer G.F.: Meshfree Approximation Methods with MATLAB. World Scientific Publishing Company, 2007, ISPN-10: 9812706348
- A. Iske: Multiresolution Methods in Scattered Data Modelling, Lecture Notes in Computational Science and Engineering, 2004, Vol. 37,, Springer-Verlag, Heidelberg.
- R.K. Beatson, M.J.D. Powell, and A.M. Tan: Fast evaluation of polyharmonic splines in three dimensions. IMA Journal of Numerical Analysis, 2007, 27, pp. 427–450.
- Some algorithms for thin plate spline interpolation to functions of two variables . Powell. M. J. D.. 1993. Cambridge University Dept. Of Applied Mathematics and Theoretical Physics Technical Report. https://web.archive.org/web/20160125162849/http://www.univie.ac.at/nuhag-php/bibtex/open_files/po94_M%20J%20D%20Powell%2003%2093.pdf . January 7, 2016. 2016-01-25 .
- Book: Evans, Lawrence. Partial Differential Equations. limited. American Mathematical Society. 1998. 0-8218-0772-2. Providence. 450−452.