Random Feature Methods Overview

• Random Feature Methods are approximation schemes that use a frozen random nonlinear basis to convert complex nonlinear problems into linear ones for both kernel learning and PDE solvers.
• They enable efficient kernel approximations and statistically controlled regression by replacing high-dimensional kernels with Monte Carlo feature maps and data-dependent compression techniques.
• In PDE applications, RFMs utilize mesh-free collocation, weak-form formulations, and adaptive localized bases to achieve robust and stable solutions while addressing conditioning challenges.

Searching arXiv for recent and foundational papers on Random Feature Methods to ground the article in the literature. Random Feature Method (RFM) denotes a class of approximation schemes in which nonlinear feature functions are sampled once and then frozen, while only linear output coefficients are determined from data or from a differential-equation residual. In the kernel literature, this appears as random feature mapping: a positive-definite kernel is replaced by a finite Monte Carlo feature map, so kernel learning becomes linear learning in feature space (Agrawal et al., 2018). In the PDE literature, the same acronym denotes a mesh-free solver that represents the unknown solution in a randomized shallow trial space and computes the coefficients by least squares, often with partition-of-unity localization and penalty enforcement of boundary or initial conditions (2207.13380). Across these usages, the defining structural property is identical: nonlinear basis generation is randomized and fixed, whereas training or solve-time optimization is linear and therefore convex in the coefficients (Wang, 2019, Kuvakin et al., 1 May 2025).

1. Terminology and common mathematical structure

A standard kernel-side formulation begins from an integral representation

$$\kappa(x,x')=\int \psi(x;\omega)\psi(x';\omega)p(\omega)\,d\omega,$$

and replaces the infinite-dimensional representation by a finite random map $z(x;S_s)\in\mathbb{R}^s$ obtained from sampled $\omega_1,\dots,\omega_s$ (Wang, 2019). For shift-invariant kernels, Bochner’s theorem yields the familiar random Fourier feature form

$$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$

with $\omega_j$ sampled from the spectral density and $b_j$ sampled from $\mathrm{Unif}[0,2\pi]$ (Agrawal et al., 2018).

A standard PDE-side formulation uses an ansatz of the form

$u_M(x)=\sum_{j=1}^{M} a_j\phi_j(x),$

or, in the localized partition-of-unity setting,

$$u_M(x)=\sum_{n=1}^{M_p}\psi_n(x)\sum_{j=1}^{J_n}u_{nj}\phi_{nj}(x),$$

with random features such as $\sigma(\omega^\top x+b)$, $z(x;S_s)\in\mathbb{R}^s$0, $z(x;S_s)\in\mathbb{R}^s$1, $z(x;S_s)\in\mathbb{R}^s$2, ReLU, sigmoids, Gaussians, or compactly supported local bases (2207.13380, Deng et al., 17 Jul 2025). Because the hidden parameters are fixed once sampled, the unknowns are only the output coefficients, and the resulting algebraic problem is linear least squares for linear PDEs and a structured nonlinear least-squares problem for nonlinear PDEs (Beek et al., 21 Jun 2025, Tan, 5 Oct 2025).

The literature therefore uses “RFM” in two closely related senses. “Random Feature Mapping” emphasizes kernel approximation and statistical learning, whereas “Random Feature Method” emphasizes mesh-free collocation, weak-form, or time-stepping solvers for differential equations. The unifying mechanism is a randomized nonlinear trial dictionary plus a linear readout (Agrawal et al., 2018, 2207.13380).

2. Kernel approximation and statistical learning

In large-scale kernel learning, RFMs are introduced to avoid the $z(x;S_s)\in\mathbb{R}^s$3 memory needed to store a kernel matrix and the $z(x;S_s)\in\mathbb{R}^s$4 time typically required for exact kernel training (Agrawal et al., 2018). Replacing the kernel matrix $z(x;S_s)\in\mathbb{R}^s$5 by $z(x;S_s)\in\mathbb{R}^s$6 changes downstream costs to those of linear methods on the feature matrix $z(x;S_s)\in\mathbb{R}^s$7, while preserving the approximation target $z(x;S_s)\in\mathbb{R}^s$8 (Agrawal et al., 2018). For random Fourier features, the Monte Carlo approximation yields uniform error $z(x;S_s)\in\mathbb{R}^s$9 on compact sets (Agrawal et al., 2018).

