Radial Basis Function Interpolation

• Radial basis function interpolation is a mesh-free method that reconstructs smooth functions by linearly combining radially symmetric kernels centered at scattered data points.
• It employs various kernels such as Gaussian, multiquadric, and thin plate splines to balance accuracy and stability while often using polynomial augmentation for exactness.
• Adaptive strategies and scalable solvers, including local methods and fast multipole techniques, mitigate high computational costs, making it effective for complex, high-dimensional problems.

Radial basis function (RBF) interpolation is a mesh-free, dimension-independent methodology for scattered data approximation: it seeks to reconstruct a smooth function, often in high-dimensional or irregular geometries, from its values at a discrete set of spatial nodes. The core principle is to express the interpolant as a linear combination of radially symmetric basis functions centered at the data sites, with coefficients determined so as to satisfy exact interpolation or approximation constraints. RBF interpolation frameworks underpin a diverse range of computational techniques, including meshless discretizations of PDEs, scientific data fitting, machine learning surrogates, and more.

1. Mathematical Framework

Let $\{x_j\}_{j=1}^N \subset \mathbb{R}^d$ be distinct scattered centers, with corresponding data values $f_j$. The canonical RBF interpolant $s(x)$ is

$$s(x) = \sum_{j=1}^N \lambda_j \, \phi(\|x - x_j\|),$$

where $\phi(r)$ is a chosen radial kernel and $\{\lambda_j\}$ are coefficients to be determined. The interpolation conditions $s(x_i) = f_i$ for $i = 1,\dots,N$ yield the dense linear system

$$A \lambda = f, \qquad A_{ij} = \phi(\|x_i - x_j\|).$$

The matrix $A$ is symmetric for most kernels, and its invertibility is characterized by the conditional positive definiteness of $\phi$.

Polynomial augmentation is frequently employed, especially with conditionally positive definite kernels (e.g., polyharmonic splines), via

$$s(x) = \sum_{j=1}^N \lambda_j \phi(\|x - x_j\|) + \sum_{k=1}^P d_k p_k(x)$$

subject to generalized moment (side) conditions, forming the block system

$$\begin{pmatrix} A & P \ P^T & 0 \end{pmatrix} \begin{pmatrix} \lambda \ d \end{pmatrix} = \begin{pmatrix} f \ 0 \end{pmatrix},$$

where $P_{jk} = p_k(x_j)$ for a polynomial basis $\{p_k\}_{k=1}^P$.

2. Kernel Selection and Conditioning

A wide variety of RBF kernels are used, with the most common including:

• Gaussian: $\phi(r) = \exp[-(\varepsilon r)^2]$ (positive definite, infinitely smooth; exhibits severe ill-conditioning for small $\varepsilon$).
• Multiquadric: $\phi(r) = \sqrt{r^2 + c^2}$ (positive definite for $c > 0$; shape parameter $c$ influences flatness and conditioning).
• Cubic: $\phi(r) = r^3$ (conditionally positive definite of order 2; no shape parameter, but may be poorly conditioned for certain node arrangements).
• Thin Plate Spline: $\phi(r) = r^2 \log r$ (conditionally positive definite in 2D; reproduces linear polynomials).
• Wendland's Compactly Supported RBFs: e.g. $\phi(r) = (1 - r)_+^q P(r)$ (compact support; positive definite up to a certain dimension depending on $q$; supports sparse assembly for large data sets).

Hybrid kernels, such as the blend of Gaussian and cubic, $$\phi_H(r) = \alpha \exp[-(\varepsilon r)^2] + (1-\alpha) r^3$$, are used to combine the accuracy of flat Gaussians with the numerical stability of scale-free kernels. At small $\varepsilon$, even a small admixture $(1-\alpha) \ll 1$ substantially lowers the system's condition number, expanding the range of admissible shape parameters for stable interpolation. The optimal $(\varepsilon, \alpha)$ for hybrid schemes are typically selected via global optimization—particle swarm optimization is effective for minimizing cross-validation or $L_2$ validation errors (Mishra et al., 2015).

3. Local, Adaptive, and Geometrically-Informed Algorithms

Global RBF interpolation suffers from $O(N^3)$ computational cost and, for many kernel choices, severe conditioning problems as $N$ grows or as the domain geometry becomes more complex. Several strategies mitigate these challenges:

