Nicolas Béreux
email: nicolas.bereux@gmail.com

Aurélien Decelle

Cyril Furtlehner

Beatriz Seoane

Learning a Restricted Boltzmann Machine using biased Monte Carlo sampling

Restricted Boltzmann Machines are simple and powerful generative models that can encode any complex dataset. Despite all their advantages, in practice the trainings are often unstable and it is difficult to assess their quality because the dynamics are affected by extremely slow time dependencies. This situation becomes critical when dealing with low-dimensional clustered datasets, where the time required to sample ergodically the trained models becomes computationally prohibitive. In this work, we show that this divergence of Monte Carlo mixing times is related to a phenomenon of phase coexistence, similar to that which occurs in physics near a first-order phase transition. We show that sampling the equilibrium distribution using the Markov chain Monte Carlo method can be dramatically accelerated when using biased sampling techniques, in particular the Tethered Monte Carlo (TMC) method. This sampling technique efficiently solves the problem of evaluating the quality of a given trained model and generating new samples in a reasonable amount of time. Moreover, we show that this sampling technique can also be used to improve the computation of the log-likelihood gradient during training, leading to dramatic improvements in training RBMs with artificial clustered datasets. On real lowdimensional datasets, this new training method fits RBM models with significantly faster relaxation dynamics than those obtained with standard PCD recipes. We also show that TMC sampling can be used to recover the free-energy profile of the RBM. This proves to be extremely useful to compute the probability distribution of a given model and to improve the generation of new decorrelated samples in slow PCD-trained models. The main limitations of this method are, first, the restriction to effective low-dimensional datasets and, second, the fact that the Tethered MC method breaks the possibility of performing parallel alternative Monte Carlo updates, which limits the size of the systems we can consider in practice.

Introduction

