Fast training and sampling of Restricted Boltzmann Machines

Weaknesses

1. We appreciate the reviewer’s feedback on the accessibility of our paper. We know how important it is to make our work as understandable as possible for a wide audience. Therefore, we have included an appendix to explain the definitions and concepts from phase transition theory. We are happy to provide further explanation if there are particular terms or sections that could benefit from additional context. Although our work contains elements from the field of statistical physics, we believe that it also makes an important contribution to the fields of computer science and machine learning by providing a practical algorithm for training, evaluating, and sampling Restricted Boltzmann Machines that has numerous applications.

2. We appreciate the reviewer’s comments on the term "highly structured dataset". While there is no precise mathematical definition, we would like to clarify what we mean by this. By "highly structured" we refer to datasets that exhibit certain notable characteristics: (i) the presence of visible and well-separated clusters in PCA projections and (ii) the difficulty Monte Carlo methods have in jumping between isolated clusters, often due to excessively long mixing times. As for "unstructured" datasets, we show that our method is effective on these as well, although it does not offer much advantage for long training times, since PCD has difficulty thermalizing both with and without pre-training, rendering pre-training mostly useless for long training times. Fig. 15 shows the training of the full MNIST dataset, where the pre-training becomes useless (in terms of matching log-likelihoods) very early in the training, from $10^3$ updates. Fig. 15 illustrates the training of RBMs on the full MNIST dataset, where the advantage of pre-training diminishes as early as $10^3$ updates, making it ineffective in terms of matching log-likelihoods beyond this point. In contrast, for the MNIST 0-1 subset shown in Figure 2 (top), the pretrained model consistently outperforms the standard PCD model in terms of log-likelihood and in a proper balance of the different peaks of the histogram of the projected generated data. For the CelebA dataset (Figure 2), the pre-trained RBM also reaches a higher log-likelihood than the standard model, but both approaches seem to converge toward similar values over time.

3. Our paper focuses on the sampling problems that highly structured datasets pose when training and evaluating RBMs, and suggests strategies to overcome these problems. The long introduction focuses on both reviewing previous work and explaining the physical reason why many of these sampling strategies may fail at the different training stages. Based on these conclusions, we propose a more appropriate training and sampling strategy as well as a new method to compute log-likelihood during training. We want to stress that while Monte Carlo is an easy method to implement, in many cases, like in RBMs, it is extremely difficult to implement and control it correctly. While we do not believe that pre-training is the main contribution of our work (we rather think that it is the trajectory AIS measure of log-likelihood or the PTT), we have made innovations that are necessary to make the mapping of D&F 2021 usable in real datasets, which we present on page 5. We can try to make them clearer.

4. The work D&F 2021 propose a rather theoretical setting to learn a low-rank RBM and it is tested only on very simple, low-dimensional synthetic data sets with specific modes and regular features to cover the low-dimensional space, and mainly focuses on theoretical aspects. The first added value of the present work is to move from theory to practice, where many details need to be specified to make the technique work for arbitrary data. In this previous work, the algorithm was tested up to only 2 constrained directions. In our construction, we reach 4 constrained dimension thanks to including the possibility of adding a trainable bias, which also turns out to be crucial to train decent low-rank RBMs with image data. Moreover, D&F's work never controlled the quality of the low-rank RBM generated samples but only that of the low-rank Coulomb machines. It turned out that when dealing with real data, one needs to carefully correct the entropy term to ensure that true RBM equilibrium configurations can be obtained by the static Monte Carlo procedure, which is crucial to sample fast the trained machines, but also to properly train the low-rank RBMs. We will present these improvements in more detail in the final version.

We thank the referee for all the interesting question posed. We answer all them one by one below.

The distinction between theoretical claims and empirical findings is not clear. It would be beneficial for the authors to clarify which parts of the study are based on theoretical analysis and which are supported by numerical experiments, particularly in the context of related work. For instance, the first- and second-order phase transition claims pertain to equilibrium properties. However, it is unclear how these phase transitions are justified when updating parameters with limited samples.

