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. A special case is the thin plate spline.

Definition

A polyharmonic spline is a linear combination of polyharmonic radial basis functions (RBFs) (denoted by $\phi$) plus an optional linear term:

$f(\mathbf{x}) \, = \, \sum_{i=1}^N w_i \, \phi(|\mathbf{x} - \mathbf{c}_i|) + \mathbf{v}^T \, \begin{bmatrix} 1 \\ \mathbf{x} \end{bmatrix}$

(1)

where

File:Polyharmonic-splines-basic-functions.png

Polyharmonic basis functions

• $\mathbf{x} = [x_1, x_2, \cdots, x_{d}]^T$ is a real-valued vector of $d$ independent variables,
• $\mathbf{c}_i = [c_{i,1}, c_{i,2}, \cdots, c_{i,d}]^T$ are $N$ vectors of the same size as $\mathbf{x}$ (often called centers) that the curve or surface must interpolate
• $\mathbf{w} = [w_1, w_2, \cdots, w_N]^T$ are the $N$ weights of the RBFs.
• $\mathbf{v} = [v_1, v_2, \cdots, v_{d+1}]^T$ are the $d+1$ weights of the polynomial.
• The linear polynomial with the weighting factors $\mathbf{v}$ improves the interpolation close to the "boundary" and especially the extrapolation "outside" of the centers $\mathbf{c}_i$. If this is not desired, this term can also be removed (see also figure below).

The polyharmonic RBFs are of the form:

$\begin{matrix} \phi(r) = \begin{cases} r^k & \mbox{with } k=1,3,5,\dots, \\ r^k \ln(r) & \mbox{with } k=2,4,6,\dots \end{cases} \\[5mm] r = |\mathbf{x} - \mathbf{c}_i| = \sqrt{ (\mathbf{x} - \mathbf{c}_i)^T \, (\mathbf{x} - \mathbf{c}_i) } \end{matrix}$

Other values of the exponent $k$ are not useful (such as $\phi(r) = r^2$), because a solution of the interpolation problem might not exist. To avoid problems at $r=0$ (since $\log(0) = -\infty$), the polyharmonic RBFs with the natural logarithm might be implemented as:

$\phi(r) = \begin{cases} r^{k-1} \ln(r^r) & \mbox{for } r < 1 \\ r^k \ln(r) & \mbox{for } r \ge 1 \end{cases}$

The weights $w_i$ and $v_j$ are determined such that the function interpolates $N$ given points $(\mathbf{c}_i, f_i)$ (for $i=1,2,\dots,N$) and fulfills the $d+1$ orthogonality conditions:

$\sum_{i=1}^N w_i=0, \;\; \sum_{i=1}^N w_i \, c_{j,i}=0 \;\;\; (\text{for} \ j=1,2,...,d)$

To compute the weights, a symmetric linear system of equations has to be solved:

$\begin{bmatrix} \mathbf{A} & \mathbf{V}^T \\ \mathbf{V} & \mathbf{0} \end{bmatrix} \; \begin{bmatrix} \mathbf{w} \\ \mathbf{v} \end{bmatrix} \; = \; \begin{bmatrix} \mathbf{f} \\ \mathbf{0} \end{bmatrix}\;\;\;\;$

(2)

where

$A_{i,j} = \phi(|\mathbf{c}_i - \mathbf{c}_j|), \;\;\; \mathbf{V} = \begin{bmatrix} 1 & 1 & \cdots & 1 \\ \mathbf{c}_1 & \mathbf{c}_2 & \cdots & \mathbf{c}_{N} \end{bmatrix}, \;\;\; \mathbf{f} = [f_1, f_2, \cdots, f_N]^T$

Under very mild conditions (essentially, that at least $d+1$ points are not in a subspace; e.g. for $d=2$ that at least 3 points are not on a straight line), the matrix of the linear system of equations is nonsingular and therefore a unique solution of the equation system exists. The computed weights allow evaluation of the spline for any $\mathbf{x}\in\mathbb{R}^d$ using equation (1).

