Fast training and sampling of Restricted Boltzmann Machines

## Abstract

Abstract Restricted Boltzmann Machines (RBMs) are effective tools for modeling complex systems and deriving insights from data. However, training these models with highly structured data presents significant challenges due to the slow mixing characteristics of Markov Chain Monte Carlo (MCMC) processes. In this study, we build upon recent theoretical advancements in RBM training, focusing on the gradual encoding of data patterns into singular vectors of the coupling matrix, to significantly reduce the computational cost of training (in very clustered datasets) and evaluating and sampling in RBMs in general. The learning process is analogous to thermodynamic continuous phase transitions observed in ferromagnetic models, where new modes in the probability measure emerge in a continuous manner. Such continuous transitions are associated with the critical slowdown effect, which adversely affects the accuracy of gradient estimates, particularly during the initial stages of training with clustered data. To mitigate this issue, we propose a pre-training phase that encodes the principal components into a low-rank RBM through a convex optimization process. This approach facilitates efficient static Monte Carlo sampling and accurate computation of the partition function. Furthermore, we exploit the continuous and smooth nature of the parameter annealing trajectory to achieve reliable and computationally efficient log-likelihood estimations, enabling online assessment during the training process, and proposing a novel sampling strategy termed parallel trajectory tempering that outperforms previously optimized MCMC methods.
Our results demonstrate that this innovative training strategy enables RBMs to effectively address highly structured datasets that conventional methods struggle with. Additionally, we provide evidence that our log-likelihood estimation is more accurate than traditional, more computationally intensive approaches in controlled scenarios. Moreover, the parallel trajectory tempering algorithm significantly accelerates MCMC processes compared to existing and conventional methods.

### 1 Introduction

Energy-based models (EBMs) are a long-established approach to generative modeling, with the Boltzmann Machine (BM)  (Hinton & Sejnowski, ) and Restricted Boltzmann Machine (RBM) (Smolensky, ; Ackley et al., ) being among the earliest examples. They offer a clear framework for capturing complex data interactions and, with simple energy functions, can reveal underlying patterns, making them especially valuable for scientific applications. Pairwise-interacting models (Nguyen et al., ) have been widely used for decades for inference applications in fields like neuroscience (Hertz et al., ) and computational biology (Cocco et al., ). More recent research has shown that RBMs can also be reinterpreted as many-body physical models (Decelle et al., ; Bulso & Roudi, ), analyzed to identify hierarchical relationships in the data (Decelle et al., ), or used to uncover biologically relevant constituent features (Tubiana et al., ; Fernandez-de Cossio-Diaz, ). In essence, RBMs present a great compromise between expressibility and interpretability: they are universal approximators capable of modeling complex datasets and simple enough to allow a direct extraction of insights from data. This interpretative task is particularly challenging when dealing with more powerful models such as generative convolutional networks EBMs (Xie et al., ; Song & Kingma, ) and other state-of-the-art models in general. At this point, it is also important to stress that RBM are not only interpretable, but they perform particularly well when dealing with tabular datasets with discrete variables, and where data is not abundant, which is often the case in genomics/proteomics or neural recording datasets. For instance, RBMs have shown to be among the most reliable approaches to generate synthetic human genomes (Yelmen et al., ; ) or to model inmunological response (Bravi et al., ; ; ) or neural activity (van der Plas et al., ; Quiroz Monnens et al., ).

Despite the appealing modeling properties of RBMs, they are notoriously difficult to train, a challenge common to EBMs in general. The main difficulty arises from the computation of the log-likelihood gradient, which requires an ergodic exploration of a dynamically evolving and potentially complex free energy landscape using Markov Chain Monte Carlo (MCMC) processes. Recent studies have shown that models trained with non-convergent MCMC processes suffer from out-of-equilibrium dynamic memory effects (Nijkamp et al., ; ; Decelle et al., ). This dynamical behavior can be explained analytically using moment-matching arguments  (Nijkamp et al., ; Agoritsas et al., ).
While exploiting these effects can yield fast and accurate generative models, even for highly structured data (Carbone et al., ) or high-quality images with RBMs (Liao et al., ), this approach results in a sharp separation between the model’s Gibbs-Boltzmann distribution and the dataset distribution, thereby undermining the interpretability of the model parameters as investigated in Agoritsas et al. (); Decelle et al. (). Thus, to extract meaningful information from datasets using RBMs, it is essential to ensure proper mixing of the chains during training, in short, one needs equilibrium models.

Both the ability to train an RBM in equilibrium and to generate convincing new samples from its equilibrium measure strongly depend on the dataset in question. For typical image datasets such as MNIST or CIFAR-10, good RBMs can be obtained by increasing the number of MCMC steps. However, this approach is no longer feasible for highly structured datasets (Béreux et al., ). Natural datasets from which one seeks scientific insights are often highly structured, such as genomics/proteomics data, neural-recordings or low-temperature configurations in physical systems. These datasets typically exhibit distinct clusters, identifiable via principal component analysis (PCA) which form distant groups of similar entries often associated to different labels. We show an example of the PCA of the 5 datasets we will be studying in this work in Fig. 1. The first 4 datasets are very clustered, while CelebA is not, as it presents a rather compact PCA. A very similar (non-clustered) situation is obtained if the entire MNIST dataset is considered.
Details about these datasets are given in the caption of figure 1. During training, the model must evolve from an initial normal distribution to an increasingly multimodal distribution. Sampling from multimodal distributions is particularly challenging because the mixing times are controlled by the transition times between modes. But this is not the only difficulty. In RBMs, these distant modes can be shown analytically(Decelle et al., ; ; Bachtis et al., ) to gradually appear in the model’s equilibrium measure during training by a process formally analogous to the ferromagnetic phase transition in simple spin models, the process by which the magnetic spins of a material spontaneously align with a fixed pattern when the temperature goes below a certain transition temperature. The RBM training is characterized as a sequence of such phase transitions in which the specialization of the model parameters plays the role of temperature reduction.
These pattern-encoding phase transitions are continuous phase transitions (also known as second-order transitions or critical transitions), and these types of transitions are associated with a sudden arrest in the characteristic times of the visible variables, a phenomenon known as critical slowdown, which causes the MCMC mixing time to diverge with a power of the system size when the transition point is approached. The critical slowdown effects that occur in the training trajectory were recently studied in Bachtis et al. (). A brief introduction to phase transitions and the associated dynamical effects can be found in the SI G.

This sampling challenge not only hinders the training process, but also limits the model’s ability to generate new samples. Obtaining new and independent configurations would require an impractically large number of sampling steps.

Figure: Figure 1: Datasets. Panels A-E display 5 distinct datasets projected onto their first two PCA components. In some instances, the dots are color-coded to indicate different labels. In A, the MNIST 01 dataset, featuring images of the digits 0 and 1 from the complete MNIST collection, along with a few sample images. In B, the “Mickey” dataset, an artificial dataset whose PCA representation forms the shape of Mickey Mouse’s face. In C, the Human Genome Dataset (HGD), which consists of binary vectors representing mutations or non-mutations for individuals compared to a reference genome across selected genes. In D, the Ising dataset, showcasing equilibrium configurations of the 2D ferromagnetic Ising model at an inverse temperature of $\beta=0.44$. In E, the CelebA dataset in black and white, resized to 32x32 pixels. For more details on these datasets, please refer to the SI.
Refer to caption: x1.png

### 2 Related work

Training EBMs by maximizing log-likelihood has long been a challenge in the community (LeCun et al., ; Song & Kingma, ). EBMs gained popularity with the introduction of the contrastive divergence algorithm (Hinton, ), in which a set of parallel chains is initialized on independent examples in the minibatch and the MCMC process iterates for a few steps. Despite its widespread use, this algorithm yields models with poor equilibrium properties that are ineffective as generative models (Salakhutdinov & Murray, ; Desjardins et al., ; Decelle et al., ). An improvement is the persistent contrastive divergence (PCD) algorithm (Tieleman, ), which maintains a permanent chain in which the last configurations used to estimate the previous gradient update are reused. PCD acts like a slow annealing process improving gradient estimation quality. However, it often fails on clustered data as the statistical properties of the permanent chain quickly move away from the equilibrium measure and degrade the model (Béreux et al., ). This problem, which is primarily related to phase coexistence, can be addressed with constrained MCMC methods if appropriate order parameters are identified. For RBMs, these order parameters are related to the singular value decomposition of the model coupling matrix, which enables efficient reconstruction of multimodal distributions (Béreux et al., ). Although this method is effective for evaluating model quality, it is too computationally intensive to be used in training, even if it leads to models with good equilibrium properties. Other optimized MCMC methods, such as the Parallel Tempering (PT) algorithm  (Hukushima & Nemoto, ) (in which the model is simulated in parallel at different temperatures which are exchanged once in a while using the Metropolis rule), have been proposed to improve the training process (Salakhutdinov, ; Desjardins et al., ) and show important improvements in certain datasets (Krause et al., ).

However, PT is expensive because it requires sampling at multiple temperatures while using samples from only one model. 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 (). This leads to sharp drops in acceptance rates for exchanges between temperatures near the transition. For this reason, first-order transitions typically hinder PT, as many intermediate temperatures are needed to allow sufficient acceptance rates in the exchanges. In contrast, continuous transitions, as seen during the training trajectory, cause modes to gradually separate as the model parameters are changed.

The population annealing algorithm, which adjusts the statistical weights of parallel chains during learning each time the parameters are updated, has been proposed as an alternative approach (Krause et al., ). Reweighting chains using non-equilibrium physics concepts like the Jarzynski equation has also been explored (Carbone et al., ). However, both methods struggle with highly structured data, as frequent reweighting (which involves duplicating only a few samples) biases the equilibrium measure due
to an insufficient number of effective independent chains. To prevent chains from becoming too correlated around transition regions during training, either more sampling steps or a lower learning rate is required, which leads to very long training times to ensure equilibrium. In our case, we tried to used both methods with our datasets but the training becomes too unstable or two slow if learning rate and number of Gibbs steps are auto adjusted. Another approach uses EBMs as corrections for flow-based models that are easier to sample (Nijkamp et al., ), simplifying learning but sacrificing the interpretability of the energy function, a key goal of this work. Alternatively, evolving flow models can serve as fast sampling proposers for the EBM (Grenioux et al., ), but this requires training two separate networks and the acceptance of the proposed moves is significantly reduced as the EBM becomes more specialized.

A recent method for RBMs, called “stacked tempering” (Roussel et al., ), significantly accelerates sampling by training smaller RBMs using latent variables from previous models. This allows for fast updates through a PT-like algorithm, and the authors demonstrated that it is much faster than standard PT. However, while this method is effective for sampling a trained model, it is not suitable for learning, as it requires sequentially training multiple models stacked on top of each other. Also for RBMs, it has recently been shown in a theoretical work that it should be possible to train a low-rank RBM that accurately reproduces the statistics of the data projected along the $d$ first data principal directions through a convex and very fast optimization process (see Decelle & Furtlehner () and the discussion below) using a mapping to another model, the so-called Restricted Coulomb Machine. This low-rank model can be seen as a good approximation to the correct RBM needed to describe the data, and has the nice property that it can be efficiently sampled via a static Monte Carlo process. While this work showed the efficiency of this training process to describe simple synthetic datasets embedded in only one or two dimensions, the use of this algorithm to obtain reliable low-rank RBMs with real data requires further developments that will be discussed in this paper, in particular: (i) the training of a bias, which is crucial to gain almost for free an additional dimension, and to deal with image datasets, and (ii) a tuning of the RBM fast sampling process, which was not properly investigated in the original paper.

In this paper, we show how starting the RBM training with a low-rank RBM in very clustered datasets can drastically improve the quality of the models by circumventing the initial dynamical arrest effects that occur during the training process, at a very similar computational cost. We also show that the training trajectory can be leveraged to estimate the log-likelihood at a very little computational effort, so that both the training and test set log-likelihood can be easily and reliably computed online during training. Furthermore, we propose an efficient sampling method, which we call Parallel Trajectory Tempering (PTT). This method outperforms optimized MCMC methods such as the traditional Parallel Tempering (PT) or the “Stacked Tempering” (Roussel et al., ) only requiring storing the parameters of the model at a few stages during training. By combining these strategies, we are able to train equilibrium models that accurately capture different modes in datasets prone to mode collapse using standard methods, easily evaluate RBMs during the training process, and easily sample from multimodal distribution.

### 3 The Restricted Boltzmann Machine

