↑ arvau

Wikipedia / Polyharmonic spline

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

Definition

A polyharmonic spline is a linear combination of polyharmonic radial basis functions plus a polynomial term:

where

  • independent variables,

  • (often called centers) that the curve or surface must interpolate,

  • weights of the RBFs,

  • weights of the polynomial.

See figure below for comparison of splines with polynomial term and without polynomial term.

The polyharmonic RBFs are of the form:

), the polyharmonic RBFs with the natural logarithm might be implemented as:

orthogonality conditions

All together, these constraints are equivalent to the symmetric linear system of equations

where

to the null space is full rank, the system of equations ( 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.

Discussion

The main advantage of polyharmonic spline interpolation to achieve a good interpolation result is difficult or impossible.

Main disadvantages are:

  • and if both numbers are large, this is not practical.

Fast construction and evaluation methods

. Such an approach is advocated by some software libraries, the most notable being scipy.interpolate.RBFInterpolator. Recently, methods have been developed to overcome the aforementioned difficulties without sacrificing main advantages of polyharmonic splines.

evaluation were proposed:

  • Beatson et al. being a basis function at one point in 3 dimensions or less

  • Cherrie et al. 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. 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, 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 is the Laplace operator. For example, the biharmonic equation. 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. either by hand or using a computer algebra system approaches 0. This indicates that the right hand side is a Dirac delta function. A computer algebra system will show that

.

again indicates that the right hand side of the PDEs for the biharmonic and triharmonic RBFs are Dirac delta functions. Since

.

Polyharmonic smoothing splines

Polyharmonic splines minimize

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

is a multi-index 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

of () satisfies

A weak solution of equation () will still minimize () while getting rid of the delta function through integration.

be a polyharmonic spline as defined by equation ( satisfies ( operator to equation () yields

So () is equivalent to

The only possible solution to ( is

in equation () with equation () results in almost the same linear system as equation (

Explanation of additional constraints

In (

So on the given line,

will make the second type of term vanish.

as

Integration by parts shows that

Then

Integrating the first term of this by parts once yields

So integrating by parts twice for each term of () yields

() shows that

The top half of ( defined by

The equations

and

Just as in the previous smoothing spline coefficient derivation, the top half of (

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 adaptk.

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

References

Computer Code