Many practical details to implement and use polyharmonic splines are given 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 $\Delta^m f = 0$ for any natural number $m$. For example, the biharmonic equation is $\Delta^2 f = 0$ and the triharmonic equation is $\Delta^3 f = 0$. 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. Applying the 2D Laplace operator ($\Delta = \partial_{xx} + \partial_{yy}$) to the thin plate radial basis function $f_{tp}(x,y) = (x^2+y^2) \log \sqrt{x^2+y^2}$ either by hand or using a computer algebra system shows that $\Delta f_{tp} = 4 + 4\log r$. Applying the Laplace operator to $\Delta f_{tp}$ (this is $\Delta^2 f_{tp}$) yields 0. But 0 is not exactly correct. To see this, replace $r^2=x^2+y^2$ with $\rho^2 = x^2+y^2+h^2$ (where $h$ is some small number tending to 0). The Laplace operator applied to $4 \log \rho$ yields $\Delta^2 f_{tp} = 8h^2 / \rho^4$. For $(x,y)=(0,0)$, the right hand side of this equation approaches infinity as $h$ approaches 0. For any other $(x,y)$, the right hand side approaches 0 as $h$ approaches 0. This indicates that the right hand side is a Dirac delta function. A computer algebra system will show that

$\lim_{h \to 0}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} 8h^2/(x^2+y^2+h^2)^2 \operatorname{d}\!x \operatorname{d}\!y = 8\pi.$

So the thin plate radial basis function is a solution of the equation $\Delta^2 f_{tp} = 8\pi\delta(x,y)$.

Applying the 3D Laplacian ($\Delta = \partial_{xx} + \partial_{yy} + \partial_{zz}$) to the biharmonic RBF $f_{bi}(x,y,z)=\sqrt{x^2+y^2+z^2}$ yields $\Delta f_{bi} = 2/r$ and applying the 3D $\Delta^2$ operator to the triharmonic RBF $f_{tri}(x,y,z) = (x^2+y^2+z^2)^{3/2}$ yields $\Delta^2 f_{tri} = 24/r$. Letting $\rho^2 = x^2+y^2+z^2+h^2$ and computing $\Delta(1/\rho) = -3h^2 / \rho^5$ again indicates that the right hand side of the PDEs for the biharmonic and triharmonic RBFs are Dirac delta functions. Since

$\lim_{h \to 0}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}-3h^2/(x^2+y^2+z^2+h^2)^{5/2}\operatorname{d}\!x\operatorname{d}\!y\operatorname{d}\!z = -4\pi,$

the exact PDEs satisfied by the biharmonic and triharmonic RBFs are $\Delta^2 f_{bi} = -8\pi\delta(x,y,z)$ and $\Delta^3 f_{tri} = -96\pi\delta(x,y,z)$.

Polyharmonic smoothing splines

Polyharmonic splines solve the minimization problem

$\sum_{i=1}^N (f(\mathbf{c}_i) - f_i)^2 + \lambda \int\limits_{\mathbb{R}^d} |\nabla^m f|^2 \operatorname{d}\!\mathbf{x}$

(3)

where $\nabla^m f$ is the vector of all $m$th order partial derivatives of $f$. For example, in 2D $\nabla^1 f = (f_x\ f_y)$ and $\nabla^2 f = (f_{xx} \ f_{xy} \ f_{yx} \ f_{yy})$ and in 3D $\nabla^2 f = (f_{xx} \ f_{xy} \ f_{xz} \ f_{yx} \ f_{yy} \ f_{yz} \ f_{zx} \ f_{zy} \ f_{zz})$. In 2D $|\nabla^2 f|^2 = f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2$, making the integral the simplified thin plate energy functional.

To show that polyharmonic splines solve equation (3), the fitting term must be transformed into an integral using the definition of the Dirac delta function:

$\sum_{i=1}^N (f(\mathbf{c}_i) - f_i)^2 = \int\limits_{\mathbb{R}^d}\sum_{i=1}^N (f(\mathbf{x}) - f_i)^2 \delta(\mathbf{x} - \mathbf{c}_i) \operatorname{d}\!\mathbf{x}.$

So equation (3) can be written as the functional

$J[f] = \int\limits_{\mathbb{R}^d} F(\mathbf{x},f, \partial^{\alpha_1}f, \partial^{\alpha_2}f, \dots \partial^{\alpha_n}f ) \operatorname{d}\!\mathbf{x} = \int\limits_{\mathbb{R}^d} \left[ \sum_{i=1}^N (f(\mathbf{x}) - f_i)^2 \delta(\mathbf{x} - \mathbf{c}_i) + \lambda |\nabla^m f|^2 \right] \operatorname{d}\!\mathbf{x}.$

