Gaussian Radial Basis Functions

• Gaussian radial basis functions are infinitely differentiable, positive definite functions that provide spectral (exponential) accuracy for analytic targets.
• They are applied in interpolation and numerical PDE solutions using efficient algorithms, preconditioning techniques, and Hermite polynomial expansions to address ill-conditioning.
• GRBFs play a vital role in machine learning, surrogate modeling, and kernel methods by enhancing model expressivity and facilitating fast numerical computations.

Gaussian radial basis functions (GRBFs) are a fundamental class of infinitely smooth, rapidly decaying functions used extensively in interpolation, approximation theory, machine learning, numerical solutions of partial differential equations (PDEs), and surrogate modeling. A GRBF centered at $c$ with shape parameter $\varepsilon$ is defined as

$$\varphi(x) = \exp \left( -\varepsilon^2 \|x - c\|^2 \right),$$

where $\| \cdot \|$ is usually the Euclidean norm. GRBFs possess spectral convergence for analytic targets, closed-form Fourier transforms, and favorable analytical properties for weak and strong form discretizations. Despite their versatility and accuracy, practical deployment requires careful attention to issues such as ill-conditioning in flat regimes, optimal parameter selection, efficient algorithms for large-scale systems, and theoretical understanding of model expressivity.

1. Mathematical Properties and Representations

GRBFs are members of the class of infinitely differentiable, positive definite radial functions. Their key properties for approximation include:

• Translation invariance: On regular grids, the interpolation or distance matrix inherits a (block) Toeplitz structure, enabling spectral analysis and fast algorithms (Baxter, 2010).
• Spectral (exponential) accuracy: Interpolation of analytic functions with GRBFs converges exponentially in the fill distance, particularly when the function is analytic in a sufficiently large domain (Yarotsky, 2012, Adcock et al., 2022).
• Native Hilbert space: The Gaussian kernel’s native space is a reproducing kernel Hilbert space (RKHS) whose inner product can be explicitly characterized. The standard RKHS for the Gaussian kernel $K_\sigma(x, z) = \exp(-\sigma^2 \|x-z\|^2)$ is equipped with an explicit orthonormal basis and reproducing property (Singh, 2023).

Important analytical results include the existence of closed-form expressions for derivatives and polynomial moments: $$\frac{d^k}{dx^k} \exp(-\varepsilon^2 (x-c)^2) = P_k(x-c) \exp(-\varepsilon^2 (x-c)^2),$$ where $P_k$ is a degree-$k$ polynomial, underpinning exact quadrature and Galerkin formulations (Actor et al., 8 Oct 2024).

Fractional derivatives and integrals of GRBFs can also be expressed in closed form via generalized hypergeometric functions (Mohammadi et al., 2016), enabling high-order meshfree numerical methods for fractional PDEs.

2. Interpolation, Approximation, and Error Analysis

The GRBF interpolant for $f$ at nodes $\{x_j\}$ takes the form

$$s(x) = \sum_{j=1}^N c_j \exp\left( -\varepsilon^2 \|x-x_j\|^2 \right).$$

Error and convergence:

• For functions analytic in a complex domain extending a distance $\rho>0$ beyond $[a,b]$, the univariate interpolation error satisfies

$$\sup_{x \in [a,b]} |f(x) - I_g f(x)| \le c \cdot \left( \frac{b-a}{\rho} \right)^n$$

for some $c$ and $n$ the number of nodes (Yarotsky, 2012).

• Discrete least squares with oversampling (more sample points than centers) and appropriate regularization can achieve accuracy near machine precision, even in regimes where the interpolation matrix is severely ill-conditioned. For 1D, an optimal scaling of the shape parameter with degrees of freedom is $\varepsilon \sim cN$, with $c$ explicitly given (Adcock et al., 2022).

Multivariate and scattered data:

• In high-dimensional or non-tensor settings, the positive definiteness and rapid decay of the Gaussian kernel ensure nonsingularity of the interpolation matrix on distinct centers (Baxter, 2010).
• For large scattered datasets, block partitioning and exploitation of matrix symmetry are critical for memory and computational efficiency (Majdisova et al., 2018).

Specialized interpolants:

• Modified Hermite radial basis function (MHRBF) formulations augment the standard HRBF with polynomial scaling terms, mitigating ill-conditioning and improving accuracy at lower computational cost by requiring only first derivatives and reducing the need for high-degree polynomial augmentation (Fashamiha et al., 21 Feb 2025). For instance,