Concerning the 1st and 2nd order phase transitions, we agree that the existence of such transitions is a purely equilibrium phenomenon. However, our study focuses on the associated dynamical effects in the sampling process, which should still manifest even if the transition is a crossover due to limited sample sizes or other out-of-equilibrium phenomena. At the dynamical level (where the temperature is changed, or the parameters adjusted), the effect of the transition takes place not only at the critical point by also in its neighborhood (for instance in 2nd order phase transition the increase in the correlation length manifests itself as we approach the transition). It is also important to note that the existence of at least the first learning transitions in finite minibatch training has already been confirmed in (Bachtis et al., 2024) using finite-size scaling techniques.

In Section 4, the paper introduces pre-training for low-rank RBMs with singular value decomposition (SVD)--based weights, aiming to avoid continuous phase transitions (second-order transitions) as structural patterns gradually emerge. It is further claimed that training can proceed quickly using the PCD method after post-pre-training. Could the authors provide a more detailed explanation for this intuition? Even if second-order transitions are avoided, if there are multiple stable clustered states, capturing multiple modes with the PCD method may be challenging and could introduce bias in the estimation. However, the paper claims, "Once the main directions are incorporated, training can efficiently continue with standard algorithms like PCD, as the mixing times of pre-trained machines tend to be much shorter than at the transitions." I believe that simulating clustered models with simple PCD often results in impractically long mixing times. Indeed, in Section 5.2, it is argued that mixing is very slow for AGS in clustered data.

The reviewer's intuition is correct. The effectiveness of pretraining strongly depends on the properties of the dataset. In many cases, even after learning the first modes, the mixing time following the transitions remains too long for PCD to function properly. This is why the benefits of pretraining cannot be extended throughout the entire training process. However, we typically observe that after crossing the first transitions, the mixing time decreases by several orders of magnitude compared to the values observed during the transitions, which, for some datasets, makes it feasible to safely continue using PCD for a portion of the training. This was for instance measured in (Béreux 2023). That said, the mixing time generally increases as training progresses, eventually driving the system out of equilibrium, requiring alternative algorithms. We have now tried to clarify this point in the revised version of the manuscript.

The statement "It’s also often ineffective with highly clustered data due to first-order phase transitions in EBMs, where modes disappear abruptly at certain temperatures, as discussed by Decelle & Furtlehner (2021a)" suggests that using PT becomes challenging because the learned RBM exhibits a first-order transition at specific temperatures. However, does the existence of a first-order transition in the learned RBM typically occur regardless of the statistical model being learned? For example, if learning a model without a first-order transition, such as the Ising model without a local field, does a first-order transition still arise in the learned RBM? This seems somewhat nontrivial.

Our claim is based on the work from D&F 2021 for which it is quite clear that a change of temperature (a global parameter $\beta$ in front of the energy) can unbalance the minima when we have a bias (see app. G line 1328 to 1348). From a mathematical point of view, it is possible to provide the further analysis that is now included in the new SI-H "FIRST ORDER TRANSITIONS IN RBMS" Section.

