Power-Law Spectrum of the Random Feature Model

Elliot Paquette; Ke Liang Xiao; Yizhe Zhu

stat.MLMar 15, 2026

Power-Law Spectrum of the Random Feature Model

Authors: Elliot Paquette, Ke Liang Xiao, Yizhe Zhu

Organizations: Department of Mathematics and Statistics, McGill University · Department of Mathematics, University of Southern California

Abstract

Scaling laws for neural networks, in which the loss decays as a power-law in the number of parameters, data, and compute, depend fundamentally on the spectral structure of the data covariance, with power-law eigenvalue decay appearing ubiquitously in vision and language tasks. A central question is whether this spectral structure is preserved or destroyed when data passes through the basic building block of a neural network: a random linear projection followed by a nonlinear activation. We study this question for the random feature model: given data $x \sim N(0,H)\in \mathbb{R}^v$ where $H$ has $α$ -power-law spectrum ( $λ_j(H ) \asymp j^{-α}$ , $α> 1$ ), a Gaussian sketch matrix $W \in \mathbb{R}^{v\times d}$ , and an entrywise monomial $f(y) = y^{p}$ , we characterize the eigenvalues of the population random-feature covariance $\mathbb{E}_{x }[\frac{1}{d}f(W^\top x )^{\otimes 2}]$ . We prove matching upper and lower bounds: for all $1 \leq j \leq c_1 d \log^{-(p+1)}(d)$ , the $j$ -th eigenvalue is of order $\left(\log^{p-1}(j+1)/j\right)^α$ . For $c_1 d \log^{-(p+1)}(d)\leq j\leq d$ , the $j$ -th eigenvalue is of order $j^{-α}$ up to a polylog factor. That is, the power-law exponent $α$ is inherited exactly from the input covariance, modified only by a logarithmic correction that depends on the monomial degree $p$ . The proof combines a dyadic head-tail decomposition with Wick chaos expansions for higher-order monomials and random matrix concentration inequalities.

Figures & tables

Figure 1: Power-law spectral preservation. Eigenvalue spectra (normalized by $\lambda_{1}$ ) on log-log axes. (a) Iterated linear sketches: starting from a population covariance $D=\mathrm{diag}(j^{-1.31})$ with $v=10{,}000$ (red line), we apply three successive Gaussian sketches ( $10{,}000\to 3{,}000\to 1{,}000\to 1{,}000$ ). All three sketched spectra collapse onto the population power-law; only the tails peel off where the sketch dimension limits the rank. (b) A $4$ -layer $\tanh$ MLP (width $1024$ , random initialization) applied to CIFAR-10 data. All five spectra—input and four hidden layers—collapse onto the same $j^{-1.31}$ power-law, confirming that the spectral exponent is preserved across multiple nonlinear layers on real data.

Figure 2: Median postactivation eigenvalue (normalized so the top eigenvalue equals

1

) across 100 independent trials for Gaussian data with

\alpha=1.31

, using

v=d=10{,}000

and

m=100{,}000

Monte Carlo samples. The dashed line shows the reference

j^{-\alpha}

. Left: monomial activations

f(y)=y^{\text{p}}

for

\text{p}=1,\ldots,6

; all exhibit power-law decay. Right: non-monomial activations (ReLU,

\tanh

, Heaviside,

z^{2}e^{-z^{2}}

) also display power-law spectra, with all but

z^{2}e^{-z^{2}}

tracking the reference slope closely.

Figure 3: Eigenvalue spectra of the random-feature covariance for four data distributions (Gaussian, Rademacher, Student- $t(4)$ , and CIFAR-10) with $\alpha=1.31$ , $d=3{,}000$ , and $m=30{,}000$ Monte Carlo samples. Each panel shows a different activation function. The spectral slopes are consistent across distributions, providing strong evidence of universality. For $x^{2}$ and $x^{3}$ , the grey curve shows a zero-free-parameter theory prediction obtained by combining the lattice-point counts of all composition types in the Wick decomposition; this curve is asymptotically equivalent to the $(\log^{\text{p}-1}(j+1)/j)^{\alpha}$ rate of Theorem 1 but captures the subleading corrections. See Remark 2 for details.

Panel of Figure 4 — Figure 4: If $\pi=(4,2)$ , then there are $4$ copies of $i_{1}$ and $2$ copies of $i_{2}$ . The Feynman diagrams $\gamma_{1},\cdots,\gamma_{6}$ are all mapped to $\eta=(1,1)$ . Which is to say that they all yield the same Wick product $\mathbf{H}_{i_{1}}\mathbf{H}_{i_{2}}\left\langle\mathbf{y}^{(a),2}_{i_{1}}\right\rangle$ . Thus, $N_{(1,1)}=6$ .

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

