Fast training and sampling of Restricted Boltzmann Machines
Authors: Nicolas BEREUX
Aurélien Decelle
Cyril Furtlehner
Lorenzo Rosset
Beatriz Seoane
Venue: ICLR 2025 Poster (2025)
Keywords: Restricted Boltzmann Machine
Fast Sampling
structured data learning
training algorithm
Primary area: generative models
TL;DR: We design a new training pipeline for the Restricted Boltzmann Machine, capable of training quickly and accurately on hard multimodal dataset.

Abstract: Restricted Boltzmann Machines (RBMs) are powerful tools for modeling complex systems and extracting insights from data, but their training is hindered by the slow mixing of Markov Chain Monte Carlo (MCMC) processes, especially with highly structured datasets. In this study, we build on recent theoretical advances in RBM training and focus on the stepwise encoding of data patterns into singular vectors of the coupling matrix, significantly reducing the cost of generating new samples and evaluating the quality of the model, as well as the training cost in highly clustered datasets.  The learning process is analogous to the thermodynamic continuous phase transitions observed in ferromagnetic models, where new modes in the probability measure emerge in a continuous manner. 
We leverage the continuous transitions in the training process to define a smooth annealing trajectory that enables reliable and computationally efficient log-likelihood estimates. This approach enables online assessment during training and introduces a novel sampling strategy called Parallel Trajectory Tempering (PTT) that outperforms previously optimized MCMC methods.
To mitigate the critical slowdown effect in the early stages of training, we propose a pre-training phase. In this phase, the principal components are encoded into a low-rank RBM through a convex optimization process, facilitating efficient static Monte Carlo sampling and accurate computation of the partition function.
Our results demonstrate that this pre-training strategy allows RBMs to efficiently handle highly structured datasets where conventional methods fail. Additionally, our log-likelihood estimation outperforms computationally intensive approaches in controlled scenarios, while the PTT algorithm significantly accelerates MCMC processes compared to conventional methods.

Published as a conference paper at ICLR 2025
FAST TRAINING AND SAMPLING OF RESTRICTED BOLTZMANN MACHINES
Nicolas Béreux INRIA-Saclay, LISN Paris-Saclas University
Aurélien Decelle
Escuela Técnica Superior de Ingenieros Industriales Universidad Politécnica de Madrid
Departamento de Física Teórica
Universidad Complutense de Madrid
Cyril Furtlehner INRIA-Saclay, LISN Paris-Sae Paris-Saclay University
Lorenzo Rosset
LCQB, Sorbonne Université,
Laboratoire de Physique Théorique
École Normale Supérieure, Paris
Beatriz Scoane
Departamento de Física Teórica Universidad Complutense de Madrid
ABSTRACT
Restricted Boltzmann Machines (RBMs) are powerful tools for modeling complex systems and extracting insights from data, but their training is hindered by the slow mixing of Markov Chain Monte Carlo (MCMC) processes, especially with highly structured datasets. In this study, build recent theoretical advances in RBM training and focus on the stepwise encoding of data patterns into singular vectors of the coupling matrix, significantly reducing the cost of generating new samples and evaluating the quality of the model, as well as the training cost in highly clustered datasets. The learning process analogous to the thermodynamic continuous phase transitions observed in ferromagnetic models, where new modes in the probability measure emerge continuous manner. We leverage the continuous transitions the training process to define smooth annealing trajectory that enables reliable and computationally efficient log ligeelikhooo estimates. This approach enables online assessment during training and introduces novel sampling strategy called Parallel Trajectory Tempering (PTT) that outperforms previously optimized MCMC methods. To mitigate the critical slowdown effect in the early stages of training, we propose a pre-training phase. In this phase, the principal components are encoded into low-rank RBM through a convex optimization process, facilitating efficient static Monte Carlo sampling and accurate computation of the partition function. Our results demonstrate that this pre-training strategy allows RBMs to efficiently handle highly structured datasets where conventional methods fail. Additionally, our log-likelihood estimation outperforms computationally intensive approaches in controlled scenarios, while the PTT algorithm significantly accelerates MCMC processes compared to conventional methods.
INTRODUCTION
Energy-based models (EBMs) are long-established approach to generative modeling, with the Boltzmann Machine (BM) Hinton & Sejnowskil 1983) and Restricted Boltzmann Machine (RBM) (Smolensky) 1986 Ackley et al.) 1985) 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. (2077) have been widely used for decades for inference applications in fields like neuroscience (Hertz al. 2011 and computational biology (Cocco et al. 2018).

---

Published as a conference paper at ICLR 2025
More recent research has shown that RBMs can also be reinterpreted many-body physical mod- els (Decelle et al. 2024; Bulso & Roudi 2021), analyzed to identify hierarchical relationships in the data (Decelle et al. 2002) or used to uncover biologically relevant constituent features (Tubiana et al., 2019; Fernandez-e Cossio-Diaz, 2025). In essence, RBMs present a great compromise between expressibility and interpretabilityy: they universal approximators capable of modeling complex datasets and simple enough to allow a direct extraction of insights from data. This interpretative task particularly challenging when dealing with more powerful models, such as generative convolutional network EBMs (Xie et al. 2016; SSng & Kingma 2021) 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 are not abundant, which is often the case in genomics/proteomiis or neural recording datasets. For example, RBMs have been shown to be among the most reliable approaches to generate synthetic human genomes (Yelmen et al., 2021; 2023) or to model the immune response (Bravi et al. 2021a 2023 or neural activity (van der Plas et al. 2023 Quiroz Monnens et al. 2024).
RBMs, like other EBMs, are challenging to train due to the difficulty of computing the log-likelihood gradient, which requires ergodic exploration of complex free-energy landscape via Markov Chain Monte Carlo (MCMC). Training with non- non-convergent MCMC introduces out-of of---ufliibbiiu memory effects, as shown in recent studies (Nijkamp al. 2019; 2020 Decelle al. (202T) which can be analytically explained using moment-matching arguments (Nijkamp et al., 2019] Agoritsas et al.
)023). While these effects enable fast and accurate generative models for structured data (Carbone) et al. 2024a) and high-quality images with RBMs <Liao et al. (022)), they create sharp disconnect between the model's Gibbs-Boltzmann and dataset distributions, undermining parameter interpretabil-
ity (Agoritsas et al. 2023 Decelle et al. 2024). For meaningful insights, RBMs require proper chain mixing during training, emphasizing the need for equilibrium models.
Working with equilibrium RBMs poses three primary challenges: and sampling from the trained model. Both training and sampling rely on MCMC methods to draw equilibrium samples, but their efficiency depends heavily on the dataset's structure. While RBMs can be effectively trained on image datasets like MNIST or CIFAR- 10 with sufficient MCMC steps, this approach falters for highly structured datasets (Béreux et al. 2023), such as genomics, neural recordings, or low-temperature physical systems. These datasets often exhibit multimodal distributions with distinct clusters, as revealed by PCA (Fig. [D. which hinder mixing. Clustering is prominent in the first four datasets studied here, in contrast to the more compact PCA pattern seen in CelebA full MNIST. Dataset details are provided in Fig. [] Multimodal distributions exacerbate sampling challenges since mixing times depend on the MCMC chains' ability to transition between clusters. In RBMs, these distant modes emerge during training through continuous phase transitions (c.g., second-order or critical transitions) that encode patterns (Decelle et al. 2017 2018 Bachtis al. 2024). These transitions lead to critical slowdown (Hohenberg & Halperin 1977). where mixing times diverge with system size as the transition point nears. Critical slowdown along the training trajectory was recently analyzed in (Bachtis et al., 20024). For brief introduction to phase transitions and associated dynamical effects, see Supplemental Information (SD All this hampers both the training process and the model's ability to generate independent samples, as achieving adequate mixing requires prohibitively long simulations. Model quality evaluation adds to the challenge, particularly in unsupervised learning. A common qualitative approach compares equilibrium samples to data via PCA projections. Quantitatively, moment comparisons between data and generated distributions offer insights but fail to detect mode collapse. Log-likelihood computation provides single-score metric and tracks overfitting during training, but exact computation is infeasible due to the intractable partition function. Instead, methods such as Annealed Importance Sampling (AIS) (Neal 2001) are used for estimation but become unreliable for highly multimodal distributions and are computationally impractical for online evaluations (see Section 4.1p.
training, evaluating model quality,
In this paper, we build on recent advances in understanding the evolution of the landscape during the training of RBM models to propose 3 innovative approaches for tackling the challenges of evaluating. generating samples and training RBMs on highly structured datasets:
1. We propose a novel method for estimating log-likelihood (LL) by leveraging the softness of the training trajectory. This approach enables efficient online LL computation during training at mini- mal cost or offline reconstruction using trajectory annealing or parallel tempering. Moreover, we demonstrate that this method is not only less computationally expensive than existing techniques but also more accurate in controlled experiments. The method is introduced in Sect. 4.1
2

---

Published as a conference paper at ICLR 2025
A MNISTO1
C Human Genome
D - lsing 2D
- CelebA
Mickey
B-0.44
32x32
Ezropon
South Asio
Eas Astas
Auerican
Ailoan
Juf
0.S-
0:25-
0,00
0.0
0,00
2,26
0.25
-.50
022
6.00
0:25
0.0 0.6
PCJ
PC]
PCT
PCJ
PC1
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 I 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 '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 3 = 0.44. In E, the CelebA dataset in black and white, resized 32x32 pixels. For more details on these datasets, please refer to the SI.
2. We propose variation of the standard parallel tempering algorithm, where parameters from different stages of the learning process are exchanged instead of temperatures. This new sampling approach significantly outperforms previous optimized methods and is detailed in Section 4.2
3. We introduce a pretraining strategy based on the low-rank RBM framework of (Decelle & Furtlehner 2021a) to mitigate gradient degradation in early training, caused by mixing time divergence when encoding distant modes in structured datasets. Our approach extends this framework to real datasets by increasing constrained directions and incorporating trainable bias, essential for image data. We also refine static Monte Carlo sampling for efficient equilibrium sampling. See Section5land SI [I details.
We also provide an in-depth discussion in Sect. 2 on the physical reasons behind the failure of previous algorithms when handling clustered datasets.
2 RELATED WORK
Training EBMs by maximizing log-likelihood has long been challenge (LeCun et al. 2006 Song & Kingma| 202T). The introduction of contrastive divergence (Hinton) 2002) popularized EBMs by using short MCMC runs initialized minibatch examples, but this approach often produces models with poor equilibrium properties, limiting their effectiveness as generative models (Salakhutdinov & Murray 2008; Desjardins et al., 2010 Decelle et al.. 2021). Persistent contrastive divergence (PCD) (Teleman, 2008) improves gradient estimation by maintaining a continuous MCMC chain between updates, acting like slow annealing process. However, PCD struggles with clustered data, as the chain's 's properties can drift from the equilibrium measure, degrading the model (Béreux et al.) 2023). This issue, tied phase coexistence, can be mitigated with constrained MCMC methods, leveraging order parameters from the singular value decomposition of the RBM coupling matrix to reconstruct multimodal distributions (Béreux al. (002).. While effective for model evaluation, this approach is too computationally intensive for practical training, despite producing models with strong equilibrium properties.
The population annealing algorithm, which adjusts the statistical weights of parallel chains after each parameter update, has been proposed as an alternative (Krause et al. 2018). Similarly, reweighting chains using non-equilibrium physics concepts like the Jarzynski equation has been explored (Carbone) et al. 2024b). However, both methods face challenges with highly structured data, where frequent reweighting (duplicating only few samples) skews the equilibrium measure due to insufficient number of independent chains. To mitigate this, more sampling steps or lower learning rates are required, leading to excessively long training times to maintain equilibrium.
3

