Models of Heavy-Tailed Mechanistic Universality
Authors: Liam Hodgkinson
Zhichao Wang
Michael W. Mahoney
Venue: ICML 2025 poster (2025)
Keywords: heavy tails
random matrix theory
Bayesian learning theory
Primary area: theory->learning_theory
TL;DR: A general family of random matrix models are proposed to explore attributes which give rise to heavy-tailed spectral behavior.

Abstract: Recent theoretical and empirical successes in deep learning, including the celebrated neural scaling laws, are punctuated by the observation that many objects of interest tend to exhibit some form of heavy-tailed or power law behavior. In particular, the prevalence of heavy-tailed spectral densities in Jacobians, Hessians, and weight matrices has led to the introduction of the concept of *heavy-tailed mechanistic universality* (HT-MU). Multiple lines of empirical evidence suggest a robust correlation between heavy-tailed metrics and model performance, indicating that HT-MU may be a fundamental aspect of deep learning efficacy. Here, we propose a general family of random matrix models---the *high-temperature Marchenko-Pastur (HTMP) ensemble*---to explore attributes that give rise to heavy-tailed behavior in trained neural networks. Under this model, spectral densities with power laws on (upper and lower) tails arise through a combination of three independent factors (complex correlation structures in the data; reduced temperatures during training; and reduced eigenvector entropy), appearing as an implicit bias in the model structure, and they can be controlled with an "eigenvalue repulsion'' parameter. Implications of our model on other appearances of heavy tails, including neural scaling laws, optimizer trajectories, and the five-plus-one phases of neural network training, are discussed.

Models of Hlaay-Taallss Mechanistic Universality
Liam Hodgkinson' Zhichao Wang? 23 Michael W. Mahoney 234
Abstract
success of deep learning. Prominent examples include de- velopments in overparameterized models (Jacot et al., 2018; Soltanolkotabi et al., 2018; Belkin et al., 2019; Nakkiran et al., 2021; Bartlett et al., 2021). nonvacuous generalization bounds (Bartlett et al., 2017; Dziugaite & Roy, 2017; Cao & Gu, 2019; Lotfi et al., 2022), statistical mechanics of learn- ing (Martin & Mahoney, 2017; 2021b; Goldt et al., 2019; Bordelon et al., 2024a), random matrix theory (RMT) (Pen- nington & Worah, 2017; Liao et al., 2020; Mei & Montanari, 2022; Hu & Lu, 2022; Wei et al., 2022; Arous et al., 2024; Wang et al., 2024; Atanasov et al., 2025), and robust metrics to assess model quality (Jiang et al., 2020; Martin et al., 2021; Martin & Mahoney, 2021a). In general, however, most "predictions" for model performance provided by the- ory still tend to be ad hoc, providing limited practical utility beyond heuristics. As part of the effort develop effective theoretical explanations, it is important to recognize, char- acterize, and explain empirical properties of deep learning models that are not shared by classical statistical models.
Recent theoretical and empirical successes in deep learning, including the celebrated neural scaling laws, are punctuated by the observation that many objects of interest tend to exhibit some form of heavy-tailed or power law behavior. In particular, the prevalence of heavy-tailed spectral densities in Jacobians, Hessians, and weight matrices has led to the introduction of the concept of heavy-tailed mechanistic universality (HT-MU), Multiple lines of empirical evidence suggest a robust correlation between heavy-tailed metrics and model perfor- mance, indicating that HT-MU may be a funda- mental aspect of deep learning efficacy. Here, we propose a general family of random matrix models- the high-temperature Marchenko-Paie Pastuu (HTMP) ensemble- to explore attributes that give rise to heavy-tailed behavior in trained neural networks. Under this model, spectral densities with power laws on (upper and lower) tails arise through a combination of three independent fac- tors (complex correlation structures in the data; re- duced temperatures during training; and reduced eigenvector entropy), appearing an implicit bias in the model structure, and they can be controlled with an "eigenvalue repulsion" parameter. Im- plications of our model on other appearances of heavy tails, including neural scaling laws, opti- mizer trajectories, and the five-plus-one phases of neural network training, are discussed.
One of the more prominent and ubiquitous of such prop- erties is the frequent appearance of heavy-tailed distribu- tions for various objects of interest. This includes gradient norms (Simsekli et al., 2019) and loss curves (Hestness et al., 2017; Kaplan et al., 2020; Hoffmann et al., 2022). It also includes the singular values (eigenvalues of the product of a matrix with its conjugate transpose) of various matrices, in- cluding data covariance (Sorscher et al., 2022; Zhang et al., 2023), activation (conjugate kernel) (Pillaud- Vivien et al., 2018; Agrawal al., 2022: Wang et al., 2023), Hessian (Xie et al., 2023), Jacobian (Wang al., 2023), and weight matri- ces (Martin & Mahoney, 2021b). In classical settings, these objects typically exhibit (or are assumed to exhibit) Gaus- sian universality (Tao & Vu, 201 Lu & Yau, 2022; Dandi et al., 2023). Indeed, much of classical statistical theory is built around the concentration associated with Gaussian universality. Our objective is to identify a new class of RMT models to describe heavy-tailed spectral behavior in each of these objects. Summarizing contributions:
1. Introduction
Recent years have witnessed remarkable efforts directed to- ward providing a theoretical underpinning to "explain" the
'School Mathematics and Statistics, University Mel- bourne, Australia 'Department of Statistics, University of Cal- ifornia, Berkeley CA. USA è "International Computer Science In- stitute, Berkeley CA. USA 'Lawrence Berkeley National Labora- tory, Berkeley CA. USA. Correspondence to: Liam Hodgkinson <Ihodgkiisoo @wwimmlb.cdu..uuu
Modeling Framework. We propose a general modeling framework for probabilistic analyses of trained neural network feature matrices, including the derivation of their spectral densities.
Proceedings of the ndnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
Beyond the Marchenko- Pastur Law. We examine the High-TTmperature Marchenko-Paan Passtur (HTMP) distribu-

---

Models of Heavy-Tailed Mechanistic Universality
tion (Dung& & Duy, 2021), recently-introduued RMT distribution that generalizes the Marchenko-Pastu (MP) law, and we propose its consideration in the origins of HT-
Consequently, the frequency of heavy-tailed distributions ap- pearing in modern ML suggests new Heavy-Tailed Mech- anistic Universality (HT-MU), a term coined in Martin & Mahoney (2020). Universality in this context describes sys- tems near critical points, whose observables exhibit heavy- tailed, power-law behavior with shared exponents. We build our HT-MU models upon the work of Martin & Mahoney (2020); Martin et al. (2021); Martin & Mahoney (2021b), who identified heavy-tailed spectral distributions in pre- trained weight matrices, and posed the question: what con- stitutes universality in neural network weights? In RMT, universality denotes the emergence of system-independent properties derivable from a few global parameters defining an ensemble (Edelman & Rao, 2005; Edelman & Wang, 2013). In statistical physics, on the other hand, universality arises in systems with very strong correlations, at near a critical point phase transition; and it is characterized by measuring experimentally certain "observables" that display heavy-tailed behavior, with common-or universal-power law exponents. Although trained weight matrices are not truly random, but rather strongly correlated through training, RMT nonetheless provides a useful descriptive framework.
MU. We argue that HTMP arises in the spectral density of feature matrices with non-trivial structure, influenc- ing eigenvalue spacings according to a new parameter k. As more structure imposed, k. decreases, resulting in heavier-tailed spectra. This phenomenon occurs indepen- dently of the behavior of matrix elements- the elements of HTMP models need not have heavy-tailed behavior.
Connections with Existing Heavy-Tailed Properties. We apply the HTMP model to derive neural scaling laws (Kaplan et al., 2020), explain the mysterious appearance of lower power laws in optimizer trajectories (Hodgkin- son al., 2022), and explain the five-plus-one phases of training of Martin & Mahoney (2019).
We begin in Section 2 with brief overview of current ap- proaches to heavy-tailed spectral behavior and, in particular, heavy-tailed mechanistic universality (HT--MU). We intro- duce our modeling framework in Section 3, proposing that trained feature matrices can be modeled by a density of the inverse- inverse-Wishart type. Then, in Section 4, we derive the cor- responding limiting spectral density, identifying three inde- pendent influences on heavy-tailed spectral behavior. Given these results, in Section 5, we apply our model to broader consequences of heavy tails in machine learning (ML). In Section 6, we provide a brief conclusion. Additional mate- rial, including proofs, may be found in appendices.
HT-MU has several important practical consequences. Most famously, spectral power laws in the activation matrices of deep neural networks give rise to neural scaling laws that are obeyed by the test loss with respect to both the number of parameters and training time (Bahri et al., 2024; Maloney et al., 2022). These laws reveal powerful selection criteria for designing compute-optimal models under a given budget (Kaplan et al., 2020; Hoffmann et al., 2022; Muennighoff et al., 2023; Paquette et al., 2024). In the absence of signifi- cant compute or the associated training/testing data, power laws encountered in the eigenspectra of weight matrices (and in the reciprocal of gradient norms; see Hodgkinson et al. (2022)) have been found to be unusually strong in- dicators of model quality. In particular. HT-MU underlies Heavy-Tailed Seff-Regularization (HT-SR) Theory (Martin & Mahoney, 2021b), providing framework for predict- ing trends in the quality of state-of-the-art neural network mithls-without access to training or testing data (Martin et al., 2021; Yang et al., 2023). Following these principles, practitioners can diagnose and improve werroveee--thhoe models down to the
2. Heavy-Tailed Mechanistic Universality
Heavy-tailed distributions refer to broad range of dis- tributions (Resnick, 2007; Nair et al., 2022), most often with densities that decay slower than exponential, but the phrase is often (informally) used interchangeably with the prescription of power laws, which decay at polynomial rate, f(z) ~ cr as x + oo." Importantly, as an em- pirical matter, many other densities also exist that can ex- hibit tails that are empirically indistinguishable from power laws (Clauset et al., 2009). These include log-normals (e.g., in layer norms (Hanin& & Nica, 2020)) and exponentially- truncated power laws of the form fxx ~ cz cx- BB for small 3 >0 0 (which can provide more informative fits for the spectra of weight matrices (Yang et al., 2023)).
-with very modest
amounts of compute- cumpuuteuundrrpeffrrmni level of individual layers (Zhou al., 2023; Lu et al., 2024).
Approaches to HT-MU. Physical explanations for HT- MU range from the phase boundary of spin glasses. to di- rected percolation (Bouchaud & Potters, 2003), to self orga- nized criticality (SOC) at the edge of chaos (Bertschinger et al., 2004; Cohen et al., 2021), to jamming transitions crossing from underparameterizzed to overparameterized models (Geiger et al., 2019; Spigler et al., 2019). It is chal- lenging to translate these into a rigorous statistical model. On the other hand, RMT practitioners tend to examine the
In statistical physics (Sornette, 2006; Bouchaud & Potters, 2003) and the statistical physics of learning (Seung al., 1992; Watkin et al., 1993; Engel & den Broeck, 2001), the presence of power laws is particularly special, suggesting an underlying universal mechanism driving their appearance.
PPwerr laws can also occur at the bottom part of the density, taking the form f(x) ou ex. x -> 07, e.g.. for gradient norms, as in Hodgkinson al. (2022), and this will be relevant for us.
2

---

Models of Heavy-Tailed Mechanistic Universality
Power Law
Inverse Gamma
feature matrices are not independent (indeed, that was the original motivation for HT-SR theory and the introduc- tion of HT-MU); and second, empirical results demon- strate that, while eigenvalues of weight matrices modern state-of-the-ar models are heavy-tailed, elements are not (Martin et al., 2021; Martin & Mahoney, 2021b).
Mechanism
Elements
Spectrum
jid Heavy-Tailed Elements
Kesten Phenomenon
dxx
Population Covariance
V1x
d1x
Structured Matrices (Thm. 4.1)
X
Empirical Observations (Features)
V
Empirical Observations (Weights)
X
Kesten Phenomenon. In natural systems, the ubiquity of power laws is often attributed the SOC hypothesis, which asserts that dynamical systems in the neighbor- hood of a critical point exhibit power law behavior (Bak et al., 1987). In probability theory, this phenomenon arises from a mechanism discovered by Kesten (1973), where recursive systems on the edge of stability exhibit heavy-tailed stationary fluctuations. Considered with respect to gradient descent steps, Kesten's 's mechanism can explain the origins of heavy-tailed size fluctuations in the stochastic optimizer (Hodgkinson & Mahoney, 2021; Gurbuzbalaban et al., 2021); and analyses treat- ing the neural network architecture itself as recursive reveal similar findings (Vladimirova et al., 2018; Hanin & Nica, 2020; Zavatone- Veth & Pehlevan, 2021). Here, too, though, theoretical results disagree with empirical results: it is the distribution of eigenvalues that is heavy- tailed, not the elements (Martin et al., 2021; Martin & Mahoney, 2021b). Indeed, the Kesten phenomenon in its currently-studied form primarily seems to occur for chaotic training behavior (Yang et al., 2023). It is an open question how to extend the results of Hodgkinson Mahoney (2021); Gurbuzbalaban et al. (2021) to perform Kesten-like iteration in the spectrum domain using free probability (which is of interest since it may reveal heavy tails in eigenvalues but not elements).
Table 1. Comparison of various mechanisms, and their capacity to yield power laws, in feature matrix elements and feature matrix spectral densities, as well as their capacity to yield an inverse Gamma law for the spectral density in neighborhood of zero. Empirical observations on feature matrices appear in Figure 1. Empirical observations weight matrices appear in Martin et al. (2021); Martin & Mahoney (2021b).
interactions of individual models (Li & Sompolinsky, 2021; Wang et al., 2023), even under heavy-tailed data (Maloney al., 2022; Bordelon et al., 2024b). Such analyses directly consider the effects of certain model architecture choices (e.g., depth) but they often do not apply to trained models, appear intractable for complex architectures, and obscure the underlying "reasons" for HT-MU we seek to identify.
To develop models of HT-MU, and to account for the broad range of objects in deep learning that exhibit heavy-tailed behavior, we focus on a class of feature matrices, encom- passing activation, neural tangent kernel (NTK), and Hes- sian matrices. These are fundamental theoretical objects upon which other secondary ("observable") objects, includ- ing weight matrices, gradient norms, and loss curves, all depend. As in Martin & Mahoney (2021b), we use the ma- chinery of RMT as a "stand-in" for generative model of these feature matrices in state-of-the-art neural networks.
Our search for theories to underpin HT-MU revealed three primary categories of observable phenomena, which we outline and compare below and in Table 1. Most promi- nently, we consider the weight matrices of trained neural networks: such weight matrices are well-known to have spectral distributions that are strongly heavy-tailed, while having elements that not heavy-tailed (Martin et al., 2021; Martin & Mahoney, 2021b). We also identify (in the right-most column of Table 1) further universal property we observe in the spectra of NTK feature matrices: the left edge displays an inverse Gamma law. 2 See Figure 1 for a representative summary of the results. Additional details are provided in Appendix A.
Heavy-Tailed Population Covariance. A popular hy- pothesis is that the complex correlations exist in the data, are expressed as heavy-tailed spectra, and are passed onto features during training. This shifts the origin of HT-MU from the model to the data. Covariances of large datasets often exhibit power law spectra, and analyses centered around this approach have had predictive suc- cess for simple models (Paquette et al., 2024; Li et al., 2023). Correlations in the data clearly play a significant role, but implicit model biases should also play role. Otherwise, different model architectures trained on the same dataset should exhibit similar power laws (see Sec- tion 4.3), which empirical results do not display (Martin et al., 2021; Yang et al., 2023).
Independent Heavy-Tailed Elements. It is known that independent matrix elements exhibiting power laws with small tail exponents give rise to heavy-tailed spectral den- sities (Arous& & Guionnet, 2008). This is less relevant for our discussion for two reasons: first, the elements of real
A New Approach using Structured Matrices. We pro- pose an alternative approach, one which considers the effect of implicit model bias exhibiting structure in the feature matrices. Our claims regarding the observed tail behavior of feature and weight matrices (as in Table 1) are summarized in the following metatheorem.
"The reciprocal of the exponentially truncated power law is the inverse Gomma law: fx) ~ ex e A/e as x + 0t.
3

---

Models of Heavy-Tailed Mechanistic Universality
VGG11
ResNet9 dks =
ResNet18
0.0065
dks 0.0255
10a
1:85 B - x000 10
0.75
a= 1.76
8 1.5
8 --244 10
10-1
-1
10
10 4
10-1_
10
10"
107
10-2
10
10a
10
10
10f
10
10a
10a
10
10
10a
10*
10a
10a
(a) Initialization: CIFAR-I0 & VGG11
(b) Initialization: CIFAR-10 & ResNet9
(c) Initialization: CIFAR-10& & ResNet18
VGG11 des = 0.0761
ResNer9
ResNet18 des = 0,0378
0.0464
10
10
10-1
8-1.7x 10-
10
10-1
10-3
10-3
10"4.
10-5
10-1
10-7
10
10
10"
101
10a
10
10'
10'
(d) Post-training: CIFAR-I0 VGG11
(e) Post-training: CIFAR-IO & ResNe19
(f) Post-training: CIFAR-10 & ResNet18
Figure 1. Distributions of spectral values (and inverse Gamma fits near zero) of the NTK matrix for VGG11, ResNet 9, and ResNet18 models trained on 1000 sundoml--sanpled datapoints from the CIFAR- 10 dataset, at initialization and post-training.
Stochastic optimization models of the form (1) have been considered previously (Mandt et al., 2016), and they have strong links to Bayesian inference (see Germain et al. (2016) and Appendix B), which itself has strong links to the statis- tical physics of generalization (Mezard & Montanari, 2009).
Metatheorem (Theorem 4.3 and Proposition 5.3).
Spectral densities of trained feature matrices can be modeled the free multiplicative convolution the label covariance spectrum and the reciprocal of an HTMP distribution, exhibiting both an inverse- inverse-Gamma law at O and a power law at oc.
In particular, applying (1) to the training loss optimizes a PAC-Bayes bound on the test error (see Appendix B). Equa- tion (1) generalizes the approach of Xie et al. (2023), is equivalent to early-stopped graph heat-kemel diffusion (Ma- honey & Orecchia, 2011), and (in ridge regression) co- incides with known regularizers for early stopping with stochastic gradient descent (Sonthalia et al., 2024). As T de- creases during training, an optimizer following (1) smoothly interpolates between Te and the final "optimal" density.
Spectral densities trained weight matrices can be modeled with the HTMP distribution, exhibiting an sxponentially-truncaaee power law at cc.
3. Modeling Framework
Tracking the precise spectral evolution of neural networks throughout training is generally challenging. both due to the- oretical complexity and practical limitations (e.g., limited access to large-scale training dynamics). Instead, we pro- pose a general modeling approach grounded an entropic regularization perspective on stochastic optimization.
Since the features are of primary interest in our study, split an arbitrary learning task into subproblem of finding optimal model coefficients for a prescribed set of features, and a main problem of identifying those features. For ex- ample, as in Golub & Pereyra (1973), neural network can be separated into its final layer and its final activations (rep- resenting trained features). Let L be loss function that depends collection of model coefficients θ and features i with initial densities, тo and πg, respectively, where me can depend conditionally on Φ. In general, the optimally- trained features satisfy x' = argming mine L(O, Φ), if a unique minimizer of L in exists for each 4. Replac- ing deterministic minimization in this expression with our stochastic optimizer (1), the probability density of "optimal features" x" for L(O, b) is computed in Proposition 3.1. We present its proof in Appendix C.I.
3.1. Entropic Regularization Setup
Let θ be model coefficients with an initial density me To allow for possible early stopping in the optimization, we consider a stochastic optimizer that minimizes the loss by monotonically reducing the Bullback-Leibbe divergence to distribution of loss minimizers, representing "optimal" behavior of stochastic gradient descent (Chaudhari & Soatto, 2018). Define the stochastic minimization operation for any function f()), sming T with temperature 7 > 0 as:
FR.7 smin f(0) = min E (())) + TPP
(1)
where P is the set of probability densities on the support of ne. We let argsmin se." denote the corresponding mini- mizing distribution ((),, provided it exists and is unique.
Proposition 3.1 (OPTIMAL FEATURE DENSITY), De- note the minimizer of the feature density as g()) := xa.7 J0.t L(O, ②), where smin is defined by argsmin smin