We consider a Curie-Weiss model over $s_i = \pm 1$ with $E= -\sum_{i<j} s_i s_j$ at a given temperature $\beta_{CW}$. We proceed with learning using Bernoulli $\sigma_i = \{0,1\}$ visible variables and one hidden gaussian node $\tau$ of variance $1/N$ (for simplicity) and with a field on the visible nodes. The Hamiltonian is given by $E = -\sum_{i} \sigma_i w_i \tau - \sum_i s_i \eta_i$. In such setting, the free energy is given by     $

        -f &= \frac{m_{\tau}^2}{2} - N^{-1} \sum_i\log\left[ 1 + \exp(w_i m_{\tau} + \eta_i)\right] \text{ with } m_{\tau} = \frac{1}{N}\sum_i w_i {\rm sigm}\left[1 + \exp(w_i m_{\tau} + \eta_i)\right]         %m_{\tau} &= \frac{1}{N}\sum_i w_i {\rm sigmoid}\left[1 + \exp(w_i m_{\tau} + \eta_i)\right]     $

    and we can identify the optimal learned parameters of the RBM     $

        w_i = 2 \sqrt{\beta_{CW}} \text{ and } \eta_i = -2 \beta_{CW}     $

    Starting from the optimal RBM, we can multiply the energy of the system by a factor $\beta$ and compute the free energy. We show on Fig. 9 of the SI H, we illustrate the result on ten different values of $\beta$ from $[0.8,1.05]$ and we clearly observe the presence of a first order transition: when $\beta$ is lowered the local minima with the large magnetization is gradually destroyed: first it is subdominant and then disappear. When $\beta$ is increased, we observe the exact contrary behavior, thus showing a clear first order transition at $\beta=1$.

We also demonstrate the same phenomenon using an RBM trained on a real dataset (HGD). Since this dataset is intrinsically low-dimensional, we efficiently apply the Tethered method from (Béreux 2023) to compute the potential and probability distribution as a function of $m$. In Fig. 11, we show that performing an annealing experiment by multiplying the energy by a factor $\beta$ reveals a first-order transition. For smaller values of $\beta$, the cluster (top-right) becomes sub-dominant and eventually disappears below $\beta=0.7$. For larger $\beta$ values, the central cluster vanishes, providing a clear indication of another first-order transition.

In the phase diagram of A. Decelle’s Thermodynamics of Restricted Boltzmann Machine and Related Learning Dynamics does not appear to be a first-order transition, and the AT line may suggests continuous phase transitions dominated by Full-step RSB. Thus, the claim regarding first-order transitions requires further elaboration. If a first-order transition is present, it would be essential to validate this by examining the free energy from the equilibrium state of the learned model, which could likely be accomplished by evaluating the partition function using the proposed method.

The phase diagram in this reference is presented within the proper context, where the absence of biases generally results in only second-order phase transition lines. For the occurrence of first-order phase transitions, we refer the reviewer to our response to their previous comment and the new SI-H.

Concerning the 1st and 2nd order phase transitions, we agree that the existence of such transitions is a purely equilibrium phenomenon.

In other words, is it accurate to interpret that, instead of focusing on the dynamic phase transitions arising in non-equilibrium processes, your algorithmic improvements are grounded in the properties of equilibrium states and their associated critical phenomena?

The mixing time decreases by several orders of magnitude compared to the values observed during the transitions, which, for some datasets, makes it feasible to safely continue using PCD for a portion of the training.

I am not entirely sure I understand. Once clusters form, the RBM's block Gibbs sampling becomes confined to those clusters. For example, with MNIST, starting with an image of the digit "1" and running the RBM rarely results in transitions to images of other digits. With respect to your comment that "the mixing time is reduced by several orders of magnitude compared to the values observed during the transition," could it be that pretraining biases the states of the learned clusters, leading to transitions that remain confined within those clusters? This wouldn’t typically be described as faster mixing of the Markov chain. In PCD and CD, Gibbs sampling is conducted in parallel using diverse initial values, with the randomness from these initializations enabling the practical generation of high-quality images.

We thank the reviewer for the careful reading and useful comments.

Weaknesses (we answer following the order of the bullet points)

I find it a bit challenging to identify the two main contributions in the paper as those are totally disentangled in their presentation between Sec. 4 and Sec. 5.2. I strongly recommend adding a list of bullet points at the end of section 1 to clearly list the contributions of work and crossref to the corresponding point in the paper. This would substantially help navigate the paper.

We thank the reviewer for this excellent suggestion. We have now added a list of bullet points at the end of Section 1 to clearly outline the contributions and included cross-references to the relevant sections for easier navigation. We have made a common answer about the novelties of this work in the official comments section.

