Community Hub0 Subscribers
Community hub
Radial basis function interpolation
View on Wikipediafrom Wikipedia
Radial basis function (RBF) interpolation is an advanced method in approximation theory for constructing high-order accurate interpolants of unstructured data, possibly in high-dimensional spaces. The interpolant takes the form of a weighted sum of radial basis functions.[1][2] RBF interpolation is a mesh-free method, meaning the nodes (points in the domain) need not lie on a structured grid, and does not require the formation of a mesh. It is often spectrally accurate[3] and stable for large numbers of nodes even in high dimensions.
Many interpolation methods can be used as the theoretical foundation of algorithms for approximating linear operators, and RBF interpolation is no exception. RBF interpolation has been used to approximate differential operators, integral operators, and surface differential operators.
Examples
[edit]Let and let be 15 equally spaced points on the interval . We will form where is a radial basis function, and choose such that ( interpolates at the chosen points). In matrix notation this can be written as
![{\displaystyle [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/738f7d23bb2d9642bab520020873cccbef49768d)
Choosing , the Gaussian, with a shape parameter of , we can then solve the matrix equation for the weights and plot the interpolant. Plotting the interpolating function below, we see that it is visually the same everywhere except near the left boundary (an example of Runge's phenomenon), where it is still a very close approximation. More precisely the maximum error is roughly at .
The function sampled at 15 uniform nodes between 0 and 1, interpolated using the Gaussian RBF with a shape parameter of .
The interpolation error, , for the plot to the left.
Motivation
[edit]The Mairhuber–Curtis theorem says that for any open set in with , and linearly independent functions on , there exists a set of points in the domain such that the interpolation matrix
This means that if one wishes to have a general interpolation algorithm, one must choose the basis functions to depend on the interpolation points. In 1971, Rolland Hardy developed a method of interpolating scattered data using interpolants of the form . This is interpolation using a basis of shifted multiquadric functions, now more commonly written as , and is the first instance of radial basis function interpolation.[5] It has been shown that the resulting interpolation matrix will always be non-singular. This does not violate the Mairhuber–Curtis theorem since the basis functions depend on the points of interpolation. Choosing a radial kernel such that the interpolation matrix is non-singular is exactly the definition of a strictly positive definite function. Such functions, including the Gaussian, inverse quadratic, and inverse multiquadric are often used as radial basis functions for this reason.[6]
Shape-parameter tuning
[edit] |
Many radial basis functions have a parameter that controls their relative flatness or peakedness. This parameter is usually represented by the symbol with the function becoming increasingly flat as . For example, Rolland Hardy used the formula for the multiquadric, however nowadays the formula is used instead. These formulas are equivalent up to a scale factor. This factor is inconsequential since the basis vectors have the same span and the interpolation weights will compensate. By convention, the basis function is scaled such that as seen in the plots of the Gaussian functions and the bump functions.
-
A Gaussian function for several choices of -
A plot of the scaled bump function with several choices of shape parameter
A consequence of this choice is that the interpolation matrix approaches the identity matrix as leading to stability when solving the matrix system. The resulting interpolant will in general be a poor approximation to the function since it will be near zero everywhere, except near the interpolation points where it will sharply peak – the so-called "bed-of-nails interpolant" (as seen in the plot to the right).
On the opposite side of the spectrum, the condition number of the interpolation matrix will diverge to infinity as leading to ill-conditioning of the system. In practice, one chooses a shape parameter so that the interpolation matrix is "on the edge of ill-conditioning" (eg. with a condition number of roughly for double-precision floating point).
There are sometimes other factors to consider when choosing a shape-parameter. For example the bump function has a compact support (it is zero everywhere except when ) leading to a sparse interpolation matrix.
Some radial basis functions such as the polyharmonic splines have no shape-parameter.
See also
[edit]References
[edit]- ^ Hardy, Rolland (March 1971). "Multiquadric equations of topography and other irregular surfaces". Journal of Geophysical Research. 76 (8): 1905–1915. Bibcode:1971JGR....76.1905H. doi:10.1029/JB076i008p01905.
- ^ Richard, Franke (January 1982). "Scattered Data Interpolation: Tests of Some Methods". Mathematics of Computation. 38 (157): 181–200. doi:10.1090/S0025-5718-1982-0637296-4. hdl:10945/40152.
- ^ Buhmann, Martin; Nira, Dyn (June 1993). "Spectral convergence of multiquadric interpolation". Proceedings of the Edinburgh Mathematical Society. 36 (2): 319–333. doi:10.1017/S0013091500018411.
- ^ Mairhuber, John C. (1956). "On Haar's Theorem Concerning Chebychev Approximation Problems Having Unique Solutions". Proceedings of the American Mathematical Society. 7 (4): 609–615. doi:10.2307/2033359. JSTOR 2033359.
- ^ Hardy, Rolland L. (1971). "Multiquadric equations of topography and other irregular surfaces". Journal of Geophysical Research. 7 (8): 1905–1915. Bibcode:1971JGR....76.1905H. doi:10.1029/JB076i008p01905.
- ^ Fasshaur, Greg (2007). Meshfree Approximation Methods with MATLAB. World Scientific Publishing. ISBN 978-981-270-633-1.
Radial basis function interpolation
View on Grokipediafrom Grokipedia
Fundamentals
Definition and Motivation
Radial basis function (RBF) interpolation is a numerical technique for approximating a function from a set of scattered data points in multiple dimensions. Given data points where and , the interpolant is constructed as a linear combination of radial basis functions centered at the data points:
where is a univariate radial basis function, denotes the Euclidean norm, are coefficients, is a low-degree polynomial (often zero for positive definite , but required for conditionally positive definite cases to ensure uniqueness), and the coefficients are chosen to satisfy the interpolation conditions for .[7] This formulation leverages the radial symmetry of , meaning the basis functions depend only on the distance from the center , enabling isotropic approximation regardless of the data's spatial orientation.[7]
The primary motivation for RBF interpolation arises in scenarios involving irregularly spaced or scattered data, where traditional methods like polynomial interpolation or tensor-product splines become impractical or inefficient, particularly in high dimensions. Unlike grid-based approaches that require structured meshes and suffer from the curse of dimensionality—needing exponentially more points as the dimension increases—RBF methods operate directly on unstructured data without necessitating triangulation or tessellation.[7] For instance, spline interpolation in multiple dimensions often demands complex domain partitioning, whereas RBFs provide a meshfree alternative that scales flexibly with data distribution and dimensionality.[7]
Key benefits of RBF interpolation include its ability to produce smooth approximants, as most common are infinitely differentiable, and the guarantee of unique solutions when is positive definite or conditionally positive definite, ensuring the interpolation matrix is nonsingular under mild conditions on the data points.[7][8] This uniqueness, established for classes of RBFs like multiquadrics and thin-plate splines, makes the method robust for applications in computer graphics, geophysics, and scientific computing where data irregularity is common.[8]
Historical Background
The theoretical groundwork for radial basis function (RBF) interpolation was established by I. J. Schoenberg in 1938 through his analysis of positive definite functions in metric spaces, which provided the mathematical foundation for functions that ensure stable interpolation properties.[9] Radial basis function interpolation originated in the early 1970s with Rolland L. Hardy's introduction of the multiquadric method in 1971, specifically designed for interpolating scattered meteorological and topographic data in geophysics. Hardy's approach emphasized smooth surface fitting without requiring a mesh, which addressed limitations in traditional grid-based methods for irregular data. This innovation marked the practical beginning of RBF techniques, initially applied in cartography and environmental modeling.[10] In 1982, Richard Franke evaluated multiple techniques for scattered data interpolation, identifying RBF methods as particularly effective for producing smooth and accurate surfaces from irregularly spaced points. C. A. Micchelli provided a rigorous theoretical framework in 1986, proving that certain RBFs generate invertible interpolation matrices via conditionally positive definite properties, enabling unique solutions under polynomial constraints.[8] In the 1990s, RBF interpolation saw significant popularization for scattered data problems, driven by J. C. Mason's editorial work compiling key algorithms and theoretical insights in volumes like Algorithms for Approximation (1990), which highlighted RBFs' versatility. These developments, along with Edward J. Kansa's 1990 extension to solving partial differential equations, facilitated integration into computer graphics for surface reconstruction and geostatistics for spatial prediction, broadening RBFs beyond geophysics. Following the turn of the millennium, computational advances in matrix solving and parallel processing enabled larger-scale RBF applications, spurring further theoretical refinements and widespread adoption in scientific computing.Mathematical Framework
Interpolation Problem Setup
Radial basis function (RBF) interpolation addresses the problem of reconstructing a function from a set of scattered data points in multiple dimensions. Given a set of distinct points $ \mathbf{x}_1, \dots, \mathbf{x}_N \in \mathbb{R}^d $ and corresponding function values $ f_1, \dots, f_N \in \mathbb{R} $, the goal is to find an interpolant $ s: \mathbb{R}^d \to \mathbb{R} $ such that $ s(\mathbf{x}_j) = f_j $ for $ j = 1, \dots, N $. The interpolant takes the form
where $ \phi: [0, \infty) \to \mathbb{R} $ is a radial basis function, $ \lambda_1, \dots, \lambda_N \in \mathbb{R} $ are coefficients to be determined, and $ p $ is a polynomial of low degree (typically total degree at most $ m-1 $) to ensure conditional positive definiteness when required. This formulation is particularly suited for scattered data, as it imposes no grid or mesh structure on the points, unlike finite element or spline methods that rely on connectivity.
The interpolation conditions lead to a linear system of equations. Substituting the form of $ s $ into the conditions yields $ A \boldsymbol{\lambda} + P \boldsymbol{\mu} = \mathbf{f} $ and $ P^T \boldsymbol{\lambda} = \mathbf{0} $, where $ A $ is the $ N \times N $ interpolation matrix with entries $ A_{ij} = \phi(|\mathbf{x}_i - \mathbf{x}_j|) $, $ \mathbf{f} = (f_1, \dots, f_N)^T $, $ P $ is the $ N \times q $ matrix whose columns are the basis functions of the polynomial space evaluated at the $ \mathbf{x}_i $, $ \boldsymbol{\lambda} = (\lambda_1, \dots, \lambda_N)^T $, $ \boldsymbol{\mu} $ are the polynomial coefficients, and $ q = \binom{d + m - 1}{m - 1} $ is the dimension of the polynomial space. For strictly positive definite $ \phi p \equiv 0 $, $ m=1 $), simplifying to the square system $ A \boldsymbol{\lambda} = \mathbf{f} $. Uniqueness of the solution is guaranteed under the assumption that the points are distinct and $ \phi $ is strictly positive definite (or conditionally positive definite of order $ m $ with the polynomial augmentation), ensuring the augmented system is invertible.
This setup operates effectively in any dimension $ d \geq 1 $, providing flexibility for high-dimensional problems where traditional mesh-based interpolation fails due to the curse of dimensionality or irregular data distribution. The method's reliance on pairwise distances via the radial kernel $ \phi $ allows it to handle arbitrarily scattered configurations without preprocessing for structure.
Radial Basis Functions
Radial basis functions (RBFs) are univariate functions that depend only on the radial distance $ r = | \mathbf{x} - \mathbf{x}_i | $ from a center point, exhibiting radial symmetry such that $ \phi(r) = \phi(| \mathbf{x} | ) $ for $ r \geq 0 $.[7] This symmetry ensures translation and rotation invariance in the interpolation process. Many RBFs are positive definite or conditionally positive definite, meaning the associated interpolation matrix is invertible (possibly after augmenting with polynomials for conditional cases) to guarantee unique solutions.[7][1] Several common RBFs are employed in interpolation, each parameterized by a shape factor $ \epsilon > 0 $ (except for certain splines) that controls the function's decay and flatness. The Gaussian RBF is defined as
which is infinitely differentiable and positive definite, providing smooth interpolants with infinite support.[7][1] The multiquadric RBF takes the form
which is conditionally positive definite of order 1 (requiring linear polynomial augmentation for uniqueness) and also infinitely smooth with infinite support, often used for its robustness in scattered data scenarios.[7][1] The inverse multiquadric is given by
positive definite and infinitely smooth like the Gaussian, but decaying more slowly, suitable for approximating functions with moderate variation.[7][1]
For specific dimensions, the thin-plate spline RBF, commonly used in two dimensions ($ d=2 $), is
which lacks a shape parameter and is conditionally positive definite of order 2 (requiring quadratic polynomial augmentation); it is smooth away from the origin but minimizes bending energy in surface fitting.[7][1] Selection among these RBFs involves balancing smoothness against support: for instance, the Gaussian offers superior smoothness but global influence due to infinite support, whereas conditionally positive definite types like the multiquadric may require augmentation yet provide flexibility in handling irregular data.[7] The shape parameter $ \epsilon $ influences this trade-off and matrix conditioning, with optimal values often tuned separately.[1]
Solution Techniques
Direct Methods
Direct collocation methods for radial basis function (RBF) interpolation involve solving the full linear system $ A \boldsymbol{\lambda} = \mathbf{f} $, where $ A $ is the $ N \times N $ interpolation matrix with entries $ A_{ij} = \phi(| \mathbf{x}_i - \mathbf{x}_j |) $, $ \boldsymbol{\lambda} $ are the unknown coefficients, and $ \mathbf{f} $ contains the data values at the interpolation points $ \mathbf{x}_i $. For strictly positive definite radial functions $ \phi $, the matrix $ A $ is symmetric positive definite, ensuring a unique solution that can be obtained via direct linear solvers such as Gaussian elimination or LU decomposition.[7] These direct methods exhibit a computational complexity of $ O(N^3) $ time and $ O(N^2) $ space due to the dense nature of $ A $, making them suitable for small to moderate datasets with $ N $ up to approximately 1000 points, beyond which storage and solution times become prohibitive without specialized hardware.[7] RBF interpolation matrices are often severely ill-conditioned, particularly for globally supported functions like multiquadrics or thin-plate splines, with condition numbers growing exponentially with $ N $ and depending on the node distribution and shape parameter. To mitigate this, Gaussian elimination with partial pivoting is commonly employed to enhance numerical stability during factorization, while preconditioners can further improve conditioning for larger $ N $, though exact solutions remain feasible only for smaller problems.[7] For conditionally positive definite radial functions $ \phi $ of order $ m $, which reproduce polynomials up to degree $ m-1 $, the basic system is augmented with a polynomial term $ p(\mathbf{x}) = \sum_{k=1}^M \mu_k p_k(\mathbf{x}) $ of degree at most $ m-1 $, where $ {p_k} $ form a basis for the polynomial space, and $ M = \binom{d + m - 1}{m - 1} $ in $ d $ dimensions. The coefficients satisfy orthogonality conditions $ \sum_{i=1}^N \lambda_i p_k(\mathbf{x}_i) = 0 $ for $ k = 1, \dots, M $, leading to the block linear system$$ \begin{bmatrix} A & P \ P^T & O \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \ \boldsymbol{\mu} \end{bmatrix}
\begin{bmatrix} \mathbf{f} \ \mathbf{0} \end{bmatrix}, $$ where $ P $ is the $ N \times M $ polynomial matrix with entries $ P_{ik} = p_k(\mathbf{x}_i) $, and $ O $ is the $ M \times M $ zero matrix. This augmented system, of size $ (N + M) \times (N + M) $, preserves positive definiteness under the side conditions and is solved using the same direct techniques as the unaugmented case, with $ M $ typically small compared to $ N $.Iterative and Approximate Methods
Iterative solvers address the challenges of solving the large, dense linear systems arising in radial basis function (RBF) interpolation, particularly for datasets with thousands or millions of points. For symmetric positive definite systems, such as those from multiquadric or thin-plate spline RBFs, the preconditioned conjugate gradient (PCG) method is commonly employed to achieve faster convergence. This approach was first proposed by Dyn, Levin, and Rippa for thin-plate spline interpolation, motivated by variational principles in reproducing kernel Hilbert spaces. For indefinite or nonsymmetric systems, the generalized minimal residual (GMRES) method serves as an effective Krylov subspace solver, often combined with flexible preconditioning to handle variable preconditioners at each iteration. Preconditioning strategies, including diagonal scaling or more advanced techniques like approximate cardinal basis functions and multigrid methods, are crucial to mitigate the ill-conditioning exacerbated by small shape parameters ε. Approximate methods reduce the prohibitive O(N²) computational complexity of exact matrix operations, where N is the number of data points, by exploiting structure in the RBF kernel. The fast multipole method (FMM) accelerates matrix-vector multiplications in iterative solvers through hierarchical expansions and quadrature approximations, achieving O(N log N) or better complexity; a kernel-independent variant using band-limited RBF approximations has been particularly effective for scattered data interpolation. Hierarchical matrices (H-matrices) further enable this by providing blockwise low-rank approximations of the dense interpolation matrix, supporting fast arithmetic operations and the construction of preconditioners such as approximate LU factorizations or inverses, with overall complexity reduced to O(N log N). The radial basis partition of unity (RBF-PU) method decomposes the domain into overlapping subdomains with local RBF supports, solving compact local interpolation systems and blending solutions using partition of unity weights like Shepard's method, which scales well for large, unstructured datasets by limiting each local solve to a small number of points. Domain decomposition techniques partition the computational domain into non-overlapping or overlapping subdomains, solving RBF interpolation restricted to each subset independently and then combining the local approximants to form a global solution. This approach, formulated in the reproducing kernel Hilbert space via von Neumann's alternating projection theorem, converges linearly under mild overlap conditions and enables parallel computation; for instance, it has solved interpolation problems with up to 5 million centers efficiently on standard hardware. Local solves can employ direct factorization for conditionally positive definite RBFs after reformulation into symmetric positive definite systems, avoiding the full global matrix assembly. Kansa's method, originally developed for scattered data approximation, uses unsymmetric collocation with RBFs to form the interpolation system, resulting in a nonsymmetric matrix that is typically solved via iterative methods like GMRES. This technique simplifies implementation by directly collocating the RBF approximant at data sites without symmetry enforcement, and while general nonsingularity proofs are unavailable, singularities are rare in practice for typical configurations.Properties and Analysis
Approximation Theory
Radial basis function (RBF) interpolation is underpinned by a reproducing kernel Hilbert space framework, where the native space serves as the natural function space for the analysis. For a positive definite RBF , the native space is the completion of the span of translates under the inner product induced by the kernel, making it a Hilbert space of functions with finite norm . The RBF interpolant to a function at data points is the unique element in that satisfies the interpolation conditions for all and minimizes the native space norm among all such interpolants. This minimization property ensures that provides the optimal recovery of in the native space norm, as derived from the representer theorem for kernel methods.[11] Error estimates for RBF interpolation rely on the power function, which quantifies the worst-case approximation error. For , the pointwise error bound is given by , where the power function measures the distance from to the data set in the native space geometry, specifically , with the interpolation matrix and the vector of kernel evaluations at . In the uniform norm, the global error satisfies for some constant independent of , assuming lies in the native space. This bound highlights that the approximation quality improves as the power function decreases with denser data sites. Convergence of RBF interpolants to the target function is analyzed in terms of the fill distance , which characterizes the density of the data points in the domain . For RBFs of smoothness order (e.g., Matérn kernels with parameter ), as and , the error converges at rate for , with the constant depending on the separation radius and domain properties. For infinitely smooth RBFs like the Gaussian , the native space embeds into the space of analytic functions, yielding exponential convergence for sufficiently smooth , though practical rates may saturate due to the fixed shape parameter. These rates assume quasi-uniform data distribution and hold in bounded domains.[7] Saturation effects in RBF approximation occur when the error convergence rate does not exceed the intrinsic smoothness order of the RBF, even for smoother target functions , limiting the method to optimal rates within the native space. To achieve higher-order approximation for polynomials, augmented RBF systems incorporate polynomial subspaces: the interpolant takes the form , where is a polynomial of degree at most , ensuring exact reproduction of polynomials up to degree if the data sites unisolvently determine such polynomials. This augmentation enhances saturation bounds, allowing error rates up to for smooth in suitable Sobolev spaces, as seen in polyharmonic spline interpolants.[12]Shape Parameter Tuning
The shape parameter ε governs the flatness and width of the radial basis function φ, directly influencing the trade-off between approximation accuracy and numerical stability in radial basis function (RBF) interpolation. Small values of ε produce wide, flat basis functions that cause the interpolant to resemble a constant approximation and thereby degrading accuracy, particularly for functions with varying scales or high-frequency components, while resulting in an ill-conditioned interpolation matrix , often nearly singular, which amplifies rounding errors and leads to instability. In contrast, large ε values produce narrow, peaked basis functions that enable high-fidelity interpolation with improved numerical stability.[13][14] A primary strategy for tuning ε involves leave-one-out cross-validation (LOOCV), which selects ε by minimizing the sum of squared prediction errors computed by iteratively excluding each data point and evaluating the interpolant from the remaining points at the excluded location. Rippa's algorithm efficiently implements this by deriving a closed-form expression for the leave-one-out errors using the inverse of , avoiding repeated solves of reduced systems and achieving overall complexity, where is the number of points; it has been shown effective for multiquadric, inverse multiquadric, and Gaussian RBFs across two-dimensional datasets. Generalized cross-validation (GCV) provides a robust approximation to LOOCV for exact interpolation, minimizing the score
which uses the diagonals of the inverse interpolation matrix for an unbiased estimate of predictive error without explicit leave-one-out computations.[15][6]
Additional tuning approaches include a priori heuristics based on data geometry, such as ε ≈ 1/(0.815 d) where d is the average nearest-neighbor distance, or more generally ε ∼ 1/h with h the fill distance (the radius of the largest empty ball in the domain), which scales the basis functions to the data resolution. Bayesian methods model the interpolation error as a Gaussian process and employ optimization to find ε that minimizes expected error, offering efficiency for high-dimensional searches. Multi-scale decompositions address varying local features by assigning different ε values across hierarchical levels, with coarser scales using larger ε for global trends and finer scales employing smaller ε for details, enhancing overall stability and accuracy in hierarchical RBF networks.[13][16][17]
For the Gaussian RBF, φ(r) = exp(-ε² r²), the limit ε → ∞ yields narrow, peaked basis functions enabling high-fidelity interpolation in theory, while the limit ε → 0 leads to extreme numerical instability, as the condition number κ(A) grows exponentially with 1/ε, often visualized as a sharp increase for small ε ≈ 0.1 (scaled to data) for moderate N, necessitating careful tuning to balance accuracy gains against conditioning degradation.[13][18]