Local RBF Interpolation Overview

• Local radial basis function interpolation is a meshless technique that constructs approximants from restricted neighborhoods using compact support, partition of unity, or pointwise methods.
• It reduces dense global coupling by employing localized strategies—such as compactly supported kernels and adaptive stencils—to lower computational costs for large, scattered data sets.
• While enhancing numerical scalability and accuracy, these methods rely on careful parameter tuning and domain-specific adjustments to handle stability and conditioning challenges.

Searching arXiv for recent and foundational work on local radial basis function interpolation. Local radial basis function interpolation denotes a family of meshless interpolation schemes in which the approximation at a target location is constructed from a restricted neighborhood, subdomain, or stencil rather than from all data sites globally. The common motivation is computational and numerical: classical RBF interpolation uses the ansatz

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

or, in polynomially augmented form,

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

and enforcing interpolation yields a dense linear system that is expensive to solve directly (Baxter, 2010, Skala, 2017). Local formulations alter this picture in different ways: by using compactly supported radial basis functions so that $\phi(r_{ij})=0$ beyond a finite radius (Skala, 2017), by decomposing the domain into overlapping patches and solving only small local systems combined through partition of unity weights (Cavoretto et al., 2018, Cavoretto, 2014), by forming pointwise neighborhoods for each query location (Gao et al., 2019), or by using local stencils as in RBF-generated finite differences (Londoño-Arboleda. et al., 2018). Although these variants differ in algebraic structure, they share the aim of replacing the dense, globally coupled interpolation problem by sparse, localized, or multilevel computations while retaining the meshless character of RBF approximation.

1. Classical formulation and the case for localization

The standard interpolation problem begins with scattered or gridded samples $\{(x_i,h_i)\}_{i=1}^N$ or $\{x_j,f_j\}$, and seeks an interpolant expressed as a linear combination of radial kernels evaluated at pairwise distances (Baxter, 2010, Skala, 2017). In one common block formulation with polynomial augmentation,

$$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$

where $\Phi_{ij}=\phi(\|x_i-x_j\|)$ and $P$ contains polynomial basis terms (Skala, 2017). This framework is meshless and depends on distances rather than on a mesh or triangulation, which is one reason RBF methods are used for scattered $d$-dimensional interpolation (Skala, 2017).

The principal difficulty is that global kernels such as Gaussian, multiquadric, inverse multiquadric, and thin-plate spline induce full dense matrices (Zhao et al., 2016, Skala, 2017). For large point sets, the storage and solution cost of these matrices is prohibitive. One paper states that for $N \approx 10^6$, a dense $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$0 matrix is clearly infeasible, while the dense/global RBF costs are roughly $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$1 in memory and $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$2 for a direct solve (Skala, 2017). Another paper emphasizes that even matrix-vector multiplication for dense RBF interpolation without fast summation methods costs $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$3 (Zhao et al., 2016).

Local RBF interpolation is therefore best understood as a collection of strategies for reducing global coupling. In some approaches, locality is imposed directly through the kernel support; in others, it is imposed through the computational domain decomposition, neighborhood selection, or transformed basis. This suggests that “local” in the RBF literature is not a single construction but a structural principle: the influence of each basis function, patch, or stencil is restricted in space or in algebraic support.

2. Compact support, sparsity, and neighborhood-limited influence

A direct route to locality uses compactly supported radial basis functions. In that setting, the basis is defined on a normalized support interval by

$$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$4

with a scaling factor $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$5 used in practice so that the support radius is not fixed to $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$6 (Skala, 2017). The decisive local property is explicit: if two points are farther apart than the support radius, then

$$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$7

Each basis function therefore influences only nearby points, and the interpolation matrix becomes sparse rather than dense (Skala, 2017).

This sparse-versus-dense distinction is central to the large-data literature. Global RBFs such as thin-plate splines produce a full dense matrix $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$8, whereas CSRBFs produce a sparse matrix because many entries are exactly zero (Skala, 2017). The practical consequence is not merely reduced memory; sparse direct or iterative solvers become viable, and matrix-vector operations are cheaper. The same source emphasizes that CSRBFs are especially advantageous when the data set is large, the points are scattered or unstructured, a mesh or triangulation would be expensive or impractical, the problem is higher-dimensional, memory limits make dense RBF matrices infeasible, or local interpolation behavior is desired (Skala, 2017).