where $\alpha_i$ is a multi-index that ranges over all partial derivatives of order $m$ for $\mathbb{R}^d$. In order to apply the Euler-Lagrange equation for a single function of multiple variables and higher order derivatives, the quantities

${\partial F \over\partial f} = 2\sum_{i=1}^N (f(\mathbf{x}) - f_i) \delta(\mathbf{x} - x_i)$

and

$\sum_{i=1}^n \partial^{\alpha_i} {\partial F \over\partial (\partial^{\alpha_i}f)} = 2\lambda \Delta^m f$

are needed. Inserting these quantities into the E-L equation shows that

$\left( \sum_{i=1}^N (f(\mathbf{x}) - f_i) \delta(\mathbf{x} - \mathbf{c}_i) \right) + (-1)^{m} \lambda \Delta^m f = 0.$

(4)

Let $f$ be a polyharmonic spline as defined by equation (1). The following calculations will show that $f$ satisfies (4). Applying the $\Delta^m$ operator to equation (1) yields

$\Delta^m f = \sum_{i=1}^N w_i C_{m,d} \delta(\mathbf{x} - \mathbf{c}_i)$

where $C_{2,2} = 8\pi,$ $C_{2,3}=-8\pi$, and $C_{3,3}=-96\pi$. So

$\sum_{i=1}^N \delta(\mathbf{x} - \mathbf{c}_i) (f(\mathbf{x}) - f_i + (-1)^m\lambda C_{m,d} w_i) = 0$

(5)

is equivalent to (4). Setting

$f(\mathbf{c}_j) - f_j + (-1)^m\lambda C_{m,d} w_j = 0 \quad \text{for} \ j=1,2,\dots,N$

(6)

(which implies interpolation if $\lambda=0$) solves equation (5) and hence solves (4). Combining the definition of $f$ in equation (1) with equation (6) results in almost the same linear system as equation (2) except that the matrix $A$ is replaced with $A + (-1)^m C_{m,d}\lambda I$ where $I$ is the $N\times N$ identity matrix. For example, for the 3D triharmonic RBFs, $A$ is replaced with $A + 288\pi\lambda I$.

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 phi = exp(-r²) which gives a good interpolation as well. Finally, the figure includes also the non-polyharmonic spline phi = r² 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).

File:Polyharmonic-splines-example1.png

Interpolation with different polyharmonic splines that shall pass the 4 predefined points marked by a circle (the interpolation with phi = r² is not useful, since the linear equation system of the interpolation problem has no solution; it is solved in a least squares sense, but then does not pass the centers)

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 = r² is no longer included). Since phi = (scale*r)^k = (scale^k)*r^k, the factor (scale^k) 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 phi = exp(-k*r²) with k=1, the interpolation is no longer reasonable and it would be necessary to adapt k.

File:Polyharmonic-splines-example1-scale100.png

The same interpolation as in the first figure, but the points to be interpolated are scaled by 100

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 = r² 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.

File:Polyharmonic-splines-example1-no-polynomial.png

The same interpolation as in the first figure, but without the polynomial term

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 $e^{-k\cdot r^2}$ needs to be tuned, so that $k$ is selected according to the underlying grid of the independent variables. If this grid is non-uniform, a proper selection of $k$ to achieve a good interpolation result is difficult or impossible.

Main disadvantages are:

• To determine the weights, a linear system of equations must be solved, which is non-sparse. The solution of a non-sparse linear system becomes no longer practical if the dimension n is larger as about 1000 (since the storage requirements are O(n²) and the number of operations to solve the linear system is O(n³). For example n=10000 requires about 100 Mbyte of storage and 1000 Gflops of operations).

• To perform the interpolation of M data points requires operations in the order of O(M*N). In many applications, like image processing, M is much larger than N, and if both numbers are large, this is no longer practical.

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 O(log(N)) instead of O(N).

See also

• Inverse distance weighting
• Radial basis function
• Subdivision surface (emerging alternative to spline-based surfaces)
• Spline