A Theoretical Framework For Overfitting In Energy-based Modeling

Methods and evaluation criteria

We answer this comment about the choice of the spectrum in the answer to Rev. MN9R.

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.

Other Comments Or Suggestions We refer to the general reply (posted on Reviewer SPJf's section) about extensions to general EBMs. We 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.

Questions for authors About 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. Additional 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.

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 case even though the ratio #samples/#parameters ($=M/N^2$) is typically below 1 in our experiments but 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 the interpolation threshold from above ($\rho>1$) there is a departure of the test LL from the train LL corresponding to overfitting in this proportional scaling limit. Then when 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$ is 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.

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: We 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.

Methods and evaluation criteria: The 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).

Experimental Designs Or Analyses: We 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. The 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).

Strengths and weaknesses In 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.

The 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.

About 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).

Questions For Authors This 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.

Methods and evaluation

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.
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. 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. 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.

Theoretical claims The 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. It 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.

Experimental Designs Or Analyses

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.
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.

Relevance for more complex non-linear models See the common answer written in the rebuttal of rev. SJPf.

1. Common reply about model relevance and applicability to more complex EBMs

We thank all the reviewers for their comments.

All 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

$$\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)$$

where 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.

[Hastie2022] Hastie, Montanari, Rosset, Tibshirani Ann. of Stat. (2022)

[Hyvarinen2005] Hyvärinen, Dayan, JMLR (2005)

[Jacot2018] Jacot, Gabriel, Hongler, NeurIPS (2018)

[Shizat2019] Chizat, Oyallon, Bach, Neurips(2019)

2. Specific answer to the reviewer:

Methods and evaluation criteria: The 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$.

Regarding 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.

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. In 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.