I find that the structure of sections 5.2 and 5.2.1 can be improved. In particular, I find it confusing that Parallel Trajectory Tempering is introduced in section 5.2, and Parallel Tempering approaches are discussed in section 5.2.1. I find this logically inefficient as I believe that a more natural yet easier-to-follow flow would be to first introduce Parallel Tempering approaches and then explain what makes PTT different compared to existing approaches from the literature. As this is one fundamental contribution of this work I believe it is crucial to rework these sections such that the actual novelty emerges more clearly from the discussion.

We sincerely thank the reviewer for this valuable suggestion. In the revised manuscript, we have reorganized the discussion to make it more linear and better aligned with the proposed structure.

The discussion around eq. (4) is rather crucial for the paper as it represents one of the main contributions of this work. Currently, the novelty with respect to Decelle and Furtlehner (2021a) is not very clear to me, and I would appreciate it if the authors could elaborate more on this. Moreover, what's the intuition behind the "magnetizations" along each of the singular vectors? Is there any correspondence with the magnetization as a physical observable? As far as I understand, those should be the projections along the unitary vectors of the visible variable. Is that correct? If all my understanding is correct, then the new contribution of this work is to use a bias initialization along a direction , which augments the dimensionality of the system by one in the bias direction. If all above is still all correct, I wonder the following:

How beneficial is to have such an augmented direction for the bias compared to the naive approach proposed in Decelle and Furtlehner (2021a)?

Have the authors conducted any ablation studies to compare the differences in performances between Decelle and Furtlehner (2021a) and their new approach from an empirical standpoint?

The benefit of having an augmented direction is two-fold:

• The method requires evaluating a discretized multidimensional integral. Having the augmented direction allows one to learn an extra dimension without the computational overhead of discretizing it since this direction is independent of the rest of the model.
• This extra dimension is not learned using the hidden features of the RBM but only with the bias, resulting in RBMs with far fewer hidden nodes when trained with this additional direction.

The work by D&F never attempted to train the low-rank RBM on real datasets, focusing solely on simple artificial clustered data. Without incorporating the learning of biases and additional directions, their approach fails to produce low-rank RBMs capable of adequately generating images, even for a basic dataset like MNIST01. Moreover, the D&F method faces issues in generating samples from the low-rank RBM, as the generated samples do not align with the statistics of the model's equilibrium distribution. This discrepancy can be verified by running a standard MCMC simulation on the samples generated using the static Monte Carlo method proposed by D&F. To resolve this issue, it is necessary to correct the entropy computation to address the mismatch caused by the fact that the generated samples do not have the same exact magnetization as those sampled from $p(m)$.

Regarding the second question, we have not yet conducted a systematic study on the performance impact of adding more directions. Preliminary results suggest that the improvement strongly depends on the dataset. For MNIST01, having at least 4 directions was crucial to achieve better performance than PCD training, whereas for the HGD dataset, 2 or 3 directions appeared to yield similar results. We plan to include an ablation study in an updated version of the paper. Concerning the term magnetization, the referee is right it comes from the study of spins systems in physics, we can here see it simply as the projection of the spins on a particular direction as the referee's said.

Please state what the free parameters to learn are in equation 3. If u and $\bar{u}$ are the singular directions, then the free parameters would be $w_\alpha$? In general I found the description of the low-rank approach unclear and this important section needs work to make it simpler and more clear to the reader.

In Eq. (3), for the low-rank method, the matrix $\boldsymbol{u}$ is determined using the PCA while $\bar{\boldsymbol{u}}$ and $\eta_a$ are randomly initialized in the Coulomb Machine and not learned. So only $w_\alpha$ and $\theta_\alpha$ are indeed adjusted through gradient ascent. We have included some sentences explaining this in the new version. Nevertheless, we will try to make an effort to improve the explanation of the pre-training method.

For figure 14 it would be useful to show the distribution of the PCA projected data to see how well the RBM matches the projected data distribution.

Thank you for the comment. The new figure is now in the new version of the paper.

