Fast, accurate training and sampling of Restricted Boltzmann Machines

We thank the referee for his comments and questions.

Q1: The main reasons for using RBMs are twofold: they are easy to interpret and are particularly well suited for tabular datasets with discrete variables such as DNA or protein sequences, both in terms of generation power and sample efficiency, which makes them particularly suitable for inference applications in physics, bioinformatics or neuroscience. Transformer-based generative methods require a large amount of training data. Most applications in science cannot afford this luxury, as experimental annotations are expensive or private. Moreover, once these models are trained, apart from their impressive generating power, it is almost impossible to extract understandable information from their parameters, which is the typical goal of data-driven scientific studies.

We copy below a list of recent references published in high-impact journals using the RBMs in biology:

• [Bravi, B., et al. Probing T-cell response by sequence-based probabilistic modeling. PLOS Computational Biology, (2021)]
• [Di Gioacchino et al. Generative and interpretable machine learning for aptamer design and analysis of in vitro sequence selection. PLoS computational biology, (2022)]
• [Bravi, B. et al. RBM-MHC: a semi-supervised machine-learning method for sample-specific prediction of antigen presentation by HLA-I alleles. Cell systems, (2021)]
• [S. Quiroz Monnens et al, The Recurrent Temporal Restricted Boltzmann Machine Captures Neural Assembly Dynamics in Whole-brain Activitye Life (2024)]
• [Bravi, B. A transfer-learning approach to predict antigen immunogenicity and T-cell receptor specificity. ELife, (2023)]
• [Tubiana, J.et al. Funneling modulatory peptide design with generative models: Discovery and characterization of disruptors of calcineurin protein-protein interactions. PLoS computational biology (2023)]
• [Bravi, B. Development and use of machine learning algorithms in vaccine target selection. npj Vaccines, (2024)]

We say that RBMs are interpretable precisely because their architecture is simple enough to allow a high level of analytical treatment. For example, their energy function can be rewritten in terms of an explicit many-body spin interaction system. This analogy makes it a very powerful inference engine that overcomes the standard Ising inverse approaches restricted to pairwise interactions. The structure of the probability distribution can also be explored with perturbative tools that allow, for example, to hierarchically cluster data for phylogenetic reconstruction (see [8]) or to extract possible mutational paths between data points [Mauri, E., Cocco, S., & Monasson, R. (2023). Mutational paths with sequence-based models of proteins: from sampling to mean-field characterization. Physical Review Letters, 130(15), 158402.].

Q2: the question of the novelty is answered in the Author rebuttal section

Weaknesses:

• memory costly PTT: It is true that the PTT sampling requires more memory than a standard sampling procedure, but it allows the generation of samples in models where it is simply impossible to generate samples in a reasonable amount of time otherwise. Memory is much easier to handle than time, especially in our case the models are moved to GPU memory one at a time, and only when we need to generate samples of them. Otherwise, they are kept in RAM. Even though the method requires more memory than standard Gibbs sampling, it is a very limited amount when compared to sample large models. Moreover, pre-training allows keeping this memory cost very low (in our case, only 7 times more memory was needed than standard sampling). We will mention this in the limitations section.

Minor comments:

• We will correct the problem with the auto-reference.
• We apologize for the non-anonymized GitHub link, we were not very sure how to do it otherwise.

On the PT In the attached PDF, we include a sketch explaining the difference between first and second order first transitions in the attached Fig. 1. In Fig. 2 D-E, we illustrate the effect of the disappearance of a peak upon lowering the temperature in the HDG dataset and the failure of the PT strategy of this dataset (Fig.2C). We also compare these results with the data sampled using 10^7 ordinary Gibbs MCMC steps (that should be long enough to guarantee equilibrium given the measures of Fig.4D).

The reviewer suggests using an improved version of PT using a non-uniform reference configuration. We were not aware of this modified algorithm, and we will discuss it in the new version of the manuscript and this reference, and compare it with the PTT. We agree with the referee that this reference probability will add a bump at the center of the probability distribution, which should to mitigate (not making them disappear) the effect of the first order transitions. Yet, our guess is that the efficiency of this algorithm must depend a lot on the properties of the dataset and the spatial location of the barriers, and it should struggle also with datasets with multiple clusters. We could not do a comparison during the rebuttal week due to the lack of available resources, but we will do it for the final version of the paper.