---

Models of Heavy-Tailed Mechanistic Universality
(I) and 7 are the temperature parameters for O and 4,
This is in the form of the Master Model Ansatz, with
respectively. The probability density function q( satisfies
= pm/2 and / = poo /2.
Neural Tangent Kernel (NTK), Consider the NTK
(( x (.)
Gram matrix ( E R Xmn, where each m x m block is given by J(D); = DSe, (x;) ( De (xz), where Dfe E Rdx is the Jacobian of the model fe.a : X + Y C =Rm with parameters @ E R Coined by Jacot et al. (2018), the NTK is central to the analyses of generalization performance in neural net- works (Huang & Yau, 2020). The NTK approxima- tion, treated in Rudner et al. (2023) and Wilson et al. (2025), is the linearized model fe, (2) ~ fer, + (z) (0 - -²), under the loss
where (D) is marginal (Gibbs) likelihood with prior re:
(Ø) = Ee-x, exp (- -L,,
Proposition 3.1 suggests that stochastic optimizer concen- trates on regions with high marginal likelihood (also called model evidence), with the degree of concentration deter- mined by the ratio T 7/n.. Importantly, in this result, feature learning (controlled by 77) allowed to occur at different rate than coefficient learning (controlled by T).
3.2. Master Model
L(O, 4) = fl (X) - -I
(4)
We now apply Proposition 3.1 to reveal densities for three popular types of feature matrices: activation matrices; NTK matrices; and Hessian matrices. Of particular interest are the late stages of training, where 7, 77 + 0+, and we will also assume that r/n - p 0. In each case we present below, we observe that the trained feature matrix generally follows an önverse-Ve mmvers-Wishart--yee Wishart-type density (Mardia et al., 2024, $3.8), of the form (2). Hence, we propose a Master Model Ansatz for the probability densities of trained feature matrices, with some parameters o. B and initial density m:
By the same arguments used to derive (3). =velY- fe- (X)) E Rmm, then
27(2) x exp(- (A0((Y(022(00000 + (0-)(((0))))))-~>))) 1)))
det(o2](()+ + 11)
Applying Proposition 3.1 for M = ={()), we obtain
((M) x (det -/2 exp (- ( (M). 2
po²
(5)
Master Model Ansatz. For feature matrices M, g(M) (d) (e--)-)-------- (MAEEEMMM))(())))) (M).. (2)
which is also in the form of the Master Model Ansatz, with a = p/2 and 3 = =po²/2 Given that the Fisher in- formation matrix possesses the same nonzero eigenvalues as the corresponding NTK, we can derive a comparable result for the Fisher information matrix well; we omit the details for brevity. The NTK approximation is known to be effective for models in the late stages of training (Fort et al., 2020). A question is: can we justify (5) with- out appealing to a linear approximation? To answer this affirmatively, consider the setting where d >> mn, where Hodgkinson al. (2023) have developed asymptotic ap- proximations to the marginal likelihood under the loss (4) for very general models in the interpolating regime. Their arguments lead to similar finding to (5), when the model is trained to implicitly regularize a lower bound on its variance. See Appendix D for details.
Activation Matrices. Following Golub & Pereyra (1973); Pillaud-Vivis et al. (2018), consider a multilayer neural network with m outputs. ff : X + y C R" which is parameterized in terms of its final linear layer, W (uyyk) E d and a feature vector y : X + RRd of the last hidden layer. That is, the k-th output of f is defined by f((())= = E Wyx4p)()) (1:) for l E [m]. For a dataset D - [1i, yi J= X x V, we can investigate the ridge regression problem of minimizing L(W, b) = | - 11 + M [2:, where μ > 0, kj = (2,) and Y = ( E R= For тe = (N I) and
G²_ a²
1+ z the marginal likelihood takes the form
Hessian Matrix. For more general losses T.(O, ()) com- puting the marginal likelihood becomes intractable, but can be well-estimated by Laplace approximation in the regime where 7 + 0+. Letting & denote appropriate hyperparameters or architecture choices, Simon (2015, Theorem 15.2.2) immediately gives
exp tr T + @5 det + 31)m22
+Y))
(+( x
(3)
See Appendix C.2 for a derivation of (3). We now apply Proposition 3.1 to this model to find the minimizing density for क. Let Σ- = YYT, and consider change of variables to .M- = (1+ 2nr0² - -Φ" T + IZ 7
2, x det L(O" -0r - (LL() 11 + O()).
I. Then, given the density of M before training as π, we can conclude that the density of M after training is
Applying Proposition 3.1, up to an (TT) term,
(M) x (det My-pm/2 exp( (opp-3pp'u(EEM----
(()) x det det .)) 0/2 10
no (T).
5

---

Models of Heavy-Tailed Mechanistic Universality
This expression need not satisfy the Master Model Ansatz. However, for losses of the form (D.1) (see Appendix D), ViL((Or, .) = ش (f(x)) when () = 0. so the spectrum of the Hessian is equivalent (up to zeros) to that of the NTK. Thus, the version of the Maserr Model Ansatz (5) (see Appendix D) also applies for the Hessian for small training loss.
explicitly and examined in greater detail in Appendix E. The proof of this theorem is deferred to Appendix E.2.
Theorem 4.1 (HTMP). Consider a sequence of matrices
MN obeying the high-temperature inverse- Wirseer ensem- ble (6) with x = (NN)). Assume ((N) == --</2-1 + for some constant y E (0, 1) N + The empirical spectral distributions of 2n(N) My converge to:
k/2
4. Heavy- Taay-TTailddae Spectral Behavior
(a) MP, (Marchenko-Pastur) if ((N) + oc; or
With the Master Model Ansatz (2) for trained feature matri- ces in hand, we can explore potential theoretical explana- tions for HT-MU that fit within this structure. We begin by outlining a proposed family of random matrices that satis- fies (2) and that exhibits spectral behavior (from Table 1) reminiscent of observations of HT-MU.
(b) HTMP Ty.e
(high-temperaturee Marchenko-Pusee if
<(N)
HTMP densities generalize the celebrated MP law, with an additional shape parameter, K, that permits interpolation between heavy-tailed spectra (N + 0+) and MP spectra (x - 20). We highlight the visual differences between MP and HTMP densities in Figure 2, with the key discrepancy lying in the tail behavior: HTMP spectral densities exhibit heavier tail behavior and have infinite support (even though the elements of HTMP need not have heavy-tailed behavior).
4.1. The HTMP Spectral Density
Obtaining a spectral density from the Master Model Ansatz: (2) requires diagonalization: for M a symmetric matrix, we can perform the change of variables M + QAQT, where Q is an orthogonal matrix of eigenvectors and A is a diagonal matrix of eigenvalues. The covariance matrix E can be removed from (2) by change of variables, so for now let Σ 7 (we return the general E scenario later in Section 4.3). With the effect of D removed, only the influence of the density π of feature matrices M at initialization remains to be characterized, to completely determine the spectral density from (2).
4.2. Structured Feature Matrices
Our next objective is to investigate the nature of the param- eter x in the HTMP distribution and motivate the use of the high-temperature beta-ensemble (6) in the context of the Master Model (2). To do so, we shall first consider the joint density of eigenvalues when π is uniform over different structured matrix classes. Our motivating feature matrix is the NTK matrix, defined as an N x N matrix comprised of n X n blocks, each of size m x m. so N = mn. We consider cases where N is uniform over: (a) Diagonal matrices; (b) Commuting block diagooaa matrices, where each m x m block on the diagonal commutes with each other, and all other blocks are zero; (c) Symmetric block-diagonal matri- ces, where each m x m block on the diagonal is an arbitrary symmetric -ositive-definite matrix, and all other blocks are zero; (d) Kronecker-like matrices, where the eigenvector ma- trix Q is a Kronecker product Q @ Q2. where Q1 E mm and Q2 E Ra*"; and (e) arbitrary (symmetric) matrices. Activation and NTK matrices for neural networks have been hypothesized exhibit Kronecker-like structures, as in case (d) (Martens & Grosse, 2015). In Appendix F, we derive the joint eigenvalue density for each of these classes.
Although the joint density of eigenvalues can be compli- cated, depending on RT, we argue (see Section 4.2) that much of the behavior is captured by the extent of the eigenvalue repulsions. To isolate this effect, we consider the family of beta-ensembles (Dumitriu & Edelman, 2002), parameter- ized by 0 S K S N, for the joint density of eigenvalues:
(A AND II
V(A))
II 1........
AA,
(6)
For (6) to follow (2), we must take V(A) = --A Beta-ensembles are well-studied objects in RMT, but they are typically considered theoretical curiosities, rather than "physical" models of matrices (Dumitriu & Edelman, 2002; Forrester, 2010). The 1/N "high temperature" scaling has also heen examined (Forrester& & Mazzuca, 2021), but with- out application. We argue (in Section 4.2) that x determines the "degree of randomness" in the matrix model, or con- versely, the rigidity of the matrix structure.
The key observation, noted in Remark F.6, is that the joint densities for (a)-() all exhibit absolute differences of eigen- values, |\;- - Ayl. differing mainly in how many such terms are multiplied together. This variation reflects differing degrees of eigenvalue repulsion, consequence of the re- duced randomness in the eigenvectors, as formalized in the ceweneeecoo-eigevalle identity (Denton et al., 2022, Theo- rem 1). Because structured matrices lack uniformly random eigenvectors, their joint eigenvalue densities are not sym- metric, complicating the determination of limiting spectral
Fortunately, Dung & Duy (2021) have derived the spec- tral density for high-temperature beta-ensembles (6). We use their result to introduce a novel RMT class of high- temperature Marchenko-Pase Pasttr (HTMP) densities, param- eterized by an aspect ratio ? and a structure parameter K, in Theorem 4. 1.1. The HTMP class of densities is defined
6

---

Models of Heavy-Tailed Mechanistic Universality
Structure
varying function L(z), and Eir(MF+ < +00. Then the spectral measure paay of My satisfies pMy +000) Ca L(z as 0c for constant Co.
(a) Diagonal
(b) Commuting block diagonal
(c) Symmetric block diagonal
(m
1)
mam
Proposition 4.2 follows from Kotodziejek & Szpojankowski (2022, Lemma 7.2). A similar argument shows that popula- tion covariance can also generate power laws for M
(d) Kronecker-like matrix
(c) Symmetric (no structure)
mn
Table 2. Behavior of K* across matrix structures for matrices with
While spectral distributions of many datasets, including CIFAR-10 and large language datasets, tend to exhibit power law behavior (Clauset et al., 2009; Zhang et al., 2023), Proposition 4.2 "predicts" that if the population covariance is the only mechanism by which heavy-tailed spectra occur, then the tail exponents should not differ between model architectures. This is an example of the power-law in, power- law out (PIPO) principle. More recent analyses have shown how other decisions, including architectural decisions, can alter the power law (Maloney et al., 2022). However, these results hold only for specialized models. One advantage of our approach is that the direct influence of the data can be separated, so our conclusions depend only on the resulting structure of the feature matrix, independently of the data.
n x n blocks of size m x m. Assuming n > m. the first term in each expression provides the dominant behavior.
density. To make these densities symmetric, one can ap- proximate them by symmetric polynomials: isolating the leading-order order behavior leads to (6). Therefore, by proposing (6) as a variational family of approximations the spectral density, for arbitrary joint density of eigenvalues g, the temperature x² can be identified by
argmin (gt.(9<<lg)) (g))
(7)
We consider the behavior of x² for each of the five differ- ent matrix structures described above. See Appendix G for details. In particular, in Proposition G.I, we prove that k" counts the number of eigenvalue repulsions in q, when it exhibits simple product form, as in cases (a), (c). (e). The more complex structures, (b) and (d), require numerical methods. In Algorithm I of Appendix G, an efficient numer- ical procedure for estimating R= is provided (satisfying a central limit theorem, see Proposition G.5). Using symbolic regression, the relationship between x* and m, n can be ascertained. Our findings are summarized in Table 2. We observe a direct correlation between the size of x* and how restrictive the matrix structures are, with the most restrictive (diagonal) and least restrictive (arbitrary) cases occupying the smallest and largest possible values of k*
4.4. Tail Behavior
We now put everything together to state our main theorem. Theorem 4.3. Let Mx denote the N x N Gram matrix of trained feature matrix obeying the Master Model (2). For PN. the empirical spectral distribution of MN and E, the label covariance matrix with spectral measure jus, we huve: (x(A) + (uz & ()) as N + c0, where X denotes free multiplicatiie convolution. Under the ensemble (6), the limiting density p(>) = -=> PPTTMP(>") if is finite; p(A) = 2 (X-1) if oo, where PMP and рHTMP are the probubility density functions of MP,, and HTMP,, defined in Theorems E.1 and E.2, respec- rively. Additionally,
(8)
4.3. Including Population Covariance
The most popular hypothesis behind the origin of heavy- tailed spectral behavior is the occurrence of this behavior in the d 3 We can understand this in terms of the Master Model Ansatz as follows. For a fixed prior π, let My be distributed according to the (6). Assuming a homogeneous prior, Ms. is equivalent in distribution to 21/2 Thus, the spectral density of My can be computed from Mj and E using Voiculescu's S-transform (Voieulescu, 1987). We then have the following proposition, which highlights how heavy-tailed spectra in My can arise from heavy-tailed spectra in the population covariance Σ, assum- ing that the feature matrix M, for isotropic data exhibits lighter-tailed spectra than the data. This true whenever regardless of whether (6) holds.
Power Law at oot For some constant C+ >> 0.
bounded support 3/2
if = 70,77 #
p(x) c4x c+r -
if = 20, 1 as J. - no.
otherwise.
Inverse Gamma Law at O: If = oo, the support of p is bounded away from zero; otherwise, for some 3_ > 0,
3_ as +
p(x) ~ xc - exp
Remark 4.4. The power law for the limiting density p con- tains a tail exponent that gets heavier as k decreases, i.c. as the structure of the underlying matrix becomes more rigid. Interpreting this as increasing implicit model bias, our findings are in line with conjectures of Martin & Mahoney (2021b) and Simsekli et al. (2019), claiming that heavier tails imply stronger model biases, and therefore better model quality and generalization ability. Importantly, the elements of these models need not have heavy-tailed behavior.
Ms =
Proposition 4.2. Assume that the spectral measure μi of satisfies p=((2, ~ 2 - [(2) as + oc, for slowly
'Even if this is less interesting from the perspective of HT-MU and this paper, we include this for completeness.
7

---

Models of Heavy-Tailed Mechanistic Universality
To the best of our knowledge, our HTMP model represents the first RMT ensemble that captures key empirical proper- ties of (strongly-correlate)) modern state-of-the-art neural networks (Martin & Mahoney, 2021b; Martin et al., 2021).
epoch T, then CN ~ cT a as T + oo (Agrawal et al., 2022). Velikanov & Yarotsky (2024) provide a comprehen- sive theoretical framework establishing power law rates of convergence across multiple optimizers (see also Velikanov & Yarotsky (2021)), under the assumption that the spectral density of J satisfies a power law near zero.
5. Applications
Here, we describe how our theory relates to several heavy- tailed results observed empirically for neural networks.
5.2. Optimizer Trajectories
Contrary to the popular belief that stochastic optimizers exhibit Gaussian fluctuations, Mandt et al. (2016); Simsekli et al. (2019) observed that the distribution of gradient norms of large neural networks during training exhibits a power law. This behavior manifests in heavy-tailed fluctuations during training, and it enables rapid escape from basins of poor-performing models in the loss landscape (Nguyen et al., 2019). We use the terms lower and upper power law to correspond to polynomial behavior in the distribution of gradient norms about zero and infinity, respectively:
5.1. Neural Scaling Laws
Perhaps the most famous appearances of power laws in deep learning are the neural scaling laws, which assert that the test loss at the end of training scales as power law with respect to both the number of parameters, d, and the size of the dataset, N (Hestness et al., 2017; Kaplan et al., 2020). In its modern form given by Hoffmmnn et al. (2022), using our notation, the test loss is observed to behave L~ Lo + A + N²: for some constants a, B, A, B, Lo > O. Since the total cost C' of training satisfies C' x Nd, for fixed computational budget, one can identify the optimal choices of N and d to achieve minimal test loss.
B
P(I | S x) C_ x.
as x - 0+
(10a)
P(I > x) C4. as co.
(10b)
We present a neural scaling law for ridge regression on the activation matrix satisfying the spectral density function in Theorem 4.1. Unlike previous scaling law work, instead of assuming power law in the dataset (c.g., done by Wei et al. (2022); Defilippis et al. (2024); Paquette et al. (2024); Lin et al. (2024)), we show how the scaling limit depends on the power law in the feature matrix L discussed in Section 3. The following proposition is proven in Appendix C.3.
Both behaviors have been observed in practice, with gen- eralization performance tied to the lower tail exponent a (Hodgkinson et al., 2022). There has been significant theo- retical justification for the upper power law (10b) in terms of the Kesten mechanism (Hodgkinson & Mahoney, 2021; Gurbuzbalaban al., 2021; 2022), but there has been little justification for the lower power law (10a). Our analysis can partially explain this too. The following proposition is proven in Appendix C.4*.
Proposition 5.1 (NEURAL SCALING LAW), Consider the activation matrix scenario in Section 3 with m = 1 and p X + R.d. Suppose that each label y; _ sp(x)),, E [n] and & satisfies the conditions of Theorem 4.1, with Σ = I and parameters k and q, For the solution vo argmin L(w,9) D), consider the generalization error
Proposition 5.2. Consider the loss function (4). Assuming the residuals fe (X) - Y are normally distributed and independent of an inverse Wishart-distributed NTK matrix J, the stochastic gradients L satisfy (10).
5.3. The 5+1 Phases of Training
L := E s.w. ~C( 21,,,,)),,)0)) -
The more difficult power laws to explain are those appear- ing in the weight matrices themselves (Martin & Mahoney, 2019; 2020; 2021b). Power law exponents in the spectrum of weight matrices are simple to compute and are strongly predictive of model performance (Martin et al., 2021; Yang et al., 2023; Zhou et al., 2023). In a sequence of papers, Mar- tin & Mahoney observed six classes of behaviors in trained weight matrices, with a smooth transition from random- like MP law to a heavy-tailed density, before experiencing a "rank collapse." " Excluding the rank collapse phase, the five primary phases are: (a) Random-Like: pure noise, modeled some spikes occurring outside the bulk of density; (c) Bulk+Spikes: spikes are dis- tinet and separate from the MP bulk; (d) Bulk-Decay. tails extend so that the support of the density is no longer finite;
Then, for p =t with E (0, 1), with high probability,
L 12+ - $) as +
We compare this result with previous work in Appendix H. Several works have outlined theories to explain the appear- ance of neural scaling laws (Bahri et al., 2024; Maloney et al., 2022). Each assumes the eigenvalues of a covariance matrix (either the final layer of activations or J) exhibit power law decay ' ~ ck- We remark that the shape of the spectrum of J after training is key, as the spectrum at initialization cannot predict neural scaling laws (Vyas et al., 2023; Bordelon et al., 2024b) for feature learning.
by a MP density; (b) Bleeding-
Similar to neural scaling laws is the appearance of power law decay in the rate of convergence to zero in the train- ing loss. That is, letting fr denote the model obtained at
'Proposition 5.2 for the beta-ensemble (6) remains an open problem, although we conjecture that it holds universally
8

---

Models of Heavy-Tailed Mechanistic Universality
1.0
K=
0.8
0.6
0.4
0.0
0.0
0.0
0.0
0.5
1.0
1.5
2.0
2.5
0.0
0.5
100
1.5
2.0
2.5
3
5
(a) Random-Like: batch size 1000
(b) Bleeding- Outt batch size 800
(c) Spik+Spikes:: batch size 250
0.8
0.10
0.8
5.5
1.9
n= 10
0.08
0.6
0.06
0.4
0.4
0.04
0.2
0.2
0.02
0.0
0.00
(d) Bulk-Decay: batch size 100
(c) Heavy-aaalllee Tailed: batch size 50
(1) Rank Collapse: batch size 5
Figure 2. The "S+1 phases of training" in weight matrices, estimated by Theorem 4.1. Compare with Figure 12 of Martin & Mahoney (2021b) (which is Figure 14 the technical report version their paper). All spectral densities (black) are compared to a corresponding MP density with an aspect ratio 7 = 0.3255. The red dashed lines are density functions in Theorem 4.1 with different K. The top row comprise cases with x = oc; the last row involves k= = 5.5, x 1.9 and x = 10-3 from left to right. See Section 5.3 for details.
and (e) Hlavy-Tailled: the tails become more heavy-tailed, exhibiting the behavior of a (possibly truncated) power law. The transition from (a) to (e) is also seen in Thamm et al. (2022), as well some non-uniformity in the eigenvectors, indicative of matrix structure. This smooth transition be- tween multiple phases is primary motivation of this work. We find that this behavior is displayed by a combination of a nontrivial covariance matrix to capture the spikes and the HTMP class with decreasing K. Indeed, we have the following, Proposition 5.3, proved in Appendix C.5.
structure and/or dynamics, we present a modeling frame- work for spectral densities of trained feature matrices, in- cluding activation matrices, Hessian matrices, and NTK matrices, culminating in Master Model Ansatz (2). Based this, we explored several factors contributing to the ori- gins of HT-MU, In line with the PIPO hypothesis, if the covariance of the data exhibits heavy-tailed spectra, trained feature matrices should be expected to adopt this behavior. Moreover, if a model exhibits bias on the distribution of eigenvectors, reducing eigenvector entropy, the joint density of eigenvalues can be approximated by a temperaumee Wiitret ensemble. The local eigenvalue spacings can be characterized with new parameter, K, and we pro- pose the use of the HTMP y.u model for the limiting spectral density. This distribution becomes more heavy-tailed as x decreases (more implicit bias, see Theorem 4. I). This model was used to explain recent examples of heavy-tailed phenom- ena: obtaining new scaling laws; explaining the presence of lower power laws in optimizer trajectories; and describing the five phases of training from Martin & Mahoney (2019).
Proposition 5.3 (WEIGHT MATRICES). Consider the acti- vation matrix scenario in Section 3. Let A = WTW, where W = argminy L(W, 4). Assume M क5 T satisfies (2) and (2MM) has joint eigenvalue density (6) with V As m. d + 0o, for n in Theorem 4.1, the empirical spectral distribution of 3
-- A converges to HTMP, 7,K-
In Figure? 2, we show that the HTMP family, HTMP, 7. provides good fits to observed spectral distributions across the five phases of training, compared with Martin & Ma- honey (2021b, Figure 2). In Figure 2, we plot the eigenvalue distributions of the fully connected layer of AleiAlexNet (a smaller AlexNet used by Zhang et al. (2021)) after training, using different batch sizes. In all the cases, we ensure the training losses are smaller than 0.01. Figure 2 indicates how training dynamics affect the heaviness of the tail, related to k, for the trained weight matrices.
We anticipate that the HTMP density can serve as valu- able starting point for investigating heavy-tailed spectral phenomena. Nevertheless, we emphasize that our primary objective was to establish a theoretical framework capable of elucidating the underlying causes of heavy-tailed spectral behavior (HT-MU). It is premature to precisely link model generalization to our hyperparameters, and / in (2) and k in (6). In future work, we can examine more closely this connection and extend the scaling law in Proposition 5. to more general cases.
6. Conclusions
Motivated by complementary lines of work deseribing heavy-tailed properties in state-of-the-ee neural networks

---

Models of Heave-Tailed Mechanistic Universality
Acknowledgements
Bak, P., Tang, C., and Wiesenfeld, K. Self-organized criti- cality: explanation of the 1/f noise. Physical Review Letters, 59(4);381, 1987.
The authors thank Dr. Chris van der Heide for reading our manuscript and providing valuable feedback. LH is supported by the Australian Research Council through a Discovery Early Career Researcher Award (DE240100144), ZW is supported by the National Science Foundation under Grant No. DMS-1928930, as resident at the Simons Laufer Mathematical Sciences Institute in Berkeley, California, in Spring 2025. MWM acknowledges partial support from DARPA, DOE, NSF, and ONR.
Bartlett, P. L., Foster, D. J., and Telgarsky, M. J. Spectrally- normalized margin bounds for neural networks. Advances in neural information processing systems, 30, 2017.
Bartlett, P. L.,., Montanari, A., and Rakhlin, A. Deep learning: a statistical viewpoint. Acta numerica, 30:87-201, 2021.
Bégin, L., Germain, P., Laviolette, F., and Roy, J.-F. PAC-
Bayesian bounds based on the Rényi divergence. In Arti- ficial Intelligence and Statistics, pp. 435 044, 2016.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be
Belkin, M., Hsu, D., Ma, S., and Mandal, S. Reconciling modern machine-learning practice and the classical bias- variance trade-off. Proceedings of the National Academy of Sciences, 116632)15849-1555 2019.
specifically highlighted here.
Bertschinger, N., Natschläger, T., and Legenstein. R. At the edge of chaos: Real-time computations and self- organized criticality in recurrent neural networks. Ad- vances in neural information processing systems, 17, 2004.
References
Absil, P.-A., Mahony, R., and Sepulchre, R. Optimization algorithms on matrix manifolds. Princeton University Press, 2008
Bordelon, B., Atanasov, A., and Pehlevan, C. A dy- namical model of neural scaling laws. arXiv preprint arXiv:2402.01092, 2024a.
Agrawal, K. K., Mondal, A. K., Ghosh, A., and Richards, B. a-ReQ: Assessing representation quality in self- supervised learning by measuring eigenspectrum decay. Advances in Neural Information Processing Systems, 35: 17626-17638, 2022.
Bordelon, B., Atanasov, A., and Pehlevan, C. How feature learning can improve neural scaling laws, arXiv preprint arXiv:2409. 17858, 2024b.
Ameli, S., van der Heide, C., Hodgkinson, L., Roosta, F., and Mahoney, M. W. Determinant estimation under mem- ory constraints and neural scaling laws. arXiv preprint arXiv:2503.00 04424, 2025.
Bouchaud, J. P. and Potters, M. Theory of Financial Risk and Derivative Pricing: From Statistical Physics to Risk Management. Cambridge University Press, 2003.
Anderson, G. W., Guionnet, A., and Zeitouni, O. An Intro- duction to Random Matrices. Cambridge Studies in Ad- vanced Mathematics. Cambridge University Press, 2009.
Boucheron, S., Lugosi, G., and Massart, P. Concentration Inequalities: Nonasymptotic Theory of Independence. Oxford University Press, 02 2013.
Arous, G. B. and Guionnet, A. The spectrum of heavy tailed random matrices. Communications in Mathemati- cal Physics, 278(3):715-75) 2008.
Cacoullos, T. On upper and lower bounds for the variance a function of random variable. The Annals of Proba- bility, 10(3):799--8)) 8099 1982.
Arous, G. B., Gheissari, R., Huang, J., and Jagannath, A. High-dimensional SGD aligns with emerging out- lier eigenspaces. In The Twelfth International Conference Learning Representutions, 2024.
Cao, Y. and Gu, Q. Generalization bounds of stochas- tic gradient descent for wide and deep neural networks. Advances in neural information processing systems, 32, 2019.
Atanasov, A., Bordelon, B., Zavatone- Veth, J. A., Paquette, C., and Pehlevan, C. Two-point deterministic equivalence for stochastic gradient dynamics in linear models. arXiv preprint ar
Chaudhari, P. and Soatto, S. Stochastic gradient descent performs variational inference, converges to limit cycles for deep networks. 2018 Information Theory and Applications Workshop (ITA), pp. 1-10. IEEE, 2018.
2025.
Bahri, Y., Dyer, E., Kaplan, J., Lee, J., and Sharma, U. Ex- plaining neural scaling laws. Proceedings of the National Academy of Sciences, (2((2))e2311187112) 2024.
Clauset, A., Shalizi, C. R., and Newman, M. E. Power- law distributions in empirical data. SIAM Review, 51(4): 661-703, 2009.
10

---

Models of Heavy-Tailed Mechanistic Universality
Cohen, J. M., Kaur, S., Li, Y., Kolter, J. Z.,, and Talwalkar, A. Gradient descent neural networks typically occurs at the edge of stability, arXiv preprint arXi::2103.0000 2021.
(eds.), Advances in Applied Mathematics, Modeling, and Computational Science. Springer, 2013.
Engel, A. and den Broeck, P. V. Statistical mechanics of learning. Cambridge University Press, New York, NY, USA, 2001.
Cranmer, M. Interpretable machine learning for science with PySR and SymbolicRegression, jl. arXiv preprint arXiv:2305.01 00.55.05882 2023.
Forrester, P. J. Log-gases and random matrices. Princeton University Press, 2010.
Dandi, Y. Stephan, L.., Krzakala, F., Loureiro, B., and Zde- borová, L. Universality laws for Gaussian mixtures in generalized linear models. Advances Neural Informa- tion Processing Systems, 36:54754-54768, 2023.
Forrester, P. J. and Mazzuca, G. The classical 3-ensembles with 3 proportional 1/N: from loop equations Dyson's disordered chain. Journal of Mathematical Physics, 62(7), 2021.
Defilippis, L., Loureiro, B., and Misiakiewicz, T. Dimension-free deterministic equivalents and scaling laws for random feature regression. The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
Fort, S., Dziugaite, G. K., Paul, M., Kharaghani, S., Roy, D. M., and Ganguli, S. Deep learning versus kemel learning: an empirical study of loss landscape geometry and the time evolution of the neural tangent kemel. CoRR, abs/2010. 151511151,, 10, 2020.
Denton, P., Parke, S., Tao, T., and Zhang, X. Eigenvectors from eigenvalues: A survey of a basic identity in linear algebra. Bulletin of the American Mathematical Society, 59(1):31-58, 2022.
Geiger, M., Spigler, S., d' Ascoli, S., Sagun, L.., Baity-Jesi, M., Biroli, G., and Wyart, M. Jamming transition as paradigm to understand the loss landscape of deep neural networks. Physical Review E, 100(1):012115 2019.
DIME. NIST Digital Library of Mathematical Functions. https:/ /tl::/ddii /dlmff....////////////////////// nist. gov/ / Release 1.2.4 of 2025- 03-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. Cohl, and M. A. McClain, eds.
Germain, P. Lacasse, A., Laviolette, F., and Marchand, M. AAC-Bayesian learning of linear classifiers. In Pro- ceedings of the 26th Annual International Conference on Machine Learning, pp. 353-360, 2009.
Germain, P., Bach, F., Lacoste, A., and Lacoste-Julien, PAC-Bayesian theory meets Bayesian inference. Ad- vances Neural Information Processing Systems, 29, 2016.
Donsker, M. D. and Varadhan, S. R. S. Asymptotic evalua- tion of certain Markov process expectations for large time. IV. Communications on Pure and Applied Mathematics, 36(2):183-25) 1983.
Goldt, S., Advani, M., Saxe, A. M., Krzakala, F., and Zde- borová, L. Dynamics of stochastic gradient descent for two-layer neural networks in the teacher-student setup. Advances in neural information processing systems, 32, 2019.
Dumitriu, I. and Edelman, A. Matrix models for beta sembles. Journal of Mathematical Physics, 43(11):58300 2002.
Dumitriu, and Edelman, Global spectrum fluctuations for the 8-hermite and B-laguerre ensembles via matrix models. Journal of Mathematical Physics, 47(6), 2006.
Golub, G. H. and Pereyra, V. The differentiation of pseudo- inverses and nonlinear least squares problems whose vari- ables separate. SIAM Journal on numerical analysis, 10 12):413-432, 1973.
Dung, T. H. and Duy, T. K. Beta laguerre ensembles in global regime. Osaka Journal of Mathematics, 58(2): 435-450, 2021.
Gurbuzbalaban, M., Simsekli, and Zhu, L. The heavy- tail phenomenon in SGD. In International Conference on Machine Learning, pp. 3964-3975, 2021.
Dziugaite, G. K. and Roy, D. M. Computing nonvacuous generalization bounds for deep (stochastic) neural net- works with many more parameters than training data. arXiv preprint arXiv: /70.///00./008 2017.
Gurbuzbalaban, M., Hu, Y., Simsekli, U.. Yuan, K., and Zhu, I. Heavy-tail phenomenon in decentralized sgd. arXiv preprint arXiv:2205. 06680, 2022.
Edelman, A. and Rao, N. R. Random matrix theory. Acta Numerica, 14:233-297, 2005.
Hanin, B. and Nica, M. Products of many large random matrices and gradients in deep neural networks. Com- munications in Mathematical Physics, 376(1):287-32) 2020.
Edelman, A. and Wang, Y. Random matrix theory and its innovative applications. In Melnik, and Kotsireas,