In Physics, Restricted Boltzmann Machines (RBMs) are only a spin system or an Ising model in which the coupling matrix between spins has a bipartite topology. In Machine Learning, it is a famous generative model introduced many decades ago by Smolensky [START_REF] Smolensky | Information processing in dynamical systems: Foundations of harmony theory[END_REF] as one of the first generative models capable of successfully learning complex/real datasets. It was later popularized by Hinton [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF][START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF], who (at the time) proposed a viable algorithm for training these machines. After a brief period of notoriety in which RBMs were used as a pretraining method for larger neural networks [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF], they fell into disuse after 2010 with the introduction of new generative models such as Generative Adversarial Network (GAN) [START_REF] Goodfellow | Generative adversarial nets[END_REF] and Variational Autoencoder (VAE) [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. Indeed, these models have a clear advantage: despite eventual instabilities of the training process, they do not rely on complex and costly Monte Carlo simulations. Nevertheless, RBMs have made a strong comeback in the field of Computational Biology in recent years, especially for purposes of pattern extraction [START_REF] Tubiana | Emergence of compositional representations in restricted boltzmann machines[END_REF][START_REF] Shimagaki | Selection of sequence motifs and generative hopfield-potts models for protein families[END_REF][START_REF] Tubiana | Learning protein constitutive motifs from sequence data[END_REF][START_REF] Tubiana | Learning compositional representations of interacting systems with restricted boltzmann machines: Comparative study of lattice proteins[END_REF][START_REF] Van Der Plas | Compositional restricted boltzmann machines unveil the brain-wide organization of neural assemblies[END_REF][START_REF] Bravi | Rbm-mhc: a semisupervised machine-learning method for sample-specific prediction of antigen presentation by hla-i alleles[END_REF][START_REF] Yelmen | Creating artificial human genomes using generative neural networks[END_REF].

RBMs remain an interesting generative model for many reasons. On the ML side, RBMs are simpler models than GAN or VAE. They have only one hidden layer, resulting in simple and potentially interpretable learned features. They can also learn complex datasets and are practical even when the number of samples of that dataset is limited. During the last decade, they have also regained the interest of physicists [START_REF] Tubiana | Emergence of compositional representations in restricted boltzmann machines[END_REF][START_REF] Huang | Statistical mechanics of unsupervised feature learning in a restricted boltzmann machine with binary synapses[END_REF][START_REF] Decelle | Thermodynamics of restricted boltzmann machines and related learning dynamics[END_REF][START_REF] Barra | Phase diagram of restricted boltzmann machines and generalized hopfield networks with arbitrary priors[END_REF][START_REF] Decelle | Restricted boltzmann machine: Recent advances and meanfield theory[END_REF], given their closeness to the wellknown Ising model, for which computing the phase diagram in its many variations has been a long-standing hobby in the Statistical Physics community. In the case of RBMs, it is more difficult to obtain a relevant phase diagram, since the coupling constants of the trained model must be correlated due to the learning process. Nevertheless, the phase diagram obtained in simpler cases can at least be used as a guide to understand the macroscopic behavior of RBMs in the high--temperature oe linear regime of the model. Lately, some recent work also has attempted to analyze the learning dynamics of the models [START_REF] Harsh | place-cell'emergence and learning of invariant data with restricted boltzmann machines: breaking and dynamical restoration of continuous symmetries in the weight space[END_REF].

Despite all these efforts, training RBMs remains a difficult task, especially because of the instability of the training procedure, which is mainly related to the use of nonconvergent Markov Chain Monte Carlo (MCMC) processes to estimate the gradient during training, as stated in a recent work [START_REF] Decelle | Equilibrium and non-equilibrium regimes in the learning of restricted boltzmann machines[END_REF]. In this work, we present a new learning method that uses biased sampling Monte Carlo methods [START_REF] Landau | A guide to Monte Carlo simulations in statistical physics[END_REF]. We will show that the classical alternative sampling MCMC method combined with Glauber Heat-Bath dynamics used in training RBMs is impractical for some datasets: the mixing time of the chain is simply too long. By using biased MC sampling, we can avoid this situation, especially at the beginning of learning. The drawbacks of this method is first that it breaks the conditional independence of the two layers of the RBM. Because of that, parallelization is less efficient and therefore restrain the size of the dataset one can be considered. Second, it needs to discretize an order parameter into possibly different directions, reducing its use when only few directions are needed.

This paper is organized as follows. We start with a definition of RBMs and their classical training procedures in Section 2. We then motivate our work by discussing the break of ergodicity phenomena encountered in RBMs trained with data following multimodal distributions. Then, in Section 3, we introduce the Tethered MCMC (TMC) sampling method and in Section 4, and its combination with the stochastic gradient ascent during the learning process in Section 4. [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF]. In Section 6, we compare the new TMC training strategy with the standard training scheme in several artificial and real datasets.

Definitions

An RBM is an Ising spin model with pairwise interactions defined on a bipartite graph of two noninteracting layers of variables: the visible nodes, v = {v i , i = 1, ..., N v }, represent the data, while the hidden nodes, h = {h a , a = 1, ..., N h }, are the latent representations of this data and attempt to establish arbitrary dependencies between the visible units. Typically, the nodes are binary-valued in {0, 1}, yet Gaussian or arbitrary distributions on real-valued bounded support are also used, ultimately making RBMs adaptable to more heterogeneous datasets. Here we are concerned only with binary {0, 1} variables for both the visible and hidden nodes. Other approaches that use truncated Gaussian hidden units [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] and provide an activation function of type ReLu for the hidden layer work well, but our experiments should not be affected by this choice. The energy function of an RBM is taken as

E[v , h; w , b, c] = - ia v i w ia h a - i b i v i - a c a h a , (1)
with w being the weight matrix, and b, c the visible, and hidden biases (or the magnetic fields in the Physics language), respectively. The Boltzmann distribution of such a model is then given by

p[v , h|w , b, c] = exp(-E[v, h; w , b, c]) Z , ( 2
)
where Z is the partition function. RBMs are usually trained using gradient ascent of the loglikelihood function, , of the training dataset, = {v (1) , • • • , v (M ) }, being defined as

(w , b, c| ) = M -1 M m=1 ln p(v = v (m) |w , b, c) = M -1 M m=1 ln {h} e -E[v (m) ,h;w ,b,c] -ln Z.
The gradient of is then composed of two terms: the first one accounts for the interactions between the RBM's response and the training set, and the same for the second, but using the samples drawn by the machine itself. The expression of the gradient w.r.t. all the parameters is given by

∂ ∂ w ia = 〈v i h a 〉 -〈v i h a 〉 , (3)
∂ ∂ b i = 〈v i 〉 -〈v i 〉 , (4)
∂ ∂ c a = 〈h a 〉 -〈h a 〉 , (5)
where

〈 f (v , h)〉 = M -1 m {h} f (v (m) , h)p(h|v (m) ),
denotes an average of an arbitrary function f (v, h) over the dataset, and 〈 f (v , h)〉 , the average over the Boltzmann measure in Eq. [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF].

A large part of the literature on RBMs focuses on schemes for approximating the r.h.s. of Eqs. 3-5, usually referred to as negative term of the gradient. Theoretically, these terms can be computed to arbitrary accuracy using parallel MCMC simulations, as long as the number of MC steps is large enough to ensure a proper sampling of the equilibrium configuration space. The most popular approximation schemes for computing the negative part are: (i) the Contrastive Divergence (CD) introduced decades ago by Hinton [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF], or its later refinement the (ii) Persistent CD (PCD) [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]. In both cases, the negative part of the gradient is estimated after a few k sampling steps, typically k ∼ (1), and the differences are related only to different clever ways of initializing the Markov chains. In the CD recipe, the parallel MCMC simulations are initialized from the same samples of the mini-batch that are used to compute the positive part of the gradient, while in the PCD recipe, the final configurations of each Markov chain are stored from one parameter update to the next to initialize the subsequent chain. In practice, the consequences of the lack of convergence of all these MCMC processes under these approximations during learning have very rarely been studied in the Literature, and it has recently been shown to have dramatic and nonmonotonic effects on the generation performance of RBMs [START_REF] Decelle | Equilibrium and non-equilibrium regimes in the learning of restricted boltzmann machines[END_REF]. In fact, this work has shown that the CD method is generally a very poor training method because it fits models whose equilibrium distribution is drastically different from the dataset distribution. In contrast, the PCD recipe appeared to fit RBMs with good equilibrium properties in image datasets, but the quality of the samples generated by these machines was limited by the number of sampling steps k used during learning.

A closer examination of the PCD scheme shows us examples where this strategy also fails completely in learning good model parameters. For example, it was shown in [START_REF] Decelle | Exact training of restricted boltzmann machines on intrinsically low dimensional data[END_REF] that simple artificial datasets with low effective dimension, such as well-separated point clusters, cannot be learned correctly with PCD. One might naively think that such a simple dataset should be straightforward to learn with an RBM, considering that this problem is easily solved with a Gaussian Mixture model. On the contrary, the main problem of such "clustered" datasets is that it is very difficult to ergodically sample the configuration space of multimodal distributions. In a one-dimensional case, i.e., clusters separated along a spatial direction, the typical probability density profile along this direction should correspond to well-fitted Gaussian peaks centered at each cluster center and separated by regions with essentially zero probability. It is well known that in such a case the MCMC chains can take an extremely long time to jump from one cluster to another. In Physics, we would speak of large free energy barriers to surmount, and it is exactly the same phenomenon encountered in the vicinity of a first order phase transition. This shows that a learning procedure that relies on parallel MCMC samplings to estimate the correlations of the model should be extremely costly or even fail miserably when trying to learn clusterized data.

Many interesting datasets form compact clusters after dimensional reduction, which is the typical case of genome or protein sequence data, for example. The previous discussion anticipates that it must not only be difficult to obtain reliable and effective models for such datasets using standard learning schemes (such as CD /PCD), but also to evaluate the quality of a given trained model. Indeed, standard sampling techniques are unable to equilibrate, which means that a convergent generation of samples may be prohibitively long. As a consequence, approximate techniques typically used to compute , such as Importance Sampling [START_REF] Krause | Algorithms for estimating the partition function of restricted boltzmann machines[END_REF], must be completely inaccurate on such datasets.

The numerical problems associated with the apparition of multiple coexisting phases have long been known in Physics or Chemistry. One of the possible solutions is to break the metastability with external constraints. This idea can be easily implemented using biased MC or reweighting methods, such as the Umbrella Sampling method [START_REF] Torrie | Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling[END_REF], microcanonical methods [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF][START_REF] Martin-Mayor | Microcanonical approach to the simulation of first-order phase transitions[END_REF], or the Tethered MC (TMC) method [START_REF] Fernandez | Tethered monte carlo: Computing the effective potential without critical slowing down[END_REF][START_REF] Fernández | Equilibrium fluid-solid coexistence of hard spheres[END_REF][START_REF] Martin-Mayor | Tethered monte carlo: Managing rugged freeenergy landscapes with a helmholtz-potential formalism[END_REF], which is adapted here to training and generation sampling of RBMs.

Breakdown of ergodicity

Multimodal distributions are notoriously hard to sample because simulations remain trapped in one or several metastable free energy minima for a very long time and fail to visit the entire configuration space. We illustrate this difficulty in Fig. 1. For this purpose, we considered an artificial dataset with 5 point clusters in 1000 dimensions, where all 5 are aligned in a single dimension. For this dataset, it is possible to train an RBM exactly as described in Ref. [START_REF] Decelle | Exact training of restricted boltzmann machines on intrinsically low dimensional data[END_REF]. We use this perfectly trained RBM to generate new data using either standard MCMC moves (using Glauber-Heath-Bath dynamics), -a sampling procedure that will henceforth be referred to as SMC -, or the TMC method that we will introduce later. In both cases, the term MC sweeps refers to an attempt to update all N v and N h variables of the model.

Before proceeding with the discussion of Fig. 1, let us briefly discuss the projections we use to compare the original and generated datasets. We begin by extracting the principal components of the dataset ω α with α = 0, • • • N v -1 according to standard Principal Component Analysis (PCA). Thus, ω 0 marks the direction of maximum variance of the dataset, ω 1 the next, and so on (with the constraint that all {ω α } are orthonomal). Now we use these directions to project each data point, v (m) ,

m (m) α = v (m) • ω α N v , (6)
which defines the datum's magnetization along the α principal direction. In Fig. 1, we show the first two magnetizations, m 0 ) for the different sets of data. Note that one needs just one direction to split up the 5 clusters (they are aligned by construction), and it is straightforward to capture such a direction with the PCA. The choice of Figure 1: We illustrate the difficulties of sampling a multimodal distribution with standard Glauber MC moves of a dataset of 5 aligned point clusters in high dimensions. On the top, we show the projection of the dataset/generated samples along the first (m 0 ) and the second (m 1 ) dataset's principal components. Each dot is one sample. In the bottom, we show the histogram of all these m 0 . The original dataset is displayed in black, and the generated datasets are shown in blue or orange, depending on the algorithm used to sample the configuration space of the exactly trained RBM from Ref. [START_REF] Decelle | Exact training of restricted boltzmann machines on intrinsically low dimensional data[END_REF], from random initializations. Orange samples are obtained after 10 5 standard Glauber MC (SMC) steps and blue samples, after 30 TMC sweeps. using the PCA in general will be justified in Sect. 5. Clearly, from a random initialization, local SMC moves get trapped in just one of the 5 data clusters even after performing 10 5 MC sweeps, thus failing to visit (and generate) data in the rest of the clusters. The breakdown of ergodicity is easy to understand here: in order to jump from one cluster to another with local moves, the trajectory must visit configurations with an intermediary m 0 , and such configurations have extremely low probability. In the TMC sampling, one can enforce the simulations to visit these rare but crucial intermediate states by introducing external constraints. As a result, this sampling method is able to generate samples in all the five clusters with the correct associated weights.

Sampling with the tethered MC method

We use the TMC sampling method [START_REF] Fernandez | Tethered monte carlo: Computing the effective potential without critical slowing down[END_REF][START_REF] Martin-Mayor | Tethered monte carlo: Managing rugged freeenergy landscapes with a helmholtz-potential formalism[END_REF] to sample efficiently the configuration space of clustered models. TMC is a refinement over the popular Umbrella Sampling method [START_REF] Torrie | Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling[END_REF] that permits simplifying notably the reconstruction of canonical expected values. This simplification will be crucial to combine this sampling technique with the stochastic training dynamics. In the TMC method, a soft constraint m(v ) ∼ m1 is introduced in the partition function Z via the identity

α/2π ∞ -∞ d m e -α
where α controls the strength of the constraint. In this way, Z can be written in terms of the TMC effective potential, Ω( m),

Z = ∞ -∞ d m e -αΩ( m) , (7)
with

e -αΩ( m) = α 2π {v ,h} e -E(h,v )-α 2 [m(v )-m] 2 . ( 8
)
Then, the probability of m is given by

p( m) = e -αΩ( m) Z . (9)
Clearly, when α → ∞, p( m) recovers the canonical probability of finding

m(v ) = m, that is, 〈δ (m(v ) -m)〉.
Finite values of α soften this constraint, which is necessary to facilitate sampling. We can use this probability to define a new ensemble, the tethered ensemble, where m(v )∼ m. Then, the expected value of an observable O(h, v ) in this ensemble is

〈O〉 m = {h,v } O(h, v ) ω(h, v , m) {h,v } ω(h, v , m) , ( 10
) with ω(h, v , m) = e -E(h,v)-α 2 [m(v )-m] 2 . ( 11
)
Now, the principal innovation of the TMC method is realizing that

Ω ( m) = dΩ d m = 〈 m -m(v )〉 m , (12)
which means that Ω( m) can be recovered up to any desired precision, from a numerical integration of Monte Carlo averages of several independent simulations at fixed m

Ω( m) -Ω( mmin ) = m mmin dΩ d m d m = m mmin m -m(v ) m d m . ( 13
) So does p( m), e -αΩ( m) mmax mmin e -αΩ( m ) d m → p( m), (14)
as long as, mmin and mmax , the two extrema of the integration interval, are safely chosen so that e -αΩ( m) is essentially zero beyond that range. Also, p( m) ≈ 〈δ (m(v ) -m)〉 for large α.

In this work, we use typically α ∼ 10 4 , which means that the approximation is rather good.

In addition, we can exploit the TMC strategy to estimate the free energy barrier between two metastable states. Interestingly, it is straightforward to show that the canonical average of O, 〈O〉 , that is, the average with respect to the Boltzmann measure of Eq. ( 2) (the quantity one wants to compute in practice), can be recovered using the following identity:

〈O〉 = ∞ -∞ 〈O〉 m p( m). ( 15
)
In other words, model averages can be recovered up to any desired precision just by integrating TMC averages (with no approximation involved). Such a turnover between ensembles is only useful if sampling ergodically the TMC measure [START_REF] Tubiana | Learning compositional representations of interacting systems with restricted boltzmann machines: Comparative study of lattice proteins[END_REF] is much easier than in the original unconstrained problem (2). This is particularly true when only one phase is stable at a given m. In the absence of metastabilities, thermalization is extremely fast. This means that one can use the TMC weight of Eq. ( 11) to avoid phase coexistence (and long time dependencies) in your MC simulation.

Example: TMC Sampling of the 5 aligned cluster dataset

Let us come back to the 5 clusters dataset example used to illustrate the ergodicity problem (discussed around Fig. 1). In the light of the preceding discussion, it is now clear that we can efficiently sample the whole configuration space running TMC simulations with m running in the range of observed m 0 . We have summarized the TMC sampling procedure for this example in Fig. 2. The first step is to discretize the range of m 0 needed to integrate numerically Ω ( m) in [START_REF] Bravi | Rbm-mhc: a semisupervised machine-learning method for sample-specific prediction of antigen presentation by hla-i alleles[END_REF] and to obtain p( m) [START_REF] Yelmen | Creating artificial human genomes using generative neural networks[END_REF]. To do this, we need to select a good interval of m, that is, a mmin and mmax . We can do this using the minimum and maximum projection of the dataset along ω 0 , plus a safety border to ensure zero probability beyond these points. In order to avoid scaling problems between RBMs models when defining m, we have always normalized the data projections conveniently, so that they lie in the interval [0, 1], and considered extra borders of ±0.1 or ±0.2 to define the range of m. Then, we discretize uniformly the interval [ mmin , mmax ] using N m points (for the sampling in Fig. 2 we considered 250 discretization points, even though we show much less in the sketch to enhance the comprehension). Then, we run independent TMC simulations at each of these m values. In general, we consider several parallel runs for each m to estimate Ω ( m) in ( 12). This is much faster than iterating TMC many times because they can run in parallel. Note that, the TMC weight of Eq. ( 11) breaks the independence of the visible units provided by the bipartite lattice structure, which means that we lose the possibility of updating all visible variables at the same time (the socalled alternating Gibbs sampling) which makes the TMC simulations much slower than the unconstrained model. Once we have computed Ω ( m) for all the discretized values of m, we integrate this function numerically to extract Ω and p( m). We compare the p( m) extracted following this pipeline2 with the histogram of m 0 obtained from the dataset in Fig. 3, showing a perfect match between both distributions. Finally, one can use this p( m) to improve the estimation of the gradients via Eq. ( 15) during the learning process, as we will discuss below, or just to generate samples from our model. In this second case, the procedure is simple, draw a set of msampling values from the p( m) distribution, and run independent and short TMC runs (initialized at random or other non-informative configuration) at these ms. The distribution of the projections of the samples generated using this procedure were shown in blue in Fig. 1, showing again a perfect match with the distribution of the original dataset.

Generalization to multiple order parameters

The method can easily be generalized to k constraints. Let m = (m 1 , • • • , m k ) represent these parameters. The method remains the same: a constraint is imposed on the parameters at the same time, hence the weights of the TMC ensemble are now given by :

ω(h, v, m) = e -E(h,v)-α 2 k l=1 ( ml -m l ) 2 . ( 16
)
The gradient of the potential is now

∇Ω( m) = 〈 m -m〉 m , (17)
and the potential can be recovered at fixed m through the numerical integration of ∇Ω( m):

Ω( m) = Ω( mmin ) + ∇Ω( m )d m , ( 18
)
where represents any path leading from mmin = ( m1 min , • • • , mk min ) to m. Finally, p( m) can be recovered analogously to the single order parameter case [START_REF] Yelmen | Creating artificial human genomes using generative neural networks[END_REF]. In black, we show the Tethered probability p( m) obtained from the numerical integration of Ω( m) plotted against 〈m〉 m. We compare this probability with the histogram of the projection of the dataset along ω 0 . An almost perfect match is seen, even though the two probabilities should only match in the limit α → ∞.

Training RBM with the Tethered Monte Carlo method

The TMC method is therefore very efficient to sample quickly isolated clusters. It is then straightforward to use it to compute the negative term of the gradient eqs. [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF][START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. For this goal, we need to run several independent simulations at fixed m and estimate 〈s i 〉 m, 〈h a 〉 m and 〈s i h a 〉 m. We then recover the unconstrained negative term of the gradient using a numerical integration as in [START_REF] Huang | Statistical mechanics of unsupervised feature learning in a restricted boltzmann machine with binary synapses[END_REF]. In order to improve statistics in the gradient computation, we average each observable O, not only using M parallel samples (running each at fixed m), but also in time T thermal history. This means that each 〈O〉 m is estimated by

O m = 1 M M k=1 1 T T t=1 O (k) (t), x (19)
during the simulation. The code we used for the article can be freely downloaded on GitHub.

Learning dynamics and principal components of the dataset

RBMs have been thoroughly studied during the last decade in the context of statistical physics (see Ref. [START_REF] Decelle | Restricted boltzmann machine: Recent advances and meanfield theory[END_REF] for a recent review on the topic) and very important progresses have been made related to the understanding of their thermodynamics and learning mechanisms [START_REF] Tubiana | Emergence of compositional representations in restricted boltzmann machines[END_REF][START_REF] Huang | Statistical mechanics of unsupervised feature learning in a restricted boltzmann machine with binary synapses[END_REF][START_REF] Barra | Phase diagram of restricted boltzmann machines and generalized hopfield networks with arbitrary priors[END_REF][START_REF] Harsh | place-cell'emergence and learning of invariant data with restricted boltzmann machines: breaking and dynamical restoration of continuous symmetries in the weight space[END_REF].

In particular, we now know that the RBMs pattern encoding process is triggered by learning the PCA of the dataset, i.e. the directions {ω α } introduced in Sect. 3. At the beginning of the learning, when the machine is in a paramagnetic state [START_REF] Decelle | Thermodynamics of restricted boltzmann machines and related learning dynamics[END_REF][START_REF] Decelle | Spectral dynamics of learning in restricted boltzmann machines[END_REF][START_REF] Ichikawa | Statistical-mechanical study of deep boltzmann machine given weight parameters after training by singular value decomposition[END_REF], the learning dynamics can be easily understood when decomposing the RBM's weight matrix onto its singular value decomposition (SVD), w = UΣV T , and projecting the gradient equations 3-5 onto each of these modes. In the regime where all the elements of w are small, it can be shown that the first columns of V progressively align with the principal components of the dataset one by one, a process that can be easily followed through the emergence of new eigenvalues in the SVD. At larger learning times, the linear approximation of the gradient is no longer valid and the relation between the dataset's principal directions and the ones of the weight matrix is no longer true. This means that in the early times of the learning process, the RBM passes from being in a paramagnetic phase to a ferromagnetic phase, with non-zero magnetizations m α (recall (6)) along a growing number of directions α, related to the number of the dataset's principal components encoded in the SVD rotation matrix V and the number of hidden nodes of the RBM. We also know that the relevant directions to project the data are for some time {ω α }, and therefore related to the eigenmodes of the weight matrix. These theoretical results suggest that m α must be good order parameters to control thermalisation (and to break metastabilities) at least for some time. The main drawback here, is that the longer the learning, the higher the number of constraints needed to sample ergodically the phase space. In the following, we will use the TMC method to constraint the value of the magnetization along 1 or 2 fixed PCA directions.

Results

In this section, we will compare the quality and the dynamics of RBMs trained using either the standard Glauber MC updates (SMC), or our new TMC sampling strategy to estimate the negative part of the gradient. In both cases, we will always keep the PCD initialisation scheme to enhance as much as possible the thermalisation during the training and perform the same number of Gibbs steps between each parameter update. We further designate SMC-RBM to machines trained with the usual procedure (SMC), and TMC-RMB, to machines trained with TMC. We begin with results on artificial datasets to end with applications in real ones.

Artificial datasets

To illustrate the TMC learning strategy, we will first analyze the learning process of artificial datasets. At a second stage, we will move to real clustered datasets. We will consider 2 different high-dimensional datasets having N v = 1000 visible variables both. The first one, consists of two clusters supported on a one-dimensional subspace given by the vector u = 1/ N v . The second, is made of three separated clusters supported on a two-dimensional subspace. In both cases not only one (or two directions) is sufficient to separate the clusters, but they also lie in a low-dimensional subspace: their extension to the other dimensions is only due to the noise in the creation of the dataset. Both datasets were designed to exaggerate a clustered nature, this means that, points are generated so they form compact clusters separated by large empty regions. In both cases, points can be easily split up in clusters by projecting the dataset on the first (ω 0 , for the 1D case) or the first two (ω 0 and ω 1 for the 2D case) principal components of the dataset.

For this type of datasets, we expect a simple spin-flip type MCMC algorithm to fail dramatically when it comes to sampling ergodically a correctly trained RBM, because jumping from one cluster to another requires surpassing a very large free energy barrier (note that the region between clusters has essentially zero probability by construction). Such a problem is likely to make learning a good RBM on this type of dataset impossible in practice. Indeed, one would need a prohibitive number of MC sweeps to compute the negative term of the gradient, and the chains would never reach equilibrium even with standard methods such as PCD.

One-dimensional dataset

As anticipated, the 2 trained RBMs cannot be sampled ergodically using SMC moves as we illustrate in Fig. 4-A: simulations get trapped in one of the two metastable states from the very first MC sweeps and remain there at least until 10 5 MC sweeps. With this dataset, the learning of the SMC-RBM is very unstable, therefore we used different learning rate for the two RBMs: 10 -2 for the TMC-RBM and a much slower one, 10 -4 , for SMC-RBM. At this point, it is important to stress that one manages to generate points belonging to two different clusters only thanks to the random initialization of the Markov chains. In particular, we initialize each v i as {0, 1} with 50% of probability, which means that initial configurations have in average m 0 ∼ 0.5, which depending on the profile of the potential Ω, tend trap the chains in the nearest free-energy minima, independently on its relative depth. Despite both RBMs generating data in two separate clusters, such clusters are not the same for the 2 machines. This is much clearer if we compare the distribution of the magnetizations m 0 of the generated samples after 10 5 SMC sweeps with each RBM, with those of the dataset, see Fig. 4-B2. Clearly the SMC-RBM generates samples in the wrong region (with wrong values of m 0 ) and this situation does not improve at all if we continue the learning for much longer. On the contrary, the samples generated with the TMC-RBM are at the right position. On Fig. 4-B2, we see only one (red) peak since the SMC is incapable of jumping from one cluster to the other, still the second peak would match the position of the dataset as can be checked on the panel A of the same figure.

As justified above, we can efficiently sample the phase-space of these two trained machines using the TMC sampling method (using simulations at fixed m 0 ). We show the results in plain lines in Fig. 4-B In A we show the evolution of m 0 (t) for 1000 independent simulations initialized at random (in blue we show the samplings of SMC-RBM, in red of the TMC-RBM). Fig. A shows that ergodicity is completely broken: Markov chains remain trapped in the state closest to their random initialization, and not a single jump from one state to the other is ever observed. This means that thermalization is not even close, and it is impossible to properly estimate the relative height of the two peaks from SMC simulations. In B, we instead take both machines under the microscope using the TMC algorithm to properly explore the equilibrium measure.

The TMC results are shown in thick solid lines. In B-1 we show the TMC effective potential Ω( m), which has 2 relative minima for both models. We mark the position of these minima with large dots. We can easily see that they coincide with the mean projection of the two states we sampled with SMC in A. However, in the case of the SMC-RBM, the left minima is much lower, which translates into a much higher TMC probability once p( m) is computed in B-2. It is clear that none of the reconstructed equilibrium probabilities match the histograms obtained with the non-ergodic SMC sampling. Fig. 4-A, but their relative depth is significantly different for the SMC-RBM model. In the TMC-RBM case, the two minima have the same potential's value, meaning that they correspond to two states equally probable: the equilibrium distribution of the TMC-RBM match remarkably well with the dataset's distribution. Instead, in the case of the SMC-RBM, one minima is much higher than the other, and in fact, it is essentially suppressed in the equilibrium distribution p( m0 ), see Fig. 4-B2. If we could wait long enough, all generated samples of the SMC-RBM would fall to the lower-m 0 minimum.

The conclusions of this simple experiment are rather worrying. First, when the clusters are well-separated in the learned RBM, the persistent chains of the classical training method remain trapped in one of the two clusters. The problem, here, is that it is impossible that the parallel chains are correctly balanced between the two minima of the potential. Still, since they can not escape from their local minimum, it gives the impression that the RBM is learning correctly the dataset while in reality the true equilibrium probability is not well-adjusted. This also means that the gradients cannot correct this problem. In fact, at longer training times, we observe that the SMC-RBM is not adjusting the location of the peaks correctly because of this effect.

We can illustrate this effect studying several intermediate models in the training of the SMC-RBM model, see Fig. 5-C where t age refers to the number of parameters updates performed. In this figure, we compare: (i) the equilibrium constrained measure obtained with the TMC method, and the histogram of the m 0 of the samples either in (ii) the persistent chain, or (iii) obtained with a long SMC sampling. We also include, next to each distribution, the trajectory followed by several independent Markov chains during the SMC sampling. We can see that at, and below t age = 91, each Markov chain jumps several times from one cluster to another in 5000 MC sweeps, and because of that the TMC p( m) matches very well the distribution of the persistent and generated samples, even when the distribution is bimodal. As the training progresses, for instance at t age = 96 these jumps get rarer and rarer, which means that the persistent chains get trapped in each of the states with the proportion reached the last time the standard sampling was able of thermalising. Such proportion remains fixed for all the rest of the training, even if the equilibrium distribution of the model is actually unimodal. Neither the samples generated with standard MC moves are reliable.

Fig. 5-C let us illustrate another troublesome effect. Since the distribution of the persistent or generated chains differ strongly from that of the equilibrium, the negative term of the gradient is very badly computed and the learning trajectory is strongly affected by this. In fact, since the persistent chain is almost correctly located, the gradient is probably only adjusting little details of the parameters, while in reality, the equilibrium samples of the true model do not correspond to these samples and therefore the RBM end up very badly trained. Indeed the normal evolution of the training moves towards increasing the relative probability of the low m 0 state (the persistent chain contained an over-representation of high m 0 samples), which only contributes to suppress more and more the probability of the second state.

A sign of this training problem can be easily detected in this simple dataset by monitoring the eigenvalues of the weight matrix during the learning, λ α . In this simple dataset, we observe that only the first eigenvalue λ 0 is turned on during the training. This λ 0 is associated to an eigenvector, ω W 0 , that in both kinds of machines, the TMC-RBM and the SMC-RBM, aligns perfectly with the first principal component of the dataset, ω PCA 0 in Fig. 5-A. Yet, the value of λ 0 does differ from one training to another. Indeed, for the TMC-RBM, the learning stops at a given moment and the first eigenvalue remains constant from then and on. While in the SMC-RBM scheme, the λ 0 keeps growing as the gradient is not computed correctly. This is translated in a decrease of the effective temperature of the model (and a sharpening of the equilibrium peak that is appreciated in the compact trajectories of the t age = 481 RBM in Fig. 5-C).

Two-dimensional dataset

We demonstrate again the feasibility of our approach for a case where the clusters span a twodimensional subspace (see orange dots in Fig. 6). As in the 1D case, a SMC sampling fails to jump from one cluster to another one.

In this case, we rely on a 2D version of the TMC method. We need to constrain the magnetization along two principal directions because a single constraint is not able to split up the 3 clusters because they superimpose when projected on a single direction. A TMC-1D strategy would then suffer from convergence problems too (as the SMC). In this case, we show directly the results of the TMC approach, since it is clear from the previous example that the usual training algorithm cannot work. We show on Fig. 6 the TMC effective potential Ω extracted at various steps of the learning together with the projection of the dataset along the two chosen directions. On the last panel of Fig. 6 (bottom-right) we plot the learning curves of the eigenvalues of the weight matrix. We see that, once the first two eigenvalues are learned, the machine stops evolving on linear timescale. Still, some more eigenvalues are increasing at large number of epoch indicating that we should probably lower the learning rate to stabilize the learning. We also show the evolution of the effective potential at different instants during the training process. We can appreciate how the potential adjust to fit the 3 different minima, and also how the final TMC reconstructed probability matches perfectly the distribution of the dataset in the t age = 1001 figure. This shows how this strategy works perfectly on these difficult datasets.

Real datasets

We now apply our TMC training method on real datasets. It is clear that big improvements are only expected in datasets that are clustered in some way. We therefore focus on two cases. First on the MNIST dataset reduced to the digits 0 and 1. It has been observed recently, that SMC samplings struggle to jump from one digit to another, making convergent generation samplings extremely long [START_REF] Roussel | Barriers and dynamical paths in alternating gibbs sampling of restricted boltzmann machines[END_REF]. The second example is taken from the 1000 Human Genome Project [START_REF] Consortium | A global reference for human genetic variation[END_REF], already studied in a generative context recently [START_REF] Yelmen | Creating artificial human genomes using generative neural networks[END_REF][START_REF] Decelle | Equilibrium and non-equilibrium regimes in the learning of restricted boltzmann machines[END_REF]. In this last dataset, a strong clustered structure is observed when the data are projected along the first principal components of the dataset, {ω α }.

MNIST 0-1 dataset

This dataset was created by taking only the 0 and 1 digits from the MNIST dataset [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF]. We obtain slightly more than 10 4 images. Unlike the full dataset, this reduced dataset splits into two well-separated clusters when projected along ω 0 (the first principal direction of the dataset centered on zero), see Fig. 7.

In this example, we will again compare the performance of a TMC-RBM with that of a SMC-RBM. A critical moment in training this dataset is when the first principal direction is learned. In this regime, the phase space of the RBM consists of two separate modes where an SMC dynamics has trouble jumping from one to the other, while TMC should have no problem. Later in training, as more and more directions are learned, we will see that the strong acceleration of TMC is damped because the trained RBM seems to create alternative dynamical paths to jump from 0s to 1s that do not have to overcome this large barrier in m 0 . We can characterize these dynamical effects by computing the integrated autocorrelation times [START_REF]Monte carlo methods in statistical mechanics: foundations and new algorithms[END_REF][START_REF] Amit | Field theory, the renormalization group, and critical phenomena: graphs to computers[END_REF], τ int , of the evolution of m 0 (t) in the two different sampling dynamics:

τ int = 1 2 + κτ int t=1 ρ m (t) ρ m (0) , with ρ m (t) = 1 M -t M -t s [m 0 (s) -m 0 ] [m 0 (s + t) -m 0 ] ,
where κ is a small number, in our case κ = 6 as in [START_REF]Monte carlo methods in statistical mechanics: foundations and new algorithms[END_REF][START_REF] Amit | Field theory, the renormalization group, and critical phenomena: graphs to computers[END_REF], and m 0 is the empirical estimator of the equilibrium 〈m 0 〉 using M independent samples.

In practice, we find that the relaxation time of the SMC dynamics is enormous (not even measurable in "reasonable" times and higher than 10 7 ) when only one direction was learned, while the TMC relaxes in only a few sweeps. We estimate the SMC relaxation time in each RBM in Fig. 8-B. This divergence of times is related to the apparition of two very distant metastable states, see Fig. 8-A, just as we discussed in the context of the artificial 1D dataset. Then, in a second learning phase, the SMC relaxation time falls back to measurable values using the transverse directions learned by the RBMs. Even though SMC thermalization times remain quite high, SMC sampling becomes competitive with TMC sampling in terms of physical computation times because TMC relaxation times grow rapidly when m 0 is no longer able to break metastabilities, and because TMC sampling is much slower than SMC because the update of visible variables cannot be parallelized. Nevertheless, quite unexpectedly, we find that TMC-RBMs have much faster SMC relaxational dynamics than SMC-RBMs and that the factor between the two increases with the number of parameter updates, see Fig. 8-B.

When we compare the samples generated with the two learned RBMs at t age = 1001 parameter updates, we find that those obtained with the TMC-RBM are qualitatively better than those from the SMC-RBM. However, we emphasize that this comparison is difficult to see from a direct visual inspection, see Fig. 9A, since the difference is not related to the overall quality or definition of the digits, but rather to the ratio of 1s to 0s and the variability of the shape of

A B

Figure 8: A: Histogram of the dataset along the first direction ω 0 , together with the tethered potential Ω( m0 ) and the corresponding reconstructed TMC probability p( m0 ) of both RBMs in the early stages of learning (here t age = 61). We clearly observe the formation of two metastable states in both machines. In practice, we observe that Markov chains using SMC dynamics are trapped in one of these metastable states since initialization and cannot visit the rest of the phase space. Exactly as described in the synthetic 1D dataset. B: We estimate the SMC integrated autocorrelation time τ SMC int of m 0 (t) as a function of the number of training updates used to train the two machines: TMC-RBM and SMC-RBM. The shaded area marks the region where τ SMC int exceeds 10 7 MC sweeps, while τ TMC int (for the slower m constraints ) remains below 15 MC sweeps. When t age 10 2 , the SMC relaxation time drops to feasible values in both machines. Nevertheless, the relaxation dynamics of TMC-RBM is significantly faster than that of SMC-RBM. the generated digits. The comparison becomes clearer if we project the generated data onto the first 4 principal components and compare the histograms obtained with those obtained with the MNIST projections, see Fig. 9C. It is clear that the equilibrium distribution of TMC-RBM reproduces much better the original. Another interesting effect is that the typical time needed to jump from one digit to another (from one cluster to another) during the generation sampling is faster for the TMC-RBM, as illustrated in Fig. 9B.

The reason why TMC-RBMs present faster relaxations than SMC-RBMs is not clear, but we think it is related to the effects of the nonequilibrium regime that were encoded in the model during training when the mixing times of the model were decades above the number of Gibbs steps used to train the machine, as estimated in Fig. 8B. Indeed, in Ref. [START_REF] Decelle | Equilibrium and non-equilibrium regimes in the learning of restricted boltzmann machines[END_REF] we observed that one of the consequences of RBMs trained deeply in non-equilibrium was precisely their extremely slow relaxation. This effect is also probably related to the degradation phase observed in Ref. [START_REF] Dabelow | Three learning stages and accuracy-efficiency tradeoff of restricted boltzmann machines[END_REF].

Human genome dataset

We considered the population genetics' dataset of Ref. [START_REF] Consortium | A global reference for human genetic variation[END_REF][START_REF] Colonna | Human genomic regions with exceptionally high levels of population differentiation identified from 911 whole-genome sequences[END_REF]. This dataset corresponds to a sub-part of the genome of a population of 2504 individuals (each providing two strands of DNA, yielding a dataset of 5008 samples) sampled from 26 populations in Africa, East Asia, South Asia, Europe, and the Americas. Each variable is binary {0, 1} and indicates whether or not a particular gene is altered compared to the human reference genome. A total of 805 alterations are reported, reflecting a high proportion of the population structure present in the whole genome dataset. The dataset is therefore made of 805 binary variables and 5008 independent samples. For our study, the dataset was made of 4500 samples for the training set. This dataset has a very clustered nature, as can be appreciated in the dataset's projection along the second and third principal directions shown in Fig. 10 where each dataset sample is represented as a black dot. One can also see that SMC dynamics get completely trapped in these different clusters in a well-trained machine, as we illustrate with different color trajectories each corresponding to an independent sampling of the RBM. We therefore constrain these directions for the training of the RBM using the TMC method. For this case, we will exhibit two clear advantages of this method. First, we show the feasibility of training the RBM on a real dataset using a set of two constrains. In this setting, the TMC method is not only useful to estimate better the negative term of the gradient, but it also gives the possibility to monitor how-well the equilibrium measure of the RBM is matching the data through the TMC effective potential. It is worth of mentioning that evaluating the quality of an RBM (that is, ensuring a proper adjustment of the density peaks with those of the dataset) with traditional SMC sampling methods is extremely hard, because jumps between the different metastable clusters are very rare and there are four of such clusters involved.

In fact, on Fig. 11 we show the TMC effective potential and the reconstructed distribution at various stages of the learning. It is interesting to see that it seems quite difficult to adjust correctly the 4 clusters to their correct weight. Thanks to the TMC method, we can easily track the position of the RBM's clusters and diminish the learning rate to fine-tune the learning. On Fig. 11, we can see that after t age ≈ 1200, the TMC-RBM seems to be well-learned.

The second crucial improvement is that, only with the TMC method we can guarantee that all the clusters are visited regularly. Indeed, we already showed in Fig. 10 that we suffered strong break of ergocity issues even after 10 5 SMC sweeps. This fact tells us that, on this dataset, each persistent chain of the PCD will remain trapped in certain regions of the phase space and will not be representative of the equilibrium measure of the RBM. As a consequence, the persistent chains will not estimate correctly the negative term of the gradient.

Finally, it is possible to try to train this dataset using the classical PCD -SMC algorithm. However, what we see is that the learning dynamics has a lot of trouble adjusting correctly the different minima, and that in addition, even if one can "cherry pick" the exact number of updates where the trained SMC-RBM has all the good dataset minima, such a process requires using the TMC method to track the presence of the good minima in the effective potential because SMC samplings are completely stuck. Even in that case, we observe that again TMC-RBMs relax faster than the equivalent SMC-RBM (for the same t age ). Yet, SMC dynamics are that slow for both kinds of machines, that computing τ SMC int , as done in Fig. 8 for MNIST 0-1s, is beyond our numerical capacities. Still we can compare for both machines, the relaxation time in the TMC runs (which are much shorter). We obtain that τ

TMC int [SMC -RBM] ∼ 4 • τ TMC int [TMC -RBM],
showing the TMC-RBM is again faster that the PCD one.

Runtime comparison

The TMC method has two major disadvantages. First, the global constraint imposed on the measure breaks the conditional independence of the two layers: If the constraint applies to the visible nodes, the hidden nodes are still independent given a set of visible nodes, but the converse is no longer true. Because of the constraint, the visible nodes are no longer independent once the hidden nodes are fixed, since there is an explicit interaction between the visible nodes. Therefore, we cannot create a new state for the visible nodes in parallel, which means that part of the advantage (speed increase) of using GPUs is lost. Another disadvantage of the TMC method is that we no longer have to calculate the correlations between the visible and hidden nodes for a single RBM, but for as many discretization step RBMs as we used to constrain the order parameters. However, once the machine is learned, the time required to generate samples with TMC is significantly reduced compared to SMC. For each generated for RBMs trained on either CPU/GPU and using the TMC or SMC method during learning. The curves were generated with a Geforce RTX 3090 as GPU and an AMD Ryzen 9 5950X as CPU. data point, the process is as follows: first, a TMC constraint m is randomly drawn from the reconstructed probability p( m), and then a TMC sampling process is performed at fixed m, a dynamical MC process that converges to equilibrium in a very short time thanks to the constrained ensemble. In Fig. 12 we show how the runtimes compare when training the RBM with TMC or SMC and on CPU or GPU. The SBM-GPU is much faster than all other implementations, making it the perfect candidate when the phase space if not clustered. Still, when an energy barrier is present, the only mean to reach thermalization in a large system is to use the TMC approach. As future improvement, it shall be possible to tether the hidden magnetization (rather than the visible one). In that case, the performance bottleneck would be reduced since the number of hidden nodes is generally much smaller than the number of visible ones, in particular in low-dimensional datasets.

Conclusion

We have demonstrated the feasibility and benefits of using a biased Monte Carlo sampling strategy for both training -i.e., to better estimate the negative terms of the gradient in order to update the parameters of the model -and for generating new samples in trained RBMs. The largest improvements are observed for highly clustered datasets, where classical local MCMC moves get stuck in one or several of these clusters and fail to sample ergodically the entire phase space (in reasonable time). We show that the Tethered Monte Carlo method is extremely efficient for sampling the equilibrium measure of these types of models, as long as one can find good observables to uniquely identify each cluster. We show that the projection along the first PCA directions in the early stages of training are good order parameters as motivated from the mean-field results [START_REF] Decelle | Thermodynamics of restricted boltzmann machines and related learning dynamics[END_REF][START_REF] Decelle | Spectral dynamics of learning in restricted boltzmann machines[END_REF]. We show that this TMC strategy is crucial to correctly learn artificial low-dimensional clustered datasets. On these datasets, the standard learning approach that combines parallel local MCMC updates with PCD completely fails to obtain reliable models for the data. On real datasets, we find that as learning progresses, more ordering parameters are needed to avoid the presence of metastabilities during sampling, which means that the TMC strategy with only one or two constraints also starts to suffer from ergodicity problems. We note, however, that even if TMC does not thermalize easily, the models trained with TMC exhibit significantly faster relaxation dynamics than the models obtained with the standard approach.

Moreover, we find that even for RBMs trained with standard sampling techniques, the TMC sampling approach is extremely useful to assess the quality of a given trained machine, since it allows a direct calculation of the probability distribution of the model, projected along a set of selected directions, thus providing the possibility to estimate how well the distribution of the model matches the empirical one. The developed approach is thus a practical alternative to the most commonly used Gibbs sampling algorithm when the properties of the dataset drive learning toward a model with well-separated modes.

Despite the impressive performance of the TMC model on clustered datasets, it has two important drawbacks: First, it is much slower than the traditional alternative sampling method in terms of computation time because the global constraint introduces correlations between the visible (or hidden) variables, thus removing the conditional independence of the two layers that favors parallelization. For this reason, alternative approaches involving only local constraints should be considered in the future to restore parallel performance. A possible improvement could be to constrain the value of some hidden nodes with a condition that does not remove the conditional independence. However, a clever mechanism must be found to select the hidden nodes to be fixed. Second, the dimensionality of the constraints. Since we need to integrate the potential exactly, we need to discretize the space into as many dimensions as the number of constraints allows. This makes the method cumbersome even with three constraints and probably impossible with more.

1 ,

 1 and the distribution p(m (m)

1 )Figure 2 :

 12 Figure 2: We describe the basic steps of the TMC sampling procedure, which can be used either to generate equilibrium samples of the model or to estimate the loglikelihood gradient during training.

Figure 3 :

 3 Figure 3: Reconstruction of the probability distribution profile of the dataset using TMC sampling. In red we show the derivative of the TMC effective potential, Ω ( m), obtained from TMC simulations at fixed values of m and α = 2 • 10 4 .In black, we show the Tethered probability p( m) obtained from the numerical integration of Ω( m) plotted against 〈m〉 m. We compare this probability with the histogram of the projection of the dataset along ω 0 . An almost perfect match is seen, even though the two probabilities should only match in the limit α → ∞.

- 1 B- 2 Figure 4 :

 124 Figure4: We compare RBMs trained with two different sampling procedures. We analyze models trained for 481 parameter updates and a learning rate of 10 -2 for the TMC-RBM and 10 -4 for the SMC-RBM. The SMC-RBM is trained with PCD and 10 MC sweeps of Glauber dynamics (the standard procedure, SMC), and the TMC-RBM is trained with PCD and 10 TMC1D sweeps fixing the projection m 0 . Once both machines are trained, we can sample their configuration space with SMC moves in A or TMC moves in B. In A we show the evolution of m 0 (t) for 1000 independent simulations initialized at random (in blue we show the samplings of SMC-RBM, in red of the TMC-RBM). Fig.Ashows that ergodicity is completely broken: Markov chains remain trapped in the state closest to their random initialization, and not a single jump from one state to the other is ever observed. This means that thermalization is not even close, and it is impossible to properly estimate the relative height of the two peaks from SMC simulations. In B, we instead take both machines under the microscope using the TMC algorithm to properly explore the equilibrium measure. The TMC results are shown in thick solid lines. In B-1 we show the TMC effective potential Ω( m), which has 2 relative minima for both models. We mark the position of these minima with large dots. We can easily see that they coincide with the mean projection of the two states we sampled with SMC in A. However, in the case of the SMC-RBM, the left minima is much lower, which translates into a much higher TMC probability once p( m) is computed in B-2. It is clear that none of the reconstructed equilibrium probabilities match the histograms obtained with the non-ergodic SMC sampling.

Figure 5 :

 5 Figure 5: In A we show the projection of the first eigenvector of the RBM weight matrix W , ω W 0 , onto the direction of the first principal component of the dataset ω PCA 0 , as a function of the number of parameter updates t age , for a TMC-RBM and a SMC-RBM. In B we show the evolution of the first 10 eigenvalues of W for the two machines. Only the first eigenvalue is expressed, since only one direction is needed for this dataset: the others are almost zero. In C we examine the SMC-RBM at 7 different stages of training (t age = 71, 81, 91, 96, 111, 221, 481 updates). In particular, in the right panel of C we compare the equilibrium distribution of m 0 obtained with the TMC method p( m), the histogram of the dataset, the histogram of the persistent chain during learning, and the histogram of the samples generated with 50000 SMC sweeps. To illustrate the trapping phenomena, we show in the left panel of C the trajectory of several Markov chains (in the m 0 -m 1 space) obtained during the SMC generation process (in color) and compare them with the scatter plot of the dataset projections (black dots).

Figure 6 :

 6 Figure 6: A-C: In orange dots, we show the m 0 and m 1 projections of the samples of a synthetic dataset of 3 clusters. In a contour map, we show the TMC effective potential Ω reconstructed at different moments of the training of a TMC-RBM (t age refers to the number of updates). At t age = 1001, we also show the TMC probability distribution (in a heatmap), which correctly reflects the density of the three clusters. In D we show the evolution of the eigenvalues of the W matrix as a function of t age .

Figure 7 :

 7 Figure 7: Histogram of the projection m 0 of all samples of the dataset MNIST 0-1 along its first principal component. We show in blue and orange the histogram of the magnetizations of only the 0 and 1 digits. It is clear that m 0 efficiently separates both digits.

Figure 9 :

 9 Figure 9: For a TMC-RBM and a SMC-RBM, both trained for 1001 parameter updates with the same learning rate γ = 0.01 and the same k = 10, we show: A: 50 equilibrium samples generated (with each model) and compared to 50 samples of the original MNIST dataset. In order to ensure equilibration, we performed 10 7 SMC sweeps prior to imaging B: For both machines, we show the scatter plot of the SMC dynamical trajectories projected along the first two directions: m 1 (t), m 0 (t), all 10 3 sweeps, for 5 • 10 5 SMC sweeps. It can be clearly seen that the sampling of the SMC-RBM is slower, since there are fewer jumps (red curve) between digits 0 and 1. C: The histogram of the data generated by the two RBMs (after 10 7 SMC steps), in red by the SMC-RBM and in blue by the TMC-RBM, projected along the first four principal directions of the dataset. The black line corresponds to the histogram of the dataset. It is clear that the SMC-RBM has many more anomalies than the TMC.

Figure 10 :

 10 Figure 10: We show the training dataset projected along the second and third principal directions (namely m 1 and m 2 ) in black dots. The samples are grouped in separate clusters. In the colored solid lines, we show different SMC sampling trajectories {m 1 (t), m 2 (t)} t of 10 5 MC sweeps extracted from a TMC-RBM trained for t age = 1000 updates, λ = 0.01 and k = 10. One can clearly appreciate that there is no jump between the high and low m 1 sectors during this period.

Figure 11 :

 11 Figure11: Evolution of the potential along the second and third principal directions during the learning of the TMC-RBM. We clearly observe how the second direction is learned by the RBM, followed by the third. After one thousand parameter updates, the reconstructed probability distribution covers the dataset.

Figure 12 :

 12 Figure 12: We show the average wall time per epoch (averaged over 20 epochs)for RBMs trained on either CPU/GPU and using the TMC or SMC method during learning. The curves were generated with a Geforce RTX 3090 as GPU and an AMD Ryzen 9 5950X as CPU.

(m(v )m) 2=1,[START_REF] Smolensky | Information processing in dynamical systems: Foundations of harmony theory[END_REF] technically the global constraint could perfectly concern both the visible and the hidden units, m(h, v )∼ m, or only the hidden m(h)∼ m too.

actually p( m) versus 〈m 0 (v )〉 m for the whole range of m simulated

Acknowledgments

This work was supported by the Comunidad de Madrid and the Complutense University of Madrid (Spain) through the Atracción de Talento programs (Refs. 2019-T1/TIC-13298 for A.D. and 2019-T1/TIC-12776 for B.S.), the Banco Santander and the UCM (grant PR44/21-29937); and Ministerio de Economía y Competitividad, Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER) (Spain and European Union) through the grants. PID2021-125506NA-I00 and PGC2018-094684-B-C21. N.B. would like to thank the Department of Theoretical Physics of the Complutense University of Madrid for the kind reception and support during his stay.


Ricardo Amils

Muriel Gargaud

José Cernicharo Quintanilla

Henderson James Cleaves

William M Irvine

Daniele Pinti

Michel Viso

Franck Poitrasson
email: franck.poitrasson@cnrs.fr

Iron Isotopes

Keywords:

Iron has four naturally occurring stable isotopes: 54 Fe, 56 Fe, 57 Fe, and 58 Fe of respective abundances 5.80 %, 91.72 %, 2.20 %, and 0.28 %. They are increasingly used to determine the source of the iron present in geological or biological materials and/or to characterize the chemical reactions that may have lead to the form of iron under study. The isotopic signature of Fe is particularly sensitive to redox effects, including certain biological reactions such as bacterial dissimilatory iron reduction (DIR).

Overview

Although less variable in abundances than carbon or oxygen stable isotopes, iron isotopes show analytically significant mass-dependent isotopic variations in nature. The geochemical community has investigated the range of Fe isotope compositions in both high-and low-temperature geological environments. Contrary to early expectations, there are no specific biological-type iron isotope signatures. On the other hand, the largest isotopic fractionation was observed for low-temperature redox reactions, and biological dissimilatory iron reduction (DIR) is one of those biological reactions that affect significantly iron isotope compositions. Another requirement for inorganic or organic reactions to leave an isotopic imprint is that the reaction exchange is not complete, that is, not all of the original pool of iron has been consumed. Nature-based and experimentally-based astrobiology studies have shown that iron isotopes are particularly useful to study traces of life and/or their environment when used in conjunction with other mineralogical, biological, and chemical indicators.

Basic Methodology

The isotopic composition of iron in natural samples has been measured since the 1940s by electron impact ionization and, more recently, by thermal ionization mass spectrometry (TIMS) after sample dissolution and chemical purification using anion exchange chromatography. The level of precision reached relative to natural variations rendered the study of Fe isotope signatures of little interest until the end of the twentieth century. The advent of the second generation of multiple collector inductively coupled plasma mass spectrometers (MC-ICP-MS) led to an order of magnitude improvement in analytical precision [START_REF] Belshaw | High precision measurement of iron isotopes by plasma source mass spectrometry[END_REF]. This opened the way to an exponential development of these "non-traditional isotopes," such as iron isotopes. Although the purification chemistry is straightforward and has been know for some time [START_REF] Strelow | Improved separation of iron from copper and other elements by anion-exchange chromatography on a 4% cross-linked resin with high concentrations of hydrochloric acid[END_REF], MC-ICP-MS techniques vary among laboratories. Mass bias corrections range from the simple sample-standard bracketing approach to more sophisticated Cu-or Ni-doping methods. A few groups have also adopted the double spike technique, which remains the only way to determine precise Fe isotope composition by TIMS. Isotopic results are conventionally reported using the delta notation with reference to the IRMM 14 reference material (unfortunately no longer available from the European Institute for Reference Materials and Measurements, Geel, Belgium, but it may be conveniently replaced by the IRMM 524A reference material) that has a composition close to that of chondrites. For example:

δ 57 Fe (‰) € = 57 Fe/ 54 Fe sample 57 Fe/ 54 Fe IRMM 14 -1 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ×10 3
and this notation can be equally defined for δ 56 Fe. Some groups still report their Fe isotope compositions with reference to average of igneous rocks. Although being isotopically homogeneous, this convention should be abandoned since some igneous rocks depart significantly from this average value (e.g., some granites from the Earth's continental crust; [START_REF] Poitrasson | Heavy iron isotope composition of granites determined by high resolution MC-ICP-MS[END_REF], and therefore, this may generate some confusion among non-specialists.

Key Research Findings

A first important finding is that the iron isotope compositions of igneous rocks that form the Earth's crust vary very little. As a result, they make a baseline against which low-temperature, superficial biogeochemical processes that display larger isotope variations may be compared to [START_REF] Beard | Application of Fe isotopes to tracing the geochemical and biological cycling of Fe[END_REF][START_REF] Poitrasson | On the iron isotope homogeneity level of the continental crust[END_REF]. Various field-based, experimental, and especially theoretical studies conducted over the past two decades revealed that the main driving factors generating Fe isotope variations are redox and coordination chemistry changes. For example, the largest equilibrium Fe isotope fractionation among two coexisting iron reservoirs (~4.5 ‰ in δ 57 Fe) was observed between aqueous Fe 2+ and Fe 3+ species [START_REF] Welch | Kinetic and equilibrium Fe isotope fractionation between aqueous Fe(II) and Fe(III)[END_REF][START_REF] Anbar | Theoretical investigation of iron isotope fractionation between Fe(H 2 O) 6 3+ and Fe(H 2 O) 6 2+ : implications for iron stable isotope geochemistry[END_REF]. The importance of the bonding environment on Fe isotope signatures is best illustrated by theoretical prediction of Fe isotope fractionation between mineral phases holding iron in the same redox state, but with different network-forming anions, such as in sulfide and carbonate minerals [START_REF] Polyakov | The use of Mössbauer spectroscopy in stable isotope geochemistry[END_REF][START_REF] Blanchard | Iron isotope fractionation between pyrite (FeS 2 ), hematite (Fe 2 O 3 ) and siderite (FeCO 3 ): a first-principles density functional theory study[END_REF]. Ab initio computations also illustrate the importance of iron coordination chemistry and nature of aqueous ligands to interpret the Fe isotope composition of environmental samples (Fuji et al., 2014). The finding that mass-dependent Fe isotope variation occurs in nature is interesting in itself. Even more important however is the discovery that for the same rocks, the stable isotopes of iron show contrasted fractionation when compared to those of carbon [START_REF] Matthews | Controls on ironisotope fractionation in organic-rich sediments (Kimmeridge Clay, Upper Jurassic, southern England)[END_REF]. Although not unexpected given the contrasted interatomic bonding environments of Fe and C, this illustrates that the stable isotopes of different elements can trace different processes and/or different properties of the studied reservoirs. As such, combining different isotopic tracers is certainly a promising approach for future research, notably to detect biological activity [START_REF] Johnson | The iron isotope fingerprints of redox and biogeochemical cycling in the modern and ancient Earth[END_REF][START_REF] Severmann | Reconstructing paleoredox conditions through a Multitracer approach: the key to the past is the present[END_REF].

Applications

Since astrobiology aims at looking for traces of life on the early Earth and on planets other than the Earth, one would ideally do precise Fe isotope measurements of extraterrestrial samples to determine signatures that could be tracking life or at least environmental conditions suitable for to life. However, the need for high-precision mass spectrometry methods that were established in terrestrial laboratories over the past decade makes such delicate measurements on robotic missions on other planets impossible at present. Aside for the Moon that has never been prone to the development of life since the beginning of its history given its extreme temperatures, lack of liquid water, and lack of atmosphere, we should wait for sample return missions from other planetary bodies. In the meantime, we have to rely on meteorites thought to come from planets where life may have existed, such as Mars. However, Martian meteorites that belong to the Shergotty, Nakhla, and Chassigny (SNC) class are essentially igneous rocks. Hence, looking for traces of life on those rocks using iron isotopes was unsuccessful [START_REF] Anand | Searching for signatures of life on Mars: an Fe-isotope perspective[END_REF]. This was not unexpected given that even if traces of life were searched on bulk terrestrial fresh igneous rocks, the igneous rock Fe isotope signature would have likely overwhelmed the biologically processed Fe. Another limitation is that despite initial expectation [START_REF] Beard | High precision iron isotope measurements of terrestrial and lunar materials[END_REF], there are, like for carbon isotopes, no unique Fe isotope biological signatures [START_REF] Anbar | Nonbiological fractionation of iron isotopes[END_REF]. Indeed, inorganic processes may also produce strong Fe isotope fractionation leading to light products, notably when kinetic effects are involved. Hence, the use of Fe isotopes to search for traces of biological activity should necessarily be conducted in conjunction with other mineralogical, chemical, and isotopic indicators.

Astrobiology also looks at extreme terrestrial environments where life may occur despite hostile conditions. This is the case of acid mine drainage waters from the Rio-Tinto-Odiel sulfide deposit basin, southern Spain, that is characterized by anomalously acidic waters (1.5< pH <3.9) containing elevated concentrations of dissolved Fe (>20 000 ppm). Despite these conditions, bacterial activity occurs. This geological setting is of interest to astrobiologists because it forms specific minerals including jarosite that also occur at the Martian Meridiani Planum site. Iron isotopes were used by [START_REF] Egal | Iron isotopes in acid mine waters and iron-rich solids from the Tinto-Odiel Basin (Iberian Pyrite Belt, Southwest Spain)[END_REF] to study this Spanish site. Although the biological activity imprint could not be detected through Fe isotope compositions, this study combined mineralogical, aqueous geochemistry, and iron isotope determinations to put constraints on the mechanisms that led to the formation of such an extreme geochemical environment on Earth.

Archean apex cherts from Marble Bar, Western Australia, is another interesting geological case study for astrobiologists since it has been claimed to contain the oldest fossils of microbes found yet on Earth, ca. 3.5 Ga years ago, although this finding generated a lot of debate. [START_REF] Pinti | Biogenic nitrogen and carbon in Fe-Mn-oxyhydroxides from an Archean chert, Marble Bar, Western Australia[END_REF] used iron isotopes to study Fe-Mn-oxyhydroxides containing C and N with a biological-like isotopic signature that occur in these cherts. They found δ 57 Fe values compatible with oceanic hydrothermal precipitation, thereby strengthening the inference that C and N isotope compositions mark 3.5 Ga old biological activity.

Archean banded iron formations (BIFs) are also obvious geological records to look for ancient traces of life. Combining a detailed petrographic study with C, O, and Fe isotope investigations of carbonates from these sedimentary formations led [START_REF] Heimann | Fe, C, and O isotope compositions of banded iron formation carbonates demonstrate a major role for dissimilatory iron reduction in similar to 2.5 Ga marine environments[END_REF] to highlight the importance of bacterial dissimilatory iron reduction back to ca. 2.5 Ga ago. Studying again Fe and C isotopes in magnetites and carbonates from BIFs of ca. 2.5 and 3.8 Ga, [START_REF] Craddock | Iron and carbon isotope evidence for microbial iron respiration throughout the Archean[END_REF] concluded that this bacterial process is in fact recorded in the oldest available sedimentary rocks we have on Earth.

Experimental investigations may also help to study processes that may have occurred on the early Earth. For example, [START_REF] Gronstal | Laboratory experiments on the weathering of iron meteorites and carbonaceous chondrites by ironoxidising bacteria[END_REF] examined experimentally how the first living organisms on Earth, that is, prokaryotes, could have extracted energy by oxidizing metallic iron from meteoritic material. Such a form of Fe, which is nearly nonexistent on the present Earth's surface, was supposedly abundant at that time corresponding to the end of the Late Heavy Bombardment and there was an anoxic atmosphere. Interestingly, these authors found that whereas metal weathering occurred regardless of the presence of iron-oxidizing bacteria, the weathering products clearly showed isotopically contrasted iron pathways during oxidation with and without bacteria. This suggests that the biological action in meteorite weathering could be highlighted with Fe isotopes (Fig. 1). Among subsequent experimental studies, [START_REF] Amor | Mass-dependent and -independent signature of Fe isotopes in magnetotactic bacteria[END_REF] investigated the Fe isotope fractionation during the growth of magnetotactic bacteria and found a unique mass independent iron isotope fractionation between the culture medium and the bacterially-produced magnetites. The next challenge is to find geological formations where this effect can be recorded despite obvious risks of dilution of this isotopic signature in the rock matrix [START_REF] Amor | Iron-biomineralizing organelle in magnetotactic bacteria: function, synthesis and preservation in ancient rock samples[END_REF].

Future Directions

Iron isotopes have not yet been applied frequently to astrobiology problems. More case studies therefore need to be investigated, either experimentally, and from natural examples of extreme terrestrial environments, or of the oldest remains of the early Earth. There are however more fundamental issues that need to be addressed first to be able to fully exploit the potential of iron isotope signatures in this endeavor. In a number of studies, equilibrium isotopic fractionation between phases is taken as a reference against which the observed natural variations may be interpreted. This is especially important in low-temperature environments where kinetic isotope fractionations are more likely to be observed given the sluggishness of chemical reactions under such conditions and therefore the large amount of time required to reach equilibrium. For this, more experimental and numerical determinations of the fractionation factors are required. It is important to conduct both approaches simultaneously since cross-validation of experimental, theoretical, and natural determinations is clearly needed [START_REF] Hill | Effects of changing solution chemistry on Fe3+/Fe2+ isotope fractionation in aqueous Fe-Cl solutions[END_REF]. Another important development is the need for in situ data. Meteorites and very old terrestrial samples have had a protracted history that generated chemical and isotopic heterogeneity. Producing bulk sample data unveils only a portion of the message. Doing precise and accurate in situ Fe isotope determinations [START_REF] Horn | In situ iron isotope ratio determination using UV-femtosecond laser ablation with application to hydrothermal ore formation processes[END_REF][START_REF] Steinhoefel | Deciphering formation processes of banded iron formations from the Transvaal and the Hamersley successions by combined Si and Fe isotope analysis using UV femtosecond laser ablation[END_REF][START_REF] Czaja | Biological Fe oxidation controlled deposition of banded iron formation in the ca. 3770 Ma Isua Supracrustal Belt (West Greenland)[END_REF] is more likely to fully unfold the history that iron isotopes can tell.

See Also

. Archean Traces of Life

. Banded Iron Formation .

Biomarkers

.

Isotope

. Isotopic Fractionation (Interstellar Medium)

.

Jarosite

. Late Heavy Bombardment

. Rio Tinto

Fig. 1 :

 1 Fig.1: Products of iron meteorite weathering in aqueous solutions with (Biomet) and without (Abiomet) iron-oxidizing bacteria. Both the weathered layers at the surface of the meteorites (scraping) and precipitates elsewhere in the experimental cell (precipitate) were sampled. Iron isotope composition relative to the starting meteorite shows contrasted Fe isotope signatures for these weathering products with and without bacteria, suggesting contrasted reaction pathways (Error bars within symbols; From[START_REF] Gronstal | Laboratory experiments on the weathering of iron meteorites and carbonaceous chondrites by ironoxidising bacteria[END_REF]