• I find it hard to follow why the authors are considering different sampling schemes and therefore what the aim of this section is. I presume this is considering alternative sampling approaches after the low-rank pre-training has been applied. However, I struggle to follow a clear recommendation or conclusion as to which method might be more suitable*

We believe the reviewer is referring to Section 5, where we compare various sampling methods—Alternating Gibbs Sampling (AGS), Parallel Tempering (PT), Stacked Tempering (ST)—with the algorithm proposed in this paper: Parallel Trajectory Tempering (PTT).

In this section, we assume a well-trained model (obtained using the pretraining + PCD procedure described earlier) and focus on sampling independent equilibrium configurations. For highly structured (i.e., clustered) datasets, this task is inherently challenging. In standard local Monte Carlo methods like AGS, the mixing time is ruled by the typical times for chains to jump between clusters, resulting in very long equilibration times. In these kinds of models, neither PT performs well because of the existence of first order transitions in temperature: the fact that some modes of the probability measure just disappear at neighboring temperatures. This last phenomenon is now explained in detail in a new SI H section about first order transitions.

To address this problem, we propose the PTT, an algorithm designed to drastically reduce equilibration time by leveraging the training trajectory of the model. To evaluate the efficiency of sampling algorithms, we track the average number of cluster jumps made by the Markov chains in average. In Figure 4, we compare PTT to previously developed methods, including AGS, ST, PT, and its variations. Our results demonstrate that PTT is significantly more efficient—often by several times or even orders of magnitude—compared to other methods.

We have revised and reorganized the entire section in the new version of the manuscript to clarify the contributions and enhance the discussion.

In the conclusion the authors claim to have introduced a method that enables "precise computation of log-likelihood". I cannot see anything in the main text that relates to this. There is no experiment I can see that measures the quality of the log-likelihood approximation. Please give some evidence to support this assertion.

We are quite puzzled by this comment. Section 5.1 in the previous version of the manuscript was devoted to using the training trajectory as an annealing process to estimate the log-likelihood (LL). We explained how this approach allows the LL to be computed online during training, as well as offline or through thermodynamic integration of configurations obtained with the PTT sampling.

Figs. 3 and 11–14 in the SI (Figs. 3 and 15–18 in the current version) systematically compare in all datasets the trajectory-based LL estimates with the exact values in RBMs with a small number of hidden nodes, where all states can be enumerated. These figures also contrast our trajectory-based method with the standard temperature AIS approach, which is not only less precise but also significantly more computationally expensive than the trajectory-based method.

We have included a set of bullet points in Section 1 of the paper to highlight that this is one of the important contributions of this work.

We again thank the reviewer for their careful reading and constructive comments. We will attempt to respond to all comments individually below.

The paper suffers from a lack of clarity of presentation and lack of clarity of novelty.

We apologize for any lack of clarity in the paper, which may stem from differences in language between communities. We have tried to reorganize and rephrase some parts of the manuscript to make it clearer overall. However, if the reviewer could point out certain unclear sections, we would be happy to revise them to make them more accessible.

Regarding the novelty, as discussed in our general response, we respectfully disagree with the reviewer’s assessment. While the reviewer’s focus appears to be on the pretraining proposal, we consider our primary contributions to be the new algorithms for estimating the log-likelihood and accelerating the sampling process. These innovations provide substantial improvements over the state-of-the-art and are broadly applicable to energy-based models beyond RBMs. To make these diverse contributions clear and accessible, we will highlight them in bullet points at the beginning of the paper.

The paper mentions that the idea of a low-rank approach has already been used by others and it's unclear to me what novelty there is in any of the sampling approaches used after the pre-training phase.

To our knowledge, the low-rank initialization approach has only been proposed theoretically in Decelle-Furtlehner (2021) and has not been applied to real data or used as a pre-training method. Additionally, leveraging the training trajectory to compute the log-likelihood and perform parallel tempering is a novel contribution of this work. This approach significantly enhances the performance of previously optimized methods across all datasets considered, marking an important new contribution of this work.

