The Underlying Universal Statistical Structure of Natural Datasets
Authors: Noam Itzhak Levi
Yaron Oz
Venue: ICML 2025 poster (2025)
Keywords: Random Matrix Theory
Data Structure
Universality
Gaussian data Empirical Data Estimation
Power Law Scaling
Primary area: theory->everything_else
TL;DR: We analyze natural and synthetic datasets using RMT tools, finding that some universal properties are related to the strength of correlations in the feature distribution and connect these with the number of samples required to reach ergodicity.

Abstract: We study universal properties in real-world complex and synthetically generated datasets. Our approach is to analogize data to a physical system and employ tools from statistical physics and Random Matrix Theory (RMT) to reveal their underlying structure.
Examining the local and global eigenvalue statistics of feature-feature covariance matrices, we find: (i) bulk eigenvalue power-law scaling vastly differs between uncorrelated Gaussian and real-world data, (ii) this power law behavior is reproducible using Gaussian data with long-range correlations, (iii) all dataset types exhibit chaotic RMT universality, (iv) RMT statistics emerge at smaller dataset sizes than typical training sets, correlating with power-law convergence, (v) Shannon entropy correlates with RMT structure and requires fewer samples in strongly correlated datasets. These results suggest natural image Gram matrices can be approximated by Wishart random matrices with simple covariance structure, enabling rigorous analysis of neural network behavior.

The Underlying Universal Statistical Structure of Natural Datasets
Noam Levi' Yaron Oz?
Abstract
regarding the distribution. The second moment, however, cannot be set to O, and holds valuable information. This observation motivates the study of the empirical covariance (Gram) matrix, ΣM M XXT.
We study universal properties in real-world com- plex and synthetically generated datasets. Our approach is to analogize data to a physical sys- tem and employ tools from statistical physics and Random Matrix Theory (RMT) to reveal their underlying structure. Examinnig the local and global eigenvalue statisties of feature-feature co- variance matrices, we find: (i) bulk eigenvalue power-law scaling vastly differs between uncorre- lated Gaussian and real-world data, (ii) this power law behavior is reproducible using Gaussian data with long-range correlations, (iii) all dataset types exhibit chaotic RMT universality, (iv) RMT statis- ties emerge at smaller dataset sizes than typical training sets, correlating with power-law conver- gence, (v) Shannon entropy correlates with RMT structure and requires fewer samples in strongly correlated datasets. These results suggest natu- ral image Gram matrices can be approximated by Wishart random matrices with simple covariance structure, enabling rigorous analysis of neural net- work behavior.
The properties of ΣM in real world data are entirely un- known a priori, as we do not know how to parameterize the process which generated natural images. Nevertheless, in- teresting observations have been made. Empirical evidence shows that the spectrum of EM for various datasets can be separated into set of large eigenvalues (O(1)), a bulk of eigenvalues which decay power law X i I-o (Ruderman, 1997 Caponnetto and De Vito 2007) and large tail of small eigenvalues which terminates at some fi- nite index n. Since the top eigenvalues represent the largest overlapping properties across different samples, these are not simply interpreted without more information on the un- derlying distribution. The bulk of the eigenvalues, however, can be understood representing the correlation structure of different features amongst themselves, and has been key to understanding the emergence of neural scaling laws (Ka- plan et al. 2020: Maloney et al. 2022)
In this work, we study both the power law behaviors present in natural datasets, and their spectral statistics, with the goal of obtaining universal, analytically tractable model for real world Gram matrices, regardless of their origins. While this may not be feasible for any MM. fortunately, the standard datasets used today are high dimensional and contain many samples, a ubiquitous regime found in complex systems, and typically studied using Random Matrix Theory (RMT).
1. Introduction
Natural, or real-world, images are expected to follow some underlying distribution, which can be arbitrarily complex, and to which we have no direct access to. This distribution could have infinitely many nonzero moments, with varying relative importance compared to one another. In practice, we only have access to a very small subset of samples from the underlying distribution, which can be parameterized as X dd where d is the dimension of each image vector and M is the number of samples. The first moment of the data can always be set to 0, since we can remove the mean from each sample, without losing information
RMT a powerful tool for describing the spectral statisties of complex systems. It is particularly useful for systems that are chaotic but also have certain coherent structures. The theory predicts universal statistical properties, provided that the underlying matrix ensemble is large enough to suffi- ciently fill the space of all matrices with a given symmetry, a property known as ergodicity (Guhr et al. 1998). Ergod- icity has been observed in a variety of systems, including chaotic quantum systems (Bohigas et al. 1984 Mehta 1991] Pandey 1983), financial markets, nuclear physics and many others (Plerou et al. 1999; Brody [988 Efetov] 1997). To demonstrate that a similar universal structure is also ob- served for correlation matrices resulting from datasets, we will employ several diagnostic tools widely used in the field
'École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland : Raymond and Beverly Sackler School of Physics and Astronomy, Tel-Aviv University. Tel-Aviv 69978, Isracl. Cor- respondence to: Noam Levi tnoam.levi@cpfl....
Proceedings of the nnnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

The Underlying Universal Statistical Structure of Natural Datasets
of quantum chaos. We will analyze the global and local spectral statistics of empirical covariance matrices gener- ated from three classes of datasets: (i) Data generated by sampling from a normal distribution with specific corre- lation structure for its features, (ii) Uncorrelated Gaussian Data (UGD), (iii) Real-world datasets composed of images, at varying levels of complexity and resolution. Our research aims answer the following questions:
2. Background and Related Work
Neural Scaling Laws Neural scaling laws are a set of em- pirical observations that describe the relationship between the size of a neural network, dataset, compute power, and its performance. These laws were first proposed by Kaplan et al. (220) and have since been confirmed by number of other studies (Maloney et al.) 2022 Hernandez et al. 2022) and studied further in (Ivgi et al.
2022; Alabdulmohsin el al.
2022 Sharma and Kaplan 2022, Sorscher et al. 2022 De
Is a power-law spectrum a universal property across real-world datasets?; what determines the scaling exponent and what properties should an analytical model of the dataset have in order to follow the same behavior?
howski 2023] Fernandes et al. 2023). The main finding of neural scaling laws is that the test loss of neural network scales power-law with the number of parameters in the network. This means that doubling the number of param- eters roughly reduces the test loss by 2 However, this relationship does not persist indefinitely, and there is point of diminishing returns beyond which increasing the number of parameters does not lead to significant improvements in performance. One of the key challenges in understanding neural scaling laws is the complex nature of the networks themselves. The behavior of a neural network (NN) is gov- erned by large number of interacting parameters, making it difficult to identify the underlying mechanisms that give rise to the observed scaling behavior, and many advances have been made by appealing to the RMT framework.
What universal properties of datasets can be gleaned from the empirical covariance matrix and how are they related to local and global statistical properties of RMT?
How to quantify the extent to which complex data is well characterized by its Gram matrix?
What, if any, are the relations between datasets power laws, entropy and statistical chaos diagnosties?
Our primary contributions are:
We find that power-law spectra appears across various datasets. It is governed by a single scaling exponent a, and its origin is the strength of correlations in the underlying population matrix !. We accurately recover the behavior of the eigenvalue bulk of real-world datasets using Wishart matrices with the singular values of a Toeplitz matrix (Gray) (006) as its covariance. We dub these Correlated Gaussian Datasets (CGDs).
Random Matrix Theory RMT is a branch of mathemat- ics that was originally developed to study the properties of large matrices with random entries. It is particularly suited to studying numerous realizations of the same system, where the number of realizations M co, the dimensions of the system d + oo, and the ratio between the two tends to a constant d/M + S 1, E R Results from RMT calculations have been applied to a wide range of problems in Machine Learning (ML), beyond the scope of neural scaling laws, including the study of nonlinear ridge regres- sion (Pennington and Worah 2017]. random Fourier feature regression (Liao et al. 2021 ) the Hessian spectrum (Liao) and Mahoney| 2021), and weight statistics (Martin and Ma- honey 2019: Thamm et al.) [022). For a review of some of the recent developments, we refer the reader to Couillet and 1ao (022) and references therein.
2. We show that generically, the bulk of eigenvalues' dis- tribution and spacings well described by RMT pre- dictions, verified by diagnostic tools typically used for quantum chaotic systems. This means that the CGD model is a correct proxy for real-world data covariances.
3. We find that the effective convergence of the empirical covariance matrix as a function of the number of sam- ples correlates with the corresponding RMT description becoming good description of the statistics and the eigenvalues power law decay.
Universality Considerable work has been dedicated to the concept of universality, i.e. that certain features are shared between seemingly disparate systems, when the systems are sufficiently large. For instance, spectra that are generated by different dynamical processes may have similar distribu- tions (Bao et al. 2013; Baik et al.) 2004 Hu and Lu| 2022 Bai and Silverstein 2010). Universality is powerful since it often happens that System A's complex structure is diffi- cult to analyze, and can be explained by system B, which lies in the same universality class, and is much easier to study. In our work, system A represents weal-world datasets
4.
The Shannon entropy is correlated with the local RMT structure and the eigenvalues behavior, and is substan- tially smaller in strongly correlated datasets compared to uncorrelated data. Additionally, it requires fewer samples reach the distribution entropy.
There are systems which display multiple correlation scales, showing several bulk exponents (Levi and Oz 2023).
2

---

The Underlying Universal Statistical Structure of Natural Datasets
10a
CeNeb 10.5]
ImageNutt4 5) Shapes (0.4)
CIFAR
10a
Fashion 25)
(:1)
MN ST (0.25
CGO (0.5)
101
10a
10"
10a
<
10:2-
10-2
100%
corruption
80%
corruption
10"
10
40%
corruption
corruption
FMNIST no comuptiin
10
10a
10:
10a 10a 10*
Figure 1. Left: Scree plot of Egj. for several different vision datasets, as well as for UGD and CGD with fixed a. Here, the number of samples is taken to be the entire dataset for each real-world dataset, and M = 50k for the gaussian data, where we set c = We see a clear power law for the cigenvalue bulk as X; x i i-1-a where all real-world datascts display S 1/2. Right: The power-law exponent value can be tuned from a = 1/4 to a = -1 by corrupting the FMNIST dataset with a varying amount of normally distributed noise.
with unknown statistics generated from a complex process, while system B is our CGD, whose Gram matrix is a simple Wishart matrix. The fact that real world datasets fall in the same universality class as CGD allows us to replace its complex covariance matrix by the simple CGD one, while retaining the information encoded in its spectrum.
transformations. In order to study Eq. ① using the RMT approach, we define Σij, as a single sample realization of the population random matrix ensemble Eyy. and thus Eij.M is the empirical ensemble average, i.e. Σm = (Ea)aeM = MI >n Di approximating the limits of . + 0o, d + co. If M and d are sufficiently large, the statistical properties of Mas will be determined entirely by the underlying symmetry of the ensemble. We refer this case the "RMT regime".
3. Correlations and power-law spectra
In this section, we analyze the feature-feature covariance matrix for datasets of varying size, complexity, and origin. We consider real-world as well as correlated and uncorre- lated gaussian datasets, establish a power-law decay of their eigenvalues, and relate it to a correlation length.
3.2. Data Exploration
We study the following real-world datasets: MNIST (LeCun) et al. 2010), FMNIST (Xiao et al.) 2017). CIFARIO (cif)
IinyIMMAGENE (Torralba et al.) 2008), and CelebA (Liu
et al. 2015) (downsampeld to 109 89 in grayscale). We proceed to center and normalize all the datasets in the pre-processing stage, to remove the uninformative mean contribution. The uncorrelated gaussian data is repre- sented by a data matrix whose elements X10 E R, where each column is drawn from a jointly normal distribution (0, Iaxa). We then construct the empirical covariance matrix Σn = 1 ΣM XiaXaj E Rdxd To generate cor- 32
3.1. Seature-Feature Empirical Covariance Matrix
We consider the data matrix X Xia E RdxM, con- structed of M columns, each corresponding to single sample, composed of d features. In this work, we focus on the empirical feature-feature covariance matrix, defined as
M
Eijas
Σ M
1 XiaX E Rdxd
(1)
related gaussian data, we repeat the same process, changing the sample distribution (0 Eaxa), where we choose a specific form for Σ which produces feature-feature correla- tions and includes natural cut-off scale, as
Intuitively, the correlations between the different input fea- tures, Xia. should be the leading order characteristic of the dataset. For instance, if the Xio are pixels of an image, we may expect that different pixels will vary similarly across similar images. Conversely, the mean value of an input fea- ture is uninformative, and so we will assume that our data is centered in pre-processing stage.
Toe ij = S.j. Ts = Lij + +cj - 11" = (U*SV); (2)
where a, c E R. The matrix D is a positive semi definite diagonal matrix of singular values S constructed from T, a full-band Toeplitz matrix. The sign of a dictates whether correlations decay (negative) or intensify (positive) with distance along a one-dimensional feature space?)
random matrix ensemble is a probability distribution on the set of x d matrices that satisfy specific symmetry properties, such as invariance under rotations or unitary
CCorrelation strength which grows with distance hallmark of some one-dimensional physical systems, such as the Coulomb
3

