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

Scaling Laws from Sequential Feature Recovery: A Solvable Hierarchical Model
May 14, 2026Arie Wortsman-Zurich, Hugo Tabanelli, Yatin Dandi +2Hierarchical Representation Learning Scaling Laws
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)
Random Matrix Theory for Deep Learning: Beyond Eigenvalues of Linear Models
Jun 16, 2025Zhenyu Liao, Michael W. MahoneyRandom Matrix Theory Neural Network Generalization
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
Spectral Convergence of Random Feature Method in Multiple Dimensions
Sep 3, 2026Pingbing Ming, Hao YuPDE Solving Random Feature Methods

May 14, 2026J/K move · Enter open · S save

Scaling Laws from Sequential Feature Recovery: A Solvable Hierarchical Model

Arie Wortsman-Zurich, Hugo Tabanelli, Yatin Dandi, Florent Krzakala+1

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

Open paper AskPDF arXiv

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.

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