A Theoretical Framework For Overfitting In Energy-based Modeling
Authors: Giovanni Catania
Aurélien Decelle
Cyril Furtlehner
Beatriz Seoane
Venue: ICML 2025 poster (2025)
Keywords: energy-based models
overfitting
random-matrix theory
inverse problems
early-stopping
training dynamics
Bolzmann Machine
Primary area: general_machine_learning->unsupervised_and_semisupervised_learning
TL;DR: We analyze the impact of limited data in the training of energy-based models, focusing on eigenmodes' dynamics, giving a theoretical perspective of early-stopping and data-correction protocols to improve the quality of the inferred model

Abstract: We investigate the impact of limited data on training pairwise energy-based models for inverse problems aimed at identifying interaction networks. Utilizing the Gaussian model as testbed, we dissect training trajectories across the eigenbasis of the coupling matrix, exploiting the independent evolution of eigenmodes and revealing that the learning timescales are tied to the spectral decomposition of the empirical covariance matrix. We see that optimal points for early stopping arise from the interplay between these timescales and the initial conditions of training. Moreover, we show that finite data corrections can be accurately modeled through asymptotic random matrix theory calculations and provide the counterpart of generalized cross-validation in the energy based model context. Our analytical framework extends to binary-variable maximum-entropy pairwise models with minimal variations.
These findings offer strategies to control overfitting in discrete-variable models through empirical shrinkage corrections, improving the management of overfitting in energy-based generative models.
Finally, we propose a generalization to arbitrary energy-based models by deriving the neural tangent kernel dynamics of the score function under the score-matching algorithm.