---

Published as a conference paper at ICLR 2025
In our case, applying these methods to our datasets resulted in training that either too unstable or prohibitively slow when auto-adjusting learning rates and Gibbs steps. Another approach uses EBMs as corrections for flow-based models (Nijkamp et al. 2022), , simplifying sampling but sacrificing the interpretability of the energy function, which is central to our work. Alternatively, evolving flow
models act as fast sampling proposers for HBMs (Grenioux et al. 2023), but this requires training two separate networks and leads to lower acceptance rates as the EBM becomes more specialized.
The Parallel Tempering (PT) algorithm (Hukushima & Nemoto 1996), in which the model is simulated in parallel at different temperatures which are exchanged once in a while using the Metropolis rule, has been proposed both speed up the sampling and to improve the training process (Salakhutdinov [2099] Desjardins et al. (010)) showing important improvements in certain datasets Krause et al. 2020).
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 2021a). 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 the model parameters are changed. The "Stacked Tempering" algorithm (Roussel et al. 2003). recently proposed for RBMs, accelerates sampling by training smaller RBMs with latent variables from previous models, enabling fast updates via a PT-like algorithm. While significantly faster than PT, it unsuitable for learning, as requires sequential training of multiple stacked models.
Separately, a theoretical study (Decelle& & Furtlehner 2021a/ demonstrated that low-rank RBM can be trained via a convex, fast optimization process to reproduce the statisties of the data projected along the d first data principal directions through a convex and very fast optimization process (see (Decelle) & Furtlchner 2021a) and the discussion below) using a mapping Restricted Coulomb Machine. This low-rank model can be seen as good approximation to the correct RBM needed to describe the data, and has the nice property that 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) tuning of the RBM fast sampling process, which was not properly investigated in the original paper.
another model, the so-called
In this paper, show that the training trajectory can be exploited to estimate log-likelihood with minimal computational cost, enabling reliable online computation. We also introduce Parallel Trajectory Tempering (PTT), an efficient sampling method that outperforms optimized MCMC approaches like Parallel Tempering (PT) and Stacked Tempering (Roussel et al. 2023), requiring only a few model parameters stored at different stages of the training. Additionally, we demonstrate that initializing RBM training with low-rank RBM in highly clustered datasets mitigates the effect of initial dynamical arrests, significantly improving model quality at a comparable cost. Combining these strategies allows us to train equilibrium models that capture diverse modes in datasets prone to mode collapse, efficiently evaluate RBMs during training, and sample from multimodal distributions.
THE RESTRICTED BOLTZMANN MACHINE
The RBM consists of N, visible nodes and N₂, hidden nodes. In our study, we primarily use binary variables {0, 1} or 11 for both layers. The visible and hidden layers interact through weight matrix W, with no intra-layer couplings. Each variable further influenced by local biases, 0 for visible units and n for hidden units. The Boobs-Boltzzan distribution for this model is given by
(I) 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, (v))., and the model's marginal distribution over the visible variables, p(v) = En exp [-H(e,, h)]/2 = exp(-H(()) /2 Minimizing the Kullback-Leibler divergence is equivalent to maximizing the
p(u, h) = 1/2 exp (H(
hi] where H(v, h) = -iE ii Wigh - Σ, 0,v - E., nahi.
4

---

Published as a conference paper at ICLR 2025
likelihood of observing the training dataset in the model. Thus, the hog-likelhhod L (HH log Z can be maximized using the classical stochastic gradient ascent. For a training dataset = (m(m) mm=1...,M. the log-likelihood gradient is given by
ac ac (mi) - (M)RBMM Dna = (ha)D (ha)REM,
-wigha) (v.hh))BMM JW..
ac
00,
where (f(u, h))p M Em Eyns f(olm) h)p(hp/(m) denotes the average with respect to the entries in the dataset, and (f(,,))) respect to p(v, h).
Since Z is intractable, gradient estimates typically rely on N, independent MCMC processes, replacing observable averages (o(v, h)) with Σ o(v(r), // where (v(+), h()) are the final states of R parallel chains. Reliable estimates require proper mixing before each parameter update, but enforcing equilibrium at every step is impractical and slow. The widespread use of non- convergent MCMC processes underlies most training difficulties and anomalous dynamics in RBMs, as discussed in (Decelle et al. 2021). In order to minimize cut-of-equilibrinn to keep R permanent (or persistent) chains, which means that the last configurations reached with the
effects, it is often useful
MCMC process used to estimate the gradient at given parameter update t, Pt = {(o)") N ((r)) ) are used to initialize the chains of the subsequent update t + 1. This algorithm is typically referred to as PCD. In this scheme, the process of training can be mimicked as a slow cooling process, only that instead of varying a single parameter, e.g. the temperature, whole set of parameters O (we, O,, 77:) is updated at every step to 0 θ + 7 VL, with VL, being the gradient in equation 2 estimated using the configurations in Pt. and y being the learning rate.
Typical MCMC mixing times in RBMs are initially short but increase as training progresses (Decelle ctt al. 2021). Studies reveal sharp mixing time peaks when the model encodes new modes (Bachtis
al. 2024). phenomenon known critical slowing down Hohenberg & Halperin (197)) associated with critical phase transitions. In RBMs, these transitions are related with the alignment of the singular vectors of W with the dataset's principal components (Decelle al.. 2017 2018 Decelle & Furtlehner] [2021b Bachtis et al.] 2024). However, these dominant modes can be identified via PCA before training and efficiently encoded through convex optimization, circumventing standard maximum likelihood training. This approach prevents gradient degradation due to diverging mixing times the first second-order transitions. The advantages of using the low-rank construction from Decelle & Furtlehner (2021a) as pretraining method will be discussed in Section S
4 TRAJECTORY ANNEALING PARADIGM
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 such as log-likelihood (LL.) estimation through Annealing Importance Sampling (AIS) (Krause et al., 2020), making them inaccurate and meaningless. The evolution of parameters during training can be seen smoother anncaling process compared to the standard temperature annealing used in algorithms like annealing importance sampling (AIS) for LL estimation in sampling algorithms like simulated annealing or PT. This is because it avoids the first-order transitions that typically occur with temperature in clustered datasets, see SI for more detailed explanations of this phenomenon in clustered datasets. However, these methods can still be applied to trajectory annealing using the parameters of different stages of training, rather than scaling the model parameters by a common factor 3,, as done in temperature- protocols.
4.1 TRAJECTORY ANNEALING IMPORTANCE SAMPLING
The standard method for log- LL estimation in the literature is AIS Neal 2001): see (Krause et al.) 2020) for a review on its application to RBMs. A detailed algorithm description is provided in S AIS estimates the partition function of model with parameters O { 0, n) by gradually increasing the inverse temperature 8 from 0 to over N₂ steps. At 3 = O, where Po = Pref- the partition function Zo is tractable, allowing estimation of Z; at neighboring 3; values via path sampling. The reference distribution Prer is typically uniform, but improved performance has been reported when using the independent site distribution centered on the dataset, see SI for details.
5

---