---

The Underlying Universal Statistical Structure of Natural Datasets
3.3. Correlations Determine the Spectral Noise to Data Transition
of results for the specific matrix ensemble which both real-world and Gaussian datasets converge.
We begin by reproducing and extending some of the results from Maloney et al. (2022). In Fig. [ we show the M eigenvalue power law decay for the different classes of data (i.e. real-world, UGD and CGDs). We find that for all datasets, the eigenvalues bulk scales as a power-law
We interpret Mar for real world data as single realization, drawn from the space of all possible Gram matrices which could be constructed from sampling the underlying pop- ulation distribution. In that sense, Ex itself is a random matrix with an unknown distribution. For such a random matrix, there are several universality classes, which depend the strength of correlations in the underlying distribu- tion. These range from extremely strong correlations, which over-constrain the system and lead to the so called Poisson ensemble (Atas et al., 2013). to the case of no correlations, which is equivalent to sampling independent elements from normal distribution, represented by the Gaussian Orthogo- nal Ensemble (GOE) (Mehtal 2004). These classes are the only ones allowed by the symmetry of the matrix XXT.. provided that the number of samples and the number of features are both large. Since the onset of the RMT regime depends on the population statistics, it is a priori unknown. Determining if real data covariances converge to an RMT class, and to which one they converge to at finite sample size would inform the correct way model real-world co- variances.
A x xio a E R. = doulk.
(3)
where i - 10 is approximately where the power law be- havior begins and dbulk is the effective bulk size, where the power-law abruptly terminates. We stress that this be- havior repeats across all datasets, regardless of origin and complexity.
The value of a can be readily explained in terms of correla- tions within our CGD model. Taking the Laplace Transform of the second term in Eq. (2) the bulk spectrum is given by Appendix B as
1+a (
hulk = c. r(a+1)
(4)
where [(x) is the Gamma function. This implies that the value of a determines the strength of correlations in the original data covariance matrix. For real-world data, we consistently find that 0, which corresponds to increas- ing correlations between different features. In contrast, for UGD, the value of a -1, and the power-law behavior vanishes. Interpolating between UGD, and real-world-data, the CGD produces power-law scaling, which can be tuned from -1 < a S 0, the case of decaying long range correlations, or S a < oc for increasing correlations, to match any real-world dataset we examined. Lastly, we can extend this statement further and verify the transition from correlated to uncorrelated features by corrupting real-world dataset (FMNIST) and observing the continuous deterioration of the power-law from a ~ 1/4 to a = 11, implying that the CGD can mimic the bulk behavior of both clean and corrupted data.
Below we review the tools used in our analysis. While we provide an overview of each diagnostic, we refer the reader to Tao 2012); Kim et al. (2023) for a more comprehen- sivereview. We then apply these tools to gain insights into the statistical structure of the datasets.
Spectral Density: The empirical spectral density of a matrix Σ is defined as,
1 - Σ 8(A (2)),
PS(A)
(5)
where 5 is the Dirac delta function, and the A,(E), 1,..., n. denote the n cigenvalues of E, including multi- plicity. The limiting spectral density is defined the limit of Eq. as + D0.
Level Spacing Distribution and r-statistics: The level spacing distribution measures the probability density for two adjacent eigenvalues to be in the spectral distance s, in units of the mean level spacing A. The procedure for normaliz- ing all distances in terms of the local mean level spacing is often referred to as unfolding. We unfold the spectrum of the empirical covariance matrix M (p) by standard meth- ods (Kim et al. 2023). reviewed in Appendix A. Ultimately, the transformation X; + = p();) is performed such that e, shows an approximately uniform distribution with unit mean level spacing. Once unfolded, the level spacing is given by s; = c++1 - e, and its probability density function p(s) is measured.
4. Global and Local Statistical Structure
4.1. Random Matrix Theory
In this section, we move on from the eigenvalue scaling be- havior, to their statistical properties. We begin by describing the RMT diagnostic tools, often used to characterize RMT ensembles, with which we obtain our main results. We de- fine the matrix ensemble under investigation, then provide an overview of each diagnostic, concluding with a summary
(p)
and Riesz gases (Lee and Yang 1966 (Smorodinsky 1953), which display an inverse pawer-law repulsion, while decaying correla- tions are common in the 1-d Ising model [Ising 1925).
4

---

The Underlying Universal Statistical Structure of Natural Datasets
10
10a
10
10
10a
10'
10'
10
10"4
10
10
10-1
10-1
10
10-1
27
114
10-a
3T
2/5
10:
y7w
FPRCST
10:31
CGD CIFAR10
10
ImiogeNet
PMNST (30000 samp eal
152000 samplea)
10"
Imageliet 150000 samples)
10
10
10
10a
10
10
10
10
10?
10
10a
20a
FMN ST
UED
USD
CCD RMNIST, (KL=6.001)
HP
CIFARLO
CGD CIPARIO, (KL-0.200)
102
magetlet
Imagehot, (SL=2.01)
10a
10a
P())
P(A)
P(A)
100
100
10c
10-1
10-1
0.0
0.2
0.4
0.6
0.8
1.0
1.2
10
0.0
0.2
0.4
0.6
0,8
1.0
12
14
1.6
0.0
0.2
0.4
0.6
0.5
10
1.2
1.4
1.6
Figure 2. Top row: Scree plot of ExgMr for several different configurations and datasets. We show the cigenvalues of the population covariance matrix Io the eigenvalues for the empirical covariance of the full real-world dataset with M = 50k and finally the eigenvalues the empirical covariance using the same yToe, with M = 50k. The datasets used here are (left to right): FMNIST, CIFARIO, ImageNet. Bottom row: Spectral density for the bulk of eigenvalues for the same datasets, as well as comparison against UGD of the same dimensions. The À indicates normalization over the maximal eigenvalue among the bulk. We also provide the KI. divergence between the CGDs and the real-world data distributions.
The distribution p(s) captures information about the short- range spectral correlations, demonstrating the presence of level repulsion, i.e., whether p(s) + 0 as s - O, which is a common trait of the GOE ensemble, as the probability of two eigenvalues being exactly degenerate zero. Further- more, the level spacing distribution p(s) for certain systems known. For integrable systems, it follows the Poisson distribution p(s) e while for chaotic systems (GOE), it is given by the Wigner surmise
The expectation value of the ratios ry takes very specific values if the energy levels are the eigenvalues of random matrices: for matrices in the GOE the ratio is (r) x 0.536.
The value is smaller for integrable systems, approaching (r) Ro 0.386 for a Poisson process (Atas et al. 2013)
Spectral Form Factor: The spectral form factor (SFF) is a long-range observable that probes the agreement of a given unfolded spectrum with RMT at energy scales much larger than the mean level spacing. It can be used to detect spectral rigidity, which is a signature of the RMT regime.
ps(s)= ee b28"
where 3, Zz. and by depend on which universality class of random matrices the covariance matrix belongs to (Mehta 2004). In this work, we focus on matrices that fall under the universality class of the GOE, for which | 1, as we show that both real-world data and CGD covariances belong to.
The SFF is defined as the Fourier transform of the spectral
2- 2-point correlation function (Cotler et al. 2017 Liu 2018)
EΣ ples)e ( 2 (8) K(r) = [(+)) Z
While the level spacing distribution depends on unfolding the eigenspectrum, which is only heuristically defined and has some arbitrariness, is useful to have additional diagnos- tics of chaotic behavior that bypass the unfolding procedure. The r-statisties, first introduced in Oganesyan and Huse (2007). is such diagnostic tool for short-range correlations, defined without the need to unfold the spectrum.
where 7(r) = Tre- Tre- The second equality is the nu- merically evaluated SFF (Juntajs et al. 2020), where e, is the unfolded spectrum, and Z = Ip(e))) is chosen to ensure that K(T + oo) R 1.
The SFF has been computed analytically for the GOE en- semble, and it reads
Given the level spacings 8;, defined as the differences be- tween adjacent eigenvalues-.. ... < A; < A1+1 < ... without unfolding, one defines the following ratios:
Koor(T) (T) = 2r rIn(1 + 27) for 0 < 1, KGog(T) = 2 - rln + 1) - 1)] for T
(9)
Several universal features occur in chaotic RMT ensem- bles, manifesting Eq. (9) and discussed in detail in Liu
r Min(s;, s 1)/ax(ss, s1+1), 05 (7)
5

---

The Underlying Universal Statistical Structure of Natural Datasets
2.5
GGE
UGD THNIST CIFARIO
1.0
Wigner Surmise
spectral form factor
2.0
CCD CIFARIO
CCD PHNIST
0.8
CGD CIFAR10
CGD FHNIST
10'
1.5
UGD FMNIST CIFARIO
CGD
PANIST
109
OGn
CIPAR30
1.0
USD
20-1
0.5
0.2
FMN
-0-22
GOE
0.0
0.2
0.8 1.0 0.0 0.5
1. 3.0 3.5 4.0
10-4
10 100
10a
r
unfolded
7
Figure 3. The probability density (left), the unfolded level spacing distribution (center) and the spectral form factor (right) of EM for FMNIST, CIFARIO, their CGDs, and UGD, obtained with M = 50000, Black curves indicate the RMT predictions for the GOB distribution from Eq. (II). These results indicate that the bulk of real-world data cigenvalues belongs to the GOE universality class, and that system has enough statistics to converge to the RMT predictions.
(2018); Kim et al. (2023). We mention here only two: (i) The constancy of K ((T) for r 2 1 is simply a consequence of the discreteness of the spectrum. (ii) The existence of a timescale that characterizes the ergodicity of a dynami- cal system. is defined as the time when the SFF of the dynamical system converges the universal RMT compu- tation. More concretely, it is indicated by the onset of the universal linear ramp 2r as in equation 9 which is absent in non-ergodic systems.
In Fig. we show that the CGDs capture not only the power law decay of the eigenvalue bulk, but also the spectral density and the distribution of eigenvalues, for ImageNet, CIFARIO, and FMNIST, measured by the Kullback-Leibler divergence (KL) (Kullback and Leibler [195]). We further emphasize this point by contrasting the distributions with the MP distribution, which accurately captures the spectral density of the UGD datasets. This measurement alone is in- sufficient to determine that the system is well approximated by RMT, and we must study other statistical diagnostics.
4.2. Insights from the Global and Local Statistics
4.2.2. LEVEL SPACING DIAGNOSTICS
4.2.1. EIGENVALUES DISTRIBUTIONS
RMT predicts that certain local and global statistical prop- erties are determined uniquely by symmetry. Therefore, the empirical covariance matrix must lie either in the GOE ensemble if it is akin to a quantum chaotic system or in the Poisson ensemble, if it corresponds to an integrable system.
While the power law behavior of the bulk of eigenvalues is certainly meaningful, it is not the only piece of information that can he extracted from the empirical covariance matrix. Particularly, is natural to inquire whether the origin of the power-law scaling determines also the degeneracy of each eigenvalue. We can test this hypothesis by comparing the global and local statistics of the bulk between real-world data and their CGD counterparts.
Both the level spacing and r statistics (the ratio of adjacent level spacings) probability distribution functions and SFF for a Wishart matrix in the limit of d, M + oo and d/M 7/ are given by the GOE universality class:
For the gaussian datasets we generate, there are known pre- dictions for the spectral density, level spacing distribution, r-statistics and spectral form factor. In these special cases, the empirical covariance matrix in Eq. (T) is known as a M).
PGOE (s) = se Is? (F) 44-223 2/3, (11)
PGOE(r) = 27 (r +2)
Wishart matrix (Wishart 1928): : Eig.M
4 (1 (2) O - -r),
The spectral density p()) of Wishart matrix is given by the generalized Marčenko- Pasttur (MP) law (Silverstein and
In Fig. we demonstrate that the bulk of eigenvalues for various real-world datasets behaves as the energy eigenval- ues of a quantum chaotic system described by the GOE universality class. This result matched by both the UGD and the CGDs, as expected of a Wishart matrix. Here, the dataset size is taken to be M 50000 samples, and the results show that this sample size is sufficient to provide proper sampling of the underlying ensemble.
Bai 1995; Couillet and Liao. 2022), which depends on the details of and specified in App. B for certain limits. For x² o²Ta Ta, the spectral density is given explicitly by the MP distribution as
(Amax 2)(A Amin)
p(A) 2TO2
(10)
7t
for A E [Amin Amax] and 0 otherwise. Here, a E
Large random real symmetric matrices helong in the orthogo- nally invariant class.
Amax/min = 21 vr)², = d/M and d, M -> cc.
6

---

