Cascade of phase transitions in the training of energy-based models

Detailed anwers to the questions:

1. The referee is right, we should have been more careful when defining the quantity. The distribution of the hidden nodes is indeed of zero mean and variance of order $1/N_v$. The value is indeed $N=N_v$. This is the correct scaling to have a large dimensional (large $N_v$) limit.
   Let us clarify better this point. In a general formulation, we can consider that the probability distribution of the hidden nodes follows a Gaussian distribution with a variance $\sigma_h^2$. In this case, the effective model on the visible nodes is given by the following Hamiltonian $H=-\sigma_h^2 ( \sum_i s_i w_i )^2$, where we assume that the weights are of order $\sim O(1)$. This model can be analyzed as a function of \sigma_h. In practice, if $\sigma_h^2$ scales as $N^\alpha$, where $\alpha<-1$, the system is dominated by its entropy and remains in the paramagnetic phase where the average of each variable $\langle s_i \rangle = 0$ since the interaction is too weak. If \alpha>-1 the interaction is very strong and the system is frozen in configurations where $\langle w_i s_i \rangle = \pm 1$. Therefore the correct scaling is when the weights are of order $\sim O(1)$ and the variance of the hidden $\sim O(1/N_v)$ in order to observe a transition (as the magnitude of the weights change) of a disordered phase without particular orientation, and an ordered phase where the variables are parallel to $w_i$. For practical experiments, it is not necessary to impose the scaling since the learning process bring the system into the interesting regime where information is stored in the model parameters, as illustrated by our experiments.

2. The sentence at lines 161-162, refers to the type of phase transition occurring in the system. In the case with multiple features, the first phase transition is very similar to the one with one feature, in the sense that it is of the same “nature”. In this regime, corresponding to the first phase transition, everything happens as if there were just one main direction, which is $\xi^1 + \xi^2$. We will correct this sentence to make it clearer in the final version.

3. The theory is written for the binary-Gaussian RBM in the main text and for the binary-binary RBM in the Appendix. In the latter case, we show that in order to observe an interesting behavior we need to have an extensive (large) number of hidden nodes. On this aspect the gap is not strong. The main difference between theoretical and experimental parts relies on the dataset which is much simpler in the former case. Based on the article [8], we expect that even if the dataset present various modes (small compared to the number of visible/hidden nods), the system can still be described in a mean-field setting and, consequently, the phenomenology observed for one or two modes can be extended to many modes (but only a small number). This is confirmed by our experiments. Now we do not claim that this description is valid all along the training process, and we can see that at large training time the mean-field theory description is not as accurate anymore. We will comment on this in the final version.

4. This sentence might be indeed unclear and thus will be rewritten in the final version. By this sentence we intended to clarify that as the system is learning the first direction to split the data, the dataset projected on the remaining direction is mostly Gaussian distributed. The reason is that in the other direction, the projection is essentially given by the sum of random i.i.d. weights (non sparse) multiplied by the spins $s_i$. Furthermore, we can provide in the Appendix, if the reviewer thinks it would be useful, a figure showing that the projections on these directions corresponds to random noise as long as the other directions are not learned.

We thank the reviewer for reading our paper and for their comments. We think that they are instrumental in improving the quality of our work. We think that the gap between the theoretical part and the experiments highlighted by the reviewer is partly due to a gap in the presentation (which we will revise accordingly) and not in the analysis. Indeed, for simplicity, we have chosen to present the binary-Gaussian RBM in the main text, but we do analyze the binary-binary RBM in the Appendix. There is a clear relationship between the projections and projected weights in the theoretical section and the singular value decomposition described in the empirical evaluation. We think that this has not been clarified in the current version but it can be easily explained. We will discuss this and make it explicit in the revised version. Apart from that, the discrepancy between the theoretical and the experimental part is due to simplifications in the theoretical part, e.g. all hidden nodes have the same weights, which was done to keep the analysis tractable.