---

Models of Heavy-Tailed Mechanistic Universality
Hastie, T., Montanari, A., Rosset, S., and Tibshirani, R. Surprises in high-dimensional ridgeless least squares in- terpolation. Annals of statistics, 50(2):949, 2022.
Advances in Neural Information Processing Systems, 31, 2018.
Jiang, Y., Neyshabur, B., Mobahi, H., Krishnan, D., and Bengio, S. Fantastic generalization measures and where find them. In International Conference on Learning Representations, 2020.
He, K., Zhang, X., Ren, S., and Sun, J. Deep residual leam- ing for image recognition. In Proceedings of the IEEE
conference on computer vision and pattern recognition, pp. 770-778, 2016.
Kaplan, J., McCandlish, S., Henighan, T., Brown, T. B., Chess, B., Child, R., Gray, S., Radford, A., Wu, J., and Amodei, D. Scaling laws for neural language models.
Hestness, J., Narang, S., Ardalani, N., Diamos, G., Jun, H., Kianinejad, H., Patwary, M. M. A., Yang, Y., and Zhou, Y. Deep learning scaling is predictable, empirically, arXiv preprint arXiv:1712.004400 2017.
arXiv preprint arXiv:2001 08361, 2020.
Kesten, H. Random difference equations and renewal theory for products of random matrices. Acta Mathematica, 131 00:207-248, 1973.
Hodgkinson, L. and Mahoney, M. W. Multiplicative noise and heavy tails in stochastic optimization. In Interna- tional Conference on Machine Learning, pp. 2262-42744 2021.
Kimura, T. and Mazenc, E. A. The Schur expansion of characteristic polynomials and random matrices. arXiv preprint arXiv:21 2021.
Hodgkinson, I.,, Simsekli, U., Khanna, R., and Mahoney, M. W. Generalization bounds using lower tail exponents in stochastic optimizers. In International Conference on Machine Learning, pp. 8795,, 2022.
Kotodziejek, B. and Szpojankowski, K. A phase transition for tails of the free multiplicative convolution powers. Advances in Mathematics, 403:108398, 2022.
Hodgkinson, L., van der Heide, C., Salomone, R., Roosta, F,, and Mahoney, M. W. The interpolating information criterion for overparameterized models. arXiv preprint arXiv:2307.07785v1 2023.
Krizhevsky, A. Learning multiple layers of features from tiny images. Master's thesis, University of Toronto, 2009.
Krizhevsky, A., Sutskever, I., and Hinton, G. E. Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems, 25, 2012.
Hoffmann, J., Borgeaud, S., Mensch, A., Buchatskaya, E., Cai, T., Rutherford, E., Casas, D. d. L.. Hendricks, L. A., Welbl, J., Clark, A., et al. Training compute-optinal large language models. arXiv preprint arXiv:2203. 2022.
Li, Q. and Sompolinsky, H. Statistical mechanics of deep linear neural networks: The backpropagating kernel renor- malization. Physical Review X, 11(3):031059 2021.
Holzmüller, D. On the universality of the double descent peak in ridgeless regression. In International Conference on Learning Representations, 2020.
Li, Y., Zhang, H., and Lin, Q. On the asymptotic learning curves kernel ridge regression under power-law decay. Advances in Neural Information Processing Systems, 36: 49341 143664, 2023.
Hotelling, H. The generalization of Student's ratio. In Breakthroughs in statisties: Foundations and basic theory, pp. 54-65. Springer, 1992.
Liao, Z., Couillet, R., and Mahoney, M. W. A random matrix analysis random fourier features: beyond the gaussian kernel, precise phase transition, and the corresponding double descent. Advances in Neural Information Process- ing Systems, 33: 33:13939-13955 2020.
Hu, H. and Lu, Y. M. Universality laws for high- high--iigenne dimensional learning with random features. JEEE Transactions on Information Theory, 193--1993,,,,4,, 2022.
Huang, J. and Yau, H.-T. Dynamics of deep neural networks and neural tangent hierarchy. In International conference on machine learning, pp. 4542-4551, 2020.
Lin, L., Wu, J., Kakade, S. M., Bartlett, P., and Lee, J. D. Scaling laws in linear regression: Compute, parameters, and data. The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
Ismail, M. F., Letessier, J., and Valent, G. Linear birth and death models and associated laguerre and meixner polynomials. Journal of Approximation Theory, 55(3): 337-348, 1988.
Lotfi, S., Finzi, M., Kapoor, S., Potapczynski, A., Goldblum, M., and Wilson, A. G. PAC-Bayes compression bounds so tight that they can explain generalization. Advances in Neural Information Processing Systems, 35:31459- 31473, 2022.
Jacot, A., Gabriel, F., and Hongler, C. Neural tangent ker- nel: Convergence and generalization in neural networks.
12

---

