Polyharmonic spline

From COSSAN Wiki
Jump to: navigation, search

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}


  • \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:

  \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) }

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:

  \mathbf{A} & \mathbf{V}^T \\
  \mathbf{V} & \mathbf{0} \end{bmatrix}
  \mathbf{w} \\
\end{bmatrix} \; = \; 
  \mathbf{y} \\


  A_{i,j} =  \phi(||\mathbf{c}_i - \mathbf{c}_j||), \;\;\;
  \mathbf{V} =  
         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}.

See Also


  1. Duchon, J., Splines minimizing rotation-invariant seminorms in Sobolev spaces, Constructive Theory of Functions of Several Variables, Lecture Notes in Mathematics 571 (W. Schempp, K. Zeller, eds.), Springer-Verlag, Berlin, 1977, pp. 85-100