Random Feature Model Insights

• Random feature models are finite-dimensional approximations of kernel methods that use a randomized, frozen first layer to generate nonlinear features followed by a linear readout.
• They are crucial in understanding deep learning dynamics near the infinite-width limit and are applied in operator learning, latent-variable modeling, and diffusion-based generative tasks.
• The method leverages spectral concentration, Monte Carlo kernel approximations, and both sparse and adaptive feature designs to enhance performance and regularization.

A random feature model is a parametric approximation to kernel interpolation or regression methods in which a randomized, frozen first layer induces nonlinear features and only a linear readout is trained. In the standard setting, the predictor takes the form $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp\,\phi(x,\omega_k)$$, with random parameters $\omega_k$ sampled independently of the data; equivalently, it is a two-layer network with random frozen first-layer weights and learned second-layer coefficients. Random feature models also play a distinguished role in the theory of deep learning, describing the behavior of neural networks close to their infinite-width limit, but the same construction now appears in operator learning, latent-variable modeling, privacy-preserving interpolation, and diffusion-inspired generative modeling (Nelsen et al., 2020, Liao et al., 2024, Aguirre-López et al., 2024).

1. Canonical formulation and kernel interpretation

For supervised learning with samples $(x_j,y_j)$, a standard random feature regression model uses randomized basis functions and solves a linear problem in the coefficients. In the Fourier-type case, one writes

$$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$

with feature map

$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$

and design matrix

$A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$

In the over-parameterized regime $N\ge m$, the non-private model can be trained by the min-norm interpolator

$$c^\sharp \in \arg\min_{Ac=y}\|c\|_2, \qquad c^\sharp=A^\dagger y=A^*(AA^*)^{-1}y$$

when $AA^*$ is invertible (Liao et al., 2024).

The kernel viewpoint is equally central. By Bochner’s theorem, a continuous shift-invariant positive definite kernel can be written as