Panel of Figure 5 — Figure 5: For the set $g_{a_{1}},g_{a_{2}},g_{a_{3}},g_{a_{4}}$ , there exists $10$ distinct Feynman diagrams.

Panel of Figure 6 — Figure 6: The complete Feynman diagrams for $\mathbb{E}\left[\left\langle g_{a_{1}}g_{a_{2}}\right\rangle\left\langle g_{b_{1}}g_{b_{2}}\right\rangle\left\langle g_{c_{1}}g_{c_{2}}\right\rangle\right]$ .

Panel of Figure 7 — Figure 7: The left and right Feynman diagrams correspond respectively to the Wick product given by $\mathbb{E}[g_{a_{1}}g_{a_{9}}]\mathbb{E}[g_{a_{2}}g_{a_{11}}]\mathbb{E}[g_{a_{3}}g_{a_{10}}]\mathbb{E}[g_{a_{7}}g_{a_{8}}]\mathbb{E}[g_{a_{4}}g_{a_{1}2}]\left\langle g_{a_{5}}g_{a_{6}}\right\rangle$ and $\mathbb{E}[g_{a_{5}}g_{a_{9}}]\mathbb{E}[g_{{a_{6}}}g_{a_{11}}]\mathbb{E}[g_{a_{2}}g_{a_{3}}]\left\langle g_{a_{1}}g_{a_{4}}g_{a_{7}}g_{a_{8}}g_{a_{10}}g_{a_{12}}\right\rangle.$

Figure 8: Power-law spectral exponent across network layers for CIFAR-10 data. Each curve shows the fitted slope of the activation covariance eigenvalues at each layer (fit range: indices

1

–

100

), plotted against relative depth. The dashed gray line marks the input data slope (

-1.29

). Dashed curves: random initialization; solid curves: trained networks. For the trained ResNet-20, the slope recovers to the input value at the final layer. The trained CNN shows flattening at middle layers followed by steepening at the final layer.

Architecture	Layer	Dim	Slope	$R^{2}$
Input (CIFAR-10)		3072	$-1.29$	0.998
MLP, ReLU (random)	L1 / L2 / L3 / L4	1024	$-1.14$ / $-1.06$ / $-0.99$ / $-0.95$	$>0.990$
MLP, ReLU (trained)	L1 / L2 / L3 / L4	1024	$-1.22$ / $-1.00$ / $-0.91$ / $-1.08$	$>0.984$
MLP, $\tanh$ (random)	L1 / L2 / L3 / L4	1024	$-1.28$ / $-1.28$ / $-1.27$ / $-1.28$	$>0.998$
MLP, $\mathrm{He}_{3}$ (random)	L1 / L2 / L3 / L4	1024	$-1.27$ / $-1.21$ / $-1.07$ / $-0.83$	$>0.990$
MLP, ReLU+RMSNorm (random)	L1 / L2 / L3 / L4	1024	$-1.08$ / $-0.99$ / $-0.90$ / $-0.85$	$>0.998$

Appendix

Table 1: Power-law slopes and

R^{2}

for activation covariance eigenvalues across layers on CIFAR-10, fitted over indices

1

–

100

. The MLP rows compare three activation functions and two normalization schemes at random initialization, plus a trained ReLU MLP. The randomly initialized ResNet-20 shows remarkably stable slopes (

-1.17

-1.25

) across all layers. For the trained ResNet-20, the final-layer slope (

-1.35

) approximately recovers the input slope.

Figure 9: Power-law slope across layers of a

4\times 1024

MLP at random initialization for three activation functions. The drift rate increases with the effective Hermite degree of the activation:

\tanh

(

\text{p}=1

dominant) shows near-perfect preservation, ReLU (mixed p ) drifts moderately, and

\mathrm{He}_{3}

(

\text{p}=3

) drifts fastest.

Figure 10: Power-law slope across layers of a

4\times 1024

ReLU MLP at random initialization with different normalization schemes. Both RMSNorm and LayerNorm accelerate the drift toward flatter spectra compared to no normalization, with LayerNorm having the strongest whitening effect.

Figure 11: Eigenvalue spectra of the activation covariance at selected layers of a trained ResNet-20 on CIFAR-10. Red lines show power-law fits. The input slope of $-1.29$ is recovered at the final layer (G3-B3).

Figure 12: Final-layer covariance spectra for each architecture, compared with the input. The trained ResNet-20 final layer (rightmost) has slope $-1.29$ , matching the input (leftmost).

Random Feature Methods

Random feature methods approximate kernels and neural networks with finite randomized embeddings to reduce computation while retaining useful accuracy.

See all 46 papers