We answer the question into details in the official comments while commenting them very briefly:

Questions:

1. The referee is right, we should have been more careful when defining the quantity. The distribution of the hidden nodes is indeed of zero mean and variance of order $1/N_v$. The value is indeed $N=N_v$. This is the correct scaling to have a large dimensional (large $N_v$) limit. We provide a detailed explanation in the official comment.
2. The sentence at lines 161-162, refers to the type of phase transition occurring in the system. In the case with multiple features, the first phase transition is very similar to the one with one feature, in the sense that it is of the same “nature”. In this regime, corresponding to the first phase transition, everything happens as if there were just one main direction, which is $\xi^1 + \xi^2$. We will correct this sentence to make it clearer in the final version.
3. The theory is written for the binary-Gaussian RBM in the main text and for the binary-binary RBM in the Appendix. In the latter case, we show that in order to observe an interesting behavior we need to have an extensive (large) number of hidden nodes. On this aspect the gap is not strong. The main difference between theoretical and experimental parts relies on the dataset which is much simpler in the former case. Based on the article [8], we expect that even if the dataset present various modes (small compared to the number of visible/hidden nods), the system can still be described in a mean-field setting and, consequently, the phenomenology observed for one or two modes can be extended to many modes (but only a small number). This is confirmed by our experiments. Now we do not claim that this description is valid all along the training process, and we can see that at large training time the mean-field theory description is not as accurate anymore. We will comment on this in the final version.
4. This sentence might be indeed unclear and thus will be rewritten in the final version. More details are provided in the official comments.
5. We thank the referee for pointing out that we should clarify further the panel E of Fig. 4. The aim of this figure was to show the “hysteresis” phenomenon which is a clear signature, devised and observed in statistical physics, to unveil that a high-dimensional probability measure had a phase transition at which it split in two different lumps. A complete explanation is provided in the official comments.
6. A detailed comparison between ref [8] and our work is provided in the official comments and will be better discussed in our final work. Concerning the relationship between sec 4 and 5, the reviewer’s comments helped us to understand that we have to clarify the relationship between the SVD decomposition used in Sec. 5 with the projections used in Sec. 4. There is a direct connection between the PCA and the preferred directions ($\xi^1 + \xi^2$ and $\xi^1 -\xi^2$) used in the theoretical analysis, but it was not clear enough in the current version. We put a detailed explanation again in the official comments and we will try to make it clearer in the final version of our work.

Minor:

1. In the Mattis model $\xi$ correspond to a preferred direction (pattern) of the model, which is a frozen variable. In the Mattis model, $\xi$ can be any fixed binary vector. We will add a section in the appendix to describe with more detail the statistics of the model, as well as the conditions for $\xi$.
2. We thank the referee for pointing out these issues with Fig 1, we will correct them and define properly the quantities.
3. The referee is right, we will correct the typographical error.
4. The learning rate comes in front of the r.h.s. of eq. 3. We will add it in the final version. The learning rate naturally impacts the learning behavior, the crucial ingredient is to keep the ratio learning_rate/N small.
5. We plan, as mentioned before, to add first a section describing the Mattis model, and then another one more specific explaining the BG-RBM in a more general context.
6. We agree with the referee, we will carefully read the whole work to avoid confusion in the definitions and correct the typographical errors.
7. We thank the referee, we will check for this problem in the final version.

Limitations: We will add a new section Limitations, in which we will discuss the limitations of the current analysis: (i) the extent to which these results are applicable to more complicated models, (ii) the limitations of the analytical setup, and (iii) the limitations of the numerical setup (problems with insufficient MCMC sampling).

Thank you for the authors' careful responses and explanations. After reading the responses, I have resolved most of my concerns regarding the paper. However, I remain uncertain about the readability of the manuscript, as I am unable to review the revised version, which is expected to differ significantly from the current draft. I understand that this limitation is due to the conference's policy, and while the authors may be disappointed, it makes it difficult for me to increase the presentation score. Nonetheless, I acknowledge the significance of this work and have updated my score accordingly (4 -> 6). I hope the authors will take care to clarify the notations and statements, and include the limitations of this work in the final version.