This trapping phenomena associated to the existence of first order transitions in temperature is reproduced when trying to compute the partition function using a temperature annealing, and it is beyond the failures observed in the computation of AIS. Instead, the trajectory tempering only crosses second order transitions, which means that a high temperature mode smoothly splits on two when lowering the temperature, which avoids the disappearance of modes.

We thank the reviewer for its comments and suggestions. 1. PT: The reviewer disagrees with us on the overall performance of standard PT in RBMs. We have several comments on this. First, the reviewer disagrees with the statement that PT is costly and assesses that the reliability of PT in RBMs has been well established in previous references. It is more or less textbook knowledge in statistical physics that the PT algorithm does not work in the presence of a 1st order phase transition (Ph-Tr). The reason is that if new distant modes simply emerge when the temperature decreases, the chains coming from high temperatures will never nucleate them if the number of MCMC steps performed between temperature-swaps is not long enough because the barriers between modes are large. As discussed theoretically in [41], this is the case when dealing with datasets clustered in a low-dimensional space (which is the case for genetic or protein data). Schematically, for these datasets, one tries to encode several distant modes at $\beta=1$, which are usually not symmetric around the zero magnetization. Since the temperature is only coupled to the energy and not to the entropy, a change in temperature leads to “1st order” Ph-Tr (because the modes are more sensitive to temperature changes the farther they are from the origin, see [41]). As a result, there are modes that simply disappear from the equilibrium measure when the temperature is increased. On the other hand, the chains located at these vanishing clusters at low temperatures melt into other clusters when they visit the high temperatures. When iterating the standard PT process, the chains occupying these clusters will gradually disappear, giving the false impression that the configurations belonging to these clusters are not typical at $\beta=1$, when they are. In contrast, the “2nd order” transitions are characterized by a merging of the modes. If the “2nd order" Ph-Tr occurs at temperature, this would mean that one hump gradually splits into two (or more) as the temperature decreases. The transitions that occur along the training trajectory are 2nd order transitions, and the transitions encountered when the temperature is increased (in clustered data sets) are 1st order Ph-Tr. Both statements are analytically proven in RBMs (see [26,27,28] for the trajectory, and [41] for the temperature).

We further elaborate on this point in the official comments.

2. LL: The LL was shown in Fig. 3 for MNIST and in Fig.5 for the HGD. We will include the LL of the remaining dataset in the Supplemental. We will add error bars to the curves by repeating several times the estimation. Following the reviewer’s suggestion, we will compare the different strategies to compute the LL in a controlled set-up where the partition function can be exactly computed. In particular, we will consider RBMs trained with the datasets discussed in the paper but with only a few hidden nodes (we can compute the exact LL as long as $N_h\le20$). We will also illustrate these cases to discuss the differences between PT and PTT.

3. It is not clear to us which artifacts the reviewer refers to. The samples in Fig. 3B are real samples from the learned RBM and not the probabilities for each pixel. The white dots that sometimes appear in the black background are rare events due to fluctuations at T=1. Further away, one starts to observe overfitting effects. If necessary, we can use the mean values of the pixels instead of samples, which smooths the images into gray levels.

4. The learning rate we used seemed to work well in our case and was maintained for all experiments and methods we used, for the sake of comparison. While we agree that keeping the learning rate constant is generally not the best strategy, we also believe that it is an easier setting to compare different methods and avoid adjusting more and more hyperparameters that would introduce hard-to-control effects. Moreover, adjusting the learning rate would similarly affect both strategies using PCD (i.e., with or without pre-training). We propose to show experiments with smaller learning rates in the Appendix, showing that the overall behavior between training schemes remains unchanged. We will also include error bars in the LL figures to illustrate the typical fluctuations of this measure.

5. We will include in the appendix the comparison between the samples obtained with PTT and those obtained after sampling $10^7$ “classical” MCMC steps (several hours long) the HGD model. We show this data in the attached PDF. They both perfectly match. We also show that the PT sampling fail to reproduce the upper cluster. We could not compare with the modified PT during this rebuttal week, but we will do it in the final version.