The RBM is composed by $N_{\mathrm{v}}$ visible nodes and $N_{\mathrm{h}}$ hidden nodes. In our study, we primarily use binary variables $\{0,1\}$ or $\pm 1$ for both layers. The two layers (visible and hidden) interact via a weight matrix $\bm{W}$, with no direct couplings within a given layer. Variables are also adjusted by visible and hidden local biases, $\bm{\theta}$ and $\bm{\eta}$, respectively. The Gibbs-Boltzmann distribution for this model is expressed as

$$ $p(\bm{v},\bm{h})=\frac{1}{Z}\exp\left[-\mathcal{H}(\bm{v},\bm{h})\right]\text{ where }\mathcal{H}(\bm{v},\bm{h})=-\sum_{ia}v_{i}w_{ia}h_{a}-\sum_{i}\theta_{ i}v_{i}-\sum_{a}\eta_{a}h_{a},$ (1) $$

where $Z$ is the partition function of the system. As with other models containing hidden variables, the training objective is to minimize the distance between the empirical distribution of the data, $p_{\mathcal{D}}(\bm{v})$, and the model’s marginal distribution over the visible variables, $p(\bm{v})=\sum_{\bm{h}}\exp\left[-\mathcal{H}(\bm{v},\bm{h})\right]/Z=\exp
\left[-H(\bm{v})\right]/Z$. Minimizing the Kullback-Leibler divergence is equivalent to maximizing the likelihood of observing the training dataset in the model. Thus, the log-likelihood $\mathcal{L}=\langle-H(\bm{v})\rangle_{\mathcal{D}}-\log Z$ can be maximized using the classical stochastic gradient ascent. For a training dataset $\mathcal{D}=\{{\bm{v}}^{(m)}\}_{m=1,\dots,M}$, the log-likelihood gradient is given by

$$ $\displaystyle\frac{\partial\mathcal{L}}{\partial w_{ia}}$ $\displaystyle=\langle v_{i}h_{a}\rangle_{\mathcal{D}}-\langle v_{i}h_{a} \rangle_{\rm RBM},\;\frac{\partial\mathcal{L}}{\partial\theta_{i}}=\langle v_{ i}\rangle_{\mathcal{D}}-\langle v_{i}\rangle_{\rm RBM},\;\frac{\partial \mathcal{L}}{\partial\eta_{a}}=\langle h_{a}\rangle_{\mathcal{D}}-\langle h_{a }\rangle_{\rm RBM},$ (2) $$

where $\left\langle{\cdot}\right\rangle_{\mathcal{D}}$ denotes the average with respect to the entries in the dataset, and $\left\langle{\cdot}\right\rangle_{\mathrm{RBM}}$ with respect to $p(\bm{v},\bm{h})$. Since $Z$ is intractable, the model averages in the gradient are typically estimated using $N_{\mathrm{s}}$ independent MCMC processes, and observable averages $\left\langle{o(\bm{v},\bm{h})}\right\rangle_{\mathrm{RBM}}$ are replaced by $\sum_{r=1}^{R}o(\bm{v}^{(r)},\bm{h}^{(r)})/R$, with $(\bm{v}^{(r)},\bm{h}^{(r)})$ being the last configurations reached with each of the $r=1,\ldots,R$ parallel chains. To obtain reliable estimates, it should be ensured that each of the Markov chains mix well before each parameter update. However, ensuring equilibrium at each update is impractical, slow and tedious. The common use of non-convergent MCMC processes is the cause of most difficulties and weird dynamical behaviors encountered in training RBMs as discussed in Decelle et al. (). In order to minimize out-of-equilibrium effects, it is often useful to keep $R$ permanent (or persistent) chains, which means that the last configurations reached with the MCMC process used to estimate the gradient at a given parameter update $t$, $\bm{P}_{t}\equiv\{(\bm{v}^{(r)}_{t},\bm{h}^{(r)}_{t})\}_{r=1}^{R}$, are used to initialize the chains of the subsequent update $t+1$. This algorithm is typically referred as PCD. In this scheme, the process of training can be mimicked to a slow cooling process, only that instead of varying a single parameter, e.g. the temperature, a whole set of parameters $\bm{\Theta}_{t}=(\bm{w}_{t},\bm{\theta}_{t},\bm{\eta}_{t})$ are updated at every step to $\bm{\Theta}_{t+1}=\bm{\Theta}_{t}+\gamma\bm{\nabla}\mathcal{L}_{t}$ with $\bm{\nabla}\mathcal{L}_{t}$ being the gradient in equation 2 estimated using the configurations in $\bm{P}_{t}$, and $\gamma$ being the learning rate.

Typical MCMC mixing times in RBMs are short at the beginning of training but increase as it progresses (Decelle et al., ). Recent analytical and numerical studies show that the mixing times exhibit sharp, localized peaks when the model encodes a new mode during training (Bachtis et al., ). This phenomenon, known in physics as critical slowing down, is related to the system going through critical phase transitions. For RBMs, the first phase transitions during maximum likelihood training are associated with the encoding of the principal components of the dataset into the singular vectors of the weight matrix $\bm{W}$ (Decelle et al., ; ; Decelle & Furtlehner, ; Bachtis et al., ). Since these directions can be determined before training by performing PCA on the data, one can pre-train an RBM to encode them accurately. This avoids crossing the initial phase transitions encountered during the standard training and the associated degradation of the gradient calculation.

### 4 The low-rank RBM pretrained

In Decelle & Furtlehner (), it was shown that it is possible to train exactly an RBM containing a reduced number of modes in the weight matrix $\bm{W}$ by exploiting a mapping between the RBM and a Restricted Coulomb Machine and solving a convex optimization problem, see the SI A.2. In other words, it is possible to train an RBM whose coupling matrix has this simplified form

$$ $\bm{W}=\sum_{\alpha=1}^{d}w_{\alpha}\bar{\bm{u}}_{\alpha}\bm{u}_{\alpha}^{\top },\qquad\text{with}\qquad(\bm{u}_{\alpha},\bar{\bm{u}}_{\alpha})\in{\mathbb{R} }^{N_{v}}\times{\mathbb{R}}^{N_{h}},$ (3) $$

where $\{w_{\alpha}\}$ are the singular values of the coupling matrix and the right singular vectors $\{\bm{u}_{\alpha}\}_{\alpha=1}^{d}$ correspond exactly to the first $d$ principal directions of the data set. Under this assumption, $p(\bm{v})$ a function of only $d$ order parameters, given by the
magnetizations along each of the $\boldsymbol{u}_{\alpha}$ components, i.e. $m_{\alpha}(\bm{v})=\bm{u}_{\alpha}\cdot\bm{v}/\sqrt{{N_{\mathrm{v}}}}$. In particular,

$$ $H(\bm{v})=-\sum_{\alpha=0}^{d}\theta_{\alpha}m_{\alpha}-\sum_{a}\log\cosh\left (\sqrt{N_{\mathrm{v}}}\bar{u}_{a}\sum_{\alpha=1}^{d}w_{\alpha}m_{\alpha}+\eta_ {a}\right)=\mathcal{H}(\bm{m}(\bm{v})),$ (4) $$

where $\bm{m}=\left(m_{1},\ldots,m_{\alpha}\right)$ and $\theta_{\alpha}$ is the projection of the visible bias on the extended $\bm{u}$ matrix. In order to initialize the visible bias as suggested in Montavon & Müller (), a direction $\bm{u}_{0}$ is used for the magnetization $m_{0}$ that is only present in the bias term. This is an improvement of this work with respect to Decelle & Furtlehner (), as it allows us to gain almost for free one additional direction.

The optimal parameters can be determined by solving a convex optimization problem, detailed in SI A.2. Once the model is trained, we obtain a probability distribution $p(\bm{m})$ in a much lower-dimensional space than the original $p(\bm{v})$, which can be easily sampled globally using the inverse transform sampling method. However, since this approach requires discretizing the $\bm{m}$-space for both training and generation, it becomes impractical for intrinsic space dimensions $d>4$ due to memory constraints. These low-rank RBMs are trained to replicate the statistics of the dataset projected onto its first $d$ principal components. Despite their simplicity, the low-rank models can still generate an approximate version of the dataset, as shown in Supplemental Figure 15 for the five datasets previously presented.

In the initial stage of the standard RBM learning process, the model encodes the strongest PCA components of the data through multiple critical transitions (Decelle et al., ; ; Bachtis et al., ). Pre-training with the low-rank construction allows us to bypass these transitions and avoid out-of-equilibrium effects caused by critical slowing down associated to these transitions, which becomes prohibitive in the case of very clustered datasets. 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 that at the transitions.

Figure: Figure 2: We compare the equilibrium samples generated by RBMs trained on the MNIST01, Mickey, HGD, Ising2D and CelebA datasets using two different training schemes: PCD, and PCD initialized on low-rank RBMs. To assess the fitting of the modes, we show a scatter plot of the projections of the samples in the first two principal directions of each dataset, together with the histograms along each direction. Black dots and histograms refer to the original data distribution, while the red color is used to represent the generated data. In D we show the evolution of the log-likelihood in time for the two training methods computed with the online AIS trajectory. In continuous lines we show the train and in dashed, the test. Results are obtained by averaging with respect to 10 trainings initialized at the same low-rank RBM, and the standard deviation between all is shown in a shade around the line, which is hardly seen due to the small variation from run to run.
Refer to caption: x2.png

Overcoming these transitions has dramatic implications for the quality of the models we can train and how accurately they reproduce the statistics of the data, especially for short training times.
In Fig. 2, we show for 5 datasets the equilibrium samples drawn from 2 RBMs trained with identical number of samples, minibatch size, $k=100$ Gibbs steps, and learning rate $\gamma=0.01$, but different training strategies. In particular, we consider an RBMs trained from scratch with the standard PCD (Tieleman, ), and a final machine trained with PCD and pre-trained with a low-rank RBM. We also considered machines trained with the recently proposed Jarzynski reweighing method (Carbone et al., ) (see SI E for our specific implementation in the RBM), but we do not include the analysis here because in most of the datasets the training method eventually becomes very unstable reaching very poor results. In all datasets, the quality of the generated samples and the test log-likelihood are significantly better when pre-training is used. We see that the pre-training+PCD training perfectly balances the different modes (here we show the projections on the first 2 PCA directions, but the agreement extends further to several other principal directions). We can also compare the log-likelihood of all 5 datasets and find that the pre-trained RBM achieves higher values. In all cases, the log-likelihood is measured using the annealed importance sampling trajectory (Tr-AIS) method, and the samples are drawn using the trajectory PT (PTT) algorithm. Both algorithms will be described in the next section 5.

Already from the scatter plots, we see that the pre-training is key in obtaining models where all modes are properly balanced in clustered datasets, but it’s also beneficial in the maximum test log-likelihood we can achieve. We observe that the pretraining is not particularly useful when the dataset is not much clustered, as it is the case of CelebA and the entire MNIST dataset as shown in Fig. 16 of the SI I. This result can be easily understood as the encoding of the first modes of the data do not involve splitting the distribution into two isolated modes, something that happens in the case of clustered data.

### 5 Trajectory annealing

The evolution of parameters during training can be seen as a smoother annealing process compared to the standard temperature annealing used in algorithms like annealing importance sampling (AIS) for log-likelihood estimation or the PT algorithm for sampling. This is because it avoids the first-order transitions that typically occur with temperature in clustered datasets. However, these methods can still be applied to trajectory annealing by using the parameters from different stages of training, rather than scaling the model parameters by a common factor $\beta_{i}$, as is done in temperature-based protocols.

#### 5.1 Trajectory Annealing importance sampling

The standard method for computing the log-likelihood (LL) in the literature is Annealed Importance Sampling (AIS) (Neal, ). For a detailed review of its application in RBMs, see Krause et al. (). In this approach, the partition function of a model with parameters $\bm{\Theta}=\{\bm{w},\bm{\theta},\bm{\eta}\}$ is estimated by gradually increasing the inverse temperature $\beta$ from 0 to 1 over $N_{\beta}$ steps. At $\beta=0$ (where $p_{0}=p_{\mathrm{ref}}$), the partition function $Z_{0}$ is tractable, enabling the estimation of $Z_{i}$ at neighboring $\beta_{i}$ values through path sampling methods. The reference configuration $p_{\mathrm{ref}}$ is commonly taken uniform, but improvements have been reported in RBM when fixing it to the independent site distribution that describes the empirical center of the dataset. For a detailed description of the algorithm, see SI D.