A Theoretical Framework For Overfitting In Energy- Modeling
Giovanni Catania' Aurélien Decelle 12
Cyril Furtlehner Beatriz Seoane'
Abstract
architecture and optimization, employing various implicit regularization mechanisms (Belkin, 2021). Finding the op- timal balance between maximizing data utility, preserving generalization, and ensuring the privacy of training data rep- resents a critical trade-off that can be challenging to pinpoint. In supervised learning tasks like classification, overfitting is readily identified using standard practices. Metries like test-set accuracy, particularly when augmented by cross- validation in data-scarce scenarios, clearly signal overfitting, enabling strategies like early stopping, regularization, and hyperparameter tuning to mitigate it. Furthermore, training and generalization performance in regression and classifi- cation tasks are now well understood certain simplified regimes, such as high-dimensional ridge (Atanasov et al., 2024; Advani et al., 2020; Saxe et al., 2014; Tomasini et al., 2022) logistic (Mai al., 2019; Loffredo et al., 2024) regression or numerous more complex setting of non-linear regression in various scaling regime (see for instance (Mei et al., 2018; Arnaboldi et al., 2023; Saad & Solla, 1995) among many other recent works). This gives the possi- bility to assess some simple indicator like the generalized cross-validation (GCV) (Golub al., 1979), an exact rela- tion between train/ /est errors valid for the ridge regression that can be derived using leave-one-out argument (see e.g. (Furtlehner, 2023)). This methodology is also relevant in deep learning contexts (Wei et al., 2022), particularly in over wver-parametriieed regimes where it aligns with observed stochastic gradient descent behaviors (Patil et al., 2024).
We investigate the impact of limited data on train- ing pairwise energy-based models for inverse problems aimed at identifying interaction net- works. Utilizing the Gaussian model as testbed, we dissect training trajectories across the eigen- basis of the coupling matrix, exploiting the in- dependent evolution of eigenmodes and reveal- ing that the learning timescales are tied to the spectral decomposition of the empirical covari- ance matrix. We see that optimal points for early stopping arise from the interplay between these timescales and the initial conditions of training. Moreover, show that finite data corrections can be accurately modeled through asymptotic random matrix theory calculations and provide the counterpart of generalized cross-validation in the energy based model context. Our analytical framework extends to binary-variable maximum- entropy pairwise models with minimal variations. These findings offer strategies to control overfit- ting in discrete-variable models through empirical shrinkage corrections, improving the management of overfitting in energy-based generative models. Finally, we propose a generalization to arbitrary energy-based models by deriving the neural tan- gent kernel dynamics of the score function under the score-matchin algorithm.
Recent advancements in training and architecture have greatly enhanced the generative capabilities of neural net- work models across various fields (Bengesi et al., 2024), en- abling the creation of photorealistic images, credible speech synthesis, and biologically functional synthetic proteins (Wu et al., 2021). Despite this progress, selecting optimal mod- els from a pool remains challenging due to noisy training data often leading to undetected overfitting. Yet, detecting overfitting in unsupervised learning settings, particularly for generative modeling, is elusive but crucial, especially with sensitive datasets like human genomic data (Yelmen et al., 2021; 2023) and copyrighted content. Unlike supervised learning, unsupervised learning lacks clear overfitting in- dicators, complicating model development and validation. While some theoretical insights on optimal regularization tuning for simple energy-based models exist (Fanthomme
1. Introduction
Controlling overfitting is basic in machine learning, par- ticularly as modern, wneetparameterized architectures en- hance learning capabilities. To prevent learning noise or irrelevant patterns, numerous empirical solutions have been proposed that modulate the model's implicit bias through its
'Departamento de Física Teórica, Universidad Complutense de Madrid, Spain. "Escuela Técnica Superior de Ingenieros Indus- triales. Universidad Politécnica de Madrid, Spain "Inria-Saclay, Université Paris- SSclay, LISN, Gif Yvette, France.. Correspon- dence to: Giovanni Catania <gcatania ucm.es>.
Proceedings of the nrnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

A thcoretical framework for overfitting in energy-based modeling
et al., 2022), practical indicators such as early stopping points during training dynamies remain undefined. More- over, estimating log-likelihood for model selection poses significant computational challenges (Béreux et al., 2022: 2024). Consequently, there is an urgent need for methods to detect and mitigate overfitting in these contexts.
2. Gaussian Model
The Gaussian Energy-Based Model (GEBM) specifies a multivariate Gaussian distribution for real-valued variables æ E RN, characterized by 2-body interactions encoded within symmetric, positive-definite coupling matrix, J e RNNN, The GEBM is the simplest model that effectively captures the first and second-order statistics of set of data. For the purposes of this analysis, we assume O means for the data components, thus simplifying the initial model by excluding the learning of external biases. Nevertheless, the theoretical framework presented below can be readily extended to the above to accommodate non-zero means. The probability distribution of configuration V is then:
This paper focuses on energy-based models (EBMs) (Ack- ley et al., 1985). which encode the empirical distributions of various data types -scch neural recordings (Roudi et al., 2009), images (Du & Mordatch, 2019), and ge- nomic (Yelmen et al., 2021) or proteomic (Morcos et al., 201 sequences- inte a probability framework rooted in Boltzmann's 's law. By adopting a Bayesian approach, EBMs aim to maximize the likelihood function, enabling the gen- eration of new data that closely resembles the training set and facilitates the extraction of detailed microscopic in- sights. EBMs range from simple Boltzmann Machines (BMs) and Restricted Boltzmann Machines (RBMs) to more complex architectures like convolutional neural networks, making them versatile in statistical physics for solving in- verse problems like deducing Hamiltonian parameters from observed data. The interpretability of simple EBMs en- ables uncover underlying rules within datasets: this ca- pability has proven highly effective in fields ranging from neuroscience (Roudi et al., 2009) to bio- bio-molecular struc- ture prediction (Cocco et al., 2018) predominantly through pairwise maximum maximum-entropy models. Recent advancements extend these applications, using complex EBMs to infer high-order interactions (Decelle et al., 2024; 2025; Feinauer et al., 2022; Feinauer & Lucibello, 2022) constitutive patterns (Tubiana et al., 2019; Decelle et al., 2023), signifi- cantly deepening our comprehension of data structures.
p (x J) = (2x) N/2 deet Je 1 Jm
Jo
(1)
It is straightforward to check that the population covariance
matrix of such distribution is C = E, [x> = J with
Eअ [1] denoting the average with respect to (1).
Inference problem. Consider dataset D = {x*) J.M
J=1 with
M entries generated with GEBM model with coupling matrix J Our objective is then to find the parameters J that best approximate the empirical distribution of the data formally (a)) = M- M- == 5(x x²/) with the probabilistic model (1). Without prior information about the model parameters J, the maximum likelihood (ML)
estimator jML.M is calculated as jMLLM = (Car)
where ÕM is the empirical covariance matrix from M data points, provided it is invertible (MacKay, 2003). Denoting with N the number of data components (data dimensions), this condition requires that M 2 N, assuming samples to be independent. Clcarly, when M- + co, JML.,M recovers the true set of parameters used to generated the data, J
This work develops a theoretical framework for understand- ing and mitigating overfitting in EBMs. We begin with a simple Gaussian model as a fundamental non-trivial exam- ple, using to quantitatively analyze overfitting through synthetic experiments with predefined ground truths. We examine eigenvalue dynamics using artificial covariance ma- trices that simulate real datasets, exploring how overfitting arises from different learning timescales associated with various eigenmodes of the empirical covariance matrix. We address inaccuracies in learned eigenvalues with corrections based on random matrix theory (RMT), showing that the quality of model generation in EBMs is less affected by the lower modes of the covariance matrix, while the accu- racy of inferred couplings is significantly impacted. We demonstrate that regularization techniques like shrinkage corrections are crucial to counteract overfitting, providing a robust framework to refine EBM training by considering umite-sample-siee effects. This approach also informs our analysis of more complex models like the BM, underscoring the importance of regularization strategics to enhance model reliability and predictive accuracy.
Training dynamics. The GEBM stands out as one of the few high-dimensional inference problems where an analyt- ical expression for the ML estimator is available, indepen- dent of both M and N. However, our focus here is on the training dynamics associated with an iterative maximization of the likelihood function through gradient ascent dynamics, as is typical in EBMs. This approach allows us to explore the adaptive process of parameters' estimation over time.
In the GEBM, the log-likelihood (LL) of the param- eters J depends only C. and it reads L(J) - En, Jiy CM JiyC + 1 log det J. This quantifies how well J matches the observed data. In a standard gradient ascent algorithm, the update rule for the parameters reads:
BC J J' + 7
where 7 is the learning rate. Assuming non-symmetric per- turbations on the parameters Jej- the gradient in (2) reads:
ac Cyyy + +Es [m,xa]= - Cyy + {J +J
(3)
JJ.y
2

---

A thcoretical framework for overfitting in energy-based modeling
The second equality comes from the exact expression of 2-point correlations of the Gaussian model in terms of its coupling matrix J: this is equivalent to assume that we perfectly sample the model with an infinite amount of con- figurations at any I. For a more generic EBM, another source of noise should be added due to the finite number of samples (or chains) used to estimate the empirical correlations used to compute the gradient.
related the concept of spectral bias (Rahaman et al., 2019) although here the decomposition is spectral rather than in a Fourier basis and it has been characterized theoreti- cally in the case of linear regression (Advani et al., 2020), The interaction between these varying timescales can result in an initial phase where the strongest components of the dataset's PCA are effectively captured, followed by a phase where training begins to adjust noise-dominated directions, potentially leading to overfitting.
To analyze the learning dynamics, use the spectral de- composition of J and project the gradient onto its eigen- basis. denoted by V = [va)a- From (3), the gradient projected on modes a and /3 is expressed as:
We' We've shown that GEBMs' learning dynamics are governed by the spectral decomposition of C.M, with sinite-sample effects arising from changes in the spectrum due to finite M. This falls within the realm of random matrix theory (RMT) (Potters & Bouchaud, 2020), as we detail shortly.
ac 0J
2M +
6us J.
(4)
0.
where cas is the projection of C.M. Generally, for any matrix M, we write map de
Asymptotic RMT analysis. In our simplified setting of GEBMs. the parameters of the model J (t) along the learn- ing trajectory are an explicit function of the empirical covariance matrix C.M upon choosing the same constant B, = B > -/cc in (6) for the initialization (B = -1/c cor- responds to the initialization J(] = 0) and assuming that (t) is aligned with C.M at t = 0. This choice simplifies considerably the analysis, and as explained in Appendix G. Using RMT, all relevant quantities can be derived in closed forms based solely on the population spectrum v and the aspect ratio p = M /N, under the asymptotic proportional scaling where M, NN + co with p held constant.
v2 Mus.
This approach enables us formulate set of evolution equations for the eigenvalues {Ja} and the rotation of the cigenvectors. By assuming an infinitesimal learning rate, we can transform the timerete-time update equation (2) into a continuous set of differential equations (see Appendix A for further derivation details):
dJa dt Ja
do
e.
MM
a
for a 8, (5)
dt
Ja
where F is a timescale set by the learning rate, T = 1/7. From Eq. (5) we see that eigenvectors of J stop rotating when they align with the eigenvectors of C.M to cancel out the numerator aM Eq. (5) can be integrated analytically, and MM can be replaced by the matrix eigenvalue 2Mr ac
We are interested in the train and test energies:
Eurain = - and Erest = N -T[CC,
the coupling error E def N- JI with || 11?. the Frobenius norm, and the LL (train and test)
The solution of (5) can be expressed in an explicit (although not closed) form using Lambert Wo function, namely:
L. def 2N log det
2 Ptrain,test
where C def lim CM is the population matrix. In addition, we will also explore the behavior of the maxi- mizer of the log-likelihood with regularization, i.e. C[J] = I I.rrann [J] - AJ)) focusing on A(J) = Tr[J2] for L2 ridge regularization, and A(J) Tr[J] for L, lasso reg- ularization on the spectrum. T. is applicable since J is symmetric and remains positive definite throughout the tra- jectory, ensured by the logarithmic barrier.
GM Wo (e0)
Ja (t) G.M +
(6)
with the constant B3 is fixed by the initial condition at t t=0 O. Eq. (6) delineates the evolution of each eigenvalue, which progresses independently once the eigenvectors of J align with those of CM, A crucial aspect of this equa- tion is that the relaxation time it takes for an eigenvalue to reach its steady-stan state value )((²) = =lim, Ja (t)==//a) is inversely proportional to the square of the corresponding eigenmode in the covariance matrix: indeed, for + oo Eq. (6) describes an exponential relaxation to the fixed point with timescale x (
We simply quote here the result of the asymptotic lim- its (for p >>> 1, more details in Appendix G), based on RMT (Marčenko & Pastur, 1967; Ledoit & Péché. 2011). First, the spectral density D of AM reads in this limit:
This relationship shows that the evolution of each eigen- value is closely linked to the significance of the correspond- ing cigenvector in representing the data, so that stronger modes in the covariance matrix are learned more quickly than weaker ones. The idea that information is learnt pro- gressively starting from strong PCA's directions is closely
0xx)= pA.(z)
((())
(8)
1 (-),(-) +r.
where A(z) = A, (x) + iA and T(=) = r,) (T)))))
for 2 = x + =0 obey the self-consistent equations
A(z) = -
[(2) 1 ((dz)z
--A(a]=
3

---

A thcoretical framework for overfitting in energy-based modeling
(a)
MNIST
(b)
Synthetic
ground truth coupling matrix J' is constructed, through
100
M 109
10
. -02
=²= Using this setup, we generate multivariate Gaussian distribution and extract M data points from it. These data points are then used to train new GEBM us- ing the empirical covariance matrix, CM,, derived from the M of these samples, with the goal of inferring the original model parameters.
M=4-100
2 .M 10
M 1116--00
10
10
10-1
M = 1.2- 10
M 00
10-6
10-2
0.0
0.5
1.0
0.0
0.5 a/N
1.0
a/N N
As previously discussed, the training dynamics of each mode of J are directly linked to the eigenvalues of C.M, To enhance this analysis, we have developed a synthetic model for the spectrum of C*, which influences the spectrum of CM in scenarios with finite datasets. This model closely mimics the eigenvalue spectra of real datasets, as illustrated in Figure 1-(a), which shows the eigenvalue spectrum (in descending order) of covariance matrices from MNIST for several sizes M (more examples
Figure 1. (a): Eigenvalue spectra of the empirical covariance ma- trices for MNIST dataset (Deng, 2012). Black lines show spectra using the full dataset size (M") ). while scatter colored points rep- resent subsets (M < M*) (b): Black line shows synthetic pop- ulation eigenvalue spectrum based (10) for N 100, r r=0.9 0.9, 0.9, 7 = 1.1. Ø1 = 10 22 = 10; colored points show the eigenvalues from CM calculated by sampling different M config- urations from GEBM model with J =²= (Eq. (I)).
given in Appendix B). Our analysis reveals that the spectrum of C.M remains rel- atively stable w.r.t. M for significant number of modes, indicating aa x c = = However, smaller eigenvalues fluctuate markedly with M; they tend to be underestimated as M decreases, suggesting aM < c for small 6%+ and are slightly overestimated for the larger eigenvalues. This behavior is rigorously characterized using RMT tools in simplified data models (Baik& & Silverstein, 2006; Ledoit & Péché, 2011) 1). Additional insights into the conservation of eigenvectors across modes are detailed in Appendix B.
in terms of the population spectral density v(dx). In turn we obtain
Errain = J J0 dyj(t) [A + (w))
Erest = + 10 (w))
(0 1),
while the coupling error takes the form
v(dz) Eg = Jo at + S0 v(dr)j(z) [j(x) - 2 2(1 (1-p)+20A.(o))
Inspired by these findings, we will characterize our synthetic population matrix C* by an eigenvalue spectrum {c² N generated according a mixture of power laws. The cumu- lative distribution is defined as follows:
where j(r) is one of the analytical functions Jr. and jtz corresponding respectively the time dependent, Li and Le regularized forms of J. Remarkably, for the Lı regularized coupling matrix we get a deterministic relation between the train and test energies (Appendix G)
P[A<x] =r D-
(2 + +
1 (1.22)
(10)
where I1 (1.b) denotes the indicator function in the interval (). This setup distinguishes between "strong" modes with c > 1 and "weak" modes with c < 1. with their prevalence controlled by parameter r. Parameters B and y represent the power-law exponents for these two cate- gories, with I1 and z2 the respective lower and upper cut- off'. Model (10) is chosen i) mimic the eigenvalue distri- bution of a realistic dataset's i covariance matrix, though our numerical results are robust to specific spectral details, and ti) to extract asymptotic quantities in the mantinuous-densiig limit N, M - 0o (with p p= M/N finite) using RMT.
Exest (1-p-'Eurate
Eusinn
which is the counterpart of GCV for GEBMs which might be usable in practice for arbitrary EBM (in the same way as GCV can be used for deep regression models), as it allows one to get an estimation of the test LL. This concept remains have focused
topic for future research. Instead, strategies for data cleaning and regularization, specifically employing shrinkage techniques (Bun et al., 2017). Using model of the data defined in the next section, allows us to specify v(x) in order to assess these strategies by comparing with the expected optimal performances given by RMT.
Figure 1-(b) shows an example spectrum of population eigenvalues ct from Eq. (10) with N 100 (black line), alongside empirical estimates of finite-data eigenvalues aM (scatter points) for various M values, demonstrating that strong modes remain stable despite finite-M noise, while weak modes are consistently underestimated. The full ma- trix C' is finally assembled by projecting the diagonal ma- trix of these eigenvalues onto a random orthogonal matrix U {ux ) resulting in C* =ea C u u
3. Modeling realistic data covariances
To effectively study the impact of finite number of data on the learning process of a GEBM in a controlled setting, we need to define a synthetic model that facilitates the analysis of different learning timescales. The first step is to artifi- cially create a population covariance matrix C", from which
4

---

A thcoretical framework for overfitting in energy- -based modeling
analytic
numeric
(c))
gir", (911
10
10"
10
10
10P
10
10-3
CC
(c)
Ja(t)
100
109
10
10
10P
10a
10
1
10
-20
(b)
10
test- I.,
cooso 039S3
10
(a)
(b)
Lapt
10-1
10a
10
10
-40-
10P
44444
E stable
7
7
Ec
10
10
10"
109
10
10F
Figure 2. Training dynamics of the GEBM from a population matrix C" (in (a), with system size and parameters matching those in Fig. 1-(b)), and from an empirical covariance matrix C.M (generated from C through (2). with م == 2.11 (in (b)). (a)-(b) display the analytic evolution of eigenvalues Ja toward the steady- state (lines), and a comparison with numerical training (points, in (a)). In all cases the initial condition is identity matrix.
P
Figure 3. Results for GEBM training with finite data. (a)-(b)-(c) display respectively the reconstruction error &,, the test-LL and the generation error Ec, all plotted vs time, for various sample sizes M (indicated by color gradient from blue to red for increasing p= M/N. Dashed black lines refer to a training from C" (i.e M oo). (d): comparison between time of minimum reconstructon (circles), maximum test LL (diamonds) and time at which the generation error converges its steady-state value. These quantities are also shown in the related panels for better clarity.
4. Training Dynamies on Synthetic Data
The introduced synthetic model enables analysis of training dynamics in two scenarios: i) an ideal setting with an in- finite amount of samples using C* as the data covariance matrix, and ii) a more realistic situation with a finite dataset, represented by the empirical covariance matrix C.M
trained J (t) defined in the previous section. Fig 3-(a) illus- trates the error's 's evolution over training time. Beginning from an identity matrix, at low p = M/N values, the er- ror displays marked non-monotonic behavior, peaking at specifie mmin (p) before stabilizing at the training's fixed point. At higher p. the error decreases monotonically until stabilization, following a trend consistent with the M o scenario (i.e. using C",, in black dashed line). This behavior, also noted in complex EBMs (Decelle et al., 2024; Agoritsas et al., 2023), underscores the GEBM's utility as a simple model yet capturing complex phenomena in EBMs.
Training Dynamics with Infinite Data. For clarity, we be- gin by training our GEBM using the population covariance matrix C Fig. 2-(a) illustrates the evolution of the cou- pling matrix eigenvalues Ja, comparing analytical solutions from Eq. (6) and numerical iterative training using Eq. (2), both starting from the same initial condition (J) = 1). The analytical and numerical results match perfectly, demon- strating the expected time-scale separation: eigenvalues Ja corresponding to stronger covariance modes converge faster to their fixed point y((²) = = = =1/C,, which are the small- est in the coupling matrix, while weaker modes converge slower. Starting from an initial condition (0) that does not commute with C the coupling matrix must initially align its eigenvectors with those of C", process detailed in Appendix C and guided by Eq. (5). Following this align- ment, the eigenvalues evolve independently according to Eq. (12), supporting our analytical approach. Notably, train- ing directly from the population matrix achieves perfect reconstruction of the original model, thereby avoiding any discrepancies or generation errors as expected.
This analysis shows that with limited data, there is opti- mal training duration beyond which model inference accu- racy declines, highlighting sweet point for early stopping, However, detecting this point without ground truth is chal- lenging: it does not coincide with the peak of test LL (as in (b)), a phenomenon also noted in RBMs (Decelle et al., 2024). Moreover, the generation's 's quality, given error be- tween C(t) = Jt) and the population matrix C* in (c) (computed as Ec I CC) (l),,, stabilizes well before mii (p) and remains flat afterwards. This suggests that the generation quality of the GEBM isn't 't solely dependent on the model itself, as evidenced by consistent generation er- rors at both the minimum-erre point tmin and the training's fixed point, indicating this metric fails to capture the deteri- oration of model parameters over time. These optimal times are shown against p in Fig 3-dd). Additional evaluation metrics are discussed in Appendix F.
Impact of Finite Datasets: Interplay of Initialization and Time Scales Favoring Early Stopping Strategies. We explore the training dynamics using finite-data estimates of the population covariance matrix, CM,, With any finite M. the GEBM trained with ŌM will show discrepancies from the true model J'. We track these discrepancies by computing the reconstruction error & between J and the
Fig. 2-(b) illustrates the evolution of eigenvalue Ja over time for a sample size with 2=1.11. Initially, the stronger
5

---

A thcoretical framework for overfitting in energy-based modeling
N 100
N
200
N-300
RMT
Asymptotic analysis. Our findings so far have been es- tablished by numerically integrating the gradient ascent dynamics (i.e. using Eq. (2) with a slow learning rate), or with the analytical expression for the eigenvalue evolution (Fiq. (6)) In both cases, we utilized empirical covariance ma- trices extracted from a finite number of samples M, sampled from the distribution (I) with finite N. These results are almost insensible to the choice of the population spectrum as discussed in Appendix F.
6x 10
1.5
10a
100
Togt/ Jan/N
E(J
10
10
(a)
101
max--LL
10 1
10P
10
min 5
E(J", , Jnnseree)
6 10
p 1.5
10-1
maa
m 8
x 109
10
2
10"
2 10
(c)
10
(d)
We demonstrate that the phenomena of overfitting and finite- M corrections can be accurately modeled using RMT to predict the N, M + co limit, thereby removing the need for empirical data. Detailed methodologies are provided in Appendix G. For a constant = M/N = 1.5, Fig. 4-(a) compares the reconstruction error of (t) over the train- ing period for various N values (colored lines) against the asymptotic RMT prediction (dashed black lines), showing strong consistency as N increases. This agreement extends to the evolution of the test LL (not shown) and the timing of the minimum error and peak test LL as functions of p (see Fig. 4-bb)). Notably, the optimal stopping times for the two estimators do not coincide, yet finite (N, M) trainings align precisely with the asymptotic predictions.
10
X
10
10"
10
10=
P
Figure 4. (a)-(e): for p= M/N- 1.5 we plot the reconstruction enror during training (in (a), () and the final reconstruction obtained using a L2-norm regularization (in (c). vs A). (b): train- ing time achieving optimal reconstruction error (points) and time of maximum test I.I. (squares), plotted vs p. (d): optimal value of regularization prior À vs p. again selecting the optimum w.r.t. reconstruction eιror and w.r.t the test LL. All panels show com- parisons between numerical results for various N (colored lines) against asymptotic results from RMT (black line).
modes (deep red curves), which are less influenced by low- M induced noise, quickly stabilize, aligning their eigen- values aM close to the population values c marked by black crosses. This early alignment to very small and rela- tively accurate values makes the error curve (Ei) for finite M closely resemble that of the M + scenario. However, after period around fiin (p). the model starts encoding the weaker modes (blue lines), which are systematically under- stated relative to the population, causing J(x0) =1 1/ca to significantly exceed the ground-truth J = -///-- If train- ing begins from small Jn values, there's a critical point where the cigenvalues temporarily align more closely with their ground-truth than at the fixed point, effectively cre- ating an optimal time t% where Jo (42) RV Já. This align- ment markedly decreases discrepancies between the trained model's eigenvalues and those of the true model, highlight- ing the significance of initial conditions in training dynamics. Yet, the specific initial values of (0) are less critical, as long as they are substantially smaller than 1/2M for the weaker modes of ÕM (see Appendix D for further details).
5. Protocols to mitigate overfitting
In the GEBM, non-monotonic behavior stems from adjust- ments to the eigenvalues of C.M compared to the popu- lation covariance matrix. In fact, one can easily check that replacing the population eigenvectors while retaining M M-dependent eigenvectors to form an optimally corrected del
matrix, Õ. non-mone m--ootoooiccc effects on model quality and overfitting, as shown in Fig. (a) (green), and in other datasets or p values we see the bump completely disappear. While effective, this approach is useless for real experiments where the popula- tion matrix is unknown. Nonetheless, this idealized scenario informs the design of protocols aimed at minimizing overfit- ting and reducing reliance on uncontrollable early-stopping strategies. We now explore common strategies to mitigate overfitting within our framework, focusing on regularization and shrinkage corrections. We also introduce versatile downsampling-guide mode-fitting scheme that allows cir- cumvent the traditional limitations of RMT strategies, and design corrections that should be valid beyond GEBMs.
val-pop = Σ. ca u M u almost eliminates the
Now, we can also explain the stable generation performance of the GEBM, shown in Fig. 3-(C), using scale separation ar- guments. Generation error mainly depends on the strongest e values, which are learned early on, whereas the over- all model quality is controlled by the weakest aM (where ((>>) = 1//2/>, ), which minimally affects generation error due to their small value. While this phenomenon appears unique to the GEBM, a similar effect is observed in binary pairwise EBMs (cf. Sec. 6).
Regularization. In machine learning, regularization pri- ors are standard for preventing overfitting. In the GEBM, they constrain the growth of eigenvalues Ja. avoiding sub- optimal fixed points affected by mode fluctuations in CM For training dynamics, L2 regularization is applied to the coupling matrix J, and similar outcomes are achieved with projected L Li-reuularization on J's eigenbasis, (a protocol that facilitating asymptotic RMT analysis). The impact of
6

---

A thcoretical framework for overfitting in energy-based modeling
MI
GM
(x(99))
shrinkage fit
(a)
10
10*
ElJ Smiterea)
(c)
10-1
10-1
E(J" E(J", J(C)
2.0
2.2
2.8
E E(J, Jufemma)
10.0
(b)
--
RMT
EIJ",
1.8
1.6
10
10-2
1.4
1.0
10P
108
10
10
10P
107
1.2
num. updates
100-
2.8
mumerri
analytic
10"
10
10"
10
10a
(b1)
(b2)
(b3)
Juft)
Figure 5. Effect of data-correction protocols on the training GEBM (in (a), for p 2.8) and on the final model's quality as function of p (in (b)): comparison of the reconstruction error E between training from an empirical covariance matrix CM (blue), optimal L2 lz-reuularization (W.rt. reconstruction, red), shrinkage formula (cyan). The settings are the same as in Fig.3.
5.2
141.38
10
10
10 10"
10P
10
10
10a
10
Dum. updates
mim. updates
mum. updates
Figure 6. Results on the BM for the inverse Ising problem. (a)-(b): Training dynamics. (a): Reconstruction error (Frobenius norm) vs time (number of updates), for different values of dataset's size. (b1) (bb))-)))))3): evolution of eigenmodes Ja during training for 3 values of p: comparison between numeries (blue lines) and analytic curve (red). Black crosses indicate the true models' eigen- values BJa. (c): effect between data-correction strategies on the inferred model. Comparison between optimal Izzregularizztion (red), standard ML-training (bluc), modes-fitting (orange) and best reconstruction computed at fain (green).
regularization at the fixed point is studied for finite N and in the N + oc limit via R.T. Fig. 4 (c) shows the reconstruc- tion error function of X, with empirical results aligning closely with RMT predictions. An optimal Aopt minimizes the error but does not match to the value that maximizes the test LL (see (d)). complicating Aopt 's identification without knowing the population parameters, akin to identifying opti- mal early stopping. Further details on regularized training and RMT are provided in Appendices H.1 and G. Red line in Fig. 5-(a) illustrates the error over time for a L2 regularized training using the optimal parameter xopt
We then use the extrapolated eigenvalues to clean our co- variance matrix and run a new training. The evolution of the reconstruction error is shown in Fig. 5-(a) (yellow).
Shrinkage correction protocols are pivotal in statistical learning and signal processing for estimating covariance matrices, particularly when the sample size is small relative to data dimensionality (Bun al., 2017). Some of these pro- tocols use rotationally invariant estimators (RIEs) to adjust eigenvalues distorted by sampling noise (Ledoit & Wolf, 2004; 2020), while preserving eigenvectors, ensuring cor- rections are independent of the coordinate system. Based on RMT, RIEs align the eigenvalues of finite-sample covari- ance matrices to minimize the deviation of the covariance matrix from the population one. Using the optimal RIE from (Bun et al., 2017), we correct our CM matrices, and use them for training our GEBMs. Depicted in light blue in Fig. 5-(a). this new training shows significant improvements in model inference quality, although some non-monotonic behaviors persist. The main drawback of this approach is that it is specific to the GEBM case.
Comparison of strategies. The effectiveness of various strategies to counteract overfitting is depicted in Fig. 5-(b), presenting the reconstruction error across different p values. Notably, the optimal L2 regularization, the C.M val- strategy, and the performance at the optimal early stopping point derived from RMT all follow similar trajectories, with a 1//p scaling for large p as expected. While the high per- formance of these strategies stems from knowing the true model to optimize parameters not usable in practice we demonstrate that similar performance can be achieved with RMT-based shrinkage corrections or empirical polynomial fits. These methods do not require prior knowledge of the model, making them especially suitable for real-world infer- ence applications.
6. Boltzmann Machine for inverse Ising
Polynomial fit of eigenmodes. To overcome the limitation of RIEs, we introduce a simple strategy to correct empiri- cally the eigenvalues of C.: the idea is to downsample our dataset to obtain 8 with m < M and use the correspond- ing eigenvalues 2h values to extrapolate the m + oo limit from a linear fit in 1/m (as expected from (Baik & Silver- stein, 2006)). Additional details are given in Appendix H.2.
We extend our analysis to the Boltzmann Machine (BM) the so-called inverse Ising problem (Nguyen et al., 2017), adapting our approach to binary variables æ x= ( This model is able to capture multimodal distributions through its pairwise energy function E(x) = =-Etey Jyyaa, Σ, her, with parameters 0 = (J,h). Due to the lack
7

---

A thcoretical framework for overfitting in energy-based modeling
of closed-form solutions for the correlation functions and likelihood in BMs, employ a mean-field approximation suitable e.g. at high temperatures. This approximation allows for an analytic, albeit not exact, expression linking the model's correlation matrix C' to the coupling matrix J as C' = (Lv- - J)- facilitating an analytical exploration of the training dynamies shown in (Agoritsas et al., 2023) and further elaborated in Appendix I.
merical results (eigenspectrum of J()) with the analytic curve from Eq. (12). While the trends align qualitatively, particularly in capturing the separation of time scales, the analytic curves consistently underestimate the actual eigen- value evolution due to the overlooked diagonal constraint in the BM model. Nonetheless, this timescale separation is similar to that observed the GEBM: stronger covari- ances (ea > 1 ++ Ja > 0) are learned faster, while weaker covariances (Ga < ++ Ja < <0) take longer. This pat- tern indicates that the training dynamics are dominated by the convergence of weak e In sparse Ising models, this involves learning the negative spectrum of J. Unlike in GEBMs where negative eigenvalues do not exist, in BMs, these later-encoded modes significantly impact the overall reconstructed J due to their large absolute eigenvalues, even though they have negligible effect in sampling quality, Ac- curately inferring sparse Ising models hinges on effectively learning weaker covariances, heavily influenced by finite- data noise. However, training good generative models is considerably faster, see Appendix I.
Similar to GEBMs, analytical description of spectral dynamics can be applied to BMs, though with certain lim- itations due to two main factors: a) the ML estimator for the coupling matrix (t) may not strictly preserve the same eigendecomposition as ČiM, despite typically observing nice alignment for the strongest modes; and b) the binary nature of the variables is not accounted for in the diagonals of the covariance matrices. We must remove the diagonal constraint to allow independent evolution of the modes, sim- ilar to the GEBM scenario, since a spherical constraint as in (Fanthomme al., 2022) would introduce mode coupling. This decision is crucial for deriving approximate analytical expressions for training dynamics; without it, the problem becomes intractable in time. However, the proper fixed point alone can still be effectively analyzed using mean field techniques (Kappen & Rodríguez, 1998).
Strategies similar to those used in GEBMs can be employed to mitigate overfitting in BMs, with comparable outcomes as shown in Fig. 6-(c). Shrinkage formulas are not applica- ble to BMs, yet the empirical polynomial fit correction for the eigenvalues proves still effective at reducing overfitting effects. This is a very good outcome as it is the only non- informative correction (as the identification of tmin or Aapt requires knowing J*).
As detailed in Appendix I, we can project the gradient on the spectral basis of CM and obtain approximate analytic expression for the evolution of the eigenvalues of J:
dJa dt
(11)
I-Jo
7. Theoretic extension to generic EBM learning
Wo BBe Bee- (cM))) Z
(1) 2
(12)
whose fixed point is J ~ 1 (22)) This fixed point is shifted due to our unconstrained diagonal, and neither is J traceless, as compared to the complete treatment. However, Eq. (12) still qualitatively captures the training dynamics of BMs, revealing significant differences from the GEBM sce- nario where Ja and ca are no longer inversely proportional. Notably, the smallest ca values are associated with negative Ja values which are not necessarily small in absolute terms, which significantly contribute the reconstruction of J. Additionally, for positive Ja. Ja increases when 2M does.
The analysis of overfitting for simple models such as the GEBM or the BM constitutes a preparatory attempt to ad- dress this question in the broader context of EBMs. Let us see now to which extent the analyses of overfitting carried out so far is also relevant for more general EBMs. The line of arguments bears some similarity with the one justi- fying that high-dimensional linear regressions are relevant to analyze deep learning (Belkin et al., 2018; Hastie et al., 2022). To this end let us consider the score-matching al- gorithm (Hyvärinen & Dayan, 2005) theoretical proxy for the analyses of overfitting in EBM. Even though this approach might appear sub-optimal in many circumstances we postulate that the mechanisms leading to overfitting have similar origin as in more sophisticated methods. Consider a generic EBM of the form p( p(x(0) = Z z- (0) e -E( 0) where 0 E RR. is the vector of parameters and having a train set D {æ,, i = I,... M) } of size M. Defining the score function as (xx10) def -V, (0)) E(x())))
Results. We conducted numerical experiments training an Ising-BM on equilibrium data sampled from a 2D Ising model (i.e. defining J on a 2D nearest neighbors lattice) with N 8 x 8 spins at high temperature (3- = 0.1). Fig- ure 6 presents the results: Panel (a) shows the reconstruction error between the trained model and the ground truth 3J" for different p values, revealing non-monotonic trend at low P that mirrors observations made with GEBMs (a be- havior which is robust w.r.t. the system size, see Fig. 20 in Appendix I). Panels (b1)-(b3) track the eigenvalue evo- lution during training for three P values, comparing nu-
the score matching loss is given by (0 = 2 _ (a)))) - V log p(xx) [*]
27
8