For kernel ridge regression with random feature mapping, the out-of-sample prediction gap relative to full KRR is controlled under an “unbiased and bounded random features” assumption (Wang, 2019). The central upper bound decays as $\omega_1,\dots,\omega_s$0 in the number of random features, and the corresponding lower bound shows that the $\omega_1,\dots,\omega_s$1 rate is essentially optimal without stronger structural assumptions (Wang, 2019). This establishes that the approximation is not merely heuristic: the gap between RFM-KRR and KRR is quantitatively controlled by the number of sampled features, the regularization level, and the empirical kernel geometry (Wang, 2019).

A major development is data-dependent compression of an initially large random feature pool. The compression procedure of “Data-dependent compression of random features for large-scale kernel approximation” begins from an up-projection of size $\omega_1,\dots,\omega_s$2, constructs a subsampled coreset objective on sampled index pairs, and then applies Frank–Wolfe or GIGA to obtain a sparse weight vector with only $\omega_1,\dots,\omega_s$3 nonzeros (Agrawal et al., 2018). The resulting compressed kernel

$\omega_1,\dots,\omega_s$4

retains the statistical guarantees of the larger random feature model while reducing feature count, memory, and downstream training cost (Agrawal et al., 2018). In the paper’s experiments, the method attains lower kernel approximation error and higher classification accuracy than baseline RFM and JL compression for a fixed compressed feature count, including on a dataset with more than 50 million observations (Agrawal et al., 2018).

Two further lines of work refine the statistical picture. First, “Avoiding The Double Descent Phenomenon of Random Feature Models Using Hybrid Regularization” interprets RFMs as ill-posed linear inverse problems and shows that a hybrid of early stopping and weight decay, with both chosen by generalized cross-validation, removes the spike near the interpolation threshold without a validation set (Kan et al., 2020). Second, “Random Features Outperform Linear Models: Effect of Strong Input-Label Correlation in Spiked Covariance Data” shows that under anisotropic spiked covariance, strong input-label correlation can make RFMs outperform linear models, and the asymptotic behavior becomes equivalent to that of noisy polynomial models whose effective degree depends on the correlation strength (Demir et al., 2024). Together, these papers show that random features are not only a computational surrogate for kernels; their behavior is strongly shaped by regularization and data geometry (Wang, 2019, Kan et al., 2020, Demir et al., 2024).

3. PDE formulations

The PDE literature recasts RFM as a mesh-free collocation or residual minimization method. In the strong-form setting, one samples interior and boundary collocation points, evaluates the PDE operator and boundary operator on the frozen feature ansatz, and solves a weighted least-squares problem for the coefficients (2207.13380). The 2022 bridge paper makes three ingredients explicit: representation of the approximate solution using random feature functions, collocation to treat the PDE, and the penalty method to treat boundary conditions “on the same footing” (2207.13380). The paper emphasizes that multi-scale representation and rescaling of the loss are crucial in practice (2207.13380).

Localized PDE RFMs frequently use overlapping subdomains and partition-of-unity windows. In one formulation,

$\omega_1,\dots,\omega_s$5

where $\omega_1,\dots,\omega_s$6 are local windows and only the coefficients $\omega_1,\dots,\omega_s$7 are trainable (Beek et al., 21 Jun 2025). In another, each subdomain is mapped to a reference box and features are written as $\omega_1,\dots,\omega_s$8, with “partition hyperplanes” defined by the preactivation zero sets (Deng et al., 17 Jul 2025). These localized constructions are central in later adaptive and preconditioned variants (Beek et al., 21 Jun 2025, Deng et al., 17 Jul 2025).

Weak Random Feature Method extends this program to PDEs without strong solutions by replacing pointwise strong residuals with weak-form identities against a finite family of localized test functions (Kuvakin et al., 1 May 2025). The weak formulation is written as

$\omega_1,\dots,\omega_s$9

and the test space is built from windowed sine families on each subdomain (Kuvakin et al., 1 May 2025). In contrast to strong-form RFM with interface derivative matching, WRFM enforces continuity across subdomains but does not impose differentiability conditions (Kuvakin et al., 1 May 2025). This broadens the admissible solution class and is specifically motivated by weak solutions and non-smooth data (Kuvakin et al., 1 May 2025).

