dataesr/scientific-paragraphs-categorization huggingface dataset:
Nevertheless, RBMs have made a strong comeback in the field of Computational Biology in recent years, especially for purposes of pattern extraction  [8] [9] [10] [11] [12] [13] [14]. 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  [8, [15] [16] [17] [18] , 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.
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. 3. 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  [22]  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. (2). 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  [3] , or its later refinement the (ii) Persistent CD (PCD)  [23].
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  [24]  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  [25] , 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.
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. [24]. 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. [24] , 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.
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  (13)  and to obtain p( m)  (14). 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.
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 pipeline 2  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.
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. (3) (4) (5). 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  (15). 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. [18]  for a recent review on the topic) and very important progresses have been made related to the understanding of their thermodynamics and learning mechanisms  [8, 15, 17, 19]. 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  [16, 32, 33] , 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.
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 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.
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.
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  [34]. The second example is taken from the 1000 Human Genome Project  [35] , already studied in a generative context recently  [14, 20]. 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  [36]. 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.