---

A thcoretical framework for overfitting in energy-based modeling
which, thanks to a by part integration rewrites
where ((r) ( (a), 1,... M} is the vector
of empirical features spanning HE. Additionally 0g is di- rectly read off from (b(=)00)) at first order in O, in the lazy regime
(AEE(=))))(AA((((l)))) (8) Em | [[E([10))) AF(2))) + Cst. (13)
Notice first that for the GEBM, this leads to learning dynamics of the coupling matrix corresponding to dJ(t) dt
(x)10)) R V (s))0 (0: - O2o)
(20)
-(CJ(() + J(I)C) I leading to the solution
Using the parameter-sample duality eventually leads to
0: 0 + Je Jc(C Com) ge0
(x)
(14)
(21)
when assuming that the initial condition commutes with C. More generally, the dynamies of the score function is governed by a neural tangent kernel (NTK) (Jacot et al., 2018). We have dxp( (b)) di -Exr (r( x²) (1b((()))))) - V log p(x'))
where
M
C
V Vorl (0r Vor
(22)
M
(15)
gf def
Vou Tx(xs)
(23)
with Ke x²) = (x: (00(2(00)000 (z) Integrating by parts the second term we obtain
In GEBM case we recover (14) by letting (/x = Ox.
(0.))) -Exr Ke( x') (v(2')6))
K(x, x') 2 (a:x' xx x²) and ((z)) x. leading to ((M) = C with CCM)- = 1 M x
(z)) (16)
dt
M
with (x) -iy Ogr (ar))a())))))) (x010)
8. Discussion
EB Var (i( x²)] As for supervised learning we expect kernel regime for large enough network's width (Chizat et al., 2019). Then K becomes deterministic, the dynamics is linear, with (xx t) an explicit function of the kernel matrix K(x,, on the training set. Indeed, the NTK dynamics takes place on a reproducing kemnel Hilbert space (RKHS) of finite dimension corresponding either to Hp Span (x) q = I,... P). or to the HM del Span { (xx, v,), = I,. depending respec- tively on whether we are the under or over vee-pparameteriee regime. In the latter case K.s" def l K (æ., )x)) is full rank and we have
This work presents theoretical framework to understand overfitting in simple energy-based models, using the eigen- decomposition of the data covariance matrix to analyze training dynamics. We illustrate how the principal com- ponents control a timescale separation, where information progressively encoded from the strongest to the weakest data modes. Due to varying impacts of finite-size noise on different components. this results in an early-stopping point dictated by their interplay. Furthermore, we show that finite sample corrections can be very accurately described analytically using asymptotic RMT analyses. This analysis provide us with an analogous of the GCV in the context of EBM, which deserves further empirical investigations. This analysis exact for Gaussian EBMs and approximate for Ising-BMa -Ms at high temperatures. capturing similar phe- nomena observed in more complex EBMs like RBMs. We discuss data-correction protocols typically used to mitigate overfitting and propose to extend these strategies to more complex models leveraging higher- order data correlations (e.g. by exploiting the SVD decomposition). Further inves- tigations into RMT may clarify how early-stopping points relate to the asymptotic properties of the covariance ma- trix's spectrum or which should be the proper observables to pinpoint them without prior knowledge of the data model. Finally, an extension of the theory to EBM via a neural tan- gent kernel dynamics of the score function deserves further experimental investigations to find relevant hypothesis for the spectrum of population covariance matrices of tangent features.
M7
1 k((-to)
((1)
B ((oo)
(17)
K
where $(t) and 5 are respectively the vectors (1(2x)),), 1,... M) and (ő(x.), = 1....M) .) and assuming (to) = O. In any case we can consider only the projection of 1p the RKHS, its transverse part being assumed to be zero at L = lo- As result, the dynamics takes place in the "empirical" RKHS and we have
1 M M
11 Σ
K
(18)
where the vector 3 is obtained from (17) yiclding finally
0xx)0,)) -xx) [k)
3 (20
(19)
K

---

A thcoretical framework for overfitting in energy-based modeling
Acknowledgments
Belkin, M., Ma, S., and Mandal, S. To understand deep learning need understand kernel learning. In proc. of ICML, pp. 541-549, PMLR, 2018.
Authors acknowledge financial support by the Comunidad de Madrid and the Complutense University of Madrid through the Atracción de Talento program (Refs. 2019- T1/TIC-13298 & Refs. 2023- 5A/TIC-28934), the project PID2021-125506NA- financed by the "Ministerio de Economía Competitividad, Agencia Estatal de Investi- gación" (MICIU/AEI/1O. 13039/501 100001100000 10001 1033), the Fondo Europeo de Desarrollo Regional (FEDER, UE) and the French ANR grant Scalp (ANR-24- CE23-1320).
Bengesi, S., El-Sayed, H., Sarker, M. K., Houkpati, Y., Irungu, J., and Oladunni, T. Advancements in generative ai: A comprehensive review of gans, gpt, autoencoders, diffusion model, and transformers. IEEE Access, 2024.
Béreux, N., Decelle, A., Furtlehner, C., and Seoane, B. Learning a restricted Boltzmann machine using biased monte carlo sampling. arXiv preprint arXiv:2206.01310 2022.
Impact Statement
This paper aims to advance the field of Machine Learning by deepening our understanding of generative models under data scarcity. While our findings may have broad societal implications, we do not identify any that require specific emphasis this stage.
Béreux, N., Decelle, A., Furtlehner, C., Rosset, I.., and Scoane, B. Fast, accurate training and sampling of restricted Boltzmann machines. arXiv preprint S.Kiv:2400.15376, 2024.
Bun, J., Bouchaud, J.-P., and Potters, M. Cleaning large correlation matrices: Tools from random matrix theory, Physics Reports, 666:1-109, 2017. ISSN 0370-1573, doi: https:/loo.nr////00001/.//sysspp
References
Ackley, D. H., Hinton, E., and Sejnowski, T. J. A learning algorithm for Boltzmann machines. Cogmitive science, 9 (1):147-169, 1985.
URL https://www.cieencele com/ cocieene//rrticle/pp/////33/
63033377 Cleaning large correlation matrices: tools from random matrix theory.
Advani, M. S., Saxe, A. M., and Sompolinsky, High- dimensional dynamics of generalization error in neural networks. Neural Networks, 132:428 4468,44662226 2020.
Bun, J., Bouchaud, J.-P., and Potters, M. Overlaps between eigenvectors of correlated random matrices. Phys. Rev. E, 98:052145, Nov 2018. doi: 10.1 10.1103/PhysRevE.08 052145. URL https: / //t.///ink..///////// aps orppoorg//o//1// 1103 P Reys/ 98 052145.
Agoritsas, E., Catania, G., Decelle, A., and Scoane, B. Ex- plaining the effects of non-convergent MCMC in the train- ing of energy-based models. In Krause, A., Brunskill, E., Cho, K., Engelhardt, B., Sabato, S., and Scarlett, (eds.), Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Ma- chine Learning Research, pp. 322-336. PMLR, 23-29 Jul 2023. URL https: /ttp::/ss////pceepppeeee digs mlr. .tles// html.
Chizat, L., Oyallon, E., and Bach, F. On lazy training in differentiable programming. In proc. of NeurlPS, 32, 2019.
v202 10.02/agoritssa2222
Cocco, S., Feinauer, C., Figliuzzi, M., Monasson, R., and Weigt, M. Inverse statistical physics protein sequences: a key issues review. Reports on Progress in Physics, 81 (3):032601, 2018
Arnaboldi, L., Stephan, L., Krzakala, F,, and Loureiro, B. From high- figh-dimensional & mean-field dynamics to dimen- sionless odes: A unifying approach to sgd in two-layers networks. In The Thirty Sixth Annual Conference on Learning Theory, pp. 1199-1227. PMLR, 2023.
Consortium, G. et al. A global reference for human genetic variation. Nature, 526(7571):68, 2015.
Atanasov, A., Zavatone- Veth, J., and Pehlevan, C. Scal- ing and renormalization in hign-dimensional regression. arXiv preprint arXiv:2405.000 000900 2024.
Decelle, A., Fissore, G., and Furtlehner, C. Thermody- namics of restricted Boltzmann machines and related learning dynamics. Journal of Statistical Physics, 172 (6):1576-1088, 2018. doi: https://doi.orr/11..00 s -10955-018-2105-y.
Baik, J. and Silverstein, J. W. Eigenvalues large sample covariance matrices of spiked population models. Journal multivariat analysis, 97(6):3382 1408, 2006.
Decelle, A., Scoane, B., and Rosset, L. Unsupervised hierar- chical clustering using the learning dynamics of restricted Boltzmann machines. Physical Review E, 108(1):014110 2023.
Belkin, M. Fit without fear: remarkable mathematical phe- nomena of deep learning through the prism of interpola- tion. Acta Numerica, 30:203-246 2021.
10

---

A thcoretical framework for overfitting in energy-based modeling
Decelle, A., Furtlehner, C., Gómez, A. D. N., and Seoane, B. Inferring effective couplings with restricted Boltzmann machines. SciPost Phys. 16:095, 2024, doi: 10.21468/ SciPostPhys.16.44. 16.4.995 URL https: scipost. org/10. 2.1 / Sci Phys 4. 095.
Jacot, A., Gabriel, F,, and Hongler, C. Neural tangent kernel: Convergence and generalization in neural networks. In In proc. of NeurIPS, volume 31, 2018.
Kappen, H. J. and Rodríguez, F. B. Efficient Learning in Boltzmann Machines Using Linear Re- sponse Theory. Neural Computation, 10(5):1137- 1156, 07 1998. ISSN 0899-7667.
Decelle, A., de Jesús Navas Gómez, A., and Scoane, B. In- ferring high-order couplings with neural networks. 2025. URL hwtps:t///aaaatta org/ abs abs/2501.... 06108.
doi: 10.1162/
0899766963000113366 URL https:
Delon, J., Desolneux, A., and Salmona, A. Gro- was-wasserttiin distances between gaussian distribu- tions. Journal of Applied Probabiliry, 5984) 2022. doi: 10. 10.10177jp:22222.16.
Kiwata, H. Estimation of quenched random fields in the inverse ising problem using a diagonal matching method.
Phys. Rev. E, 89:062135, Jun 2014. doi: 10.1103/ PhysRewE.89.0621335. URL https: //link aps. http///...110///////eaaae . 062135.
Deng, L.. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine, (9(6):141-1) 2012.
Krizhevsky, A., Hinton, G., et al. Learning multiple layers of features from tiny images. 2009.
Du, and Mordatch, Implicit generation and modeling with energy based models. Advances in Neural Informa- tion Processing Systems, 32, 2019.
Ledoit, O. and Péché, S. Eigenvectors of some large sample covariance matrix ensembles. Probability Theory and Related Fields, 151(1):233-266 2011.
Fanthomme, A., Rizzato, F., Cocco, S., and Monasson, R. Optimal regularizations for data generation with prob- abilistic graphical models. Journal of Statistical Me- chanics: Theory and Experiment, 202255) may 2022. doi: 10. (0.1188/1742-54468888 //dx org
Ledoit, O. and Wolf, M. A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multi- variate Analysis, 88(2):665411, 2004. ISSN 0047-259X. doi: http:///o..ar///ooooolsss0//2x20 URL htttp::/www.scienccctt.. scfeenee/atticee//////0222200000
URL https: ac650c.
Feinauer, C. and Lucibello, C. Reconstruction of pairwise interactions using energy-based models. In Mathematical and Scientific Machine Learning, pp. 291-313. PMLR, 2022.
Ledoit, O. and Wolf, M. Analytical nonlinear shrinkage of large-dimensional covariance matrices. The Annals of Statistics, 48(5):3043 - 3065, 2020. doi: 10.1214/
19-A0S1921. URL https: 19 19-A01192.
Feinauer, C., Meynard- Piganeau, B., and Lucibello, C. In- terpretable pairwise distillations for generative protein sequence models. PLOS Computational Biology, 18(6): c1010219, 2022.
Loffredo, E., Pastore, M., Cocco, S., and Monasson, R. Restoring balance: principled under/oversampling of data for optimal classification. In Proceedings of the 41st In- ternational Conference on Machine Learning, 1CML/24. JMLR.org, 2024.
Furtlchner, C. Free dynamics of feature learning processes. J.Stat.Phy Phys, 190(3):51, 2023.
Golub, G., Heath, M., and Wahba, G. Generalized cross- validation as a method for choosing a good ridge parame- ter. Technometrics, 21(2):215 223- 1979.
MacKay, D. J. C. Information Theory, Inference, and Learn- ing Algorithms. Copyright Cambridge University Press, 2003.
Hachem, W., Loubaton, P., and Najim, J. Deterministic equivalents for certain functionals of large random matri- ces. 2007.
Magnus, J. R. and Neudecker, H. Matrix Differential Cal- culus with Applications in Statistics and Econometrics. John Wiley, second edition, 1999. ISBN 0471986321 9780471986324 047198633X 9780471986331.
Hastie, T., Montanari, A., Rosset, S., and Tibshirani, R. Sur- prises in high-dimensional ridgeless least squares inter- polation. The Annals of Statistics, 50(2):949 9866. 2022.
Mai, X., Liao, Z., and Couiller, R. A large scale analy- sis of logistic regression: Asymptotic performance and new insights. In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3357-3361. IEFF, 2019.
Hyvärinen, A. and Dayan, P. Estimation of non-normalized statistical models by score matching. JMLR, 6(4), 2005.

---

A thcoretical framework for overfitting in energy-based modeling
Marčenko, V. and Pastur, L. Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik, 1(4):457-483, 1967.
Saxe, A. M., McClelland, J. L.,, and Ganguli, S. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. In Bengio, Y. and LeCun, Y. (eds.), ICLR, 2014. URL http: dblp. uni-trier. de/ 4, html #SaxeMG13. 3.
Mei, S., Montanari, A., and Nguyen, P.-M. A mean field view of the landscape of two-layer neural networks. Pro- ceedings of the National Academy of Sciences, 1 115(33): E7665-E7671, 2018.
db/ con
Suzuki, M. and Kubo, R. Dynamics of the Ising model near the critical point. i. Journal of the Physical Society of Japan, 24(1):51 --560, 1968. doi: 10. 143/SSS1.24.51
Morcos, F., Pagnani, A., Lunt, B., Bertolino, A., Marks, D. S., Sander, C., Zecchina, R., Onuchic, J. N., Hwa, T., and Weigt, M. Direct-coupling analysis of residue coevolution captures native contacts across many pro- tein families. Proceedings of the National Academy of Sciences, (89449):E1233-11330 2011.
Tomasini. U. M., Sclocchi, A., and Wyart, M. Failure and success of the spectral bias prediction for Laplace kernel ridge regression: the case of low-dimensional data. In Chaudhuri, K., Jegelka, S., Song, L.,, Szepes- vari, C.. Niu, G., and Sabato, S. (eds.), Proceedings of the 39th International Conference on Machine Learn- ing, volume 162 of Proceedings of Machine Learn- ing Research, pp. 21548 21583, PMLR, 17-23 Jul 2022. URL www.t::p/tt......... htt.htttsss:/prrceeeee /ttp:///ccccceeesseeeeeee ml prrsss// v162 htsl//tmmsiini20 html.
Nguyen, C., Zecchina, R., and Berg. J. Inverse statistical problems: from the inverse Ising problem to data science. Advances in Physics, 66(3):197-261, 2017. doi: 10. 1080/ 00018733.2017.1341660 URL https: Reoi.org
Tubiana, J., Cocco, S., and Monasson, R. Learning protein constitutive motifs from sequence data. Elife, 8:e39397, 2019.
Patil, P., Wu, Y., and Tibshirani, R. Failures and successes of cross- validation for early-stopped gradient descent. In International Conference on Artificial Intelligence and Statistics, pp. 2260-2268. PMLR, 2024.
Wei, A., Hu, W., and Steinhardt, More than a toy: Random matrix models predict how real-world neural representa- tions generalize. In International Conference on Machine Learning, pp. 23549-23558 PMLR, 2022.
Potters, M. and Bouchaud, J.-P. A first course in random ma- trix theory: for physicists, engineers and data scientists. Cambridge University Press, 2020.
Wu, Z., Johnston, K. E., Arood, H., and Yang, K. K. Protein sequence design with deep generative models. Current opinion in chemical biology, 65:18-27, 2021.
Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F, Bengio, Y., and Courville, A. On the spectral bias of neural networks. In Chaud- huri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Machine Learn- ing, volume 97 of Proceedings of Machine Learning Research, pp. 5301-5310. PMLR, 09-15 Jun 2019. URL htt.s::p/// http://procccee http///oooceennnneneeeee mlr, press/v97/ rahaman19a 9a. html.
Yasuda, M. and Tanaka, K. Susceptibility propaga- tion by using diagonal consistency. Phys. Rev. E, 87:012134, Jan 2013. doi: 10.1103/PhyseevE.88 wtw.appp...ooooooooo/..oooooooop
012134. URL 1103/PhysRevE. 0120143
Yelmen, B., Decelle, A., Ongaro, L., Marnetto, D., Tallec, C., Montinaro, F., Furtlehner, C., Pagani, L.,, and Jay, F. Creating artificial human genomes using generative neural networks, PLoS genetics. 17(2)ee1093303, 2021.
Ricci-Tersenghi, F. The bethe approximation for solving the inverse Ising problem: a comparison with other inference methods. Journal of Statistical Mechanics: Theory and Experiment, 2012(08) aug 2012. doi: 1081088/ 1742-5468820012//8/P80000 URL https: http//xxxooo///oooo doi. org /0/ 108880088880
Yelmen, B., Decelle, A., Boulos, L. L., Szatkownik, A.. Furtlehner, C., Charpiat, G., and Jay, F. Deep convo- lutional and conditional neural networks for large-scale genomic data generation. PLoS Computational Biology, 19(10): 1584, 2023.
64 /tte/////////011////0000000
Roudi, Y., Aurell, E., and Hertz, A. Statistical physics of pairwise probability models. Frontiers computational neuroscience, 3:652, 2009.
Saad, D. and Solla, S. Dynamics of on-line gradient descent learning for multilayer neural networks. Advances in neural information processing systems, 8, 1995.
12

---

A thcoretical framework for overfitting in energy-based modeling
A. Derivation of projected gradient equations
Starting from the log-likelihood's derivative w.r.t a parameter Jıj, we can assume that in the limit of an infinitely small learning rate ? + 0 we can replace the discrete-time update equation for the parameter (2) into a differential equation for the evolution of each parameter Jay:
0C
dJy dt
OC
Jy (t (t+ J1y (t) +
(24)
0Jky j(()
0Jú J(()
We now decompose the rhs of the above expression in terms of time-evolution of eigenvalues and eigenvectors of J at time Given the eigendecomposition Jis- Σn of Jr we have
dJij d EojJuos Σ di
dop dl
d. di
de}
Jyu
+ r.
+ 1 Jr
(25)
di
dl
We now project this on the eigenbasis of the eigenvectors of J, which after simple algebraic manipulations leads to
dJs = d0d dJa
Σ duf
dus
Σ 1
+ (1- - 600) (1-6a0)
3g
+ Σ
Jau
(26)
dt
dt
dt
dt
dJa
du
+ (1 11--8ah San) (Jg Ja) Σ dad dt
(27)
dt
where we used the property (u = = 0 because they are vectors of an orthonormal basis. Finally, combining Eqs. (4) and (27) separating the contributions for a = and a + 3 we get to Eq. (5) the main text.
final note on the log-likelihoods' gradient: in the first expression (3) we have assumed that perturbations are not symmetric, that is when taking the derivative w.r.t. Jıj for i + j we assume that Jij Jji. Assuming instead symmetric perturbations one would get slight different form of the Log-likelhhoo''s 's gradient w.r.t. (3), given by:
OC JJ.g
Aij [- +ON+((1))) (J
(28)
with Aij 1 From the point of view of the training fixed point this is an issue, the ML estimator is exactly the same in both cases. However, the modified gradient (28) leads to slight different dynamics: in particular, it is not anymore true that the dynamics can be exactly decomposed into a separate evolution for the different eigenvalues of J by following the above steps. One could either include symmetric constraints on the eigendecomposition when computing the projected gradient (a more cumbersome process from a mathematical point of view, see e.g. (Magnus & Neudecker, 1999)) or simply double the learning rate on the diagonal terms = j to compensate for the factor Ary. Nevertheless, the difference between the analytic evolution (i.e. Eq. (6), obtained assuming mom-yymmttiic perturbation) and a numerical training performed using the gradient (28) are almost coincident as seen from Figure 7. None of the results presented the manuscript is affected by such difference in the gradient computation.
B. Finite-M fluctuations of eigenbasis
We detail additional eigenvalue spectra of covariance matrices from various datasets in Fig. 8, complementing those in Fig. of the main text. The spectra for CIFAR-I0 (Krizhevsky et al., 2009) and the Human Genome Dataset (Consortium al., 2015) are displayed in (a) and (b), respectively. Panel (e) illustrates the empirical covariance matrix from equilibrium configurations sampled from a 2D Ising model with periodic boundaries high temperature (3 0.1), i.e. in a paramagnetic phase. Panel (d) presents another synthetic spectrum, generated through mixture of power-laws from Eq. (10), but using different set of parameters than those used for Fig. 1-(b) (and used for the trainings in Figs. 2 and 3). This new spectrum was used for the figures involving comparisons with RMT (Figs. 4, 5), for numerical stability issues with the integration of the RMT equations. Nonetheless, all the results presented in the main text about the training dynamics of the GEBM can be perfectly reproduced on a wide range of the parameters defining the population eigenvalues, so that the qualitative picture that emerges from our analysis is extremely robust with respect to specific details of the spectrum.
13

---

A thcoretical framework for overfitting in energy-based modeling
zaxlytic symmetric perarbations)
mumetii (symmetric peraubatons)
[4]]
Jolt)
-0'
10
10
10
100 10
10a
10 10
72
Figure Difference in eigenvalues' evolution the training of GEBM when imposing symmetry or allowing asymmetry in the perturbation of Jej. The points correspond to numerical results obtained by enforcing symmetry on Jıy after each update during training (i.e. using Eq. (28)), while the lines represent analytical expressions derived for the case of non-symmetric perturbations (i.e. Eq. (6)). The setting is the same as Figure 2 in the main text.
(a)
CIFAR10
(b)
HGD
(c)
Ising 2D
(d)
Synthetic B
10
10a
M 4-
M 10
M 5-10P
10
M 102
M 10
10a-
M M=2-11²
M = 102
M =2 2210"
M 1.6 10
M- 3-0
M 3V-T8 7-10
M 10a
3 10"2_
M 3.2--00
M 102
M MM-22 1.5
10"
15 10
10-2
10P
10
10-4
10-1
M 5-10
M - 4.5-10
M ->.00
M -0
10-6
10-8.
10-1
10
0.5
1.0
0.0
0.5
1.0
0.5 a/N
1.0
0.0
0.5
a/N N
a/N
a:/ N
Figure 8. (a)-(b)-(e): Eigenvalue spectra of the empirical covariance matrix of real datasets, respectively CIFAR- 10 (in (a)). Human Genome Dataset (b)), and dataset made of equilibirum configurations of 2-d Ising model of size N 162 at 0.1 (in (c)) Black lines represent the spectrum computed with the full set of available data (of size M"), while scatter colored points show the result for subset of data M < MMM"" (d): the black line shows synthetic population eigenvalue spectrum generated according to (10), with N = 100, and r = 0.5, 8 = 1.0, = 0.5, x1 10
10; colored points display the eigenvalues of the empirical covariance matrix CM computed by sampling configurations from a GEBM with J = C- (from (1)) for different values M.
Fig. 9 illustrates how the eigenbasis of the covariance matrix for real datasets remains consistent against downsampling. Starting with the eigenbasis decomposition the covariance matrix for the largest available dataset M consideedd our closest approximation to the population matrix CM"" C C ww denote its eigenvector matrix as U = {u Ja=z These eigenvectors are arranged columnwise and sorted in descending order by their corresponding eigenvalue. For each reduced sample size M M we perform similar decomposition on the resultant empirical covariance matrix C. with MI
its basis represented as UM M = 2u J a=1- To evaluate the preservation of eigenvectors, we calculate the norm of the matrix product between a projection operator P defined as p Ui. (incorporating the first n eigenvectors of C*) and the a-th eigenvector of CM U.M. This measurement determines whether u falls within the subspace spanned by the first n eigenvectors of C", thereby helping to mitigate eigenvector oscillations due to exchanges between the ordering of the associated eigenvalues. Fig. shows the norm
pm T .u plotted versus n and each value of a, for the same datasets of Fig. 8. We can observe that for high values of M most eigenvectors eigenvectors of C are contained in the subspace spanned by C", as the norm of such matrix product raises sharply to 1 when n R a. On the other hand, when M is lowered (lower panels on each subplot) the conservation starts to deteriorate. especially in the middle-lower part of the spectrum. Interestingly, we observe that the most conserved directions (at least in (a)-(b) are both the strongest covariance modes and the lowest ones, a phenomenon already highlighted in (Bun et al., 2018).
well preserved: this means that most
14

---

A thcoretical framework for overfitting in energy-based modeling
MINSTT
HGD
764
M - 6300
indes
ca0do indess
PM
|
0.0
1.0
1.0
305
1000
moode Baddes
M Tm
|
||
0.0
800
projector index
ponjector indls 71
(c)
Ising
(d)
Gaussian Model
1.0-
250
1.0+
300
|p 1
26530
index
|| -u
30002
Booke indeoo
M
M1
0.0
1.0
incles
moode indec
|
M
M
150 projecsee inclec 29
250
projoctor Index
Figure 9. M fluctuations of eigenvectors in the covariance matrix of datasets. The four panels show the norm of the matrix product between the n-th projection operator P", containing the first n eigenvectors of the population matrix C" (for a real dataset, we just take the covariance matrix with the full available data M") and the a-th eigenvector of the covariance matrix CM with M < M (a)-(b)-(c) respectively refer to MNIST dataset (Deng, 2012), Human Genome dataset, and equilibrium configurations drawn from a 2d Ising model (same setting in Fig. 8-(c). (d) refers to synthetic Gaussian Model generated as discussed in the main text with the same settings in Fig. 1. All panels show the results for two values of M, a larger one at the top and a lower one at the bottom. Results are plotted w.r.t. the projector index and each line correspond to different a.
C. Training dynamics in GEBM with non-commutative initialization
This section provides a brief follow-up to what discussed in the first part of Section 4, concerning the training dynamics of a GEBM. For simplicity we focus here only on the infinite-sample scenario (i.e. when training from C*), although the same reasoning holds also for finite data. We are also interested in describing the training dynamics for a generic initialization of the matrix J, which in general will not commute with C In this scenario, the model has also to learn the eigenvectors of C Fig. 8 shows the evolution of the coupling matrix eigenvalues Ja according to Eq. (6) (shown with solid lines), in comparison
numerical training done iteratively maximizing the likelihood as in Eq. (2) (points). The initial condition here is a matrix (0) constructed from random population of modes (0) U10, 1] and projected on to random orthogonal matrix. In this way, (0) and C do not commute. At the beginning of the training, there is indeed a discrepancy between theory and simulations, because of the wrong assumption of independence of eigenvalues. Once eigenvectors align, the evolution proceeds independently for each eigenvalue and perfectly follows Eq. (6).
Note that the initial oscillations of the modee-to-mode eigenvector overlap in Fig. 10-(b) due to the fact that eigenvalue learning non monotonic at the beginning, so that there is an initial exchange in the ordering of the eigenvectors. Nonetheless, after initial transient all the eigenvectors align to their counterparts in the covariance matrix. This alignment process and with much faster timescale w.r.t. the learning of eigenvalues (especially the ones associated to weaker covariances): for this reason, the assumption leading to our analytic description about independency on the eigenvalues' evolution remains justified for practical purposes. This reasoning about eigenvector alignment holds for any input covariance matrix: what determines the non trivial dynamics of the reconstruction error (explained in Sec. 4) is fully determined by the noise in the eigenvalues.
15

---

thcoretical framework for overfitting in energy-based modeling
(a) 10!
100
10
Jo10
10-1
analytic
mumerie (e)))
(b)
0.8-
0.6
0.2
mummeric
10
10-1
10
7
Figure 10. Training dynamics of the GEBM from population matrix C The system size and the parameters defining C' are the as in Fig. 8-(d). (a): Evolution of eigenvalues, comparison between analytic solution (full line) and numerical training (points). The initial condition (0) constructed from random distribution of modes and projecting it back to random orthogonal matrix which differs from the eigenbasis of C so that the two matrices do not commute. (b): Aligmment of eigenvectors, computed as the mode-to-mode overlap between eigenvectors of the population matrix u,, and eigenvectors of J. i.e. wa. Red-ish (resp. blue-ish) colors correspond strong (resp. weak) covariances c. The learning rate is set to 10-3
D. Robustness of results w.r.t. initialization at finite M
We show in Fig. 11 some analogous results w.r.t. Fig. 2-(b) for the training dynamics of the GEBM in the case of finite M (here set p MM/N 2.11) by varying the initialization. In this case, do not care about cigenvectors' alignment as in the previous section: we only consider different initializations for the eigenmodes of the coupling matrix, i.e. [Ja) We can observe how the non-monotonic behavior of the reconstruction error (plotted in the bottom rows for each initialization) robust against different standard, uninformative and small initializations, and it keeps appearing as long as Jn (0) < 1/2 for the majority of the cigenvalues (especially the ones corresponding to weak covariances). We can also observe how, increasing the initial conditions higher values than the fixed point (i.e. moving from columns (1)-(2)(33) to the right most ones (4)- the non-monotonic behavior disappears, indicating that the early-stopping break-point longer exists and that there is way mitigate the overfitting effects with this strategy. However, is common practice start with small values. Moreover, in more complex EBMs where sampling is required to estimate the correlations of the model in the LL gradient (e.g., BMs or RBMs), it may be very bad idea to assume extreme initializations (i.e. far from an uninformative initialization where the parameters the model are small): this could lead to ergodicity problems in sampling, as the model may get stuck in spin-glass-lite phases, a phenomenon that has been well studied in several EBMs (see e.g. (Decelle et al., 2018)).
16

---

A thcoretical framework for overfitting in energy-based modeling
(0)
1.0
(,) 0.01
J,
U110
J,
3[1,30
J.(0)
50.0
(a2)
a33)
10
(.())
10
10-1
(al)
(a4)
(a5)
(b1)
(b2)
(b3)
(b4)
E(J J(B))
10
(b5)
10
10F
10
10a
10 10
10
Figure 11. Results on the training dynamics GEBM finite amount of data by varying the initial conditions. Panels (a)'s (top row) show eigenvalues' evolution (according Eq. (12)), while panels (b)'s (bottom row) show the corresponding reconstruction error w.r. the ground truth model J"; all quantities are shown versus time. Each column corresponds instead to different initial condition. From left to right: an identity-like initialization (J, (0) 1) in column (1), in Fig. 2-(b); small coupling initialization (J)(0) 100 in column (2); two random initialization of modes (resp. in the boundaries [10 " 10] and [1; [1:30] in column (3)-(4)); constant initialization very high values larger than fixed point, i.c. (Ja(0) = > 1/66 column (5). trainings are performed analytically, with an
empirical covariance matrix GM generated with the same settings as in Fig. 2 with p = M/N = 2.11.
E. Additional generation quality metrics
In this section consider an additional metric to compute the discrepancy between the trained model and the true one, namely the Wasserstrein distance (Delon et al., 2022). Figure 12 shows the same data as in Figure of the main text, now including the evolution of the Wasserstrein distance between the trained model and the true one w.r.t. training time. Also this quantity shows non-monotonic behavior in f especially for low M, with a clear carly-stopping point. the rightmost panel, we compare the locations of the minima of each error estimator (and the maximum of the log-likelihod)) as functions of p. We observe that the time point corresponding to the minimum Wasserstein distance follows trend very similar that of the maximum likkeliihood.
F. GEBM analysis with various eigenvalue spectra
This section presents additional results on the GEBM, analogous to Fig. in the main text, obtained using alternative spectra for the population covariance matrix. To assess the robustness of our findings with respect to spectral choice, we replicate the analysis using both synthetic and empirical spectra.
First, in Fig. 13, consider synthetic spectrum distinct from Eq. (10): for N 100, generate 10 dominant modes with amplitudes uniformly distributed in [2, 10), and a bulk of N k 90 noisy modes with amplitudes in [10-1, 1]. The qualitative behavior, including overfitting effects, remains unchanged.
Next, we repeat the analysis using empirical spectra: specifically, the eigenvalues of the sample covariance matrices from the MNIST and Human Genome Dataset (HGD), shown in Figs. 14 and 15, are used as population spectra. The resulting dynamics, analogous to Fig. 12, again show no qualitative deviations. all cases, the non non-monotonic temporal behavior of key metrics and the early-stopping times preserved.
17

---

thcoretical framework for overfitting in energy-based modeling
(a) Reconstruction error
(b)
test log-likelihood
(e)
10"
10
10a
p
M/N
early stopping times t
102
10
10
10a
10
10
10
10a
10
104
training time
training time
(c)
Generation error
(d) Wasserstein distance
100
minmmmm reomeiurciio erner
sttady generation erear
10
10
10a
10
10
10
10a
10
100
10
10a
training time
training time
P M/N
Figure 12. We compare the results shown in Fig. 3, obtained using various generation quality measures, with the corresponding curves computed using the Wasserstein distance. (a)-(b)-(c) display respectively the reconstruction error &, the test-LL, the Wasserstein distance and the generation егтоr &c, plotted time, for various sample sizes M (indicated by a color gradient from blue for increasing p= M/N. Dashed black lines refer to training from C* (i.e. M -> oc). (c): comparison between time of minimum reconstruction (circles), maximum test LL (diamonds), minimum Wasserstein distance (squares) and time at which generation emoor converges to its steady-state value. These quantities also shown in the related panels for better clarity, Apart panel (d), this figure contains same information and quantities as Fig.
(a) Reconstruction error
(b) test log-likelihood
(e)
10"
10'
102
p
M/N
early stooo dooppig times t*
102
10
10
10a
10
10
10
102
10
104
training time
training time
(c)
Generation error
(d) Wasserstein distance
-00
mimeumm rroonstruction cener
simudy gunnttion
10-2
10
102
10
10
10
10a
10
100
10
102
training time
training time
p == M/N
Figure Same plots in Fig. this time obtained by training GEBM starting from synthetic population covariance matrix spectrum of dimension N = 100. with set of 10 dominant modes with amplitudes uniformly distributed in the interval [2, 10), and bulk 90 noisy modes with amplitudes uniformly distributed between 10 and
18

---

thcoretical framework for overfitting in energy-based modeling
Reconstruction error
(h)
test log-likelihood
early stopping times t
10
10'
10a
P
M/N
109
10a
10
10
10
10
10
reonsiruction
training time
training time
maddums tos
steady generation
(c)
Generation error
(d) Wasserstein distance
107
mastarun Wizsesstrcin distance
105.
10 training time
10
10
10
10a
10
10
102
training time
p = M/N
Figure 14. Same plots as Fig. 12, this time obtained by training a GEBM starting from the eigenvalue spectrum of the empirical covariance matrix computed from the MNIST datasct, with cutoff 10-6 to filter out weak modes.
a)
Reconstruction error
(b) test log- lilkiiiiioodd
(e)
108
early stopping times t*
10
10"
10
10a
0
M/N
106
10a
105
108
10a
10a
108
reometmution maxionms test
training time
training time
10s
stedd genenation moncarum Warneestre distance
(c)
Generation error
(d) Wasserstein distance
104
103
102
10a
105
108
10a
10a
108
10t
102
training time
training time
p -= D-M/N
Figure 15. Same plots as Fig. 12, this time obtained by training a GEBM starting from the eigenvalue spectrum empirical covariance matrix computed from the HGD dataset, with cutoff 10 to filter weak modes.
19

---

A thcoretical framework for overfitting in energy-based modeling
G. Asymptotic analysis through Random Matrix Theory
G.1. General case
Various quantities appearing in the core of the manuscript are explicit function of the empirical covariance matrix C. and as such are amenable to asymptotic analysis thanks random matrix theory (RMT). These quantities are respectively the train, test energy (associated to the EBM), the coupling error and the LL (train and test). For the sake of clarity, we repeat them here:
1 Tr(J
Errain
(29)
N
1 1 Tr(JO*),
Erest
(30)
N
1 ||J J 11
Eg
(31)
N
L
log cle([]]
Eirain,tests
(32)
2N
where C limm C. is the population matrix, || }| the Frobenius norm while J is the estimation of the coupling
matrix from the train samples x assumed to be of the form Fz, with E[zz =I, FFT C.M ||ar|| distributed
w.r.t some density o(T). Depending on the setting (dynamical, spectral Ê, or L2) J may appear in three different explicit
(Ln) and j(رر of C.M. We have
functional form jt. Jn
(x) = (1+ Wol Wol-o- training dynamics,
(33)
((,1),
a
(L, (spectral) regularization),
(34)
1 ax
)((2) a
(L2 regularization).
(35)
derivation Eqs. (34)-(35) given in Appendix H.I. The j: corresponds to the situation where all eigenvalues Ja have the initial condition J, (0) 0 and follow the time evolution Eq. (6). Let call generically j the functions given above.
Using the resolvant
G (3) def
zI
we can express the various quantities of interest with help of Cauchy integrals
Exrain 1
2in ff (3)))
Frest 2im ff [+
2i f. dz (((3)))(((M)))) ("0
1
E, =
1 of dz logli(=)]
Tr Tr(G)
Etrain,test
LLtrain,
2im Jc
2
where C is a contour of integration around the real axis. Next, from RMT, the proportional asymptotic limit M, N - 00 with fixed M/N = p. GMM) has a deterministic equivalent (Hachem et al., 2007) G. defined as
1
G(2)
-I- -(=)C*
20

---

thcoretical framework for overfitting in energy-based modeling
with C* the population matrix and (=) is given implicitly by the following self-consistent equations (Marčenko & Pastur, 1967)
A(2) = f - [())T'
a(T)
(36)
v(dz)r T(2) 2 A(z)x'
(37)
with vddz) the spectral density of the population matrix and where
lim a N.
TT[ M
(38)
For sake of clarity we do not consider the Marchenko-Pastie equations in full-generality and actually assume the fluctuation of z to be negligible i.e. take (77) 5(T - 1). Letting p(dx) the asymptotic limit of the empirical spectrum in the proportional regime, its Stieltjes transform is given by the trace of the resolvent:
/ v(dx)
((2)
2- x
Then the bulk spectrum is given by the Stieltjes transform
g + ie) = gr) + in {((1),
{
which rewrites (disregarding the pole at = 0 for p < 1)
pA:(v) p
(,(0)
f(y)
(39)
Ty
Ty
[1-r, T2)2²2
Along the contour integrate y ie with infinitesimal. In the limit O, both A and I' may acquire finite imaginary part which we write as
lim A(z) = A, (v) = Avv)
lim -(=)=T,(0 + (v))
c-0
In terms of these quantities we obtain the following equations for the train and test energies:
| dyj(w) (.(r)) + A,())
Etrain =
P π /0 dyj()) + (p<))((0)c()))
Frest
where c(p) given implicitly by
v(dx)x
(40)
+ c(p)
Errain may also be written
Exrain = =f dybi(y)7(y).
with F(y) given in (39).
The coupling error takes the form for any > 0
Eg
v
p(dy))²(y)
2
( v(dy)
(1- p) p) + 2p)r (w))
p p
y
(p<1) (1 p)) ())
1
(pp)
0
but practice we consider only the unram--arameterrieed regime corresponding to p 1.
21

