Radial basis functions: Theory¶
Interpolation¶
Let \(f\) be a function which has been sampled at distinct points \(\mathbf{x}_1, \mathbf{x}_2, \ldots\),
\[f_i = f(\mathbf{x}_i) \quad (i = 1, 2, \ldots)\]
and we wish to create an interpolant (or approximant) of \(f\) from the data \(f_i\), \(\mathbf{x}_i\). If \(\phi\) is a non-negative function defined on \([0, \infty)\) then
\[\phi_i(\mathbf{x}) = \ \phi\left(\left\| \mathbf{x} - \mathbf{x}_i \right\| \right)\]
is a radially-symmetric function centred on \(\mathbf{x}_i\). A radial basis function interpolant \(s\) for \(f\) is a linear combination of such \(\phi_i\), i.e.,
\[s(\mathbf{x}) = \sum_{i=1}^{n} \alpha_i \phi_i(\mathbf{x})\]
where the coefficients \(\alpha_i\) are determined by the requirement that \(s\) be an interpolant:
\[f_j = s(\mathbf{x}_j) = \sum_{i=1}^{n} \alpha_i \phi_i(\mathbf{x}_j)\]
which is a linear system which can be written in matrix form
(1)¶\[\mathbf{f} = \Phi \mathbf{a}\]
where \(\mathbf{f}\) is the vector of \(f_i\), \(\mathbf{a}\) is the vector of \(\alpha_i\), and \(\Phi\) is the matrix with entries \(\phi\left(\left\| \mathbf{x}_j - \mathbf{x}_i \right\| \right)\): one solves the linear system (1) to find the \(\alpha_i\).
This interpolation method was introduced by Hardy [2] for applications in geophysics, and has since become one of the most popular methods for multidimensional scattered data interpolation. The texts by Buhmann [1] or Wendland [4] should be consulted for a detailed treatment of the theory.
Approximation¶
In the case that the \(f_i\) are contaminated by noise, one typically seeks an approximation rather than an interpolation, and this can be achieved by penalising the size of the function \(s\) by finding
\[\min\left\{ \ \sum_{j=1}^{N}\left| s(\mathbf{x}_j) - f_j \right| + \ \lambda\left\| s \right\| \right\}\]
where \(\lambda\) is a smoothing parameter which balances the smoothness and accuracy of the approximation. It can be shown that the minimiser of this problem is the solution to the linear system
(2)¶\[\mathbf{f} = (\Phi + \lambda I)\mathbf{a}\]
where \(I\) is the identity matrix of the same order as \(\Phi\). Clearly the system (2) reduces to the interpolation system (1) as \(\lambda\rightarrow 0\).
See Wendland’s paper [3] for a fuller discussion and references on the practicalites of radial basis function approximation.
The polynomial part¶
The quality of the interpolation (or approximation) is often improved by combining a low-order polynomial interpolant with that of the radial basis functions; in other words one seeks an interpolant \(s\) of the form
\[s(\mathbf{x}) = \ \sum_{i=1}^{n} \alpha_i \phi_i(\mathbf{x}) + p(\mathbf{x}).\]
In this module, the polynomial part is implemented with MVPoly
polynomials.