Languages:

1. We thank the referee: we should have defined better these quantities, and we will do it in final version.

2. RBMs are much easier to interpret than other models based on DNN. The main reasons for this are, first, that the weights of the latent variables interact directly with the visible nodes and therefore offer many possibilities to analyze them: SVD, looking at their value as a function of the visible nodes, etc. For example, Ref [15] shows that the protein functionality can be inferred from the RBM's latent space. Second, Ref. [7] showed that the RBM can be mapped to a system of interacting visible nodes, from which not only pairwise couplings between nodes can be inferred, but also higher-order couplings. Third, ref. [8] showed that the minima of the approximated free energy can be used to cluster data hierarchically. One can accurately approximate the free energy of the RBM because it has a very simple structure. There are also many recent papers that use the RBM for interpretive applications in biology. A list can be found in the rebuttal of reviewer SD84.

3. We will add the citation in the final version.

4. We will add more precise citation to the figures for the final version.

About the learning rate: As discussed with other referees, we believe that it would be better to reduce the learning rate for practical applications. To allow comparison between methods, we apply the same scheme to all training methods. A high learning rate exacerbates the problems associated with the lack of thermalization during PCD training, and not the other way around. Furthermore, the learning rate is not an absolute value. The optimal value for each machine strongly depends on the number of visible and hidden nodes used. Also on the number of MCMC steps. The machines discussed in Tielemann’s paper do not have the same number of nodes as ours, so the direct comparison of lr does not make sense.

Our samples show no artifacts. Tielemann's work shows the probabilities, not the variables. This is obvious because the numerical images are displayed in grayscale and not in black and white pixels. The display of probabilities is common in the literature for binary variables and image data. When we display the probabilities, we get a perfect black in the background, as in Tielemann, and our digits are also better shaped and more diverse when we work with the full MNIST and this lr. Finally, the (random) presence of white pixels could be more related to the value of the bias, which is set at the beginning of the training. The bias is initialized to reproduce the frequencies without weights before it starts learning. However, for pixels that do not fluctuate, a threhold is used to avoid infinite bias. The flipping of some spins is probably related to a “weak” threshold.

We thank the reviewer for its reply. Our paper and also ref. [41] focus on binary-binary RBMs in describing highly clustered datasets, where the different modes are separated by large barriers that are already visible at the PCA level, where typical training approaches drastically fail.

On the comparison of PT over PCD. The cases discussed in the rebuttal are good examples of situations where a PCD training can work better than a PT training (the version of PT where a temperature swap is followed by only one or a few Gibbs steps). In these data sets, the equilibrium measure can be sampled by performing many accumulated Gibbs steps, but not by performing many PT steps. This is because visiting the higher temperatures makes the sampling less ergodic than at T=1 because each time the chain visits high temperatures, a cluster is suppressed and standard sampling does not have enough time to re-nucleate it before visiting high temperatures.

You can solve this problem by significantly increasing the number of MCMC steps $k$ between temperature exchanges and/or the number of temperatures in the ladder, but in this case the overall acceleration of PT becomes less clear because you are forced to sample $N_\beta$ systems in parallel that are not particularly useful for anything else. The total cost for $k$ sampling steps at $\beta=1$ is then $N_\beta\times k$. For this reason, we have said that PT is computationally very demanding. Ref [40], for example, shows an order of magnitude improvement over AGS in the jumps' timescale in the case of MNIST 0/1, but at the cost of simulating 10 systems in parallel which is much slower as reported in Table S4. They also report no further acceleration when $N_\beta$ is increased from 10 to 1000. For proteins and the 2D Ising model, the situation is even worse.

Finally, we would like to mention that the ref https://www.sciencedirect.com/science/article/pii/S0004370219301948 cited by the reviewer analyzes machines trained with CD-25 and CD-1. The CD training produces very pathological models with a high LL but with a very strange dynamic behavior and, more importantly, models that are completely unable to generate proper samples (see ref. [21] or ref. [r2] suggested by the reviewer). We are not sure to what extent the results of this work are affected by this.