---

A thcoretical framework for overfitting in energy-based modeling
G.2. Special of spectral 4 Regularization
The case corresponding the form (34) can be treated more directly without use of Cauchy integrals. In that case, considering instead the resolvant
G
I aCM
with the inverse penalty
Ā
introduced here for convenience. This leads to the following form of the various quantities of interest
Etrain
Tr(Gian (M)
(41)
N
Erest
Tr[GiMO
(42)
Beouplings
|aG - J]
(43)
I.Irainnttstt
2N Trilog o
Ferain,test
(44)
2
In the scaling limit we again have a deterministic equivalent (Hachem et al., 2007) of the resolvent the form
G I+AC*
where the fixed point equations now read (a(() = 5(T 1))
r v(dx)
(45)
0
+ Ax
A a
(46)
T+T
with again I' given by (38). The expression for Etrain,test are straightforward in the scaling limit:
Errain = 1 - / v(dx)
1 1+ Ar
r
Erest
p
Remarkably, thanks to leave-one out argument there is a deterministic relationship between the train and test energy. For 8 E Irrain, we have a leave-one out relation of the form
GLM
G (M)
1 + A oz (M)
where (M)) is the resolvant obtained after removing s from the train set. Assuming the samples to be the form Ds = F 2z 1. with FF C with z, -O, MTI) such that E E(z,z,) then for large N we have the concentration property
o G (an)
G M + +0
x.
M
As result for large N, M we immediately obtain
Exest
E Etrain
1 Frest
22