AIS is quite costly because it requires a large number of annealing steps $N_{\beta}$ for reliable estimations, making it impractical to compute the log-likelihood online during training. As a result, it is often calculated offline using a set of models at different temperatures. However, since the training process itself follows an annealing trajectory, it makes more sense to use the evolution of parameters during training as annealing steps rather than relying on temperature. Additionally, at the start of training, parameters are set to zero, making $Z_{0}$ tractable, which also holds when initializing with a low-rank RBM. Since updating the partition function only involves calculating the energy difference between the new and old parameters, this can be efficiently computed online during each parameter update, allowing for very small integration steps. We refer to this method as *online* Tr-AIS. Alternatively, similar to standard AIS, multiple models can be saved during training for an *offline* Tr-AIS estimation.

To evaluate the reliability of our new approach, we trained several RBMs using datasets with a maximum of 20 hidden nodes, allowing for the exact computation of the partition function through direct enumeration. We then calculated the difference between the true partition function and the estimations obtained from various methods. The results for the HGD dataset are presented in Fig. 3, while similar qualitative results for the other datasets can be found in the SI, see Figs. 11, 12, 13 and 14. In Fig. 3, we display the estimation errors in the LL using three different methods: standard temperature Annealed Importance Sampling (AIS) (panel A), AIS with a reference distribution (panel B), and Tr-AIS (panel C). For the temperature AIS estimates, we show results for different values of $N_{\beta}$ in distinct colors. In panel C, we compare the online estimation in black with the offline Tr-AIS estimation, using the same models analyzed in panels A and B. Additionally, we present the error in LL estimation using a few machines and the trajectory parallel tempering (PTT), which will be detailed in the next section. Overall, the online Tr-AIS estimation consistently outperforms the other methods, and only the PTT estimation yields similar results in offline analyses.

Figure: Figure 3: Comparison of the error in (LL) estimation (with respect to the exact value) using different methods: AIS (left), AIS with a reference distribution (middle), and Tr-AIS (ours, right). For Tr-AIS, we present three approaches to computing it: online during training, offline using the models saved during training, and with PTT. The RBM was pretrained and trained with a total of 20 hidden nodes (allowing for an exact computation of the LL) for 10,000 gradient steps. PTT uses a subset of saved models such that the acceptance rate between two consecutive models is 0.25.
Refer to caption: x3.png

#### 5.2 Standard Gibbs sampling vs. Parallel Trajectory Tempering

One major challenge with structured datasets is quantifying the model’s quality, since sampling the equilibrium measure of a well-trained model is often too time-consuming. This affects the reliability of generated samples and indirect measures as log-likelihood’s estimation through Annealing Importance Sampling (AIS) (Krause et al., ), making them inaccurate and meaningless.

Figure: Figure 4: Comparison of the performance of different MCMC sampling methods on RBMs trained with the pretraining+PCD procedure on the MNIST01 (row-1), the HGD (row-2) and Ising 2D datasets (row-3). In A and B columns, we show the trajectory of two independent Markov chains (red and orange) iterated for $10^{4}$ MCMC steps using either PTT or alternate Gibbs sampling (AGS), projected onto the first two principal components of the dataset. The position of the chains is plotted every 10 MCMC steps. The black contour represents the density profile of the dataset. In column C, we show the average number of jumps between the two regions separated by the dashed grew line in the PCA plots using different MCMC methods: Alternative Gibbs Sampling (AGS), Parallel Tempering (PT) (Hukushima & Nemoto, ), without or with reference configuration (PTref) (Krause et al., ), the Stacked Tempering (Roussel et al., ), and the Trajectory Parallel Tempering (PTT) proposed in this work. The average is calculated over a population of 1000 chains and the shadow around the lines indicates the error of the mean.
Refer to caption: x4.png

To illustrate this problem, let us consider the MNIST01, the HGD and the Ising 2D datasets. MNIST01 and Ising 2D are strongly bimodal in the PCA and the HGD trimodal, as shown in Fig. 1. Let us consider that we want to sample the equilibrium measure of the RBMs trained using low-rank RBM pretraining. In order to draw new samples from these models, one would typically run MCMC processes from random initialization and iterate them until convergence. The mixing time in such case is controlled by the jumping time between clusters. To accurately estimate the relative weight between modes, the MCMC processes must be ergodic, requiring many back-and-forth jumps. However, as shown in Figs. 4–A1-3, Alternate Gibbs sampling (AGS) dynamics are extremely slow, rarely producing jumps even after $10^{4}$ MCMC steps. The red curves in Figs. 4–C1-3 show the mean number of jumps over 1000 independent chains as a function of MCMC steps, indicating that a proper equilibrium generation would require at least $10^{6}-10^{7}$ MCMC steps. In order to speed up sampling, we consider other traditional optimized MCMC methods such as the PT (with or without reference), the Stack Tempering (Roussel et al., ) and our Parallel Trajectory Tempering (PTT) that will be explained below. The evolution of the number of jumps between modes observed with each of these algorithms, as a function of the total number of AGS steps performed (i.e. the number of parallel models multiplied by the number of AGS steps performed at each model) is shown in Fig. 4–C1-3. Before discussing our PTT algorithm and the results, let us briefly introduce the PT algorithm.

##### 5.2.1 Parallel Tempering approaches

Parallel tempering (PT) is commonly used to simulate multiple systems in parallel at different inverse temperatures $\beta$ (Hukushima & Nemoto, ; Marinari & Parisi, ), with occasional temperature swaps between models based on the Metropolis rule. Since this move satisfies detailed balance, the dynamics will converge to the Boltzmann distribution $p_{\beta}(\bm{x})\propto\exp(-\beta\mathcal{H}(\bm{x}))$ at each $\beta$. However, a more general approach involves simulating multiple systems, each with its own set of parameters. We can define a system as $p_{\bm{\Theta}^{t}}(\bm{x})$, where $\bm{\Theta}^{t}$ denotes the parameters of the model indexed by $t$. This includes the classic PT case, where temperature scaling is applied, $\bm{\Theta}^{t}=\beta^{(t)}\bm{\Theta}$, with $\bm{\Theta}$ being the parameters of the model to be sampled.

In this work, we introduce a new algorithm that exchanges parameters of models saved at different stages of learning, using the learning trajectory to define the parameter set $\bm{\Theta}^{t}$, where $t$ refers to a specific point in the training process. We call this method Parallel Trajectory Tempering (PTT). The rule for swapping configurations $\bm{x}=(\bm{v},\bm{h})$ between neighboring machines indexed by $t$ and $t-1$ is given by

$$ $p_{\mathrm{acc}}(\bm{x}^{t}\leftrightarrow\bm{x}^{t-1})=\min\left(1,\exp\left( \Delta\mathcal{H}^{t}(\bm{x}^{t})-\Delta\mathcal{H}^{t}(\bm{x}^{t-1})\right) \right),$ (5) $$

where again the index $t$ indicates that the Hamiltonian is evaluated using the parameters $\bm{\Theta}^{t}$. This move satisfies detailed balance with our target equilibrium distribution $p(\bm{x})=\exp(-\mathcal{H}(\bm{x}))/Z$, ensuring that the moves lead to the same equilibrium measure. The reason for using the trajectory is that PT algorithms work very well when a system undergoes a continuous phase transition. As discussed in the appendix G, adjusting the inverse temperature can lead to discontinuous transitions for which PT is not known to perform very well, while following the trajectory, the system undergoes several continuous transitions as discussed in Decelle et al. (); Bachtis et al. ().

