Radial basis functions: Theory
==============================

Interpolation
-------------

Let :math:`f` be a function which has been sampled at distinct
points :math:`\mathbf{x}_1, \mathbf{x}_2, \ldots`,

 .. math::

     f_i = f(\mathbf{x}_i) \quad (i = 1, 2, \ldots)

and we wish to create an interpolant (or approximant) of :math:`f`
from the data :math:`f_i`, :math:`\mathbf{x}_i`.  If :math:`\phi`
is a non-negative function defined on :math:`[0, \infty)` then

 .. math::

     \phi_i(\mathbf{x}) = \
     \phi\left(\left\| \mathbf{x} - \mathbf{x}_i \right\| \right)

is a radially-symmetric function centred on :math:`\mathbf{x}_i`.  A
*radial basis function interpolant* :math:`s` for :math:`f` is a linear
combination of such :math:`\phi_i`, i.e.,

 .. math::

     s(\mathbf{x}) = \sum_{i=1}^{n} \alpha_i \phi_i(\mathbf{x})

where the coefficients :math:`\alpha_i` are determined by the requirement
that :math:`s` be an interpolant:

 .. math::

     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

 .. math::
     :label: interp

     \mathbf{f} = \Phi \mathbf{a}

where :math:`\mathbf{f}` is the vector of :math:`f_i`,
:math:`\mathbf{a}` is the vector of :math:`\alpha_i`, and :math:`\Phi` is
the matrix with entries
:math:`\phi\left(\left\| \mathbf{x}_j - \mathbf{x}_i \right\| \right)`:
one solves the linear system :eq:`interp` to find the :math:`\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 :math:`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
:math:`s` by finding

 .. math::

     \min\left\{ \
     \sum_{j=1}^{N}\left| s(\mathbf{x}_j) - f_j \right| + \
     \lambda\left\| s \right\|
     \right\}

where :math:`\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

 .. math::
     :label: approx

     \mathbf{f} = (\Phi + \lambda I)\mathbf{a}

where :math:`I` is the identity matrix of the same order as
:math:`\Phi`. Clearly the system :eq:`approx` reduces to
the interpolation system :eq:`interp` as
:math:`\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
:math:`s` of the form

 .. math::

     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 :class:`MVPoly`
polynomials.

References
----------

.. [1] Martin D. Buhmann,
       "Radial Basis Functions: Theory and Implementations",
       Cambridge University Press, 2003.


.. [2] R. L. Hardy,
       "Multiquadric equations of topography and other irregular surfaces",
       Journal of Geophysical Research, 76(8):1905--1915, 1971.

.. [3] Holger Wendland,
       "Computational Aspects of Radial Basis Function Approximation",
       in K. Jetter *et al.* (eds.)
       Topics in Multivariate Approximation and Interpolation,
       Elsevier B.V., 2005, 231--256.

.. [4] Holger Wendland,
       "Scattered Data Approximation"
       Cambridge University Press, 2004.