---

thcoretical framework for overfitting in energy-based modeling
which can be reverted
Etrain
Erest
(47)
Prrain
Concerning the error on the couplings, we have
1 N
E Tr -
1 [TTT(Grmm) + Tr[J 20Tr(G
= N
From Ledoit-Péchet the last term simply reads (up (1/vMM) corrections):
M Tr[
J-
Tr J M I+ AC
] M J++
Tr
J
For the first term use the following identity:
G GLM2 = =G(M) + +a- d
(
de
As result, asymptotically we have:
Tr G Tr
ta d
Tr
II AC
da
1 AC
1
N(a)TT
AC
(I+ AC
For this we need to compute A A((a) which can be done from the self-consistent equation (46,45):
A P
A'(a)
a² p Q
with
Q = 1 Tr
A²C
M
(T+ AC*)2
Ultimately we obtain:
E, 1 Tr (1-a0
M 1 AC c -1
1
a²(1 - N)
AC
Tr
M
+ AC*
So we have
- 02
a 1
/ u(dx))
Ar
E, = u(dx)) [ Ax
(1 + +xr)2
Finally, concerning I.rrainnttsst we don't see how to avoid the Cauchy integral, but the train-test relationship (47) has an important consequence, because it allows us to get a very precise estimation of the test likelihood when M becomes large:
Errain
LLrest (J) log dct(()
1 - p Euain
1
long as J is the function (34) of C(a)), For general EBM models have LL of the form
Litrainttet [J] log Z[J] tesain,test [J]
so by analogy with GCV, not excluded that we can use this train-test relation in practice.
23

---

A thcoretical framework for overfitting in energy-based modeling
H. Details on data-correction protocols
In this section, we analyze the effect of different ways to improve the estimation of the covariance matrix's eigenvalues in order to avoid diminish the effect of overfitting during the training dynamics.
H.I. Training dynamies with regularization prior for finite N
We first discuss what happens to the training in the presence of a regularization. We employ two regularization protocols: a standard norm, and projected L1 -ormm.. The choice of the second regularization is justified because it allows have aasdimmum---osteeiiii coupling matrix which commutes with the original covariance matrix C.M, it happens in the absence of regularization, thus facilitating the asymptotic analysis through RMT discussed in Appendix G.
1.1.1.1 L2 REGULARIZATION
The log-posterior now reads
1 M
- T (J²) (TP) log (J (JI LD (J)()()))))))))))))(
(48)
where Ả is the reguavlarization strength. The derivative w.r.t. the parameters now reads
1 Plogp( D) = [- + (J (J - NWo
(49)
M
OJij
Notice that the new term commutes with the second one (J and J-1 are diagonal in the same basis), so even in this case the maximum-at a-osterrioo matrix JMPP will share the same basis as C.M, as it happens in the absence of regularization. Therefore, we can apply the same reasoning discussed in the main text and project the log-posterior's gradient on the basis of J. The evolution equation of each eigenvalues reads:
dJa
1-2M
-
1
2 AJa.
(50)
dt
J.
Although there exist no closed expression for the full time-dependent solution of Eq. (50), it is possible least compute analytically its fixed point:
(( 2A
(2)
M
V (eM) + +AA
(51)
The full coupling matrix corresponding to the above fixed point is finally computed projecting back Eqs. (51) onto the cigenbasis of C that is Lox) (A) E. -!00 (A) u u
H.1.2. SPECTRAL Ż NORM
This regularization schemes utilizes a L1 -nomm but on the projected basis of the coupling matrix J. This construction still allows to employ a similar formula to Eq. (12) to describe the evolution of eigenvalues, each one independently on the others. The original differential equation describing the evolution of Ja is now modified as
1 dJa
1
GM
X
(52)
7 dt
Jn
whose fixed point reads
((00) Ea (A)
(53)
e. + X
Note that equations (51)- (53) just derived are the same ones Eqs. (35)-(34) in Appendix 6, respectively.
1.3. RESULTS ON THE EFFECT OF REGULARIZATION
We can check the performances of either type of regularization by looking at the training fixed point, and at how it modifies the quality of the inferred model. In Fig. 16-(a), we show the reconstruction error E, computed between the ground truth and the inferred model the presence of a regularization prior with strength A, both for the L2 (solid lines) and for the spectral-T., norm (dashed lines), for a given value of number of samples: here we have = M/N 1.5. Comparison
24

---

A thcoretical framework for overfitting in energy-based modeling
(a)
1
L
L N = 100
m1 zun E
mate-L min-S
6 10
L N = 200
p= 1.5
10-1
E(J, Jalured)
L = 300
4 x 10
T RMT
3 10
10
12- N 100
2x 10
L-- 200
Lg- 300
10-3
Lg RMT
10
10-3
10
10
10
10" 10'
10 10
10
10
A
P
Figure 16. Effect of the regularization priors on the inferred model's quality. (a): the plot shows the reconstruction error & computed between the ground truth and inferred model the presence a regularization prior with strength A. Solid lines refer to the L2 norm, while dashed lines to the spectral-L, norm, both discussed in Section 1. For each prior, we compare finite-size results (colored lines) RMT asymptotic estimations. (b)-(c): we show the values the regularization strength App that achieves optimal reconstruction the model (lines with scatter points), and the optimal value of the regularization that maximizes the test log-likelihood (lines with scatter diamonds). Panels (b) (resp. (c)) refers to the optimal values when using the I-- (resp. the spectral T -nomm. Note that (b)-(c) are on the same y-scale. Settings the same in Fig. 4.
shown between finite-size trainings (colored lines) and RMT estimation (black lines). All the quantities are plotted versus the regularization strength A: we can observe how there is clear non-monotonic behavior with a minimum developing at a certain Appt- As one might expect the optima value differs between the two regularization priors, i.e. X * X Results on both priors are shown in the same panel to highlight that the two regularization schemes have qualitatively the same effect on the final reconstruction error: that is, the quality of the inferred model at the optimal value is the same for both regularizations. Panels (b)-(c) show instead the optimal value of À computed either by minimizing the reconstruction error (shown with points) or by maximizing the test log-likelihood (diamonds). Panel (b) is actually repetition of Fig. 4-(d) and refers to the Ly-prior, while (c) refers to the L1 spectral prior. Again, we can observe a similar behavior of the two norms, with the only difference that Laa X independently on the chosen criterion. Finally, all the optimal values go to 0 when P + 00, as expected.
What is the effect of the regularization on the training dynamics? Considering that the standard training has non-monotonic behavior w.r.t. the training time, we would expect that, since the regularization strongly improves on the models' quality w.r.t. the standard case (at least at the optimal optimizing regularization strength Aopt). such a non-monotone behavior is diminished. This is indeed the case, as shown by Fig. 17-(a), displaying different training curves for different regularization strengths: the closest the regularization to its optimal value (highlighted in red), the smoother the model's quality is w.r.t. training time. At the optimal point the model's 's quality completely non-monotonic and approaches the fixed point the same reconstruction error the minimum w.r.t. time.
Actually, due to the simplicity of the GEBM it is even possible interpret the L2 regularization shrinkage correction protocol. Consider indeed the training fixed point given by Eq. (50). We stress again that the maximum- maximum-posterior matrix J has the same basis decomposition as the empirical covariance matrix, because the regularization term commutes with the other two terms in Eq. (50). By the dualism between covariance matrix and coupling matrix in the Gaussian EBM, we can think at the reciprocal values of Eq. (50) as eigenvalues a corrected covariance matrix w.r.t. CM, depending on X. We can therefore define another aggevalue-correeter covariance matrix, by using the analytic fixed point on the training dynamics obtained through the regularization:
C.M -al- La(
Σ
M MT u u
(54)
((00)) La (2)
By definition, the training fixed point obtained using i) training dynamics with the un-touched empirical covariance matrix EM plus the regularization term úi) regularization-free dynamics using matrix (54) are the same. Fig. 17-(a) shows indeed how the regularization modifies the eigenvalues of CM when interpreting the fixed point of the training dynamics
25

---

A thcoretical framework for overfitting in energy-based modeling
Model's quality
Corrected Eigenvalues
10
Training with GM
10
100
(b)
GM +2g(A
6 10
10-1
10
10
Eigenvalues of CM
x 10
10-1
10--
Gu
2 x 10
Gu
x
10
10-2
10
102 10 10-1
10
10
Figure 17. Effect of different L2-norm regularization strengths on the GEBM's learning dynamics. Panel (a) shows the reconstruction error vs time. The dotted blue line corresponds to the standard training over C.M All the other full lines correspond to training with a certain value of regularization strength obtained by numerically solving Eq. (50) for all modes. The regularization strength increases from yellowish colors to blueish (see the colorbar at the right). The curve corresponding to the optimal regularization that minimizes the reconstruction error after training is highlighted red. Panel (b): plot of the equivalent covariances modes corrected by the regularization. For each curve, scatter plot these values Eq. (54) against the population eigenvalues. Here we set p = 1.66.
as a shrinkage correction, Each set of points shows the quantities 1/J(0)) (A) (i.e. the eigenvalues of (54) ) vs the population ones c It intuitive to notice that the optimal regularization (red points) is the one that makes such corrected eigenvalues close as possible to the population ones. This entire reasoning holds analogously with the spectral Lı norm.
H.2. Empirical shrinkage correction through modes fitting
simple way to perform a heuristic shrinkage correction is to down-sample the empirical covariance matrix and estimate the asymptotic eigenvalues through fitting procedure. Starting from the available dataset with M samples whose covariance matrix is CM we randomly extract subsets of N < -M samples and estimate the eigenvalues of the size-reduced covariance matrices CrM. Every time a down-sampling procedure of this kind is performed, both the eigenvalues and the cigenvectors will be different from the original CM,, however, we here suppose to account for the eigenvalues, keeping the basis fixed to the one of Č.M After applying this computation to different values of N < m < M, for each eigenvalue o we can fit the resulting data {e M) according to
è (m)
m" A + BA
(55)
with v being an exponent of choice. The coefficients B will represent the asymptotic estimate of the a-th eigenvalue of the population matrix, c, corresponding the + oo extrapolation. After fitting each mode separately, we can construct an eigrrvalue-corrected covariance matrix as
Ō = Σ
B
M MT u
u
(56)
This procedure can in principle be generalized e.g by using a combination of powers in the fitting function, although in this work we only restricted to the functional form (56); secondly, the estimation of the fitting coefficients can be improved by collecting mean values of eigenvalues { } () by evaluating the -ampledd covariance matrix Ot multiple times with differrntt subsets of the M original data. In the experiments presented in the main text for the GEBM (orange line in Fig. 5-(b),, used simple linear fitting in 1 1/m (v 1) and 10 random resampling for each value of m, from which the mean eigenvalue is extracted to perform the fit. This linear scaling (in 1 //m) for the finite-size fluctuation of the eigenvalues in the GEBM is also justified by theoretical evidence (see e.g. (Ledoit & Péché, 201 1)). In the experiments for the Ising-BM BM instead, (orange line in Fig. 6-(c)) we find that best reconstruction is achieved with v = 1/2, while a linear fitting has always very bad performances. Also in this case we used 10 resampling steps. An example of such a fitting procedure is shown in Figure 18 for both the GEBM (in (a)) and for the BM (in (b)), each case for a given value of M.
26

---

thcoretical framework for overfitting in energy-based modeling
(a)
1/M
1.75-
2.0-
*
AM
B,
1.50-
23F
1.5
1.66
1.25-
1.00
1.0
18
0.30
0.5
0.25
p=7.81
1/M
0.0
0.00-
0.000
0,0012
0.0014
0.006
0.0088
0.010
0.02
0.04
0.06
0.88
0.10
0.12
1/m
1/ym
Figure 18. Examples of eigenmode fitting procedures for the GEBM and the Ising-BM. Each panel illustrates the procedure used to fit covariance matrix CM by downsampling to m M samples, order to extrapolate their behavior as m -> co. following Eq. (55). Panel (a) refers to the GEBM used in the main text (e.g.. Fig. 5). For fixed value of M such that p M/N = 1.66. we show a subset of eigenvalues {2, (denoted by markers), along with their downsampled counterparts {c(" for several values of m M (small circles). These downsampled eigenvalues are obtained by randomly selecting samples from the full dataset and computing the eigenvalues of the resulting covariance matrix, averaged over 10 independent instances. Dashed black lines correspond to fits using Eq. (55) with 1, and colored diamonds indicate the extrapolated intercepts B. i.c., the estimated eigenvalues Crosses mark the population eigenvalues {<a], showing that the fitted extrapolations are significantly closer to the true population spectrum than the empirical eigenvalues obtained from M samples. Panel (b) shows the analogous procedure for the Ising-BM, using = 7.81. In this case, the best fits are obtained with v = 1/2, and the horizontal axis is accordingly rescaled as 1 1//m.
eigenmodes
+ co.
Eigendecomposiioon of training dynamics on Boltzmann Machine
We consider Ising-liit like Boltzmann Machine for the inference of binary-valued data. The probability of configuration where E [-1,1] at given parameters is expressed as
1 7 e Jiyar sch.
p (w J.
(57)
We suppose to generate equilibrium configurations from known model with 0" (J", h²) (eventually these parameters are rescaled by external factor /3 that plays the role of an external inverse temperature), and we want to infer back the original model through likelihood maximization procedure. The LL of certain set of parameters O = =J,,h) is given by
M
Co (J) = M logg (x" J,h. Σ ED [xyæ5] +
Σ h.Ep [x:] log Z
(58)
where Ep !:] denotes the average w.r.t. the dataset D = {xx, | In what follows, for the analytic treatment of the training dynamics the ML procedure we neglect the problem of learning the external fields h;. This assumption consistent with a scenario where the data have null magnetizations. Therefore, the gradient the LL w.r.. the couplings J4p reads
OLp OJ,
ED [xv0z] - Ey [mmay] GM Es [2,2,)
(59)
Then, the couplings are updated as
Joy(+ J1(())))(())))) @Cv(t)
(60)
TJ
where y is the learning rate. From here, we can assume an ideal training with an infinitesimal learning rate to recast the evolution equation the matrix in time in the following matrix form
dJi
0LD
dJ
C. (xxx ))
(61)
dt
0Jgy J()
dt
27

---

A thcoretical framework for overfitting in energy-based modeling
where Ey [ denotes the average w.r.t. to model (57), and 7 1/7- Note that, since we assumed to neglect local magnetizations, the r.h.s (61) computed in its diagonal entries is 0. This is consistent with the fact that self-couplings Ja do not evolve in time, because they correspond to constant energy terms in the energy function. In order to make the models' correlation E, [x=²"] analytically treatable, we implement now a meenn-field approximation. We can exploit the following exact expression (self-consistent) for the correlator (Suzuki& & Kubo, 1968), which we further expand for high temperatures
Ey [mnarg]]= 8ig + (1-
- 8ez) B,
T2 Lanh Σ Jaatk Rs dig + (1 - -ds)) Σ JuHEs [Igaa]),
(62)
Then, in matrix form (Cij = EE (1) we get
CN + JC diag [JC]
(63)
where IN is the identity matrix of size N. and the operator diag [M] extracts the diagonal part of a matrix M. The last term is introduced to correct the wrong estimation of the diagonal entries of C, which should be equal to 1. The problem is that the above equation does not admit simple analytical solution for the model's correlation matrix C, which was the original goal. Instead, is typical implemented in the literature is the following expression of the linear response correlations
C- (JJ = (IN - J)
(64)
which gives a reliable estimate of the correlation matrix and is at the core of well-studied mean- mean-field like expression for the inferred couplings (Kappen & Rodríguez, 1998; Ficci-Tersengh, 2012). However, the diagonal entries of (64) are not in general equal to 1. Normally, this is not an issue because one is interested in correlations for * j (i.e. off-diagonal entries). However, in our approach such diagonal mismatch creates non-null gradient on the diagonal entries of J. Indeed, by plugging Eq. (64) into the LL's gradient, we get
dJ C (IN J)-1
(65)
dt
Now, the diagonal part of the r.h.s. of Eq. (65) not null anymore. In order to circumvent this additional issue, a possible solution could be to modify the gradient by remove the diagonal terms which would results in an nonphysical evolution of the self-couplings - "by hand":
dJ = C -0 - J) + diag [TN - (IN - -J)
(IN--]
(66)
dt
The last term in the above expression correctly fixes the diagonal problem for the matrix J in the gradient. This strategy is similar to what carried out in (Fanthomme et al., 2022) where the authors impose add spherical constraint on the gradient in the form of Lagrange multiplier. However, this leads to a complicated expression even for the training fixed point because the evolution of all the eigenvalues is now coupled. For this reason, since here we are interested in the dynamics of training, adding this constraint would result into a system of coupled differential equations for the eigenvalues, which has computationally the same complexity of the original problem, so there would be no gain in that. The simplest strategy is therefore to use the approximate expression for the correlator and avoid adding the constraint, so to use the gradient (65) as it is. Although it might seem crude approximation, it still allows us to decompose the dynamics in the same way as for the GEBM. Before going on, it is worth noticing that the diagonal matching problem is at the core of some refined mean-field like approximations for binary (lsing-like)
rikimum-enttrop models which exploit e.g. iterative diagonal consistency tricks (see e.g (Yasuda& & Tanaka, 2013; Kiwata, 2014) ). Therefore, as explained in Sec. 4 for the GEBM we can project the log-likelihood's gradient onto the eigenbasis of J, leading to an expression for the rotation of eigenvectors (same one as in (5)) and another one for the evolution of its eigenvalues, which is given below:
dJa
MM Caa
1
(67)
dt
IN Ja
which is the same equation shown in the main text (Eq. (11)). As explained in the main text, the solution of the above equation describes an independent evolution of eigenvalues which is not quantitative accurate with respect to the numerical results, but still captures the qualitative trend: in particular, it perfectly describes the separation of timescales in terms of PCA's modes during the training dynamics, an effect which is observed also in the numerical results (see panels (b) in Fig. 6).
28

---

thcoretical framework for overfitting in energy-based modeling
(a)
Ec(CMM C ())
10a
Ec(Caet C
10a
10*
105
10
10"
10
10
100
10a
10
t
t
Figure Supplementary results on the Boltzmann Machine for the inverse Ising problem. The model, dataset and training setting are the same in Fig. 6. (a): show the error between the covariance matrix of generated configurations from the model (along the training trajectory) and the covariance matrix of the training Ö., so E = O (t) (b): plot of the generation error computed w.r.t test again between the covariance matrices, i.c. Exee = || Ct. Both quantities plotted versus training time (number of updates) for different values of M shown in the colorbar. The learning rate set 7 10
10
10"
10
10
p
E(J", J(O))
10-2.
10
10
10a
10a
10
updates
Figure 20. Supplementary results on the Boltzmann Machine for the inverse Ising problem. The model, dataset and training setting are the same in Fig. 6, except for the system size which here is equal lattice size of L 32, that N I.2 1024. We show the Frobenius norm of the error between the true model and the trained one versus training time (number of updates) for different values of p M/N shown in the colorbar. The learning rate is 7 10-2
29


Reviews:
[{"review_id": "9kkBprV3ZP", "reviewer": "Reviewer_SJPf", "summary": "The authors present an analysis of training dynamics and overfitting in different settings (infinite data, limited data, continuous domain, binary domain) for a specific class of EBMs. The basic idea is to project the training dynamics to the principal components of the coupling matrix, which allows (in the class of models the authors study) for strong analytic claims.The authors also study common methods for mitigating overfitting (like regularization) in their framework.\n\nUPDATE AFTER REBUTTAL:\n\nThe paper is limited in that they are studying GEBM, which limits its immediate applicability, but the work is nonetheless interesting, novel and thoroughly conducted. The authors also addressed my concerns. I therefore change my rating to 'accept'.", "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": "## **1. Common reply about model relevance and applicability to more complex EBMs**\nWe thank all the reviewers for their comments. \n\nAll Reviewers have raised important concerns regarding the generality of the results presented in the paper, particularly their applicability to more complex or realistic generative models. While we fully acknowledge that our analysis is far from the architectures typically employed in contemporary machine learning, we would like to emphasize that, to the best of our knowledge, this is the first work  to tackle the problem of overfitting in energy based generative models. To start with something concrete and tractable we consider here Gaussian multivariate models which might indeed look over-simplistic and irrelevant regarding modern ML. The logic here is a bit similar as the one considering linear regression for understanding deep learning [Belkin2018,Hastie2022]. Actually we found a concrete parallel to that and **plan to add a new section in the main-text and supplementary material discussing it in detail**. The argument goes as follows: if we consider score-matching [Hyvarinen2005] as a proxy for studying theoretically the learning of EBMs, we arrive at a description of the learning dynamics in terms of a neural tangent kernel dynamics [Jacot2018] of the score function ($\\psi(x,\\theta) =-\\nabla_x E(x,\\theta)$, i.e. the gradient of the energy function of the EBM) of the form\n$$\n\\frac{d\\psi(x\\vert\\theta_t)}{dt} = -\\hat{\\mathbb E}_{x'}\\Bigl[K_t(x,x')\\psi(x'\\vert\\theta_t)\\Bigr] + \\hat \\phi_t(x)\n$$\nwhere the kernel $K$ and the source $\\hat\\phi$ is build on a tangent space corresponding to the derivative of the score function w.r.t the parameters. Considering then the dynamics corresponding to a lazy training regime [Chizat2019] where the NTK (and also $\\hat\\phi$) can be assumed deterministic we end up with a similar dynamics as for Gaussian EBM, driven now by the empirical covariance matrix of score features instead of the input features. The empirical covariance matrix on the input feature is now replaced by a covariance of tangent score functions, which are as well expected to be in the random matrix regime and overfitting will as well occur when weak modes of this covariance starts to be learned. So in the end we have potentially a much broader theory, which encompasses exponential models and more generally the EBM in the kernel regime.\n\n\n[Hastie2022] Hastie, Montanari, Rosset, Tibshirani Ann. of Stat. (2022)\n\n[Hyvarinen2005] Hyvärinen, Dayan,  JMLR (2005)\n\n[Jacot2018] Jacot, Gabriel, Hongler, NeurIPS (2018)\n\n[Shizat2019] Chizat, Oyallon, Bach, Neurips(2019)\n\n## **2. Specific answer to the reviewer:**\n\n**Methods and evaluation criteria:** \nThe experiments we have conducted on the GEBM, although performed at relatively small system sizes — with most results in the main text corresponding to $N = 100$ — exhibit excellent agreement with the asymptotic predictions derived from Random Matrix Theory (RMT). As shown in Fig.~4 (notably panels a and c), the empirical results obtained at $N = 100$ already display perfect overlap with the theoretical curves expected in the limit $N, M \\to \\infty$.\n\nRegarding the Boltzmann Machine (BM), it is true that the system size employed ($N = 64$ spins) is smaller than what is typically used in standard machine learning applications. Nevertheless, the qualitative behavior we observe is robust with respect to system size. These BM experiments are intended to demonstrate that the key phenomenology observed in the GEBM carries over to the BM case. In particular, we highlight the emergence of negative eigenvalues in the coupling matrix and the role of finite-sampling noise, which primarily affects the learning dynamics associated with the smallest principal components of the data. To strengthen the robustness of these findings, we plan to perform additional numerical experiments at larger system sizes for the final version of the manuscript.\n\n**Using real datasets --** Our results do not depend on the specific dataset used, as they are based on a given arbitrary covariance matrix, which can be derived from any dataset. The use of synthetic data is motivated by the need to control the asymptotic limit required for the Random Matrix Theory (RMT) analysis, and to ensure that the eigenvalues are not too small at this limit — a condition that would otherwise significantly increase the computational cost, since small eigenvalues necessitate extremely fine discretizations for solving the numerical RMT equations. Nevertheless, we can always consider a cutoff which should only have effects far away from the early stopping point.\nIn the final version of the paper, we will include an additional section in the Appendix where we reproduce the results using  different eigenvalue spectra taken from real datasets and different parameters for the synthetic model, explicitly showing that the qualitative behavior remains unchanged.", "claims_and_evidence": "The claims are very well supported: The authors base them on theoretical analysis and show excellent agreement with experiment. The projection of the training dynamics into eigenspace is well done and interesting.", "methods_and_evaluation_criteria": "The methods and evaluation criteria show quite exactly what the authors want to show. However, in order for the paper to be relevant for the wider field the evaluations are not sufficient: The experiments are small scale and the data very far from what ML methods are typically applied to. It is not clear how this analysis would work in more realistic scenario.", "theoretical_claims": "I did not do the math but it seems to be relatively straight-forward and the results are sensible.", "experimental_designs_or_analyses": "I am quite confident that the experiments are valid. They coincide with the theoretical predictions and the setting is relatively small-scale and controlled.", "relation_to_broader_scientific_literature": "The relate to random matrix theory and inverse Ising problems. The paper is well-embedded in this field.", "essential_references_not_discussed": "Not that I know of, but some of the topics touched have a vast literature spanning several decades, so I am not entirely sure.", "other_strengths_and_weaknesses": "Strengths:\n\n1) Setup & Execution: The paper is beautifully written and the analysis is clean and done well, the visualization and experiments are convincing\n2) Novel results: The results are non-trivial and the agreement with experiment is very good\n3) Interesting to some groups: The paper is very interesting to anyone working in the field of statistical methods applied to ML (and related fields like inverse Ising/Potts etc.). However (read below under \"weaknesses\"), these groups do not correspond to the core audience of ICML.\n\n\nWeakness:\n\n1) Model relevance: The type of models analyzed are toy models at best. I am not aware of anyone using these kind of models in moden ML (please correct me if I'm wrong). The most relevant case I know of where this analysis might apply are inverse methods applied to protein sequences, but even that seems to have been largely abandoned. The applied papers that the authors cite are either old or quite niche.\n2) Data relevance: The authors use either synthetic data or MNIST.\n\n\nI think the paper is worth publishing if the authors argue that their analysis could be extended to a more realistic setting. I am quite confident that if pursued consistently this type of analysis could be valuable for the wider community.", "other_comments_or_suggestions": "Very minor: It's \"Frobenius\", without the umlaut.", "questions_for_authors": "Could you argue that this type of analysis could be extended to more realistic settings? For now my rating is \"Weak Accept\" but I am happy to increase the rating if you can argue 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": "xfDbd6uITM", "reviewer": "Reviewer_MN9R", "summary": "This paper proposes a theoretical framework for analyzing overfitting in energy-based models (EBMs). The framework is built upon two special cases of EBMs (namely, Gaussian EBM and Boltzmann Machine for inverse Ising model), which admit analytical (or partially/asymptotically analytical) solutions for the learning dynamics and stationary points (optimal solutions). This framework is later used to analyze the driving factors behind overfitting, suggest new overfitting mitigation strategies and reinterpret existing approaches.\n\nIt is claimed that the main reason behind the overfitting in the EGMs considered is the interplay between initialization and different learning timescales associated with eigendecomposition of the weight matrix. It is also claimed that the analysis provided is relevant to more complex models.\n\n## update after rebuttal\n\nDuring the rebuttal, most of my questions have been addressed.\nThe remaining are:\n- Minor experimental analysis concern regarding the confidence intervals.\n- Ablation on spectra lacks details (e.g., the expression used to define the spectra for Figure R2).\n- I am still not convinced that one should mimic the spectral properties of the empirical datasets considered in the paper, because this data is of non-linear structure, and covariance matrix is a poor statistic to describe the properties of such datasets.\n\nAlthough complex NN-based models are not directly targeted, the analysis is thorough and rigorous, which is a decent start. I have increased my score to 4. For more information, please refer to my last rebuttal reply.", "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": "**Methods and evaluation**\n\n  1. It is true that the choice of the generation quality measure with the Frobenius norm does not reflect a measure of distance between two probability distributions in general, but it is inspired by real experiments where actually such a metric is widely used, for instance in inverse Ising/Potts approaches in order to compare e.g. correlations matrices between the original dataset and generated configurations. In the GEBM it is also true that the metrics suggested by the Reviewer is a better measure of discrepancy between distributions. We have computed in particular the Wasserstein distance between the true GEBM and the inferred one along the training dynamics: also this quantity turns out to display a non-monotonic behavior w.r.t. the training time, and moreover the optimal time computed using this metric seems to be the closest one to the optimal time obtained by maximizing the test-LL. We will discuss these other measures in the final version of the manuscript.  \n  2. **Choice of the spectrum.** For the GEBM analysis, we do not employ real datasets directly. Instead, we construct synthetic covariance matrices designed to mimic the spectral properties typically observed in empirical data, with the aim of explaining early stopping effects reported in real experiments. Specifically, we chose a synthetic spectrum with two distinct branches to visually and analytically distinguish between dominant (strong) and subdominant (weak) modes in the population covariance matrix. We would like to emphasize that this two regimes spectral structure is supported by previous works identifying similar two-regime behavior in real datasets (see, e.g., [Yang, R. et al., ICML 2022)]. In the final version of the manuscript we will display the real and synthetic eigenvalue spectra in log-log scale where this common trait can be more easily visualized. That said, the theoretical framework we develop does not rely on the specific details of the eigenvalue spectrum, provided the population covariance matrix is non-degenerate. We tested a variety of synthetic spectra and observed no qualitative differences in the results. For example, modifying the parameters in Eq.~(10) leads to the same learning dynamics as discussed in the manuscript.\n  The choice of synthetic spectrum is guided by numerical constraints: the RMT equations require discretization of the spectral density, and very small eigenvalues significantly increase computational cost due to the need for extremely fine resolution in the integration.\n  We understand the concern and will include a dedicated appendix in the revised version of the manuscript to explicitly illustrate the robustness of our results with respect to changes in the eigenvalue spectrum.\n\n**Theoretical claims**\nThe presence/absence of the term $\\Lambda_{ij}$ arises whether one assumes that the perturbation of the log-likelihood w.r.t. to one interaction $J_{ij}$ to be symmetric or not. Of course this holds only when $i\\neq j$, that is why the term $\\Lambda_{ij}$ appears. \nIt is actually true that in the following steps of the derivation we did not use this additional term when projecting the gradient, nor in the numerical gradient ascent procedure: this is basically equivalent to assume that perturbations instead are non-symmetric ( a more detailed discussion about the two different ways of taking likelihood's derivatives is given e.g. in [Magnus, J. R., Neudecker, H. (1999)]). A practical solution to this problem is to absorb it in the learning rate, meaning that the diagonal terms $J_{ij}$ will evolve in time with a learning rate doubled w.r.t. to the off-diagonal terms to compensate that factor.  \n\n**Experimental Designs Or Analyses**\n\n  1. About the mean values on Figures 1,7: both figures show the mean eigenvalue spectra of the empirical covariance matrix of a dataset at a given number of samples $M$: here the average is intended w.r.t. a certain number $n$ of different downsampling of $M$ samples from the original full dataset with $M^*$ samples. We have chosen a large enough $n$ (specifically $n=1000$) such that the standard error of any mean eigenvalue is not appreciable in the plot's scale. Having said this, we will add errorbars or shadows to the figure to highlight the standard deviation between realizations, reduce $n$ in order to have disjoint subset of samples to compute the covariance matrix.\n  2. This is correct. In the left panel, we wish to compare the theoretical behavior (continuous time limit) with the one obtain by maximizing the likelihood using GD. The continuous limit is correct for very small learning rate, which is why $\\gamma$ is very small. On the right panel, we only consider the resolution of the evolution equations, hence the learning rate didn't matter here as it only renormalize the time.\n\n**Relevance for more complex non-linear models** See the common answer written in the rebuttal of rev. SJPf.", "claims_and_evidence": "Overall, I find the claims related to GEBMs and BMs for the inverse Ising model to be well-supported by the theoretical and empirical evidence. The paper provided a comprehensive analysis of these models and corresponding overfitting mitigation strategies.\n\nThe only major claim I see problematic is the connection to non-toy-ish EBMs. Although the manuscript can provide some intuition regarding the general EBM case, the whole analysis revolves around spectral properties of the coupling/correlation matrices. The same goes for the protocols to mitigate overfitting. It is unclear how to extend these excellent results to a general, non-linear case.", "methods_and_evaluation_criteria": "I am mostly satisfied with the proposed method and evaluation criteria. However, there are several concerns which I would like to raise.\n\n1. On lines 262-263 col.2 and in Figure 3, $\\mathcal{E}_C = \\Vert \\hat{C} - C\\Vert_F$ is referred to as \"generation quality\". However, to my knowledge, this quantity is not connected to any of the widely accepted divergencies used to assess the generation quality: neither $f$-divergencies [3a], nor Wasserstein distances [3b]. Perhaps, the authors should refer to $ \\mathcal{E}_C $ as just \"covariance matrix reconstruction error\". I also suggest adding the closed-form expressions from [3a,3b] to track how well the learned $J$ reproduces the original distribution. This will enable overfitting analysis from the distribution matching perspective, which is more interesting in the context of generative learning. These additional metrics might also provide new information, thus complementing the analysis based on log-likelihood, energies and matrix errors.\n2. I feel uncertain about using MNIST, CIFAR10, HGD and Ising 2D datasets for correlation matrix spectral analysis. These datasets are of highly non-linear structure, which makes me question the relevance of the spectra reported in Figures 1,7 for the research provided. It is unclear why one should try to reproduce the spectral properties of covariance matrices of datasets which are far from being Gaussian.\n3. Still, taking Figures 1,7 into consideration, it is unclear why it was decided to introduce nonsmoothness at $x=1$ in (10). From FIgures 1,7 it is clear that real spectrum is smooth. Perhaps, the authors should provide additional experiments showing that the results (at least, the key results) are robust to the choice of the spectral density.\n\n[3a] Frank Nielsen, Kazuki Okamura. \"On the $f$-divergences between densities of a multivariate location or scale family\". arXiv:2204.10952\n\n[3b] Salmona et al. \"Gromov-Wasserstein Distances between Gaussian Distributions\". Journal of Applied Probability, 2022, 59 (4). hal-03197398v2", "theoretical_claims": "I satisfied with the theoretical part of the work. I find the corresponding theoretical claims convincing and backed up with not only proofs, but also the experimental results. Below I provide my (mostly, minor) concerns regarding this part.\n\nI can not follow the derivation on lines 93-109 col.2. Particularly, I am unable to reproduce the gradient provided in (3). My own calculations yield $ [-C^M_{ij} + (J^{-1})_{ij}] $. Additionally, if we assume $C^M=0$ (which, of course, does not correspond to any real case, but can be used for the sake of analysis), we should get $\\partial L / \\partial J{i j} =  \\partial \\log\\det J / \\partial J{ij} = (J^{-1}){ij}$, which does not hold when using (3). Finally, in (4) the $\\Lambda$ matrix seems to disappear. Please, clarify.", "experimental_designs_or_analyses": "I have minor concerns regarding the experimental setups and analyses.\n\n1. Perhaps, using several seeds and reporting mean values and error bounds in Figures 1,7 will be more convincing to the latter analysis (it is hard to verify the under-/overestimation observations judging by several samples only).\n2. From the text, I understand that Figure 2 is provided to highlight the difference of $J_\\alpha(t)$ dynamics when (a) ground-truth or (b) empirical covariance matrix is used. However, it is therefore unclear why different learning rate $\\gamma$ was used.", "relation_to_broader_scientific_literature": "I can not pinpoint any specific prior work that targets the problem of overfitting in EBMs.", "essential_references_not_discussed": "I can not recall any specific work which is essential for understanding the context of the submission and is not cited.", "other_strengths_and_weaknesses": "The paper is comprehensive, well-written and was enjoyable to read.", "other_comments_or_suggestions": "I understand that the topic is deeply rooted in statistical mechanics. I also acknowledge the freedom of choosing the notation which is more convenient for the authors. Still, I find the latter to be a little confusing.\n\n1. From my experience, in papers on machine learning or statistics, $\\langle \\dots \\rangle$ is rarely used to denote expectation; rather, it can be used for empirical averaging, whereas the expectation is denoted by $\\mathbb{E}$.\n2. Additionally, denoting the ground-truth covariance matrix by $\\hat{C}$ and the empirical estimate by $C^M$ is also very unusual. It is usually expected that `\\hat` is used for estimates, and no accent marks are used for the ground-truth values.\n\nTypos:\n1. Lines 57, 60 col.1: inconsistent usage of spaces before and after \"—\".\n1. Similar for lines 636-638.\n1. Lines 74-75 col.2: \"-\" is used instead of \"—\".\n1. Line 101 col.2: \"the gradient of (2)\" - perhaps, \"the gradient in (2)\" (as it is the gradient of $\\mathcal{L}$) was intended.", "questions_for_authors": "1. How the framework proposed can be extended to complex EBMs?\n2. Have you tried other spectral densities instead of (10)?", "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": "kaTQEg14kE", "reviewer": "Reviewer_JaLg", "summary": "The given paper provides a theoretical analyses of overfitting in Energy-based models. Particularly, the scope of the paper is restricted to the analyses of Gaussian Energy-based Model (GEBM), wherein the authors show that the maximum likelihood (ML) training dynamics of GEBM is decomposed into different timescales. Specifically, the dominant mode features (corresponding to higher eigenvalue) are learnt early on whereas non-dominant mode features (corresponding to lower eigenvalue) are learnt later on during training. \n\nThe provided analyses shows that for finite number of training samples, the test log-likelihood (LL) improves upto certain time and then starts deteriorating. The authors provide analytical expression for this optimal stopping time ($t_{opt}$) using methods from RMT. Furthermore, the authors also provide several other results on optimal scaling factor, etc using RMT for a regularized training of GEBM. The authors conclude that overfitting observed in GEBM is mainly caused by the variation/noise in non-dominant modes in finite training sample setting, while this problem can be solved exactly when one has access to true model parameters. To this end, the authors propose to overcome the phenomenon of overfitting by using RMT to predict true model parameter via asymptotic results. However, most of the proposed solutions rely on quantities that rely on true model parameters. Lastly, the authors propose to fit a linear model to learn the eigenvalues of the data covariance, then use the model to extrapolate and use the extrapolated eigenvalues to determine the quantities of interest and avoid overfitting in finite sample cases.", "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": "**Claims and evidence:**\nWe thank the Reviewer for the suggestion, we will put in the appendix additional plots highlighting the behavior in $m$ of the down-sampled eigenvalues for both the GEBM and the BM.\n\n**Methods and evaluation criteria:**\nThe eigenvalue spectra of the covariance matrices of real datasets as MNIST/CIFAR10 is depicted to guide us in a specific choice of a functional form of a a continuous spectral density used to derive the asymptotic results through RMT. The analysis on the GEBM can however be done using directly the eigenvalue spectrum of a real dataset. The only critical thing that should be considered is the addition of a threshold in the eigenvalues to avoid the existence of very small ones that would make the convergence extremely slow, and that have no effect on the timescales where the early stop point occurs. We will add some results on this regard in the appendix (see also last reply to Reviewer SJPf).\n\n\n**Experimental Designs Or Analyses:**\nWe would like to stress that the experiments on real datasets are not very meaningful here, as the only thing that matters is their covariance matrix, i.e. the highest-order sufficient statistics of both the GEBM and BM. Our setting is not meant to reproduce with any level of quality samples from real datasets (e.g. MNIST/CIFAR10 images). In other terms, the GEBM is not used here to learn a \"good\" EBM to generate good-looking samples: that would be basically impossible because the GEBM only encodes for the first two moment statistics of the data with a single multivariate Gaussian, so that any sample generated according to the learnt GEBM would be far from a realistic image in the datasets cited above. The same holds for the BM: the binary BM does perform well as generative model, in the sense that it generates good Ising model samples, but this model, without hidden nodes, is known to perform very badly with images even at the level of MNIST dataset, as discussed e.g. in [Decelle et al, SciPost Physics, 2024], because high-order interactions are important.\nThe reason why we mention real datasets in the paper is to justify our choice of the synthetic spectrum but this procedure only looks at the eigenvalue spectrum of the real datasets' covariance matrix. Still, as mentioned in the reply to Reviewer MN9R, even such a reduced information about the spectral structure of real datasets is important to determine generalization properties even in deep networks (Yang, Mao, Chaudhari, ICML 2022). \n\n**Strengths and weaknesses**\nIn the manuscript, we discuss both GEBMs and BMs. For the latter, the analytical treatment is restricted to the high-temperature phase and remains approximate; nonetheless, we demonstrate that it successfully captures the qualitative behavior observed in real experiments. In particular, we show that early stopping effects in BMs are analogous to those identified in GEBMs, and that they can be mitigated by applying the same regularization recipe introduced in the GEBM setting.\n\nThe models considered in this work serve as a controlled playground to analyze overfitting in a highly tractable setting. As discussed in our general response, we believe that the insights obtained here can be extended to more complex and realistic setups; however, such generalizations require dedicated studies, which we plan to pursue in future work. Our current goal is to establish a solid theoretical foundation upon which these future developments can be built.\n\nAbout point 2: as explained before it is not meaningful here to compare single generated samples as it would be with datasets of images, but we will compare e.g. the generation error (i.e. the error between the covariances of generated samples and the population one) or other distance measures sugggested by Reviewer MN9R in the regularized and the un-regularized case. As stressed in the paper, it is true that in the GEBM the most dominant contribution to the error in the coupling matrix comes from the weak PCA directions, which are the least dominant in measures of generation quality, so we do not expect an improvement as net as in Fig.5(b).\n\n\n**Questions For Authors**\nThis is correct: from Eq.1 the population covariance matrix is defined as the inverse of the GEBM's coupling matrix: that would also correspond to the empirical covariance matrix obtained computed with an infinite set of independent samples from Eq.1.", "claims_and_evidence": "The authors have provided empirical and experimental proofs for most of their claims. However, I would like to see the following results as well:\n\n1. Can you provide a plots of $\\{c_\\alpha^m\\}$ against different $m$ so that one can verify that a linear fit would be good for extrapolation?\n2. A similar plot for the Boltzman Machine learning would be helpful.", "methods_and_evaluation_criteria": "**Method**: The authors haven't provided any new method as such. Rather they give a rigorous theoretical insights into overfitting of EBMs (although this is limited to GEBM).\n\n**Evaluation**: I understand that the provided analyses is limited to GEBM. However, the authors have shown spectral density on several datasets such as MNIST and CIFAR10. In that case, can the authors comment on how does their analyses apply to these datasets?", "theoretical_claims": "I checked the theoretical correctness of the provided proofs and I cannot find any obvious mistake. However, I might have missed any tricky mistake as I am not well versed with RMT.", "experimental_designs_or_analyses": "1. The authors show spectral density of complex datasets like MNIST and CIFAR10, however, they don't show visual results of sampled datapoints obtained after training.\n2. The authors should sample datapoints from MNIST and datasets mentioned in Fig. 7. Then compare how well these sampled datapoints are with few baselines like IGEBM (Du & Mordatch, 2019).\n3. One should use Negative log-likelihood (NLL) and FID for the above comparison.", "relation_to_broader_scientific_literature": "N/A", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "#### Strengths\n1. The paper is fairly well written and presented.\n2. To my knowledge, this is one of the first work that discusses the phenomenon of overfitting in EBMs. However, the scope of the provided analyses is limited to GEBM.\n3. The conclusion drawn from the paper seems fair: the training dynamics first learns the dominant modes (which should correspond to low-frequency components) and then learns non-dominant modes (which should correspond to high-frequency components). This is in line with empirical observations as well.\n\n#### Weaknesses\n1. The paper seems weak empirically/experiment-wise. Although the paper provides a 'theoretical framework', the claims need to be verified on real world datasets.\n2. The authors should show the result of sampling after training the model to verify the correctness of the method. E.g., a visual example of samples obtained using overfitted model against samples obtained using regularized model. \n3. Identification of few key quantities like - $\\lambda_{opt}$, $t_{opt}$ are not feasible since the true model parameters are unknown in practice.\n4. The scope of the proposed analyses is limited to GEBM. However, somehow the observations are consistent with non-GEBMs like BM and datasets like CIFAR10.", "other_comments_or_suggestions": "1. Line 40, second col: generative -> generative modelling\n2. Line 316, second col: circunvent -> circumvent\n3. Line 373, second col: ~~Ref~~ \n4. Line 396, second col: exit -> exist", "questions_for_authors": "1. $\\hat{C}$ has been defined two times - after Eq. 1 and after Eq. 7. Are both the definitions same due to LLN?", "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": "47l9HyB0uN", "reviewer": "Reviewer_BiKA", "summary": "This paper analyses the training dynamics of learning a multi-dimensional Gaussian distribution from data. The training dynamics considered here is a continuous time gradient ascent optimizing the maximum likelihood objective. Under some assumptions on the starting point of the learning dynamic, the authors use techniques from random matrix theory to study the effect of finite training samples on the learning. They use this analysis to study effect of regularization and early-stopping. The analysis is also extended to Ising type models with the learning done in the high temperature/mean-field setting.", "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": "**Methods and evaluation criteria**\n\nWe answer this comment about the choice of the spectrum in the answer to Rev. MN9R.\n\n**Relation To Broader Scientific Literature:** Our work is focusing on the full-batch case since it is the typical case for the BM since you only need to compute the covariance matrix once.. We will add a few comments on the effect of considering minibatches in our framework. The second point answered below.\n\n**Other Comments Or Suggestions**\nWe refer to the general reply (posted on Reviewer SPJf's section) about extensions to general EBMs. \nWe are unsure why the Reviewer suggests removing or altering the purpose of the discrete-variable analysis. The aim of this section is to demonstrate that the main features and insights obtained from the GEBM analysis—particularly the interpretation of the early stopping point—also apply to the case of discrete variables, within an otherwise identical setup. In both cases, the analysis is carried out on the standard gradient ascent algorithm for the log-likelihood maximization, which is the most generic way to train EBMs of arbitrary complexity. We do not claim that the proposed \"cleaning approach\" is sample-optimal or superior to other existing methods in the specific case of the inverse Ising problem; for instance, we are not discussing any comparison with other training schemes (e.g. pseudo-likelihood maximization or the interaction screening method), nor comparison with other mean-field like techniques suited for the Ising-BM. However, we will incorporate a discussion of the suggested references in the revised manuscript, and we would be glad to refer to a more recent review if the Reviewer could kindly provide a reference.\n\n**Questions for authors**\nAbout the alignment of $J(t)$ with $C^M$ at $t=0$: yes the dynamics remains aligned for all time in that case. If the initial condition $\\boldsymbol{J}(0) $ is aligned with $\\boldsymbol{C}^{M}$, it means that the two matrices commute, namely $\\left[\\boldsymbol{C}^{M}, \\boldsymbol{J}(0)\\right]=0$. From here, any gradient ascent step will add a term proportional to $-\\boldsymbol{C}^{M}+ \\boldsymbol{J}^{-1}(t)$: since both terms commute with which keeps the two matrices commute with each other. Alternatively, by looking at Eq. 5 (right) one immediately sees that if the initial condition is aligned with the eigenvector basis of $\\boldsymbol{C}^M$, then $c_{\\alpha, \\beta}=0$ and the two matrices remain aligned at all times.\nAdditional numerical details about the alignment of the eigenvectors at the first stage of the training are discussed in Appendix B; it could be interesting to have a characterization of this transient time, but we don't have it yet and should follow from the more general theory mentioned in the common answer, that we plan to only sketch in this paper.\n\n**About overfitting:** We thank the Reviewer for this question, it is indeed a point that might be worth discussing in the core of the paper. Our interpretation is the following: In the context of GEBM the weak modes of the empirical covariance that are eventually learned by the model correspond to directions of variance poorly estimated by the data and lead to overestimated coupling eigenvalues. As noticed by the Reviewer we are in the under-parameterized regime in this \ncase even though the ratio #samples/#parameters ($=M/N^2$) is typically below 1 in our experiments\nbut what matters is the ratio $\\rho = M/N$ which becomes critical when equal to one because the covariance matrix ceases then to be invertible at that point. Nevertheless, when approaching \nthe interpolation threshold from above ($\\rho>1$) there is a departure of the test LL from the train LL\ncorresponding to overfitting in this proportional scaling limit. Then \nwhen looking at the more general formulation mentioned above corresponding to a kernel regime of score matching, this appears as a special case where the regression factorizes into $N$ independent problems (once the coupling matrix is aligned with the covariance matrix), which leads to consider $M/N$ scaling instead as $=M/N^2$.  If we now consider the kernel setting with general EBM, there is no such factorization in general and we recover the standard picture of overfitting with $\\rho = M/P$ where $P$\nis the number of parameters. The over-parameterized regime then corresponds to learning the weak modes of the Gram matrix of the score features, i.e. typically high frequency modes  able to define a localized energy function on each sample point.\n\n**On the derivation of Eq. 5** We will add the details of the derivation in the appendix. The derivation first consider the decomposition $J_{ij} = \\sum_\\alpha v_i^\\alpha J_\\alpha v_j^\\alpha$ before projecting the gradient on this new basis. Then we identify the diagonal terms leading to the dynamics of the eigenvalues and the off-diagonal ones leading to rotations in the space of the eigenvectors.", "claims_and_evidence": "The claims made in the paper are not clearly supported by analysis. This does not give any theoretical framework to analyse the training of EBMs. The analysis only applies to gaussian fitting where the explicit training dynamic can be analyzed. The general case is much more involved with a non-convex optimization dynamic. Also, applying the proposed methodology to the visible Botlzmann machine case will just give results that are not very interesting. In that case, the mean-field approach fails in the low-temperature regime where interesting multi-modality of the model emerges.", "methods_and_evaluation_criteria": "The use of synthetic data is probelmatic. Fig 1, says that the synthetic model closely mimics the real data sets. This does not seem to be true. The range of eigenvalues is completely different. Also the $M \\rightarrow \\infty$ limit of the synthetic model shows non-smooth behavior. It is unclear, if the results of the synthetic data model actually closely mimic that of actual datasets. This can be addressed by numerical comparison of the training dynamics of both the cases.", "theoretical_claims": "(See questions)", "experimental_designs_or_analyses": "(see questions)", "relation_to_broader_scientific_literature": "The authors have not clarified how their work connects to existing literature in analyzing gradient descent for convex problems under stochastic noise. Also, for models on discrete variables there exists other works that show that the learning problem can be solved sample optimally using methods based on learning conditionals. (See questions)", "essential_references_not_discussed": "(see questions)", "other_strengths_and_weaknesses": "(see questions)", "other_comments_or_suggestions": "> The paper can be improved if the initial claims are toned down. As mentioned before, the analysis here does not extend to general EBMs. If the authors think that practical training protocols for generals EBMs can be studied using the gaussian approximation, then they should support that claim strongly with numerical or theoretical evidence in the paper.\n\n> The discrete model analysis has to be removed or reworked to a more useful approach analyzing algorithms like interaction screening [VMLC16] or logistic regression [WSD19]. The review by Nguyen et.al. [CZB17] cited in the paper is unfortunately outdated for the inverse ising problem and does not consider the substantial advances in developing efficient algorithms for this problem in the last 7-8 years. \n\n[VMLC16] Vuffray, Marc, et al. \"Interaction screening: Efficient and sample-optimal learning of Ising models.\" Advances in neural information processing systems 29 (2016).\n\n[WSD19] Wu, Shanshan, Sujay Sanghavi, and Alexandros G. Dimakis. \"Sparse logistic regression learns all discrete pairwise graphical models.\" Advances in Neural Information Processing Systems 32 (2019).\n\n[CZB17] Nguyen, H. Chau, Riccardo Zecchina, and Johannes Berg. \"Inverse statistical problems: from the inverse Ising problem to data science.\" Advances in Physics 66.3 (2017): 197-261.", "questions_for_authors": "> Could the authors explain here or add a note to the Supp. Mat. showing the intermediate steps on how Eq (5) is derived? There some clarity lacking here. For instance the LHS of the second equation should be for the inner product between d(v^\\beta)/dt and v^\\alpha.\n\n> Do the authors assume that J(t) and C^M are aligned for all t in the training dynamic? If these are aligned at t= 0, does this imply alignment at all times?\n\n> Is the overfitting reported in this case really a case of overfitting in the traditional sense? In my view, the hypothesis model is not over-parametrized here. The fact that the test metrics achieve their optimum at a different time than the training metric is just an artifact of the different noise instantiations in the test and train dataset. That is, they are literally different metrics whose difference decreases in infinite sample limit.  This does not seem to be an example of a complex model essentially memorizing the training samples and then failing to generalize on the test dataset.", "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 presents a theoretical framework for analyzing overfitting in energy-based models (EBMs). The framework is grounded in two specific cases—Gaussian EBMs and Boltzmann Machines applied to the inverse Ising model—which allow for analytical solutions to the learning dynamics and stationary points. This foundation is then used to investigate the underlying causes of overfitting, propose new mitigation strategies, and offer reinterpretations of existing approaches.

The paper was reviewed by four reviewers. After the rebuttal, it continues to receive mixed ratings. The primary point of disagreement centers on whether the paper retains value despite its limited generalizability and applicability.

Reviewer BiKA (Weak Reject) expressed concern that the analysis is limited to Gaussian EBMs, where explicit training dynamics can be analyzed. The reviewer notes that extending this to more general cases—characterized by non-convex optimization dynamics—is significantly more complex, which may limit the broader applicability of the findings. Reviewer JaLg (Weak Reject) acknowledged that the paper is among the first to explore overfitting in EBMs, but noted that the analysis is confined to GEBMs, further limiting its scope.

In contrast, the other two reviewers were more favorable. Reviewer MN9R (Accept) felt that most concerns were addressed during the rebuttal and described the paper as comprehensive, well-written, and enjoyable to read. Reviewer SJPf (Accept) agreed that the focus on GEBMs limits the immediate applicability of the results, but found the work to be novel, interesting, and thoroughly conducted. The reviewer also confirmed that the concerns were adequately addressed in the rebuttal.

The AC facilitated a discussion among the reviewers. Overall, while there is a shared concern about the limited applicability of the work—particularly that the EBMs studied are somewhat “toy” models and not well-suited for modeling complex, real-world data—there is also recognition of the paper’s strengths. From both theoretical and empirical perspectives, the submission is solid.

Reviewers MN9R and SJPf actively championed the paper during the discussion. Although the scope is limited to specific instances of EBMs, the AC believes the theoretical insights offered by the authors represent a meaningful contribution. Such theoretical advancements are as valuable as empirical breakthroughs.  The AC also finds that the rebuttal adequately addressed the technical concerns raised, and that the paper is well-written, clearly organized, and makes a worthwhile contribution to the field, thus recommending an acceptance. 

The AC encourages the authors to revise the paper by incorporating the reviewers' suggestions to further strengthen the work.
Decision: Accept (poster)