Time-dependent problems admit two main RFM formulations. “The Random Feature Method for Time-dependent Problems” introduces a space-time RFM with a space-time partition of unity, two feature constructions—space-time coupled (STC) and separation-of-variables (SoV)—and two solution strategies: a global space-time least-squares solve and a block time-marching variant (Chen et al., 2023). A distinct line, “A Discrete-Time Random Feature Method for Nonlinear Evolution Equations with Implicit-Explicit Runge--Kutta Time Stepping,” keeps time out of the feature input and instead represents the spatial solution at each time level in the random feature trial space

$$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$0

while an IMEX-RK(4,3) scheme advances the state through stage-wise linear least-squares solves (Zhou et al., 28 Apr 2026). These two constructions differ philosophically: one learns a global space-time surrogate; the other advances causally in time while keeping the spatial approximation randomized and linear (Chen et al., 2023, Zhou et al., 28 Apr 2026).

4. Convergence, conditioning, and stability

Several papers now provide explicit convergence statements for RFM. For one-dimensional second-order elliptic problems, “Spectral connvergece of random feature method in one dimension” proves spectral convergence when the solution lies in Gevrey classes and algebraic convergence for Sobolev regularity (Ming et al., 10 Jul 2025). The paper distinguishes super-exponential, exponential, stretched-exponential, and algebraic regimes, and also establishes rates for partition-of-unity-enhanced RFM in terms of patch size (Ming et al., 10 Jul 2025). This is one of the clearest formal analyses connecting smoothness class to RFM approximation rate (Ming et al., 10 Jul 2025).

The same paper gives a sharp warning about conditioning: the singular values of the RFM matrix decay exponentially, while the condition number grows exponentially as the number of features increases (Ming et al., 10 Jul 2025). It further proves that partition of unity can mitigate this excessive singular-value decay by effectively slowing the index-wise decay of the global spectrum (Ming et al., 10 Jul 2025). This makes conditioning a structural issue rather than a purely implementation-level inconvenience (Ming et al., 10 Jul 2025).

Conditioning has become a major theme in PDE RFMs. “Local feature filtering for scalable and well-conditioned Random Feature Methods” applies strong RRQR factorization to local design blocks, discards nonexpressive features, and uses the retained triangular factors to construct a right preconditioner for LSQR on the original least-squares system (Beek et al., 21 Jun 2025). The method operates directly on the rectangular design matrix rather than on normal equations, and its deterministic bound relates the condition number of the patched orthonormalized matrix to overlap interactions between neighboring subdomains (Beek et al., 21 Jun 2025). The same paper reports that in a 1D harmonic oscillator test the preconditioned operator reduces the condition number dramatically relative to the unpreconditioned system and improves both iteration counts and relative $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$1 error (Beek et al., 21 Jun 2025).

For nonlinear PDEs, the least-squares system itself becomes nonlinear and often severely ill-conditioned at scale. “Robust and efficient solvers for nonlinear partial differential equations based on random feature method” introduces an inexact Newton method with right preconditioning (IPN) and an adaptive multi-step variant (AMIPN) based on randomized Jacobian compression, QR-based preconditioning, LSQR inner solves, and derivative-free line search (Tan, 5 Oct 2025). The paper presents these methods as a response to extremely large and ill-conditioned nonlinear least-squares systems in three-dimensional and time-dependent RFM discretizations (Tan, 5 Oct 2025). In that setting, preconditioning is not merely beneficial; it is the mechanism that makes large-scale nonlinear RFM solves practical (Tan, 5 Oct 2025).

Time-dependent theory adds a distinct stability dimension. The space-time RFM paper proves that block time-marching can exhibit an exponentially growing factor with respect to the number of time blocks, whereas the global ST-RFM has an upper bound with a sublinearly growing factor in the number of time subdomains (Chen et al., 2023). The discrete-time IMEX-RFM paper derives a fully discrete error estimate that separates three sources: coefficient perturbations, best-approximation error of the RFM trial space, and third-order temporal discretization error (Zhou et al., 28 Apr 2026). These results make clear that “RFM accuracy” in evolution problems is an interaction between randomized spatial approximation and the chosen temporal formulation (Chen et al., 2023, Zhou et al., 28 Apr 2026).