Our PTT approach relies on the ease of sampling configurations from an initial RBM. In classical PCD training, the early “non-specialized” models thermalize quickly because they behave like Gaussian distributions. For RCM+PCD, independent configurations can be efficiently sampled from the low-rank RBM using a static Monte Carlo algorithm, as described in Section A.1 of the SI. The trajectory flow then significantly speeds up convergence to equilibrium, especially in the presence of a multimodal structure. The time interval between successive models is chosen to keep the acceptance probability for exchanges between neighboring models (Eq. 5) around 0.25 (we observed that the total time (number of AGS multiplied by number of machines) needed to thermalize a PTT chain did not depend on the acceptance fixed in the interval [0.2,0.5].
In practice, effective models can be saved during training by monitoring the acceptance rate between the last saved model and the current one. Pre-trained models require far fewer saved models, as many are naturally positioned near major phase transitions. The number of models used for each sampling process, along with a detailed description of the algorithm, is provided in the SI.

In Fig. 4–C1-3, we present the number of jumps between clusters as a function of the total number of AGS steps, where in PT-like schemes this refers to 1 (AGS step + one swap proposal) multiplied by the total number of models run in parallel. In all cases, PTT sampling outperforms the other methods, showing particularly remarkable results for HGD. Additionally, our method has the advantage of not requiring extra training like Stacked Tempering or simulating unnecessary temperatures.

Moreover, the equilibrium samples obtained from different models during sampling can be used to compute the log-likelihood, similar to the standard PT approach Krause et al. (). In Fig. 3–right we show that the quality of this estimation is comparable with the online trajectory AIS. Another desirable advantage of the PT algorithms is that one can easily check the thermalisation of the MCMC process by investigating the ergodicity of the random walk in models Banos et al. (). In addition, we evaluated other indicators of sample quality and potential overfitting. As detailed analysis is in Sect. H of the SI. Our method consistently provides high-quality samples across all datasets and, thanks to precise likelihood computation, can effectively detect overfitting.

### 6 Conclusions

In this work, we tackle the challenges of training high-quality equilibrium RBMs on highly structured datasets. We focused on improving the sampling process and accurately evaluating test log-likelihood throughout training. Our findings reveal that for datasets that are tightly clustered in their principal component analysis representation, introducing a pre-training phase can enhance training quality. This phase involves encoding the principal directions of the data into the model using a convex optimization approach. However, this pre-training method, that is specific to RBMs, appears to have limited benefits for less structured datasets or for long trainings, as chains get out of equilibrium, and should be combined in the future with optimized MCMC methods to achieve its maximum capacity.

Furthermore, we developed a new approach that exploits the progressive learning of features during RBM training to enhance evaluation and sampling. This approach enables online and precise computation of log-likelihood and leads to an efficient sampling strategy akin to the standard parallel tempering algorithm.
Our results demonstrate that both our proposed methods outperform traditional algorithms without the need for simulating models at unnecessary temperatures or training additional models to facilitate sampling. Furthermore, both the trajectory AIS method for log-likelihood estimation and the trajectory parallel tempering algorithm are highly versatile, making them applicable to more complex energy-based models. Additionally, the low-rank RBM model can serve as an effective pre-initialization step for deeper architectures, and the trajectory PT algorithm is suitable for direct use in any EBM, regardless of complexity.

### 7 Code availability

The code and datasets are available at
[https://github.com/anonymousrbm/iclr-2025](https://github.com/anonymousrbm/iclr-2025)

## Supplemental Information

### Appendix A Details of the pre-training of a low rank RBM

#### A.1 The low-rank RBM and its sampling procedure

Our goal is to pre-train an RBM to directly encode the first $d$ principal modes of the dataset in the model’s coupling matrix. This approach avoids the standard procedure of progressively encoding these modes through a series of second-order phase transitions, which negatively impact the quality of gradient estimates during standard training. It also helps prevent critical relaxation slowdown of MCMC dynamics in the presence of many separated clusters.

Given a dataset, we want to find a good set of model parameters ($\bm{w}$, $\bm{\theta}$ and $\bm{\eta}$) for which the statistics of the generated samples exactly match the statistics of the data projected onto the first $d$ directions of the PCA decomposition of the training set. Let us call each of these $\alpha=1,\ldots,d$ projections $m_{\alpha}=\bm{u}_{\alpha}\cdot\bm{v}/\sqrt{N_{\mathrm{v}}}$ the magnetizations along the mode $\alpha$, where $\bm{u}_{\alpha}$ is the $\alpha$-th mode of the PCA decomposition of the dataset. A simple way to encode these $d$-modes is to parameterize
the $w$-matrix as:

$$ $\bm{w}=\sum_{\alpha=1}^{d}w_{\alpha}\bar{\bm{u}}_{\alpha}\bm{u}_{\alpha}^{\top },\qquad\text{with}\qquad(\bm{u}_{\alpha},\bar{\bm{u}}_{\alpha})\in{\mathbb{R} }^{N_{v}}\times{\mathbb{R}}^{N_{h}},$ (6) $$

where $\bm{u}$ and $\hat{\bm{u}}$ are respectively the right-hand and left-hand singular vectors of $\bm{w}$, the former being directly given by the PCA, while
$w_{\alpha}$ are the singular values of $\bm{w}$. Using this decomposition, the marginal energy on the visible variables, $\mathcal{H}(\bm{v})=\log\sum_{\bm{h}}\exp\mathcal{H}(\bm{v},\bm{h})$ can be rewritten in terms of these magnetizations $\bm{m}\equiv(m_{1},\ldots,m_{d})$

$$ $\mathcal{H}(\bm{v})=-\sum_{a}\log\cosh\left(\sqrt{N_{\mathrm{v}}}\bar{u}_{a} \sum_{\alpha=1}^{d}w_{\alpha}m_{\alpha}+\eta_{a}\right)=\mathcal{H}(\bm{m}(\bm {v})).$ (7) $$

Now, the goal of our pre-training is not to match the entire statistics of the data set, but only the marginal probability of these magnetizations. In other words, we want to model the marginal distribution

$$ $p_{\mathrm{emp}}(\bm{m})\equiv\sum_{\bm{v}}p_{\mathrm{emp}}(\bm{v})\prod_{ \alpha=1}^{d}\delta\left(m_{\alpha}-\frac{1}{\sqrt{N_{\mathrm{v}}}}\bm{u}_{ \alpha}^{T}\bm{v}\right),$ (8) $$

where $\delta$ is the Dirac $\delta$-distristribution. In this formulation, the distribution of the model over the magnetization $\bm{m}$ can be easily characterized

$$ $\displaystyle p(\bm{m})$ $\displaystyle=\frac{1}{Z}\sum_{\bm{v}}e^{-\mathcal{H}(\bm{v})}\prod_{\alpha=1} ^{d}\delta\left(m_{\alpha}-\frac{1}{\sqrt{N_{\mathrm{v}}}}\bm{u}_{\alpha}^{T} \bm{v}\right)$ (9) $\displaystyle=\frac{1}{Z}\mathcal{N}(\bm{m})\exp{\sum_{a}\log\cosh\left(\bar{u }_{a}\sum_{\alpha=1}^{d}w_{\alpha}m_{\alpha}+\eta_{a}\right)}$ (10) $\displaystyle=\frac{1}{Z}e^{-\mathcal{H}(\bm{m})+N_{\mathrm{v}}s(\bm{m})}= \frac{1}{Z}e^{-N_{\mathrm{v}}f(\bm{m})}$ (11) $$

where $\mathcal{N}(\bm{m})=\sum_{\bm{v}}\prod_{\alpha=1}^{d}\delta\left(m_{\alpha}-
\frac{1}{\sqrt{N_{\mathrm{v}}}}\bm{u}_{\alpha}^{T}\bm{v}\right)$ is the number of configurations with magnetizations $\bm{m}$, and thus $S(\bm{m})=\log N(\bm{m})/N_{\mathrm{v}}$ is the associated entropy. Now, for large $N_{\mathrm{v}}$ the entropic term can be determined using large deviation theory, and in particular the Gärtner-Ellis theorem:

$$ $p_{\mathrm{prior}}(\bm{m})=\frac{e^{N_{\mathrm{v}}s(\bm{m})}}{2^{N_{\mathrm{v} }}}\approx\exp\left(-N_{\mathrm{v}}\mathcal{I}(\bm{m})\right),$ (12) $$

with the rate function

$$ $\mathcal{I}(\bm{m})=\sup_{\bm{\mu}}\left[\bm{m}^{T}\bm{\mu}-\phi(\bm{\mu}) \right]=\bm{m}^{T}\bm{\mu}^{*}-\phi(\bm{\mu}^{*}),$ (13) $$

and

$$ $\displaystyle\phi(\bm{\mu})$ $\displaystyle=\lim_{N_{\mathrm{v}}\to\infty}\frac{1}{N_{\mathrm{v}}}\log\left \langle e^{N_{\mathrm{v}}\bm{m}^{T}\bm{\mu}}\right\rangle=\lim_{N_{\mathrm{v}} \to\infty}\frac{1}{N_{\mathrm{v}}}\log\frac{1}{2^{N_{\mathrm{v}}}}\sum_{\bm{v} }e^{\sqrt{N_{\mathrm{v}}}\sum_{\alpha=1}^{d}\mu_{\alpha}\sum_{i}u_{\alpha,i}v_ {i}}$ (14) $\displaystyle=\lim_{N_{\mathrm{v}}\to\infty}\frac{1}{N_{\mathrm{v}}}\sum_{i=1} ^{N_{\mathrm{v}}}\log\cosh\left(\sqrt{N_{\mathrm{v}}}\sum_{\alpha=1}^{d}\mu_{ \alpha}u_{\alpha,i}\right).$ (15) $$

Then, given a magnetization $\bm{m}$, we can compute the minimizer $\bm{\mu}^{*}(\bm{m})$ of $\phi(\mu)-\bm{m}^{T}\mu$ which is convex, using e.g. Newton method which converge really fast since we are in small dimension. Note that in practice we will obviously use finite estimates of $\phi$, assuming $N_{v}$ is large enough.
As a result we get $\bm{\mu}^{*}(\bm{m})$ satisfying implicit equations given by the constraints given at given $N_{v}$:

$$ $m_{\alpha}=\frac{1}{\sqrt{{N_{\mathrm{v}}}}}\sum_{i=1}^{{N_{\mathrm{v}}}}u_{i} ^{\alpha}\tanh\left({\sqrt{{N_{\mathrm{v}}}}\sum_{\beta=1}^{d}u_{i}^{\beta}\mu ^{*}_{\beta}}\right).$ (16) $$

It is then straightforward to check that spins distributed as

$$ $p_{\mathrm{prior}}(\bm{v}|\bm{m})\propto e^{{N_{\mathrm{v}}}{\bm{\mu}^{*}}^{T} \bm{m}(\bm{v})}$ (17) $$

fulfill well the requirement, as $\left\langle\bm{u}_{\alpha}^{T}\bm{v}/\sqrt{{N_{\mathrm{v}}}}\right\rangle_{p_
{\mathrm{prior}}}=m_{\alpha}$. In other words, we can generate samples having mean magnetization $m_{\alpha}$ just by choosing $v_{i}$ as

$$ $p_{\mathrm{prior}}(v_{i}=1|\bm{m})=\text{sigmoid}\left({2\sqrt{{N_{\mathrm{v}} }}\sum_{\alpha=1}^{d}u_{\alpha,i}\mu^{*}_{\alpha}(\bm{m})}\right)$ (18) $$

The training can therefore be done directly in the subspace of dimension $d$. In Ref. Decelle & Furtlehner (), it has been shown that such RBM can be trained by mean of the Restricted Coulomb Machine, where the gradient is actually convex in the parameter’s space. It is then possible to do a mapping from the RCM to the RBM to recover the RBM’s parameters. In brief, the training of the low-dimensional RBM is performed by the RCM, and then the parameters are obtrained via a direct relation between the RCM and the RBM’s parameters. The detail of the definition and of the training of the RCM is detailed in the appendix A.2.

#### A.2 The Restricted Coulomb Machine

As introduced in Decelle & Furtlehner (), it is possible to exactly train a surrogate model for the RBM, called the Restricted Coulomb Machine (RCM), on a low dimensional dataset without explicitly sampling the machine allowing to learn even heavily clustered datasets.

The RCM is an approximation of the marginal distribution of the RBM with $\{-1,1\}$ binary variables:

$$ $\mathcal{H}(\bm{v})=-\sum_{i}v_{i}\theta_{i}-\sum_{a}\log\cosh\left(\sum_{i}w_ {ia}v_{i}+\eta_{a}\right).$ (19) $$

We then project both the parameters and variables of the RBM on the first $d$ principal components of the dataset:

$$ $m_{\alpha}\coloneqq\frac{1}{\sqrt{N_{v}}}\sum_{i=1}^{N_{v}}v_{i}u_{i\alpha}, \quad w_{\alpha a}\coloneqq\sum_{i=1}^{N_{v}}w_{ia}u_{i\alpha},\quad\theta_{ \alpha}\coloneqq\frac{1}{\sqrt{N_{v}}}\sum_{i=1}^{N_{v}}\theta_{i}u_{i\alpha}$ (20) $$

with $\alpha\in\{0,\ldots,d\}$ and $\bm{u}$ the projection matrix of the PCA.

The projected distribution of the model is then given by

$$ $p_{\text{RBM}}(\bm{m})=\frac{\exp\left(N_{v}\left[\mathcal{S}(\bm{m})+\sum_{ \alpha=1}^{d}\theta_{\alpha}m_{\alpha}+\frac{1}{N_{v}}\sum_{a=1}^{N_{h}}\log \cosh\left(\sqrt{N_{v}}\sum_{\alpha=1}^{d}m_{\alpha}w_{\alpha a}+\eta_{a} \right)\right]\right)}{Z}$ (21) $$

where we ignore the fluctuations related to the transverse directions and $\mathcal{S}[\bm{m}]$ accounts for the non-uniform prior on $\bm{m}$ due to the projection of the uniform prior on $\bm{s}$ for the way to compute it.

The RCM is then built by approximating

$$ $\log\cosh(x)\simeq|x|-\log 2,$ (22) $$

which is valid for $x$ large enough. The probability of the RCM is thus given by:

$$ $p_{\text{RCM}}(\bm{m})=\frac{\exp\left(N_{v}\left[\mathcal{S}(\bm{m})+\sum_{ \alpha=1}^{d}\theta_{\alpha}m_{\alpha}+\sum_{a=1}^{N_{h}}q_{a}\left|\sum_{ \alpha=1}^{d}n_{\alpha}m_{\alpha}+z_{a}\right|\right]\right)}{Z}$ (23) $$

where

$$ $q_{a}=\sqrt{N_{v}\sum_{\alpha=1}^{d}w_{\alpha a}^{2}},\quad n_{a}=\frac{w_{ \alpha a}}{\sqrt{\sum_{\alpha=1}^{d}w_{\alpha a}^{2}}},\quad z_{a}=\frac{\eta_ {a}}{\sqrt{N_{v}\sum_{\alpha=1}^{d}w_{\alpha a}^{2}}}.$ (24) $$

This can be easily inverted as

$$ $w_{\alpha a}=\frac{1}{\sqrt{N_{v}}}q_{a}n_{a}\qquad\text{and}\qquad\eta_{a}=q_ {a}z_{a},$ $$

in order to obtain the RBM from the RCM.
The model is then trained through log-likelihood maximization over its parameters. However, this objective is non-convex if all the parameters are trained through gradient ascent.
To relax the problem, since we’re in low dimension, we can define a family of hyperplanes $(\bm{n},\bm{z})$ covering the space and let the model only learn the weights of each to the hyperplane. We can then discard the ones with a weight low enough for the approximation equation 22 to be bad.

The gradients are given by

$$ $\frac{\partial J(\bm{\Theta})}{\partial q_{a}}=\mathbb{E}_{\bm{m}\sim p_{ \mathcal{D}}(\bm{m})}\left[|\bm{n}_{a}^{T}\bm{m}+z_{a}|\right]-\mathbb{E}_{\bm {m}\sim p_{\text{RCM}}(\bm{m})}\left[|\bm{n}_{a}^{T}\bm{m}+z_{a}|\right],$ (25) $$

$$ $\frac{\partial J(\bm{\Theta})}{\partial\theta_{\alpha}}=E_{\bm{m}\sim p_{ \mathcal{D}}(\bm{m})}\left[m_{\alpha}\right]-\mathbb{E}_{\bm{m}\sim p_{\text{ RCM}}(\bm{m})}\left[m_{\alpha}\right].$ (26) $$

The positive term is straightforward to compute. For the negative term, we rely on a discretization of the longitudinal space to estimate the probability density of the model and compute the averages.

#### A.3 Learning the bias

To initialize the model with its bias set to the dataset’s mean, we define the 0 0-th direction of the longitudinal space as the direction obtained by projecting the configuration space onto the normalized mean vector of the dataset.

The RBM will only contribute to this direction through the visible bias:

$$ $H(\bm{m}(v))=-\sum_{\alpha=0}^{N_{v}}\theta_{\alpha}m_{\alpha}-\frac{1}{N_{v}} \sum_{a}\log\cosh\left(\sqrt{N_{v}}\sum_{\alpha=1}^{d}m_{\alpha}w_{\alpha a}+ \eta_{a}\right).$ (27) $$

This implies that this direction is independent of the others within the model. Consequently, we can create an independent mesh for this direction to learn it. The rest of the training is performed using the same procedure as before.

### Appendix B Sampling via Parallel Tempering using the learning trajectory

Figure: Figure 5: Scheme of PTT. We Initialize the chains of the models by starting from a configuration $\bm{x}_{0}^{(0)}$ and passing it through the machines along the training trajectory, each time performing $\tilde{k}$ mcmc steps. For pre-train+PCD, $\bm{x}_{0}^{(0)}$ is a sampling from the RCM, otherwise it is a uniform random initialization. The sampling consists of alternating one mcmc step for each model with a swap attempt between adjacent machines. For pre-train+PCD, at each step we sample a new independent configuration for $\mathrm{RBM}_{0}$ using the RCM.
Refer to caption: x5.png

Assuming we have successfully trained a robust equilibrium model, there remains the challenge of efficiently generating equilibrium configurations from this model. Although models trained at equilibrium exhibit faster and more ergodic dynamics compared to poorly trained models, the sampling time can still be excessively long when navigating a highly rugged data landscape. Consequently, we devised a novel method for sampling equilibrium configurations that draws inspiration from the well-established parallel tempering approach. In this traditional method, multiple simulations are conducted in parallel at various temperatures, and configurations are exchanged among them using the Metropolis rule. Unlike this conventional technique, our method involves simultaneously simulating different models that are selected from various points along the training trajectory. This approach is motivated by the perspective that learning represents an annealing process for the model, encountering second-order type phase transitions during training. In contrast, annealing related to temperature changes involves first-order phase transitions, making traditional parallel tempering less effective for sampling from clustered multimodal distributions.

A sketch of the Parallel Trajectory Tempering (PTT) is represented in fig. 5. Specifically, we save $t_{\mathrm{f}}$ models at checkpoints $t=1,\dots,t_{\mathrm{f}}$ along the training trajectory. We denote the Hamiltonian of the model at checkpoint $t$ as $\mathcal{H}_{t}$, and refer to the Hamiltonian of the RCM model as $\mathcal{H}_{0}$. We define $\mathrm{GibbsSampling}(\mathcal{H},\bm{x},k)$ as the operation of performing $k$ Gibbs sampling updates using the model $\mathcal{H}$ starting from the state $\bm{x}$. In all our sampling simulations we used $k=1$.

The first step is to initialize the models’ configurations efficiently. This involves sampling $N$ chains from the RCM model, $\bm{x}_{0}^{(0)}\sim\mathrm{RCMSampling}(\mathcal{H}^{0})$, and then passing the chains through all the models from $t=1$ to $t=t_{\mathrm{f}}$, performing $k$ Gibbs steps at each stage: $\bm{x}_{t}^{(0)}\sim\mathrm{GibbsSampling}(\mathcal{H}_{t},\bm{x}_{t-1}^{(0)},k)$.

The sampling process proceeds in steps where we update the configuration of each model except $\mathcal{H}_{0}$ with $k$ Gibbs steps, and sample a completely new configuration for the RCM model $\mathcal{H}_{0}$. Following this update step, we propose swapping chains between adjacent models with an acceptance probability given by:

$$ $p_{\mathrm{acc}}(\bm{x}_{t}\leftrightarrow\bm{x}_{t-1})=\min\left(1,\exp\left( \Delta\mathcal{H}_{t}(\bm{x}_{t})-\Delta\mathcal{H}_{t}(\bm{x}_{t-1})\right) \right),$ (28) $$

where $\Delta\mathcal{H}_{t}(\bm{x})=\mathcal{H}_{t}(\bm{x})-\mathcal{H}_{t-1}(\bm{x})$.

We continue alternating between the update step and the swap step until a total of $N_{\mathrm{mcmc}}$ steps is reached.
The sampling procedure is illustrated in the following pseudo-code:

$$ $p_{\mathrm{acc}}(\bm{x}_{t}^{(n)}\leftrightarrow\bm{x}_{t-1}^{(n)})=\min\left( 1,\exp\left(\Delta\mathcal{H}_{t}(\bm{x}_{t}^{(n)})-\Delta\mathcal{H}_{t}(\bm{ x}_{t-1}^{(n)})\right)\right)$ $$

A comparison of performances between PTT and standard Gibbs sampling is reported in Fig.4.

#### B.1 Controlling equilibration times with PTT

One of the advantages of PT-like approaches is that they allow easy control over thermalization times by studying the random walk in temperatures (or models for the PTT) visited by each run over time, as proposed in Ref. Banos et al. (). From now on, we will refer only to the number of models used in the PTT sampling, denoted as $N_{\mathrm{m}}$. In an ergodic sampling process, each run should spend roughly the same amount of time at each of the $N_{\mathrm{m}}$ different model indices, as all models have the same probability.

The mean temperature/model index is known to be $\left\langle n\right\rangle=(N_{\mathrm{m}}-1)/2$ (assuming the model indices run from 0 0 to $N_{\mathrm{m}}-1$). This simplifies the computation of the time-autocorrelation function of the model indices $n_{i}(t)$ visited by each sample $i$ at time $t$:

$$ $C(t)=\frac{\overline{(n_{i}(t+t_{0})-\left\langle n\right\rangle)(n_{i}(t_{0}) -\left\langle n\right\rangle)}}{\overline{(n_{i}(t_{0})-\left\langle n\right \rangle)^{2}}}$ $$

Here, the symbol $\overline{(\cdots)}$ refers to the average over all samples $i$ and initial times $t_{0}$ obtained from our PTT runs. This allows us to compute the exponential autocorrelation time, $\tau_{\mathrm{exp}}$, from a fit to

$$ $C(t)\sim Ae^{-t/\tau_{\mathrm{exp}}}\quad\text{for large }t,$ $$

as well as the integrated autocorrelation time, $\tau_{\mathrm{int}}$, using the self-consistent equation Sokal (); Amit & Martin-Mayor ():

$$ $\tau_{\mathrm{int}}=\frac{1}{2}+\sum_{t=0}^{6\tau_{\mathrm{int}}}C(t).$ $$

These two times provide different insights: $\tau_{\mathrm{exp}}$ gives insight into the time needed to thermalize the simulations (setting the length of the runs to at least $20\tau_{\mathrm{exp}}$ is generally a safe (and very conservative) choice to ensure thermalization Sokal ()), while $\tau_{\mathrm{int}}$ gives the time needed to generate independent samples. To keep track of the fluctuations, we compute a different $C^{m}(t)$ for each $m=1,\ldots,N_{\mathrm{m}}$ system index simulated in parallel in the PTT, where each $C^{m}(t)$ is averaged over 1000 independent realizations (with all runs initialized to the same model index). We then compute the average over all parallel runs $C(t)=\overline{C^{m}(t)}$ and use this curve to extract $\tau_{\mathrm{exp}}$ and $\tau_{\mathrm{int}}$. The errors in $C$ and the times are then extracted from the error of the mean of the values extracted with each of the $C^{m}(t)$ curves.

For example, in Fig. 6 (left), we show the averaged curves $C(t)$ obtained with the PTT sampling of the machines discussed in the main text in the MNIST01 (top) and HGD (bottom) datasets. The extracted times are $\tau_{\mathrm{exp}}=258(29)$ and $\tau_{\mathrm{int}}=9.0(1.5)$ for MNIST01, and $\tau_{\mathrm{exp}}=48(10)$ and $\tau_{\mathrm{int}}=3.8(7)$
for HGD. These numbers guarantee us that the $10^{4}$ steps used to generate the scatterplots in Fig. 2 were long enough to ensure equilibration in this dataset, but also in the others, as shown in Table 1.

Figure: Figure 6: (Left) We show the autocorrelation function $C(t)$ obtained as in Eq. equation B.1 when we analyze the random walk in the model indices during PTT sampling with the same RBMs analyzed in the main text (trained with pretraining+PCD) for the MNIST01 (top) and HGD (bottom) datasets. Both curves were obtained with a PTT run of $5\cdot 10^{4}$ MCMC steps, where the models used were automatically selected during training using the PTT moves acceptance criterion of 0.25. (Right) In blue we show the $\tau_{\mathrm{exp}}$ as a function of the target acceptance rate for MNIST01 (top) and HGD (bottom). The error bar results from the error of the mean of the values extracted from the different $C_{m}(t)$ and the shading marks the range between the largest and the smallest value of the $N_{\mathrm{m}}$ curves. In orange, we show the total simulation costs ($\tau_{\mathrm{exp}}$ multiplied by the number of models to be simulated in parallel) as a function of the target acceptance rate. This number of models saved, as a function of the acceptance rate, is shown in green in the inset.
Refer to caption: extracted/5987275/fig/fcorr_PTT_MNIST01.png

**Table 1: Autocorrelation times associated to the random walk dynamics between models in the PTT sampling obtained with the models trained with pre-training+PCD on different datasets. Models used in the PTT are all chosen with the $a=0.25$ criterion.**
| dataset | $\tau_{\mathrm{exp}}$ | $\tau_{\mathrm{int}}$ |
| --- | --- | --- |
| MNIST01 | 258(29) | 9.0(15) |
| HGD | 48(10) | 3.8(7) |
| Mickey | 2(0) | 1(0) |
| Ising | 18(3) | 2.4(5) |

With these measures, we can try to optimize the number of models that we use for the PTT and that are automatically selected during training. The criterion here is that the acceptance of PTT exchanges between neighboring models must not go below a target acceptance $a$. In the figures shown so far, $a=0.25$. Let us justify this choice. The total physical time we need to simulate $N_{\mathrm{m}}$ models grows proportionally with $N_{\mathrm{m}}$ (unless one manages to properly parallelize the sampling, which is feasible but we did not do it), while the time needed to flow between very well trained and poorly trained models controlling the mixing time scales as $\sqrt{N_{\mathrm{m}}}/a$, where $a$ is the average acceptance of the exchanges between the models. Since $N_{\mathrm{m}}$ grows with $a$ in practice (as shown in the inset of Fig. 6–right), it is expected that the total time can be optimized for a given value of $a$. In Fig. 6–right we show $\tau_{\mathrm{exp}}$ in blue (and $\tau_{\mathrm{exp}}N_{\mathrm{m}}$ in orange) as a function of $a$ for the data sets MNIST01 (top) and HGD (bottom).
These figures show that the choice of $a$ is not particularly crucial as long as it is between 0.2 and 0.5.

### Appendix C Training details

We describe in Tables 2 and 3 the datasets and hyperparameters used during training. The test set was used to evaluate the metrics.

**Table 2: Details of the datasets used during training.**
| Name | #Samples | #Dimensions | Train size | Test size |
| --- | --- | --- | --- | --- |
| CelebA | 30 000 | 1024 | 60% | 40% |
| Human Genome Dataset (HGD) | 4500 | 805 | 60% | 40% |
| Ising | 20 000 | 64 | 60% | 40% |
| Mickey | 16 000 | 1000 | 60% | 40% |
| MNIST-01 | 10 610 | 784 | 60% | 40% |
| MNIST | 50 000 | 784 | 60% | 40% |

**Table 3: Hyperparameters used for the training of RBMs.**
| Name | Batch size | #Chains | #Epochs | Learning rate | #MCMC steps | #Hidden nodes |
| --- | --- | --- | --- | --- | --- | --- |
| CelebA |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 500 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 500 |
| HGD |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Ising |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Mickey |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| MNIST-01 |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| MNIST |  |  |  |  |  |  |
| PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |
| Pre-train+PCD | 2000 | 2000 | 10 000 | 0.01 | 100 | 200 |

**Table 4: Training time of the RBMs**
| Name | Training time (s) |
| --- | --- |
| CelebA |  |
| PCD | $119.2\pm 0.8$ |
| Pre-train+PCD | $3134.0\pm 0.6$ |
| HGD |  |
| PCD | $106.2\pm 0.3$ |
| Pre-train+PCD | $266\pm 3$ |
| Ising |  |
| PCD | $62\pm 1$ |
| Pre-train+PCD | $104\pm 2$ |
| Mickey |  |
| PCD | $118.3\pm 3$ |
| Pre-train+PCD | $246.0\pm 0.7$ |
| MNIST-01 |  |
| PCD | $100.6\pm 0.3$ |
| Pre-train+PCD | $258.0\pm 0.9$ |
| MNIST |  |
| PCD | $105.7\pm 0.2$ |
| Pre-train+PCD | $233.9\pm 0.2$ |

**Table 5: RCM Hyperparameters**
| Dataset name | # Dimensions | # Hidden nodes | Training time (s) |
| --- | --- | --- | --- |
| CelebA | 4 | 100 | 3006 |
| HGD | 3 | 100 | 150 |
| Ising | 3 | 100 | 32 |
| Mickey | 2 | 100 | 118 |
| MNIST-01 | 3 | 100 | 148 |
| MNIST | 3 | 100 | 119 |

All experiments were run on a RTX 4090 with an AMD Ryzen 9 5950X.

### Appendix D Computation of the likelihood using AIS

To compute the likelihood of the RBMs we used the Annealed Importance Sampling method introduced by Neal () and recently studied in the case of the RBM Krause et al. ().

#### D.1 Offline AIS

We follow the method described in Krause et al. (). The AIS estimate is based on creating a sequence of machines, at one end a machine from which we can compute the partition function and at the other end the machine from which we want to compute the partition function. Let’s consider that we have $i=0,\dots,N_{\rm rbm}$ machines described by

$$ $p_{i}(\bm{x})=\frac{e^{-\mathcal{H}_{i}(\bm{x})}}{Z}=\frac{p_{i}^{*}(\bm{x})}{ Z},$ $$

where $p_{0}=p_{\rm ref}$ denote the RBM from which we can compute the partition function (typically a RBM where the weight matrix is zero) and $p_{N_{\rm rbm}}=p_{\rm target}$ the machine from which we want to estimate the partition function. For each machine, we denote its set of parameters $\bm{\Theta}_{i}=\{\bm{w}_{i},\bm{\theta}_{i},\bm{\eta}_{i}\}$. Following Krause et al. (), defining a transition rule $T_{i}(\bm{x}_{i-1}\to\bm{x}_{i})$ in order to pass from a configuration $\bm{x}_{i-1}$ of the system $i-1$ to $\bm{x}_{i}$ of the system $i$, that we have the following identity

$$ $\frac{Z_{N_{\rm rbm}}}{Z_{0}}=\left\langle\exp\left(-\sum_{i=1}^{N_{\rm rbm}} \mathcal{H}_{i+1}(\bm{x}_{i})-\mathcal{H}_{i}(\bm{x}_{i})\right)\right\rangle_ {\rm ForwardPath}.$ (29) $$

There, if one can compute $Z_{0}$, the estimate of $Z_{N_{\rm rbm}}$ is obtained by averaging the r.h.s. of Eq. 29 over the forward trajectory starting from a configuration $\bm{x}_{0}$ using the operators $T_{i}$. While the method is exact, the variance of the estimation relies on on close are two successive distribution $p_{i}$, $p_{i+1}$.

###### Temperature AIS:

a commonly used implementation of AIS, is to multiply the Hamiltonian of the system by a temperature, and to the different systems indexed by $i$, by a set of $\beta_{i}$ where $\beta_{0}=0$, $\beta_{i}<\beta_{i+1}$ and $\beta_{N_{\rm rbm}}=1$. Therefore, it is equivalent to choosing $\Theta_{i}=\beta_{i}\{\bm{w},\bm{\theta},\bm{\eta}\}$. Another version of Temperature AIS consist in choosing a slightly different reference distribution $p_{0}$ that take into account the bias of the model. In this setting, the intermediate distribution are taken as

$$ $p_{i}(\bm{x})=\frac{p_{\rm ref}^{(1-\beta_{i})}e^{-\beta_{i}\mathcal{H}(\bm{x} )}}{Z_{i}}.$ $$

The main difference is that a prior distribution is used as reference $p_{0}=p_{\rm ref}$ and the temperatures interpolate between this prior and the target distribution. It is quite common to take

$$ $p_{\rm ref}(\bm{x})=\frac{\exp\left(\sum_{i}\theta_{i}v_{i}+\sum_{a}\eta_{a}h_ {a}\right)}{\prod_{i}2\cosh(\theta_{i})\prod_{a}2\cosh(\eta_{a})}.$ $$

###### Trajectory AIS:

the general derivation of AIS allow in principle to choose the set forward probability and/or intermediate system as one sees fit. Since again, a change of temperature in the system can induce discontinuous transition for which, successive distribution will be far away. While, when considering system that are neighbors in the learning trajectory, the system should pass through continuous transition which do not suffer this problem. Hence, we consider a set of RBMs along the learning trajectory. The RBMs are chosen following the prescription that the exchange rate between configuration (see Eq. 5) is about $\sim 0.3$). Then the formula Eq. 29 is used to estimate the partition function.

#### D.2 Online AIS

The online AIS (Annealed Importance Sampling) is the trajectory AIS evaluated at every gradient update during training. A set of permanent chains is initialized at the beginning of the training process. Every time the parameters are updated, a step is performed on these chains before using them to compute the current importance weights. This approach is computationally efficient relative to the cost of sampling (1 step vs. 100 steps) and allows for a very precise estimation.

### Appendix E Training of the RBM using the Jarzynski equation

In this section, we describe the procedure introduced in Carbone et al. () for training energy-based models by leveraging the Jarzynki equation, and we adapt it to the specific case of the RBM. At variance with Carbone et al. (), we introduce an additional assumption on the transition probability of the parameters that will lead us to the the reweighting implemented in population annealing methods Weigel et al. ().

In one of its formulations, the Jarzynski equation states that we can relate the ensemble average of an observable $\mathcal{O}$ with the average obtained through many repetitions of an out-of-equilibrium dynamical process. Let us consider a path $\boldsymbol{y}=\{(\boldsymbol{x}^{0},\boldsymbol{\theta}^{0}),(\boldsymbol{x}^
{1},\boldsymbol{\theta}^{1}),\dots,(\boldsymbol{x}^{t},\boldsymbol{\theta}^{t})\}$, where $\bm{x}^{i}=(\bm{v}^{i},\bm{h}^{t})$ refers to the model variables and $\bm{\theta}^{i}=(w^{i},\eta^{i},\theta^{i})$, the model parameters, both at time $i$.

If we consider the training trajectory of an RBM, $p_{0}\rightarrow p_{1}\rightarrow\dots\rightarrow p_{t-1}\rightarrow p_{t}$, we can write the average of an observable $\mathcal{O}$ over the last model of the trajectory $p_{t}$ as

$$ $\langle\mathcal{O}\rangle_{t}=\frac{\left\langle\mathcal{O}~{}e^{-W^{t}}\right \rangle_{\mathrm{traj}}}{\left\langle e^{-W^{t}}\right\rangle_{\mathrm{traj}}} =\frac{\sum_{r=1}^{R}\mathcal{O}(\boldsymbol{y}_{r})~{}e^{-W_{r}^{t}}}{\sum_{r =1}^{R}e^{-W_{r}^{t}}},$ (30) $$

where $W^{t}$ is a trajectory-dependent importance factor and the averages on the rhs are taken across many different trajectory realizations.
In particular, $W_{t}$ is given by the recursive formula:

$$ $\displaystyle W^{t}(\bm{y})$ $\displaystyle=$ $\displaystyle\mathcal{H}^{t}(\bm{x}^{t})-\mathcal{H}^{0}(\bm{x}^{0})-\sum_{i=1 }^{t}\log\tilde{T}^{i}(\bm{y}^{i}\rightarrow\bm{y}^{i-1})+\sum_{i=1}^{t}\log T ^{i}(\bm{y}^{i-1}\rightarrow\bm{y}^{i})$ (31) $\displaystyle=$ $\displaystyle W^{t-1}(\bm{y})+\mathcal{H}^{t}(\bm{x}^{t})-\mathcal{H}^{t-1}( \bm{x}^{t-1})-\log\frac{\tilde{T}^{t}(\bm{y}^{t}\rightarrow\bm{y}^{t-1})}{T^{t }(\bm{y}^{t-1}\rightarrow\bm{y}^{t})}$ (32) $$

Where $T$ and $\tilde{T}$ are, respectively, the forward and backward transition probabilities.

The transition matrix between two states $\bm{y}^{t-1}\rightarrow\bm{y}^{t}$ is factorized into two different subsequent movements:

- 1.
Update the variables $\bm{x}^{t-1}\to\bm{x}^{t}$ on the fixed model $\bm{\theta}^{t-1}$,
- 2.
Update the parameters $\bm{\theta}^{t-1}\to\bm{\theta}^{t}$ on fixed variables $\bm{x}^{t}$.

In other words,

$$ $T(\bm{y}^{t-1}\to\bm{y}^{t})=w_{\bm{\theta}^{t-1}}(\bm{x}^{t-1}\to\bm{x}^{t})t _{\bm{x}^{t}}(\bm{\theta}^{t-1}\to\bm{\theta}^{t}),$ (33) $$

and

$$ $\tilde{T}(\bm{y}^{t}\to\bm{y}^{t-1})=t_{\bm{x}^{t}}(\bm{\theta}^{t}\to\bm{ \theta}^{t-1})w_{\bm{\theta}^{t-1}}(\bm{x}^{t}\to\bm{x}^{t-1}),$ (34) $$

so that

$$ $\log\frac{\tilde{T}^{t}(\bm{y}^{t}\rightarrow\bm{y}^{t-1})}{T^{t}(\bm{y}^{t-1} \rightarrow\bm{y}^{t})}=\log\frac{w_{\bm{\theta}^{t-1}}(\bm{x}^{t}\to\bm{x}^{t -1})}{w_{\bm{\theta}^{t-1}}(\bm{x}^{t-1}\to\bm{x}^{t})}+\log\frac{t_{\bm{x}^{t }}(\bm{\theta}^{t}\to\bm{\theta}^{t-1})}{t_{\bm{x}^{t}}(\bm{\theta}^{t-1}\to \bm{\theta}^{t})}.$ (35) $$

We know the ratio between the transition probabilities between $\bm{x}^{t-1}\to\bm{x}^{t}$, because evolution of the variables during the simulation is performed using MCMC simulations that satisfy detailed balance. This means that:

$$ $\log\frac{w_{\bm{\theta}^{t-1}}(\bm{x}^{t}\to\bm{x}^{t-1})}{w_{\bm{\theta}^{t- 1}}(\bm{x}^{t-1}\to\bm{x}^{t})}=\mathcal{H}^{t-1}(\bm{x}^{t})-\mathcal{H}^{t-1 }(\bm{x}^{t-1}),$ (36) $$

and we are left with

$$ $\displaystyle W^{t}(\bm{y})$ $\displaystyle=$ $\displaystyle W^{t-1}(\bm{x})+\mathcal{H}^{t}(\bm{x}^{t})-\mathcal{H}^{t-1}( \bm{x}^{t})-\log\frac{t_{\bm{x}^{t}}(\bm{\theta}^{t}\to\bm{\theta}^{t-1})}{t_{ \bm{x}^{t}}(\bm{\theta}^{t-1}\to\bm{\theta}^{t})}.$ (37) $$

The standard formula that we use in the paper assumes that the changes in $\bm{\theta}$ are so small that we can consider the transition probabilities to be symmetrical and the last term in the rhs cancels out. In that case, we recover the known expressions of the population annealing.

Notice that, since Eq. equation 30 is an exact result, the importance weights should, in principle, eliminate the bias brought by the non-convergent chains used for approximating the log-likelihood gradient in the classical PCD scheme. However, after many updates of the importance weights, one finds that only a few chains carry almost all the importance mass. In other words, the vast majority of the chains we are simulating are statistically irrelevant, and we expect to get large fluctuations in the estimate of the gradient because of the small effective number of chains contributing to the statistical average. A good observable for monitoring this effect is the Effective Sample Size (ESS), defined as Carbone et al. ()

$$ $\mathrm{ESS}=\frac{\left(R^{-1}\sum_{r=1}^{R}e^{-W^{(r)}}\right)^{2}}{R^{-1} \sum_{r=1}^{R}e^{-2W^{(r)}}}\in[0,1],$ (38) $$

which measures the relative dispersion of the weights distribution. Then, a way of circumventing the weight concentration on a few chains is to resample the chain population according to the importance weights every time the ESS drops below a certain threshold, for instance 0.5. After this resampling, all the chain weights have to be set to 1 ($W^{(r)}=0\,\,\forall r=1,\dots,R$).

### Appendix F Tuning of the hyperparameters

In order to justify our choice of hyperparameters in the main text, we illustrate on Fig. 7 the value of the log-likelihood as a function of the number of hidden nodes and of the learning of the machine. In order to be fair, the comparison is done at a constant number of gradient updates of $10000$. In this regime, the machines that perform the best have a learning rate of $\gamma\sim 0.1,0.005$ and about $N_{h}=500,1000$ hidden nodes. We can expect that smaller learning rates reach the same value of the log-likelihood but for a much longer training. We clearly see the benefit of the RCM, where the log-likelihood achieves much higher values in general.

Figure: Figure 7: Comparison of the end log-likelihood reached by the RBM using different pairs of values for the learning rate and the hidden layer dimension on MNIST-01. The left panel shows these result when performing PCD training and the right one the results when pre-training the RBM before PCD.
Refer to caption: x6.png

### Appendix G Continuous phase transition

In statistical physics, and more precisely in mean-field model, a model described by the Gibbs-Boltzmann distribution is studied in the infinite size limit, $N\to\infty$, and characterize the state of the system by investigating the free energy of the system, defined by minus the log partition function

$$ $F[\bm{\theta}]=-\log(Z[\bm{\theta}])$ $$

One of the goals is thus to determine the typical configurations that should be expected given a set of model parameters. In the case of the RBM, the question is what are the typical visible (and/or hidden) variables that one expects to observe. For example, one may ask what typical energies one expects to observe for the parameters $\bm{\theta}$. For this purpose, let us rewrite the partition function as a sum over all possible values of the energy

$$ $Z=\sum_{E}e^{-E}g(E)=\sum_{E}e^{-E+S(E)}=\sum_{E}e^{-F(E)}=\sum_{E}e^{-Nf(E)}$ $$

where $g(E)$ the number of states with a given energy $E$, and thus $S(E)=\log g(E)$ is the entropy. Finally, the free energy is defined as $F(E)=E-S(E)$ and $f(E)=F(E)/N$ is the free energy per variable. All $Z$, $E$ and $F$ depend on the model parameters. $N$ refers to the number of variables.

For $N\to\infty$, the sum in Z is then dominated by the states with minimum free energy. A “first order”, or discontinuous transition, occurs in the parameters if the first derivative of $f(E)$ is discontinuous with respect to the parameters (in physics, the parameter is typically the temperature). A “second order”, or continuous transition, occurs when the first derivative is continuous, but the second derivative is discontinuous.

In the case of the RBM, the first learning transition can be mapped to a very simple model, the so-called Curie-Weiss model, as recently described in Ref. Bachtis et al. (). The Curie-Weiss model is fully solvable, and for this reason the meaning of a continuous transition is perhaps easier to understand there. Let’s give it a try. The Curie-Weiss model consists in a set of $N$ discrete variables $s_{i}=\pm 1$. We define the (exponential) distribution over these variables as

$$ $p(s)=\frac{1}{Z}\exp(\frac{\beta}{N}\sum_{i<j}s_{i}s_{j}+\beta H\sum_{i}s_{i})$ $$

where $\beta$ and $H$ are parameters of the model and $Z$ the normalization constant also called partition function. The question is to understand the “structure” of the distribution as a function of $\beta$.

###### Case H = 0 𝐻 0 H=0 italic_H = 0 :

we expect that for $\beta$ small each spin behaves as an isolated Rademacher random variables (each spin being $\pm 1$ with probability one-half), while for large value of $\beta$, the distribution of the system is dominated by configurations where the variables have (almost) all the same signs. The key mathematical aspect is to study the system in the limit where $N\to\infty$. A way to study this distribution is to investigate the moment generating function. It is possible to show (rigorously in this precise model), that the probability of $m=N^{-1}\sum_{i}s_{i}$, also called the magnetization of the system, is given by $p(m)\propto\exp(-N\Omega(m))$ (in the large $N$ limit) where $\Omega(m)$ is a large deviation function. In the CW model, it is given by

$$ $\Omega(m)=\frac{\beta m^{2}}{2}-\log\left[2\cosh(\beta m)\right].$ $$

The structure of $\Omega(m)$ is such that for $\beta<\beta_{c}=1$, it is a convex function with only one minimum in $m=0$, while for $\beta>\beta_{c}$ it has two symmetric minima $m=\pm m(\beta)\neq 0$ with $\Omega(m(\beta))=\Omega(-m(\beta))$ (the two minimum are equally probable). The values (in function of $\beta$) of the $m$ are continuous: $m=0$ for $\beta\in[0,1]$, and then positive until it saturates to $1$ as $\beta\to\infty$. In practice, it means that the distribution $p(m)$ pass from an unimodal distribution to a bimodal one as $\beta$ increases. The point where it happens is at $\beta_{c}=1$ and is called a second order (continuous) phase transition because de value of $m$ changes continuously from zero to non-zero values. A particular interesting properties is that at $\beta_{c}$, the function $\Omega(m)$ develops a non-analyticity in its second-order derivative, that can be linked with long-range fluctuations.
For a recent review Kochmański et al. (), and for more rigorous results Sinai ().

###### Case H ≠ 0 𝐻 0 H\neq 0 italic_H ≠ 0 :

in this case we expect that the system will always be biased toward configurations where the spins have the same direction(sign) as the bias (or field) $H$. Yet, for small value of $\beta$, the bias will be quite small, hence the magnetization will be small but not zero, $m\sim\beta H$. When $\beta$ is large, the system will now have a stronger and strong magnetization toward the magnetic field $m\sim{\rm sign}(H)$. Now we can again compute the large deviation function of the system

$$ $\Omega(m)=\frac{\beta m^{2}}{2}-\log\left[2\cosh(\beta m)\right]+Hm.$ $$

Again, it is possible to analyze the function $\Omega$ as a function of $\beta$ and $H$. For $\beta<1$ the function is convex and there is only one minimum located at $m>0$. When $\beta>1$ two situations can be distinguished. First, let’s denote

$$ $H_{\rm sp}(\beta)=\sqrt{\frac{\beta-1}{\beta}}-{\rm atanh}\left(\frac{\beta-1} {\beta}\right).$ $$

If $H>H_{\rm sp}(\beta)$, then the large deviation function has only one minimum for positive $m$. Now if $H<H_{\rm sp}(\beta)$, then $\Omega$ will present two minima. A first one for positive $m$ and a second one for negative $m$. The reason is that the pairwise interactions “manage” to create a metastable state for negative magnetization even in the presence of a positive field or bias. It is clear that this state will be unstable w.r.t. the one with positive magnetization which traduces in a minimum of higher value in the large deviation function. By changing the value of $H$ from large positive values toward large negative ones, we then pass from a system with only one state, then two states with one being subdominant, and finally back to one state but where the magnetization has changed its sign. In the infinite size limit, the system will always dominated by the state where $\Omega$ is at its lowest value. However, if one prepare a system close to a metastable state, the system will remain trap there until it disappears.

Figure: Figure 8: Left: The case $H=0$ for various temperature. We illustrate the symmetric shape of the large deviation function and how the two minima emerges from the centers. The minima of $m$ go continuously from zero to non-zero values. In the inset we also illustrate the shape of the probability distribution for moderately large size $N=500$. For $N\to\infty$ they converge toward a $\delta$ distribution. Right: The case $H\neq 0$. In this example, we illustrate what happen where the bias parameter $H$ is changed at a fix temperature $\beta=1.15>1$. We see that for small values of $H$, $\Omega$ has two minima, yet the dominant one suppresses exponentially the other one for large system size. For larger values of the bias the metastable minimum disappears. We show how the place of the minimum jump when the bias goes from negative to positive values.
Refer to caption: x8.png

###### How is this related to the training of the RBM ?

First, it has been shown in previous works Decelle et al. (); Bachtis et al. () that the RBM’s training undergoes a set of continuous transitions, at least at the beginning of the learning. Hence, changing the parameters of the RBM following the learning trajectories, the transitions are continuous.
Then in Decelle & Furtlehner () it is shown that the relaxed version of the RBM, namely the
RCM undergoes a first-order transition when changing the temperature. Schematically
to understand that, it is enough to see how the parameters $\beta$ acts on the system. In the case of the RCM, the energy term behaves in the standard way, i.e. linearly with the inverse temperature $\beta$ and the free energy takes the ordinary form

$$ $\log(Z)=-\beta F=-\beta E+S.$ $$

Therefore $\beta$ act as a trade-off parameters between the energy $E$ and the entropy $S$. By varying $\beta$, we change the contribution of the energy while keeping the entropy fixed. Hence, when the corresponding large deviation function of the RCM has reached a point having many minima with the same value of the free energy up to $~{}\mathcal{O}(1/N)$ corrections, all these minima realize in principle different trade-off between energy and entropy. As a result a slight change in $\beta$ will offset this sensible balance and the
respective free energies will acquire differences of order ${\mathcal{O}}(N)$ leading to favor dramatically one state against the others, i.e. the one having originally the lowest [resp. the largest] energy when $\beta$ is increased. It is very likely that by changing $\beta$ a discontinuous transition is encountered, where the minima of $\Omega$ vanishes far away from each other rather than a continuous one.
Concerning the RBM the argument for the occurence of such a phase transition, is only slightly more involved than for Coulomb machine because for RBM the energy (obtained after summing over latent variables) has a non-linear behavior w.r.t. $\beta$. A formal argument would go schematically as follows: just writing the local minima of the free energy (function of magnetization in the reduced intrinsic space) and looking at how the equilibrium is displaced with temperature shows that their energy vary indeed in a non-linear way with temperature and non-uniformly regarding their position on the
magnetization manifold, while the entropy contribution again do not change. This again necessarily implies that changing even slightly the temperature will
break the highly sensitive equilibrium obtained between these state corresponding to the multimodal distribution and will again typically favor one state among all.
This scenario is expected to occur as soon as the data are located on low dimensional space ($d=O(1)$ ) compared to a large embedding space ($d=O(N),\ N\gg 1$)
which is quite common in our point of view, at least for the type of data we are interested in.

To conclude, moving the parameters along the learning trajectory, leads the system to undergo continuous transition where $p_{\bm{\Theta(t)}}(\bm{s})$ is very close to $p_{\bm{\Theta}(t+\delta t)}(\bm{s})$. Changing the temperature tends to produces discontinuous transition where $p_{\beta(t)\bm{\Theta}}(\bm{s})$ can happen to be very far from $p_{\beta(t+\delta t)\bm{\Theta}}(\bm{s})$.

### Appendix H Overfitting and privacy loss as quality indicators

In this section, we examine the quality of the samples generated, regarding overfitting and privacy criteria which have been defined for genomic data in particular.
We look at this on the models trained with PCD with and without pre-training.
We focus on the human genome dataset, as shown in Fig. 1–C, to evaluate the ability of various state-of-the-art generative models to generate realistic fake genomes while minimizing privacy concerns (i.e., reducing overfitting). Recent studies (Yelmen et al., ; ) have thoroughly investigated this for a variety of generative models. Both studies concluded that the RBM was the most accurate method for generating high-quality and private synthetic genomes. The comparison between models relies primarily on the Nearest Neighbor Adversarial Accuracy
($AA_{\rm TS}$) and privacy loss indicators, introduced in Ref. Yale et al. (), which quantify the similarity and the level of ”privacy” of the data generated by a model w.r.t. the training set. We have
$AA_{\rm TS}=\frac{1}{2}\bigl{(}AA_{\rm True}+AA_{\rm Synth}\bigr{)}$
where $AA_{\rm True}$ [resp. $AA_{\rm Synth}$] are two quantities in $[0,1]$ obtained by merging two sets of real and synthetic data of equal size $N_{s}$ and measuring respectively the frequency that a real [rep. synthetic] has a synthetic [resp. real] as nearest neighbor.
If the generated samples are statistically indistinguishable from real samples, both frequencies $AA_{\rm True}$ and $AA_{\rm Synth}$ should converge to 0.5 at large ${N_{\mathrm{s}}}$.
$AA_{\rm TS}$ can be evaluated both with train or test samples and the privacy loss indicator is defined as ${\rm Privacy\ loss}=AA_{\rm TS}^{\rm test}-AA_{\rm TS}^{\rm train}$ and is expected to be strictly positive.
Fig. 10 shows the comparison of $AA_{\rm TS}$ and privacy loss values obtained with our two models, demonstrating that the pre-trained RBM clearly outperforms the other model, and even achieves better results ($AA_{\rm TS}$ values much closer to 0.5) than those discussed in Yelmen et al. (; ).

Figure: Figure 9: Evolution of the AATS (first row) and privacy loss (second row) indicator during the training of pretrained RBMs. The online estimation of the train and test LL is given at the third row. The overfitting is clearly detected in the HGD and MICKEY dataset and the divergence between the train and test LL coincides with a sharp rise in the privacy loss.
Refer to caption: x10.png

Figure: Figure 10: Comparison of the evolution of overfitting metrics along the log-likelihood between a pretrained RBM (A) and a RBM trained from scratch using PCD (B). (C) shows the comparison between the AASyn and AATruth between the two machines. The machines were selected in order to have a comparable privacy loss (D), right before starting to overfit.
Refer to caption: x13.png

Figure: Figure 11: Comparison of log-likelihood estimation w.r.t. the exact ll on the MNIST-01 dataset.
Refer to caption: x14.png

Figure: Figure 12: Comparison of log-likelihood estimation w.r.t. the exact ll on the MNIST dataset.
Refer to caption: x15.png

Figure: Figure 13: Comparison of log-likelihood estimation w.r.t. the exact ll on the Ising dataset.
Refer to caption: x16.png

Figure: Figure 14: Comparison of log-likelihood estimation w.r.t. the exact ll on the MICKEY dataset.
Refer to caption: x17.png

### Appendix I Extra figures

Figure: Figure 15: Samples generated with the low-rank RBM trained with 3 directions plus a bias. We show the samples generated by the pre-trained machined for the 5 datasets discussed in the paper and projected on the PCA of the dataset. For the image datasets, we also show the images generated.
Refer to caption: x18.png

Figure: Figure 16: We show the evolution of the tr-AIS log-likelihood in the training of the entire MNIST dataset (containing digits from 0 to 9) with and without the pre-training. The pre-training does not suppose a particular advantage but neither is a disadvantage.
Refer to caption: x19.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: *** introduction

Whilst the RBM is well known, it would be helpful I feel for a reader to have the definition of the model earlier in the text. It currently isn't defined until near the end of page 4. Please introduce the RBM formally earlier in the text.

Notation: inconsistent use of $N_v$ and $N_{\text{v}}$ throughout, similarly for $N_h$.

Equation 1: it might be better to write W_{i\alpha}, rather than w_{i\alpha} since w is used later for the "singular values".

*** page 2

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.

*** page 5

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

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.

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.

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.

It's unclear to me what contribution the authors are claiming to make. They state that the learning of the low rank parameterisation of W has been done before. Please clarify what the contributions of the paper are.


*** Section 5

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.

*** Section 6

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.


*** Supplementary material

The use of the term "mode" isn't very clear. The phrasing suggests that the first d modes of the maximum likelihood trained RBM should correspond to the d "modes" of the PCA solution. I'm not sure I know what this means. What are modes of a PCA solution?

The notation \hat{u} is confused with \bar{u}.

Why use $w$ here whereas $W$ is used in the main text?

The derivation is quite confusing. For example the dependence on \bar{u} in equation 7 disappears without explanation. Indeed \bar{u} seems to be never properly defined.

Please state clearly what are the parameters of the model that are being learned.

Section A.2. The claim as before of exact training is incorrectly made here.

The notation in equation 20 is confusing, such as w_{\alpha,a}=\sum_i w_{ia}u_{i\alpha} -- are arabic and latin indices meant to indicate referencing a different entity, even though both objects are labelled w?

In general I find the supplementary material confusing. I believe it is trying to fit an RBM projected to the d-dimensional subspace defined by PCA of the data to the empirical data distribution in that same subspace. However, approximations are clearly required in order to compute the projected RBM distribution. Given that, for a very low dimension d then one can easily discretise the model and carry out a simple maximum likelihood fit. If that is what is being done, it is not well explained and rather misleading (since this requires approximations itself).

An alternative (and standard) way to compute the marginal p(m) is to use the integral (Fourier) representation of the Dirac delta function. This means that the summation over v can be then carried out exactly, leaving only a d-dimensional integral to exactly compute p(m). This can also be carried out using discretisation for small d. The authors are (as I can understand) also using discretised integrals, so I'm unclear why they don't employ the standard Fourier Delta representation approach to compute p(m) -- this would seem to involve less approximations that the approach the authors consider.
soundness: 2
strengths: RBMs are an important model and finding appropriate ways to train them is a topic of significant interest. The paper highlights the phenomenon of critical slowing down and how pre-training the model with a low-rank approximation of the parameter matrix can help the model overcome some of the slowing down effects.
summary: The paper discusses approximations to train a restricted Boltzmann machine (RBM). The first is to pre-train the RBM by fitting a constrained (low-rank) form of the RBM to the low-dimensional PCA space of the data. This can help with finding a good initial solution. After this various MCMC approaches are considered to continue training.
weaknesses: The paper suffers from a lack of clarity of presentation and lack of clarity of novelty.

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.

In terms of presentation, there are notational inconsistencies and a general lack of clarity in terms of the main ideas. 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.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: ## Questions, Small comments and typos
- Would it be possible for the authors to provide a sketch and pseudocode for their PTT algorithm as a standalone and in comparison to PT? This would be very helpful to get a better understanding of the contribution of this work. 
- Is there any intuition behind the bump observed in Figure 3 at around $10^3$ gradient updates (left and middle plot).
- Layout: there's a problem with Figure 2. The x-axis is sometimes completely or partly cut. I strongly recommend carefully checking this, aligning the plots, and making sure such problems are removed. 
- In general the authors often refer to the Appendix as SI (I assume Supplement Information). I guess this acronym has not been defined anywhere. I identify its first occurrence in line 96. Perhaps the authors can define what SI is or, alternatively just all it appendix. 
- Line 235: I'd recommend adding a reference for critical slowing down. This comment applies to earlier occurrences of this concept.
- Line 459: grew -> grey
- Line 512: Banos et al. (2010) might need to be wrapped in parenthesis \cite -> \citep
soundness: 3
strengths: - The paper is well-written and easy to follow. 
- It represents a pleasant read that is accessible to a broad audience. 
- The literature review and related work section read well and are exhaustive.
- The idea of pre-training the RBM to encode the principal components is simple yet very effective. 
- Leveraging the analogy between critical slowing down and the struggle of RBM during training to be ergodic and discovering all modes of the distributions is elegant and intuitive (though I suppose this is not a novelty of this paper, it is very nicely pictured in the introduction). 
- The numerical experiments look solid and aligned with the theoretical insights given in the main text. 
- I have not thoroughly checked the mathematical details in the appendix, but at first glance, they look good.
summary: The claimed novelties of this work are twofold. 
First, this paper proposes low-ranking training of RBMs by directly encoding the principal components throughout a convex-optimization process. This pre-training component proves to be very efficient when data are particularly clustered. In such cases, target densities are highly multimodal, and the model struggles to "discover"all the modes from scratch during training without the pre-training phase. This autonomous discovery of new modes is often associated with second-order phase transitions, similar to systems from statistical mechanics, where critical slowing down prevents the discovery of all modes in finite time efficiently. 

As a second contribution, the paper also investigates how to use a variation of parallel tempering (PT) algorithms, termed parallel trajectory tempering, to sample more efficiently and obtain log-likelihoods estimates. In simple terms, parallel trajectory tempering (PTT) essentially relies on the same idea of parallel tempering of swapping between models at different temperatures using the Metropolis rule (and therefore retaining detailed balance). However, differently from PT, PTT swaps a full set of parameters $\Theta^t$ instead of the temperature $\beta$ only. In that sense, it can be thought of as a generalization of PT. 

Numerical experiments in Fig. 2 prove the pre-trained low-rank RBM to be more capable of identifying all modes in highly clustered data, while Figs. 3-4 show that PTT allows more accurate loglikelihood estimation and faster yet more efficient sampling from all modes of distribution compared to standard alternate Gibbs sampling (AGS).
weaknesses: - 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.  
- 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. 
- 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 $\boldsymbol{u}_0$, 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?

This latter point is crucial in assessing the effective novelty of this work. At the moment the reason for the lower score is primarily due to my perception of limited novelty. I am more than happy to discuss this with the authors during the rebuttal and revisit my score upon clarification of my concerns above (and below, see, e.g., the first bullet point in the **Questions** cell).

---

confidence: 5
contribution: 3
initial_rating: 3
presentation: 3
questions: - Does critical slowing down occur in the energy-based model when the hidden variables are traced out, or does it occur in the joint distribution that includes the hidden variables? If the phase transition occurs in the joint measure, does the traced-out distribution also exhibit a phase transition?
- What is the definition of $\bar{u}$?
- Could the authors provide a detailed derivation of Equation (4)? The terms $\bar{u}_{a}$ and $\eta_{a}$ are currently undefined.
- The phrase "a direction $\bar{u}_0$ is used for the magnetization $m_0$ that is only present in the bias term" is unclear. Could you explain this in more detail?
- Is it possible to learn DBM without pre-training using the pre-training with weights introduced by [1] ?

[1] Yuma Ichikawa and Koji Hukushima, Statistical-mechanical Study of Deep Boltzmann Machine Given Weight Parameters after Training by Singular Value Decomposition.
soundness: 2
strengths: - The paper offers a novel contribution by proposing a pre-training technique and a new sampling approach for RBMs inspired by their thermodynamic properties. This builds on the existing theoretical analyses of RBMs.
- To my knowledge, extending replica Monte Carlo methods to a learning trajectory is original and intriguing.
- Including a specialized physics background in the Appendix makes the paper accessible even to readers without a physics background.
summary: This research proposes an efficient training approach for structured data in RBMs by employing pre-training based on simple convex optimization, which significantly facilitates learning for structured datasets. Furthermore, the study introduces a novel sampling and log-likelihood evaluation method that leverages the model's learning process, differing from conventional Parallel Tempering.
weaknesses: 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. 

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

- 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.
- If a first-order transition does exist, then the exchange probability in PT would approach zero near the transition. Has this phenomenon been observed? Additionally, it would be helpful to evaluate the round-trip rate of PT and PTT.
- While it is argued that preparing models at different temperatures is challenging for PT, it should be noted that the proposed approach also requires storing models during the learning process.
- The CelebA data in Figure 2 appears to be truncated.

Because the high performance has been verified numerically, the score can be raised if the above statement is cleared.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Why didn't you apply this to larger problems?

2. What are situations where the pre-training fails?

3. Is PTT useful for generating samples during training, using only earlier parts of the training run?
soundness: 2
strengths: 1. The idea of low-rank pre-training is interesting and seems like it could be useful if it scaled up.

2. The idea of doing AIS across the training run is creative and clever.

3. Parallel tempering across training steps seems new.
summary: This paper studies algorithms for the training of Restricted Boltzmann Machines (RBMs).  It argues that "highly structured" data require different algorithms than those that have been successful for, e.g., image datasets.  There are three algorithmic ideas that are discussed: 1) Pre-training an RBM using an "exact" procedure that produces low-rank weight matrices; 2) Estimating log-likelihoods using annealed importance sampling across steps of a training run; and 3) Using parallel tempering for sampling, again using different steps of training. Evidence for the efficacy of these procedures is provided in the form of curves from training runs on a few small datasets.
weaknesses: 1. I think this paper has a somewhat limited audience.  It mostly builds upon work from a small group of authors, using language most familiar to that community.  (For example, one person's work is cited thirteen times in the references.)  A significant amount of jargon is used that keeps this from being a readable stand-alone paper.  This is coupled with heuristic explanations for things that appear to rely on sharing the particular statistical mechanical point of view of this subcommunity.

2. Much of the motivation for the work centers on "highly structured" data, which is not defined clearly.  The authors indicate that this corresponds to the existence of clusters.  The paper does not show examples of the methods succeeding or failing in the presence of this structure.  For example, the Celeb-A dataset is given as an example of a dataset in which there are not clusters and so it is not "highly structured".  However, Figure 2 does not seem to show us that this matters for the pre-training procedure.  Figure 15 is
similar.  Why does one conclude that the bottom row of Fig 2 and Fig 15 are significantly different from what we see in the top row of Fig 2?

3. The main text is highly verbose, with most of the actual concrete content being in the appendices.  I don't think anything novel is introduced until page six.

4. I find it difficult to appreciate precisely what the contribution of Section 4 is.  As I understand it, the insight is "do Decelle & Furtlehner (2021a) before you do PCD".  This is useful information, but between this section and Appendix A, I'm not sure where the boundary is between this and D&F (2021a).

5. While the ideas of section 5 are interesting and Figure 3 is intriguing, the empirical results are at the level of "preliminary findings" on a single small problem.  Even with the vastly smaller compute resources of 15 years ago, RBM researchers were studying larger problems.

6. The title is too broad relative to what the paper delivers.

Typos:
 - L161-162: "two slow"
 - L478: "exchanges parameters" but I think you mean "exchanges configuration".
 - L775-776: \bar{u} vs \hat{u}.
 - L836-837: "gradient is convex" -- surely you mean the training objective is convex in the parameters.

Meta-review: The paper introduces a pretraining strategy for RBMs that enables better coverage of all the modes of a target density and more accurate partition function estimation. The paper also proposes a new sampling algorithm that outperforms existing MCMC algorithms. Experimental results demonstrate the validity of these claims.

Decision: Accept (Poster)