May 14, 2026

Scaling Laws from Sequential Feature Recovery: A Solvable Hierarchical Model

We propose a simple mechanism by which scaling laws emerge from feature learning in multi-layer networks. We study a high-dimensional hierarchical target that is a globally high-degree function, but that can be represented by a combination of latent compositional features whose weights decrease as a power law. We show that a layer-wise spectral algorithm adapted to this compositional structure achieves improved scaling relative to shallow, non-adaptive methods, and recovers the latent directions sequentially: strong features become detectable at small sample sizes, while weaker features require more data. We prove sharp feature-wise recovery thresholds and show that aggregating these transitions yields an explicit power-law decay of the prediction error. Technically, the analysis relies on random matrix methods and a resolvent-based perturbation argument, which gives matching upper and lower bounds for individual eigenvector recovery beyond what standard gap-based perturbation bounds provide. Numerical experiments confirm the predicted sequential recovery, finite-size smoothing of the thresholds, and separation from non-hierarchical kernel baselines. Together, these results show how smooth scaling laws can emerge from a cascade of sharp feature-learning transitions.

Arie Wortsman-Zurich, Hugo Tabanelli, Yatin Dandi +2

Département d’Informatique, Ecole Normale Supérieure, PSL & CNRS · Information Learning and Physics Laboratory, École Polytechnique Fédérale de Lausanne (EPFL) · Statistical Physics of Computation Laboratory, École Polytechnique Fédérale de Lausanne (EPFL)

Hierarchical Representation Learning Scaling Laws Spectral Methods Statistical Learning Theory

Ask PDF

Jun 16, 2025

Random Matrix Theory for Deep Learning: Beyond Eigenvalues of Linear Models

Modern Machine Learning (ML) and Deep Neural Networks (DNNs) often operate on high-dimensional data and rely on overparameterized models, where classical low-dimensional intuitions break down. In particular, the proportional regime where the data dimension, sample size, and number of model parameters are all large and comparable, gives rise to novel and sometimes counterintuitive behaviors. This paper extends traditional Random Matrix Theory (RMT) beyond eigenvalue-based analysis of linear models to address the challenges posed by nonlinear ML models such as DNNs in this regime. We introduce the concept of High-dimensional Equivalent, which unifies and generalizes both Deterministic Equivalent and Linear Equivalent, to systematically address three technical challenges: high dimensionality, nonlinearity, and the need to analyze generic eigenspectral functionals. Leveraging this framework, we provide precise characterizations of the training and generalization performance of linear models, nonlinear shallow networks, and deep networks. Our results capture rich phenomena, including scaling laws, double descent, and nonlinear learning dynamics, offering a unified perspective on the theoretical understanding of deep learning in high dimensions.

Zhenyu Liao, Michael W. Mahoney

School of Electronic Information and Communications, Huazhong University of Science and Technology, Wuhan, China · ICSI, LBNL, and Department of Statistics, University of California, Berkeley, USA

Random Matrix Theory Neural Network Generalization Neural Network Training Dynamics Neural Scaling Laws High-Dimensional Asymptotics

Ask PDF

Sep 3, 2026

Spectral Convergence of Random Feature Method in Multiple Dimensions

We first prove spectral convergence of the random feature method (RFM) for multidimensional targets in Sobolev, Gevrey, ultra-analytic, and bandlimited classes. The analysis establishes general high-probability approximation estimates in the interpolation scale generated by a kernel integral operator. On a single event determined only by the sampled features, one random space approximates every target in a prescribed source ball; moreover, for each target, a single coefficient vector defines an approximant that attains spectral accuracy simultaneously in all admissible error norms. For both regularity-adapted frequency distributions and uniform distributions on growing frequency windows, the resulting rates range from super-exponential to algebraic, depending on the regularity of the target. Second, we establish abstract error estimates for strong- and weak-form RFM discretizations, thereby converting the preceding approximation bounds into convergence estimates for multidimensional second-order elliptic boundary value and eigenvalue problems. Finally, for random feature matrices (RFMtxs), we prove super-exponential singular-value decay with Fourier features and exponential decay with

\tanh

features, together with corresponding condition-number lower bounds. The analysis identifies a common mechanism: the same spectral approximation that yields high accuracy also drives severe ill-conditioning.

Pingbing Ming, Hao Yu

PDE Solving Random Feature Methods Neural Network Approximation Theory

Ask PDF

Abstract

Random Feature Methods

Explore similar work

Scaling Laws from Sequential Feature Recovery: A Solvable Hierarchical Model

Random Matrix Theory for Deep Learning: Beyond Eigenvalues of Linear Models

Spectral Convergence of Random Feature Method in Multiple Dimensions