Models of Heavy-Tailed Mechanistic Universality
Lu, H., Zhou, Y., Liu, S., Wang, Z., Mahoney, M. W., and Yang, Y. AlphaPruning: Using heavy-tailed self regular- ization theory for improved layer-wise pruning of large language models. In Thirty-eighth Conference on Neural Information Processing Systems, 2024.
Martin, C. H., Peng, T. S., and Mahoney, M. W. Predicting
trends in the quality of state-of-the-an neural networks without access to training testing data. Nature Commu- nications, 12(4122):1-13, 2021.
McAllester, D. A. Some PAC-Bayesian theorems. In Pro- ceedings of the eleventh annual conference on Computa- tional learning theory, pp. 230-234, 1998.
Lu, Y. M. and Yau, H.-T. An equivalence principle for the spectrum of random inner-produce kernel matrices with polynomial scalings. arXiv preprint 1rXiv:2205.06300 2022.
McCarthy, J. E. Random commuting matrices. arXiv preprint arXiv:2305.20029, 2023.
Mahoney, M. W. and Orecchia, L.. Implementing regulariza- tion implicitly via approximate cigenvector computation. In Proceedings of the 28th International Conference on Machine Learning, pp. 121-128, 2011.
Mei, S. and Montanari, A. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communicattinns on Pure and Applied Mathematics, 75(4):6677766, 2022.
Maloney, A., Roberts, D. A., and Sully, J. A solvable model of neural scaling laws. arXiv preprint arXiv:2210.16835 2022.
Menon, G. and Trogdon, T. Lectures on random matrix theory. unpublished notes (March 2018), 2015.
Mandt, S., Hoffman, M., and Blei, D. A variational anal- ysis of stochastic gradient algorithms. In International Conference on Machine Learning, pp. 354-3633 2016.
Mezard, M. and Montanari, A. Information, physics, and computation. Oxford University Press, 2009.
Mohamadi, M. A., Bac, W., and Sutherland, D. J. A fast, well-founded approximation to the empirical neural tan- gent kernel. In International Conference on Machine Learning, pp. 25021-50522201 250081, 2023.
Mardia, K. V., Kent, J. T., and Taylor, C. C. Multivariate analysis. John Wiley & Sons, 2024.
Martens, J. and Grosse, R. Optimizing neural networks with cronecker-factorce approximate curvature. In Interna- tional conference on machine learning, pp. 2408-2417, 2015.
Muennighoff, N., Rush, A., Barak, B., Le Scao, T., Tazi, N., Piktus, A., Pyysalo, S., Wolf, T., and Raffel, C. Scaling data-constrained language models. Advances in Neural Information Processing Systems, (6:50358-50377 2023,
Martin, C. H. and Mahoney, M. W. Rethinking generaliza- tion requires revisiting old ideas: statistical mechanics ap- proaches and complex learning behavior. arXiv preprint arXiv: 2017.
Nair, J., Wierman, A., and Zwart, B. The Fundamentals of Heavy Tails: Properties, Emergence, and Estimation. Cambridge Series in Statistical and Probabilistic Math- ematics. Cambridge University Press, Cambridge, UK, 2022.
Martin, C. H. and Mahoney, M. W. Traditional and heavy tailed self regularization in neural network models. In International Conference on Machine Learning, pp. 4284 4293, 2019.
Nakkiran, P., Kaplun, G., Bansal, Y., Yang, T., Barak, B., and Sutskever, I. Deep double descent: Where bigger models and more data hurt. Journal of Statistical Mechan- ics: Theory and Experiment, 2021(12)1220033 2021.
Martin, C. H. and Mahoney, M. W. Heavy-tailed universality predicts trends in test accuracies for very large pre-trained deep neural networks. In Proceedings of the 2020 SIAM International Conference on Data Mining, pp. 505-513. SIAM, 2020.
Nguyen, T. H., Simsekli, U., Gurbuzbalaban, M., and Richard, G. First exit time analysis of stochastic gradient descent under heavy-tailed gradient noise. Advances in neural information processing systems, 32, 2019.
Martin, C. H. and Mahoney, M. W. Post-mortem on deep learning contest: a simpson's 's paradox and the comple- mentary roles of scale metrics versus shape metrics. arXiv preprint arXiv:2 00734, 2021a.
Noumi, M. Macdonald Polynomials: Commuting Family of q- D-Difference Operators and Their Joint Eigenfunctions, volume 50. Springer Nature, 2023.
Martin, C. H. and Mahoney, M. W. Implicit self- regularization in deep neural networks: Evidence from random matrix theory and implications for learning. The Journal of Machine Learning Research, 2211)7747-7551, 2021b.
Paquette, E.,, Paquette, C.,, Xiao, L., and Pennington, J. 4+3 phases of compute-optimal neural scaling laws. In The eight Annual Conference on Neural Information Processing Systems, 2024.
13

---

Models of Heavy-Tailed Mechanistic Universality
Pennington, and Worah, Nonlinear random matrix theory for deep learning. Advances in neural information processing systems, 30, 2017.
Tao, T. and Vu, V. H. Random matrices: Universality of local eigenvalue statistics. Acta Mathematica, 206(1): 127-204, 2011.
Pillaud-Vivies L.,, Rudi, A., and Bach, F. Statistical op- timality of stochastic gradient descent on hard learning problems through multiple passes. Advances in Neural Information Processing Systems, 31, 2018.
Thamm, M., Staats, M., and Rosenow, B. Random matrix analysis of deep neural network weight matrices. Physical Review E, 106(5):054124) 2022.
Velikanov, M. and Yarotsky, D. Explicit loss asymptotics in the gradient descent training of neural networks. Ad- vances in Neural Information Processing Systems, 34: 2570-25822 2021.
Resnick, S. I. Heavy-Tail Phenomena: Probabilistic and Statistical Modeling. Springer-Verlag, 2007.
Rudner, T. G., Kapoor, S., Qiu, S., and Wilson, A. G. Function-space regularization in neural networks: A prob- abilistic perspective. In International Conference on Ma- chine Learning, pp.
Velikanov, M. and Yarotsky, D. Tight convergence rate bounds for optimization under power law spectral con- ditions. Journal of Machine Learning Research, 25(81): 1-78, 2024.
2023.
Seung, H. S., Sompolinsky, H.,, and Tishby, N. Statistical mechanics of learning from examples. Physical Review A, 4555):6006-609 1992.
Vladimirova, M., Arbel, J., and Mesejo, P. Bayesian neural networks become heavier-tailed with depth. In NeurIPS 201 2018-Thirty-secoon Conference on Neural Information Processing Systems, pp. 1-7, 2018.
Simon, B. Advanced complex analysis. American Mathe- matical Soc., 2015.
Voiculescu, D. Multiplication of certain non-commuting random variables. Journal of Operator Theory, pp. 223- 235, 1987.
Simonyan, K. and Zisserman, A. Very deep convolu- tional networks for large-scale image recognition. arXiv preprint 140.: 2014.
Vyas, N., Bansal, Y., and Nakkiran, P. Empirical limita- tions of the NTK for understanding scaling laws in deep learning. Transactions Machine Learning Research, 2023.
Simsekli, U. Sagun, L., and Gurbuzbalaban, M. A tail-index analysis stochastic gradient noise in deep neural net- works. In International Conference on Machine Learning, pp. 5827-5837, 2019.
Wang. Z., Engel, A., Sarwate, A. D., Dumitriu, I., and Chiang, T. Spectral evolution and invariance in linear- width neural networks. Advances in neural information processing systems, 36:20695-20725 2023.
Soltanolkotabi, M., Javanmard, A., and Lee, D. Theo- retical insights into the optimization landscape of over- parameterized shallow neural networks. IEEE Transac- fions Information Theory, 65(2):742-769, 2018.
Wang. Z., Wu, D., and Fan, Z. Nonlinear spiked covariance matrices and signal propagation in deep neural networks. In The Thirty Seventh Annual Conference Learning Theory, pp. 4891-4955 2024.
Sonthalia, R., Lok, J., and Rebrova, E. On regularization via early stopping for least squares regression. arXiv preprint arXiv:2406.00 2.XY:22006.04455, 2024.
Watkin, T. L. H., Rau, A., and Biehl, M. The statistical mechanics of learning rule. Rev. Mod. Phys., 65(2): 499-556, 1993.
Sornette, D. Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools. Springer- Verlag, Berlin, 2006.
Wei, A., Hu, W., and Steinhardt, More than a toy: Random matrix models predict how real-world neural representa- tions generalize. In International Conference on Machine Learning, pp.
Sorscher, B., Geirhos, R., Shekhar, S., Ganguli, S., and Mor- cOs, A. Beyond neural scaling laws: beating power law scaling via data pruning. Advances in Neural Information Processing Systems, 35:19523-19536, 2022.
2022.
Wilson, A. G. Deep learning is not so mysterious or different. arXiv preprint 2xXi::5503002113 2025.
Spigler, S., Geiger, M., d' Ascoli, S., Sagun, L., Biroli, G., and Wyart, M. A jamming transition from under- to over- parametrization affects generalization in deep learning. Journal of Physics A: Mathematical and Theoretical, 52 (47):474001, 2019.
Wilson, J., van der Heide, C., Hodgkinson, L., and Roosta, F. Uncertainty quantification with the empirical neural tangent kernel. arXiv preprint arXiv:2502.029 02000027007 2025.
14

---

Models of Heavy-Tailed Mechanistic Universality
Xie, Z., Tang, Q.-Y., Cai, Y., and Sun, M. On the power- law Hessian spectra in deep learning. arXiv preprint 22X1:2201113011, 2023.
Yang. Y., Theisen, R., Hodgkinson, L.., Gonzalez, E., Ramchandran, K., Martin, C. H., and Mahoney, M. W. Test accuracy generalization gap: Model selection in NLP without accessing training testing data. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 3011-302 2023.
Zavatone- Veth, J. A. and Pehlevan, C. Depth induces scale- averaging in overparameterized linear bayesian neural networks. In 202/ 55th Asilomar Conference on Signals, Systems, and Computers, pp. 600-607. IEEE, 2021.
Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O. Understanding deep learning (still) requires rethinking generalization. Communications of the ACM, 64(3):107- 115, 2021.
Zhang, L.-X. Central limit theorems of recursive stochastic algorithm with applications to adaptive designs. The Annals of Applied Probability, pp. 3630-3658, 2016.
Zhang, Y., Kang, B., Hooi, B., Yan, S., and Feng. J. Deep long-tailed learning: A survey. IEEE Transactions Puttern Analysis and Machine Intelligence, 2023.
Zhou, Y., Pang, T., Liu, K., Martin, C. H., Mahoney, M. W., and Yang. Y. Temperature balancing, layer-wise weight analysis, and neural network training. In Shirth-seventh Conference on Neural Information Processing Systems, 2023.
15

---