About the toy dataset: The dataset used in reference r1 provided by the reviewer is not a good example. It is simply too small (N=16). The divergence of times associated with bimodality increases with exp(N). For N=16, the thermalization problems associated with first-order phase transitions are much smaller than the ones we address in our work.

On first-order transitions in RBMs: First-order transitions in temperature are common when working with “clustered” data sets. This is not necessarily the case for most datasets where training algorithms are tested in RBMs.

We disagree with the reviewer's statement that first-order transitions do not occur in other versions of RBMs. Could the reviewer give us a reference stating that all transitions in the Gaussian binary RBM are second order? We would be very interested in such results for our research.

Both first and second order transitions have been reported in Gaussian Mixture Clustering, see for example [Lesieur, T., De Bacco, C., Banks, J., Krzakala, F., Moore, C., & Zdeborová, L. (2016, September). Phase transitions and optimal algorithms in high-dimensional Gaussian mixture clustering. In 2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton) (pp. 601-608). IEEE] or [J. Schneider (Phys Rev E 1998), First-order phase transitions in clustering], the latter showing precisely that first-order transitions in temperature are frequent in this problem. Moreover, using the method proposed in [25], it can be measured that an RBM trained with a simple data set consisting of a Bernouilli mixture of two clusters undergoes a first-order transition when $\beta$ is interpolated from 1 to 0.

On the other hand, the same arguments about RBMs learning PCA through second-order phase transitions in the early stages of training also apply to binary-binary, binary-Gaussian, and Gaussian-binary RBMs, as in Ref [Decelle, A., & Furtlehner, C. (2021). Restricted Boltzmann machine: Recent advances and mean-field theory. Chinese Physics B, 30(4), 040202].

On the PCA: The reviewer doubts the reliability of the PCA to best describe the multimodality of the data. We agree with the reviewer that the slowest mode could be related to a non-linear decomposition of the dataset and that, if this were the case, our method would not be particularly useful to avoid this. Our model is designed for data that is highly clustered in a low-dimensional space, which is typically already captured at the PCA level. This is the case of the datasets analyzed in this work (and the typical situation encountered in dna/protein datasets). We will try to make this distinction clearer and discuss this limitation in a new section: "Limitations".

The RBM learning is triggered by the encoding of the covariance matrix (and thus the PCA), it is only at a second stage when it starts to encode the higher-order statistics. The problem is that if the dataset is very clustered (as in the cases discussed in the work), the permanent chain used to estimate the negative part of the gradient remains very far from equilibrium from that point and on (chains typically get trapped in one of several isolated modes no matter their statistical weight), leading to very biased trainings. In these situations, our method can really make a difference. At variance, if the encoding of the first PCA components does not imply dividing the phase space in separate clusters (as it is the case for instance when dealing with the entire MNIST dataset), the encoding of these first directions do not imply insurmountable mixing times (in Fig. 3C of ref. [21] in the manuscript estimated that the first transition implied a growth of thermalisation times of order only 100 MCMC steps, which can be easily tackled by PCD), which means that a pre-training makes no particular speed-up in that case. The situation could be different for larger dimensions due to the critical-slowing down (ref.[28] explicitly computes the mixing times for training using resized versions of MNIST). For this reason, training the entire MNIST datasets is much easier than training the 0-1 MNIST one.

In our experience, however, for highly clustered datasets, the strongest time divergences appear at the beginning of training (when learning the PCA) because it is the moment where different modes are separated by very large barriers, later stages are typically much easier to sample because barriers become finite. This statement was quantified for some datasets used in this work in ref. [25] (see for instance Fig. 8 in [25] for the 0-1 MNIST). This is not the case for the entire MNIST, where mixing times softly grow along the training (see Fig. 3 [21]).

To make these statements clearer, we propose to show the explicit evolution of the exponential correlation times of the different models along the training time, which show very sharp jumps at an initial stage when the first eigenvalues of W grow, further decrease to much lower values that only very softly increase in time.

We thank the reviewer for its careful reading of the paper and for its constructive comments. We attempt to respond to the comments individually below.

The answers to the concerns about novelty of this paper have been answered in the section "Author Rebuttal", as they were shared by several reviewers.