The Underlying Universal Statistical Structure of Natural Datasets
4.2.3. HIGHER ORDER STATISTICS
CIFAR-10
Going beyond the nearest-neighbor local statistics, one can probe structure at increasingly larger scales. In particular, the neighbor local statistics can be defined as
1.75
n-1
n=2
6
1.50
1.25
As+2n Ai+n
1.00
(12)
\itn X
0.75
for which the GOE distribution is known analytically (Tekur
0.50
et al. 2018
0.25
(r + 2 PGOB (10)) 7 (1 + r + 1+32/2
0.00
3
4
(13)
p(m)
GOE data
where = =n 33 and Z is a normalization factor. -
In Fig. we show that the nearest-neighhor statistics are in good agreement with the GOE theory, but as the value of n gets larger, farther neighbors are probed and discrepan- cies between the distributions emerge. Probing these higher order local statistics in greater detail can offer system- atic path towards explaining the discrepancies between the expected learning curves for neural network performance trained on gaussian vs. real world data (Loureiro et al.] 2021), since ((") statistics explicitly measure deviations from gaussianity.
2.00
n=1
o=4
1.75
1.50
1.25
1.00
0.75
0.50
0.25
0.00
1
3
p(n)
4.2.4. EFFECTIVE CONVERGENCE
Figure 4. Top: The ((()) statistics for CIFAR1O. Bottom: The r()) statistics for GOE data. The local statisties at diffcrent scales (n) shown in different colors, while the black curves indicate the theoretical curves. For small n = 1, 2, 4 the ,(()) distributions match between real-world data and the GOE results, but at larger scales n = 6 a deviation begins to appear.
Having confirmed that CGDs provide good proxy for the bulk structure for large fixed dataset size, we may now ask how the statistical results depend on the number of samples.
As discussed in Sec. 3.1
Σs can be interpreted as an ensem- ble average over single realizations of the true population covariance matrix E. As the number of realizations M in- creases, a threshold value of Merit is expected to appear when the space of matrices that matches the effective dimen- sion of the true population matrix is fully explored.
where TGOE = 4 - 2-3 0.536 for the GOE.
Next, we compare the results obtained for Merit from 8(M) to the ones obtained from the global statistics by using spectral distance measure for the eigenvalue bulk given by
The specific value of Merit can be approximated without knowing the true effective dimension by considering two different evaluation metrics. Firstly, convergence of the lo- cal statistics of Σn,, given by the point at which its level spacing distribution and r value approximately match their respective RMT ensemble expectations. Secondly, conver- gence of the global spectral statistics, both of Ens to that of Σ and of the empirical parameter os to its population expectation
JaM al = A(M),
(15)
where aus is the measured value obtained by fitting a power- law to the bulk of eigenvalues for fixed dataset size M, while a represents the convergent value including all sam- ples from dataset.
Lastly, we compare the empirical Gram matrix EM with the convergent result Σ. obtained using the full dataset by taking
Here, we define these metrics and measure them for different datasets, obtaining analytical expectations for the CGDs, which accurately mimic their real-world counterparts.
Σs Σ <(M)EE)
(16)
where |A| is the spectral norm of A, and eMM will be our measure of the distance between the two covariances.
We deduce Meris from the local statisties by measuring the difference between the empirical average value and the theoretical one given by
In Fig. 6
we show the results for each of these metrics separately as a function of the number of samples M. We find that the 8(M) parameter, which a measure of local
ITM TRMMTI 8(M)rGO.,
(14)
7

---

The Underlying Universal Statistical Structure of Natural Datasets
1.0
MNIST
MN
0---
FMN ST
CIFARIO
0.
CIFARIO
FMNNST
FNNNSST
CCD PHIUST
A(MM/A
0.6-
-6(4))
UGD
0.4-
0.4
Mort
ST
CIFARIO
0.2
FINIIT
0.
0.0 0.5
10
0.5
0.
Muamplea
Ha/HH
HydH
Figure 5. Convergence of the various metrics in Eqs. 14 to 16 in relation to entropy for the bulk of cigenvalues. Left: The Shannon entropy H as a function of dataset size M. Center: Convergence the normalized a metric AM/ / to its asymptotic value as a function of the normalized entropy HM/ H. Right: Convergence of the normalized statisties metric dar / its asymptotic value as function of the normalized entropy H H. We show the results for CIFARIO, FMNIST, MNIST, UGD, and the FMNIST CGD
NNNIST
MNIST
0.175
12-
FARLO
CED PRNIST
CGD FHNIST
12:2
0.125
UCO
UGD
4000
AUM)
0.035
Ment
Mont
0.088
0.05
10'
10'
Mesmptes
Mesmpiex
Mcamp
The a distance metric (M) for the bulk of eigenvalues. Right: The full matrix comparison metric c(M). We show the results for CIFARIO. FMNIST, UGD, and the FMNIST CGD as a function of the number of samples. The results show that the bulk distances decrease as 1/M, where M is the number of samples, asymptoting to a constant value at similar values of Merit ~ d (black dashed), where is number of features.
Figure 6. Left: The r distance metric 8(M) for the bulk of eigenvalues. Center:
statistics, converges to the expected GOE value at roughly the same Marit as the entirely independent A (M) parameter, which measures the scaling exponent a. The combination of these metrics confirms empirically that the system has become ergodic at sample sizes roughly Merit d, which is much smaller than the typical size of the datasets.
4.2.6. ENTROPY, SCALING EXPONENT AND RMT
In Fig. 5 (left), we see that the entropy saturation is corre- lated with the effective convergence in Fig. as a function of the number of samples, while the middle and right plots show the correlation between the convergence of the entropy, the scaling exponent, and the r-statistics, respectively. We see that real data and gaussian data with the same scaling exponent exhibit similar convergence behavior.
4.2.5. DATASETS ENTROPY
The Shannon entropy Shhannon 1948) of a random variable a measure of information, inherent to the variable's possible outcomes [Rényi 1956). given by H - Pt log(p:) where pi is the probability of given outcome and is the number of possible states. For covariance matrices, we define p: given the spectrum as p: \// > D"bik X; where Vbulk is the number of bulk eigenvalues.
4.3. Eigenvector Phenomenology
The results of our work concern the eigenvalue behavior of natural dataset covariances, as these seem to be well cap- tured by a universal RMT description. The eigenvectors, however, specify the non-universal basis, which is highly dataset dependent, and should not be captured by a ran- dom orthogonal matrix. Nevertheless, for completeness, we explore the phenomenology of covariance eigenvectors in App. E Our main observations are: (i) the eigenvectors obtained from an empirical covariance Σx align, on aver- age, with the population as M1/2 for small sample sizes and then exponentially with M (Fig. [2]. and the transi-
In Fig. 5 (left) we plot the Shannon entropies of real and gaussian datasets as function of the number of samples. The entropies grow linearly and reach plateau whose value is related to the correlation strength, with strong correlation corresponding to low entropy. We see the same entropy for both the gaussian and real datasets that have the same scaling exponent, implying that they also share the same eigenvalues degeneracy.
We omit UGD from the center panel, as = -1 for any M.
8

---

through con- The Underlying Universal Statistical Structure of Natural Datasets
tion point Mye is significantly smaller than the eigenvalue threshold Merit: (ii) the eigenvectors corresponding to the largest eigenvalues, as well as the smallest ones, align faster than the bulk of eigenvectors (Figs. 13 and [4), in contrast with the GCD model (Fig. 15 a property that seems to be related to how spatially localized these eigenvectors are. These results imply that accomplishing learning tasks for which the spectrum insufficient and eigenvector properties needed, requires different scale of data samples.
eigenvectors in neural networks, as both components criti- cally shape information processing.
Finally, while we empirically observe chaotic spectral prop- erties in real-world data, the origin of this chaos remains unclear. Is it intrinsic to the strongly correlated structure of natural data, or does it arise from noise the sampling process? Resolving this trolled studies of synthetic datasets with tunable correlation and noise-wi deepen our understanding of data complex- ity.
5. Conclusions
In this work, we demonstrate that the bulk eigenvalue spec- tra of Gram matrices for real-world data can be accurately modeled by a Wishart matrix with a shift-invariant corre- lation structure and defining exponent a. Crucially, these Gram matrices universally conform to the GOE statistics across various data- gett-generating processes. This universality indicates that our approximation captures not only the scal- ing behavior of eigenvalues but also their full distribution, which can be rigorously derived using tools from RMT.
A Acknowledgements
We would like to thank Yohai Bar-Sinai, Marat Freytsis, Alex Maloney, Dan Roberts, and Jamie Sully for useful discussions and comments. N.L. would like to thank the Milner Foundation for the award of a Milner Fellowship. The work of N.L. is supported by the A14Science program at EPFL. The work of Y.O. is supported in part by Israel Science Foundation Center of Excellence. This work was performed in part at Aspen Center for Physics, which is supported by the U.S. National Science Foundation grant PHY-2210452.
Our findings bridge a critical gap between theoretical mod- els of data-often used study neural network (NN) behav- ior-and real-world datasets. While Normally distributed data has advanced our understanding of NN learning dynam- ics (Gerace et al.) 2021 Mei and Montanari 2020 dAscoli et al. 2020), parameter scaling limits (Maloney et al. (0022) and weight evolution (Arous et al. 2018), such models in- herently lack the structured correlations present in natural data. As a concrete application, accurate modeling of spectral densities-enabler by our framework-is essential for solving the dynamics of even simple teacher-studen models with correlated inputs.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal con- sequences of our work, in particular due to the large model sizes considered in this work, but we do not feel there are specific aspects of this work with broader impacts beyond the considerations relevant to all large machine learning models.
show in App. [ that
We propose an RMT-based data model that preserves ana- lytical tractability while incorporating realistic correlations. Although similar assumptions appear implicitly in neural scaling law studies, our work strengthens this paradigm with rigorous mathematical foundations, empowering practition- ers to derive predictions better aligned with real data.
References
Daniel I.. Ruderman. Origins of scaling in natural images. Vision Research, 37(23)3385-3398, 1997. ISSN 0042- 6989. doi: 10.1016/500422
Our results also constrain the statistical properties of real- world data distributions. Specifically, the convergence of Gram matrices to GOE (rather than Poisson) statistics under finite sampling imposes non-trivial moment conditions on the underlying data distribution. This insight could inform both the inference of data-gent generetenrrtnngggg processes and the de- sign of synthetic datasets that replicate natural correlations.
A. Caponnetto and E. De Vito. Optimal rates for the regu- larized least-squares algorithm. Found Comput Math, 7 (2):331-368, 2007. doi: 10. 1000000/12000000000111
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models, 2020.
Although our analysis focuses on Gram matrices, which in- herently discard spatial relationships, this simplicity enables a broad applicability. Extending our framework to language datasets, audio signals, or other modalities may reveal uni- versal power-law behaviors transcending data types. Future work should explore the interplay between eigenvalues and
Alexander Maloney, Daniel A. Roberts, and James Sully. A solvable model of neural scaling laws, 2022.
Thomas Guhr, Axel Groeling, and Hans A. Weidenmüller. Random-mant theories in quantum

---

The Underlying Universal Statistical Structure of Natural Datasets
physies: common concepts. Physics Reports, (99(4- 00088-4. URL
Utkarsh Sharma and Jared Kaplan. Scaling laws from the data manifold dimension. J. Mach. Learn. Res., 23: 9:1-9:34, 2022.
jun 1998. doi: 10. 10.1016660370--5533007
httpt::://mmxrrg///appppapee///p v23/20-1111 html
RIS0370-157332283322000000
O. Bohigas, M. J. Giannoni, and C. Schmit. Characteriza- tion of chaotic quantum spectra and universality of level fluctuation laws. Physical review letters, 52(1):1, 1984.
Ben Sorscher, Robert Geirhos, Shashank Shekhar, Surya Ganguli, and Ari Morcos. Beyond neural scaling laws: beating power law scaling via data pruning. In NeurIPS, 2022.
M. L. Mehta. Random matrices, volume 111. Academic Press, 1991.
Lukasz Debowski. A simplistic model of neural scal- ing laws: Multiperiodic santa fe processes. CoRR, abs/2302.09049, 2023. doi: 10.485550ar2xww22302.
A. Pandey. Random matrix theory and quantum chaos. Reviews of Modern Physics, 55(4))807-823, 1983.
09049. URL https:///ol/ 0ttp://///...rr1104455550
arXiv. arXiv.2302.c 09049
V. Plerou, P. Gopikrishnan, B. Rosenow, I. A. N. Amaral, H. E. Stanley, and Stanley M. S. Random matrix theory and financial markets. Physical Review E, 60(5):6519- 6532, 1999.
Patrick Fernandes, Behrooz Ghorbani. Xavier Garcia, Markus Freitag, and Orhan Firat. Scaling laws for multilingual neural machine translation. CoRR, 0oss2302.096500. 2023. doi: 10.48550VarXiv.2302.
T. A. Brody. Random matrix models in nuclear physics. Reports on Progress in Physics, 44(4):1125-1199) 1981.
09650. URL https: bttps://doi... /ow. https///doi.org//0..0 org/10. 48550/
arXiv. FrXiv.2302.0 09650
K. B. Efetov. Supersymmetry and disorder in quantum mechanics. Cambridge University Press, 1997.
Jeffrey Pennington and Pratik Worah. Nonlinear random matrix theory for deep learning. In Advances in Neural Information Processing Systems, pages 2637-2646, 2017.
Noam Levi and Yaron Oz. The universal statistical structure and scaling laws of chaos and turbulence, 2023.
Zhenyu Liao, Romain Couillet, and Michael W Mahoney. A random matrix analysis of random fourier features: beyond the gaussian kernel, a precise phase transition, and the corresponding double descent. Journal of Statistical Mechanics: Theory and Experiment, 2021 dec 2021. doi:
Robert M. Gray. Toeplitz and circulant matrices: A re- view. Foundations and Trends@ in Communications and Information Theory, 233) 2006. ISSN
1567-2190. doi: 10.1561/01000000006. URL http: 7/dx. wtw.o/.org/1..1110111111.
Zhenyu Liao and Michael W. Mahoney. Hessian eigenspec- tra of more realistic nonlinear models, 2021.
Danny Hernandez, Tom Brown, Tom Conerly, Nova Das- Sarma, Dawn Drain, Sheer El-Showk, Nelson Elhage, Zac Hatfield-Dodds, Tom Henighan, Tristan Hume, Scott Johnston, Ben Mann, Chris Olah, Catherine Olsson, Dario Amodei, Nicholas Joseph, Jared Kaplan, and Sam Mc- Candlish. Scaling laws and interpretability of learning from repeated data, 2022.
Charles H Martin and Michael W Mahoney. Traditional and heavy-tailed self-regularizzation in neural network models. arXiv preprint
2019. URL 08276]
Matthias Thamm, Max Staats, and Bernd Rosenow, Random matrix analysis of deep neural network weight matrices. Phys. Rev. E, 106:054124, Nov 2022. doi: 10.1103/ PhysRevF. 10.. URL. 10t.
Maor Ivgi, Yair Carmon, and Jonathan Berant. Scaling laws under the microscope: Predicting transformer perfor- mance from small scale experiments. In Yoav Goldberg, Zomitsaa Kozareva, and Yue Zhang, editors, Findings of the Association for Computational Linguistics: EMNLP 2022, Abu Dhabi, United Arab Emirates, December 7-11, 2022, pages 7354-7371. Association for Computational Linguistics, 2022. URL. /ettps://aclanteee ogg org/2022. findings www.nndiingseennp.....m 544
Romain Couillet and Zhenyu Liao. Random Matrix Methods for Machine Learning. Cambridge University Press, 2022. doi: 10. 0000//9911091190 10011 28490.
Zhigang Bao, Guangming Pan, and Wang Zhou. Univer- sality for the largest eigenvalue of sample covariance matrices with general population. The Annals of Statis- tics, 43(1), feb 2015. doi: 10. 10.12144144aos1238 URL https: /ttp::////oo..oor1101111111
Ibrahim M. Alabdulmohsin, Behnam Neyshabur, and Xiao- hua Zhai. Revisiting neural scaling laws in language and vision. In NewrlPS, 2022.
121
10

---

The Underlying Universal Statistical Structure of Natural Datasets
Jinho Baik, Gerard Ben Arous, and Sandrine Peche. Phase transition of the largest eigenvalue for non-null complex sample covariance matrices, 2004.
Vadim Oganesyan and David A. Huse. Localization of interacting fermions at high temperature. Phys. Rev. B, 75:155111, 2007. doi: 10.1103/PhyseevB.75 15$111. URL htt.ppp://iix.apps.ooo 1 HDD3/PhysRevB.c 75. 75.151.11111
Hong Hu and Yue M. Lu. Universality laws for high- dimensional learning with random features, 2022.
Jordan S. Cotler, Guy Gur-Ari, Masanori Hanada, Joseph Polchinski, Phil Saad, Stephen H. Shenker, Douglas Stanford, Alexandre Streicher, and Masaki Tezuka, Black holes and random matrices. J. High Energy Phys., 2017:118, 2017. doi: 10. N0.10077/HEEP05200)1188 URL. httphhttt:///itt wpr /tt com/article. 10. (200//EEE00
Zhidong Bai and Jack W Silverstein. Spectral analysis of large dimensional random matrices, volume 20, Springer, 2010.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Pierre Haffner. The mnist database of handwritten digits. 2010.
Han Xiao, Salim Rasul, and Richard S Zemel. Fashion- mnist: novel image classification benchmark based on fashion articles. arXiv preprint 170. 2017.
Junyu Liu. Spectral form factors and late time quantum chaos. Physical Review D, 98(8), oct 2018. doi: 10.1103/
physrevd.98.0860222 URL httpps://ooi org/10 110342Fphyszevd. 98 086026
Cifar- 10. URL https: www.. c8 wwwonto...roooooooooooooooo.
J. Juntajs, J. Bonca, T. Prosen, and L. Vidmar. Quan- tum chaos challenges many-body localization. Phys. Rev. E, 102:062144, 2020. doi: 10.1103/PhysRewE.100 062144. URL http::tttsss....eeeeme http:///i.ikaaas..../////////// 1
Pkriz/ www.c.ff..hheeeeee
Antonio Torralba, Andreas A Efros, and Christopher An- derson. Tiny imagenet: A benchmark for evaluation of image classification algorithms. International Journal of Computer Vision, 69(2):203-2 220000222.2222 2008.
org/dbi/10
102.
John Wishart Wishart. Generalised product moment distri- bution samples from an indefinitely large population. Biometrika, 20(1-2)33--5, 1928.
Ziwei Liu, Zihang Luo, Xiaogang Wang, and Xiaoou Tang. Celeba: A large-scale celebrity face attribute dataset. http://mmlab.ie wwt.ollla............ cww.....keee..e..e....... cetu.h//prnjeett/Ceeee..... 2015.
Jack W Silverstein and Z. D. Bai. On the empirical dis- tribution of eigenvalues of class of large dimensional random matrices. Journal of Multivariate Analysis, 54 (2):175-199, 1995.
D. Lee and C. N. Yang. Statistical mechanics of charged particles. Zeitschrift für Physik, 196:433 453, 1966. doi: 10.1007//F002550005.
M. A. Smorodinsky. On the classical motion of charged particles. Journal of Mathematical Physics, 4:1005-1011, 1953. doi: 10.1063//.1703711
Solomon Kullback and Richard A Leibler. Information the- ory and statistics. The Annals of Mathematical Statistics, 22(1):79-86, 1951.
Ernst Ising. Beitrag zur theorie des ferromagnetismus. Zeitschrifi für Physik, 31:853-865, 1925. doi: 10.1007/ BF01343133.
S. Harshini Tekur, Udaysinh T. Bhosale, and M. S. San- thanam. Higher-order spacing ratios in random matrix theory and complex quantum systems. Physical Re- view B, 98( (8(10), September 2018. ISSN 2469-9969. doi:
Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux. Distribution of the ratio of consecutive level spac- ings in random matrix ensembles. Phys. Rev. Lett. 110:084101, 2013. doi: 10. 10.1103/hhysReveett..10 084101. URL https::::///ikk aps. org.//////00./00//0///////// 1103/PhysRevLett.o 110.
10. 10.110//pyyseeebb.8..oo 10430. URL http:/ //tp::///xxx// /txp//xxooo/ o org/10. 1ttp/1..113//Phyaeeee
98 104305
Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mezard, and Lenka Zde- borová. Learning curves of generic features maps for re- alistic datasets with a teacher-student model. In Advances in Neural Information Processing Systems, volume 34, 2021.
M. L. Mchta. Random Matrices. 3 edition, 2004.
Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012.
Claude Elwood Shannon. A mathematical theory of commu- nication. Bell System Technical Journal, 27(3):379 423, 1948.
Joonho Kim, Yaron Oz, and Dario Rosa. Quantum chaos and circuit parameter optimization. Journal of Statistical Mechanics: Theory and Experiment, (2:0230044 feb 2023. doi: 1T0.008881112224444220 URLp:eeoooooooLLLLtRtaaaaat.aaa.
Alfred Rényi. On measures of information and entropy. Proceedings of the Fourth Berkeley Symposium on Math- ematical Statistics and Probability, 144):547-561 1956.
URL https:77
doi wtt.o////111111821111111.

---

The Underlying Universal Statistical Structure of Natural Datasets
Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mezard, and Lenka Zdeborova. Generalisation error in learning with random features and the hidden man- ifold model. Journal of Statistical Mechanics: Theory and Experiment, (12) 24013, dec 2021. doi: 10. 1088/1742 URL
Hyunjune Sebastian Seung, Haim Sompolinsky, and Naftali Tishby. Statistical mechanics of learning from examples. Physical review A, 45(8):6056, 1992.
Timothy LH Watkin, Albrecht Rau, and Michael Biehl. The statistical mechanics of learning a rule. Reviews of Mod- ern Physics, 65(2):499, 1993.
10.
Andreas Engel and Christian Van den Broeck. Statistical me- chanics of learning. Cambridge University Press, 2001.
Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptoties and double descent curve, 2020.
David L Donoho. Statistical modeling of natural images with wavelets. Proceedings of the National Academy of Sciences, 92(12):111--199 1995.
Stephane dAscoli, Maria Refinetti, Giulio Biroli, and Flo- rent Krzakala. Double trouble in double descent: Bias and variance (s) in the lazy regime. In International Con- ference Machine Learning, pages 2290. PMLR, 2020.
Noureddine El Karoui, Derek Bean, Peter Bickel, Chingh- way Lim, and Bin Yu. On robust regression with high- dimensional predictors. Proceedings of the National Academy of Sciences, 11033))14557-1456) 2013.
Gerard Ben Arous, Song Mei, Andrea Montanari, and Mihai Nica. The landscape of the spiked tensor model, 2018.
Andrew M Saxe, James L McClelland, and Surya Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. In International Conference on Learning Representations, 2014.
Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75(4):667-766, 2022.
Lenka Zdeborov's and Florent Krzakala. Statistical physics of inference: Thresholds and algorithms. Advances in Physics, 65(5):453-552 2016.
Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mezard, and Lenka Zdeborova. Generalisation error in learning with random features and the hidden manifold model. In Proceedings of the 37th International Confer- ence on Machine Learning, pages 3452-344 -362.. PMLR, 2020.
David Donoho and Andrea Montanari. High dimensional robust m-estimation: Asymptotic variance via approxi- mate message passing. Probability Theory and Related Fields, 166(3-4):935-966)) 2016.
Sebastian Goldt, Bruno Loureiro, Galen Reeves, Florent Krzakala, Marc Mezard, and Lenka Zdeborova. The gaussian equivalence of generative models for learning with shallow neural networks. In Proceedings of the 2nd Mathematical and Scientific Machine Learning Confer- ence, pages 426-4711 PMLR, 2022.
Song Mei and Andrea Montanari. The generalization error random features regression: Precise asymptotics and double descent curve. arXiv preprint urXiv:1908.05330 2019.
Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. arXiv preprint 1033::990 2019.
Alexander Wei, Wei Hu, and Jacob Steinhardt. More than a toy: Random matrix models predict how real-world neu- ral representations generalize, In Proceedings of the 39th International Conference on Machine Learning, pages 23549-23552 23588. PMLR, 2022.
Emmanuel J Cand'es Pragya Sur, al. The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression. The Annals of Statistics, 48(1):27-42) 2020.
Federica Gerace, Florent Krzakala, Bruno Loureiro, Lu- dovic Stephan, and Lenka Zdeborová. Gaussian universal- ity of perceptrons with random labels. 133003, 2023.
Benjamin Aubin, Florent Krzakala, Yue M Lu, and Lenka Zdeborov'a a. Generalization enror in high-dimensional perceptrons: Approaching bayes eπтоr with convex opti- mization. In Advances in Neural Information Processing Systems, volume 33, 2020.
Luca Pesce, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan. Are gaussian data all you need? extents and limits of universality in high-dimensional generalized linear estimation. In Hal Daumé III and Aarti Singh, edi- tors, Proceedings of The 40th International Conference on Machine Learning, volume 162, pages 10-15. PMLR, 2023.
Fariborz Salchi, Ehsan Abbasi, and Babak Hassibi. The per- formance analysis of generalized margin maximizers on separable data. In International Conference on Machine Learning, pages 8417-8426. PMLR, 2020.
12