$$k(x,x')=\int_{\mathbb R^d} e^{i\langle \omega,x-x'\rangle}\,d\rho(\omega) =\int_{\mathbb R^d} e^{i\langle \omega,x\rangle}\overline{e^{i\langle \omega,x'\rangle}}\,d\rho(\omega),$$

and Monte Carlo sampling with $\omega_k$0 gives

$\omega_k$1

More generally, if $\omega_k$2, then the induced random-feature function space and the RKHS coincide: $\omega_k$3 This makes the random feature model simultaneously a finite-dimensional linear model in randomized coordinates and a finite approximation to kernel learning (Ma et al., 2024).

2. Conditioning, optimization, and implicit regularization

A central theoretical object is the random feature matrix itself. For Fourier-type features

$\omega_k$4

high-dimensional concentration results show that under explicit conditions on dimension, complexity ratio, and variance product, the normalized Gram matrices satisfy

$\omega_k$5

with high probability in the under- and over-parameterized regimes, respectively. Consequently, the singular values of the normalized feature matrix concentrate near $\omega_k$6, which directly controls conditioning, numerical stability, invertibility, and the behavior of ridge or ridgeless regression (Chen et al., 2022).

Finite random feature width also changes the estimator itself. For Gaussian random feature ridge regression with $\omega_k$7 features, $\omega_k$8 data points, and explicit ridge $\omega_k$9, the average predictor is close to kernel ridge regression with an effective ridge $(x_j,y_j)$0 determined by

$(x_j,y_j)$1

where $(x_j,y_j)$2 are the eigenvalues of the training Gram matrix. The key facts are that $(x_j,y_j)$3 and $(x_j,y_j)$4 monotonically as $(x_j,y_j)$5 grows. Finite random feature sampling therefore acts as an implicit regularizer rather than merely a kernel approximation error (Jacot et al., 2020).

The training trajectory itself can also be characterized exactly. In high-dimensional random feature regression under gradient flow, the full time evolution of training and test errors can be calculated analytically through contour-integral formulas and spectral measures associated with the random feature Gram matrix. This yields explicit descriptions of model-wise, sample-wise, and epoch-wise descents, and shows how double and triple descent develop over time (Bodin et al., 2021).

3. High-dimensional spectra, polynomial structure, and depth

One line of work studies how random feature maps transform spectral structure. If the input covariance $(x_j,y_j)$6 has power-law eigenvalue decay

$(x_j,y_j)$7

and the random feature map is a Gaussian sketch followed by the monomial activation $(x_j,y_j)$8, then the population random-feature covariance

$(x_j,y_j)$9

inherits the same power-law exponent. More precisely, for

$$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$0

and for the remainder of the spectrum it remains $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$1 up to polylogarithmic factors. The nonlinearity changes the spectrum only through logarithmic corrections (Paquette et al., 15 Mar 2026).

A complementary statistical-mechanics analysis shows that, for Gaussian covariates, a random feature model is effectively equivalent to a polynomial student with an additional noise term. When $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$2 and $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$3, the random feature model learns polynomial components order by order; $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$4 controls the highest representable order, $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$5 controls how many orders the data can resolve, and interpolation peaks arise when the number of learned parameters becomes comparable to the sample size. This yields staircase generalization curves and a precise account of multiple descent phenomena beyond proportional scaling (Aguirre-López et al., 2024).

Depth does not remove the random-feature character, but it changes the spectral problem. In an $$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$6-layer deep random feature model,

$$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$7

only the last linear readout is trained. A universality theorem shows that the deep nonlinear model is asymptotically equivalent to a deep linear Gaussian model

$$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$8

matched in first and second moments at each layer. The covariance recursion

$$f^\sharp(x)=\sum_{k=1}^N c_k^\sharp \exp(i\langle \omega_k,x\rangle),$$9

demonstrates that depth changes the eigendistribution of the final features and therefore has a tangible effect on performance despite the fact that only the last layer is being trained (Bosch et al., 2023).

A sharper limitation appears in online student–teacher analysis. For any finite ratio $$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$0 between student hidden width and input dimension, the asymptotic generalization error remains strictly positive. Only when the student width is exponentially larger than the input dimension does the probability of containing a random hidden feature highly aligned with the teacher become large enough to permit an approach to perfect generalization (Worschech et al., 2023).

4. Structured, sparse, and adaptive feature design

A large branch of the literature replaces dense, generic random features by structured or sparse constructions. Sparse random feature expansions use a large random dictionary and recover a sparse coefficient vector by basis pursuit with noise tolerance,

$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$1

This gives non-asymptotic high-probability generalization guarantees in data-scarce regimes and improves when the target has coordinate sparsity, compact clusters of the spectrum, or rapid spectral decay. In low-order problems, additional gains come from sparse random weights supported on only $$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$2 coordinates (Hashemi et al., 2021).

HARFE is a related model for high-dimensional sparse additive regression. It uses random sparse connectivity in the random feature matrix to match an order-$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$3 additive target and learns the output coefficients by solving a sparse ridge regression problem with a hard-thresholding pursuit algorithm. The method is guaranteed to converge with an error bound depending on the noise and the parameters of the sparse ridge regression model, and the sparse-connectivity design is intended to match low-order variable interactions (Saha et al., 2022).

Another direction makes the feature distribution itself data dependent. Nonuniform random feature models using derivative information replace conventional uniform sampling of neuron parameters by densities that concentrate where the target’s local derivatives indicate useful neurons. For Heaviside and sigmoid constructions the densities are organized by gradient information; for ReLU and softplus they are organized by Hessian information. This suggests that the gap between uniform random features and fully trained shallow networks is largely attributable to parameter sampling rather than only to finite width (Pieper et al., 2024).

Learnable activation functions provide a second form of adaptivity. In the Random Feature model with Learnable Activation Functions, the activation is represented as a weighted sum of Gaussian radial basis functions,

$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$4

and the finite model becomes

$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$5

This construction can represent a broad class of random feature models whose activations belong in $$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$6, while requiring only about twice the parameter number compared to a traditional random feature model (Ma et al., 2024).

5. Beyond standard regression: latent variables, operators, and generative models

Random features also serve as a modeling device rather than merely a computational approximation. In random feature latent variable models, the nonlinear map from latent coordinates $$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$7 to observations is represented as

$$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$8

with $$\phi(x)=\frac{1}{\sqrt{N}} \begin{bmatrix} \exp(i\langle \omega_1,x\rangle)\ \vdots\ \exp(i\langle \omega_N,x\rangle) \end{bmatrix}\in \mathbb{C}^N,$$9 a random Fourier feature map. The model remains nonlinear in the latent variables but linear in the coefficients, which yields closed-form derivatives with respect to the latent coordinates and supports Gaussian, Poisson, logistic, negative binomial, and multinomial likelihoods. This replaces the exact GP mapping of the GPLVM by a finite random-feature approximation and makes nonlinear latent-variable inference tractable for exponential-family observations (Gundersen et al., 2020).

The same idea extends to operator learning. For an operator $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$0 between Banach spaces, the Banach-space random feature model is

$A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$1

where the random features are $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$2-valued functions on the input space $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$3. In the PDE setting, this defines a non-intrusive data-driven emulator for parameter-to-solution maps. Upon discretization, the model inherits mesh-invariant approximation error with respect to the true PDE solution map and the capability to be trained at one mesh resolution and then deployed at different mesh resolutions (Nelsen et al., 2020).

Random features have also been inserted directly into generative modeling. The Diffusion Random Feature Model replaces the usual diffusion score or noise-prediction network by a time-conditioned random feature architecture. For fixed diffusion time, the model is exactly a random feature model, and the paper derives a total-variation sampling bound whose random-feature contribution decays as $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$4. This gives a mathematically explicit bridge from random feature approximation to denoising score matching and diffusion sampling (Saha et al., 2023).

A different use of the framework appears in missing-data theory. There, a unique underlying random features model keeps the latent data-generating mechanism fixed while the number of observed features varies, allowing a unified analysis of naive zero imputation. In the general infinite-feature setting, the deterministic error from finite random feature approximation, missingness, and naive imputation satisfies

$A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$5

under the well-specified assumption $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$6 (Ayme et al., 2024).

6. Privacy, fairness, and empirical domains

In privacy-preserving learning, the random feature model is studied in the over-parameterized min-norm interpolation regime and then privatized by output perturbation. With

$A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$7

the method achieves an $A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$8-DP guarantee once the Gaussian noise is calibrated to the sensitivity

$A_{j,k}=\exp(i\langle \omega_k,x_j\rangle).$9

The same paper argues that random Fourier features may reduce disparate impact under the excessive-risk-gap metric because

$N\ge m$0

for every input, so the groupwise Hessian-trace mechanism driving disparate impact disappears in the stylized analysis (Liao et al., 2024).

Empirically, random feature models have been deployed on neural spike train recordings, images, and text data; on PDE operator maps such as viscous Burgers’ equation and a variable coefficient elliptic equation; on Fashion-MNIST and instrumental audio in diffusion-style generation; and on benchmark regression problems with sparse additive or low-order structure (Gundersen et al., 2020, Nelsen et al., 2020, Saha et al., 2023). This suggests that “random feature model” is best understood not as a single estimator, but as a family of randomized finite-dimensional surrogates for kernels, operator-valued regressors, latent-variable models, and shallow neural architectures. Across these settings, the recurring themes are explicit feature maps, linear training in the final layer, and theory driven by spectral concentration, kernel approximation, or function-space representation.