$$s(x) = \sum_{i=1}^N \left\{ w_i \prod_{j=1}^d (x_j - x_{i,j})^n + (x - x_i)^{2n} b_i \right\} \phi(\|x - x_i\|) + \sum_k \lambda_k p_k(x)$$

with $n \sim 4$ performing optimally.

3. Efficient Numerical Algorithms and Conditioning

Spectral analysis and preconditioning:

• Translation invariance and the Toeplitz structure enable Fourier-domain diagonalization. Quadratic forms $y^\top A y$ with $A_{j,k} = \varphi(|x_j-x_k|)$ can be expressed in the frequency domain via a symbol function $o(\xi)$, providing explicit eigenvalue bounds: $$\text{ess}\, \inf_\xi o(\xi) \leq \text{eigenvalues of }A \leq \text{ess}\, \sup_\xi o(\xi).$$ Preconditioners based on approximating $1/o(\xi)$ through banded Toeplitz matrices minimize the number of conjugate gradient iterations, which can become independent of problem size (Baxter, 2010, Baxter, 2010).

Stable evaluation in the flat regime:

• As $\varepsilon \to 0$, the interpolation matrix becomes increasingly ill-conditioned. Hermite polynomial–based expansions, particularly methods deriving from the generating function (HermiteGF), enable stable evaluation by transferring ill-conditioning to analytically controlled transformation matrices. The basis function is expanded as

$$\exp(-\varepsilon^2 (x-y)^2) = \sum_{n=0}^\infty a_n(y; \varepsilon, \gamma) \psi_n(x; \gamma, \varepsilon)$$

with $\psi_n$ a scaled Hermite polynomial basis (Yurova et al., 2017). This method extends naturally to tensor grids and anisotropic or high-dimensional settings via the Hagedorn generating function.

Hybridization for improved conditioning:

• Hybrid radial basis kernels, such as

$$\phi(r) = \alpha \exp(-(\varepsilon r)^2) + \beta r^3,$$

combine GRBF accuracy with the favorable conditioning of cubic splines. Parameter tuning via global particle swarm optimization further mitigates ill-conditioning while retaining spectral accuracy (Mishra et al., 2016).

4. Applications in PDEs and Numerical Integration

Galerkin and variational methods:

• A variational scheme based on trainable GRBFs exploits analytic polynomial–Gaussian moment integrals, assembling finite element–like mass and stiffness matrices analytically on both bounded and unbounded domains. The weak form inner products are written as

$$\int_{\Omega} (\nabla \phi_i \cdot \nabla \phi_j) dx,$$

which are reduced to polynomial moments against a Gaussian distribution. On $\mathbb{R}^d$ (unbounded domains), integration is exact and the method is conforming; for bounded domains, error estimates balance interior approximation against boundary or consistency error via penalization (Actor et al., 8 Oct 2024).

Boundary element and collocation methods:

• Approximating unknown boundary values in the boundary element method with GRBFs (instead of piecewise polynomials) yields substantial accuracy gains and improved conditioning, especially for high-dimensional domains and singular integral evaluation. Optimizing the source point locations via quadrature error minimization ensures accurate treatment of singularities with standard Gaussian quadrature (Hosseinzadeh et al., 2023).

Kinetic theory and fluid moments:

• In plasma kinetic theory, expanding the velocity-space distribution function in shifted Maxwellians (GRBFs) enables analytic evaluation of Rosenbluth potentials and collision operators, facilitating efficient simulation of Fokker–Planck collision dynamics in both 2D and 3D (Hirvijoki et al., 2015). The GRBF representation also aligns closely with the local thermodynamic equilibrium structure of the underlying physics.

Fourier and diffraction integral computation:

• Representing 2D functions as sums of GRBFs allows semi-analytic evaluation of their Fourier transforms and related optical diffraction integrals via rapidly converging series. The closed-form 2D Fourier transform of a shifted GRBF underpins efficient optical simulations, especially with simultaneous evaluation for multiple defocus parameters (Martinez-Finkelshtein et al., 2015).

5. Primal-Dual Optimization, Machine Learning, and Interpretability

Matrix decomposition and low-rank approximation:

• Nonlinear matrix approximation with GRBF components models a matrix $K$ as a sum of GRBF outer products: $$K_{ij}^0 = b + \sum_{k=1}^r a_k \exp\left( - (u_i^{(k)} - v_j^{(k)})^2 \right)$$ where $u^{(k)}, v^{(k)}$ are optimized vectors. Such decompositions can achieve comparable or better L2 error than SVD with fewer parameters and offer better visual fidelity (e.g., in image compression) and improved community detection in graphs due to their proximity-based structure (Rebrova et al., 2021).

