Radial Basis Function Residual LS (RBF–RLS)

• RBF–RLS is a mesh-free numerical method that uses radial basis functions combined with least-squares minimization to approximate functions and solve PDEs.
• It employs both global and localized schemes with efficient algorithms like SVD/QR regularization and FFT acceleration to achieve robust convergence and stability.
• The method connects deterministic numerical analysis with probabilistic models such as Bayesian RBF networks and LS-SVR, enabling effective uncertainty quantification and parameter recovery.

Radial Basis Function Residual Least Squares (RBF–RLS) is a class of mesh-free numerical approximation and regression methodologies involving radial basis expansions whose coefficients are determined by minimizing the residuals of functional conditions in a discrete least-squares sense. It encompasses both global and localized approaches for high-accuracy function approximation, PDE solution, and regression, providing a framework with well-founded stability and regularization properties, efficient numerical implementations, and deep connections to probabilistic and machine-learning models such as Bayesian RBF networks and LS-SVR.

1. Mathematical Formulation

RBF–RLS uses a linear combination of radial basis functions (typically Gaussians) to approximate a target function or solve a PDE. The standard global ansatz over a bounded domain $\Omega\subset\mathbb R^d$ is

$u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$

where $\phi(r)=\exp(-(\epsilon r)^2)$ is the Gaussian basis with shape parameter $\epsilon > 0$, and $\{\xi_j\}_{j=1}^N$ are the centers, often distributed both inside $\Omega$ and in a thin layer outside to mitigate edge effects and ensure coverage (Adcock et al., 2022).

Given $M\gg N$ collocation or residual points $\{x_i\}_{i=1}^M$, the residuals are defined via linear functionals $L_i$, which can encode pointwise evaluation $L_i[u]=u(x_i)$ (approximation) or differential operators $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$0 (PDEs). The least-squares system is constructed as

$u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$1

and the coefficients $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$2 solve

$u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$3

Associated normal equations $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$4 are solved using a truncated SVD or thresholded QR to regularize high condition numbers inherent to global RBF systems (Adcock et al., 2022, Zhou et al., 2023).

PDE boundary conditions are enforced by stacking additional boundary functionals, yielding an augmented overdetermined system that includes both domain and boundary residuals (Adcock et al., 2022, Larsson et al., 2017).

Localized and partition of unity variants employ overlapping subdomains with compactly supported RBFs, building global approximants from local least-squares fits weighted by Shepard-type windows (Larsson et al., 2017).

2. Theoretical Properties and Stability Analysis

The accuracy and stability of RBF–RLS depend critically on the choice of RBF shape parameter $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$5, center placement, and oversampling ratio $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$6. In the global case, stability and machine-precision accuracy are achieved by using:

• Linear scaling of the shape parameter with the number of degrees of freedom: $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$7. This yields constant overlap of neighboring Gaussians and controls the coefficient norm growth even as $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$8 (Adcock et al., 2022).
• Placement of centers both inside and outside the domain, forming a thin extension layer to prevent 'edge stagnation' and ensure error estimates that do not degrade near boundaries (Adcock et al., 2022).
• Oversampling ratio $u(x) = \sum_{j=1}^N c_j\,\phi(\|x-\xi_j\|),$9, ensuring the discrete residual $\phi(r)=\exp(-(\epsilon r)^2)$0 norm approximates the continuous $\phi(r)=\exp(-(\epsilon r)^2)$1 norm and rendering the system matrix a well-conditioned frame (Adcock et al., 2022, Zhou et al., 2023, Larsson et al., 2017).

Error bounds are established using frame theory; in 1D for Gaussians, sublinear scaling of $\phi(r)=\exp(-(\epsilon r)^2)$2 with $\phi(r)=\exp(-(\epsilon r)^2)$3 yields vanishing error as $\phi(r)=\exp(-(\epsilon r)^2)$4, while linear scaling leads to a reachable machine-precision plateau determined by the SVD cutoff $\phi(r)=\exp(-(\epsilon r)^2)$5 (Adcock et al., 2022). In partition-of-unity and local RBF–FD settings, discrete least squares enables high-order convergence and uniform stability (coercivity of the discrete bilinear form), which classical collocation schemes may lack, especially in the presence of Neumann boundaries (Tominec et al., 2020, Larsson et al., 2017).

3. Efficient Algorithms and Computational Complexity

Recent advances provide efficient solvers for the overdetermined RBF–RLS systems, especially when the RBF centers are placed on regular grids in a bounding box:

• The matrix $\phi(r)=\exp(-(\epsilon r)^2)$6 naturally splits into a 'regular' (circulant or block-circulant) part, admitting acceleration via FFTs, and a low-rank correction reflecting boundary and geometry effects (Zhou et al., 2023).
• The AZ algorithm for RBF–RLS exploits this splitting: it applies a fast solution for the regular part and corrects using a direct small-rank solver for the remainder. In 1D, this leads to overall $\phi(r)=\exp(-(\epsilon r)^2)$7 cost, with the rank of the low-rank correction bounded independently of $\phi(r)=\exp(-(\epsilon r)^2)$8. In higher dimensions, the FFT cost grows to $\phi(r)=\exp(-(\epsilon r)^2)$9 but still yields speedup vs. dense QR solvers (Zhou et al., 2023).

Practical algorithms typically involve:

• Constructing RBF expansion matrices and right-hand sides at oversampled interior and boundary points,
• Regularizing via SVD or QR with an absolute tolerance,
• Evaluating the approximant on arbitrary points via the reconstructed coefficients (Adcock et al., 2022, Zhou et al., 2023).