Compact support does not eliminate all numerical issues. The same analysis warns that large coordinate span harms numerical stability because the polynomial block $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$9 depends strongly on the absolute positions of points. If point coordinates are far from the origin or span a large interval, polynomial terms such as $\phi(r_{ij})=0$0, $\phi(r_{ij})=0$1, $\phi(r_{ij})=0$2, $\phi(r_{ij})=0$3, and $\phi(r_{ij})=0$4 can vary over many orders of magnitude, causing large ranges in entries of $\phi(r_{ij})=0$5, widely spread eigenvalues, ill-conditioning, and reduced numerical robustness (Skala, 2017). A practical implication, stated there explicitly, is that translation and scaling of the point set matter.

This localized-support view of local RBF interpolation should be distinguished from computational locality induced by fast summation. A band-limited fast multipole method for RBF interpolation introduces near-field direct interactions and far-field low-rank approximations, but the underlying interpolation remains global (Zhao et al., 2016). The locality in that case is “local in the FMM sense,” not local in the sense of compact support or patchwise independence.

3. Partition of unity and patchwise local approximants

A second major line of local RBF interpolation is the partition of unity (PU) method. Here the domain is covered by overlapping subdomains $\phi(r_{ij})=0$6, each local point set $\phi(r_{ij})=0$7 defines a small interpolation problem, and the global interpolant is assembled as

$\phi(r_{ij})=0$8

where $\phi(r_{ij})=0$9 are compactly supported PU weights satisfying

$\{(x_i,h_i)\}_{i=1}^N$0

(Cavoretto et al., 2018). This construction replaces one large dense solve by many smaller local systems $\{(x_i,h_i)\}_{i=1}^N$1 (Cavoretto et al., 2018).

The PU framework is explicitly described as effective when large data sets are considered, precisely because only linear systems of relatively small size need to be solved (Cavoretto et al., 2018). In two-dimensional interpolation on $\{(x_i,h_i)\}_{i=1}^N$2, a related cell-based PU algorithm uses local RBF interpolants

$\{(x_i,h_i)\}_{i=1}^N$3

and blends them by

$\{(x_i,h_i)\}_{i=1}^N$4

(Cavoretto, 2014). The paper emphasizes that if each local interpolant matches the data on the nodes in its subdomain, then the PU interpolant preserves interpolation because the weights sum to one (Cavoretto, 2014).

Theoretical analysis in the PU setting depends on a regular covering: bounded overlap, an interior cone condition for each subdomain, and uniform control of local fill distances by the global fill distance (Cavoretto et al., 2018, Cavoretto, 2014). The PU weights are required to be $\{(x_i,h_i)\}_{i=1}^N$5-stable in Wendland’s sense (Cavoretto et al., 2018). These assumptions allow the local errors to control the global approximation error.

An important refinement replaces fixed spherical patches by variable-shaped ellipsoidal subdomains paired with anisotropic kernels (Cavoretto et al., 2018). In that approach, each subdomain $\{(x_i,h_i)\}_{i=1}^N$6 is an ellipsoid specified by semi-axes

$\{(x_i,h_i)\}_{i=1}^N$7

and each local shape parameter becomes a diagonal positive definite matrix

$\{(x_i,h_i)\}_{i=1}^N$8

The rationale is that for track data and other anisotropically distributed point clouds, spherical patches can waste points in irrelevant directions and fail to reflect the local geometry (Cavoretto et al., 2018). Local parameter selection is then driven by a LOOCV error estimate

$\{(x_i,h_i)\}_{i=1}^N$9

optimized patch by patch via the Nelder–Mead simplex method implemented in MATLAB’s fminsearch.m (Cavoretto et al., 2018). The reported numerical evidence is strong: for artificial track data with $\{x_j,f_j\}$0, classical PU yields RMSE $\{x_j,f_j\}$1, whereas PU-LOOCV yields RMSE $\{x_j,f_j\}$2 (Cavoretto et al., 2018).