Neural networks with learnable geometry:

• GRBF neural networks equipped with a learnable precision (inverse covariance) matrix replace the usual Euclidean metric of the Gaussian kernel with a Mahalanobis distance: $$\phi( \|x - c_k\| ) = \exp\left( -\frac{1}{2} (x-c_k)^\top P (x-c_k)\right),$$ with $P$ positive definite and parameterized (often via $P=U^\top U$). The eigendecomposition of $P$ reveals the active subspace, allowing dimensionality reduction and interpretable variable importance ranking post-training. This approach is competitive with or superior to standard ML models and deep embedding methods on regression and classification tasks (D'Agostino et al., 2023).

Canonical duality in RBFNN optimization:

• Training GRBF-based neural networks with both centers and weights free leads to highly nonconvex loss landscapes. Application of canonical dual transformations recasts the nonconvex problem into an analytically tractable dual space, enabling complete classification and principled selection of local versus global solutions. Notably, the global minimum of the nonconvex objective may not yield the best generalization; cross-validation among critical points is required (Latorre et al., 2013).

Kernel methods and generalizations:

• The standard GRBF kernel is generalized by compositions with additional exponential terms, yielding a "generalized Gaussian RBF" (GGRBF) kernel with provable strict inclusion of the standard kernel case. RKHS structure, orthonormal bases, and explicit reproducing kernels are constructed, enabling kernel regression, SVM, and deep learning with improved empirical performance over standard GRBF, Sigmoid, and ReLU functions (Singh, 2023).

6. Model Expressivity, High-Dimensional Behavior, and Statistical Insights

Mean dimension and additivity:

• The mean dimension $\nu(f)$ of a GRBF product form

$$f(x) = \prod_{j=1}^d \exp\left(-\frac{(x_j - c_j)^2}{\theta^2}\right)$$

depends on the scale parameter $\theta$; tuning $\theta$ transitions the effective degree of interaction from purely additive ($\nu = 1$ as $\theta \to \infty$) to fully interactive ($\nu = d$ as $\theta \to 0$). In contrast, generalized multiquadrics have $\nu(f) = 1 + O(1/d)$, i.e., they become essentially additive with increasing dimension regardless of other parameters. GRBFs thus have a unique ability to model a continuum from additive to fully interactive structures (Hoyt et al., 2023).

Impact on high-dimensional quadrature:

• For the Keister function $f(x) = \cos(\|x\|/2)$ sampled over $\mathbb{R}^d$, despite symmetric use of all $d$ dimensions, effective (mean) dimension oscillates between 1 and 2 for large $d$, illuminating why quasi-Monte Carlo algorithms perform unusually well in such settings.

Statistical model selection and adaptive placement:

• In GP regression with RBF bases, adaptive basis placement and local hyperparameter inference using $k$-fold cross validation, together with gradient-refined basis placement in high-gradient regions, markedly enhance recovery of discontinuous or high-gradient fields and address classical "ringing" seen with global or harmonic bases (Gregg et al., 2020).

7. Extensions and Future Directions

• Structure-preserving models and Whitney forms: By "lifting" GRBF approximants to represent higher-degree forms (e.g., composed as $\phi_i\nabla \phi_j - \phi_j\nabla \phi_i$), one can build surrogate models that preserve geometric or conservation structure as required by physical law (Actor et al., 8 Oct 2024).
• Analysis of flat limits and relations to classical functions: Connections between GRBFs and classical cardinal functions (e.g., convergence to sinc as parameters go to infinity in multiquadric interpolation) have been established via rigorous Fourier and total positivity analysis (Baxter, 2010).
• Algorithmic developments: Techniques such as HermiteGF expansions, hybrid kernel approaches, banded/circulant preconditioning, and memory-optimized block-computation substantially expand the practical regime of deployable GRBF-based algorithms in high dimensions and at large scale (Yurova et al., 2017, Mishra et al., 2016, Majdisova et al., 2018).
• Intersections with machine-learnable variational algorithms: Recent advances propose schemes that integrate trainable GRBF parameters within exact analytic weak form solvers, yielding meshfree, data- and physics-informed surrogate models with favorable error guarantees, competitive with or surpassing PINNs in smooth settings (Actor et al., 8 Oct 2024).

GRBFs continue to serve as a cornerstone for meshfree numerical methods, surrogate models, statistical learning, and data-driven scientific computing, with ongoing research addressing remaining challenges in conditioning, scalability, and adaptivity across fields.