In terms of presentation, there are notational inconsistencies and a general lack of clarity in terms of the main ideas

We thank you again for your careful reading and your suggestions for improvement. We have revised the paper to address almost all of your comments. In doing so, we have ensured consistent spelling and improved the clarity of key ideas in the new, quick version. The remaining suggestions will be incorporated later this week or in the final version.

Fundamentally the approach of fitting a constrained model seems straightforward and indeed I believe there is a simple way to compute the projected distribution in the PCA space (using the Fourier integral representation of the Dirac delta function) which the authors do not discuss.

We agree with the reviewer that $p(m)$ could indeed be computed using a Fourier representation rather than large deviations. However, the Fourier transform would still require a saddle point approximation, effectively equivalent to our large $N$ approximation in the large deviation formalism. Importantly, the Coulomb machine mapping is not primarily aimed at estimating $p(m)$ for the RBM—this can be done through various methods. Instead, it serves to make the training process a convex optimization problem, as training the parameters of $p(m)$ directly in the RBM framework is often unstable and challenging. We will clarify this point in the paper.

QUESTIONS

Figure 1 isn't very easy to parse. For example, the panel on race is placed more in the Mickey column than the human genome column.

Thank you for pointing this out. We have adjusted the layout to improve clarity and ensure each panel aligns with the intended column.

Please clarify the difference between "model averages" and "observable averages" and the difference between using $N_s$ independent MCMC processes and R parallel chains. Please clarify for the reader the meaning of $<v_ih_a>_D$

Thank you for pointing out these areas needing clarification. We will ensure these terms and variables are properly defined and consistently used in the final version of the paper. We have included a sentence explaining the meaning of $<v_ih_a>_D$ below Eq.(2) in the new version.

Section 4: It is not correct that it is possible to train ``exactly" an RBM with a reduced number of modes. Approximations are required, as explained in the supplementary material.

We agree with the reviewer regarding the comment on the exactness of low-rank training. The training is performed precisely rather than exactly, with residual errors arising from the asymptotic approximation at finite $N$ and the finite resolution of the mesh. We will systematically revise the paper to remove all the references to the exact training.

We thank the reviewer for the suggestion, now we understand better. The reason why we want to encode the first PCA components in the $W$ matrix is to mimic what standard maximum likelihood training does in the first moments of training time, as shown in several previous theoretical and numerical works, which is encoding the first PCA components. We proposed to use this pre-training to simply skip this first part of the training and directly start training where one should have arrived if the gradient had been accurately estimated in a standard training. The reviewer suggests starting the training with another general low-rank model that is not necessarily associated with PCA. We do not know what would happen if a standard PCD training was restarted from there (whether or not the training dynamics would try to re-learn the initial PCA components and thus destroy the pre-trained model). The reliability of this idea should be tested in practice as we have no experience with it. Yet we expect this to work poorly as a projection of the data on small set of random directions usually don't show its structure which instead tend to appear on the first PC. In that case the pre-training become useless.

The low-rank construction method we used can be performed with any orthonormal basis. We use the PCA because we want to mimic the training.

Dear Reviewer,

We believe that we have addressed most of your questions and concerns about the weaknesses of our work. We would also like to emphasize that while your suggestions on low-rank construction approaches could be considered for alternative methods of constructing rank-low RBMs, the actual development of low rank RBMs was not a key contribution of our work.

To clarify the focus of the work, we have prepared a revised version of the manuscript. In this new version, we have placed more emphasis on sampling and log-likelihood computation to illustrate how these processes exploit the continuous nature of phase transitions during learning. We have also included a new figure (Fig. 4AB) and a new subsection discussing the limitations of the standard parallel tempering algorithm in the context of first-order transitions. In addition, we have added a new section in the SI (Section H) where we analyze in detail the landscape of probability density in a simplified environment where temperature changes.

We hope that these updates add more clarity and depth to our work, and ask you to reconsider our manuscript in light of these changes.