A distinct computational contribution within PU methods is the use of geometric search structures. In the cell-based 2D algorithm, the domain is partitioned into square cells with side length equal to the subdomain radius, so that each search inspects only a fixed $\{x_j,f_j\}$3 neighborhood of cells (Cavoretto, 2014). The paper reports CPU-time reductions relative to an earlier strip-based method; for $\{x_j,f_j\}$4 and $\{x_j,f_j\}$5, $\{x_j,f_j\}$6 versus $\{x_j,f_j\}$7 (Cavoretto, 2014). This illustrates a recurring theme in local RBF interpolation: localization is often as much about search and data structures as about the kernel itself.

4. Pointwise local interpolation and adaptive neighborhood selection

Another form of locality is pointwise neighborhood interpolation: for each target location, a small local data set is selected and a separate RBF interpolant is built only on that neighborhood. A concrete example is the adaptive local MQ-RBF method for missing-value estimation in geographical data (Gao et al., 2019). For each interpolated point, the algorithm computes its distance to all known data points, selects the 20 points with the smallest distances, and uses these $\{x_j,f_j\}$8 samples as the local interpolation data set (Gao et al., 2019). The local interpolant uses the multiquadric basis

$\{x_j,f_j\}$9

with Euclidean distance $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$0 (Gao et al., 2019).

The distinctive feature is adaptive shape-factor selection based on the distribution of the local neighbors. The method defines the expected density over the full data set,

$$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$1

the local density in the neighborhood,

$$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$2

and the density statistic

$$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$3

This quantity is normalized by a fuzzy membership function to $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$4, after which the shape factor $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$5 is chosen through a piecewise triangular membership rule with five preset levels $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$6 (Gao et al., 2019). The workflow is explicitly local and repeated independently for each missing point.

The method was evaluated on three DEM data sets from Beijing, Chongqing, and Longyan, each with randomly masked $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$7 of observed samples as missing values. The reported data sizes are: Beijing, 1,111,369 known values and 123,592 missing values; Chongqing, 1,074,379 known values and 97,525 missing values; Longyan, 1,040,670 known values and 119,050 missing values (Gao et al., 2019). Relative to $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$8NN interpolation and AIDW, the adaptive RBF method has the best accuracy in terms of RMSE, but is the least efficient computationally (Gao et al., 2019). The source further notes that because the test DEM points are relatively evenly distributed, the advantages of adaptive shape tuning may not be fully realized, and may show more advantage on unevenly distributed data (Gao et al., 2019).

This pointwise model differs from PU in a consequential way. PU constructs a globally defined interpolant by blending local approximants; the adaptive local MQ-RBF method instead performs target-by-target imputation. A plausible implication is that pointwise local RBF interpolation is especially natural for missing-value recovery and localized prediction tasks, whereas PU is more naturally suited to constructing a continuous global field over the entire domain.

5. Local stencils, finite differences, and transformed bases

Local RBF interpolation also appears in discretization schemes where the objective is not a global interpolant but the generation of local differentiation weights. In RBF-generated finite differences for the 2D Helmholtz equation, each node $$\begin{bmatrix} \Phi & P \ P^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \ a \end{bmatrix} = \begin{bmatrix} h \ 0 \end{bmatrix},$$9 is assigned a stencil

$\Phi_{ij}=\phi(\|x_i-x_j\|)$0

formed from the $\Phi_{ij}=\phi(\|x_i-x_j\|)$1 nearest neighbors (Londoño-Arboleda. et al., 2018). On this stencil, the local interpolant is

$\Phi_{ij}=\phi(\|x_i-x_j\|)$2

where $\Phi_{ij}=\phi(\|x_i-x_j\|)$3 is the Bessel function of the first kind used as an oscillatory radial basis (Londoño-Arboleda. et al., 2018). The local coefficients are found from the stencil interpolation matrix

$\Phi_{ij}=\phi(\|x_i-x_j\|)$4

and the differential-operator weights are then derived from this local interpolant (Londoño-Arboleda. et al., 2018).

The paper stresses that local interpolation can itself be ill-conditioned. For the Bessel RBF stencils, condition numbers can exceed $\Phi_{ij}=\phi(\|x_i-x_j\|)$5, especially for larger wavenumbers, symmetric stencil geometries, and regular square or hexagonal node patterns (Londoño-Arboleda. et al., 2018). The stabilization device is a diagonal perturbation,