Published as a conference paper at ICLR 2025
23
Figure 2: Comparison of LL estimation error (relative to the exact value) across different methods: AIS (A), AIS with reference distribution fixed to the independent site distribution matching the dataset's empirical center (B), and AI- (ours, C). For Tr-AIS, we evaluate three approaches: online during training, offline using saved models, and with PTT. The RBM, pretrained and trained with 20 hidden nodes (allowing exact LL computation), for 10 gradient steps on the HGD dataset. PIT selects subset of saved models ensuring 0.25 acceptance rate between consecutive models. Lines show the mean LL over 10 independent runs, with shaded areas representing one standard deviation.
AIS is quite costly because it requires large number of annealing steps Ng for reliable estimations, making it impractical to compute the log- likeeiiiihood online during training. As a result, it is often calculated offline using 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 Zo 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 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. 2
while similar qualitative results for the other datasets can be found in the SI, see Figs. [5 16
and 18. In Fig. 2. we display the estimation errors in the LL using three different methods: standard temperature AIS (panel A), AIS with non-flat reference distribution (panel B), and Tr-AIS (panel C). For the temperature AIS estimates, we show results for different values of N₂ 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 just few machines saved during the training 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.
4.2 EXPLOITING THE TRAINING TRAJECTORY FOR SAMPLING
The trajectory annealing paradigm can also be exploited to significantly accelerate the sampling of equilibrium configurations, as we show in Fig. Before presenting our parallel trajectory annealing (PTT) algorithm, we first discuss the limitations of standard methods like Alternating Gibbs Sampling (AGS), and introduce the standard parallel tempering (PT) algorithm. The AGS procedure involves iteratively alternating between two steps: conditioning on all hidden variables given fixed visible variables, and then conditioning on all visible variables given fixed hidden variables.
Problems with the standard AGS method- Let us consider the MNISTO1, HGD, and Ising 2D datasets. MNISTO1 and Ising 2D exhibit strong bimodality in their PCA representations, while HGD is trimodal, as shown in Fig. [ Suppose that we have trained an RBM on each of these datasets, for example using the low-rank RBM pretraining (see Section 5. and we aim to sample the equilibrium measure of these models. Typically, this involves running MCMC processes from random initializations and iterating them until convergence. The mixing time in such cases is determined by the jumping time between clusters. Accurate estimation of the relative weights between modes requires the MCMC processes to be ergodic, necessitating frequent back-and tork--andforrt forth jumps. However, as shown in the chains in Figs. 3 B1-3,, AGS dynamics are exceedingly slow, rarely producing jumps

---

Published as a conference paper at ICLR 2025
A-1
Traj Poralie tompeting
B-1
Alsemntt Gibis Sampling
C-1
dlee eeehanges
10P
0.1-
PC2
int
10
-0.2
-02-
425
2700
0.0
0.2
0.4
0.0
0.2
0.4
10P
10a
10
PC
PC
Nowas (AGS steps/model)
A-2
Traj Paralie tempering
B-2
Gibble Sampling
C-2
0.2
0.2
clos exchingss
10P
10-2
10
0.2
0.0
10P
104
PC
Naaita [ACUS starge (modee)
A-3
Tru Parulle xmpering
B-3
Alsemate Gibis Sampling
C-3
uodhanges
0.2
0.2
102
PC2
0.0
0.0
101
-1.2-
-02-
PC
N (AGS stepe//model)
Figure 3: Comparison of the performance of different MCMC sampling methods on RBMs trained with the pretrainin++PCD procedure on the MNISTO1 (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* MCMC steps using either PTT or 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 averaged number jumps between the two regions separated by the dashed grey line in the PCA plots using different MCMC methods: AGS, Parallel Tempering (PT) (Hukushima & Nemoto 1996). without or with reference configuration (PTref) (Krause al., 2020), the Stacked Tempering (Roussel et 2023), and the Trajectory Parallel Tempering (PTT) proposed in this work. The average is calculated over population of 1000 chains and the shadow around the lines indicates the error of the mean.
even after 10 MCMC steps. The red curves in Figs. C1-3 illustrate the mean number of jumps over 1000 independent chains as a function of MCMC steps, indicating that achieving proper equilibrium generation would require at least 10--107 MCMC steps. To accelerate sampling, we therefore explore optimized Monte Carlo methods inspired by Parallel Tempering.
Trajectory Parallel Tempering- Parallel tempering (PT) used to simulate multiple systems in parallel at different inverse temperatures 3 (Hukushima & Nemoto, 1996; Marinari & Parisi 1992). 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 (=( & cxp(-33ll systems with distincts set of parameters. For set of parameters O we associate the system Pe Pel(x) (æ)). This includes the classic PT case, with 0pr = (())@) where the temperature ladder is typically defined 3())= I/Nr, where 1 = O,..., Nr 1. and Ny represents the number of systems simulated in parallel. Alternatively, recent work (Roussel et al. 2023) introduced a variant of this algorithm called Stacked Tempering, which uses a series of nested RBMs, where the hidden configurations of one machine are exchanged with the visible configurations of the next in the stack, analogous the mechanism in deep belief networks.
at each B. However, a more general approach involves simulating multiple
We introduce a novel variant where model parameters from different training phases are exchanged. Using the learning trajectory, we define the parameter set O' with indexing models at specific training steps. We call this method Parallel Trajectory Tempering (PTT). Configuration exchanges x =0 h) between neighboring models at and L- follow the rule:
Pace (ar' (11) (a--1)=min =min ex, (AH' (2x") AH(((( with AH) =H) -4- (2x), (3)
7

---

Published as a conference paper at ICLR 2025
where again the index indicates that the Hamiltonian is evaluated using the parameters O This move satisfies detailed balance with our target equilibrium distribution p(æ) = (xp(- ensuring that the moves lead to the same equilibrium measure. The reason for using the trajectory is that PT-like algorithms work very well when a system undergoes a continuous phase transition. As 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., 2018: Bachtis et al. 2024). The pseudocodes for the standard PT and our new PTT algorithms are provided in B.
discussed in the SI[G adjusting the inverse temperature
--A13 show the evolution of two independent Markov chains on the most trained model using the PTT sampling, showing a significant increase in the number of jumps compared to Figs. - B--3 We also compare the number of mode jumps observed with PTT with those of other algorithms, as a function of the total number of AGS steps performed. For the PT algorithm. we test different values of Nr and include the variant with a reference configuration (Krause et al. (020). For Stacked Tempering. we replicate the method of (Roussel et al. (023). The total AGS steps are taken as the product of the number of parallel models and the AGS steps per model, as shown in Fig. C1-3.
PTT consistently outperforms all other algorithms in all datasets.
Our PTT approach leverages the ease of sampling configurations from an initial RBM. During classical PCD training, early "non- moonsspciallized models thermalize quickly, behaving like Gaussian distributions. The same is true for RBMs pre-trained from the low-rank RBM of Decelle & Furtlehner (2021a), independent configurations for the first machine can be efficiently sampled by a static Monte Carlo algorithm (see SI A.J. The trajectory flow greatly accelerates equilibrium convergence, particularly in multimodal structures. The time interval between successive models is set to maintain an acceptance probability for exchanges (Eq. ( around 0.25. Notably, the total time required to thermalize PTT chain (AGS steps x number of machines) is unaffected by acceptance values within [0.2. 0.5] (see SB.2 for details). In practice, effective models (i.e. the parameters of the model) 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. PTT offers the advantage of requiring training (unlike Stacked Tempering) or unnecessary temperature simulations and can be fully applied to other mergy-bbsee models. Instead, the equilibrium samples obtained from different models during sampling can be used to compute the likeliie,,oodd, similar to the standard PT approach (Krause) al., 2020). In Fig. --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. 2010), as shown in SI 2 In addition, we evaluated other indicators of sample quality and potential overfitting in ST D Our method consistently provides high- samples across all datasets and, thanks to precise likelihood computation, can effectively detect overfitting.
extra
On the failure of thermal algorithms- Thermal sampling algorithms (such as Simulated Annealing or PT) have difficulties in efficiently sampling models that describe structured data sets. This is mainly because in these types of models, remote modes can simply disappear or appear discontinuously when the temperature is slightly changed. This phenomenon is known in physics as a discontinuous or first-order phase transition. For a detailed introduction to this effect and an illustration of its occurrence when changing the temperature in RBM toy model, see S The presence of this or these temperature transitions can affect the performance of the PT algorithm, as a change between two neighboring temperatures is rarely accepted (and therefore much denser temperature ladders required), but in multimodal models the problem can be even more harmful and difficult to detect: PT can even forget configurations in certain minority clusters while maintaining high swap acceptances. We illustrate this problem in Fig. 4A where we compare the equilibrium distribution p(m) of well-trained machine trained with the HGD dataset when sampled at two different temperatures 3 = 0.9 and 8 = 1. This marginal distribution can be efficiently calculated using the tethered MC algorithm Béreux al. 2023). We can see that from 3 = to 3 = 0.9 the East Asian cluster (the upper right cluster) disappears from the measure. Since the jump to this cluster is strongly suppressed with AGS, as shown in Fig. C2, simulating the system at higher temperatures does not help to find this cluster more easily, but prevents it from being found with standard I-temperature swap + IAGS scheme (because it will be never nucleated at low temperatures). Indeed, in B we show the distribution of samples generated at 3 by a long PT sampling run (with all swap acceptances
8

---

Published as a conference paper at ICLR 2025
Failure
C-
D-pretroin+ped 101
E-LL
Equilibrium distribution
e 0.9
@ ==1
300
0.25
320
240
0.00
160
Mickey
-0.25 0.25
-0.25 0.25
PCU
POU
B- Sampling B-1
PT dataset
FTT vs datasel
M
HGD
Isino2D
PT-AGS
vs AGS
CelebA 32
30366-000
Figure 4: On the left: Panel compares the marginal distribution p(m) along the first two principal components for the pretrained+Por PCD) RBM on the HGD dataset at inverse temperatures 3= and 0.9.9. Panel B presents sampling results at 3 1 using PT, PTT, and 10 AGS steps, initialized from configurations generated by both algorithms. Right panels: Panels C and D compare equilibrium samples from RBMs trained with PCD alone vs. PCD initialized on low-rank RBMs for MNISTO1, Mickey, HGD, Ising2D, and CelebA. Scatter plots and histograms of projections onto the first two principal components show data (black) vs. generated samples (red). Panel E tracks log- ligelikeldoood evolution (train: solid, test: dashed) using online Tr-AIS. Results, averaged over 10 low-rank RBM initializations, show minimal variance, as indicated by narrow shaded regions.
above 0.8) completely forgetting the Asian cluster. One might wonder this is indeed the right balance between clusters learned by the model, so we checked it by performing a long (10% AGS steps) run initialized with these P'T-geneateed samples, finding that the missing clusters gradually resettle. Instead, the PTT samples correctly identify all the clusters and the sampled distribution statistics remain unchanged when a normal AGS run is re-run from that initialisation.
THE LOW-RANK RBM PRETRAINING
In (Decelle& & Furtlehner 2021a). was shown that an RBM with only a few modes in its W can be trained precisely through mapping to Restricted Coulomb Machine and solving convex optimization problem (see A.2 This approach enables training RBMs with simplified W:
W Dd W. u.
with
(u.: U.) E IR x R.M.
(4)
where (wa) are the singular values of the coupling matrix and the right singular vectors (ua correspond exactly to the first d principal directions of the data set. Under this assumption, p(v) is function of only d order parameters, given by the magnetizations along each of the u, components, i.c. (v) Ua v/ NN. In particular,
H ((v) - Ed ma E. log cosh (VN, Ea=1 Dd Uaa M. + na H(m(vv),, (5)
where m = (m...... m.) and f,, is the projection the visible bias onto the extended u matrix. To initialize the visible bias proposed in (Montavon & Müller 2012), add an additional direction

---

Published as a conference paper at ICLR 2025
ua associated with the magnetization my generated by the part of the visible bias that is not contained in the intrinsic space of W. This is an improvement over (Decelle & Furtlehner, 2021a) because effectively adds an additional direction at minimal cost. We also found that the addition of the bias term is crucial to obtain reliable low-rank RBMs for image data. The optimal parameters are obtained by solving convex optimization problem (SI A.2p. This yields a probability distribution p p(m) in lower-dimensional space than p(v), which can be efficiently sampled using inverse transform sampling. However, discretizing the m--sace limits practicality to intrinsic dimensions d S. 4 due to memory constraints. These low-rank RBMs replicate dataset statistics projected onto the first d principal components and, despite their simplicity, can approximate the dataset, as demonstrated in Supplemental Figure 19 for the 5 datasets.
At the start of RBM training, the model encodes the strongest PCA components through multiple critical transitions (Decelle et al. 22177, 2018 |Bachti et al. 2024). Pretraining with a low-rank construction bypasses these transitions, avoiding out-of-equilibrim effects from critical slowing down, which becomes prohibitive in highly clustered datasets. Once the main directions are in- corporated, mixing times stabilize at much lower values, as barriers between clusters shrink when modes align properly [Béreux et al. (2023), making it more reliable starting point for PCD train- ing. However, since mixing times still increase over training, pretraining primarily mitigates poor initialization effects, while issues from insufficient MCMC steps may reappear later. Bypassing these transitions significantly improves model quality and data statistics reproduction, especially for short training times. In Fig. 4 we compare equilibrium samples from two RBMs trained on five datasets with identical settings (see SI(C) but different strategies: one trained from scratch using standard PCD TTeleman 2008), the other using PCD after low-rank RBM pretraining. We also tested Jarzynski reweighting (Carbone et al., 2024a) (see SI ED. but results were unstable on most datasets and are omitted. Across all datasets, pretraining improves sample quality and test log- likelihoodd. The pre-traine++PCD RBM balances modes effectively, as confirmed by PCA projections, and consistently achieves higher log-likelihood values using online Tr-AIS. Samples are drawn with the PTT algorithm, both described in Section 4.1] Scatter plots highlight the importance of pretraining for balanced models in clustered datasets, significantly improving test log-likelihood. However, for less clustered datasets like CelebA and full MNIST, pretraining offers little benefit, splitting distributions into isolated modes is unnecessary (see Fig. in SIŌ).
6 CONCLUSIONS
In this work, we address the challenges of training, evaluating, and sampling high-quality equilibrium Restricted Boltzmann Machines (RBMs) on highly structured datasets. Our primary focus is on improving the sampling process and accurately estimating the log- likeliiod during training. We introduce a novel approach that leverages the progressive feature learning inherent in RBM training to enhance both evaluation and sampling. This method facilitates precise, online computation of log-likelihood and enables an efficient sampling strategy, overcoming the limitations of the standard parallel tempering algorithm. Additionally, our results indicate that for datasets with tightly clustered structures in their principal component analysis (PCA) representation, incorporating a pre-training phase can improve training outcomes. This pre- pre-training involves encoding the principal directions of the data into the model through a convex optimization technique. However, this RBM-specific pre- training method shows limited advantages for less structured datasets or extended training durations, where chains deviate from equilibrium. In the latter case, the presented training approach should be combined in the future with optimized MCMC methods achieve its maximum capacity.
CODE AVAILABILITY
The code and datasets are available at
ACKNOWLEDGMENTS
Authors acknowledge financial support by the Comunidad de Madrid and the Complutense University of Madrid through the Atracción de Talento program (Refs. 2019-T1/7IC-13298 & Refs. 2023- 5A/TIC- -2934), the project 1-220211255000 100021--15500A-1110 y Competitividad, Agencia Estatal de Investigación" (MICIU/AEI/10) 133/99.0033//01111111113111 10001 1033), the Fondo Europeo de Desarrollo Regional (FEDER, UE) and the French ANR grant Scalp (ANR-24, CH23-1320).
financed by the "Ministerio de Economía
10

---

Published as a conference paper at ICLR 2025
REFERENCES
David H Ackley, Geoffrey E. Hinton, and Terrence J Sejnowski. A learning algorithm for Boltzmann machines. Cognitive science, 911):147-1633 1985.
Elisabeth Agoritsas, Giovanni Catania, Aurélien Decelle, and Beatriz Scoane. Explaining the effects of non-convergent sampling in the training of energy-based models. arXiv preprint arXiv:2301.09422 2023.
Daniel J Amit and Victor Martin-Mayor. Field theory. the renormalization group, and critical phenomena: graphs computers. World Scientific Publishing Company, 2005.
Dimitrios Bachtis, Giulio Biroli, Aurélien Decelle, and Beatriz Seoane. Cascade of phase transitions in the training of energy-based models. NeurIPS (2024), arXiv:2405.14689 2024.
R Alvarez Banos, A Cruz, LA Fernandez, JM Gil-Narvion, A Gordillo-Guerrerro Marco Guidetti, A Maiorano, Filippo Mantovani, Marinari, V Mayor,, al. Nature of the spin-glass phase at experimental length scales. Journal of Statistical Mechanics: Theory and Experiment, 2010(06): P06026, 2010.
Nicolas Béreux, Aurélien Decelle, Cyril Furtlehner, and Beatriz Seoane. Learning restricted Boltzmann machine using biased monte carlo sampling. SciPost Physics, 14(3):032, 2023.
Barbara Bravi, Vinod P Balachandran, Benjamin D Greenbaum, Aleksandra M Walczak, Thierry Mora, Rémi Monasson, and Simona Cocco. Probing t-cell response by sequence-based probabilistic modeling. PLOS Computational Biology, 17(9)cc1002297, 2021a.
Barbara Bravi, Jérôme Tubiana, Simona Cocco, Remi Monasson, Thierry Mora, and Aleksandra M Walezak. Rbm-mhe: semi-supervised machine-learning method for sample-specific prediction of antigen presentation by hla-i alleles. Cell systems, 12(2):195-20) 2021b.
Barbara Bravi, Andrea Di Gioacchino, Jorge Fernandez-de Cossio-Diazz Alcksandra M Walczak, Thierry Mora, Simona Cocco, and Rémi Monasson. A cuunfer-learning approach to predict antigen immunogenicity and t-cell receptor specificity. ELife, 12:e85126, 2023.
Nicola Bulso and Yasser Roudi. Restricted Boltzmann machines as models of interacting variables. Neural Computation, 33(10))2646-266) 2021.
Alessandra Carbone, Aurélien Decelle, Lorenzo Rosset, and Beatriz Scoane. Fast and functional structured data generators rooted in out-of-cquilibruum physics. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2024a.
Davide Carbone, Mengjian Hua, Simon Coste, and Eric Vanden- Efficient training of energy- based models using jarzynski equality. Advances in Neural Information Processing Systems, 36, 2024b.
Simona Cocco, Christoph Feinauer, Matteo Figliuzzi, Rémi Monasson, and Martin Weigt. Inverse statistical physics of protein sequences: a key issues review. Reports on Progress in Physics, 81(3): 032601, 2018.
A. Decelle, G. Fissore, and C. Furtlehner. Spectral dynamics of learning in restricted Boltzmann machines. Europhysics Letters, 119(6):60001, nov 2017. doi: 00.1209/095550755/11 19/60011. URL hnttps://dx /ttp/xxx.iioo...oooooo 100p://0....00000335555511119999
Aurélien Decelle and Cyril Furtlehner. Exact training of restricted Boltzmann machines on intrinsi- cally low dimensional data. Physical Review Letters, 127(15):158303, 2021a.
Aurelien Decelle and Cyril Furtlchner. Restricted Boltzmann machine: Recent advances and mean- field theory. Chinese Physics B. 30(4):040202, 20211.
Aurélien Decelle, Giancarlo Fissore, and Cyril Furtlehner. Thermodynamics of restricted Boltzmann machines and related learning dynamics. Journal of Statistical Physics, 172:1576-1609 2018.

---

Published as a conference paper at ICLR 2025
Aurélien Decelle, Cyril Furtlehner, and Beatriz Scoane. Equilibrium and non-equilibrium regimes in the learning of restricted Boltzmann machines. Advances in Neural Information Processing Systems, 14:5345-539,, 2021.
Aurélien Decelle, Beatriz Seoane, and Lorenzo Rosset. Unsupervised hierarchical clustering using the learning dynamics of restricted Boltzmann machines. Phys. Rev. E, 108:014110, Jul 2023. doi: 10. 1103/PhysRevE. 108.0140 URL http::t:: http//iikkke.. aps.nkapppp.. wttpppp..o//ooo/o//e///e///////e http//////.....////////.eeee
108. 014110
Aurélien Decelle, Cyril Furtlehner, Alfonso de Jesús Navas Gómez, and Beatriz Seoane. Inferring effective couplings with restricted Boltzmann machines. SciPost Physics, 16(4):095, 2024.
Guillaume Desjardins, Aaron Courville, Yoshua Bengio, Pascal Vincent, and Olivier Delalleau. Tempered markov chain monte carlo for training of restricted Boltzmann machines. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 145-152.
JMLR Workshop and Conference Proceedings, 2010.
Jorge Fernandez-de Cossio-Diaz. Generative Modeling of RNA Sequence Families with Re- stricted Boltzmann Machines, pp. 163-175. Springer US, New York, NY, 2025. ISBN 978- 1-0716-4079-1. doi: 10.1007/978-100716--000 1-t.----0-4-9991111 URL https:///o...or//1..... 778-1-0716-4079-2
Louis Grenioux, éric Moulines, and Marylou Gabrié. Balanced training of energy-based models with adaptive flow sampling. In ICML 2023 Workshop on Structured Probabilistic Inference and Generative Modeling, 2023.
John Hertz, Yasser Roudi, and Joanna Tyrcha. Ising models for inferring network structure from spike data. 2011.
Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8)11711-1800, 2002.
Geoffrey E Hinton and Terrence J Sejnowski. Optimal perceptual inference. In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, volume 448, pp. 448-453 Citeseer, 1983.
Pierre C Hohenberg and Bertrand I Halperin. Theory of dynamic critical phenomena. Reviews of Modern Physics, 49(3):435, 1977.
Koji Hukushima and Koji Nemoto. Exchange monte carlo method and application to spin glass simulations. Journal of the Physical Society of Japan, 6546);1604-1608, 1996.
Martin Kochmański, Tadeusz Paszkiewicz, and Slawomir Wolski. Curie-weiss magnet-a simple model of phase transition. European Journal of Physics. 34(6):1555, 2013.
Oswin Krause, Asja Fischer, and Christian Igel. Popplation-contratiieedddrgnnn Does consistency help with rbm training? Puttern Recognition Letters, 102:1-7, 2018.
Oswin Krause, Asja Fischer, and Christian Igel. Algorithmss for estimating the partition function of restricted Boltzmann machines. Artificial Intelligennce, 278:03995, 2020. ISSN 00043702. doi: 10.10166j.arintt2209..00 p:1/80004370011330 948
URL
com/ www/retrievee www.:eewwweeeeeeeeeeeeeeeeeeeeee
Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and Fujie Huang. A tutorial on energy-based learning. Predicting structured data, 1(0), 2006.
Renjie Liao, Simon Kornblith, Mengye Ren, David Fleet, and Geoffrey Hinton. censsian-berouull rbms without tears. arXiv preprint arXiv:2210. 2tt::22/00003/// 10318, 2022.
Enzo Marinari and Giorgio Parisi. Simulated tempering: monte carlo scheme. Europhysics letters, 19(6),451, 1992.
Grégoire Montavon and Robust Müller. Deep Boltzmann machines and the centering trick. Neural Networks: Tricks of the Trade: Second Edition, pp. 621-637, 2012.
12

---

Published as a conference paper at ICLR 2025
Radford M Neal. Annealed importance sampling, Statistics and computing, 11:125 1195--11,,,10 2001.
H Chau Nguyen, Riccardo Zecchina, and Johannes Berg. Inverse statistical problems: from the inverse Ising problem to data science. Advances in Physics, 66(3):197-261, 2017.
Nijkamp, Gao, Sountsov, S Vasudevan, Pang, S-C Zhu, and YN Wu. Meme should mix: learning energy energy-based model with neural transport latent space meme. In International Conference on Learning Representations (ICLR 2022).. 2022.
Erik Nijkamp, Mitch Hill, Song-Chun Zhu, and Ying Nian Wu. Learning non-convergent non- persistent short-run run meme toward energy- based model. Advances in Neural Information Processing Systems, 32, 2019.
Erik Nijkamp, Mitch Hill, Tian Han, Song-Chur Zhu, and Ying Nian Wu. On the anatomy of meme-based maximum likelihood learning of energy-based models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp- 5272-5280, 2020.
Sebastian Quiroz Monnens, Casper Peters, Kasper Smeets, Luuk Willem Hesselink, and Bernhard Englitz. The recurrent temporal restricted Boltzmann machine captures neural assembly dynamics in whole-brain activity. June 2024. doi: 10.7554/elifc.98555 98449... URL http: /txp www. org/
10. 70..7554ee fe. 98489.1
Clément Roussel, Jorge Fernandez-de Cossio-Diaz, Simona Cocco, and Remi Monasson. Accelerated sampling with stacked restricted Boltzmann machines. In The Twelfth International Conference on Learning Representations, 2023.
Ruslan Salakhutdinov and lain Murray, On the quantitative analysis of deep belief networks. In Proceedings of the 25th international conference on Machine learning, pp. 872-879, 2008
Russ R Salakhutdinov. Learning in markov random fields using tempered transitions. Advances in neural information processing systems, 22, 2009.
Ya G Sinai. Theory of phase transitions: rigorous results, volume 108. Elsevier, 2014.
Paul Smolensky. Parallel Distributed Processing: Volume by D. Rumelhart and J. McLelland, chapter 6: Information Processing in Dynamical Systems: Foundations of Harmony Theory, 194-281. MIT Press, 1986.
Alan Sokal. Monte carlo methods statistical mechanies: foundations and new algorithms. In Functional integration, pp. 131-192. Springer, 1997.
Yang Song and Diederik P Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.03288. 2021.
Tijmen Tieleman. Training restricted Boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th international conference on Machine learning, pp. 1064-1071, 2008.
Jérôme Tubiana, Simona Cocco, and Rémi Monasson. Learning protein constitutive motifs from sequence data. Elife, 8:e39397, 2019.
Thijs L van der Plas, Jéröme Tubiana, Guillaume Le Goc, Geoffrey Migault, Michael Kunst, Herwig Baier, Volker Bormuth, Bernhard Englitz, and Georges Debrégeas. Neural assemblies uncovered by generative modeling explain whole-brain activity statistics and reflect structural connectivity. Elife, 12:e83139, 2023.
Martin Weigel, Lev Barash, Lev Shchur, and Wolfhard Janke. Understanding population annealing monte carlo simulations. Phys. Rev. E. 103:053301, May 2021. doi: 10. 103/P///ssss.EE.... 103. 053301. http:/ttt://ttti aps wtwp:/poog/o/////p//poppo Ttt///.0..11///pp/// 1111111//////// http:////..nnan...rr.....//g.
103. 053301
Jianwen Xie. Yang Lu, Song-Chun Zhu, and Yingnian Wu. A theory of generative convnet. In Maria Florina Balcan and Kilian Q. Weinberger (eds.), Proceedings of The 33rd Inter- national Conference Machine Learning, volume 48 of Proceedings of Machine Learn- ing Research, pp. 2635-2644, New York, New York, USA, 20-22 Jun 2016. PMLR. htttps://proceediinge www.sm.....p...mmmmmmmmmmmm
URL
eccieeelci html
13

---

Published as a conference paper at ICLR 2025
Andrew Yale, Saloni Dash, Ritik Dutta, Isabelle Guyon, Adrien Pavao, and Kristin P Bennett. Generation and evaluation of privacy preserving synthetic health data. Neurocomputing, 416: 244-255, 2020.
Burak Yelmen, Aurélien Decelle, Linda Ongaro, Davide Marnetto, Corentin Tallec, Francesco Montinaro, Cyril Furtlehner, Luca Pagani, and Flora Jay. Creating artíficial human genomes using generative neural networks. PLoS genetics, 17(2)ee100333,, 2021.
Burak Yelmen, Aurélien Decelle, Leila Lea Boulos, Antoine Szatkownik, Cyril Furtlehner, Guillaume Charpiat, and Flora Jay. Deep convolutional and conditional neural networks for large-scale genomic data generation. PLOS Computational Biology, 19(10)=e101 1584, 2023.
14

---

Published as a conference paper at ICLR 2025
Supplemental Information
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- 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 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 dataset, we want to find a good set of model parameters (W, 0 and v7) 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 these a = 1..... d projections ma Ua v VN, the magnetizations along the mode a, where un is the a-th mode of the PCA decomposition of the dataset. simple way to encode these d-modes is to parameterize the w-matrix
W
Wo U. Uu
with
(ua, Un) IR x R
(6)
where u and are respectively the right-hand and left-hand singular vectors of W, the former being directly given by the PCA, while marginal energy on the visible variables, H(v) = log Σn exp H(v, h) can rewritten terms of these magnetizations
the singular values of W. Using this decomposition, the
(m1,..., mad
H(v) =
Σ log cosh
v VN.u. ug
Σ mama + = H(m(vv). a=1
Now, the goal of our pre-training is not to match the entire statistics of the dataset, but only the marginal probability of these magnetizations. In other words, we want to model the marginal distribution
Pemp (m)
Σ Pemp (v)
II ma a=1
VN, u ).
(8)
where 5 is the Dirac õ-distristribution. In this formulation, the distribution of the model over the magnetization m can be easily characterized
p(m)
Σe H(())
II
8 ma
1
(9)
VN.
N (m) exp
log cosh
Ue
ma + 7J
(10)
a=1
((m))N
1
(mmm
(11)
Z
where ((m) = E. IT 8 (m. -NN u
is the number of configurations with magnetiza-
tions m, and thus S(m) = log N ((m)/N is the associated entropy. Now, for large Ny the entropic term can be determined using large deviation theory, and particular the Gärtner-Ellis theorem:
S(m))
Pprior (m)
exp (-N,I(m))).
(12)
2N.
with the rate function
I(m) = sup [m 1 {(1)]] m. 11" - 21(1))
(13)
15

---

Published as a conference paper at ICLR 2025
and
x(1))
lim
(
log (emmmu))
lim No0
log
e N/ Eozi No Ex Moatr
(14)
No0 Nv
Nv
2N,
lim N.
log cosh
N a=l
Ma Ma.i
(15)
Then, given magnetization m, compute the minimizer μ (m) of Φ(1+) m 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 assuming N, large enough. As result we get μe (m) satisfying implicit equations given by the constraints given at given N.:
N uf tanh
ma
N
HB
(16)
B=1
It is then straightforward check that spins distributed
Pprior (v/m) x m)(v)
mm(v)
(17)
fulfill well the requirement, as (u v/ VNv) Pprior
= ma. In other words, we can generate samples
having mean magnetization ma just by choosing 7a
Pprior (v: =1 mmmm) sigmoid 2/N,, N, Σ Wa (m)
(18)
The training can therefore be done directly in the subspace of dimension d. In Ref. (Decelle & Furtlehner 2021a), 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 mapping from the RCM to the RBM to recover the RBM's parameters. brief, the training of the low-dimensional RBM is performed by the RCM, and then the parameters obtrained via 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 202la), , it is possible to precisely train surrogate model for the RBM, called the Restricted Coulomb Machine (RCM), on 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:
H(v) Σv:0, u2O, Σ
log cosh
[Wwne Win na
(19)
We then project both the parameters and variables of the RBM on the first d principal components of the dataset:
N. WinWios Wias
N 0
U.Wia
Woa
0. =
0guio
(20)
==1
with a E {1,.... then given by
d} and u the projection matrix of the PCA. The projected distribution of the model
exp N S(() + Ed= Pama N, EN No log cosh (VN, Ed Wa + Ma ))
PRBM(N))
Z
(21)
where we ignore the fluctuations related to the transverse directions and S[m) accounts for the non-uniform prior due the projection of the uniform prior for the way compute it.
16

---

Published as a conference paper at ICLR 2025
The RCM is then built by approximating
log cosh(2)) = |a|| log
(22)
which is valid for large enough. The probability of the RCM is thus given by:
exp Nr (()) + Σazi Ma + Eos 9a Ea Re max ma
Za
(m)
(23)
Z
where
Waa
1/a
4a
N 2
007
n
Za
(24)
0=1
Ea=
w
NN, Nr
Ea=
w
This can be easily inverted as
aa
VN. no
and
f/
in order to obtain the RBM from the RCM. The model is then trained through log- fog-likelihood maximization over its parameters. However, this objective is non-convex all the parameters are trained through gradient ascent. To relax the problem, since we're in low dimension, we can define family of hyperplanes (n, x) covering the space and let the model only learn the weights of each to the hyperplane. We can then discard the ones with low weights to keep the approximation (equation 22)
good enough.
The gradients given by
(m-
n2a m + zal] - E [m-ppc((m) [n m + 2u
(25)
Oga
Emppp (m) [ma] Em--pa((m) [ma]-
(26)
a0,
The positive term is straightforward to compute. For the negative term, we rely on discretization of the longitudinal space to estimate the probability density the model and compute the averages.
A.3 LEARNING THE BIAS
To initialize the model with its bias the dataset's mean, we define the 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:
N.
H(m(v)) = Σ Oama
Σ log cosh N.
1
No
ma Weo + 1a
(27)
a=1
This implies that this direction independent of the others within the model. Consequently, we can create independent mesh for this direction to learn it. The rest of the training is performed using the same procedure as before.
B SAMPLING VIA PARALLEL TEMPERING USING THE LEARNING TRAJECTORY
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 equi- librium exhibit faster and more ergodic dynamics compared to poorly trained models, the sampling time can still be excessively long when navigating 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, method involves simultaneously
17

---

Published as a conference paper at ICLR 2025
RBMo
RBM
RBM
RBM,
(0)
(0)
(0)
Initialization
mcme step
Swap
(1)
(3)
(1)
RCM sampling
(only pre-train+PCD)
Sampling
1 mame step
(2)
(2)
(only JarRBM and PCD)
40
11
:
Figure 5: Scheme of PTT. We Initialize the chains the models by starting from a configuration mo and passing it through the machines along the training trajectory, each time performing & meme steps. For pre-train+PCID, ( is sampling from the RCM, otherwise it is uniform random initialization. The sampling consists of alternating one meme step for each model with swap attempt between adjacent machines. For pre-train+PCD, at each step sample new independent configuration for RBMp using the RCM.
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-orde 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.
sketch of the Parallel Trajectory Tempering (PTT) is represented in fig. Specifically, we save I models at checkpoints / 1, 1,...,4 44 along the training trajectory. We denote the Hamiltonian of the model at checkpoint t as H,, and refer to the Hamiltonian of the RCM model as Ho- We define GübbsSampling(" H. æ, k) as the operation of performing & Gibbs sampling updates using the model H starting from the state æ. In all our sampling simulations we used &
The first step is initialize the models' configurations efficiently. This involves sampling N chains from the RCM model, zo RCMSampling(?4") and then passing the chains through all the models from GibbsSampling(H:: (0)) (-11.)). k).
performing & Gibbs steps at each stage: xt (0)
1 to t
The sampling process proceeds in steps where we update the configuration of each model except Ho with h Gibbs steps, and sample a completely new configuration for the RCM model Ho- Following this update step, propose swapping chains between adjacent models with an acceptance probability given by:
(ace( <> 1t- min (1, (aH,( (x:
(28)
where (x) = (x) H,-1 We continue alternating between the update step and the swap step until a total of Nmeme steps reached.
A comparison of performances between PTT and standard Gibbs sampling is reported in Fig 3
B. PSEUDO-CODE OF PTT VS PT
We provide here the pseudo code for both Parallel Tempering (Alg. I and Parallel Trajectory Tempering (Alg. 2 sampling schemes. The difference between both algorithms lies in the way the parallel models are selected. In the case of Parallel Tempering, an annealing on the temperature is made, whereas for Parallel Trajectory Tempering several models are saved during training.
18

---

Published as a conference paper at ICLR 2025
Algorithm Parallel Tempering
Require: {3,]]-(1.... Ng). {w. m). Nizer Ninerement Nsample for j 1,.... N3 do chains]j] RandomB wawdomminarychaaaaaaaaaaaaaaaaaa wawyhhainsaaaaagaaaaaaa
end for
for i = I,... Niter do N3 do
for j chains t htt.::/eeee.iiiilaaalaaaaaaaaaar
end for
for jorj=1. N3 1 do chains chai [j+] Shacconfiguraatonnnchhhmms....,1 chainsjjj+
1], 8 )
end for
end for
Algorithm 2 Parallel Trajectory Tempering
Require: (0, w 0, 17 ( Nicern Ninerement Naample
for j = l..... Nmodel do
chainsij] RandomiSnarryCaaanaa....
end for
for i = 1, Niter
for j = 1, Meedel do chains[j] + #atereeneeeediiisaaaaaannaaaaaaa chains(j) end for for 1,.,...,,,Nae - do chains(ji, chainsjj +1] swapconfigggaattoomncchhaaiii., chaishh
1], + 1)
end for end for
19

---

Published as a conference paper at ICLR 2025
dataset
Texp
Tint
MNISTOI
258(29)
9.0(15)
HGD
48(10)
3.8(7)
Mickey
2(0)
1(0)
Ising
18(3)
2.4(5)
Table Autocorrelation times associated to the random walk dynamies between models in the PTT sampling obtained with the models trained with Dre-training+PCCD on different datasets. Models used the PTT are all chosen with the 0.25 criterion.
B.2 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 al. 2010). From now on, will refer only to the number of models used the PTT sampling, denoted as Nm. In an ergodic sampling process, each run should spend roughly the same amount of time each of the N.. different model indices, all models have the same probability.
The mean temperature/model index is known be (n) (Nm - 1)/2 (assuming the model indices run from 0 to Nm 1). This simplifies the computation of the tine-autocorrelation function of the model indices nt(t) visited by each sample at time
(n, + to) (n)) (n))
C(t)
(n:(1n) (²))²
Here, the symbol (...) refers to the average over all samples and initial times to obtained from our PTT runs. This allows us to compute the exponential autocorrelation time, Texp- from a fit to
C(I) Ae-t/raa t/Texu for large 1
as well as the integrated autocorrelation time, Tint, using the self-consistent equation (Sokal 1997, Amit & Martin- Mayor 2005):
Grint
C(I).
Tint
t=
These two times provide different insights: Texp gives insight into the time needed to thermalize the simulations (setting the length of the runs to at least 20Texp is generally a safe (and very conserva- tive) choice to ensure thermalization (Sokal 1997)). while Time gives the time needed to generate independent samples. To keep track of the fluctuations, we compute a different C" (t) for each m I,, I...., Nn system index simulated parallel in the PTT, where each CTL (t) 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) CT' and use this curve to extract Texp and Tiat- The errors C' and the times are then extracted from the error of the mean the values extracted with each of the Cm (t) curves.
For example, in Fig. (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 Texp = 258 258(29) (29) and Tine = 9.0(1.5) for MNISTO1, and Texp = 48()) and Tint 3.8(7) for HGD. These numbers guarantee us that the 10* steps used to generate the scatterplots in Fig. were long enough to ensure equilibration in this dataset, but also in the others, as shown in Table
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, 0.25. Let justify this choice. The total physical time we need to simulate Nm models grows proportionally with Nm. (unless one manages to properly parallelize the sampling, which is feasible but did not do it), while the time needed to flow between very well trained and poorly trained models controlling the mixing time scales as /N aa, where a is the average acceptance of the exchanges between the models. Since N. grows with in practice (as shown in the inset of
20

---

Published as a conference paper at ICLR 2025
MNISTO1
MNIST01
2000
Na.
Lotal Unoe Nandia
fit Texp
10-1
(())
1500
10
C(t)
0.5
10
1000
Target 8OC
8000
10-3
6000
10
500
1000
0
MCMC steps
0.2
0.4
0.6
target ace
HGD
HGD
100
200
20
N
3000
Yop Nurabi
fit Texp
10
150
C(t)
0.0
0.5
-2000
10 -2
100
Garget
10-3
50-
-100
total
200
400
MCMC steps
0.0
0.2
D-4
0.6
target
Figure 6: (Left) We show the autocorrelation function C(()) obtained in Eq. equation B.2 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 PTT run of 5 10a MCMC steps, where the models used were automatically selected during training using the PTT moves acceptance criterion of 0.25. (Right) In blue show the Texp function of the target acceptance rate for MNISTO1 (top) and HGD (bottomm. The error bar results from the error of the mean of the values extracted from the different (t) and the shading marks the range between the largest and the smallest value of the Ni curves. In orange, we show the total simulation costs (Texp multiplied by the number of models to be simulated in parallel) as function of the target acceptance rate. This number of models saved, function of the acceptance rate, is shown in green in the inset.
Fig- rigght), it is expected that the total time can be optimized for given value of In Fig. right we show Texp in blue (and Texp.Nm. in orange) as a function of a for the data sets MNISTOT (top) and HGD (bottom). These figures show that the choice of is particularly crucial long it is between 0.2 and 0.5.
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
610
784
60%
40%
MNIST
50 000
784
60%
40%
21

---

Published as a conference paper at ICLR 2025
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-tran Dee-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
PCD
2000
2000
10 000
0.01
100
200
Mickey
PCD
2000
2000
10 000
0.01
100
200
Fre-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
All experiments were run RTX 4090 with AMD Ryzen 9 5950X.
D COMPUTATION OF THE LIKELIHOOD USING AIS
To compute the likelihood of the RBMs used the Annealed Importance Sampling method intro- duced by (Neal 2001) and recently studied in the case of the RBM (Krause et al. 2020).
D.1 OFFLINE AIS
We follow the method described in (Krause al.) 2020). The AIS estimate is based on creating sequence of machines, at one end a machine from which we can compute the partition function and the other end the machine from which we want to compute the partition function. Let's consider that we have i = 0,.... Nrbm machines described by
e (#)
(.)
Pr(x)
where Po Pref denote the RBM from which we can compute the partition function (typically RBM where the weight matrix is zero) and pNa.m. = Ptarget machine from which we want to estimate the partition function. For each machine, we denote its set of parameters O {w, Og, m.)- Following (Krause et al. 2020), defining a transition rule [1 + x;) in order to pass from configuration T;-1 1 of the system to æ; of the system that have the following identity
Z Nota
Nebm
exp
Σ
Hi+1(w)) H.(
(29)
Zo
ForwardPath
There, if one compute Zo. the estimate of Z is obtained by averaging the r.h.s. of Eq. 29 over the forward trajectory starting from a configuration Dn using the operators T. While the method exact, the variance of the estimation relies close two successive distribution Pt, Pi+1-
22

---

Published as a conference paper at ICLR 2025
Table 4: Training time of the RBMs
Name
Training time (s)
CelebA
PCD
119.2 + 0.8 3134.0 0.6
Pre-train+PCD
HGD
PCD
106.2 + 0.3
Pre-train+PCD
266 + 3
Ising
PCD
62 +
Pre-train+PCD
104 + 2
Mickey
PCD
118.3 + 246.0 + 0.7
Pre-train+PCD
MNIST-01
PCD
100.6 + 0.3 258.0 + 0.9
Pre-train+PCD
MNIST
PCD
105.7 + 0.2 233.9 + 0.2
Pre-train+PCD
Table 5; RCM Hyperparameters
Dataset name
Dimensions
Hidden nodes
Training time (s)
CelebA
100
3006
HGD
3
100
150
Ising
3
100
32
Mickey
2
100
118
MNIST-01
100
148
MNIST
3
100
119
Temperature AIS:
a commonly used implementation of AIS, is to multiply the Hamiltonian of
the system temperature, and to the different systems indexed by by of 3 where 30 33 < 33+1 and 3N.... = 1. Therefore,
equivalent to choosing θ = /3 ww 0, n) Another version of Temperature AIS consist in choosing slightly different reference distribution Po that take into account the bias of the model. In this setting, the intermediate distribution are taken as
(1- Bb) Pref e 3.))(2)
(x(x)
Z,
The main difference is that a prior distribution is used reference po Pref and the temperatures interpolate between this prior and the target distribution. It is quite common to take
exp ( 0.0; + Σ. haha) II, cosh (.)) II. cosh(na)
Prof (z)
the general derivation of AIS allow in principle to choose the set forward proba- bility and/or intermediate system sees fit. Since again, change of temperature in the system
Trajectory AIS:
23

---

Published as a conference paper at ICLR 2025
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. 3) is about 0.3). Then the formula Eq. 29 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 very precise estimation.
E TRAINING OF THE RBM USING THE JARZYNSKI EQUATION
In this section, we describe the procedure introduced in (Carbone et al. 2024b) 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.) 2024b), 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. 2021).
In one of its formulations, the Jarzynski equation states that observable O with the average obtained through many repetitions of an out-of-equilibrium dynamical process. Let consider a path y {(x", 0²). (x=, ('),..... (x² 0¹), where x² (v', h²) refers to the model variables and 0 = (w', 11',, 0²), the model parameters, both at time i.
relate the ensemble average of an
If consider the training trajectory of RBM, Po P1 average of an observable O over the last model of the trajectory as
Pt-1 - Ps, write the
(0), = (OE
DR Oly,)) ER WF
(30)
Wt
traj
where 8 is many different trajectory realizations. In particular, W, is given by the recursive formula:
importance factor and the averages on the rhs are taken across
(') = H²'x') - H2()))) - Σ log T' - 1) + Σ log T' - y²) (31)
1=1
w-- + (x0)) H' (x log T' 1 -
(32)
T (1 yy)
Where T and T are, respectively, the forward and backward transition probabilities.
The transition matrix between two states y² 1 + y' is factorized into two different subsequent movements:
1. Update the variables -²-1 + x² on the fixed model 0--1
2. Update the parameters 0--1 - 0' on fixed variables x'.
In other words,
T(y - = wo = 1x (-1 - x Mm (0'-1 +0).
(33)
and
T(y² + 11= = tza (0' + (01 (x' + x 1111))
(34)
so that
T (yy - -11)
Whr-1 1 xx' + x
ta (0' + 0
log (yy - ) Wgg (1 + x
llog
+ log
(35)
La (00- 02)
24

---

Published as a conference paper at ICLR 2025
We know the ratio between the transition probabilities between x'~1 - æ', because evolution of the variables during the simulation is performed using MCMC simulations that satisfy detailed balance. This means that:
W (x' + x²' x²'1) wgr- (x'-1
log
H H1- (x'- 21'11(1'"))
(36)
1- + x²)
and we are left with
t. (0' - QC-1)
(y) (x) + - log
H1-
(37)
t (0
0t)
The standard formula that we use in the paper assumes that the changes in 9 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 simulating 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. 2024b)
statistically irrelevant, and expect to get large fluctuations in the
R D -
ESS
2W2 E [0,
(38)
R- R- 1
which measures the relative dispersion of the weights distribution. Then, way of circumventing the weight concentration on 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 (W(r) 0 Vr l...... R).
F TUNING OF THE HYPERPARAMETERS
In order to justify our choice of hyperparameters in the main text, illustrate on Fig. the value of the log-liiel 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 constant number of gradient updates of 10000. In this regime, the machines that perform the best have a learning rate of ? ~ 0.1, 0.005 and about No, 500, 1000 hidden nodes. We can expect that smaller learning rates reach the same value of the log- likeehooooo but for much longer training. We clearly see the benefit of the RCM, where the ligeeiklidood achieves much higher values in general.
G CONTINUOUS PHASE TRANSITION
In statistical physics, and more precisely in mean-field model, model described by the Gibbs- Boltzmann distribution is studied in the infinite size limit, N - oc, and characterize the state of the system by investigating the free energy of the system, defined by minus the log partition function
F(0)
log( (
One of the goals is thus determine the typical configurations that should be expected given a 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 observe for the parameters 0. For this purpose, let us rewrite the partition function as a sum over all possible values of the energy
= De- Σ
E+S(E) Σ
e -F(E)
Σ -N((E)
where (EE) the number of states with 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) FEE)N is the free energy per variable. All Z, E and F depend on the model parameters. N refers the number of variables.
25

---

Published as a conference paper at ICLR 2025
PCD
RCM+PCD
-00
-00
0.1
L.44
10L11
138
2006
67313
28325
0512.5
0.1
435
A111
118
62.47
59193
24127
322L10
382
2056
0.05
1.35
100.0
1.38
30.21
111
0.05
0.33
-0.T7
EIT
23
121
2081
Leuraing Take
0.01
136.00
177
06.12
15.11
7106
72.00
M900 m
0.01
100
1.0
41.75
70
-6.27
1990 099
120
120
0.005
410-00
1112
0.005
-11.4
-0.44
-8.00
-140
-140
160
160
0,0055
0.333
1:3
0.10
6:1
0,0006
OUT
1.0
-180
-180
0.0001
1.:2
3.10
101
111
393
0.0001
1012
0.04
105
121
1 1.0
121
-000
-200
20
860
600
180
Nomber of bidklen modes
Namber of hichien modees
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.
For N + oo, 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(F) 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 mapped to a very simple model, the so-called Curie- weies model, as recently described in Ref. (Bachtis et al. (024). The Curie-Weiss model is fully solvable, and for this reason the meaning of continuous transition is perhaps casier to understand there. Let's give it try. The Curie-Weit model consists in a set of N discrete variables 8: = 11. We define the (exponential) distribution over these variables as
p(5)
Z
exp
8,8j 3H Σ si)
where 3 and H are parameters of the model and 7 the normalization constant also called partition function. The question to understand the "structure" of the distribution as a function of 3
Case H 0: expect that for 3 small each spin behaves variables (each spin being +1 with probability one-half), while for large value of 3, 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 + co. way to study this distribution is to investigate the moment generating function. It is possible show (rigorously in this precise model), that the probability of m = N ", St, also called the magnetization of the system, (in the large N limit) where 32m) is large deviation function.
isolated Rademacher random
given by p(m) x In the CW model, it is given by
2(m) = 3m3
2 log 2 cosh(Bmm).
The structure of 0mm) such that for < 3c = 1, is a convex function with only one minimum in m O, while for 3 > Be it has two symmetric minima m =m(3) F O with ((m(3) = ((-m(3)) (the two minimum are equally probable). The values (in function of (8) of the m are continuous: m = 0 for 3 E [0, 1]. and then positive until it saturates to 1 as 3 - co. In practice, it means that the distribution p(m) pass from an unimodal distribution to a bimodal one as 3 increases. The point where it happens is at B. 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 B., the function 2(m) develops non-analyticity in its second-order derivative,
26

---

Published as a conference paper at ICLR 2025
0.0
0.75
-0.50
-062
-0.025
25
0:2
-064
000
-2.50
-1.
-05
0.5
10
1.0
-2.5
0.66
dle El 2-0.3
cim di 11 2-1
-0.68
c3m d- El 111
O1m de 2) 2w12
-0.725
-070
12ft 5 1.151 Ie 3.02
2.551 2.0.
6ett A 2:115
2.151-5el
0.72
-2.775
D20 2151
1.00
-1.00
-0.75
-0.25
0.00
0.25
0.50
1.00
-0.75
-0.50
-0.25
0.25
0.50
0.75
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
Figure 8: Left:
500. For N + oo they converge toward å distribution. Right: The case H + O. In this example, we illustrate what happen where the bias parameter / is changed a temperature
1.15 1. We see that for small values of H, S2 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.
that be linked with long-range fluctuations. For recent review (Kochmański et al. 2013). and for more rigorous results (Sinai 2014).
Case + in this case we expect that the system will always be biased toward configurations where the spins have the same direction(sign) (sign) as the bias (or field) H. Yet, for small value of 3, the bias will be quite small, hence the magnetization will be small but not zero, m ~ 3H. When is large, the system will now have a stronger and strong magnetization toward the magnetic field m (H). Now we can again compute the large deviation function of the system
3m3
((m) 2 log [2 cosh (Bm)) + Hm.
Again, it is possible to analyze the function 2 function of 3 and H. For 3 the function is convex and there is only one minimum located at m > 0. When 3 > two situations can be distinguished. First, let's denote
B
3 (8-1)- 3
(B) =
8 atanh
If H > (3), then the large deviation function has only one minimum for positive m. Now if H < <H)p(3), then Q will present two minima. first one for positive m and a second one for negative m. The reason is that the pairwise interactions "manage" to create metastable state for negative magnetization even in the presence of a positive field or bias. It clear that this state will be unstable w.r.t. the with positive magnetization which traduces in minimum of higher value in the large deviation function. By changing the value of // from large positive values toward large negative ones, we then pass from 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 Ω is at its lowest value. However, if one prepares system close to a metastable state, the system will remain trapped there until it disappears.
How is this related to the training of the RBM ? First, it has been shown in previous works (De- celle et al.) 2018 Bachtis et al.) 2024) that the RBM's training undergoes 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
continuous. Then in (Decelle & Furtlehner 2021a) is shown that the relaxed version of the RBM, namely the RCM undergoes a first-order transition when changing the temperature. Schematically understand that, it is enough to see how
27

---

Published as a conference paper at ICLR 2025
the parameters 3 acts on the system. In the case of the RCM, the energy term behaves in the standard way, i.. linearly with the inverse temperature 3 and the free energy takes the ordinary form
log(Z)) BF 3E S.
Therefore 3 act as a trade-off parameters between the energy E and the entropy S. By varying we change the contribution of the energy while keeping the entropy fixed. Hence, when the corresponding large deviation function of the RCM has reached point having many minima with the same value of the free energy up to ((1/N) corrections, all these minima realize in principle different trade-off between energy and entropy. As a result a slight change in will offset this sensible balance and the respective free energies will acquire differences of order ((N) leading to favor dramatically one state against the others, i.e. the one having originally the lowest [resp. the largest] energy when is increased. It is very likely that by changing 3 discontinuous transition is encountered, where the minima of Ω vanishes far away from each other rather than a continuous one. Concerning the RBM the argument for the occurence of such 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. 3. 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-unifon 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 = (()) ) compared to a large embedding space (d = O(N), N > 1) which is quite common in our point of view, at least for the type of data we interested in.
To conclude, moving the parameters along the learning trajectory, leads the system to undergo continuous transition where (O( (s) is very close to po(t+ (s). Changing the temperature tends to produces discontinuous transition where Pa(ce (8) can happen to be very far from Pa(c+&t)e (s).
H FIRST ORDER TRANSITIONS IN RBMS
In this section, examine the appearance of first-order transitions when adjusting the temperature of a well-trained RBM on clustered datasets- a phenomenon that significantly hinders the performance of thermal sampling algorithms. We first illustrate this issue using an analytically tractable toy model, followed by an analysis of the first-order temperature transition observed in the RBM trained on the HGD dataset, as discussed in the main text.
H.1 THEORETICAL ANALYSIS IN A SIMPLE MODEL
We propose a simple model of dataset which will show that the learned RBM is having first order transition in temperature. We first consider an artificial dataset using a Curie- model. We defined the following Hamiltonian on the variables s, #1
H
Σ
p(s) = Z exp(-Br
2N.
When the inverse temperature 3r is below one, the distribution in the No + oc limit is unimodal, with the Boltzmann average of the variables s being (s) O. For 3r > the system becomes bimodal, with (s,) = m, where m satisfies the self-consistent equation m = tanh(3rm)) We propose to model this dataset using Bernoulli-Gauss RBM with one hidden node. The visible nodes, o (0. 1}. follow Bernoulli distribution, while the hidden node, т, follows Gaussian distribution with zero mean and variance 1/N,. combined with a local visible bias. When the inverse temperature 3T is below one, the distribution in the N. + oc limit is unimodal and the average w.r.t. the Boltzmann distribution of the values of the variables s's is (s) = O. When 3y is above one, the system becomes bimodal with {Is,)) m when is defined by the self-consistent equation m = (anh(3rm), We now propose to learn this dataset using a Bernoulli- Gaulls RBM with one hidden node, where the visible nodes are (0, 1}. and the hidden one is Gaussian distribution
28

---

Published as a conference paper at ICLR 2025
0.20
=105
0.15
0.10
0.00
0.05
-0.10
0.15
-0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Figure 9: Free energy -m for e [0.8, 1.05| for the optimal RBM's parameters learned on CW model with 8T 1.4. At At 1.0 the two minima are perfectly adjusted. As change the temperature of the RBM, they get destabilized showing the presence of first order transition.
of zero mean and variance 1 1/N, and local visible bias. The Hamiltonian is then
o.W. H = -o.W..
Σ noo
where W; is the weight matrix and ?/ the biases. The leaning then consist converting the binary variable of the Curie-Weiss model into 10 and to learn the coresponding dataset. By rapid inspection (in the infinite size limit), we can easily show that the optimal learned parameters for the RBM
wi = 2/ BT and 17 = 2,37.
Now, thanks the simplicity of the model, we can actually compute the free energy of the model where the Hamiltonian has been rescaled by an annealing temperature B, we obtain
m 2
-f(mT)
N-1
Σ log [1 + exp(3w,m, + Bm:]]
where m,
N Σ
sugsim [1 + exp(Srcmm + [m1]
In Fig. -3. -35 (m))) for different values of 3 within the range [0.8, 1.05). is evident that represents the coexistence point, indicating a first-order transition in temperature. Note that the probability p(m) Re e (N3(((mr) for large N, implying that for 3 > 1, the typical configuration corresponds to 1, whereas for B < corresponds simple case we observe such transition when doing annealing.
1. This clearly shows that even
H.2 FIRST ORDER PHASE TRANSITION IN TEMPERATURE ON A TRAINED RBM ON A REAL DATASET
We now consider the RBM trained on the HGD dataset using the pretraining and 50 000 updates of PCD training. As we vary the temperature,
recover the probability p(m) learned by the model using the TMC method (Béreux et al., 2023), projected along the second principal component in Fig. and on both first and second principal components Fig. 10 While < 1.0, the probability distribution is dominated by the central mode, equilibrating between all modes B = 1.0. As goes
above 1.0, the distribution becomes dominated by the two external clusters.
The effect is clearly visible when sampling the model using the PT algorithm (Fig. 12 The sampling was performed with 200 temperatures with a minimum acceptance rate of 0.78 between two temperatures. We see dataset cluster being ignored with only a few samples in it. However, when resampling sampling the model for 10th AGS steps starting from the samples obtained through PT, we clear change in the distribution, equilibrating with the missing cluster. In contrast, performing the same experiment with PTT instead of PT, we observe no shift in the empirical distribution the samples after 106 AGS, as the upper cluster already sampled.
29

---

Published as a conference paper at ICLR 2025
g 0.6
B = 0.8
8 = 0.9
==100
β ==1J 1.1
40
0.25
32
24
(100) 900
0.00
400
0.25
320
240
(0))))
0.00
160
80
0.25 0.00 0.25
-0.25 0.00
0.25
0..25
0.00 0.25
0.255..00
0.25
0.25).000.25
PCD
PCO
PCO
PCO
PCD
Projection of the negative log probability (top) and probability (bottom) learned by the model, constrained to visible magnetizations mo and m., for various values of 3. The red dots represent the training data projected onto the first two principal components.
Figure 10:
log p(m)
p(m)
80
0.60
15-
0.90
80
1.00
1.10
10
40
20
0-
0.22
0.0
0.2
0.1
0.2
0.0
0.2
0.1
PC1
PC1
Figure 11: Negative log probability (left) and probability (right) learned by the model constrained to have given my for several values of B. The equilibrium distribution of the model is displayed as the black dotted line (3 = 1).
1 OVERFITTING AND PRIVACY LOSS AS QUALITY INDICATORS
In this section, 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 the models trained with PCD with and without pre-training We focus on the human genome dataset, as shown in Fig. [.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) al. 2021 2023) have thoroughly investigated this for 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 (AArs) and privacy loss indicators, introduced in Ref. (Yale et al.) 20020). which quantify the similarity and the level of "privacy" of the data generated by model w.r.t. the training set. We have AATS = 1 (AATrue + AAsymth) where AATru [resp. AASymth are two quantitics [o, 1] obtained by merging two sets of real and synthetic data of equal size N, and measuring respectively the frequency that real [rep. synthetic] has synthetic [resp. real] as nearest neighbor. the generated samples are statistically indistinguishable from real samples, both frequencies AATrue and AAAyyuth should converge 0.5 large N,. AArs be evaluated both with train test samples and the privacy loss indicator is defined Privacy loss AAtest -ttest-A.. AAArrin and is expected to be strictly positive. Fig. [4 shows the comparison of AATs and privacy loss values obtained with our models, demonstrating that the pre-trained RBM clearly outperforms the
30

---

Published as a conference paper at ICLR 2025
A
PC1
0.1
0.1
0.3
-0.1
0.2
0.0
PTT
PC+AG5
PTT-AGS
0.2-
02-
PC1
0.1
0.0
-0.1
-0.3-
-02
-02
0.0
PCD
PCD
Figure 12: Comparison of the sampling results of PT (left) and PTT (right). The top row corresponds to sample generated with each methods compared with the dataset distribution. The bottom row compares samples obtained by performing 10 AGS starting from the PT and PTT samples with their starting position.
HGD
MICKEY
0.55
MNISTIB
AA Tunth
0.90
AX Thuth
AA Syath
0.2
Prinucy
0.000
Princy
0.085
0.05
Precazy Ioar
-77
Tria
-50
-170
-000
Toat
100
Tick
Thilning updetos)
Training Line (gradenn updites)
Troming (gradbent (pplated)
Figure 13: 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 I.I. coincides with a sharp rise in the privacy loss.
other model, and even achieves better results (AATS values much closer to 0.5) than those discussed in (Yelmen al. 202I).
EXTRA FIGURES
31

---

Published as a conference paper at ICLR 2025
A
1.0
HGD
B
1.0
HGD
AA Truth
0.5
AA Synth
0.5
AA Truth
AA Syuth
0.0
0.0
0.2
0.2
Privacy loss
Privacy loss
0.0
(nats)
-00
Train
LL (nats)
375
Train
400
-450
Test
Test
425
100
10
102
10"
10
10a
10a
10
Training time (gradient updates)
Training time (gradient updates)
C
D
Privacy loes
0.7
AASya unThohh
0.05
0.00
Prinacy Lowe
0.5
0.83
0.4
0.3
0.02
0.2
0.0
0.00
pretem
PPCD
PCD
Model
Figure 14: Comparison of the evolution of overfitting metrics along the log-likelihood between RBM trained from scratch using PCD (A) and pretrained RBM (B). (C) shows the comparison between the AASyn and AATruth between the two machines. The machines selected in order to have a comparable privacy loss (D), right before starting to overfit.
32

---

Published as a conference paper at ICLR 2025
AIS #
AIS with refernce diteribution
Thajetory extraaiion
TeAI:: mndias
Einee w.x. raset 14k
10F
PTT: 20 madimee
coline estination
10F
000
107
11
107
317
10
10
Trainnng time jendiizz spdato))
Training linn (grabbest spilatt)
Trinsigg time (padiea spian))
Figure 15: Comparison of log-likelihood estimation w.r.t. the exact Il on the MNIST-01 dataset.
AIS /
AIS with refenace dastcilution
Thajectory exiraation
muchine
107
PTT: 24 audiine
Shar w.zk coest
unline retination
100
10
5/0
5000
1f
10
10
10
10
10P
10
Tioiting Kin: (predicst sg-iates)
Trokiing time (godient updates)
Thokinn time ip (gadicat updates)
Figure 16: Comparison of log-likellnn log-likelihood estimation w.r.t. the exact Il on the MNIST dataset.
AIS ,
AI8 with referenee dstribbtion
Thijetoey extiraaions
10-13
Blnce widk cart
30
10-3
TeAAs mudhines
10-
10
500
onllne extmmatinn
107
10
10"
Tinising Line Spadicat updates)
Traiiing Lime Spoeicca spdates)
Thiniig time Spadical updates)
Figure 17: Comparison of log-likelihood estimation w.r.t the exact Il on the Ising dataset.
AIS /
A18 with refermne distiiaattooooo
Thijetoey extirsations
Enooe wat. cluet L1
10at
50
machito mochitue
ooline osiiationn
11
10P
31F
107
10F
Thaining Uime (grddeest spdistee)
Training time fgrdleest updistee)
Thiining time (g Gradiet spiiste)
Figure 18: Comparison of log-likelihoo estimation w.r.t. the exact on the MICKEY dataset.
33

---

Published as a conference paper at ICLR 2025
Lonkrank RBM generation
A- MNIST01
B Mickey
- Human Genome
DCM
0 1 0
45
001
0 1
05
62
0 0
Ear
0 0
1 1
0
02
PC3
Ising 0-0.4
E CelebA
B/W 32x32
Figure 19: Samples generated with the low-rank RBM trained with 3 directions plus 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, also show the images generated.
MNIST
100
RCM-PCD
-12
LL (tant)
-100
-00
-200
10"
10
10a
10a
Training time faradient undates!
Figure 20: We show the evolution of the tr-AIS log- likelihoodd 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 particular advantage but neither is disadvantage.
34


Reviews:
[{"review_id": "N9ew0ZbjE8", "reviewer": "Reviewer_jkKM", "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.", "questions": "*** introduction\n\nWhilst 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.\n\nNotation: inconsistent use of $N_v$ and $N_{\\text{v}}$ throughout, similarly for $N_h$.\n\nEquation 1: it might be better to write W_{i\\alpha}, rather than w_{i\\alpha} since w is used later for the \"singular values\".\n\n*** page 2\n\nFigure 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.\n\n*** page 5\n\nPlease clarify the difference between \"model averages\" and \"observable averages\" and the difference between using N_s independent MCMC processes and R parallel chains.\n\nPlease clarify for the reader the meaning of <v_ih_a>_D\n\nSection 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.\n\nPlease 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? \n\nIn 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.\n\nFor 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.\n\nIt'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.\n\n\n*** Section 5\n\nI 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.\n\n*** Section 6\n\nIn 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.\n\n\n*** Supplementary material\n\nThe 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?\n\nThe notation \\hat{u} is confused with \\bar{u}.\n\nWhy use $w$ here whereas $W$ is used in the main text?\n\nThe 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.\n\nPlease state clearly what are the parameters of the model that are being learned.\n\nSection A.2. The claim as before of exact training is incorrectly made here.\n\nThe 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?\n\nIn 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).\n\nAn 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.", "limitations": "", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 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.", "weaknesses": "The paper suffers from a lack of clarity of presentation and lack of clarity of novelty.\n\nThe 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.\n\nIn 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "dTkt8ZpMfq", "reviewer": "Reviewer_axdy", "summary": "The claimed novelties of this work are twofold. \nFirst, 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. \n\nAs 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. \n\nNumerical 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).", "questions": "## Questions, Small comments and typos\n- 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. \n- Is there any intuition behind the bump observed in Figure 3 at around $10^3$ gradient updates (left and middle plot).\n- 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. \n- 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. \n- Line 235: I'd recommend adding a reference for critical slowing down. This comment applies to earlier occurrences of this concept.\n- Line 459: grew -> grey\n- Line 512: Banos et al. (2010) might need to be wrapped in parenthesis \\cite -> \\citep", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written and easy to follow. \n- It represents a pleasant read that is accessible to a broad audience. \n- The literature review and related work section read well and are exhaustive.\n- The idea of pre-training the RBM to encode the principal components is simple yet very effective. \n- 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). \n- The numerical experiments look solid and aligned with the theoretical insights given in the main text. \n- I have not thoroughly checked the mathematical details in the appendix, but at first glance, they look good.", "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.  \n- 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. \n- 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:\n    - How beneficial is to have such an augmented direction for the bias compared to the naive approach proposed in Decelle and Furtlehner (2021a)?\n    - 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?\n\nThis 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).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "4awmr2FgZN", "reviewer": "Reviewer_dgLM", "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.", "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?\n- What is the definition of $\\bar{u}$?\n- Could the authors provide a detailed derivation of Equation (4)? The terms $\\bar{u}_{a}$ and $\\eta_{a}$ are currently undefined.\n- 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?\n- Is it possible to learn DBM without pre-training using the pre-training with weights introduced by [1] ?\n\n[1] Yuma Ichikawa and Koji Hukushima, Statistical-mechanical Study of Deep Boltzmann Machine Given Weight Parameters after Training by Singular Value Decomposition.", "limitations": "", "rating": 5, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 3, "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.\n- To my knowledge, extending replica Monte Carlo methods to a learning trajectory is original and intriguing.\n- Including a specialized physics background in the Appendix makes the paper accessible even to readers without a physics background.", "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. \n\n- 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.\n\n- 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.\n\n- 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.\n- 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.\n- 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.\n- The CelebA data in Figure 2 appears to be truncated.\n\nBecause the high performance has been verified numerically, the score can be raised if the above statement is cleared.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yUQ1s4hU3E", "reviewer": "Reviewer_mx3V", "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.", "questions": "1. Why didn't you apply this to larger problems?\n\n2. What are situations where the pre-training fails?\n\n3. Is PTT useful for generating samples during training, using only earlier parts of the training run?", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The idea of low-rank pre-training is interesting and seems like it could be useful if it scaled up.\n\n2. The idea of doing AIS across the training run is creative and clever.\n\n3. Parallel tempering across training steps seems new.", "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.\n\n2. 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\nsimilar.  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?\n\n3. 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.\n\n4. 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).\n\n5. 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.\n\n6. The title is too broad relative to what the paper delivers.\n\nTypos:\n - L161-162: \"two slow\"\n - L478: \"exchanges parameters\" but I think you mean \"exchanges configuration\".\n - L775-776: \\bar{u} vs \\hat{u}.\n - L836-837: \"gradient is convex\" -- surely you mean the training objective is convex in the parameters.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)