Polyharmonic spline


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 and natural cubic splines in one dimension.

Definition

A polyharmonic spline is a linear combination of polyharmonic radial basis functions denoted by plus a polynomial term:
where
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:
Other values of the exponent are not useful, because a solution of the interpolation problem might not exist. To avoid problems at , the polyharmonic RBFs with the natural logarithm might be implemented as:
The weights and are determined such that the function interpolates given points and fulfills the orthogonality conditions
All together, these constraints are equivalent to the symmetric linear system of equations
where
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 the Cholesky decomposition after a suitable transformation. 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. In Iske polyharmonic splines are treated as special cases of other multiresolution methods in scattered data modelling.

Reason for the name "polyharmonic"

A polyharmonic equation is a partial differential equation of the form for any natural number. For example, the biharmonic equation is and the triharmonic equation is. All the polyharmonic radial basis functions are solutions of a polyharmonic equation. For example, the thin plate radial basis function is a solution of the modified 2-dimensional biharmonic equation. Applying the 2D Laplace operator to the thin plate radial basis function either by hand or using a computer algebra system shows that. Applying the Laplace operator to yields 0. But 0 is not exactly correct. To see this, replace with . 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
So the thin plate radial basis function is a solution of the equation.
Applying the 3D Laplacian to the biharmonic RBF yields and applying the 3D operator to the triharmonic RBF yields. Letting and computing again indicates that the right hand side of the PDEs for the biharmonic and triharmonic RBFs are Dirac delta functions. Since
the exact PDEs satisfied by the biharmonic and triharmonic RBFs are and.

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 and in 3D. 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:
So equation can be written as the functional
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
and
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.
Let be a polyharmonic spline as defined by equation. The following calculations will show that satisfies. Applying the operator to equation yields
where and So is equivalent to
The only possible solution to for all test functions is
. 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 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 . For a given center
On a line for arbitrary point and unit vector
Dividing both numerator and denominator of this by shows that a quantity independent of the center So on the given line,
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 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 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
Integration by parts shows that
For example, let and Then
Integrating the first term of this by parts once yields
since vanishes at infinity. Integrating by parts again results in
So integrating by parts twice for each term of yields
Since shows that
So if and
Now the origin of the constraints 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
The equations
and
are equivalent to the two systems of linear equations and Since is invertible, the first system is equivalent to So the first system implies the second system is equivalent to 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 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 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 is reasonable. The figure also includes the radial basis functions phi = exp 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 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. Since phi = 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. 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.

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:
Recently, methods have been developed to overcome the aforementioned difficulties. For example Beatson et al. present a method to interpolate polyharmonic splines at one point in 3 dimensions in operations instead of operations.