5. Adaptive and specialized extensions

The recent literature has diversified RFM through adaptive sampling, adaptive bases, multiscale reformulations, and specialized temporal structures.

Variant	Core mechanism	Representative papers
Compressed kernel RFM	Data-dependent coreset compression to $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$2 weighted features	(Agrawal et al., 2018)
Weak RFM	Weak-form equations with localized sine test functions	(Kuvakin et al., 1 May 2025)
Adaptive PDE RFM	Gradient-guided feature redistribution or morphology-aware basis augmentation	(Deng et al., 17 Jul 2025, Hu et al., 10 Oct 2025)
Conditioning-focused RFM	RRQR filtering, randomized Jacobian preconditioning, adaptive Newton solvers	(Beek et al., 21 Jun 2025, Tan, 5 Oct 2025)
Multiscale and asymptotic-preserving RFM	Micro–macro decomposition, adaptive quadrature, Christoffel or liquid sampling	(Chen et al., 2024, Adcock et al., 18 Mar 2026, Linghu et al., 14 Jun 2026)

“AFCM” augments PDE RFM with a monitor function $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$3, then redistributes both partition hyperplanes and collocation points toward high-gradient regions (Deng et al., 17 Jul 2025). The method is explicitly motivated by low-regularity solutions for which standard uniformly sampled features under-resolve sharp peaks, layers, and singular structures (Deng et al., 17 Jul 2025). Its role is analogous to adaptive moving meshes, but the adapted objects are hidden-feature geometry and collocation sets rather than mesh nodes (Deng et al., 17 Jul 2025).

For inverse Helmholtz source problems, “A Morphology-Adaptive Random Feature Method for Inverse Source Problem of the Helmholtz Equation” introduces a two-phase scheme: Integral Adaptive RFM localizes the source support through adaptive quadrature, and a second morphology-enhancement stage adds hybrid basis functions gated by morphology activation functions such as circles, rectangles, torus-like supports, truncated Gaussians, or level-set-derived indicators (Hu et al., 10 Oct 2025). The formulation remains strictly convex in the linear coefficients for fixed morphology activation, because the geometric information enters only through the design matrix (Hu et al., 10 Oct 2025). This extension is specifically aimed at singular solutions and discontinuous or disjoint geometries (Hu et al., 10 Oct 2025).

For multiscale transport, “A Micro-Macro Decomposition-Based Asymptotic-Preserving Random Feature Method for Multiscale Radiative Transfer Equations” decomposes the unknown into macroscopic and microscopic parts and builds separate random feature approximations for each (Chen et al., 2024). The resulting least-squares residual is asymptotic-preserving: as the Knudsen number tends to zero, the residual converges to a well-conditioned diffusion-limit formulation rather than degenerating in the stiff regime (Chen et al., 2024). Here the random feature principle is coupled directly to multiscale asymptotics (Chen et al., 2024).

A different adaptive idea appears in “Christoffel Adaptive Sampling for Sparse Random Feature Expansions,” which uses sparse random feature expansions, constructs the Christoffel function of the currently selected sparse basis, and samples new points from the corresponding Christoffel measure (Adcock et al., 18 Mar 2026). The method integrates active learning with sparse recovery and is designed to reduce sample complexity in data-scarce settings (Adcock et al., 18 Mar 2026). The same adaptive-sampling logic is applicable to parametric differential equations and scientific computing tasks where observations are expensive (Adcock et al., 18 Mar 2026).

Temporal structure itself has also become an adaptive design target. “Liquid Random Feature Methods for Time-Dependent Partial Differential Equations” replaces static temporal activations by closed-form liquid responses with sampled relaxation scales $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$4, producing features that solve a scalar ODE in time and therefore embed temporal scales directly into the frozen trial space (Linghu et al., 14 Jun 2026). The paper’s density theorem and temporal-rank calculation show that distinct sampled relaxation rates contribute independent temporal directions, and the empirical ablations identify the liquid temporal response as the main source of improved finite-feature accuracy in stiff and multiscale regimes (Linghu et al., 14 Jun 2026).