---

The Underlying Universal Statistical Structure of Natural Datasets
Michael A Nielsen.
A simple formula for the erage gate fidelity of quantum dynamical opera- tion. Physics Letters A, 303(4):249-252, October 2002. ISSN 0375-9601. doi: 10. 10.1016660335--6000022 URL htt.::tpp///doo//// doi. ort//0..1/1110.. S0375 9601 (02)01272-0 0
01272-0.
Mohamed El Amine Seddik, Cosme Louart, Mohamed Tamaazousti, and Romain Couillet. Random matrix the- ory proves that deep learning representations of gan-data behave as gaussian mixtures, 2020.
Andrew M. Saxe, James L. McClelland, and Surya Ganguli. A mathematical theory of semantic development in deep neural networks. Proceedings of the National Academy of Sciences, 116(23) May 2019. ISSN 1091-6490. doi: (0..0773/ppass URL p
//dx. //A4a..o...rrr//0..000/////....2
13

---

The Underlying Universal Statistical Structure of Natural Datasets
A. The Unfolding Procedure
Here, we provide additional details on the unfolding procedure used to produce Fig. 3 in the main text.
Care must be taken when analyzing the eigenvalues of the empirical covariance matrix MM since they exhibit unavoidable numerical errors. To control for the effect of numerical errors, adopted robust phenomenological procedure that utilizes the fact that all eigenvalues of EM must be non-vanishing by definition. To ensure we consider only eigenvalues of ΣM unimpacted by edge effects, we inspect only the bulk spectrum.
Restricting to the bulk removes many eigenvalues of E as many are zero for small M. However, for larger M when Ma 's structure is clearly visible, this is not the case. The procedure ensures the eigenvalues kept are robust and not significantly impacted by numerical precision. From the significant eigenvalues of the empirical covariance matrix EM. we compute the spectrum A.
Spectrum
Unfolded Spectrum
250
200
P()))
P(e))
0.8
150
0.6
100
0.4
50-
0.2
0.000 0.005 0.010 0.015 0,020 0.025 0.030 0.035
0.040
0.0
0.2
0.4
0.6
0.8
1.0
A;
ej
Figure 7. Bulk eigenvalue distribution for the empirical covariance matrix constructed from 50000 samples of FMNIST, before unfolding (left), and after unfolding (right) The unfolded spectrum displays approximately unit mean, and defined on the interval [o, 1].
The unfolding procedure used to derive the unfolded spectrum is as follows:
Arrange the non-degenerate eigenvalues, A of the empirical covariance matrix Σar) in ascending order.
Compute the staircase function S(A) that enumerates all eigenstates of the empirical covariance matrix (Ea) whose eigenvalues are smaller than equal to X.
3. Fit smooth curve, denoted by p(A). to the staircase function. Specifically, we used 12th-order polynomial as the smooth approximation.
4. Rescale the eigenvalues A as follows:
A + + ((A,)
(17)
5. By construction, the unfolded eigenvalues e; should show an approximately uniform distribution with mean level spacing
1. This can be used to check if the procedure was successful by plotting the unfolded levels and checking the flatness of the distribution.
In Fig. 7 we show an example of the unfolding procedure for the FMNIST dataset. Specifically, we show the eigenvalue distribution before (P():) ) and after (P(e.) unfolding. Up to the quality of the smoothing function (()s), the unfolded eigenvalue distribution displays a uniform distribution on the unit interval.
B. Spectral Density for Wishart Matrices with a Correlated Features
For (usppp((s)) the Stieltjes transform G and inverse Stieljes transform ps: are defined as
G(z) = (z)() dt 2 -
"T[TT(E - =I,) Ps(A) lim 3G(A + ic),
(18)
c-
14

---

The Underlying Universal Statistical Structure of Natural Datasets
where E| ] is taken with respect to the random variable X and (2 =I,,)= the resolvent of Σ.
For the construction, discussed in the main text, and general a, there is no closed form for the spectral density. However,
in certain limits, analytical expressions can be derived from the Stieljes transform using Eq. (18). Specifically, given
determinstic expression for Σ, the spectral density can be derived by evoking Theorem 2.6 found in Couillet and Liao (2022)
which uses the following result by Silverstein and Bai 19955)
1 G(2) G(z) (())=-(()) 7
1-7 ((2) a(z)=
(-2+
1 Ti [T(=++++()))))) -1
+
(19)
72
M
where y = d/M and d, M + oc, and we substitute C from the original theorem with Σ.
The empirical covariance matrix of the Gaussian correlated datasets discussed in the main text, Wishart matrix with a deterministic covariance, and thus fits the requirements of Theorem 2.6, where E S. S Tu Li ci jl In order to use Eq. (1),, it is useful to first find the singular values of Tig. This can be done by using the discrete Laplace transform (extension of the Fourier transform), leading to
VtTU, and
(T) = S(s = - cLi-a (e 3) ce (c -< 2c-2--0, d),
(20)
where s 1...d.. d, (2x, k a) is the Lerch transcendent, and Li(x) is the Poly-log function. Note that by the definition of S, Eq. 20) is non-negative function of s. Because the identity matrix commutes with Iooe, may substitute Eq.
20
in Eq. 19 to obtain
1
G(z) -2 + {Sa(a)
(21)
where we define the sum Sa(a) to be
Σ (s)
Sa(a) Σ
(22)
1 + (T))
Since the behavior of y. (s) consider the case of a < -1, where correlations decay very quickly. In this scenario, the covariance matrix reduces to Toe (8) = 1.
(s) is intrinsically different for positive and very negative a, separate the two cases. First,
Here, Sa(a) is given simply by the a + oo limit
Sa(a + -0) Σ
1
(23)
G(2)
1 G(2)
which is precisely the case of Σ Ia.
Solving Eq. 21) using the above result yields the following expression for G(2)
G(2) -1 -1- (--1 7/ 42
(24)
2z
Finally, substituting Eq. 24 into Eq. 18) leads to the known Marčenko-Pastur (MP) law (Couillet and Liao 2022)
(Amax 2)(A Amin)
P(A) = 2T
for X E [ Amax] and 0 otherwise,
(25)
7d
where \max//min (1 2
The other interesting limit that of -1, in which the correlations do not decay quickly, and for - oc the Laplace transform of the Toeplitz matrix simplifies to
(o)
d
Ita
Toe (s) ~ 1(11 ((1+a)
(26)
15

---

The Underlying Universal Statistical Structure of Natural Datasets
Eigenvalues Distribution for FMNIST
FMNIST
10a
CGD FMNIST
Analytics
10
p(A)
10-1
10-2
0.0
0.5
1.0
1.5 A
2.0
2.5
3.0
Figure 8. Thcoretical predictions for the bulk spectral density of CGD matrix against the empirical densities of FMNIST and the CGD, The green curve represents the generalized MP distribution, given by the solution to the inverse Steiljes transform in Eq. The CGD curve has value of e 1.14, and y 380/ 1000 since that is the approximate number of bulk eigenvalues.
Using this approximation for the population covariance in Eq. 22
we obtain
d
-
cssd)) 1-a
-1
d(=) = -2
6= cI a..
(27)
(())(((s)
In the d - co limit,we can convert the integral using the Riemann definition
lim Σ (a,))x - [dd((x), Ar
b -- 1 d d'
lim
(i/d)
(28)
allowing us write the equation for G(2) as
d
((s/)) 1-a
1 G(=)èx-1-0 dx ~ -2+7 So
cr -1-1
-1
G(2) = - +
(29)
d
+ ((z))(s/d))
==1
2Fi 1, otoiii + भो cG(+)
-1
+2
(())
where 2Fi (a, by c; 2) is the Gaussian hypergeometric function. Eq. 29 is an algebraic equation which can be solved numerically, analytically approximated in certain limits.
In Fig. show the theoretical results for the spectral density of a Wishart matrix with Too covariance, for value of and matching FMNIST, against the empirical densities for FMNIST and the matching CGD. The green curve shows the generalized MP distribution given by the inverse Stieljes transform of Eq. 29
C. Robustness of the results
Here, we discuss some details regarding the robustness of our local and global statistical analyses.
16

---

The Underlying Universal Statistical Structure of Natural Datasets
For all of our analyses, focused the full Gram matrix, consisting of every sample in given dataset. This implies that only have access a single realization of ΣM empirical Gram matrix, per dataset, thus limiting our ability perform standard statistics, for instance averaging over an ensemble of EM. and obtaining confidence bands. This is not issue in the RMT regime, as the matrix itself is thought of as ensemble on to itself, and its eigenvalues have an interesting structure due to a generalization of the Central Limit Theorem (CI.T).
r statistics for varying M
statistics, rea FMNIST vs population E
GOE prediction
12
FNMSST
CGD FMNIST
1.50
Population Singular Values
1.25
p(
0.50
0.25
10-1
0.4
0.8
10
Figure 9. Left: The statistics distribution for FashionMNIST, comparing M 50000 with M = 1000 subsets. In the first case, we obtain only single realization of the Gram matrix, and the statistics appear noisy, however, when taking 40 realizations of smaller subset, still above Mexit, we see that the fit to the GOE prediction (green) improves. We will add these figures, either in the main appendices, including gondness of fit measures on the of the datasets studied in the paper. Right: The statistics distribution for MNSST and its CGD. In red, we show the singular values of the population covariance, y Toe used in the main text. In Orange, the true FMNIST distribution, obtained by taking 40 different realizations of M 1000 subset of the full dataset, leading perfect fit to the GOE prediction (blue). In green, we show the CGD using 1000 samples as well. This figure illustrates that the deterministic population covariance does not sufficiently capture the information that resides the Gram matrix, while the CGD does.
Log-Log Plot Eigenvalues (Descending Order)
Eigenvalues Distribution FMNIST
10a
FMNIST
10.
Empirical OGD PHNIST
CGD FMNEST
Population CGD FMNIST
Population
10a
10"
10
10
10
10
10a 10a
10a 10
0.0
0.5
1.0
1.5
2.0
2.5
Figure 10. Left: Scree plot for the cigenvalues of the FashionMNIST Gram matrix (bluc), its CGD (orange) using M 1000 for 50 runs, and Toeplitz population covariance matrix (green). Here, we show the population precisely in spectral scaling. The Gram matrices for FMNIST and its analogue are obtained by first normalizing the samples (mean subtraction and dividing by the standard deviation) and population covariance is rescaled by a constant factor that depends only on the input dimension d. Right: The eigenvalue distribution for FashionMNIST, Gram matrix (blue), its CGD (orange) using M 1000 for 50 runs, and the Toeplitz population covariance matrix (green). see certain local features (such the spacing between cigenvalues) is poorly captured by the deterministic population covariance.
empirical covariance matrices match
three distributions are similar, as can expected, that
We can still attempt to persuade the reader that our results are robust posteriori, by noting that the number of samples
17

---

The Underlying Universal Statistical Structure of Natural Datasets
required reach the RMT regime is approximately Merit d. This implies that we can consider sub-samples of the full empirical Gram matrix, consisting of Marit, as MMan 1 E Da=1" Naj and average over multiple sub-sample
Meris
matrices.
In Fig. 9 we see an implementation of this process for FMNIST, demonstrating that additional sampling pushes the distribution to a perfect fit for the GOE r-statistics, while in Fig. 10 we see the same type of convergence for the eigenvalue distribution.
D. Results on Augmented Data
In Fig. we show the RMT results for the eigenvalue spectra as well as the r-statistics for augmented versions of CI and FaShionMNISTT. The augmentations include shifts, rescalings and rotations, in this case taken with a range of 50% for height and width shift, rotation angle of 40 degrees and zooming range up to 20%. The results indicate that global statisitics can change substantially, as seen in the the left scree plot of Fig. I1 while the local r-statistics is maintained even for highly augmented images.
E. Eigenvector Phenomenology
In this section, we discuss the universal eigenvector behavior of two natural datasets: CIFARIO and FashionMNIST. We compute the population (full training set size) and the empirical covariances Σ, ΣM. and define their eigen-decomposition simpye as
MM = SUS
(30)
where U is the matrix of cigenvectors and S is a diagonal matrix composed of the eigenvalues of MM.. We further employ the magnitude of the cosine similarity to measure how aligned two eigenvectors are with one another, as
[Cosine similarity = Sc.i i7. izy.
(31)
which ranges between [-1, 1], as the eigenvectors are normalized.
Analyzing the eigenvector alignment with the number of samples, we observe two interesting phenomena that seem to be universal. First, in Fig. [2] we show that for both datasets, the average magnitude Sg behaves differently for small sample sizes, below and above some threshold value Meec- Concretely, for M < Mree the similarity scales as Sc xx M while for M > Maecr it scales as Sc x M/Mpor, where Moop is the population sample size, which is Mpop 50k for both datasets.
Secondly, in Fig. 13 we see clear distinction between the eigenvector convergence behavior for a random basis (bottom rows) and natural data. The GCD model is always given in a random basis, whose eigenvectors align to the population cigenvectors in descending order, from the ones associated with the largest eigenvalues, to the ones associated with the smallest eigenvalues. A different phenomenon occurs for natural data both the top eigenvectors and the smallest eigenvectors align together, while the bulk is aligned last. We believe this is tightly related the localization properties of the eigenvectors. In Figs. [4] and 15 see that the cosine similarity changes can be qualitatively traced the the Inverse Participation Ratio (IPR), defined as
IPR: =
Σ ||us.
(32)
where indicates the index of the eigenvector, while a denotes the different elements of that particular eigenvector. A large small IPR implies that the vector is highly de-localized, since each element should scale as 1/Vd and since there are d elements, IPR 1/d, while for fully localized vector, all elements apart from one will be O, therefore we have that the IPR 1. We see in Figs. 14 and 15 that for natural data the IPR is small for the top eigenvectors, then plateaus before growing for the smallest eigenvectors, implying transition from de-localization to localization. For CGD data, which has a random basis, the exact opposite occurs, where the first top eigenvectors are highly localized, and the rest are not.
18

---

The Underlying Universal Statistical Structure of Natural Datasets
Fstatistics comparison for CIFAR-10
r-statistics comparison for Fashion MNIST
0.7
10a
Cifar10
fashion, mnist
Augmented Cifar10
0.6
Augmented fashion mnist
10a
GOE
prediction
0.5
10a
0.5
10a
P(r)
P(r)
10-1
0.3
0.3
Cifar10
0.2
10-3
Augmented Gffr10
10
fashion_nee mnist
0.1
0.1
Augmented fashion mnist
10
r-statistics
r-statistics
Figure 11. Top: Scree plot and statistics distributions for CIFARIO and FashionMNIST, well augmented variations of these datasets for = 50k samples out of the datascts. Augmentations encompass shifts, rescalings and rotations, in this case taken with range of 50% for height and width shift, rotation angle of 40 degrees and zooming range up 20%. Here,
the definition of (1422) This can easily translated to form presented the main text. Bottom: Representative images taken from the 4 types datasets, in order from top row to bottom row are: CIFARIO, Augmented CIFAR10, 10, FasHionMNIST and Augmented NaSHionMNIST. The results show that the local statistical structure is maintained for the augmented parameters chosen, as as power behavior the augmented and original data. It should be clear that if the parameters are exaggerated, for instance if the zoom range grows to the level constant gradient across the image, we not expect the statistical structure to remain the same.
s S1-1, which has known GOE prediction given by P(r) (1+7++3)877-
eigenvalue bulk, with different scaling exponent between
F. Universality in Neural Network Analysis A Toy Example
Universality laws have been employed in various ways study error universality in neural networks, for instance in (Mei and Montanari 2022 Gerace et al., 2020 Goldt al. 2022 Hu and Lu 2022). In this context, looks directly the universality of the training and generalisation errors instead of the features, taking into account the labels and the task. It has also been observed to hold for correlated Gaussian data in teacher-student settings in (Loureiro et al. 202T). Furthermore, it has been shown that for a simple regression task, the computation of the error reduces to an RMT problem (Wei et al.) 2022). which is linked to the work presented in the main text regarding the data features themselves. particular, it has been noted that some cases the structure of the bulk fully characterises the error, even for multi-modal distributions, see (Gerace
19

---

The Underlying Universal Statistical Structure of Natural Datasets
1.0
Studaa Deate
1-x
0.8
4-0M0 A
Cusine Sumiacks
0.8
Cosine Siralacity
Cosine Similaciy
0.5
0.6
0.6
0.4
0.4
0.4
0.2
0.2
00
00-
0.0
Sample Size M
Sample Size M
Sample Size M
Figure 12. Average magnitude the cosine similarity between population eigenvectors and their empirical estimate from samples. Left: The results for CIFAR 10 eigenvectors. Middle: The results for FashionMNIST eigenvectors. Right: The results for CGD (a = 0.3) eigenvectors. The results quite universally show an increase in cosine similarity that scales roughly as power law (M 1/2 for real data and M1/3 for CGD) for M below some threshold value Moee: while the increase from this value and above grows exponentially.
et al. 2023 Pesce et al., 2023
As self contained example of applying universality in neural network analysis, we call upon an exceedingly simple machine learning setup, namely, training a linear network using gradient descent to learn a teacher-student mapping. We show that even in this basic example, it is necessary to apply the results demonstrated in the main text in order to correctly analyze the system, when data correlations are present in the underlying population covariance.
Teacher-student models have been the subject of long line of works (Seung et al. 1992 Watkin al.
1993 Engel and
Van den Broeck 2001 Donoho 1995 El Karoui et al. 2013 Saxe et al. 2014 Zdeborov and Krzakala
2016 Donoho
and Montanari 2016) and have experienced resurgence of interest in recent years (Mei and Montanari 2019 Hastie
et al. 2019 Cand'e es et al. 2020 Aubin et al.) 2020 Salehi et al. 2020) as a powerful tool to study the high- sigg--dimensonall asymptotic performance of learning problems with synthetic data.
The teacher-student model can be described as follows: The teacher uses ground truth information along with a probabilistic model generate data which is then passed to the student who is supposed recover the ground truth as well as possible only from the knowledge of the data and the model.
Here, we consider linear teacher-student model, where the data inputs E IR are identical independently distributed (iid) normal variables drawn from Gaussian distribution with non-trivial population covariance z; (O, Epop). We draw Nrr training samples, and the teacher model generates output labels by computing a vector product on each input y = ==- r, where ws E Rdd., assuming a perfect, noiseless teacher. The student, which shares the same model as the teacher, generates predictions I = w x, where w E Re as well. The loss function which measures convergence of the student to the teacher outputs is the standard MSE loss. Our analysis done in the regime of large input dimension and large sample size, i.e., dm, Nur co, where the ratio X = dna/ /Nee E R+ kept constant. The student model is trained with the full batch Gradient Descent (GD) optimizer for steps with learning rate n, The training loss function is given by
N N
Cur
(33)
where we define A = w - wn as the difference between the student and teacher vectors. Here, Σn, 1
N. Σmn x. is the dn x di empirical data covariance, or Gram matrix for the training set. The elements of w² and w are drawn at initialization from normal distribution Wo, w" ~ ~,(00, 1/(2d..). We do not include biases in the student or teacher weights, as they have no effect on centrally distributed data.
The generalization loss function is defined as its expectation value over the input distribution, which can be approximated by the empirical average over Ngen randomly sampled points
Lgen ES-NN [][(== - w =- Tr Egnn
(34)
Here Egen is the covariance of the generalization distribution. Note that in practice the generalization loss is computed by a sample average over an independent set, which is not equal to the analytical expectation value. The gradient descent
20

---

The Underlying Universal Statistical Structure of Natural Datasets
N 1023
N 2000
N 300
N
5000
230
Conine Similariy
400
030
402
601
030
60-
2
290
200
1040
258
N 10000
20000
N 40000
N 60000
80
Cosime Sim Hariiey
603
000
07-
601
000
60-
803
290
700
1840
2-
1040
0.0
M=
1022
M 2.00
M 3000
50003
130-
100-
200
20
201
200
60
Craimnc Similrity
K0-
co0-
544
000
57-
720
100000
20000
10000
5/00000
130
20c
Conimm Stanuuriy
201
XX--
201
500
420
41
730
TuE-
700
100
200
40
630
201
480
9-10
200
46
60)
36
+00
N
1023
N 2000
N -000
N
5000
201-
330
202
Chosto Samiaccis
407-
4,3
60-
0.0
200
200
1000
334
3000
(50
30
75
1000
550
2:00
130
N
10000
20000
N 40000
N
50000
1.0
201
202
K12
400-
02
40-
601-
0,4
220
280
250
100
2.4
240
244
3000
(50
20
1000
520
200 150
130
0.0
Figure 13. Top rows: Covariance eigenvector cosine similarity magnitude matrices function of the number of samples M, for CIFARIO. Middle rows: Covariance eigenvector cosine similarity matrices as a function of the number of samples M. for FashionMNIST. Bottom 2 rows: Covariance eigenvector cosine similarity matrices function of the number of samples M, for GCD with a = 0.3. Clearly, a random hasis aligns by first matching the cigenvectors associated with the largest cigenvalues, while natural data aligns first both the largest and the smallest eigenvectors.
21

---

The Underlying Universal Statistical Structure of Natural Datasets
1x 10
X-346
180.0
5- 080
1007
IPR
311
X=-000
S=300
10
0.003
0.002
2
1.84
1.20
Cueiat Saniaatly
Costae Simlacty
0.6
06
0.4
0.19
0.0
000
1000
Index
Index
0.0
10"1
31
JPR:
M 2000
3000
34 10000
M 30000
31 90000
X 50000
M 9200
34
14
10000
M (000)
0.:3
10
0.0
1.0
0.5
Coelax Sualleitt
Cooitr Siualaitt
20
11.0
0022
2041
J00
12.1
700
Indlx
Inciee
Figure 14. Cosine similarity for between each population eigenvector and its empirical estimate from N samples, well IPR for each cigenvector function of number of samples. Top row: IPR in log and lincar scale for CIFARIO. Second row: Absolute cosine similarity for CIFARIO. Third row: IPR log and linear scale for FashionMNISt. Fourth row: Absolute cosine similarity for FashionMNIST.
22

---

The Underlying Universal Statistical Structure of Natural Datasets
1.0
0.5
1PR
30000
1PR
N- 2200
X- 100
0.4
0.4
0858
Coaine Sinilntt
Cosine Similarit
0.6
20m
0.2
Index
Index
Figure 15. Cosine similarity for between each population eigenvector and its empirical estimate from N samples, well as IPR for each eigenvector as function of number of samples. Top row: IPR in log and linear scale for CGD with a 0.3. Bottom row: Absolute cosine similarity log and linear scale for CGD with 0.3.
equations at training step t are
= =1 (I-1= 21) A..
(35)
where y E + is the weight decay parameter, and TRd xdm is the identity.
Eq. 35
can be solved in the gradient flow limit, setting 1) 7ndt and di 0, resulting in
Á(!) 2no (())
2n. Ao. A(I)
(36)
where Ag is simply the difference between teacher and student vectors at initialization. It follows that the empirical losses, calculated over dataset admit closed form expressions as
Cur
A An
Cgen
D/e
End
Ppop
2no Dn.
(37)
Since the directions of both and the eigenvectors of En are uniformly distributed, we make the approximation that the projection of A on all eigenvectors is the same, which transforms Eq. (37) to the simple form
Lar R 1 -E
Anpeec -nnpeibys Vi
(38)
din
while the calculation for the generalization loss amounts to
Leen
Σ den i.
1
(39)
where U is a random unitary matrix used to rotate to the basis of Exr-
Now we turn to the choice Epop and the implication for Σ,,, As demonstrated in the main text, the empirical covariance matrix of many real world data-sets can faithfully modelled by a Wishart matrix with long range correlations, where the
bulk of eigenvalues is described by the population covariance Epop [(1 + (i//d) -1-a we discussed in App. B
23

---

The Underlying Universal Statistical Structure of Natural Datasets
can utilize RMT observations to give a closed formula expression for the empirical Gram matrix eigenvalues and spectral density, in terms of generalized MP law.
Following this choice of data modelling, and focusing on the bulk eigenvalues alone, it is clear that the sums Eqs. 38 and (39) are the empirical averages over the function e 4-uut fv)), if v follows the spectral density derived in App. We
can the solve the training dynamics by approximating the sum by its respective expectation value,
(To X, a, Re Eypopppo (A.a)
[uee - 4mprrt
(40)
In order to proceed further for the generalization loss, we note that the rotation matrices which form the basis for Etr are random unitary matrices, drawn from the Haar measure. This implies that we can glean further information by averaging over training realizations, which will not change the training trajectory at all, but will provide with an average generalization loss (C)o. We utilize the following property of ensemble averaging over unitary random matrices (Nielsen 2002)
! du(U)UXUt
(xx) Eu EUXUX
Tr(XXI,
(41)
din
where du((') is the Haar measure. Since the eigenvalue distribution does not change upon this averaging, the average generalization loss can be expressed as
(Cyen)u Fo Tr(Epop) x din
Σ
Anout
(42)
which can be approximated by its expectation value
(Egen (No A, a, U) (e (2.a)) [e 4nprt)
(43)
completing the dynamical analysis of the loss curves for the model hand.
Above gave toy example of how one may use our results to obtain justified theoretical predictions. Namely, solved simple teacher-student model with power law correlated data, and showed that the training dynamics and convergence both depend on the spectral density of the Gram matrices studied in the main text. We obtained analytical expressions for the training and generalization losses.
We stress that on their own, these findings do not attempt to fully explain many aspects of neural network dynamics and generalization, which depend on additional factors beyond the bulk spectrum, such as the large outlier eigenvalues, eigenvectors and higher moments.
Analyzing the interaction between these elements and learning dynamics/genneralization remains an important open question, as recent works have started to demonstrate how outliers impact early gradient steps and network collapse.
For instance, as shown by (Seddik et al. 2020) and verified by our results, the outliers can also be described by a Gaussian model, but simply not the CGD that we presented in this work, as the largest eigenvalues are expected to describe the most shared features the entire data, and do not demonstrate the local correlation structure of the bulk. They are certainly important in classification tasks, and in particular their effect, as well the effect of the different class mean values are the most important for linear classifiers, as shown in (axe et al. 2019).
Our approach focused more on the regime where one would like to understand improved performance using more and more data, where the largest cigenvalues have long been well captured, and the only performance gain that can be achieved is squeezed out of the bulk alone. This has proven a sufficient path to construct solvable models which approximate real-world generalization curves (Kaplan et al. 2020 Maloney et al. 2022: Mei and Montanari, 2022).
24


Reviews:
[{"review_id": "hzEdnlKzwU", "reviewer": "Reviewer_LgA8", "summary": "The paper investigates the universal statistical structure underlying the covariance matrices (Gram matrices) of natural datasets. Leveraging Random Matrix Theory (RMT), it demonstrates that real-world datasets and correlated Gaussian datasets (CGDs) share universal spectral properties, including a characteristic power-law decay of eigenvalues and chaotic statistical behavior consistent with Gaussian Orthogonal Ensemble (GOE) universality. This universal behavior suggests that natural data covariance structures can be effectively modeled by simpler correlated Wishart matrices, enabling analytical tractability and better understanding of neural network learning dynamics.", "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": "Dear Reviewer LgA8,  \nThank you for your positive reading of our manuscript, deeming our paper acceptable. \n\nBelow, we address your comments/questions:\n\n**Essential References Not Discussed:**\n\nWe believe the list offered is likely a mistake, since all of these references are explicitly cited in our related work section under “Random Matrix Theory”. Please let us know if we have missed something.\n\n**Questions**\n\n1)\t*How robust is the observed universality across modalities other than vision (e.g., natural language or audio)?* - The results appear to be quite universal, also across different modalities. In particular, we studied the RMT properties of language data, and found similar results. These results will be postponed to future work.\n2)\t*Could different neural architectures affect the emergence of universal spectral structures?* - Since the essence of our work deals directly with the data itself, the notion of a neural architecture doesn’t explicitly enter into our analysis. However, it is clear (and known) that performing different nonlinear transformations on data, can dramatically change its RMT behavior, since it essentially changes the moments of the underlying data distribution, leading to different sample complexities and training dynamics.\n3)\t*What are the practical implications or potential limitations when applying these universal properties to realistic neural network training scenarios?* - The practical implications could range from different sampling techniques based on PCA (and eigenvector sampling) and obtaining scaling laws for sampling required to reach ergodicity, to a more detailed study of the structure of data required for successful learning, moment by moment (second, third etc.).\n\nWe hope that our replies are sufficient to raise your confidence in our work and accept the paper. We welcome any further questions and comments.", "claims_and_evidence": "Extensive numerical experiments across multiple datasets, rigorous statistical tests, and clear comparisons between empirical data and theoretical models from RMT support the claims.", "methods_and_evaluation_criteria": "NA", "theoretical_claims": "Yes, the theoretical claims appear sound. The authors effectively apply well-established RMT concepts and derive theoretical predictions, which closely match the empirical data.", "experimental_designs_or_analyses": "Experiments are thorough, involving various real-world datasets (MNIST, FMNIST, CIFAR-10, ImageNet), uncorrelated Gaussian data, and correlated Gaussian datasets. The methodology is clearly explained and robustly implemented. However, experiments are limited to image domain only and it would be interesting to see how the findings extend beyond image modality.", "relation_to_broader_scientific_literature": "The paper effectively situates its findings within the broader context of neural scaling laws, chaos theory, and Random Matrix Theory. It provides a novel connection between empirical observations of natural datasets and theoretical models traditionally used in quantum chaos and statistical physics.", "essential_references_not_discussed": "Key related works essential to understanding the paper include studies on neural scaling laws (Kaplan et al., 2020; Maloney et al., 2022), foundational papers in RMT (Mehta, 1991; Couillet and Liao, 2022), and previous analyses of data spectral structure (Ruderman, 1997; Pennington & Worah, 2017).", "other_strengths_and_weaknesses": "Clear identification of a universal statistical structure in real-world datasets.\n\nStrong theoretical underpinning via RMT.\n\nRigorous and comprehensive experiments.", "other_comments_or_suggestions": "NA", "questions_for_authors": "How robust is the observed universality across modalities other than vision (e.g., natural language or audio)?\n\nCould different neural architectures affect the emergence of universal spectral structures?\n\nWhat are the practical implications or potential limitations when applying these universal properties to realistic neural network training scenarios?", "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": "Oi0DBp0lhv", "reviewer": "Reviewer_FARz", "summary": "This paper aims to exploit Gaussian universality to generate synthetic data that accurately captures the statistics of natural data distributions. In particular, the paper proposes a synthetic model that reproduces the eigenvalue statistics of the covariance matrices of natural data. Since the synthetic data is described by RMT, the authors focus on demonstrating that the RMT predictions match natural datasets with increasing accuracy as the sample covariance approaches the population covariance.", "questions": "", "limitations": "", "rating": 2, "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": "Dear Reviewer FARz,  \n\nWe appreciate your thoughtful reading of our manuscript. Your main concern is the statement that natural data matrices behaving as Gaussian at sufficiently large dimensions is already known in the literature. We argue that this is **not the case** and we will attempt to explain our estimation of the distinction between what was known and our work below.\n\n**Weakness and Questions**\n\nWe believe there is a subtle misunderstanding between us and the reviewer, which may come from the focus of our work. The large body of literature which concerns Gaussian universality concerns the **learning curves of neural networks when taking gaussian data as a surrogate for real data**, or the **equivalence between neural networks in certain limits and Linear noisy networks** applied to **Gaussian data** [1]. None of the works cited by the reviewer deal directly with the statistical structure of the **data itself**. This is a key difference in our opinion. We agree that the studies of generalized linear regression rely on the properties of the covariance matrix which we study, but they are an indirect measurement. Our work then serves two purposes: (1) it is an agnostic, and much stronger validation of the conjectures used in the literature to model the data (2) the r statistics (as well as the other local statistical measures) can offer a systematic path towards explaining the discrepancies between predictions for gaussian and real world data, since $r^{(n)}$ statistics explicitly measure deviations from Gaussianity (3) the local statistics inform you of the number of samples required to reach the RMT regime, showing a non-trivial universal behavior which is not discussed in the previous works (4) the fact that we find level repulsion could have implications for the localization of eigenvectors, which certainly does have an effect for learning tasks, particularly classification tasks. Note that some of the works you cited are already discussed in our text, particularly in App. F, but we should highlight these in the main text.   \n\n**Regarding specific questions:**\n\n1.\t*”How does your analysis and measurements extend the empirically known Gaussian universality principle?*” – as we replied in the above paragraph, please see our points (1,2) and (3).\n2.\t*”Do the specific measurements you proposed shine light on any specific training phenomena? E.g., learning problems like linear regression or neural networks?”*  - We have not studied the direct effects of our results on models beyond App. F, but it is clear that this is the key question that will be studied in future work. The intuition is that by studying the higher level spacing statistics on real data we can systematically describe the deviations from Gaussianity. By focusing on the effects of these deviations for learning curves, we should be able to make better predictions for learning curves as described in (2) above. Additionally, see point (4) in the previous paragraph. \n\n**Methods**\n\nThe reason for this definition of entropy is a common practice in quantum chaotic systems – we wish to interpret the data matrix as the equivalent of a physical hamiltonian whose complexity can be analyzed using its spectrum. Then, the Renyi entropy definition corresponds to a version of entanglement entropy, measuring the structural complexity of the data.\n\n**Theoretical claims**\n\nWe use the “correlated” term simply due to the Toeplitz matrix assumption. The reason for this model is that it gives a simple “real-space” interpretation of the known spectrum. In that sense the source for anisotropy comes from correlations, offering some direction towards understanding why the power law appears at all.\n\n**Experimental design**\n\nWe thank the reviewer for bringing this work to our attention. We will incorporate the proposed methods into our work. Still, we would argue that since the errors reported in Simon et al. with direct methods are at most ~5%, our results can be trusted, since we do not focus on the precise estimation of $\\alpha$, but rather showing that it is generically between 0 and 1 for natural data. We believe our observations regarding convergence trends and spacings are also reasonable within this  5% error margin.\n\n**Relation to broader literature**\n\n*”Line 78 \"and many advances have been made by appealing to the RMT framework.\" Give examples?”* - We relegated some of the relevant citations to the appendices due to space limitations, we will revise our text to better reflect the contributions in the main text.\n\n**Comments**\n\nThank you for the comment regarding Eq.1, it will be fixed.\n\nWe hope that our replies are sufficient to re-evaluate the merit of our work and change your score of our paper. We welcome any further questions and comments", "claims_and_evidence": "Claim 1: Powerlaw spectrum of population covariance matrices. This claim is well supported by the experiments and is in agreement with much previous work.\n\nClaim 2: The powerlaw spectrum originates from underlying correlational structure. This claim may be true but it is not supported by the given evidence. See Theoretical Claims section.\n\nClaim 3: The correlated Gaussian datasets are a correct proxy for natural data. This claim is well supported by the experiments, but it follows from the Gaussian universality principle, which has been extensively studied in previous work. See Essential References Not Discussed section.\n\nClaim 4: Convergence of sample covariance to population covariance corresponds to convergence of sample statistics to RMT statistics. I believe this claim already follows from Gaussian universality and law of large numbers.\n\nClaim 5: Shannon entropy is correlated with the local RMT structure, is smaller in correlated datasets, and converges to the population entropy in much fewer samples. The meaning of this claim is unclear. The proposed definition of Shannon entropy does not appear to match the commonly accepted definition for Gaussian distributions. See Methods And Evaluation Criteria section.\n\nTo me, the main claim of this paper seems to be that Gaussian universality holds. This has been confirmed in many previous works, so it's not clear what the novel contribution is.", "methods_and_evaluation_criteria": "The definition of dataset entropy is confusing to me. I don't see the connection between the entropy as defined and the entropy of the data distribution, which seems the more natural measure of entropy. What motivates the proposed definition?", "theoretical_claims": "I don't see novel proofs or derivations in this paper. There are some conceptual issues I will point out here.\n\nIt's not clear to me why the setup is called \"correlated.\" The data vectors are drawn from a Gaussian with diagonal covariance. It seems that the uncorrelated/correlated distinction being drawn is in fact a distinction between isotropic/anisotropic. Relatedly, it's not clear to me what the Toeplitz matrix brings to the analysis. The exponent in the Toeplitz matrix must be measured from the dataset to be modeled anyways; why not simply use that measurement to directly set the diagonal of the covariance according to eq. 3? With this simplification, I believe the main claim boils down to the usual Gaussian universality assumption.", "experimental_designs_or_analyses": "It is known to be notoriously difficult to directly measure powerlaw exponents (see https://arxiv.org/pdf/cond-mat/0402322). Appendices B and C in https://arxiv.org/pdf/2311.14646  also point out that it is easy to make drastic measurement errors trying to directly fit the sample covariance eigenvalues. Due to the ease of measurement errors, I do not know whether to trust the measurements in this paper. The previously mentioned paper provides a more reliable procedure, based in RMT, for exponent estimation.", "relation_to_broader_scientific_literature": "Line 78 \"and many advances have been made by appealing to the RMT framework.\" Give examples?\n\nSee Other Strengths And Weaknesses section.", "essential_references_not_discussed": "The following papers give convincing evidence for Gaussian universality in natural datasets. Many of these concern kernel ridge regression, where the random matrix in question is the kernel matrix \\phi(X)^T \\phi(X) rather than the covariance XX^T, but the result is analogous (since the spectra of covariance and gram matrices are the same). \n\nhttps://arxiv.org/pdf/2302.08923\nhttps://arxiv.org/pdf/2102.08127\nhttps://arxiv.org/pdf/2203.06176\nhttps://arxiv.org/pdf/1905.10843\nhttps://arxiv.org/pdf/2311.14646", "other_strengths_and_weaknesses": "My primary concern with this paper is that the main claim seems to be that natural data matrices behave like Gaussian random matrix ensembles when the data dimension is large enough. This is already a well-established claim; there is a broad literature on Gaussian universality in high-dimensional statistical inference problems (see Essential References Not Discussed). This paper suggests some physics-inspired diagnostics for further verifying this claim and find agreement; this is useful but, on its own, I don't think it is a sufficient contribution for this venue.\n\nThe strengths include proposing to port over some empirical measures of RMT behavior from systems exhibiting quantum chaos. Another novel idea was to suggest a data generation process via Toeplitz matrices, which produces the anisotropic Gaussian features, but I didn't find it convincing.", "other_comments_or_suggestions": "Eq 1: I believe it should be X_{ja} or X^T_{aj} instead of X_{aj}", "questions_for_authors": "1. How does your analysis and measurements extend the empirically known Gaussian universality principle?\n2. Do the specific measurements you proposed shine light on any specific training phenomena? E.g., learning problems like linear regression or neural networks?", "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": "JKdPrZJEZh", "reviewer": "Reviewer_kWZr", "summary": "The paper examines the empirical eigenvalues of Gram matrices of natural image datasets, showing that they follow a power law. It offers a simple model that reproduces this power law behavior using Gaussian data with long-range correlations. These results suggest that natural image Gram matrices can be approximated by Wishart random matrices with a simple covariance structure, enabling rigorous analysis of neural network behavior.", "questions": "", "limitations": "", "rating": 4, "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": "Dear Reviewer kWZr,\n\nThank you for your careful review and positive appraisal of our work, deeming it acceptable at ICML. We’re glad you found the work interesting and well supported by empirical evidence.\nBelow, we address the remaining issues you raised, as well as questions and comments.\n\n**Related work**\n\n1)\tWe apologize for not bringing the connection between neural scaling laws, learning curves and the spectral properties of the data more to the forefront in this section. We provide a concrete example for this connection in Appendix F, but we will clarify these connections in the related work section as well in the revised version. \n2)\tThank you for pointing out these references, we will update the revised manuscript to include them both.\n\n**Weaknesses**\n\n1)\t*Highlighting the connection to real world DNN training on CGD* - We understand the reviewer’s suggestion and agree it has considerable merit, however, our intent in this paper was to bridge a long line of existing semi-empirical works (including ones that train networks on gaussian versions of real world data, see Ref [1,2] and similar works), and purely theoretical works that make the assumption of gaussian data. We will include Ref [1,2] and several other works that perform such experiments to better connect our results to the already large existing literature.\n\n**Comments**\n\n1.\t*Line 153 right column, above eq 2: 'with' is misplaced.* - Thank you for this comment, the typo will be fixed. \n\n2.\t*I find the index notation confusing, e.g. $X_{ia}\\in \\mathbb{R}^{d×M}$; usually $X_{ia}\\in \\mathbb{R}$ as in the matrix element.* - We accept your comment, and will augment the final version accordingly.\n\n\n**Questions**\n\n1.\t*Lines 434-436: “accomplishing learning tasks for which the spectrum is insufficient and eigenvector properties are needed, requires a different scale of data samples.” - would it not be fair to say that for most typical learning tasks eigenvector properties are needed?*  - We agree that there are many cases, such as classification tasks, the eigenvectors are crucial, and we study to what extent in a future work. However, for teacher-student settings, many other regression based settings, the eigenvectors are not necessarily important, and the spectrum is sufficient to predict performance and dynamics. \n\n2) & 3.  *What information is NOT captured by the CGD? why does the scaling law behavior begins around i=10?* - First, the CGD does \n             not capture the top ~10 eigenvalues, since these eigenvalues contain the least “universal” properties of the data, and constitute the most particular aspects of every different dataset, also corresponding to the most informative eigenvectors. Therefore, the CGD, which describes the universal part of the spectrum, is shared between datasets and is only valid for the bulk of eigenvalues. Second, as stated in the main text, the CGD does not capture the eigenvector structure, as we do not assume a special basis for the CGD model. It would be interesting to explore which combinations of bases and spectra would render the CGD model a more complete description of particular datasets.\n\n4.\t*Lines 189-191: “For real-world data, we consistently find that $\\alpha>0$” - is there a simple reason for this?* - This is a very good question, which is currently still open. One possibility is that these values of power law spectra can emerge from hierarchical data structures with varying levels of rare features [3]. Answering this question could also shed light on certain aspects of neuroscience, see [4] for the emergence of power laws in human visual cortex activity data. We would also like to note that $\\alpha>0$ seems to be prevalent for natural datasets, but not necessarily for other types of data. For instance, simulations of physical systems can give rise to negative exponents (see for instance [5] for the case of turbulence).\n\nReferences:\n\n[1]  Refinetti et al., Neural networks trained with SGD learn distributions of increasing complexity (2023)\n\n[2] Székely et al., Learning from higher-order correlations, efficiently: hypothesis tests, random features, and neural networks (2024)\n\n[3] Cagnetta et al., How Deep Neural Networks Learn Compositional Data: The Random Hierarchy Model (2024)\n\n[4] Gauthaman et al., Universal scale-free representations in human visual cortex (2024)\n\n[5] Levi et al., The Universal Statistical Structure and Scaling Laws of Chaos and Turbulence (2023)", "claims_and_evidence": "As far as I can see, all claims are clearly stated and well-supported by theoretical or empirical evidence.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria make sense for the problem or application at hand.", "theoretical_claims": "I checked the theoretical claims and found no major issues with their correctness.", "experimental_designs_or_analyses": "I found the experiments to be generally sound and valid.", "relation_to_broader_scientific_literature": "I found the subsection “Neural scaling laws” under Background and Related Work somewhat non-satisfactory. It leaves an impression that the relation between the spectral properties of the data and the learning curve is unclear.", "essential_references_not_discussed": "• Ba et al. 2022 - High-dimensional asymptotics of feature learning: how one gradient step improves the representation. \n\n• Wang et al. 2023 - Spectral evolution and invariance in linear-width neural networks.", "other_strengths_and_weaknesses": "Strengths: \n\n1. The paper is well structured, clearly written, and rather easy to follow, given the heavy math. \n\n2. Tackles a very timely subject, both interesting and relevant. \n\nWeaknesses: \n\n1. I think that the main weakness is establishing the claim in the last sentence of the abstract: “image Gram matrices can be approximated by Wishart random matrices with simple covariance structure, enabling rigorous analysis of neural network behavior.” The paper could be much improved if this point was made more central, beyond what is provided in App. F, e.g by some experiments of training a DNN on real-world data vs on the corresponding CGD.", "other_comments_or_suggestions": "Typos and other small issues: \n\n1. Line 153 right column, above eq 2: 'with' is misplaced. \n\n2. I find the index notation confusing, e.g. $X_{ia}\\in\\mathbb{R}^{d\\times M}$; usually $X_{ia}\\in\\mathbb{R}$ as in the matrix element.", "questions_for_authors": "1. Lines 434-436: “accomplishing learning tasks for which the spectrum is insufficient and eigenvector properties are needed, requires a different scale of data samples.” - would it not be fair to say that for most typical learning tasks eigenvector properties are needed? \n\n2. What information is NOT captured by the CGD? \n\n3. Can you speculate why the scaling law behaior begins around i=10\n\n4. Lines 189-191: “For real-world data, we consistently find that \\alpha>0” - is there a simple reason for 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": "EJ0eaCzve7", "reviewer": "Reviewer_ZEdS", "summary": "A model of natural dataset (at least natural images) is build based on a random matrix analysis and extensive data analysis. An interpretation is proposed for the power law spectrum of the covariance matrix of input features which generically shows up on natural data. In addition a further characterization of the statistical ensemble is provided in terms of a random matrix ensemble (Gaussian orthogonal ensemble (GOE)) based on a emprirical study of the statistics of level spacing.", "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": "Dear Reviewer ZEdS,\n\nThank you for your careful review and positive appraisal of our work, tending towards acceptance. We’re glad you found the work interesting and well supported by empirical evidence.\nBelow, we address the remaining issues you raised, as well as questions and comments.\n\n**Methods and Evaluation**\n\n*” I don't see how this can be used in practice, since the exponent which is in principle unknown, is estimated itself at finite M. This looks more like an qualitative heuristic estimation.”* - We agree that defining the asymptotic exponent itself at finite $M$ is suboptimal, and therefore its value may not be entirely trustable, but this is truly the best one can do on any real-world finite dataset. The important part is the universal convergence towards this exponent, and the interesting fact that this convergence occurs in two sharply separated phases, as seen in Fig 5 (middle).\n\n**Weaknesses and Questions**\n- *”…level-spacing statistics of the Gram matrices seems not to play (to the best of my knowledge) any role for theoretical analysis of ML.”* - This is inaccurate, please see Appendix F for more details. Briefly, the role of the level spacing metric is to confirm the spectral density distribution, which is critical for analysis of even linear regression in the empirical limit (the resolvent is required to find the optimal network solution, which depends directly on the data covariance matrix’s eigenvalue distribution).\n- *”…whether a 2d model would not be more appropriate to interpret the data?”* - In order to capture the eigenvalue statistics a 1-d model is sufficient, and any higher dimensional model can work, provided that there is a component that decays as a power law. For instance if we were to take a fourier type model, as long as the spectrum had a component that decayed with the magnitude $|\\vec{k}|$ as a power law, the same type of results would hold. Our goal was simply to provide the simplest possible model.  \n-\t“*equation (9) to check, should it not be continuous at tau=1?*” Yes it should, thank you for bringing this to our attention, we will correct the formula in the revised version.\n\nWe hope that our replies are sufficient to raise your confidence in our work and possibly raise your scoring of our paper. We welcome any further questions and comments.", "claims_and_evidence": "The main claims is that the bulk spectrum statistics of the feature-feature covariance matrix  arises from long range correlations following a power law and is in GOE universality class. This claim is backed on a precise data analysis comparing the RMT model obtained with population matrices given by a full band Toeplitz matrix corresponding to a 1d model of homogeneous interactions growing with distance raised to some exponent.", "methods_and_evaluation_criteria": "On one side the RMT model for the considered 1d model of long range interactions is solved in order to get the corresponding spectral density in the proportional scaling regime (fixed ratio N_data/N_features ) using Marchenko Pastur  equations. It is then used to fit very well the spectral density of a large selection of datasets (image datasets). Then the level-spacing statistics corresponding to these datasets is determined empirically using methods from quantum chaos, namely the unfolding method used to normalize the level spacing w.r.t the spectral density, the distribution of level-spacing statistical indicator called the r-statistics and the spectral two-point  form factor,  \nwhich have well defined signature for different RM ensembles. These are then compared with the expected result for GOE yielding a very good agreement. An estimation of the number of samples needed to be in this RMT regime is also provided based on various indicators having well defined values in the matrix ensemble once the exponent is known. I don't see how this can be used in practice, since the exponent which is in principle unknown, is estimated itself at finite M. This looks more like an qualitative heuristic  estimation.", "theoretical_claims": "No specific theoretical claims", "experimental_designs_or_analyses": "As discussed already in  the Methods and evaluation criteria section, the experimental data analysis is perfectly sound.", "relation_to_broader_scientific_literature": "I think the bibliography is fine.", "essential_references_not_discussed": "No essential reference is missing in my opinion.", "other_strengths_and_weaknesses": "- strength: the statements on the statistical properties of natural data are interesting and instructive, I found the paper clear and the methodology rather convincing.\n- weakness: the scope is limited to natural image and their texture properties. This might have limited impact in ML theory (see my questions below).", "other_comments_or_suggestions": "equation (9) to check, should it not be continuous at tau=1?", "questions_for_authors": "I find the paper quite interesting but I wonder about the real scope of this work regarding ML in particular the modelling of data for theoretical analysis. Among the claims, the fact that the real datasets fall into the GOE is well expected (Wishart matrices obtained from large data embedded in large dimensional spaces), but the empirical analysis of the level-spacing has the merit of providing a clean validation of this hypothesis. Still, the impact on the theory of ML is not clear to me, since the power-law hypothesis (with possibly some additional spike model of outliers) of the population matrix is already commonly used, while the level-spacing statistics of the Gram matrices  seems not to play (to the best of my knowledge) any role for theoretical analysis of ML. Then remains the 1 dimensional long-range correlated model of the population matrix, in the form of Toeplitz matrices, which looks quite artificial to me. I guess any d-dimensional  model with spatial power-law decay would lead to the same result, but the the 2-point function would decay with distance for d>1 instead of increasing in 1d. I believe this bulk modelling  is a simplified model of texture in this context of image dataset and wonder whether a 2d model would not be more appropriate to interpret the data?", "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: This paper investigates the universal statistical structure underlying the covariance/Gram matrices of natural datasets, and compares empirical findings against theoretical predictions derived from statistical physics and random matrix theory.

The authors have done a commendable job during the rebuttal phase, and all reviewers are now generally **convinced of the significance and merit of the work**. 
It is important to note that this paper does not aim to fully resolve the question of universal structures in natural data. 
Rather, it represents a first step toward a deeper understanding of this fundamental and challenging problem.

As highlighted by Reviewer FARz and others, there remains ample room for improvement, including providing novel theoretical analysis and further clarifying the central ideas and contributions. 
Nonetheless, given the importance of the problem and the value of this initial contribution, I believe the paper is a meaningful addition to the ICML community.

As a result, I recommend acceptance of this paper to ICML 2025.
Please ensure that the additional numerical results and clarifications provided during the rebuttal phase are incorporated into the final version of the paper.
Decision: Accept (poster)