We thank the reviewer for reading our paper and for their comments. We do not entirely agree with the reviewer regarding the target audience of this article.

First, the RBM itself is a model introduced by Hinton and Sejnowski for the purposes of machine learning. This model has been used in ML for decades, it was used as a building block of deep neural networks for a while, and is still used for ML tasks today, notably for interpretability, and in particular in low-data regimes (e.g. in neuroscience, population genetics,...). . While it is true that our analysis relies on phenomena more commonly studied in physics, our work shows that understanding these phenomena is crucial to fully understand the training of RBM, and how the weights are shaped during training. To name just a few elements: the second-order phase transition implies a sudden increase in mixing time which directly affects the quality of the training of RBMs, knowing that RBMs learn the direction of PCA helps to monitor training and to understand the structure of the weights, and we show how the exponential growth of the eigenvalues of the weight matrix is determined by the variances of the clusters.

Questions:

• It is not clear how the transitions would change in general. One should expect the first singular values to still correctly match the PCA components, because at the beginning of the learning process the mixing time is very small, which means that the estimation of the negative term by the CD-1 scheme (only one step initialized in the examples in the minibatch) is roughly sufficient. This also means that the first transition should be very similar to that of the PCD-100 studied in this article. It is not yet clear what should happen after the transitions start to follow each other, as the mixing times significantly exceed the 1 step used by the algorithm. To test this, we repeated the analysis for training with MNIST using the CD-1 algorithm. We show the equivalent of Fig.3 A,B,C in Fig.1 of the attached PDF file. We find that the image of the cascade of transitions remains unchanged, but the overall behavior is different because the susceptibility increases sharply with the growth of the corresponding eigenvalue, but then decreases again. In Fig. 2 in the attached PDF file, we compare the first 5 singular vectors with the first 5 PCA components of the data set, which show extremely similar patterns (within an overall sign). This picture is very similar to what we observe when analyzing the HGD dataset, where the mixing times become so large that even PCD-100 is no longer able to thermalize properly. We will include a subsection in the paper where we briefly discuss the impact of MCMC estimation of the gradient on the transitions.

• Our analysis suggests that similar phenomena probably occur in more complex ML models. First, in Deep Boltzmann Machine, a previous work [Ichikawa, Y., & Hukushima, K. (2022). Statistical mechanical study of deep boltzmann machines given weight parameters after training by singular value decomposition. Journal of the Physical Society of Japan, 91(11), 114001] shows that the SVD of the weight matrix is crucially important as well. Second, in more general EBMs, it is also highly probable that such phase transitions occur and that the mixing time increases. Finally, the cascade of phase transitions that are observed here are very akin to the ones recently observed in the sample generation of diffusion models [23-25]. For the latter, it is conceivable that during the training dynamics the model of the score first learns the most prominent features, which needs less data (and hence less training time, at least in an online learning schedule) and then later the finer ones.

We thank the reviewer for carefully reading the manuscript and for their interest in it.

We agree with the reviewer that the next natural update of this project will be to investigate these features in more complex generative models. Furthermore, we still believe that it is important to establish secure foundations to build on, and we feel that this work is already very dense and compact to include further analysis on a different type of models, which would require different techniques and settings.

Nevertheless, we believe that the applicability of these results goes beyond RBMs for several reasons:

• It is very likely that the same kind of transitions are also present in the early stages of deep Boltzmann machine learning, since the same perturbative analysis of refs. [8,9] performed for RBMs was recently reproduced for deep Boltzmann machines in Ref. [Ichikawa, Y., & Hukushima, K. (2022). Statistical mechanical study of deep boltzmann machine given weight parameters after training by singular value decomposition. Journal of the Physical Society of Japan, 91(11), 114001.].
• Cascades of phase transitions have also been observed in other generative models, e.g. in the training of the Gaussian mixture model [18-22] or very recently in the sample generation of diffusion models [23-25]. For the latter, it is conceivable that during the training dynamics the model of the score first learns the most prominent features, which needs less data (and hence less training time, at least in an online learning schedule) and then later the finer ones.