$\Phi_{ij}=\phi(\|x_i-x_j\|)$6

with $\Phi_{ij}=\phi(\|x_i-x_j\|)$7 chosen so that the perturbed matrix has a prescribed target condition number $\Phi_{ij}=\phi(\|x_i-x_j\|)$8 (Londoño-Arboleda. et al., 2018). The reported convergence rates are about third order and sixth order depending on stencil choice and grid structure, with approximately $\Phi_{ij}=\phi(\|x_i-x_j\|)$9 on a square grid and rates around $P$0 to $P$1 on a hexagonal grid in some cases (Londoño-Arboleda. et al., 2018).

A broader multilevel interpretation of locality appears in the hierarchical-basis solver for RBF interpolation (Castrillon-Candas et al., 2011). There, the kernel may remain globally supported, but an orthonormal hierarchical basis $P$2 is constructed so that the detail basis $P$3 annihilates all polynomials up to a chosen degree and is adapted to the node geometry via an octree-like multiresolution partition (Castrillon-Candas et al., 2011). The authors state that the transformed system behaves like a localized, multilevel approximation rather than a single dense global solve. This produces a transformed matrix $P$4 whose condition number is essentially invariant to domain scaling; for the biharmonic example with $P$5, the transformed matrix had condition number about $P$6 across multiple scales, while the original system condition number reached values like $P$7 (Castrillon-Candas et al., 2011). Combined with fast summation, the total solve is reported as roughly $P$8 under the paper’s assumptions (Castrillon-Candas et al., 2011).

This suggests a useful conceptual distinction. Some local RBF methods are local in the strong sense that only nearby nodes participate in each solve. Others remain globally supported at the kernel level but reveal effective locality after basis transformation or multilevel decomposition. Both belong to the broader efficient-localized RBF tradition, but they address dense coupling in different ways.

6. Geometry-aware locality, structured grids, and parameter tuning

Locality based on Euclidean distance can fail on domains with nontrivial geometry. A recent example is rescaled localized radial basis function interpolation based on geodesic distance, RL-RBF-G (Bucelli et al., 2024). Standard localized RBF interpolation on non-matching meshes uses compactly supported Wendland $P$9 basis functions,

$d$0

with support radius chosen locally as

$d$1

where $d$2 is the $d$3-th nearest neighbor of $d$4 (Bucelli et al., 2024). The interpolant is rescaled by interpolating the constant function $d$5 and dividing: $d$6 This RL-RBF formulation is already local because the basis functions are compactly supported and each support radius is chosen from local point spacing (Bucelli et al., 2024).

RL-RBF-G modifies the distance notion itself. Instead of using Euclidean distance alone, it employs a thresholded hybrid distance based on a discrete geodesic $d$7,

$d$8

This suppresses nonphysical interactions across cuts, cavities, thin walls, or sharply curved regions (Bucelli et al., 2024). In the ring-domain benchmark, standard RL-RBF produces large oscillations near the slit, whereas RL-RBF-G removes these oscillations (Bucelli et al., 2024). The paper reports that for coarse meshes RL-RBF-G reduces the error by about an order of magnitude relative to standard RL-RBF, and that the observed convergence rate is approximately first order in the maximum source mesh size (Bucelli et al., 2024). It also reports nearly ideal strong scalability for matrix assembly up to 3072 cores on GALILEO100 (Bucelli et al., 2024).

Structured-grid settings give yet another variant of localized computation. For gridded data, the interpolation matrix becomes Toeplitz: $d$9 in one dimension and

$N \approx 10^6$0

in $N \approx 10^6$1 dimensions (Baxter, 2010). Although this is not a local method in the compact-support sense, the preconditioner is local in construction: it is formed from a small finite section of the inverse Toeplitz operator, producing a banded symmetric Toeplitz matrix $N \approx 10^6$2 (Baxter, 2010). The reported practical payoff is that the number of PCG iterations is independent of $N \approx 10^6$3 and essentially independent of the right-hand side, and more generally that iteration counts are independent of $N \approx 10^6$4 in practice (Baxter, 2010). This is a localized inverse approximation for a global interpolation operator.

