# Polyharmonic spline

Polyharmonic splines, introduced by [1] ,are very useful for interpolation of scattered input-output data in many dimensions. Polyharmonic splines are defined as a linear combination of basis functions $\phi(\cdot)$ that depend only on distances together with a polynomial (in the following equation, only a linear polynomial is considered, please refer to Response_surface for other possible polynomials):

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

where

• $\mathbf{x} = [x_1, x_2, \cdots, x_{nx}]^T$ is a real-valued vector of nx independent variables,
• $\mathbf{c}_i = [c_{1,i}, c_{2,i}, \cdots, c_{nx,i}]^T$ are N vectors of the same size as $\mathbf{x}$ (often called centers).
• $\mathbf{w} = [w_1, w_2, \cdots, w_N]^T$ are the N weights of the basis functions.
• $\mathbf{v} = [v_1, v_2, \cdots, v_{nx+1}]^T$ are the nx+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$.

The basis functions of polyharmonic splines are functions 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||_2 = \sqrt{ (\mathbf{x} - \mathbf{c}_i)^T \, (\mathbf{x} - \mathbf{c}_i) } \end{matrix}$

Basis functions of polyharmonic spline.

## Calibration of a Polyharmonic spline

The weights $w_i$ and $v_j$ are determined such that the function passes through $N$ given points (called centers) $(\mathbf{c}_i, y_i)$ (i=1,2,...,N) and fulfill the $nx+1$ orthogonality conditions:

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

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{y} \\ \mathbf{0} \end{bmatrix}\;\;\;\;$

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{y} = [y_1, y_2, \cdots, y_N]^T$

Under very mild conditions (essentially, that at least ndim+1 points are not in a subspace; e.g. for ndim=2 that at least 3 points are not on a straight line), the system matrix of the linear system of equations is nonsingular and therefore a unique solution of the equation system exists.

Once the weights are determined, interpolation requires to just evaluate the top most formula for the provided $\mathbf{x}$.