• Compactly Supported RBFs (CSRBFs): Localize interaction, yielding sparse matrices with $O(N)$ storage and subquadratic or near-linear computational cost (Skala, 2017).
• Local RBF Approaches: For each evaluation point, form and solve a small system on a local stencil (e.g., $k=20$ neighbors), adaptively tuning kernel parameters such as the shape factor $c$ to local point density, as in adaptive MQ-RBF interpolation for imputation of missing geographic data (Gao et al., 2019).
• Geodesic-Aware RBFs: In domains with topological complexity or strong curvature, standard Euclidean distances can produce spurious interpolant oscillations; replacing the RBF argument with a discrete geodesic distance (e.g. via Dijkstra's algorithm on a mesh graph) preserves field continuity and accuracy across domain features (Bucelli et al., 2024).
• Partition of Unity (PU) and RBF-FD Methods: Decompose the domain into overlapping patches/stencils, solving small augmented local problems and blending results for global smoothness. Such approaches are crucial for scaling to millions of nodes or for PDE discretizations (Skala, 2017, Shankar et al., 2018).

4. Theoretical Guarantees and Error Estimates

The convergence rate and stability of RBF interpolation are governed by the kernel's smoothness, the density and uniformity of the data sites, and the norm in which error is measured.

• Native Space and Sobolev Norms: For strictly positive definite kernels with smoothness parameter $\tau > d/2$, the interpolation error in the native Hilbert space satisfies $\|f - s\|_n = O(h^\tau)$, where $h$ is the fill distance. Under additional regularity ($f \in H^{2\tau}$), the "doubling trick" yields improved native-norm rates $O(h^{2\tau})$. For higher-order Sobolev norms $H^s$, new results show $$\|f - s\|_{H^s} \leq C h^{2\tau - s} q^{\tau - s} \|v\|_{L^2}$$, where $q$ is the separation radius and $v$ is the convolution source (Hangelbroek et al., 2023).
• Polynomial Reproduction and the Patch Test: Kernels augmented with polynomials of appropriate degree reproduce polynomials exactly (critical for numerical PDE applications), and ensure nonsingularity of the interpolation system.
• Noise and Robustness: RBF interpolation is robust to small perturbations in node locations (so long as the sequence remains a complete interpolating set) and to additive noise, with error in $L^2$ norm bounded linearly by the noise amplitude (Bouchot et al., 2016).

5. Computational Techniques and Scalability

The primary computational bottleneck in global RBF interpolation lies in solving dense linear systems of size $N \times N$. Large-scale applications necessitate the use of:

• Hierarchical Solvers and Fast Preconditioners: The Discrete Adapted Hierarchical Basis (DAHB) constructs an orthogonal basis that decouples polynomial modes, accelerating Krylov solvers and ensuring scale-invariant conditioning and subquadratic cost (Castrillon-Candas et al., 2011).
• Fast Multipole Method (FMM) for RBFs: Band-limited kernel approximation, together with hierarchical tree structures and kernel-independent translation, reduces matrix-vector multiplication from $O(N^2)$ to $O(N \log N)$, enabling large-scale interpolation and speeding up iterative solvers (Zhao et al., 2016).
• Parallelization and Distributed-Memory Implementations: For use cases such as multiphysics coupling on unstructured meshes, RBF interpolation is implemented in scalable MPI-based frameworks, with near-linear strong scalability for matrix assembly when using geometry-aware kernels (Bucelli et al., 2024).

6. Variants and Extensions

• Hermite RBF (HRBF) and Modified Hermite RBF (MHRBF): For problems requiring interpolation of function and derivative data (such as PDEs with gradient constraints), HRBF extends the standard ansatz to use both kernel values and their derivatives, leading to block systems of higher dimension. MHRBF introduces additional polynomial monomial weights, significantly improving conditioning and accuracy across a wide range of shape parameters (Fashamiha et al., 21 Feb 2025).
• Fractional Derivative RBFs and Conditioning Strategies: Fractional powers and derivatives can emulate the logarithmic behavior of thin-plate splines and interpolate intermediate smoothness classes, with improved stability or conditioning via QR-based preconditioning or the use of "false TPS" radial functions (Torres-Hernandez et al., 2019).
• Mesh Adaptation and Deformation: RBF interpolation is central to mesh morphing methodologies in curvilinear grid generation, with adaptive or iterative shape parameter adjustment guiding element quality and mesh smoothness (Zala et al., 2018).
• Scale-Free RBFs: Cubic and polyharmonic spline kernels do not require a scale/shape parameter, exhibit stable conditioning, and avoid the saturation/ill-conditioning of Gaussian or multiquadric RBFs, especially at high node density (Monnig et al., 2013, Shankar et al., 2018).

7. Applications and Paradigms

RBF interpolation is established as a foundational technique in:

• Meshless Numerical PDEs: Used in meshfree discretizations (RBF-FD, RBF-PU) for both global and local solution of elliptic, parabolic, and hyperbolic PDEs, including on manifolds and complex domains.
• Geoscientific Data Imputation: For estimating missing values in spatial data, adaptive local RBF interpolation (with locally determined shape parameters) provides higher accuracy than traditional methods, at the cost of increased computation (Gao et al., 2019).
• Numerical Quadrature on Surfaces: Polyharmonic RBFs, augmented by polynomials, yield high-order quadrature rules for both flat and curved (manifold) domains; accuracy and stability are governed by the degree of polynomial augmentation and stencil size (Shaw et al., 17 Apr 2025).
• Multiphysics Coupling and Field Transfer: Geometry-aware (geodesic) RBF interpolation is critical for robust transfer of fields between nonmatching multiphysics subdomains, eliminating spurious oscillations near topological features (Bucelli et al., 2024).
• Inverse Nonlinear Dimensionality Reduction: Scale-free cubic RBF interpolants recover stable inverse mappings for nonlinear embeddings, outperforming Gaussian RBFs by avoiding shape parameter tuning and degeneracy (Monnig et al., 2013).

In summary, RBF interpolation is a mature and flexible framework grounded in positive definiteness and invariance, supporting advancements in stability, adaptivity, high-order accuracy, and scalability. The integration of polynomial augmentation, hybrid or adaptive kernel selection, and fast solvers ensures that RBF methods remain at the forefront of meshless numerical analysis, scientific data fitting, and manifold computation (Mishra et al., 2015, Skala, 2017, Castrillon-Candas et al., 2011, Zhao et al., 2016, Shaw et al., 17 Apr 2025).