Parameter tuning remains a cross-cutting issue in both local and global methods. Shape parameters affect both accuracy and stability (Mishra et al., 2015, Cavoretto et al., 2023). A Bayesian optimization approach models the validation error as a Gaussian process and uses Expected Improvement to select promising values of $N \approx 10^6$5, reporting a “relevant decrease in search time” relative to LOOCV (Cavoretto et al., 2023). The experiments there are global rather than local, and the paper explicitly states that it does not actually develop or test local RBF interpolation, but notes that the method is conceptually suitable for patchwise local tuning in Partition of Unity schemes (Cavoretto et al., 2023). This suggests that automated parameter selection is increasingly treated as part of the local interpolation problem itself, not merely as an external calibration step.

7. Numerical behavior, advantages, and persistent limitations

Across the literature, the main numerical advantages of local RBF interpolation are consistent. Local or localized schemes reduce dense global coupling, permit sparse storage or small local solves, and are more practical for very large scattered data sets (Skala, 2017, Cavoretto et al., 2018, Cavoretto, 2014). PU methods solve many $N \approx 10^6$6 systems instead of one $N \approx 10^6$7 system (Cavoretto et al., 2018). CSRBF formulations create sparse matrices because many entries vanish exactly (Skala, 2017). Local stencil methods yield sparse global discretizations because each row has only $N \approx 10^6$8 nonzero entries (Londoño-Arboleda. et al., 2018). Geometry-aware localized transfer can eliminate spurious oscillations on complex domains (Bucelli et al., 2024).

Accuracy behavior is more nuanced. For artificial track data and real topographic data, optimized anisotropic-ellipsoidal PU interpolation is reported to be substantially more accurate than classical PU with fixed circular patches (Cavoretto et al., 2018). For local MQ-RBF imputation on DEM data, adaptive local RBF is more accurate than $N \approx 10^6$9NN interpolation and AIDW, but slower (Gao et al., 2019). In large-data interpolation with CSRBFs, sparse support makes the computation manageable, but stability still depends on coordinate scaling and polynomial augmentation (Skala, 2017). In local RBF-FD, regularization can stabilize stencil generation while preserving high-order accuracy (Londoño-Arboleda. et al., 2018).

Several limitations recur. Many localization strategies introduce new parameters: support radius or scaling factor for CSRBFs (Skala, 2017), patch sizes and anisotropic shape parameters in PU (Cavoretto et al., 2018), neighborhood size $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$00 and density-based shape levels in adaptive local MQ-RBF (Gao et al., 2019), regularization parameter $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$01 and target condition number $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$02 in local Bessel RBF-FD (Londoño-Arboleda. et al., 2018), and curvature-detection parameter $$f(x) = \sum_{j=1}^{N} \lambda_j \, \phi(\|x - x_j\|) + P_k(x),$$03 in RL-RBF-G (Bucelli et al., 2024). The literature repeatedly shows that locality alone does not guarantee numerical robustness. Large geometric span, poor scaling, symmetric stencil geometry, and mismatch between patch geometry and data anisotropy can all degrade performance (Skala, 2017, Londoño-Arboleda. et al., 2018, Cavoretto et al., 2018).

A common misconception is that every computationally accelerated RBF method is therefore “local interpolation.” The evidence does not support that equation. Fast multipole acceleration remains a global interpolation scheme with local near-field computation (Zhao et al., 2016). Toeplitz preconditioning for gridded RBF interpolation constructs a localized inverse approximation, but the interpolation operator itself remains global (Baxter, 2010). Hierarchical basis solvers reveal localized block structure without compact support at the kernel level (Castrillon-Candas et al., 2011). By contrast, CSRBFs, PU methods, query-wise local neighborhoods, and RBF-FD stencils are local in a stronger algebraic sense.

Local radial basis function interpolation is therefore best viewed as an umbrella term for meshless interpolation methodologies that control interaction range, either explicitly through support truncation and patching or implicitly through transformed multilevel representations and localized inverse approximations. The unifying theme is the replacement of a single dense global interpolation problem by restricted interactions that are more faithful to data geometry, more scalable computationally, or more stable numerically.