If we have the opportunity, we will discuss this point in a new "Limitations" section and try to make the paper more accessible by taking advantage of the additional page.

We thank the reviewer for their comments and for their careful reading of our paper. First of all, we will of course revise the notation to keep it consistent between the two parts and also make the transition and connections between the two parts clearer. We thank the reviewer for their advice.

Weaknesses:

• we will ensure that the notations are homogeneous in the whole paper for the final version
• Yes, the projections and projected weights in the theoretical section correspond to the discussed singular value decomposition described in the empirical evaluation. We agree that the paper would clearly improve by discussing and making explicit this connection. We will do it in the revised version. The susceptibility is analyzed in the theoretical section, see line 127. We show how the susceptibility diverges when the quantity $N^{-1} \sum_i w_i^2 \sim 1$ showing explicitly that the transition is of second order. We will add in the appendix a comment on the critical exponent.
• We will add structure to the numerical part in the final version, where an additional page can be added. It is true that the negative term is where the difficulty relies on when training the RBM. In practice, at the beginning of the training the sampling is very easy (the chains mix very quickly) and very crude algorithms such as one-step CD can be used to compute the gradient. When the first phase transition occurs the mixing time suddenly jumps and at that point it is important to correctly compute the negative term. In this case, the sampling can introduce a strong bias, the theoretical analysis of which is not easy since the weights are not small and, for instance, the 1-step CD would introduce non-linearity in the terms. Yet, we can check in practice, that the picture discussed in the paper does not change much when training with CD-1, as we show in the attached PDF for MNIST. Eigenvalues still turn on one by one, the first eigenvectors are very similar to those of the PCA and the growth of the eigenvalues is associated with a growth of the susceptibility along these directions. We propose to add a small subsection or appendix in the new version of the paper discussing out-of-equilibrium effects in more detail.

Questions:

• The directions onto which we project the gradient in the theoretical section indeed correspond to the principal directions of the dataset. In the case with one mode it is trivial (since the weight matrix is only a vector). In the multiple features case, it is easy to see by looking at the equation right after line 491, decomposing the correlation matrix of the dataset onto the vectors $\eta$ and $\bar{\eta}$. We will make it clearer in the final version of the work. Furthemore, the divergence of the sucesptibility in the theoretical section is already done, see line 126-127 and the equation in-between. We will be more clear about the exponent with which it diverges and how it relates directly with a second-order phase transition.
• The first phase transitions can be easily detected by monitoring the eigenvalues of the weight matrix (the singular values). In the model where the variables are $\\{0,1\\}$, the phase transitions occur when the singular values become higher than 4, while when the variables are $\\{\pm 1\\}$, the transition occurs when the singular values become higher than 1. In that sense, it is easy to detect the first phase transition. In practice, it is then possible to adapt the number of Monte Carlo steps nearby the transitions. The next phase transitions can also be detected by monitoring the second,third, … eigenvalues of the weight matrix, although it is not clear how many of them their will be. *We do not use a Langevin sampler. Since we work with discrete variables, we rather rely on Markov Chain Monte Carlo and Gibbs parallel sampling. To answer the reviewer's question: We do not observe spikes when updating the parameters unless we set a high learning rate. In this case, we do see spikes in the eigenvalues around the transition, but not in the experiences discussed in the current version of the manuscript.

Limitations:

We will add a new section Limitations, in which we will discuss the limitations of the current analysis: (i) the extent to which these results are applicable to more complicated models, (ii) the limitations of the analytical setup, and (iii) the limitations of the numerical setup (problems with insufficient MCMC sampling).