Concerning the rest of comments:

On the PCA The reviewer doubts the reliability of the PCA to best describe the multimodality of the data. We agree with the reviewer that the slowest mode could be related to a non-linear decomposition of the dataset and that, if this were the case, our method would not be particularly useful to avoid this. Our model is designed for data that is highly clustered in a low-dimensional space, which is typically already captured at the PCA level. This is the case of the datasets analyzed in this work (and the typical situation encountered in dna/protein datasets). We will try to make this distinction clearer and discuss this limitation in a new section: "Limitations". We further elaborate our answer in the official comments.

MNIST: We did not discuss the entire MNIST dataset because it is much easier to train than the 0-1 version because barriers are much less pronounced in that case. MNIST dataset can be well-trained with a PCD-10 scheme, and in fact we do not observe any particular improvement of using the pretraining strategy (nor do we observe that it is getting worse because the quality is limited by the number of steps k used at the second stage of the training PCD-k). Being larger CelebA is a much better test for our algorithm because PCD-10 struggles to get decent images. We will include some analysis in the new version of the paper (and discuss the performance on the complete MNIST in the Supplemental).

Baseline methods: Referee wonders if we did any effort to optimize the baseline methods for a meaningful comparison. Our approach proposes only to include a pre-training before initiating the standard training methods. In this sense, we tried to optimize the baseline methods first and use the same hyperparameters later. Nevertheless, we agree with the referee that no discussion about the hyperparameters was included in the manuscript, and it is an important point. We will include some figures in the supplemental material showing the performance at different learning rates and number of added hidden nodes.

Comments:

• Interpretability Simple EBMs can be directly mapped to spin ferromagetic models in physics, for instance a binary Boltzmann machine without hidden units is formally identical to a (pairwise) disordered Ising model, and thus the trained model parameters can be interpreted as physical pairwise interactions between both variables. A strategy (typically referred as inverse Ising approach) widely used in data-driven interpretability approaches in biology, such as the direct coupling analysis (DCA) of protein families or the inference of neuron connections from spike data in neuroscience. Similarly, RBMs can be mapped to a slightly more complex spin model where interactions between spins are not only pairwise but multibody. The interactions between variables in a deep convolutional generative neural network cannot be easily written in that way, and the most prominent interactions between variables cannot be easily inferred from the analysis of the network. Furthermore, the analogy between neural networks and spin systems allows us to apply all the machinery developed in physics during the last decades to explore the equilibrium behavior of complex systems to analyze the learning or probe the probability distribution function. We will extend the paragraph to make it clear and add a section in the SM illustrating the mapping to a physical system.

• Following the same physical analogy, one can describe the encoding of features during the training as a thermal annealing process, along which each new mode in the model is introduced through a phase transition which is analogous to a paramagnetic to ferromagnetic transition in temperature. These kinds of transitions are denominated in physics as second order because they are characterized by a divergence on the second order statistics of the system's order parameter (in this case, the fluctuations of the magnetization - the variables projected on the each of the PCA components-). Second order phase transitions, also known as critical transitions, have many interesting properties such as universality (identical exponents rule the divergences in very different physical systems satisfying same basic symmetries and dimensions), the emergence of system size correlations between variables, and the apparition of also universal dynamical features, such as a sudden arrest of the dynamics, an effect known in physics as critical slowing down, which in simulations is reflected by a sharp growth of the mixing times with a power of the number of variables. These dynamical effects have recently been studied in the training of RBMs in [[D. Bachtis et al, Cascade of phase transitions in the training of Energy-based models, arXiv:2405.14689 2024]].

We acknowledge that phase transitions are not common knowledge in computer science and statistics (another referee complained about the same point), so we propose to include a whole pedagogical discussion in the Supplemental Material explaining the mapping between the RBM and a physical system and the dynamical characteristics of first and second order phase transitions as they are both important to understand the properties of the training process and the failures of the sampling one.

• We thank the referee for pointing out the interesting reference by Fischer&Igel that was unknown to us, and will be discussed in a further revision of this work.

• We also thank the good suggestion about comparing our estimation of the log-likelihood against the BAR method, which will be taken into account in the final version of the paper.