Models of Heave-Tailed Mechanistic Universality
Appendices
Contents
Experiment Details
16
A.1 Empirical Spectral Distributions of NTK Matrices
16
A.2 Empirical Spectral Distributions of Trained Weight Matrices
17
Stochastic Minimization Model Optimizes PAC-Bayes Bounds
18
Proofs of Assorted Results from the Main Text
19
C.1 Proof of Proposition 3.1
19
C.2 Proof of Equation (3): for Activation Matrices
20
C,3 Proof of Proposition 5.1
C.4 Proof of Proposition 5.2
22
C.5 Proof of Proposition 5.3
23
D NTK Matrix Densities using the Interpolating Information Criterion
23
E The Tiip--Tempeertrre Marchenko-Pastut (HTMP) Distribution
25
E.1 Beta Laguerre Ensembles for Two Regimes
25
E.2 Proof of Theorem 4.1
26
E.3 Proof of Theorem 4.3
27
E,4 Properties of High-Temperature Marchenko-Passe Pasttrr Distribution
28
F
Structured Matrix Models: Weyl Formulae and Joint Eigenvalue Densities
29
Variational Approximations with Beta-Ensembler
34
G.I Counting Eigenvalue Differences
34
G.2 General Numerical Method
36
G.3 Variational Approximations for Commuting Matrix Models
37
Comparison Neural Scaling Law for Linear Regression
38
A. Experiment Details
A.1. Empirical Spectral Distributions of NTK Matrices
In this section, describe the experimental details behind the results in Figure where plot the empirical spectral distributions of NTK Gram matrices. The (empirical) NTK is
kNTK(æ, y) (fe(() Dfe(y)
T E mm
where Dfe E ddx is the Jacobian of a function fe + y c Rm with respect to its parameters θ E Rdd. The corresponding Gram matrix over {xr] C X is the matrix J E mmaxma with m blocks JI, kNTK (1), This
16

---

Models of Heavy-Tailed Mechanistic Universality
Gram matrix is often enormous when Se is a deep neural network: in double precision, the empirical NTK Gram matrix for an ImageNet-1K IK classifier requires more memory than most international datacenters (Ameli et al., 2025). While some approximations exist to estimate certain quantities of this matrix (Mohamadi et al., 2023), none of these approximations are valid for the purposes of plotting the eigenspectrum. Instead, we consider models trained on a subsample of 1000 datapoints from the CIFAR-10 dataset with 10 classes (Krizhevsky, 2009). We remark that is important, for our purposes, to ensure that the model is only trained on data contained in the matrix, as this will inevitably impact the shape of the resulting spectrum in our analysis.
In Figure we consider three convolutional neural networks of moderate size: VGGI1 (9.2M parameters) (Simonyan & Zisserman, 2014); ResNet9 (4.8M parameters); and 8 8 (11.1M parameters) (He et al., 2016). The output layer of each model is altered from its usual ImageNet counterparts to be classified into ten classes. All networks were initialized with weights randomly chosen according to the standard He initialization scheme (He et al., 2016). The top half of Figure 1 plots the eigenspectrum of these networks at initialization over 1000 entries of the CIFAR-10 10 dataset. This is a similar scenario to that of Wang et al. (2023), but with real model architectures. The inverse Gamma law is plainly visible here across the spectrum; a maximum likelihood fit to a, 8 in the law p(>) x À -3/2 is performed, with the Kolmogoroy-Sminng statistics reported². Values of 3 are small, and the tail exponent o < 2. We remark that the corresponding weight matrices in this scenario exhibit spectral densities with the Marchenko-Passe
law.
Next, the networks are trained on this set of 1000 examples using 200 epochs of a cosine annealing learning rate schedule with a 200 epoch period, starting from learning rate of 0.05 with a batch size of 64, All models achieve near-zero loss under cross-entropy and an accuracy of 100%. The lower half of Figure presents the spectral densities after training. Here, there are two distinct "bumps", both of which exhibit inverse Gamma behavior. The right tail is not too different from that of the spectrum at initialization, so we draw special attention to the inverse Gamma behavior near zero. Once again, we perform maximum likelihood fits, with the Salmogorov-Smirnon statistics reporting good adherence to the inverse Gamma law.
A.2. Empirical Spectral Distributions of Trained Weight Matrices
In this section, we present the experimental details of Figure for 5+1 phases of trained weight matrices and their empirical spectral distributions. In these experiments, we train MiniAlexNet (Martin & Mahoney, 2019) for the CIFAR-11 10 classification task with different batch sizes. This MiniAlexNet is simplified version of AlexNet (Krizhevsky et al., 2012), which comprises two convolutional blocks-Conv(3 -> 64, 3 x 3, stride 2)+ReI.U+MaxPPe by two fully-connected layers (192 x 4 x 4 + 1000 + 10) with ReLU activations.
and Conv(64 - 192,
3 x MiLLUU--aaoPoolffolloo All weights are initialized via the parametrization (Jacot et al., 2018), i.e., the entries of trainable weights are i.i.d. (0,0 /.am.). We normalize each data point in CIFAR-10 to zero mean and unit variance per channel, and rescale it so that each flattened input vector has Euclidean norm v where d = 3 x 32 x 32 is the input dimension for each data point. We train this MiniAlexNe with SGD (momentum 0.9, weight decay 10-*), using a learning rate scaled 2/ //bacch_size for different batch sizes and early stopping once the training loss falls below 0.01 200 epochs. In Figure 2, we separately train MiniAlexNet with batch sizes 1000, 800, 250, 100, 50 and 5; and we repeat the experiments 3 times for the average. Post-training, Figure 2 shows histograms of the eigenvalues of WWT E IR 1000 1000 where W is the trained weight matrix in the first fully-connected layer with input dimension 192 x 4 x 4 and output dimension 1000, for different batch sizes. Notice that, because of the NTK parameterization at initialization, the empirical spectral distribution of WWT is close to the MP, with 7 = 1000/ (192 X 4 x 4). For large batch sizes, batch_size = 1000, after training, the empirical spectral distribution of W remains similar to a random case, as depicted in Figure 2(a). The red dashed line is computed using the density function of MP,, with y 1000/( (192 4 x 4), which is described by (E.3). This MP,, can also be interpreted oo. When the batch size decreases, we observe eigenvalue bleed out from the bulk spectrum and
as HTMP,,, T.M with the emergence of many spike eigenvalues in Figure 2(b) and (c). Despite this, the bulk spectrum remains relatively close to MP,,. y Further decreasing the batch size, Figure 2(d) and (c), the bulk spectrum of the trained weight is also changed into heavier tailed distribution which is close to HTMP,,. 7.x for some finite k. The red dashed lines in Figure 2(d) and (e) are computed by the density function of HTMP 7.M with k- = 5.5, 1.9 and y = 1000 Additionally, when the batch size is small (batch size- = 5), the trained weight matrix exhibits significant rank deficiency, corresponding to HTMP y. with approaching zero. In such scenarios, the model's ability to effectively learn the task is compromised due rer.laguultiizttion Hence, Figure 2(f) is not a particularly interesting case for further exploration.
x x x 4), defined in (E.5).
'Recall that the Kolmogorov-Smiroov statistic between two cumulative distribution functions F. and Fg is dRs (Fi, Fa) sup.cR [z)(z) - F2(z))
17

---

Models of Heavy-Tailed Mechanistic Universality
B. Stochastic Minimization Model Optimizes PAC-Bayes Bounds
In this section, we outline the connection between our stochastic minimization operator, sming of Equation (1), and the celebrated PAC-Bayesian bounds (McAllester, 1998), which have seen increasing popularity and success for explaining the success of deep learning (Wilson, 2025) and for providing non-vacuous generalization bounds (Lotfi al., 2022). In particular, will show that the stochastic minimization operator applied to training losses (as Monte-Carlo estimators of the lest loss) minimizes, by definition, PAC-Bayesian upper bound on the test loss. Our discussion mostly follows expositions by Germain et al. (2016). For the sake of brevity, instead of the two-part loss function, L( ②) from Section 3, we consider arbitrary loss functions of the form L(O). This is no less general: we can once again split the parameter space into two parts and apply our arguments here to each part independently.
Donsker-Varadhan The foundation of modern PAC-Bayes proofs, as first outlined in Bégin al. (2016), the Donsker Varadhan variational formula, first introduced in Donsker & Varadhan (1983, Eqn. (1.7)). For completeness, we state both the formula and its proof in Lemma B.I.
Lemma B.1 E and f: R.
VARIATIONAL FORMULA), For any probability measure ye measurable space sup (E)) - dkl(v|lu)) = log E v
for el E (p),, the probability measure that achieves this equality satisfies
dv- du (0) x exp(s(O)),
+ E
(B.1)
Proof. For brevity, let μ = Te. Let denote an arbitrary probability measure that is absolutely continuous with respect to μ.. By the Lebesgue decomposition, μ Mae + Msing where Mre < and Msing 1 V. Then, the Radon-Nikodym derivative du du v-almost everywhere, and deu.
dKL (v|lm) (e) log v())
-Ev log v()) dkI.(ullmn) (=c(O)
((0)
Furthermore, by Jensen's 's inequality,
e (()) dd
dv < log ey (())
(() disc (0) = (xp((((())) < Ep(e) ex)(f(9)).
Ex)
de
At the same time,
(())) dpac dy
EV(
=Ee()) dKI(vIlHme) = [e(o) dKL(((1)))
and so
sup (Bm dKI.(vllx)) S Ey(e) exp(
ven
now only remains show that the measure satisfying (B.1) saturates this inequality. Letting 2 (e cxp()(e)), so that de = 2--1 =eex((((()) there is du
dvv* du- (e) - llgg 2) Ey (e) - log 2.
dKI. dkI.(="lln) = Ey. (0) log
du
Hence,
(e((e) dxL(v"Il)) = log 2 = log E e(0))
as required.
Now we return to the case where y and are both absolutely continuous measures with respect to Lebesgue measure, and are turn represented by probability densities on C Rd, and we let = (9), in line with the setup from Section 3. Most losses to be minimized are training losses of the form T., (9) = n -1 EE ((O, 3z, V:), where each (2i, U:) is an
18

---

Models of Heavy-Tailed Mechanistic Universality
input-output pair. Training losses random; a Monte-Carlo approximation of the test loss ((O) defined as the expectation (e). Let and consider the function (n(@) / (T((@) I. (+)). Then, by Lemma B.1,
E (n))) S dRL. (xk )e) + log Sn,rI
where Sn,r is defined by
T -L(1)))x((9))
(O) (E)))) exp
( F= are independent and identically distributed, then subgaussian concentration inequalities such as the Bernstein inequality (Boucheron al., 2013, Theorem 2.10) imply that RR,,, is bounded as -> oc. particular, we may assume that
K2 = sup ES,,T < +00.
n
From Markov's inequality, for 0 < < 1, (P > K+/6)= S P(Gn,r > ESn,T/0) S 8, and so with probability at least
1 1-6, 6,
Kr
Ex(e) S dkt((alTT)) + log
Recalling the definition of fa reveals the (simplified) PAC-Bayes bound, as presented in Germain et al. (2009, Theorem 2.1). Theorem B.2 (PAC-BAYES THEOREM (Germain et al., 2009)). For any probability density q absolutely continuous with respect to me on C IRd, assuming (ri, are idd, then with probability at least 5,
T (())) < Ln (0) (gll=T)
I n
K 8
log
(B.2)
As the left-hand-side of (B.2) corresponds with the test loss, we are interested in considering densities q that assign probabilities to model parameters ④ that have been fitted to the data. The density me acts the prior in the PAC-Bayes theory, and it assigns probability mass to models that are believed to be more or less likely to exhibit small test loss, prior to looking at data. This is to be contrasted with uniform PAC bounds, which typically treat all models as equally likely. A key feature of the prior is that it is data-independent, but otherwise it is arbitrary. It is this feature that has allowed (B.2) to achieve the tightest generalization bounds to date (Lotfi et al., 2022).
now also immediately clear, by definition, that the density In which minimizes the bound of the test error on the right-hand side of (B.2) is precisely given by the argsmin operator:
To. argsmnnn L,, (e),
On (
where Tn
C. Proofs of Assorted Results from the Main Text
In this section, we present proofs of several results from the main text. Below, we first prove Proposition 3.1 and subsequently complete the proof of Equation (3) for activation matrices introduced in Section 3.2. We defer the analysis of NTK matrices and the proof of Equation (5) to Appendix for further details. Subsequently, we present the proofs of Propositions 5.1, 5.3, and 5.3. Our random matrix results, as outlined in Theorem 4.1 and Theorem 4.3, are presented Section E to ensure clarity and self-containedness.
C.1. Proof of Proposition 3.1
The primary tool in the proof of Proposition 3.1 a corollary of the Donsker- Donsker-Varadhan variational formula in Lemma B.1, which we now state in the notation of Section 3.
Corollary C.I. For C Rd and f : = + R such that (()) = exp(( is integrable over θ E, argsmin 78.7 f(0) x (().
Proof. Letting denote the set of probability densities that are absolutely continuous with respect to Lebesgue measure
19

---

Models of Heavy-Tailed Mechanistic Universality
on E,
me.7
argsmin ʃ(0) argmin
( [e))) + =dki.(ul|TT)) q(S)EP
( (-1 dKL(al1Te)
argmax
(
Applying Lemma B.I, it follows that the minimizing density is proportional to q()).
Let be as in Corollary C.I. and let Z, = Be-xe (x)(= so that ((() = 2 (()) is a probability density. As an immediate consequence of Lemma B.I, since dkt (Q] = - - BB ((O) log 2,, it follows that
me. smin f(0) = -T -T- log 27.
(C.1)
Proposition 3,1 now follows immediately from Corollary C.I and (C.1): for (+) as defined in Proposition 3.1,
argsmin smin L(O, exp +
Te-7
Te.T
(- log 2,
TT. (4( (0)) / TT((1) (T).
C.2. Proof of Equation (3): for Activation Matrices
Now we prove Equation (3). Recall that
L(W, ①)
|WW- 112 + wl
where μ > 0. Our objective is to first compute
2, ()) Ew (W exp -- - Y - [1WIp
For the sake of brevity, we ignore all constants of proportionality that do not depend on 4. First, observe that
(D) /
exp | - r12 " [Iw 202 |W] dw,
and so 2,(0) is proportional to the marginal likelihood p(Y|#) p(YY: () corresponding to the likelihood
(2 mon/2 exp (-1 l y
pY|W, +)
of the model Y 4W + Şe, and the prior
(-
1
[1
p(W) x exp
20
implying W ~ ~x(00, 821), where
2/ +
1
0²
7 0²
2
Let = vec(Y) and w: = vec(W), recalling that y = (I @ k)w normal also, with Eyy|0 = and
I under the model likelihood. Since w is normal, y
Cov(y) = Cov((1 & ())) +
= =T @ क + +
8 ( +
20

---

Models of Heavy-Tailed Mechanistic Universality
Therefore,
exp 3 (1 (o + "1)]
p(u|0)= (2 (2m)
det(7 @ (ō2TTTT T + /))
Since tr(A B) vec(A) vec(B), get
@ (o2 - y (eve((@22pp) T + +11)
= tr((Y (02 + ) 'Y)
In terms of Y,
(+) x p(Y exp (= (0200 + -1)1
det("=m T + /1/m
C.3. Proof of Proposition 5.1
Recall the solution the ridge regression problem as
é = ( Ty E R
(C.2)
and, by definition, we know that
1 (Φ ① A1)-1a1 (w.,,
for any 0. Thus, the training loss
Lo (v(x:) [w.)] w.
\3w. (T ( - + XT) -100 w..
(C.3)
Here, we do not consider label noise: each label is the form of yı =p(x;) w.. Based on the definition of the generalization error defined in Proposition 5.1, we have
L = [(((((x)) w (xxx) Tw.)
(C.4)
(u w.) " w-)
(C.5)
= w. ((0 1b + 21) 1 + - 1)E(( + x -l),,
(C.6)
= Aoo Φ -
(C.7)
where we denote that Σ := [[z(x)(((x) J. Suppose that w. ~ (0, 1)) Then, we can apply Hanson- -rright inequality get
X2 L d tr T Φ + A1) - + o(1),
with high probability as d + oc. Based on our assumptions in Proposition 5.1, we set be the identity matrix, namely the isotropic case, (For general E, we would need apply Theorem 4.3 and deterministic equivalent statement for tr(T t + A2,, for any deterministic matrix A dRdd However, currently, such deterministic equivalence for inverse Wishart matrices is unknown in RMT. We leave this general statement for future work.)
Since we assume that & T d has a limiting eigenvalue distribution with density functions described in Theorem 4.1, we can claim that
L - X (-x) =: f()),
almost surely, d + where Sy.,, (z) is the derivative of Sp (z), defined in Proposition E.3, with respect to z. Here Sp.... (2) is the Stieltjes transform of HTMP,,, 7.x introduced in Theorem 4.1. For more details, see Appendix E. From (E.7), we know that
U(s/2+1,2 24
8 Spm. (- K
A> 0.
(C.8)
x + 2 U(n/2,
8.
2'
21

---

Models of Heavy-Tailed Mechanistic Universality
Now want to study the asymptotic limit of the function f(A) À Notice that X 0 represents the ridgeless regression; and, in this case, the training loss will converge zero, because of (C.3). Hence, considering the case when the estimator interpolates all the training datasets.
From DLMF, (13.3.22), recall that
d -a) 1,2).
dz
In general, from DLMF, (13.2.16-22), also know that
[(b 1) _1= + +0 T(a)
U(a,b,=)
Nib 2, b 4
(C.9)
2(a,2,2) = r(a)* 2 + 0(n 2).
(C.10)
b() r 10) 11 T(1 -b)
T(a - b-1) +0
15 RD 2,
(C.11)
T(a)
1
1,)
(In (a) + + 0((Imn)),
(C.12)
r(a)
[' - b)
b()
(2
0< Rb < 1,
(C.13)
[(a b+
U,),0,=)
(z ln=),
(C.14)
[(a + 1) [(1 b)
U(a,b,2) = l' 1)
+ (2),
Rb 20
(C.15)
Now for simplicity, (C.8), we define that
R\
27
k/2 1, b= x/2
2
c k/2, = K
K
+1+
27
2
Thus, we now have
A(A) 22 cU 1, 1;z)U(a, aU 1,b 1;)++;
U(c, d: 2)²
aa b;) (a + 1, ( d;²) U(c, =)²
Let us consider the case when and d > Hence, using the asymptotic results (C.5)-(C.15,, we can get that
f(A) x d+1
as z + 0+. Notice that here b, d are constants and E (0. 1). Therefore, when assuming that À = n nf, we can conclude the result in Proposition 5.1. For more comparison of Proposition 5.1 with previous literature, we refer to Appendix H.
C.4. Proof of Proposition 5.2
The stochastic gradients V Ln are given by 6L = oX). Z. where Z = (e) - are assumed to be independent and standard normal. Consequently,
2 ZT T D (X) = JZ,
where J has the inverse- Wishart distribution. Due to a result of Hotelling (1992), d-N+1 V is F-distributed with (N, N + 1) degrees of freedom, and therefore satisfies (10).
dN
22

---

Models of Heavy-Tailed Mechanistic Universality
C.5. Proof of Proposition 5.3
Consider the activation matrix model in Section 3. By Golub & Pereyra (1973, Theorem 1), with the pseudo inverse, trained weights are given by W ΦtY, so A = d' W W d' 'Y ( --. Since the (nonzero) eigenvalues of the product of two matrices AB are always the same those of BA, assuming Σ = YY" is full rank, the spectrum of A is equivalent to that of B Σ M where M = ΦΦ Under the conditions of Theorem 4.1, M satisfies the Master Model (2) and so
(M) x (det M) exp(-Str(EMM"))x((
Performing the change of variables M I- we find that
(BB) x (det B) exp( exp(-Btu(B~"))x(11
and since assert that n (EEB has Weyl formula that aligns with (6) for V the eigenvalues of B satisfy
N
q(1 AN)
3/2:
II
A - 1,/
1=1
N
From Theorem 4.1, with B = B(m, d) and
7 --x/2-I-
k/2
the sequence of empirical spectral distributions satisfy
B
B(m,))-1)
m,d-000
spec
HTMP 7
a k22
Since spec(B(m, d)- spec(A(m, d)), the result follows.
D. NTK Matrix Densities using the Interpolating Information Criterion
Motivated by the overparameterization (d > vi) many deep neural network models, Hodgkinson et al. (2023) developed approximations to the marginal likelihood, for very general models, that are effective in the interpolating regime. These approximations led to the development of an interpolating information criterion to judge model quality in the interpolating regime. In this section, we will show how their techniques combined with Proposition 3.1 reveal a precise family of densities for the NTK matrices of interpolating models. We will also show how the NTK approximation discussed in Section 3 compares and its implications on implicit regularization.
Suppose that (+) x e (+) and let O argming R and (Q) = (o) Rs (Oa). Consider the neural tangent kernel matrix, J(D) E m which satisfies J(D);s iij Df(z:) Dz(z)), where Df dx is the Jacobian of model, fe.* X - y C R=, in its parameters θ E Rd. If V represents the manifold Hessian (Absil et al., 2008), then the marginal Gibbs likelihood for the loss
L(O, Φ) (1 Yll
(D.1)
for inputs X and outputs Y. satisfies (Hodgkinson et al., 2023, Theorem 1)
Ro(00)
det V2 (Gn) det V Ro((*) [1 + O(T + p)),
(D) x det J(
(D.2)
where O* is the interpolator
argmin R4 (e) subject to
L(O, b) = 0.
(D.3)
Using (D.2) with Proposition 3.1 constitutes a precise formula for minimizing the density: up to an O(T + p) factor,
det V Ro (o det V (t.
P/2
exp ( (__Ro((e) T (2).
( (())
(D.4)
While we cannot elucidate the spectral density from (D.4), it provides a more accurate picture of the complexities with ascertaining the densities of trained feature matrices at scale. Indeed, we can use (D.4) to understand our model for the NTK feature matrix (5) without appealing to linear approximations.
23

---

Models of Heavy-Tailed Mechanistic Universality
First, we ask the question: if O* should yield robust predictor, what is good choice of R&? For the sake of brevity, we drop the dependence on the features + and write Rs as simply R. In the overparameterizze regime, the test loss is primarily dictated by the variance in model predictions with respect to label noise (Holzmülleer, 2020). For Lipschitz loss function, this variance depends on that of the estimator @ itself. Consequently, to ensure that solution to (D.3) exhibits small test loss, the regularizer R should be chosen to control the variance of the induced estimator. This leads us to the following Definition D.I. Here, the variance of a random vector is considered as the sum of the variances over each coordinate, i.c., VarX E: VarX;,
Definition D.I (IDEAL REGULARIZER). An ideal regularizer (for minimizing test error) is a choice of R(() that controls the induced label noise of the corresponding estimator:
R(O) 2 VaryO'(Y).
Although the precise form of an ideal regularizer may be quite complicated, simple lower bound can be derived in the scenario where the label noise has finite covariance.
Lemma D.2. If J is full rank, then VaryO*(Y) 2 tr(EA-") + 01))) where Σ Cov(YY.
Proof. First, performing Taylor series expansion in Y, we find that
Vary (Y) = -Y>))\>>>'00>>>>>> [IDy0 + ((
(D.5)
where 1 is the symmetric square root of Σ. Letting fe- denote the underlying model, since fe- in the interpolating regime, taking derivatives of both sides in Y reveals that DF D and D is a right-inverse of DF. Consequently, Dy Dy-'r 11/2 is also right-inverse of E-1/2 DF, and
| 2 |(E-1/2DF)
=tr DEDE
=Ir(EJ-),,
where the second equality follows from the identity | tr((X'X 11) Altogether, this implies that
Vary (r(Y) 2 tr(EAJ-1) - o (|[E)12)
as required.
Remark D.3. In the case where Y is normally distributed, (Cacoullos, 1982, Proposition 3.7) implics that (D.5) can be replaced by
Vary |
and so the limit o( ([]]) can be dropped.
Choosing R(O) tr (r(E)ー |DF to saturate these lower bounds, (D.4) becomes
det (2R del V R(O*)) exp (-Lm(E)(O)))) T (T))
P/2
(()) x (()) 0/2
(D.6)
This is equivalent to the density (5) obtained under the linearized model except for the second-order factor (det (²)()))/ det (00) Formally speaking, it is unclear what the influence of this factor for general make some conjectures. For (D.6) coincide with (5) under the change of variables &b + J(D),
models, however, it is important for ((On) and det Vh, R(@*) to be minimally dependent on J()). In the former case, we could that the mutual information between R(Op) and R(O*) is small as consequence of the training procedure. This has some empirical merit according to Fort et al. (2020). The latter object is more complex, but depends primarily on the null space of DF(O*) and the higher-order derivatives of F. Again, we may conjecture that the mutual information between these quantities and the eigenspectrum of J(D) is relatively minimal, although such claims require further study.
24

---

Models of Heavy-Tailed Mechanistic Universality
E. The HighhTemperrture Marchenko-Pastie (HTMP) Distribution
In this section, we present more background on the beta ensemble RMT, we then prove Theorems 4.1 and 4.3, and we then present additional properties of the HTMP distribution introduced Theorem 4.1.
E.1. Beta Laguerre Ensembles for Two Regimes
In RMT, for N 2 and B > 0, the eeta-Hermite (Gaussian) ensemble of size x N is defined by its joint eigenvalue density function:
PH(AI..... AN) x E X
ISKJEN
This reduces to the classical random matrix ensembles: GOE, GUE, and GSE, when 3 = 1, 2, respectively. We refer to Chapter of Anderson et al. (2009) for more details on beta ensemble matrices. In our setting, replacing the quadratic confinement above with a Wishart-type matrix yields the beta-Laguerre ensemble (Dumitriu & Edelman, 2002).
We define a random matrix Ay with a betu-Laguerre ensemble of size N N via its joint eigenvalue density function:
A()
x
II iml
xpx 6
-
II
14 Xo1
(E.1)
i.j=1.....N s-]
for some / = B(N)= > O and any N. D E N. Denote 3 = K/N for some & >> 0. In this section, we present RMT for the beta-Laguerre ensemble for two different regimes, based on different limits of B, or equivalently, the limits of k. We clarify that the parameter (3 the beta-Laguerre ensemble, as defined in (E.1), distinet from the / encountered in our Master Model Ansatz, presented (3.2). With slight ambiguity, the term employed in this section for a distinct purpose compared to its usage in the main text.
We let D, N - oo with D 2 N and N/D-72 10.
Then, (E.1) is asymptotically equivalent to
N
(x ( E A
a(A AN) & II 2m1
II
1>g - 1//NN
(E.2)
i.f=1.....N
We will prove Theorem 4.1 based on the joint density function (E.2). Here, we focus on the beta-Laguerre ensemble An and its joint probability density function (E.1), with 3 ((N) and D D(N) 2 N. For each value of N, let μu denote the empirical spectral distribution of Av Specifically, pA the measure defined as μuN = N-1 1 8x((As)-
We now present the limits of the empirical spectral distribution of AN, for two distinct regimes of ((N) or, equivalently, x(N)= B N. These have been summarized by Dung & Duy (2021). Later, will use these results to prove Theorem 4.1.
Firstly, we review the asymptotic results when N ((N) + co. In this case, the Marchenko-Pasee Pastur distribution appears in the limiting spectral distribution. The Marchenko-Pasee distribution MP,, with parameter y E (0, 1) is absolutely continuous on (0, oo) with finite support only on the interval I [7- where Y1 = (1 + 1/2 The corresponding probability density function is given by
1 V(7
xx 7
(x(x)=
x E In.
(F.3)
2T
Theorem E.1 (MARCHENKO-PASEE PASSTRR THEOREM). If (NN) = S((N-)) as N - so and N/D -> (0, 1), then the empirical spectral distributions μN of beta- ensemble Ax, defined by (E.1), converge weakly to the Marchenko- Pastur distribution MP,, whose probubility density function is pr (x) with parameter ", defined by (E.3).
This Marchenko-Pastur Theorem has been established by Dumitriu & Edelman (2006), for any fixed value of B. Subsequently, Dung & Duy (2021) extended this theorem to encompass general (NN), as long as the limit of N B(N) approaches intinity.
Secondly, we consider the high temperature case where 3(N) ~ - for some fixed x: > 0. Let (a, b, 2) denote the Tricomi confluent hypergeometric function (DLMF, 13.2.42), defined for negative arguments when b is not a negative
25

---

Models of Heavy-Tailed Mechanistic Universality
integer by
T(1 - b) Ula,b, b,, ((a - b+ Fi -2) +
[ T(b-1), 1) 21 F Fi b 1,,,,,, 2> 0, ((a)
(B.4)
where I is the Gamma function and Pf is the generalized hypergeometric function. The integer taking appropriate limits. Below, we present the result of the beta- Leta-Laguerre ensemble for the high-temperature case, proved by Dung & Duy (2021) using moment methods.
be obtained by
Theorem E.2. Assume that x(N) 3(N)/N - - E (0, 0o) as N - ac and N/D + E (0, The empirical spectral distributions Pa of beta-Laguere ensemble AN defined by (E.1), converge weakly to probability distribution on (0, oo) with probability density function Prys defined as
(o)(())= 2y r(k/2+ + ()0/2H) (UKK/2,-
27
K2
1-
e
27
x 2 0.
(E.5)
x
+1+
1+
27
We call this limiting spectrum the high-temperaturee Marchenko-Pastur (HTMP) distribution, denoted by HTMP 7. with parameters K, > 0.
E.2. Proof of Theorem 4.1
To prove Theorem 4.1, let My be an N N matrix possessing the high-temperature inverse-Wishart ensemble defined in (6), that is.
N
AN) &
A7
9A
II
A - X3 MM/
1J=1...,,N
We focus the matrix My whose eigenvalues reciprocals of those of MN. By change of variables A + A, the joint eigenvalue density of My denoted by Tn. satisfies
N
6/N
1
The ...........
o
2
AA
X
Furthermore, observe that
R/N
A - A
K/N
A
X;dj
kj=1...,N K
K
X K+K/NN
II 1N M/N
i.J1I...,N
and
N kNNNN
The (A (
An) x
e -BAe
II I.....N
1y - A//NN
Another change of variables \;
2x8 di shows that the matrix VN
2nd M has joint eigenvalue density
( AN) x 2 c
A
1 - M
(E.6)
Jill.. N i<j
Notice that 2 2-K- =5 1) 1 since denote
K/2
K/2
26

---

Models of Heavy-Tailed Mechanistic Universality
which turns the joint eigenvalue density of Vy into
N
q(A AN) &
II
2 11-1)+++
27
II 1.J=....N
P3 - A,
which is asymptotically equivalent to (E.2) as N + cc. Hence, is also asymptotically equivalent to (E.1), where we take 3(N) <(V)/N and N/D + 7 E (0, in (E.1). Notice that y (0, 1) indicates that we need assume a 1. Applying Theorems E.1 and E.2 shows convergence in distribution for two different regimes of ((N). Therefore, can conclude that
when x(N) + oc, the empirical spectrum distribution of Vi converges weakly to MP,;: and
when x(N) > k. the empirical spectrum distribution of Vy converges weakly to HTMP 7.M
This completes the proof of Theorem 4.1.
E.3. Proof of Theorem 4.3
The proof of (8) in Theorem 4.3 directly based on the Master Model Ansatz (2) and free probability theory. For more details on free probability theory, the definitions of multiplicative free convolution and asymptotic freeness, and their connections with RMT, see Chapter 5 of (Anderson et al., 2009). Precisely, to prove (8), we can apply an independent N X N unitary matrix UN following the Haar measure to Σ. Then, Corollary 5,4.11 in (Anderson et al., 2009) proves that MN and EN x are asymptotically free. Thus, we can conclude that the limiting eigenvalue distribution is the multiplicative free convolution μr Ø The precise definition of the multiplicative free convolution "Ø" of two probability distributions can be found in Definition 5.3.28 of Anderson et al. (2009).
Hence, when E is identity and MN follows the high-temperature inverse-Wist Wishart ensemble (6), can apply Theorem 4.1 to conclude that empirical spectrum density PN (A) p(A) weakly, where the limiting density p(A)= = PHTTR is finite, and p(A) A PM (A-1) )) if k co. Here PuP and PurSs are the probability density functions of MP, and HTMP, 7. defined in Theorems E. and E.2, respectively.
> (A-1) if
Now, we further identify the tail asymptotics of the limiting spectral density p(z) for two different cases of k(NN.
First, as x(N) + sc. applying (E.3), we know that the limiting spectral density of 2v0
MN
pxx) = 1 - xx- 7-)
E b ,
2x
assuming that y # 1. In this case, the density has bounded support. If = 1, then
1
22 x
p(x)
E (2,
2T
and p(x) -/200-0
Scaling these densities does not influence the tail behavior.
Next, we turn to the case k(N) + K E (0, oo) as N + oc, and we rely on two asymptotic properties of the Tricomi confluent hypergeometric function defined in (E.4):
(I) For any < 0, U(a, b)
T(1-6)
(=)) as |=| (DLMF, 13.2.22); and
M(a-b41)
(II) For any a, b, (a, b, 2)] |2|- |a| with arg 2|
ST (DLMF, 13.2.6).
From (I) and (II), the probability density function (x) given in (E.5) satisfies
(x) ~ f -1-4 as x. + 0+, and (x))(1) ~ Czx
t $ as x -> oc.
27

---

Models of Heavy-Tailed Mechanistic Universality
Thus, with change of variables, the density function p of & M satisfies 27B
p(x) ~ cτ 25 -1+1 as x oo, and p(x) ~ c2
by
-1-3
exp
K
as x - 0+.
2y
We now only need to observe that rescaling MN will not change the exponent in the power law, but it will change the coefficient in the exponential term.
E.4. Properties of
Marchenko-Pastur Distribution
Recall that the Stieltjes transform of spectral density p is given by Sp(2) [( =) p(z)dd. Thus, we have the following result.
Proposition E.3. The Stieltjes transform of HTMP,, 7. defined in (E.5) is given by
S Sp (2) = 1 ())((z)
(=)= I dx
U(n/2 + 1,2 + 2:
27
27
E CECRR. I.
(E.7)
24
Unn/2,
+1+ 2:
2y
The proof of the above proposition is given by Ismail et al. (1988).
From DLMF, 13.7.3, the asymptotic expansion for large of U (a, b, 2) is given by
c = (-1) [(a + b + ++
U(a.b, 2)).b,2)~2
(a % z-
(a,5)
(E.8)
(a+a)((a-b+1)
k=0
Using (E.7) and this approximation of (a, b, =). we can compute the first few moments of the HTMP V.A distribution. These are given the following proposition.
Proposition E.4. Let m, =So ²² Pr (x)dr denote the k-th moment of the HTMP, 7.x distribution for any E N. Then and b 1 3 2.4 In general,
mg 1, m. (a- b + 1), and m2 (a (u-b+ bb 1)(2a b + 2) where
k- 1
k (a.b)
Σ
(a.b)
my
k
mk 14
1=1
where of (a.b) defined in (E.8).
Proof. We seek an expansion of the form Sp (2) = Et m = as the coefficients my are precisely the moments of HTMP 7. Let z/z and b 1 + - 27 From (E.7) and (E.8), formally,
(a+1.++1) 2-h
A=U
mgz
k=
ck
(a,b)
2-
and so
Σ 1=
2
-k
(a.b) -
=
k
(a+ 1,b+1)
k=
k 0
Using the Cauchy product,
Σ
Σ
mk ma-10
(a:b)
x
(0+1,b+1)
k=0
1=0
k=
and so for t 2 O,
(a.b) (a+1,b+1)
mg
my 10
Cx
Evaluating my for ls E {0, 1}, we find that mo
(as expected),
mi = (a (a b + 1),
28

---

Models of Heavy-Tailed Mechanistic Universality
where we the fact that [( 1) (z).. Furthermore, for k 2, have
m2 (a b+ 1) b+ 2).
F. Structured Matrix Models: Weyl Formulae and Joint Eigenvalue Densities
In this section, we derive a change-of-varrable formulae (also called Weyl formulae) for obtaining joint eigenvalue densities with structured matrix models. The models we consider are outlined in Table F.I. The results proved in this section will be applied later Appendix G solve the variational problem (7) for various matrix models outlined in Table F.I.
Type
Description
Weyl Formula
General Form
(a) Diagonal
All off-diagonal entries are zero
Trivial
A diag( (A....., And
(b) Commuting block diag.
Block-diagonal with symmetric blocks A A. plus [A,, Aj] = 0 for all i,j
Thm. F.2
A = diag( (A1,.... As). where [Ai, A₁] =0
(c) Symmetric block diag.
Block-diagonal with symmetric blocks
Cor. F.4
= diag( Al...., Ak)
A = A!
(d) Kronecker-like matrix
Full symmetric matrix with commuting subblocks
Thm. F.5
A
[Aulli where
Asj = Al. and [Ay, Aui] =
0
(e) Symmetric
General symmetric matrix
Cor. F.3
A A
Table F.I. A collection of structured matrix models with explicit joint eigenvalue densities.
Computing Jacobians for matrix-valued functions presents a notationally complex task. To address this challenge, we employ the following technique derived from the change of coordinates relationships prevalent in Riemannian geometry, as outlined in (Menon & Trogdon, 2015).
Change of Coordinates Trick: Assuming that M = F(r,,..,, xNn is a matrix-valued function, let 9M = (0M);s where (0M); Ek= Jajadæ and Jijk 8F., the elements of the corresponding Jacobian matrix. Compute Ozc
tr(OM (MMMaaM)
gudradry
where
Jki
Jegk
Then for any integrable function f.
J [[M]JM= / [(F(E)....., xN)) / det g(I1,..., IN) dæ ... dzx
As simple example of this trick, Lemma F.1 considers the task of finding the Haar measure (the uniform distribution) on the space of symmetric matrices.
Lemma F.1 (SYMMETRIC MATRICES), The Haar measure on the group of N N real-valued symmetric matrices given by
dM 2
NIN-1)
dM
II
aMGx.
1=1
L5icKEN
29

---

Models of Heavy-Tailed Mechanistic Universality
Proof. Since dMy d.Mkk it follows that
tr(OMTaM) = Σ
dMjadMaa- = Σ dM +2 jck
dM3 jk
2.k
Therefore, for M E RNXN, 9ú = and 9ij = 2 for i + j. Hence, det = 1N NIN-1) result follows.
M(N-1)
2 2
and the
We are now interested in applying the trick to obtain joint densities of eigenvalues starting from densities over Haar measure on subclasses of symmetric matrices. The most famous result in this regard is the Weyl formula for the joint density of eigenvalues over the general class of symmetrie matrices. This is presented in Corollary F.3, as a special case of a useful generalization that present next in Theorem F.2. This next result treats our commuting block diagonal matrices (type (b)), and is inspired by McCarthy (2023), involving d real-valued symmetric matrices that are simultaneously diigoooiiiyydigggaliaaabl That is, these are families of matrices MI...., MM, such that M, = QA,Q", where Q is orthogonal and each A, is diagonal. Note that N = mn, then N x N matrices which are block-diagonal with commuting m x m blocks can be treated in this way.
Theorem F2 (WEYL FORMULA: TYPE (B)), Let dM be the Haar measure on the space Mm of n-ssmuultaneously diagonalizable matrices of size m x m. Let f M + R be such that f satisfies f )(QMQT)= = fMM where Q E Rm is orthogonal and M& E M Then, there exists a constant Cm > O independent of f such that
J [(M)dM = Cm,, / S(N)
II
12, - A
1.j=1.....
where A (A1,.... An) and A, diag (11:: ...y dim
Proof. Immediately, restrict the support those matrices with distinct eigenvalues. One readily verify that the set of matrices with non-distinet eigenvalues has zero measure with respect to the Haar measure symmetric matrices, so this does not affect the density. Let A =
so that M =
where M = (M,...,, M,). Then
0M dQ((A A Q(dA ................................ d QAI............................. A, AdQ
= Q [ddA)....,dd dA,) + Q T (A,,....,AAA) + (A QT
Since QTQ== = I, this implies that dQ Q = QQ aQ, and so
QM Q [(dA 1 ..... dA,)) + Q dQ(A) A.) -AA (A......AA)) dQ] QT Q [aa + Q dQ, A]] TT
Consider the variable A such that dA - Q dQ, so that 8M Q(dA + AJA, Now.
n
tr(MM 0M) Σ (r(dA))++tr(aal)))))) + 2tr A + tridA, A,12
Observe that for a diagonal matrix P =
.Pm }. we have that
[dA, Plaj = (dAP- PdlA); = dAgp, - -padAij = (pz - A.)d
and so
m Σ
tr[[A,, P13
dA, PlaydA,
1,==1
m
Σ
(p) - p:)dAis(pp p3 dAA
1-=1
Σ (p) - dd)
-1=1
2 Σ (p3 p.) Pd
Kj
30

---

Models of Heave-Tailed Mechanistic Universality
where the last two lines follow from the fact that dA dQ anti-symmettiic, and (dA),: and (dA)a -ddA) Similarly,
tr(dP(dA,P))
Σ
dpı[dA, Plkx 0,
since [dA, Plkk 0. Consequently,
n1 m
tr(aM 0M) Σ>
d 2
|A
AxIPdA33
(F.1)
jck
Therefore, using the change of coordinates trick,
|1MMMM=2
= milo-1
I [(A, A) II 1A,
A.
III
II JCk
dA gkc
JCk
and for f independent of A, the integral over the measure A separates, and so for some Cm,n > independent of f, we have that
I S(MMMM = =Cm, / S(A) II
1x, - A.Al
JCk
From this result, we immediately obtain the following two results as corollaries.
Corollary F3 (WEYL. FORMULA: TYPE (E)). Let dM be the Haar measure on the space of real-valued symmetric matrices of size x N. Let f a matrix function that depends only on the eigenvalues. Then there exists a constant CN independent of f such that
N
I f(M))M -N f() AN)
II
1
A II ==1
d\,
J.k=1 I........N
jek
Corollary FA4 (WEYL FORMULA: TYPE (C)), Let dM be the Haar measure on the space of real-valued symmetric matrices of size N x N with N = mn, that are block-diagonal with block size x m. Let f be matrix function that depends only on the eigenvalues. Then there exists a constant Cm,,. > independent of f such that
/ {(M)dM = Cm,, [ II II 1=1 j.k=1.....
MNg - AdA
jCh
Finally, we consider the case of general block matrices with commuting blocks. This case is slightly different than the setting in Theorem F.2, as orthogonal decomposition adopts Kronecker-like structure. However, we can use the same techniques to obtain the joint eigenvalue density.
Theorem F.5 (WEYL. FORMULA: TYPE (D)). Let dM be the Haar measure on the space Mm,a of mn x mn real- valued symmetric matrices comprised of x n commuting m x m blocks. Let f Mm. + B be such that f satisfies I(QMQT) J(M), where Q E mmmxm orthogonal and M E Mm Then there exists a constant Cmn, > independent of f such that
1/2
1/2
I F(M)dM -Om,n / F(A))
"
m (De )us)2
II
Σ k11=1
(A
II
dA.
Li=l...,m
i==1....m
k,1=1
i<j
where A (A1,..., An) and A; = diag (iil ...,, dim)
31

---

Models of Heavy-Tailed Mechanistic Universality
Proof. The proof proceeds in similar fashion to the proof of Theorem F.2. Once again, the support of the density is appropriately restricted, without loss of generality. Now consider the change of coordinates
M = (Q + O) (AI................. An (Q 8 0),
where Q E IR and O E Rmxm are both orthogonal, and so U Q 2 0 is also orthogonal. Note that dU
(dQ @ 0) + (Q @ dO). As before,
0M U(AA + [U (U,U AUJA))
= =UddA + [[Q &
@ (0) -UAdA [(QTdQ [) + [(T dO),
d0. Since dA and dB are both antisymmetric, for Now, consider variables A and B such that dA dd and dB dO.
dA dB (dA I) (I dB), have that
(dA @ dB) (dA² @ 1) + (T@d (dA 1) (1 dB) = - -dA @ dB),
and it is therefore also antisymmetric. Recall that
(dA & +m((-))+tm---))++ =dA; dAi (I
dBBi:
and
[dA 69 I), Am) (Aşt Aex ) [(7 & dB), Am( = (Azz dua )d
Consequently,
tr([(AA @ 1),A1)) =2 Σ
Σ (A Ax))
dAg
kd
tr([(T @ dB), = Σ
Σ 2
dB
b
On the other hand, since dAA and dB are antisymmetric,
( & Al)
(Ax dd)) 0.
i.j.k.l
Altogether,
tr(OM T 9M) = Σ>
Srه +2
Σ Kgj
dAi
+
E - RCl 1.
Aw)2
dB
(F.2
=1j=1
R.
The result follows by the change of coordinates trick.
Remark F.6. each case, the Weyl change-of- variables factor the joint eigenvalue densities involves absolute differences of eigenvalues |Λ; - >yI. and it varies only how these differences are included. These factors determine the degree of eigenvalue repulsion, they assert that configurations with eigenvalues close together should occur with reduced probability. This plays central role in the shupe of the spectral density. An important observation for us is that the strength of these repulsions only increases as we take products of cigenvalue differences, since sums do not further reduce the probability density. Heuristically, we can count the number of multiplied eigenvalue differences as a rough approximation of the degree of eigenvalue repulsion in the spectral density. Let us call this quantity 6. From the Weyl formulae, for each structured matrix class, is as follows.
32

---

Models of Heavy-Tailed Mechanistic Universality
Diagonal: 0 0
(b) Commuting block diagonal: m(m-1)
2
(c) Symmetric block diagonal: N(m-1)
2
(d) Kronecker- like matrix: 0
m(m-1] 1)
n n(n
(c) Symmetric: A N(N-1)
2
Clearly, 0 becomes larger the matrix class becomes less structured (more degrees of freedom).
In principle, one measure the amount of matrix structure according to 0. Finding the spectral density of structured matrices from their joint eigenvalue densities can be very challenging, but this becomes much simpler under an approximation that enforces all eigenvalues to be interchangeable that is, when the joint density is symmetric. Letting .e.A......A Ax) denote the Weyl change-of-variables factor, we can approximate w using Schur functions (Noumi, 2023, $3). This process was recently proposed in the context of RMT in Kimura & Mazene (2021). To reduce complexity, we can further observe that the behavior of will be dominated by the leading-order term Ili,, [ Xyl" for some y (Noumi, 2023, pp. 21). Therefore, class of approximations that should respect the degree of eigenvalue repulsion, measured by A weights each eigenvalue difference according the ratio / ((N-1) so that
II IJEE
Ayl
becomes
II 65=1.., (-)
X Ay N(--
This is only a heuristic, however, and this principle must be made more precise. In Appendix G, we show that variational approximations using beta-ensembles replicate this behavior a principled way.
For completeness, we consider one last example of a Weyl formula which will help to make the comments in Remark F.6 more explicit. Let M be a symmetric matrix with eigendecomposition QAQ". Here, we consider the scenario where only columns of Q are variable, and the other N are fixed orthogonal vectors. We refer to this scenario as the "d free eigenvectors" case, only d of the eigenvectors of M are allowed vary in the random matrix ensemble.
Corollary F.7 (WEYL FORMULA: FREE EIGENVECTORS). Let dM be the Haar measure on the space of real-valued symmetric matrices of size N x N with d free eigenvectors. Lel J be a matrix function that depends only on the eigenvalues. Then there exists a constant CN,d > independent of f such that
N
/ (anjaM-Cx.. f(Ai...., AN)
II
An II
a>
1k==1....dd
1m1
jck
Proof. Let Q = VV where U E Nxd is variable and V R. N is fixed. Then
U QTQ=
vTV = I,
VTU vTv)
and
d(Q Q) = (dUT UU wTy)
( dT
This implies that V dU = dU 0. Consequently. since dQ = (dUlo),
dA =A=QTdo (UTdV 0) (UTaU) au
9)
dT 0
This implies that dA 0 whenever i d or > d. From (F.1),
ur@MTaM Σ dx +2
Σ
1gg Xel AdA,
1.k=I... jck
and by using the change of coordinates trick, the result follows.
33

---

Models of Heavy-Tailed Mechanistic Universality
Returning Remark F.6, we can count the number of eigenvalue differences in Corollary F.7 and find 0 d2 This provides direct relationship between the number of "degrees of freedom" in the matrix structure according to the number of free eigenvectors, and the resulting Weyl formulae. We will return to this example in Appendix G.
G. Variational Approximations with Beta-Ensembles
In this section, we use the Weyl formulae from Appendix to solve the variational problem (7) and to derive the relationships in Table for the matrix models outlined Table F.T. Our derivations begin in Section G. with cases (a), (c), and (e) using an exact formula when the coannge-f-vvriablle factor of particular form. To cover cases (b) and (d) (using (c) only to test the validity of our approach), estimate the relationships using symbolic regression applied to estimated solutions of (7). This approach is centered around a numerical method outlined and justified in Section G.2. The experiments are performed and conclusions drawn in Section G.3.
Assume that a matrix M E RNNN has elements distributed according to a density p(M) x II e e V(A;)) over some class of matrices M, where ,I,...,Ay Ay are the eigenvalues of M. Such densities typically arise the form p( M) x (det ( ( for some analytic functions f.g. After performing a change of variables, the corresponding joint density of eigenvalues adopts the form
N
p( AN) & wh1 AN)
v(A.)
(G.1)
where w depends on the underlying class of matrices M. To provide a concrete bascline for the study of models of this form, we consider approximation the joint density of eigenvalues by the Laguerre beta-ensemble
qa(A AN) 2(3)
II 1K3
(G.2)
Note that | = 1 corresponds to the case where M is the set of real symmetric matrices, and | = 2 corresponds to the case where M the set of complex Hermitian matrices. Intuitively, smaller values of 3 correspond to more restrictive class M. Tridiagonal matrix models exhibitinn joint eigenvalue densities given by (G.2) were discovered Dumitriu & Edelman (2002). To approximate (G.1) by the variational family (G.2), we consider the forward variational approximation:
8 argmin dKL (gsllg)))))
R² N
(G.3)
B20
G.1. Counting Eigenvalue Differences
We shall now confirm that the variational approximation (G.3) replicates the desired heuristic behavior outlined in Remark F.6. In Proposition G.I, we show that when the change of variable factor wo is of product form, then 3 is explicitly determined by counting the number of eigenvalue differences.
Proposition G.1. For probability density p satisfying (G.1) with w(A) IIJE - 1 letting f = B.j.
20
argmin dRl (gs (p) N(N - 1)
820
Proof. A primary tool in the proof is the explicit expressions for the derivatives of F(9) = log Z(B), where
Z(B)
exp
Σ v(A)
log |4, -111 dd. Σ Kej
from it follows immediately that
(3)
F()))
Egg Σ
log 1>, - -2g1
7(3)
i<j
(3) (3)
(3(
2
Var Σ
log |A, -2,1 1/>00 0.
Z(3)
Z(3)
34

---

Models of Heavy-Tailed Mechanistic Universality
Note that since gs symmetric in its arguments, each eigenvalue is exchangeable random variable, and
20 N(N - 1) F3)
Egga
Beg log A - -111=
Σ
Bg log A Xy1 = Egg log (11
X2
Σ
Oij
AEE
ijEE
BE
On the other hand,
gg - -
Altogether, letting Zp denote the normalizing constant for p.
dKL(gsllp) Egg log
[qa(......... ,N)))))))) (
p()1 ((A).....NN)
log Zp (
20
-1) F((3) F(B).
N(N
Taking the derivative, obtain
d dk (gg (p)) =
20
(3),
d3
N(N - 1)
and since (")3) > 0, the only critical point of B + dKL(gg)lp) is N(N-1)
20
To show this is minimizer, observe that
d²
20
dkk kKL
(q=llp) (3) N(N 1)
F*)(B),
and at the critical point, d
d² dkL (p) (3)
From this, we obtain the following result as a corollary.
Corollary G.2. For the following structured matrix models, the values of 8", K² defined (G.3) are given by:
(a) Diagonal: 8r n² O.
(c) Symmetric block diagonal: 8" m=1 K² (m 1) N-I-
(e) Symmetric: 3 1, K* N.
Proof. The diagonal and symmetric cases are trivial. For the setting of block-diagonal matrices with
w (1 Aw)
Daj
Aael
1 j, 1.k=1..m
jek
as shown in Corollary F.4, 0 == n m(m-1] 1) and so Proposition G.1 implies f am( 1)
OO
mmm-
Performing similar analysis for free eigenvectors scenario in Appendix F, obtain Theorem G.3. Theorem G.3. Consider the density (2) with N satisfying an inverse- Wishart distribution over the space of symmetric matrices with d free eigenvectors. The variational approximation (7) satisfies x² = d d-1 d2
N-1 N as d, N + oc.
Proof of Theorem G.3. Applying Proposition G.I to Corollary F.4, we find that B" = 20 (N(N 1)), with By = (i.j d} and
d(d -
d(d-1)
Σ
Ojj
and so 8 N(N -1)
2
Then the result follows.
35

---

Models of Heavy-Tailed Mechanistic Universality
G.2. General Numerical Method
At higher generality, it difficult make concrete statements about the behavior of f". Instead, we can appeal to numerical methods. Solving the optimization problem (G.3) directly is feasible using iterative solvers, but is expensive and can be unstable. Instead, Lemma G.4 shows that 3 is the solution to convenient fixed point problem, and hence it can be estimated by form of stochastic fixed point iteration.
Lemma G.4. Any local minimum 3 of dKL( satisfies 3 = R(3"), where
N R(B) = Cov (log ('i................... AN) Eij=,, ==,, log Ar Xyl)
(G.4)
Vargs Eij=liki log AA Az1
Proof. Letting Zp denote the normalizing constant for
dkI.(gsllp) = log 2pp F(3) Egg log w + (B).
Any local minimum dRL (fallp) also critical point; taking derivatives 3 yields
d
d dRL (gs (gslp)) (3) Egg logg w dg
log qs
d3
However, since
d
log qs E log Dhi - FF
d3
which is zero-mean under qa, it follows that any local minimum of dKI(ggell)) satisfies
(F) = Cov
log Elog i<
log Ai Xz1
The result follows.
Suppose that for each & = ,,...,, and fixed 3 2 0, (21 (3...))
(3a)) sid qa. Letting
23,k = Σ
log (B.)) 1
i==ı,iij
U3.k = log w/ (A)
1 (B.K)) (B.k)
we can see that consistent estimator of (G.4) is given by
E (Ig F.3 (ys, 033) ( - zz)²
(G.5)
where ã3 = Et= kg and W/3 = = UR. Since (G.5) exactly the estimator for the slope simple linear regression,
it is also unbiased. Starting from Bg N. for fixed 0 < 7 < 1, consider the iterations
8 = (1- (1--21) B, ? Σ (I B I Ez ) V8.)
(G.6)
++1 ΣL (x) Br.k )2
The following central limit theorem for 3, follows from Zhang (2016, Theorem 1.1) and the unbiasedness of (G.5). Proposition G.5. Let 8 denote a critical point of in distribution to a zero-mean normal random variable as r + cc.
Then for sufficiently small 7 > O, vr(ß, - 3r) converges
Proposition G.5 suggests that the iterations (G.6) provide a reliable means of estimating 3*. practice, we have found that 7 typically suffices. The complete pseudocode is provided Algorithm 1.
36

---

Models of Heavy-Tailed Mechanistic Universality
Algorithm Stochastic Fixed Point lteration for Estimating 8
input matrix size N N, weight function w, step size 7 (e.g. 7 = 1), number of samples p
Initialize 8g (e.g. 3p N
5.
2: for 0, 2,
until convergence do
3: for k = 1 to p do
4: Sample (AI...., AN) 93,
5: Compute:
6: Tk = iij bijj=liik iis log 1r )3>
7: 1/ = log w (11,
8: end for
9: Compute means:
10:
Ek=1 The y= E Us
11:
Compute covariance and variance:
12:
Say Σ 11 (IA (Vk 9)
13:
Sax kmn (TM- Z)²
14:
Update B:
Br+1 = (1- - 3r + न 8ex
15:
Sae
if convergence criteria met (e.g.. 13,41 3rl tolerance) then
17:
Return 3
end if
end for
G.3. Variational Approximations for Commuting Matrix Models
To complete Table 2, we now look to use Algorithm with symbolic regression (Cranmer, 2023) to estimate the behavior of with respect m, n, and the shape parameter of the Laguerre beta-ensemble. To sample from qa. rely the tridiagonal matrix model of Dumitriu & Edelman (2002) (see Algorithm 2). For eigenvalue computations, recommend an implementation of the DSTEBZ routine in LAPACK, such as the eigyalsh.tridi routine in SciPy.
Algorithm 2 Generating a Sample from the Laguerre 8 EEsemmll
input 8 0, matrix size N E N, shape parameter
I: = a a + + 3(N - 1)
2: for -= 1 to do
3: Sample d, from Xa-9(i-1)
4: end for
5: for to N do
6: Sample from Xa (N-
7: end for
8: Build diagonal entries
D, = + (i=1...., n- -1), Dn = dc.
9: Build offdiagonal entries
E, = d I.-,....,7 1).
10: Compute symmetric tridiagonal matrix eigenvalues with diagonal D, off-diagonal E
To collect estimates of nM for variety of scenarios, consider ten values of both m and logarithnically-span -pacce between 3 and 100 and 1.2....,5 5. In Algorithm I, we let p 50, 7 1 and choose our tolerance in each case to be 10 10-3/N. We then perform symbolic regression these estimates with respect to m, n, and a using PySR (Cranmer, 2023) with maximum size of terms, population size of 20, 5 iterations, addition, multiplication, and division binary operations, and logarithmic and exponential unary operations.
To verify that this process capable of reconstructing expected symbolic relationships for s*, we consider the symmetric
37

---

Models of Heavy-Tailed Mechanistic Universality
block diagonal in Figure G.I. The predicted relationship x= (m 1.1) 1.05, close the exact relationship
N n" (m ---1 1) N-I given in Corollary G.2.
symmetric block diag
SR formula: (m 1.0958514) 1.0424111
120
data
100
y-x
estimated X*
40
20
40
60
Bd
100
120
symboliic regression prediction
Figure G.I. Symbolic regression predictions vs. estimations using Algorithm for x² under the symmetric block diagonal structure over the Laguerre beta-ensemble.
Next, we turn our attention to the commutative cases where explicit expressions for x are not known. The results are presented in Figure G.2, which suggest the relationships:
(b) Commuting block diagonal. R² 1 (m 2
(d) Kronecker-type. K²
These relationships closely follow the heuristic in Remark F.6.
commuting block diag
Kronecker- type
SR formula: (m 0.49962232)
SR formula: (n / m) + (m /n)
y=x
y=x
30
30
estimated k²
estimated
20
20
0
symbollc-regression on prediction
symbolic-regression prediction
Figure G.2. Symbolic regression predictions estimations using Algorithm for x" under the commuting block diagonal, and Kronecker-type matrix structures over the beta-Laguerre ensemble.
H. Comparison of Neural Scaling Law for Linear Regression
Modern scaling laws in ML attempt to quantify how test errors depend on key resources such as model size and dataset size. While many such laws have been observed empirically deep learning, there are also rigorous proofs in simpler, more tractable settings (c.g., linear models or kernel methods) that exhibit qualitatively similar phenomena. However, most previous work focuses linear regression the training dataset (Lin et al., 2024; Bordelon et al., 2024a).
In this section, we consider concrete power law setting for the data covariance, and we review some previous results of an associated ridge regression scaling law, making comparison with Proposition 5. which we proved in Appendix C.3. This high-dimensional ridge regression has been studied in many different scenarios, e.g.. Hastie et al. (2022); Wei al. (2022); Li et al. (2023); Defilippis et al. (2024). The main difference between our Proposition 5.1 and previous results is that we considered the ridge regression on a heavy-tailed feature matrix.
38

---

Models of Heavy-Tailed Mechanistic Universality
Specifically, for linear regression, consider the dataset be composed of, for 1,...,
The feature vectors I; E Rdd have covariance matrix E, whose eigenvalues {u, { exhibit a power-law decay
C ja
My (2)
for some constants > 0 and a > 0, for large j.
(H.1)
The label is generated via linear model y; (Tz, wn) Exi with i..d. Gaussian noise Et (0,,
We show how the test mean-squared error (MSE) of the ridge estimator scales with n, d, and an appropriately chosen ridge parameter An. Under mild conditions (detailed below), the MSE follows particular power law in n whose exponents depend and on how scale A, with This illustrates how data spectral structure drive nontrivial scaling phenomena in high-dimensional learning.
Given data matrix X E R of n samples (each row is n.J) and labels y = (y1....,, Vn)", ..) the ridge regression estimator with regularization > is defined by
arg min wERd
{ [ + {||w||*} (x x Ala) 1'
We will let À = An possibly depend on n (and/or d), as is common practice for controlling varianer tradeoff's in high dimensions (Wei et al., 2022; Li et al., 2023). The generalization error is given by
(A,)( = Ex (z. T6) - (x, wn))
where ~ NN(0 is independent of the training set.
We now present the asymptotic behavior of En,((An) under the spectral assumption (H.1) and a suitable scaling of An with n, as n, + oc. For simplicity, let UAU be the spectral decomposition of Σ Ra i... U is orthonormal and A [(ig((1,,... -..,........................ .a)) contains the eigenvalues (ordered). We analyze ridge regression in this eigenbasis. If we let ã U x4, then ā; (0, A), i.e., the coordinates (ī:), independent with variance μlj. Similarly, w" can be written in the same basis, = = Uw" so w(() (w"); in that basis. Thus, when 2; (O0 E), the matrix X = [aT :..;;.,;;,
has an SVD tightly connected with 11/2, but we simply recall the well-known result that in the basis U, ridge regression amounts to a coordinate-wis wise shrinkage: where
M4 M4 Am W3 + (variance term from noise).
U
Hence, in expectation (conditioning on X or in integrated sense), each coordinate is shrunk by factor 14
Bias-variance decomposition.
For [N, Σ), the test MSE can be written as
(n,d(An)
Bii(@ w 1
My (bias} vary),
MSF the T.-iner producct
coondinate wise contributions
where
var, = E - 1(
biasy
[
Ri Wjj. and the variance part, vary, scales like 22
Ignoring constants, one finds that the bias part, biasy reflecting how the noise passes through the resolvent
tg+ (XT X + 1n
(nytA.))
As n, d + oc, we replace the sum by an integral and use that b ~ C/ja. Then, the test MSE is roughly:
Eyg ( My Mg + An - 1) (wz)²²
d²
En,a(Dn) ~ Σus
Pj (Mf A ) n
bias term
variance term
If || < xo in that basis, we can assume (1072 = Oj-0'') for some B' > 0 or at least that w= is well-defined in E's eigenbasis. For simplicity, one often takes (ız²)2 = O(1) or O(j we will highlight only the main effect of the covariance spectrum {1z} here.
39

---

Models of Heavy-Tailed Mechanistic Universality
Choice of A and the Emergent Power-Law Rate. When by Cj-a, the sums above partition into low-index terms
(where j is small and Py is large) and high-index terms (where j is large and μ4y is small). Roughly speaking:
If Xn is too large, then even the large cigenvalues μ4j get shrunk heavily, leading to big bias.
If Xn is too small, then the small eigenvalues (at large j) produce a large variance.
Hence the optimal ridge parameter A, in the presence of power-law {uz} often takes the form Xn x n 2a+I, giving a trade-off that yields final MSE scaling like
(np() x n 2a+1.
(H.2)
For more details, we refer to Li et al. (2023).
Comparison with Proposition 5.1. The scaling limit get in Proposition 5.1 is the product of the power of heavy tail and the ridge scaling, whereas in (H.2), the decay rate will be dominant, by the power law of the data covariance, when we choose the ridge parameter An based on a defined in (H.1). In this sense, our scaling limit is more general, and it really depends on the combination of the regression model structure X, and the heavy tail distribution of the feature matrix. On the other hand, the above derivation of (H.2) for linear regression considered the anisotropic dataset and label noise, which not handled in our Proposition 5.1. We leave this for future work, where we will be able combine the heavy tail distribution in the feature matrix, the power law for the dataset covariance, and the model structure (for instance A, in ridge regression). This extension may require deeper RMT for inverse Wishart matrices, such as local law results by Wei et al. (2022), deterministic equivalence for (C.7), as in Defilippis et al. (2024).
40


Reviews:
[{"review_id": "oaI27RU7G4", "reviewer": "Reviewer_RfHj", "summary": "This paper explores heavy-tailed mechanistic universality (HT-MU) in deep learning by proposing a new family of random feature matrix models based on the \"high-temperature inverse-Wishart\" ensemble. The paper reviews two mechanisms and presents a third one for the emergence of power laws in different matrices related to trained neural networks. These three mechanisms are (i) complex correlations in data, (ii) reduced training temperatures, and (iii) implicit architectural biases that affect eigenvector entropy. The authors provide theoretical results linking their model to neural scaling laws, the five phases of learning, and optimizer trajectory properties.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful examination of our work and for the thoughtful feedback. We appreciate that they find our framework to be ambitious, and understand the concern about empirical verification of individual claims within our theory. One of our fundamental assumptions---that the eigenvectors are not uniformly distributed---was verified by a reference provided by Reviewer Weri; (Thamm et al., 2022), and we will include the reference to this. However, our primary goal in this paper is the development of a single theoretical model designed explicitly to display observed asymptotic spectral tail behavior (Theorem 4.3). Our discussion is used to motivate the construction of a variational family with the desired characteristics. We study this model to compare its overall predicted tail behavior against empirical phenomena. We cannot claim that this model represents specific architectures or matches empirical phenomena identically; our analysis cannot replace other fantastic work that is being conducted in this direction. To help focus our claims, we add the following to line 53:\n\nOur central objective is to identify a plausible random matrix class that exhibits inverse Gamma law spectral behavior, and the smooth transition from Marchenko-Pastur to heavier-tailed densities observed in Martin & Mahoney (2021) and Thamm et al. (2022). \n\nand replacing 63--65:\n\n- construct a parametric family of spectral densities which includes the Marchenko-Pastur law, but allows for heavy-tailed spectral tail behavior in line with empirical observations (Figures 1 \\& 5, Appendix I)\n\nWe replace lines 125--135 with the following: \n\nHowever, PIPO also fails as a standalone learning theory, as it does not account for implicit model biases. If data _alone_ influences tail behavior, models trained on the same dataset should exhibit similar power laws (Section 4.1), which contradicts empirical findings (Yang et al., 2023). At present, theoreticians analytically examine the interactions of individual models in the presence of heavy-tailed data (Maloney et al., 2022), but such analyses are intractable at scale. Originally conjectured by Martin & Mahoney (2021), here, we look for a third alternative: a universal, model-agnostic mechanism that can give rise to different heavy-tailed spectral behaviors from the same dataset. We construct a sequence of hypotheses that lead to such a mechanism, which we refer to as \"eigenvector entropy\". Matrix models in the literature typically have eigenvectors that are Haar-uniformly distributed (maximum entropy) (Anderson et al., 2010); this is also represented by _delocalization_ (Bloemendal et al. 2014). Breaking this property, reducing entropy, our framework provides a family of variational approximations (HTMP) with the right qualitative behavior: arbitrary power law spectra, and inverse Gamma laws that can arise from model design alone.\n \nRegarding more specific comments:\n- $M$ should be positive-definite, not $\\Phi$ (typo)\n- \"Near ballistic\" has been changed to \"power law exponents $\\alpha < 3\"\n- Donsker-Varadhan now cited\n- Hanin & Nica discuss heavy-tailed log-normals; not power laws, but can appear indistinguishable, see (Clauset et al., 2007). We now make this clear.\n- \"Free eigenvectors\" replaced with \"eigenvectors $v_1,...,v_N$ with $v_1,...,v_d$ Haar-uniform and $v_{d+1},...,v_N$ fixed.\n- The approximation (6) is proposed in accordance with patterned matrix models whose spectral densities we calculate in Appendix D. It is not exact, but designed to replicate qualitative phenomena.\n- See response to Reviewer Weri regarding 5+1 phases.\n- For further citations to cover claims made about the wider literature, see the response to Reviewer Weri.\n- Section 2.5 of (Anderson et al. 2010. An Introduction to Random Matrices.) showed that eigenvectors of GOE/GUE are Haar uniformly distributed; we now include this.\n- For general Wigner or sample covariance matrices, (Bloemendal et al. 2014. Isotropic local laws for sample covariance and generalized Wigner matrices) used local law results to prove the delocalization for all eigenvectors, which means all the entries of the normalized eigenvectors are $1/\\sqrt{N}$ order. For more details, we refer to local law lecture note by (Benaych-Georges \\& Knowles 2016).\n\nTo answer your questions:\n1. Yes, but this assumption is not needed if $M = \\Phi \\Phi^\\top + \\frac{\\gamma}{2\\tau} I$, so we take this in the updated version.\n2. Thank you for the suggestion; this now reads \"For more general losses...\"\n3. Linear regression with equal numbers of parameters and datapoints is one example, since $\\min_\\Theta L = 0$.\n4. The Appendix includes $L^2$ regularization, but otherwise the two are equivalent. We now add $L^2$ regularization to the activation matrices section.\n\nPlease let us know if you have remaining questions about claims; space constraints prevent us from addressing further comments, but we are happy to provide evidence in the discussion period.", "claims_and_evidence": "While the theoretical framework is ambitious, the logical flow of the paper is extremely hard to follow. \nThe paper seems to review and reject the PIPO (Population Covariance) approach; they reject it based on \"More recent analyses have shown how architectural decisions can alter the power law, but these hold only for specialized models.\" (lines 245), however, they give no reference to such recent analysis.\nRecursive model structure is reviewed and rejected in much the same manner.\nFinally, they come to their suggestion of reduced eigenvector entropy. In this section, the manuscript makes a number of confusing statements.  \n- What do the authors mean by assuming $\\Phi$ describes a positive-definite matrix? $\\Phi$ was introduced as $n \\times m$ feature mapping, with $n$ being the number of datapoints and $m$ the number of features. Generally this is not a square matrix, nor is it positive definite. When is such an assumption expected to hold?\n- The manuscript moves on to stating \"In this change of variables, it is typically assumed that the distribution of eigenvectors (Q) is uniform.\" Uniform over what domain? (the d-sphere?) Why is this assumption typical? It would be helpful to cite relevant works. \n- The concept of eigenvector entropy is never introduced and a reference for it is never given. In fact, I was hard pressed to find references to it in the literature.\n- The paper goes on to talk about \"free eigenvectors\", these are never defined.\n- This is followed by: \"Admittedly, the link between (5) and non-uniform eigenvector distributions is nontrivial, so several points of evidence are in order\",  the status of much of the statements that follow is unclear, is it a conjecture, an interpretation, or a known result?\n\n\nThe paper follows with applications. They start with the 5+1 phases of learning, this phenomenon is not nearly as widely recognized as scaling laws, for the manuscript to be self-contained, it should include at least a brief explanation of what these are.\n\nThe authors need to overhaul the presentation to clearly delineate how their evidence supports each claim, ideally with concrete examples or references.", "methods_and_evaluation_criteria": "See above", "theoretical_claims": "The paper makes it difficult to assess its claims. It is often hard to distinguish between known results, conjectures, and rigorous results in the text.", "experimental_designs_or_analyses": "The authors show that finite temperature inverse Wishart distributions fit experimental data in several settings which they label according to a nomenclature which is not clarified in the text. Putting the latter point aside, since the theoretical basis for this distribution isn't clear, much further numerical evidence and precision tests are needed in our minds to make this claim sound.", "relation_to_broader_scientific_literature": "The manuscript cites works of Mahoney extensively, in a way that stands out compared to usual academic standards. This choice stands in sharp contrast to the relative thin reference to existing literature, e.g., citing the work of Martin & Mahoney along side a single other paper as examples of \"statistical mechanics of learning\", though it is not a review or a book, and there are much more prominent works in the field. The same repeats for \"robust metrics to assess model quality\", a work by Martin & Mahoney is cited together with a single additional paper, again, for such a rich field, a review would be more useful to the reader, or a choice of more canonical works.\nAt the same time, the paper simply lacks crucial references, for example, the paper mentions the Donsker-Varadhan variational formula (Appendix G.1) but does not cite the paper.\n\nAdditionally, the paper is full of jargon that is not common in the literature, for example, the paper cites Arous & Guionnet, 2008 for the statement \"While independent matrix elements exhibiting near-ballistic power laws can give rise to heavy-tailed spectral densities\" but the jargon \"ballistic\" does not appear in the cited work. \nTo continue this problematic line of misreference, they cite Hanin & Nica, 2020 to support the statement \"These results assert power laws in the elements of feature matrices.\" but Hanin & Nica do not mention power laws.  \n\nThese are just some examples I checked, I'm sure they are not the only ones.\n\nI recommend the authors completely revise their choice of references, as they currently give a partial and skewed view of the literature.", "essential_references_not_discussed": "The whole choice on references should be revised.", "other_strengths_and_weaknesses": "Strengths\n\nThe paper takes on an ambitious goal of establishing and explaining a universal behavior in deep neural networks. This goal is ambitious and inspiring. \n\nWeaknesses \n\nThe paper should be rewritten. At the current state, I gauge it would be inaccessible to the large majority of the ICML community, and the large majority of those whose research field is theory of DL.", "other_comments_or_suggestions": "In Table 1, the manuscript mentions \"observations\" but gives no reference to these observations.\n\nThe manuscript's notation is confusing and inconsistent.\n\n1. Eq. 3 includes $\\tau$. $\\tau$ is not presented before in the manuscript and it is only presented in Appendix G.1.\n2. $\\rho$ has at least three meanings, two appear in the \"metatheorem\", and then $\\rho$ is introduced again as the ratio $\\gamma / \\eta$.\n3. What is the definition of $N$ in section 4.3 and on? In section 5.2. it seems to be the dataset size, which was previously denoted by $n$.", "questions_for_authors": "Questions:\n1. Below Eq.3 is it assumed that $M=\\Phi \\Phi^T$ is invertible? Is that mentioned in the manuscript? Is that supposed to be a pseudo-inverse?\n2. Could the authors specify what they mean by \"For more general classes of models, computing marginal likelihood becomes intractable\" in line 195? What is being generalized? Is it the choice of loss function?\n3. Later in line 207 the manuscript focuses on the scenario \"However, if L(Θ∗, Φ) is constant in Φ , that is, any choice of Φ yields the same training loss\". What could be an example of such a scenario?\n4. Under the \"Activation Matrices\" subsection $y$ is introduced as a regression target, but in Appendix $G.1$ $y$ has a different definition, together with a new notation $\\epsilon$ which is never introduced. Could the authors clarify their setting? Do the two definitions of $y$ somehow coincide?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "oZ776XH4nZ", "reviewer": "Reviewer_Weri", "summary": "This paper argues that many phenomena observed in neural scaling laws arises from universal random matrix theory effects which the authors term heavy tailed universality. The paper introduces a theory which breaks up the deep network optimization into an optimization over features and optimization over the last layer weights. The relative strength of learning for these different components are controlled by by a hyperparameter $\\rho$. This leads to a *master model* for the feature matrices that the authors are able to compute the eigenvalue densities of such adapted kernels using techniques from random matrix theory. Depending on the hyperparameters of the model, there can be 6 types of observed feature spectra. The authors make connections to scaling law literature.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive assessment of our work, and for providing several references whose inclusion greatly enhances the working document! In particular:\n\n- Thamm et al. do an excellent job of highlighting the phenomenon we are interested in and provides further evidence for our case (their Figure 4 shows a proportion of eigenvectors are non-uniform)\n- We thank the reviewer for providing the neural scaling law paper by Bordelon et al. 2024, and we admit that this paper considered the power law assumption for the Fourier spectra of the kernels at initialization. This feature learning setup is similar to our power law assumption of the feature map and is the motivation of our neural scaling law section. They also consider a similar setup to the simultaneous training of weights and features in our Proposition 3.1. Very cool! We will modify the literature review in our neural scaling law section.\n\nThe other references seem to provide excellent examples of more fine-grained analyses with other precise models of training and architectures. We agree these provide excellent context and a better illustration of the state-of-the-art in the theory of feature learning. \n\nIn addition, we include the following references:\n- (Pillaud-Vivien et al., 2018) also assume power law features.\n- (Mezard & Montanari, 2009) for statistical mechanics of learning\n- (Yang et al., 2023) \"Test Accuracy vs. Generalization Gap\" instead of Martin & Mahoney, 2021a for review papers on robust metrics\n\nWe recognize that we did not summarize the 5+1 phases of learning. We include the following at line 363:\n\nIn a sequence of papers, Martin & Mahoney observe 6 classes of empirical behaviors in trained weight matrices, comprising a smooth transition from a random-like Marchenko-Pastur to a heavy-tailed density, before experiencing rank collapse. Excluding rank collapse, the five primary phases are:\n(a) Random-Like: Pure noise, modeled by a Marchenko-Pastur density.\n(b) Bleeding-Out: Some spikes occurring outside the bulk of density.\n(c) Bulk+Spikes: Spikes are distinct and separate from the Marchenko-Pastur bulk.\n(d) Bulk-Decay: Tails extend so that the support of the density is no longer finite.\n(e) Heavy-Tailed: The tails become more heavy-tailed, exhibiting the behavior of a (possibly truncated) power law.\nThe transition from (a) to (e) is also seen in (Thamm et al., 2022). This smooth transition between multiple phases is a primary motivation of this work. We find that this behavior is displayed by a combination of a nontrivial covariance matrix to capture the spikes, and the HTMP class with decreasing $\\kappa$.\n\nTo answer your questions:\n1. The depth and architectural effects of the structures are challenging to analyze due to the adopted general approach. These aspects can be studied empirically, necessitating further research. We admit that the model depth and other architectural details in this model will strongly affect the prior distribution of the feature matrix $\\Phi$ and also the final feature matrix. In this work, we want to provide a new random matrix model (high-temperature MP law) which may resemble the spectral properties of the feature matrix $\\Phi$. By tuning the parameters in HTMP, we can mimic the spectra of the feature matrices in different scenarios, which may be at initialization, may have been well trained, or may be associated with very complicated architectures. Although we do not have a clear picture of the relationship between the architectural details and spectra of the feature matrices, our random matrix model can potentially be viewed as an equivalent simplified model to study the feature matrices with very different architectural details.\n2. Yes, we are assuming that features equilibrate more rapidly than the loss. In this way, we are examining trajectories at the end of training (the kernel learning regime in (Fort et al., 2020) \"Deep learning vs. kernel learning\") once the features are mostly trained. We elaborate more on this in the current version.\n3. Interesting! So the input data is uniform, but this example still seems to assume a power law Fourier spectrum in the target function. Unless we are missing something, this still assumes heavy tails in the data (in the labels in this case). However, we have found that (Liao & Mahoney, 2021) \"Hessian Eigenspectra\" also shows how nonlinearities can influence the spectrum, although they do not prove scaling laws from this. In this case, conditions for the nonlinearity to exhibit heavy tails are still unclear, so we consider this a case of examining individual models.", "claims_and_evidence": "The theoretical claims are supported by both proofs and numerical experiments.", "methods_and_evaluation_criteria": "The authors focus on smaller scale vision models but this is very appropriate for a theoretical paper.", "theoretical_claims": "Yes, I read through many of the derivations in the Appendix.", "experimental_designs_or_analyses": "Yes, the experiments appear correct.", "relation_to_broader_scientific_literature": "This paper studies an important problem of the origin of neural scaling laws and provides an interesting universality hypothesis about their origin beyond the lazy learning regime.", "essential_references_not_discussed": "There are several important references that were not referenced but I think should possibly be included \n\n1. Li et al analyze how the hidden layer kernels in a Bayesian neural network change in a proportional limit where width and data diverge at fixed ratio https://journals.aps.org/prx/abstract/10.1103/PhysRevX.11.031059. \n2. Zavatone Veth et al also analyze the distribution over kernels after training in Bayesian linear networks in a variety of scaling limits https://ieeexplore.ieee.org/abstract/document/9723137\n3.  Thamm et al empirically study the spectra of trained weight matrices in deep networks https://arxiv.org/abs/2203.14661. Their results may also be relevant to the final feature matrices. \n4. Bordelon et al 2024 studied the dynamics of a random feature model concurrent with the work of Paquette et al 2024 and Lin et al 2024, also capturing compute optimal laws while also treating the effect of limited data in addition to limited training time or features https://openreview.net/pdf?id=nbOY1OmtRc.  A recent follow up from that group rederived their result using techniques from random matrix theory https://arxiv.org/abs/2502.05074 including $S$-transform methods.\n5. Bordelon et al 2024 had a follow up where they considered a simple model of rich (non-kernel) learning dynamics where the kernel is allowed to adapt that reproduced a faster scaling law in terms of source and capacity exponents https://arxiv.org/abs/2409.17858. This can be attributed to the changes to the kernel during optimization. I sense there could be a connection between their results and the theory in this work which allows the feature kernel to adapt beyond its prior.", "other_strengths_and_weaknesses": "The paper is technically strong and provides many experiments to support its claims. There are some remaining questions about how architectural details like depth and other earlier layers alter their master model (see questions below).", "other_comments_or_suggestions": "1. It could be useful to the reader to quickly define or outline each of the 5+1 phases either in the main text or in the supplementary material so that the reader would not need to consult prior works. \n2. In line 951 it says \"Appendix _\"", "questions_for_authors": "1. What is the role of model depth and other architectural details in this model (like widths of earlier layers)? Does all of this enter into the prior distribution over the feature matrix $\\pi(M)$? Would it also effect the rate at which the feature matrix can evolve during learning? What if an earlier layer had a significant bottleneck in width compared to the final feature layer. This would likely change the expressivity of the network and decrease the flexibility of the final feature matrix during optimization. \n2. In section 5.3, the authors consider optimization trajectories. However, the theory provides the distribution for the *final features.* Could the authors comment on the connection? Are they assuming that the features equilibrate more rapidly than the readout and thus the loss?\n3. The authors claim that other works rely on power law assumptions in the data space rather than feature space. This is not universally true, see for example section 5.1 here https://arxiv.org/abs/2409.17858 where the data is uniformly distributed but the nonlinearity in the network causes a power law decay which sets the scaling law on that task. Could the authors revise this?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ilmK2qfs0z", "reviewer": "Reviewer_PqrH", "summary": "Recent advancements in deep learning, including neural scaling laws, have highlighted the prevalence of heavy-tailed or power law behaviors in key network components such as the Jacobian, Hessian, and weight matrices. This phenomenon, termed heavy-tailed mechanistic universality (HT-MU), has been empirically linked to model performance, suggesting its fundamental role in deep learning success. To investigate the origins of this behavior, the study introduces a general class of random feature matrix models, the high-temperature inverse-Wishart ensemble. The model identifies three key factors contributing to heavy-tailed spectral densities: (i) complex correlation structures in data, (ii) lower training temperatures, and (iii) implicit bias in model structure leading to reduced eigenvector entropy. The study further explores the implications of HT-MU on learning dynamics, neural scaling laws, and optimizer trajectories.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive assessment of our work and for taking the time to verify the correctness of our proofs. We appreciate that the reviewer finds the experiments to be sound and suitably validate our model class. We agree about the high writing density; there are many individual components and motivations in our discussion. In line with other reviewer comments, we are considering alterations to the second section to better emphasize relevant details and move away from some of the statistical physics (see comments to Reviewer RfHj for example). We will include all experimental details in the main text by adding a separate experiment section with the following:\n\n5+1 phases of learning experimental details: For Figure 5, we train a MiniAlexNet for the CIFAR10 classification task with different batch sizes. This MiniAlexNet is a simplified version of AlexNet, which contains six layers: the first three are convolutional layers, followed by max-pooling layers, and the last three are fully connected layers. The histograms in Figure 5 are the histograms of the eigenvalues of $WW^T$ where $W$ is the trained weight matrix in the first fully-connected layer with input dimension $192 * 4 * 4$ and output dimension 1000. Here $W$ is initialized by the centered normal distribution with variance $\\sqrt{1/\\mathrm{fan_in}}$. For all the histograms in Figure 5, we trained the network using SGD with momentum, with a learning rate of 0.01 and a momentum parameter of\n0.9 for 200 epochs. For each histogram, we repeat the experiments 3 times for the average. The red dot curves are the numerical simulations of the density function of the HTMP with different $\\kappa$ and $\\gamma=1000/(192 * 4 * 4)$.\n\nTo obtain NTK spectral densities, we consider larger neural networks that are trained to near-zero loss (all $>99.8\\%$) on a subsampled dataset of 1000 entries through 200 epochs of a cosine annealing learning rate schedule with 200 epoch period, starting from a learning rate of 0.05 with a batch size of 64. Each model is comprised of the following number of parameters:\n- resnet9 (4.8M parameters)\n- resnet18 (11.1M parameters)\n- vgg11 (9.2M parameters)\n- vgg13 (9.4M parameters)\n- lenet (62K parameters)\n- logistic (30K parameters)\n- densenet121 (7.0M parameters)\nThe output layer of each model is altered from their ImageNet counterparts to classify with ten classes (for CIFAR-10, SVHN, and MNIST datasets). \n\nWe also recognize the relationship of the amount of training data to be a limitation. It is possible to complete the scaling law by including the power law dependence on the number of model features, although this differs from model size in a strict sense. Unfortunately, without a precise parameterization, it is difficult to establish such a law, but along with dependencies on individual model properties (e.g. depth, see comments to Reviewer Weri), we consider this of prominent interest in follow-up work.", "claims_and_evidence": "Yes, the claims made in the submission are supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "No, this paper doesn't propose any method or evaluation criteria.", "theoretical_claims": "Yes, I checked the correctness of the proofs for theoretical claims. Proposition 3.1 provides the density of 'optimal features,' which suggests that a stochastic optimizer concentrates on regions with high marginal likelihood (also called model evidence).", "experimental_designs_or_analyses": "Yes, I checked the soundness/validity and experimental designs or analyses. The numerical experiments of the paper effectively validates the results of its theoretical analysis.", "relation_to_broader_scientific_literature": "To uniformly analyze and understand the power-law phenomenon in deep learning, the paper proposes the high-temperature inverse-Wishart ensemble model.", "essential_references_not_discussed": "No, there are not related works that are essential to understanding the (context for) key contributions of the paper, but are not currently cited/discussed in the paper.", "other_strengths_and_weaknesses": "**Strengths:**\n1. This work proposes a novel high-temperature inverse-Wishart ensemble model, providing a unified theoretical framework to explain the emergence of the power-law phenomenon in deep learning.  \n\n2. The paper presents comprehensive numerical experiments, effectively empirically validating the effectiveness of its theory.  \n\n3. Overall, the paper is well-written and easy to follow.\n\n**Weaknesses:**\n1. The overall writing density of the paper is somewhat high. It is recommended to move some less important discussions on statistical physics to the appendix.  \n\n2. It is suggested that the authors add a separate experiment section in the main text to provide details on the experimental setup, conclusions, and their analysis.  \n\n3. Regarding the theoretical analysis of the neural scaling law, the authors seem to have only derived the power-law relationship concerning the amount of training data, without considering factors such as model size. It is recommended that the authors provide a more detailed discussion on the relationship between the high-temperature inverse-Wishart ensemble model and the neural scaling law.", "other_comments_or_suggestions": "See weaknesses.", "questions_for_authors": "See weaknesses.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The submission focuses on the origins of power-law spectra in deep learning by introducing a new framework based on the so-called high-temperature inverse-Wishart ensemble as a model for random features. The authors classify the possible sources of heavy-tailed spectra and characterise the resulting spectral properties within the proposed model. The manuscript is found to be technically strong and possibly an advancement in the study of the relevance heavy tailed feature matrices in machine learning. Although the received feedback is overall positive, reviewers complained about its dense mathematical content and not clear expostion. I strongly recommend the authors to tackle a review of the text by fixing the typo and address the concerns of Reviewer ```RfHj``` about the bibliography and the large number of citations of works by Mahoney and coworkers by considering their relevance for the text. I suggest *Weak accept*.
Decision: Accept (poster)