Partition-of-unity RBF–RLS achieves further efficiency by reusing factorizations across patches and allows control of the stability norm via oversampling, enabling 5–10× speedup over direct collocation at a given accuracy (Larsson et al., 2017).

4. Extensions: Local and Partition-of-Unity RBF–RLS

RBF–RLS generalizes smoothly to several key frameworks:

• Local RBF–FD (finite difference) methods replace global expansions with localized stencils and construct the discrete least-squares system for the cardinal basis. This variant demonstrates superior stability and consistent high-order convergence for elliptic PDEs, particularly on domains with Neumann and mixed boundaries where collocation is unstable (Tominec et al., 2020).
• Partition-of-unity methods construct local RBF–RLS approximants on overlapping patches, blending them with smooth weight functions. Least-squares enforcement at an oversampled set of test points yields error estimates driven by local accuracy and an oversampling-controlled stability norm. This approach is robust against node layout and geometric irregularities, with demonstrated algebraic and spectral convergence in 2D and 3D Poisson problems (Larsson et al., 2017).

5. Applications: Function Approximation, PDEs, and Beyond

RBF–RLS is applied to:

• High-fidelity function approximation on bounded and irregular domains, achieving geometric convergence for smooth targets (Adcock et al., 2022, Zhou et al., 2023).
• Strong-form PDE solution, including elliptic (Poisson, Helmholtz), parabolic, and transport equations, via collocation or weak-form/L$\epsilon > 0$0 residual minimization (Adcock et al., 2022, Larsson et al., 2017, Tominec et al., 2020).
• Handling of singular sources: The RBF–RLS (single-layer RBF network) is well-suited for PDEs with Dirac delta terms in forcing or boundary/initial data. By direct analytic integration in the weak form, residuals associated with singular sources can be driven to zero, unlike PINNs, whose global neural tangent kernel (NTK) inhibits this decoupling. RBF–RLS maintains high accuracy in forward and inverse problems for advection–dispersion with delta sources (Reyna et al., 10 Jun 2026).

The approach extends to inverse problems, where model parameters (velocity, diffusivity, reaction rates) are estimated by embedding the RBF–RLS solve within an outer optimization loop, using analytic gradients available from the linear system structure (Reyna et al., 10 Jun 2026).

6. Statistical Foundations and Relationship to Kernel Methods

RBF–RLS and its variant, Least Squares Support Vector Regression (LS-SVR) with RBF kernels, exhibit precise equivalence to Maximum a Posteriori (MAP) inference in Bayesian RBF networks with Gaussian prior on regression weights (Mesquita et al., 2019): $\epsilon > 0$1 with the regularization parameter $\epsilon > 0$2 corresponding to the prior-to-noise variance ratio. The LS-SVR dual is identical to the equation solved by the RBF–RLS method; the predictive mean and variance of the Bayesian model match the LS-SVR predictor and provide explicit posterior uncertainty quantification (Mesquita et al., 2019).

This connection provides theoretical underpinning for the regularization and generalization behavior of RBF–RLS, enables hyperparameter optimization via evidence maximization, and directly relates the design of RBF kernels in deterministic numerical analysis to those in probabilistic inference.

7. Numerical Results and Implementation Guidelines

Empirical studies confirm the theoretical predictions:

• Global RBF–RLS systems, with optimal shape parameter scaling and sufficient oversampling, achieve geometric convergence to machine precision with coefficient norms remaining $\epsilon > 0$3 (Adcock et al., 2022).
• FFT-accelerated RBF–RLS achieves uniform accuracy and moderate solution time for both function approximation and PDEs, surpassing standard dense solvers beyond moderate $\epsilon > 0$4 (Zhou et al., 2023).
• Local and PUM RBF–RLS solvers sustain high-order and spectral rates with speedups up to $\epsilon > 0$5 over collocation on the same geometry/patch configuration, and maintain stability under node/h-refinement (Larsson et al., 2017).
• In singular-source PDEs, RBF–RLS outperforms PINNs by an order of magnitude in RMSE and achieves robust parameter recovery in inverse modeling (Reyna et al., 10 Jun 2026).
• Key implementation best practices include using an oversampling ratio $\epsilon > 0$6 for global methods, placing an extension layer of centers, using QR/SVD regularization at threshold $\epsilon > 0$7, and (for partition/FD methods) careful control of node placement and patch overlap (Adcock et al., 2022, Larsson et al., 2017, Tominec et al., 2020).

Summary Table: Principal RBF–RLS Paradigms

Variant	Domain/Structure	Regularization/Acceleration	Stability Mechanisms
Global RBF–RLS	Bounded/unstructured	SVD/QR, oversampling, FFT	Shape scaling, halo layers
Local RBF–FD–RLS	Irregular, stencils	Local QR/SVD, sparse normal	Oversampling, local polynomials
RBF–PUM–RLS	Patches, PUM structure	Local QR, partitioned solve	Shepard weights, overlap

These frameworks collectively underpin a versatile, stable, and theoretically justified suite of mesh-free least-squares solvers for approximation, deterministic PDEs, and probabilistic modeling, with robust performance in both classical and modern scientific computing applications (Adcock et al., 2022, Zhou et al., 2023, Tominec et al., 2020, Larsson et al., 2017, Reyna et al., 10 Jun 2026, Mesquita et al., 2019).