Beyond classical numerical analysis, the acronym now reaches other computational settings. “Quantum Random Feature Method for Solving Partial Differential Equations” keeps the linear least-squares structure but proposes quantum-generated random features and a QLSA-based solution pipeline (Hu et al., 9 Oct 2025). “Random Feature Spiking Neural Networks” adapts the fixed-feature principle to Spike Response Model SNNs and uses a data-driven algorithm, S-SWIM, to construct temporal random features without surrogate gradients (Gollwitzer et al., 1 Oct 2025). These extensions preserve the defining RFM pattern—freeze nonlinear features, fit a linear readout—even though the underlying state space and hardware assumptions differ (Hu et al., 9 Oct 2025, Gollwitzer et al., 1 Oct 2025).

6. Comparative position, applications, and open issues

RFM occupies an intermediate position between classical numerical methods, kernel methods, and deep-learning-based solvers. In kernel approximation, it competes with Nyström and JL-style sketches: Nyström is data-dependent but can scale poorly in high ambient dimension, whereas compressed random features preserve data dependence without Nyström’s $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$5-scaling; JL compression is data-independent and typically yields only constant-factor reductions in feature count, whereas coreset-based compression can reach $$z(x)=\frac{1}{\sqrt{J}}[\cos(\omega_1^\top x+b_1),\ldots,\cos(\omega_J^\top x+b_J)]^\top,$$6 features (Agrawal et al., 2018). In PDEs, RFM is repeatedly contrasted with PINNs: PINNs optimize all network parameters through a nonconvex residual loss, whereas RFM fixes the hidden layer and solves a linear least-squares problem, often with stronger numerical stability and lower training cost for low-dimensional scientific computing problems (2207.13380, Hu et al., 10 Oct 2025, Zhou et al., 28 Apr 2026).

A recurring misconception is that RFM is synonymous with random Fourier features for kernels. The literature surveyed here is broader. PDE RFMs use trigonometric, hyperbolic, ReLU, Gaussian, compactly supported, partition-of-unity-weighted, morphology-gated, micro–macro, and liquid temporal features (2207.13380, Deng et al., 17 Jul 2025, Hu et al., 10 Oct 2025, Linghu et al., 14 Jun 2026). A second misconception is that RFM is intrinsically data-independent. That is inaccurate for compression, adaptive sampling, weak-form test-space construction, Christoffel sampling, local filtering, and morphology-adaptive inverse solvers, all of which use data or residual structure to choose features, points, or preconditioners (Agrawal et al., 2018, Kuvakin et al., 1 May 2025, Adcock et al., 18 Mar 2026, Beek et al., 21 Jun 2025).

Applications now span kernel approximation and classification, kernel ridge regression, support vector machines, kernel PCA, elliptic and Helmholtz equations, radiative transfer, reaction–diffusion, Burgers, KdV, Cahn–Hilliard, Schrödinger, Navier–Stokes, inverse source problems, and parametric differential equations (Agrawal et al., 2018, 2207.13380, Chen et al., 2024, Zhou et al., 28 Apr 2026, Tan, 5 Oct 2025). In multiple PDE papers, RFM is reported to achieve spectral or spectral-like accuracy for smooth solutions, to remain mesh-free on complex geometries, and to compare favorably with PINNs and several classical baselines in the reported benchmarks (2207.13380, Ming et al., 10 Jul 2025, Kuvakin et al., 1 May 2025).

Open issues remain prominent. Conditioning and solver design are central in both kernel and PDE RFMs, especially as feature counts grow or nonlinear residuals are introduced (Ming et al., 10 Jul 2025, Beek et al., 21 Jun 2025, Tan, 5 Oct 2025). Sampling design is still an active area, with Christoffel adaptive sampling, gradient-driven adaptation, and data-dependent compression offering distinct but incomplete answers (Agrawal et al., 2018, Deng et al., 17 Jul 2025, Adcock et al., 18 Mar 2026). Theoretical understanding is also uneven: some regimes now have explicit error and compression theorems, whereas nonlinear, weak-form, multiscale, and adaptive PDE RFMs still rely heavily on empirical evidence or problem-specific analysis (Wang, 2019, Ming et al., 10 Jul 2025, Kuvakin et al., 1 May 2025, Chen et al., 2024). A plausible implication is that the long-term development of RFM will depend less on inventing new random features than on integrating approximation theory, sampling, and numerical linear algebra into feature-selection and solve-time algorithms.