Aurélien Decelle
email: adecelle@ucm.es

Cyril Furtlehner

Restricted Boltzmann Machine, recent advances and mean-field theory *

Keywords: RBM, machine learning, statistical physics PACS: 02.50.-r, 02.30.Z, 05.70.F

This review deals with Restricted Boltzmann Machine (RBM) under the light of statistical physics. The RBM is a classical family of Machine learning (ML) models which played a central role in the development of deep learning. Viewing it as a Spin Glass model and exhibiting various links with other models of statistical physics, we gather recent results dealing with mean-field theory in this context. First the functioning of the RBM can be analyzed via the phase diagrams obtained for various statistical ensembles of RBM leading in particular to identify a compositional phase where a small number of features or modes are combined to form complex patterns. Then we discuss recent works either able to devise mean-field based learning algorithms; either able to reproduce generic aspects of the learning process from some ensemble dynamics equations or/and from linear stability arguments.

Introduction

During the last decade, machine learning has experienced a rapid development, both in everyday life with the incredible success of image recognition used in various applications, and in research [START_REF] Goodfellow | Deep learning[END_REF][START_REF] Mehta | [END_REF] where many different communities are now involved. This common effort involves fundamental aspects such as why it works or how to build new architectures and at the same time a search for new applications of machine learning to other fields, like for instance improving biomedical images segmentation [START_REF] Ronneberger | International Conference on Medical image computing and computer-assisted intervention[END_REF] or detecting automatically phase transitions in physical system. [START_REF] Carrasquilla | [END_REF] Machine learning classical tasks are divided into at least two big categories: supervised and unsupervised learning (putting aside reinforcement learning and the more recently introduced approach of self-supervised learning). Supervised learning consists in learning a specific task -for instance recognizing an object on an image or a word in a speech-by giving the machine a set of samples together with the correct answer and correcting the prediction of the machine by minimizing a well-design and easy computable loss function. Unsupervised learning consists in learning a representation of the data given an explicit or implicit probability distribution, hence adjusting a likelihood function on the data. In this latter case, no label is assigned to the data and the result depends thus solely on the structure of the considered model and of the dataset.

In this review, we are interested in a particular model: the Restricted Boltzmann Machine (RBM). Originally called Harmonium [START_REF] Smolensky | Information Processing in Dynamical Systems: Foundations of Harmony Theory[END_REF] or product of experts, [START_REF] Hinton | [END_REF] RBMs were designed [7] to perform unsupervised tasks even though they can also be used to accomplish supervised learning in some sense. RBMs are part of what is called generative models which aim to learn a latent representation of the data in order to later be used to generate statistically similar new data -but different from those of the training set. There are Markov Random Fields (or Ising model for physicists), that were designed as a way to automatically interpret an image using a parallel architecture including a direct encoding of the probability of each "hypothesis" (latent description of a small portion of an image). Later on, RBMs started to take an important role in the Machine Learning community, when a simple learning algorithm introduced by Hinton et al., [START_REF] Hinton | [END_REF] the contrastive divergence (CD), managed to learn a non trivial dataset such as MNIST. [8] It was in the same period that RBMs became very popular in the ML community for its capability to pre-train deep neural networks (for instance deep auto-encoder), in a layer wise style. And, it was then showed that RBMs are universal approximator [9] of discrete distributions, that is, an arbitrary large RBM can approximate arbitrarily well any discrete distribution (which led to many rigorous results about the modelization mechanism of RBMs [START_REF] Montúfar | Restricted boltzmann machines: Introduction and review[END_REF] ). In addition, RBMs offer the possibility to be stacked to form a multi-layer generative model known as a deep Boltzmann machine (DBM). [START_REF] Salakhutdinov | Deep Boltzmann machines[END_REF] In the more recent years, RBMs continued to attract scientific interest. Firstly because it can be used on continuous or discrete variable very easily. [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF][START_REF] Yasuda | [END_REF][START_REF] Cho | Improved learning of Gaussian-Bernoulli restricted Boltzmann machines[END_REF][START_REF] Yamashita | To be Bernoulli or to be Gaussian, for a restricted Boltzmann machine[END_REF] Secondly, because the possible interpretations of the hidden nodes can be very useful. [START_REF] Hjelm | [END_REF][START_REF] Hu | Human brain mapping[END_REF] Interestingly, in some cases, more elaborate methods such as GAN [START_REF] Goodfellow | Generative adversarial nets[END_REF] are not working better. [START_REF] Yelmen | Creating artificial human genomes using generative models[END_REF] Finally it can be used for other tasks as well such as classification or representation learning. [20] Besides all these positive aspects, the learning process itself of the RBM remains poorly understood. The reasons are twofold: firstly, the gradient can be computed only in an approximated way as we will see; secondly, simple changes may have terrible impact on the learning or, messed up completely with the other meta-parameters. For instance making a naive change of variable in the MNIST dataset [START_REF] Cho | Enhanced gradient and adaptive learning rate for training restricted Boltzmann machines[END_REF][START_REF] Tang | Data normalization in the learning of restricted Boltzmann machines[END_REF] can affect importantly the training performance. Another example, when varying the number of hidden nodes, keeping the other mete-parameters fixed, will affect not only the representational power of the rbm but also the learning dynamics itself.

The statistical physics community, on its side, has a long tradition of studying inference and learning process with its own tools. Using idealized inference problems, it has managed in the past to shed light on the learning process of many ML models. For instance, in the Hopfield model, [START_REF] John | [END_REF][24][25][26] a retrieval phase was characterized where the maximum number of patterns that can be retrieved can be expressed as a function of the temperature. Another example is the computation of the storage capacity of the Perceptron [27] on synthetic datasets. [28,29] In these approaches, the formalism of statistical physics explains the macroscopic behavior of the model in term of its position on a phase diagram in the large size limit.

From a purely technical point of view, the RBM can be seen for a physicist as a disordered Ising model on a bipartite graph. Yet, the difference with respect to the usual models that are studied in statistical physics is that the phase diagram of a trained RBM involves a highly non-trivial coupling matrix where the components are correlated as a result of the learning process. These dependencies make it non-trivial to adapt classical tools from statistical mechanics, such as the replica theory. [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF] We will illustrate in this article how methods from statistical physics still have helped to characterize both the equilibrium phase of an idealized RBM where the coupling matrix has a structured spectrum, and how the learning dynamics can be analyzed in some specific regimes, both results being obtained with traditional mean-field approaches.

The paper is organized as follows. We will first give the definition of the RBM and review the typical learning algorithm used to train the model in Section 2. Then, in Section 3. , we will review different types of RBMs by changing the prior on its variables and show explicit links with other models. In Section 4. , we will review two approaches that characterize the phase diagram of the RBM and in particular its compositional phase, based on two different hypothesis over the structure of the parameters of the model. Finally, in Section 5. we will show some theoretical development helping to understand the formation of patterns inside the machine and how we can use the mean-field or TAP equations to learn the model.

Definition of the model and learning equations 2.1. Definition of the RBM

The RBM is an Ising model (or equivalently, a Markov random field), defined on a bipartite graph structure over two layers of variables: the visible nodes s i , for i = 1, . . . , N v and the hidden nodes τ a = 1, . . . , N h , denoting N v and N h the number of visible and hidden nodes respectively. In the following, we will use i, j, k, . . . to enumerate the visible variables and a, b, c, . . . for the hidden ones. No connection between any pair of visible or hidden nodes occurs . Hence, we will call 𝑤 the coupling or weight matrix and denote its elements as w ia since no other interactions are present (such as w i j or w ab ). In addition to the pairwise coupling matrix 𝑤, each visible and hidden node can have a local magnetic field, or local bias (we will refer to it as bias in the rest of the article), respectively named θ i and η a . We can introduce the following Hamiltonian

ℋ[𝑠, 𝜏 ] = -∑ ia s i w ia τ a -∑ i θ i s i -∑ a η a τ a , (1)
from which we define a Boltzmann distribution

p(𝑠, 𝜏 ) = 1 Z exp(-ℋ[𝑠, 𝜏 ]).
where Z is given by

Z = ∑ {𝑠},{𝜏 } exp(-ℋ[𝑠, 𝜏 ]).
The structure of the RBM is represented on Figure 1 where the visible nodes are represented by black dots, the hidden nodes by red dots and the weight matrix by blue dotted lines. The benefit of having a bipartite structure is that, when keeping fixed an entire layer, in our case all the visible or all the hidden nodes, the variables of the other layer become statistically independent. In other words, the measure p(𝑠|𝜏 ) and p(𝜏 |𝑠)

In MNIST, it is usual to consider binary variable {0, 1} to describe the dataset. Taking instead {±1} naively will affect dramatically the learning of the RBM.

xxxx-2

A c c e p t e d M a n u s c r i p t factorizes over the visible/hidden nodes respectively. This is an important property to keep in mind since it will be used in the learning procedure of the model. We will see that this property is widely used during the learning in order to draw new samples using a Monte-Carlo Markov Chain (MCMC) by Gibbs sampling. Historically, the RBM was first defined with binary {0, 1} variables for both the visible and the hidden nodes in line with the sigmoid activation function of the perceptron, hence being directly intepretable as spin-glass model of statistical mechanics. A more general definition is considered here by introducing a prior distribution function for both the visible and hidden variables, allowing us to consider discrete or continuous variables. This generalization will allow us to see the links between RBMs and other well-known models of machine learning. From now on we will write all the equations for the generic case using the notation q v (σ ) and q h (τ) to indicate an arbitrary choice of "prior" distribution. Averaging over the RBM measure corresponding to Hamiltonian (1) will then be denoted by

⟨ f (𝑠, 𝜏 )⟩ ℋ = ∑ {s,τ} p(𝑠, 𝜏 ) f (𝑠, 𝜏 ) (2)
where here Σ can represent both discrete sums or integrals and with the RBM distribution defined from now on as

p(𝑠, 𝜏 ) = 1 Z q v (𝑠)q h (𝜏 ) exp(-ℋ[𝑠, 𝜏 ]). (3)
It is worth mentioning that, the choice of the prior distribution can be rephrased in terms of an activation function on the conditioned distribution over the visible or hidden variables. Therefore, when specifying a prior distribution, we will systematically indicate the corresponding activation function for the hidden layer, that is p(𝜏 |𝑠), which is obtained using the Bayes theorem

p(𝜏 |𝑠) = p(𝑠, 𝜏 ) ∑ τ p(𝑠, 𝜏 ) = q h (𝜏 ) exp(-ℋ[𝑠, 𝜏 ]) ∑ {𝜏 } q h (𝜏 ) exp(-ℋ[𝑠, 𝜏 ])
Before entering more into the technical details about the RBM, it is important to recall that it has been designed as a "learnable" generative model in practice. In that sense, the usual procedure is to feed the RBM with a dataset, tune its parameter w, θ and η such that the equilibrium properties of the learned RBM reproduce faithfully the correlations (or the patterns) present in the dataset. In other words, it is expected that the learned model is able to produce new data statistically similar but distinct from the training set. To do so, the classical procedure is to proceed with a stochastic gradient ascent (to be explained in Section 2.2.) of the likelihood function that can be easily expressed. Usually the learning of ML models involves the minimization of a loss function which happens here to be minus the log likelihood, thus in the following we will refer to Stochastic Gradient Descent (SGD) instead. First, consider a set of datapoints {s

(d)
i }, where d = 1, . . . , M is the index of the data. The log-likelihood is given by

ℒ = 1 M M ∑ d=1 log ∑ {𝜏 } p(𝑠 (d) , 𝜏 ) = 1 M M ∑ d=1 log p(𝑠 (d) ) = 1 M M ∑ d=1 log ∑ τ q v (𝑠 (d) )q h (𝜏 ) exp -ℋ[𝑠 (d) , 𝜏 ] -log(Z) = 1 M M ∑ d=1 ∑ i θ i s (d) i + log q v (𝑠 (d) ) + ∑ a log ∑ τ a q h (τ a ) exp ∑ i s (d) i w ia τ a + η a τ a -log(Z) xxxx-3 A c c e p t
∂ ℒ ∂ w ia = 1 M M ∑ d=1 ∑ τ a q h (τ a )s (d) i τ a e ∑ i s (d) i w ia τ a +η a τ a ∑ τ a q h (τ a )e ∑ i s (d) i w ia τ a +η a τ a -⟨s i τ a ⟩ ℋ = 1 M M ∑ d=1 s (d) i ∑ τ a τ a p(τ a |𝑠 (d) ) -⟨s i τ a ⟩ ℋ = ⟨s i τ a ⟩ data -⟨s i τ a ⟩ ℋ (4)
where we used the following notation

⟨ f (𝑠, 𝜏 )⟩ data = 1 M M ∑ d=1 ∑ {𝜏 } f (𝑠 (d) , 𝜏 )p(𝜏 |𝑠 (d) ). (5)
The gradients for the biases (or magnetic fields) are

∂ ℒ ∂ θ i = ⟨s i ⟩ data -⟨s i ⟩ ℋ (6) ∂ ℒ ∂ η a = ⟨τ a ⟩ data -⟨η a ⟩ ℋ (7)
It is interesting to note that, in expression (4), the gradient is very similar to the one obtained in the traditional inverse Ising problem with the difference that in the inverse Ising problem the first term (sometimes coined "positive term") depends only on the data, while for the RBM, we have a dependence on the model (yet simple to compute). Once the gradient is computed, the parameters of the model are updated in the following way

w (t+1) ia = w (t) ia + γ ∂ ℒ ∂ w ia w (t) ia ,θ (t) i ,η (t) a (8)
θ (t+1) i = θ (t) i + γ ∂ ℒ ∂ θ i w (t) ia ,θ (t)
i ,η

(t) a (9) η (t+1) a = η (t) a + γ ∂ ℒ ∂ η a w (t) ia ,θ (t)
i ,η

(t) a (10) where γ called the learning rate tunes the speed at which the parameters are updated in a given direction, the superscript t being the index of iteration. A continuous limit of the learning process can be formally defined by considering t real and replacing t + 1 by t + dt, γ by γdt and letting dt → 0. The difficulty to train an RBM resides in the difficulty to compute the second term of the gradient, the so-called "negative term", which represents , in the gradient over the weight matrix, the correlation between a visible node i and a hidden node a under the RBM distribution. Similarly, the gradient over the biases is difficult to compute, where here the negative term is given by the mean value over the visible/hidden nodes. Depending on the value of the parameters of the model (the couplings and the biases), we can either be in a phase where it is easy to sample configurations from p(𝑠, 𝜏 ), (usually called paramagnetic phase); either be (if unlucky) in a spin glass phase, where it is exponentially hard to escape from the spurious free energy minima; either be (if lucky) in a "recall" phase where the dominant states correspond to data-like configurations. But even in the latter case, it might be difficult to transit from one state to another one with random jumps if these states are separated by large energy or free energy barriers, as in the Hopfield model for instance.

Stochastic Gradient Descent

Considering the difficulty to use the Eq. ( 4) to learn the model (the computation of the negative term scales exponentially with the system size, and Monte Caro Markov chains (MCMC) can be very slow to converge), an efficient approximative scheme name contrastive divergence [START_REF] Hinton | [END_REF] (CD) has been developed in order to approximate this term. First of all, the dataset is partitioned into small subsets called minibatches, and the gradient ascent is performed sequentially over all these minibatches in a random order. As a result each gradient step is performed only over a small subset of the whole dataset at a time. In order to estimate the negative term, the principle of CD is to start many Monte-Carlo chains in parallel, as many as the number of samples in a minibatch, and to use each sample of the minibatch as an initial condition for the chain. The idea being that starting from xxxx-4 A c c e p t e d M a n u s c r i p t desired equilibrium configurations and making k steps -the number of MC steps is coined in the method : CD-k-, we expect to explore nearby configurations representative of the dataset when the machine is learned; if otherwise the chains flow away they will "teach" the RBM how to adjust the parameters. The interpretation of CD is that it tends to create a basin of "attraction" centered on the datapoints where nearby configurations will be attractive to these datapoint under the Gibbs dynamics. In practice, starting from a datapoint 𝑠 d a random configuration of the hidden layer is sampled; in turn given this a configuration of the visible layer is sampled and so on for k steps. For this we take advantage of the bipartite structure of the model to draw a whole visible or hidden layer at once thanks to the factorization of the conditional distribution p(𝑠|𝜏 ) and p(𝜏 |𝑠):

𝑠 d → 𝜏 0 ∼ p(𝜏 |𝑠 d ) → 𝑠 1 ∼ p(𝑠|𝜏 0 ) → • • • → 𝑠 k ∼ p(𝑠|𝜏 k-1 ) → 𝜏 k ∼ p(𝜏 |𝑠 k ) (11)
finally 𝑠 k and 𝜏 k are used to estimate the negative term. It is clear that the CD-k is not directly minimizing the likelihood, or equivalently the Kullback Leibler (KL) divergence between the data distribution p 0 (𝑠) and the Boltzmann one p(𝑠). In reality it minimizes the KL divergence D KL (p 0 ||p k ) between the data distribution p 0 and the distribution obtained after k MC steps p k that is defined as

D KL (p 0 ||p k ) = ∑ {𝑠} p 0 (𝑠) log p k (𝑠) p 0 (𝑠) p k (𝑠 k ) = ∑ {s 0 ,...,s k-1 } ∑ {τ 0 ,...,τ k-1 } k ∏ l=1 p(𝑠 l |𝜏 l-1 )p(𝜏 l-1 |𝑠 l-1 ) p 0 (𝑠 0 )
In Ref. [31] it is argued that this procedure is roughly equivalent to minimizing the following KL difference

ℒ CDk = D KL (p 0 ||p) -D KL (p k ||p),
up to an extra term considered to be small without much theoretical guaranty. The major drawback of this method is that the phase space of the learned RBM is never explored since we limit ourselves to k MC steps around the data configurations, therefore it can lead to estimate very poorly the probability distribution for configurations that lie "far away" from the dataset. A simple modification has been proposed to deal with this issue in Ref. [32]. The new algorithm is called persistent-CD (pCD) and consists of having again a set of parallel MC chains, but instead of using the dataset as initial condition, they are first initialized from random initial conditions and then the state of the chains is saved from one update of the parameters to the next one. In other words, the chains are initialized one time at the beginning of the learning and are then constantly updated a few MC steps further at each update of the parameters. In that case, it is not longer needed to have as many chains as the number of samples in the mini-batch even though in order to keep the statistical error comparable between the positive and the negative term it should be of the same order. More details can be found in Ref. [32] about PCD and in Ref. [33] for a more general introduction to the learning behavior using MC. In Section 5. we will intend to understand some theoretical and numerical aspect of the RBMs learning process.

Overview of various RBM settings

Before investigating the learning behavior of RBMs, let us have a glimpse at various RBM settings and their relation to other models, by looking at common possible priors used for the visible and hidden nodes.

Gaussian-Gaussian RBM

The most elementary setting is the linear RBM, where both visible and hidden nodes have Gaussian priors:

q v (s i ) = 1 2πσ 2 v exp - s 2 i 2σ 2 v q h (τ a ) = 1 2πσ 2 h exp - τ 2 a 2σ 2 h
with σ v and σ h the intrinsic variance of the visible and hidden variables respectively. After summing over hidden variables we get a multi-variate Gaussian distribution over the visible ones. If not very sophisticated, the model is yet interesting because xxxx-5

A c c e p t e d M a n u s c r i p t it presents a non-trivial learning dynamics that can be written exactly. [START_REF] Karakida | Analyzing feature extraction by contrastive divergence learning in rbms[END_REF][START_REF] Karakida | [END_REF][36][37] When using Gaussian prior, the corresponding activation function p(𝜏 |𝑠) are Gaussian centered on σ 2 h ∑ i w ia s i :

p(𝜏 |𝑠) ∝ ∏ a exp - τ 2 a 2σ 2 h + τ a ∑ i w ia s i .
Let us write the marginal distribution over the visible nodes p(𝑠) (we omit the hidden bias since it can be canceled by a redefinition of the visible one), starting from Eq. ( 3) and integrating over the hidden variables we get

p(𝑠) = 1 Z ∏ i e - s 2 i 2σ 2 v +s i θ i ∏ a dτ a exp - τ 2 a 2σ 2 h + ∑ i s i w ia τ a = 1 Z ∏ i e - s 2 i 2σ 2 v +s i θ i ∏ a exp σ 2 h 2 ∑ i j s i w ia w ja s j = 1 Z exp -𝑠 T 1 2σ 2 v - σ 2 h 2 𝑤𝑤 T 𝑠 + 𝑠 T 𝜃 = 1 Z exp -𝑠 T 𝐴𝑠 + 𝑠 T 𝜃 ( 12
)
where we define the precision matrix

𝐴 ≡ 1 2σ 2 v - σ 2 h
2 𝑤𝑤 T . Now we can also identify the conditions for the existence of the measure p(𝑠). We need the matrix 𝐴 to be strictly positive definite, hence that the highest eigenvalue of 𝑤𝑤 T remains strictly below 1/(σ 2 v σ 2 h ). More interestingly, the Gaussian prior let us write in closed form the stochastic gradients (in fact we solve the deterministic equation, not the stochastic one), hence given us some hints on the nature of the learning dynamics of non-linear RBMs, since in any case we expect a linear regime to take place at the beginning of the learning process. In the present case, we can rewrite Eq. ( 4) as

∂ ℒ ∂ w ia = 1 M ∑ d s (d) i σ 2 h ∑ j s (d) j w ja -σ 2 h ⟨s i ∑ j s j ⟩w ja = σ 2 h ∑ j C i j w ja -∑ j ⟨s i s j ⟩w ja = σ 2 h ∑ j C i j w ja -∑ j A -1 i j w ja (13)
where

C i j = ⟨s i s j ⟩ data = M -1 ∑ d s (d) i s (d)
j is the correlation between the nodes i and j in the dataset, and 𝐴 -1 the inverse of the precision matrix. At this point, following, [36] it is convenient to use the singular value decomposition (SVD) of 𝑤. We note w ia = ∑ α u α i w α v α a the eigen-decomposition of the rectangular weight matrix 𝑤, where the matrix 𝑢 and 𝑣 correspond to the left (resp. right) eigenvectors of 𝑤 associated to the visible (resp. hidden) variables and w α the eigenvalue associated to the mode α. As can be seen in Eq. (12), this transformation will diagonalize the interaction term of the Hamiltonian of the system. We can now make the following change of variables

ŝα = ∑ i u α i s i and τα = ∑ a v α a τ a
under this change of variable, the Gaussian measure factorizes where ∑ i, j,a s i w ia w ja s j = ∑ α ŝα w 2 α ŝα and therefore

-𝑠 T 𝐴𝑠 = - 1 2 ∑ α ŝα 1 -σ 2 v σ 2 h w 2 α σ 2 v ŝα
Writing the distribution in this new basis we obtain

p( ŝ) ∝ ∏ α exp - ŝ2 α 2 1 -σ 2 v σ 2 h w 2 α σ 2 v
Hence, we can obtain an exact equation for the gradient in the basis of the SVD of the weight matrix 𝑤. First, we project the equation Eq. ( 13) on the modes αβ of the SVD of 𝑤 Now to simplify we discard the fluctuations associated to the stochastic gradient, by considering instead the full gradient and an infinitesimal learning rate such that we can consider the iteration time to be continuous and identify ∂ ℒ ∂ w ia ∼ dw ia dt . As a result we obtain the time derivative of the matrix 𝑤 decomposed over its eigenmodes

∂ ℒ ∂ 𝑤 αβ = ∑ ia u α i ∂ ℒ ∂ w ia v β a = ∑ ia u α i [⟨s i τ a ⟩ data -⟨s i τ a ⟩ ℋ ] v β a xxxx-6 A c c
d𝑤 dt αβ = ∑ ia u α i d dt ∑ γ u γ i w γ v γ a v β a = ∑ iaγ u α i u γ i dw γ dt v γ a v β a + u α i du γ i dt w γ v γ a v β a + u α i u γ i w γ dv γ a dt v β a = δ αβ dw α dt + (1 -δ αβ ) 𝑢 α d𝑢 β dt w α + w β d𝑣 α dt 𝑣 β
This equation shows that, the gradient update of 𝑤 can be decomposed when projected on the SVD basis of 𝑤 into a gradient over the mode w α and a rotation of the matrices 𝑢 α and 𝑣 α . Noticing first that ⟨ ŝα τα ⟩ = σ 2 h w α ⟨ ŝ2 α ⟩, we therefore end up with the following dynamics for the singular values w α :

dw α dt = d𝑤 dt αα = ⟨ ŝα τα ⟩ data -⟨ ŝα τα ⟩ ℋ = σ 2 h w α ⟨ ŝ2 α ⟩ data -⟨ ŝ2 α ⟩ ℋ = σ 2 h w α ⟨ ŝ2 α ⟩ data - σ 2 v 1 -σ 2 v σ 2 h w 2 α ( 14
)
where in Eq. ( 14) ⟨ ŝ2 α ⟩ data denotes the variance of the components of the data on the mode α:

⟨ ŝ2 α ⟩ data = ∑ ia u α i 1 M ∑ i j s (d) i s (d) j u α j
This first result tells us that when keeping the matrices 𝑢 and 𝑣 fixed, the SGD on the mode w α will adjust the value of w α such that the r.h.s matches the variance in the direction given by 𝑢 α , giving the following limit values:

w 2 α = ⟨ ŝ2 α ⟩ data -σ 2 v σ 2 v σ 2 h ⟨ ŝ2 α ⟩ data if ⟨ ŝ2 α ⟩ data > σ 2 v 0 if ⟨ ŝ2 α ⟩ data < σ 2 v ( 15
)
We remark that, if the empirical variance given by the data is smaller than the prior variance of the visible variables the corresponding mode is filtered out. The evolution of the matrices 𝑢 α and 𝑣 α can also be obtained [37] from the following expression in the present case

Ω u αβ ≡ d𝑢 α dt T 𝑢 β = -(1 -δ αβ )σ 2 h w β -w α w α + w β - w β + w α w α -w β ⟨s α s β ⟩ data ( 16
)

Ω v αβ ≡ d𝑣 α dt T 𝑣 β = -(1 -δ αβ )σ 2 h w β -w α w α + w β + w β + w α w α -w β ⟨s α s β ⟩ data ( 17
)
of the infinitesimal rotations of the vectors 𝑢 α and 𝑣 α . In the particular case of the Gaussian-Gaussian RBM, we can note the absence of term averaged over the model ⟨.⟩ ℋ . This is due to the fact that the SVD corresponds to the eigendecomposition of the RBM measure (that is, the Gaussian measure factorizes over the singular modes) and that the Eqs. (16-17) involve correlation between modes α ̸ = β which are zero here. From Eq. (16-17), we see that a steady state is found when a direction 𝑢 α is found that diagonalizes the empirical covariance matrix of the dataset.

In short, the Gaussian-Gaussian RBM learns the principal components of the dataset and for each principal axes the weight matrix is adjusted until the strength of the corresponding modes w α reaches the value given by Eq. ( 15). Of course, modes above threshold acquire a variance which matches the variance of the dataset in this direction ⟨s 2 α ⟩ ℋ = ⟨s 2 α ⟩ data . We can somehow say that the Gaussian-Gaussian RBM is performing a sort of SVD of the dataset, keeping only the modes above a given threshold. It Actually these equations are given with a wrong sign in Refs. [37] which is corrected here.

xxxx-7

A c c e p t e d M a n u s c r i p t is worth noting that an analysis has been done in Ref. [35] where it is shown that updating the parameters of the model using the kCD approximation converges toward the same solution as the one obtained by maximizing the likelihood of the model.

We can illustrate the learning mechanism in simple cases where it is possible to solve explicitly the dynamics. First assume that the RBM has found the principal axes, i.e. consider the matrices 𝑢 and 𝑣 to be fixed. In this case the quantity ⟨ ŝ2 α ⟩ data remains constant. Letting

x α = σ 2 v σ 2 h w 2 α and δ α = ⟨ ŝ2 α ⟩ Data -σ 2 v σ 2 v ,
and rescaling time as tσ 2 v σ 2 h → t, equation ( 14) then rewrites as

ẋα = 2x α δ α - x α 1 -x α
and we obtain a solution of the form

x α (t) = f -1 α (δ α t), with f α (x) = log x x α (0) - 1 1 + δ α log γ α -x γ α -x α (0)
, and

γ α = δ α 1 + δ α .
For δ α ≪ 1 we get a sigmoid type behaviour

x α (t) x α (0) = δ α e δ α t δ α + x α (0) e δ α t -1)
.

To illustrate the rotation of the modes, consider now the situation where there are 2 modes u α , α = 1, 2 which are a linear combination of two dominant modes of the data { û1 , û2 } with identical orientation taken in this order, all other modes considered to be already properly aligned with the data. Let then θ represent the angle between u 1 and û1 (and also between u 2 and û2 see Figure 2. Equation ( 16) for this pair of modes rewrites then as

dθ dt = -σ 2 h w 2 α + w 2 β w 2 α -w 2 β ⟨s 1 s 2 ⟩ Data (t), with ⟨s 1 s 2 ⟩ Data (t) = cos θ sin θ ⟨s 2 2 ⟩ Data -⟨s 2 1 ⟩ Data , ⟨s 2 1 ⟩ Data (t) = cos 2 θ ⟨s 2 1 ⟩ Data + sin 2 θ ⟨s 2 2 ⟩ Data ⟨s 2 2 ⟩ Data (t) = sin 2 θ ⟨s 2 1 ⟩ Data + cos 2 θ ⟨s 2 2 ⟩ Data
so that finally we get a dynamical system of the form

ẋ1 = 2x 1 δ 1 cos 2 θ + δ 2 sin 2 θ - x 1 1 -x 1 (18) ẋ2 = 2x 2 δ 1 sin 2 θ + δ 2 cos 2 θ - x 2 1 -x 2 (19) θ = - 1 2 (δ 1 -δ 2 ) x 1 + x 2 x 1 -x 2 sin(2θ ) (20)
Note that at fixed x 1 and x 2 the dynamics of θ corresponds to the motion of a pendulum w.r.t the variable θ ′ = 4θ shown on Figure 2.

xxxx-8

A c c e p t e d M a n u s c r i p t Fig. 2. Angle between the reference basis given by the data and the moving one given by the RBM shown on the up left panel. Equivalence with the motion of a pendulum is indicated on the left bottom panel. Solution of (18,19,20) of two coupled modes in the linear RBM (right panel).

Gaussian-Spherical

The Gaussian-Gaussian case is interesting as a solvable model of RBM but of limited scope, since p(𝑠) reduces in the end to a multivariate Gaussian. Next, a simple non-linear RBM which remains exactly solvable is based on the so-called spherical model. [38,39] For this model, it is possible to compute the phase diagram and the equilibrium states once the coupling matrix is given -more precisely, when the spectral density of the coupling matrix is given. Here we chose the following priors to impose a spherical constraint on the hidden nodes:

q v (s i ) = 1 2πσ 2 v exp - s 2 i 2σ 2 v q h (𝜏 ) = δ ∑ a τ 2 a -σ N h N v
where σ is a parameter of the model. [40] The interest of such an RBM is first that the spherical constraint can be dealt with analytically. [40,41] Secondly the model can exhibit a phase transition unlike the Gaussian-Gaussian case. Absorbing the parameter σ 2 v in the definition of the weight matrix, to follow the computation of Ref. [40], a simple analysis in the thermodynamic limit tells us that the phase transition takes place when the parameter σ exceeds the value σ c , where σ c depends on the value of the highest mode w max and of the form of the spectrum of 𝑤 (typically, σ 2 c ∝ 1/w 2 max , where the pre-factor depends on the form of the spectrum). The condensation along this mode of the visible (resp. hidden) magnetization is then given by

m α = 1 √ L ∑ i u α i ⟨s i ⟩ ℋ = w max σ σ 2 -σ 2 c m2 α = 1 √ L ∑ a v α a ⟨τ a ⟩ ℋ = σ 2 -σ 2 c
where we defined L = √ N v N h . This type of RBM is again of limited scope to represent data. In the thermodynamic limit a finite number n = 𝒪(1) of modes can condense. They necessarily accumulate at the top of the spectrum of the weight matrix and represent a distribution concentrated on an n-dimensional sphere in absence of external fields while other non-condensed modes are responsible for transverse Gaussian fluctuations. The dynamical aspect of this model will be discussed in Section 5. .

To end up this section let us also mention that the finite size regime is amenable to an exact analysis when restricting the weight matrix spectrum to have the property of being doubly degenerated (see Ref. [40] for details).

Gaussian-Softmax

The case of the Gaussian mixture if rarely viewed like that, fits actually perfectly the RBM architecture. Consider here the case of Gaussian visible nodes and a set of discrete {0, 1} hidden variables with a constraint corresponding to the softmax activation function: [42]

q v (s i ) = 1 2πσ 2 v exp - s 2 i 2σ 2 v , xxxx-9
A c c e p t e d M a n u s c r i p t

Chin. Phys. B Vol. 30, No. x (2021) xxxx

q h (𝜏 ) = ∏ a (δ τ a ,0 + δ τ a ,1 ) δ ∑ b τ b ,1 .
With this formulation, we indeed see that the conditional probability of activating a hidden node is a softmax function

p(τ a = 1|𝑠) = exp (∑ i w ia s i + η a ) ∑ b exp (∑ i w ib s i + η b )
.

It is easy from this expression to recognize the equations of the Gaussian mixture model (GMM), [START_REF] David | Information theory, inference and learning algorithms[END_REF][START_REF] Christopher | Pattern recognition and machine learning[END_REF] where the latent variable τ a indicates if a sample belong or not to the center a. The position of the associated center is given by the vector 𝑤 a . It is even clearer when writing the marginal over the visible nodes after summing over the hidden nodes in Eq. ( 3)

p(𝑠) = 1 Z ∑ a exp η a + ∑ i - s 2 i 2σ 2 v + θ i s i + s i w ia = 1 Z ∑ a exp η a + ∑ i - 1 2σ 2 v s i -σ 2 v [w ia + θ i ] 2 2 + 1 2 σ 2 v [w ia + θ i ] = 1 Z ′ ∑ a ρ a exp ∑ i - 1 2σ 2 v s i -σ 2 v [w ia + θ i ] 2 2
by identifying

ρ a ≡ exp η a + ∑ i σ 2 v w ia + σ 2 v θ i 2 ∑ b exp η b + ∑ i (σ 2 v w ib + σ 2 v θ i ) 2 (21)
the weight of the mode a in the Gaussian mixture centered in 𝑤 a . Now we can see that the extra parameter θ i can be absorbed in the definition of the weight matrix w ′ ia = w ia + θ i . It turns out that the positive term of the gradient in equation ( 4) (ignoring ρ a ) corresponds to the gradient that is obtained in the GMM. This can be reformulated into the Expectation Maximization (EM) update by considering that p(τ a |𝑠) do not depend on w ia , hence doing the "expectation" step:

⟨s i τ a ⟩ data = 1 M ∑ d s (d) i -σ 2 v w ia p(τ a |𝑠 (d) ) (22)
If we impose that the gradient is zero, doing now the "maximization" step, we obtain

w (t+1) ia = ∑ d s (d) i p(τ a |𝑠 (d) ) σ 2 v ∑ d p(τ a |𝑠 (d) ) (23)
where the l.h.s. is to be understood as the new values for the parameters w ia while the conditional distribution on the r.h.s. depends on w (t) ia . For an RBM, one would also compute the negative term of the gradient, involving the derivative of ρ a w.r.t. w ia . We obtain the negative term

⟨s i τ a ⟩ ℋ = 1 M ∑ d σ 2 v w ia p(τ a |𝑠 (d) ) -ρ a (24)
Again, we can recover with equation ( 24) the EM update for the density of the Gaussian mode a in the GMM, by first considering that the conditioned distribution p(τ a |𝑠 (d) ) does not depend on w ia (expectation step) and by putting the l.h.s. to zero (maximization step). The fact that when using the RBM formalism we do not obtain directly the same EM equations as in the GMM is due to the different parametrization of the parameters. In the GMM, the density of each Gaussian is defined right from the beginning as an independent parameter while when using the RBM, the density of the Gaussian depends on other parameters such as the weight matrix 𝑤.

Phase transition in the learning processan interesting phenomena occurs in this model when learning position of the centers of the Gaussian while submitting the variances σ v of the Gaussian to an annealing process. [START_REF] Rose | [END_REF] First of all, starting from a very high variance (equivalently, very high temperature), we can convince ourselves that the learning will end up finding the center of mass of the dataset. Let us therefore consider that we centered the dataset beforehand: ∑ d s 23) around this point with η a = 0 and w ia ≈ 0 + ε ia , where the 𝜖 are small perturbations, we can derive the threshold where the linear perturbations get amplified. The linear stability analysis leads to the following equations for the perturbation ε:

ε (t+1) ia ≈ ∑ d s (d) i (1 + ∑ j s (d) j ε (t) ja -1 N h ∑ jb s (d) j ε (t) jb ) σ 2 v ∑ d (1 + ∑ j s (d) j ε (t) ja -1 N h ∑ jb s (d) j ε (t) jb ) = 1 σ 2 v ∑ j c i j ε (t) ja - 1 N h ∑ b ε (t) jb
where c i j is the covariance matrix of the dataset. From this expression, one sees that when the variance is higher that the largest eigenvalue Λ C of 𝑐, i.e. σ 2 v > Λ C , the solution w ia = 0 is stable. Then, when σ 2 v < Λ C , the solution is unstable and the system starts to learn something more about the dataset besides its center of mass. It is interesting to note that this threshold is very similar to the one obtained in equation ( 15) for the Gaussian-Gaussian RBM. In this model, it is then possible to study the cascade of phase transition, occurring in a hierarchical way on structured datasets. [46,47] We stress here that, even if it is possible to project the learning equations on the SVD of the weight matrix as in the two previous analysis, it does not provide much more insight since this case cannot be solved exactly by this transformation.

It is also interesting to investigate the behavior of the exact gradient (not using EM) in the presence of a learning rate γ. When using the gradient, the update equations are given by w

(t+1) ia = w (t)
ia + γ∆w ia . In that case we obtain the following equation for the linear stability

ε (t+1) ia = (1 -γ)ε (t) ia + γ σ 2 v ∑ j c i j ε (t) ja - 1 N h ∑ b ε (t) jb .
Interestingly, the threshold does not depend on the value of γ in that case, meaning that the instability is a generic properties of the learning dynamics. The only change is the speed with which the instabilities will develop.

Bernoulli-Gaussian RBM

The next case is the Bernoulli-Gaussian RBM where we consider the following prior

q v (s i ) = 1 2 (δ s i ,0 + δ s i ,1 ) , q h (τ a ) = 1 2πσ 2 h exp - τ 2 a 2σ 2 h .
Again, a Gaussian prior implies that the activation function is Gaussian. It is interesting to consider this version of the RBM through its relation with the Hopfield model [START_REF] John | [END_REF] realized in Ref. [48]. Since the hidden variables are Gaussian they can be integrated out which leads to a simple analytical form for the marginals of the visible variables. In some recent works, the opposite approach has been done, starting with an Hopfield model and expressing it as a RBM using the Hubbard-Stratonovitch (HS) transformation (expressing the exponential of a square as Gaussian integral) to decouple the interactions between each spin. [49,50] After integrating over the hidden nodes in Eq. ( 3), we end up with the following distribution

p(𝑠) = 1 Z exp σ 2 h 2 ∑ i j s i s j ∑ a w ia w ja
We recognize a Hopfield model where the patterns are given by the weights w ia of the RBM and the effective coupling between two variables i and j is J i j = ∑ a w ia w ja . We can also consider that the variances of the hidden nodes is related to the temperature of the model. Some experiments have been conducted in Ref. [51] in order to compare the learning process of the Hopfield model versus the Bernoulli-Gaussian RBM on artificial data generated from an Hopfield model with dicrete patterns. It is interesting to note that, when assuming discrete patterns, the inverse procedure can be formulated in terms of an approximated Hopfield model. Thus, the inference of the pattern can be done directly using a set of TAP equations of the Hopfield model, and it has been shown that the artificial patterns were inferred exacltly. When using the RBM's formulation, in the absence of information over the patterns, only the subspace covered by the patterns were retrieved with a weak overlap with the true patterns. In fact, in that case the marginal over the visible nodes is a function of 𝑤𝑤 T , which is invariant by rotation of the 𝑣 matrix. It explains why the learned weight matrix in the RBM context does not overlap with the true patterns. With this machine, it is also possible to impose a maximum rank in order to reduce the number of parameters needed to describe the dataset giving the possibility of a trade-off between a good description of the dataset and the number of parameters. This properties has been used in Ref. [50] to find global patterns in protein foldings, using the RBM version of the Hopfield model with q discrete states.

Gaussian-Bernoulli RBM

At this point we now focus on models where the hidden layers will have a stronger impact. The integration of the hidden layer will not end up in a simple analytical form and therefore will make it difficult to understand the effect of the features and to characterize properly the learning dynamics. We first mention the Gaussian-Bernoulli case dealing with the following priors:

q v (s i ) = 1 2πσ 2 v exp - s 2 i 2σ 2 v , (25)
q h (τ a ) = 1 2 (δ τ a ,0 + δ τ a ,1 ) . ( 26
)
When using the discrete {0, 1} variables, we obtain the sigmoid activation function for the hidden nodes

p(τ a = 1|𝑠) = 1 1 + exp(-∑ i w ia s i + η a )
.

With this parameterization, it is natural to interpret a hidden node τ as an active feature when τ = 1 and an inactive one if τ = 0. When responding to a given input through the conditional probability p(𝜏 |𝑠), the machine is turning on the hidden nodes corresponding to overlapping features with the input. Therefore, the input undergoes a non-linear decomposition on the learned features. Saying it that way, it is somewhat reminiscent of the Independent Component Analysis [START_REF] Hyvärinen | [END_REF] where a matrix X is factorized on a set of independent sources or components 𝑦: 𝑥 = 𝐴𝑦. The sources here are independent in the sense that they are independently distributed. In the context of ICA, the goal is to find the inverse of the mixing matrix in order to recover the sources from the received signal. Concerning this particular RBM, it is proven [START_REF] Karakida | [END_REF] that under some assumptions -(i) having the same number of visible and hidden nodes, (ii) that the signal comes from a set of independent sources and (iii) that the variance of the visible variables is much smaller than the mean of the signal-there exists a stable solution for the learning dynamics where the learned weight matrix corresponds to the un-mixing matrix of the signal. In this regime, the RBM acts as an ICA. In other words, if the signal 𝑠 d used as an input for the RBM can be written as a mixture of sources: 𝑠 = 𝐴𝑦, a stable solution of the learning process consists in recovering the inverse mixing matrix in the weight matrix: 𝑤 = 𝐴 -1 .

To end up with this variant of the RBM, it is interesting to note that the prior variance of the visible variables here is in principle a fixed parameter. It has been noted that when using the prior (25-26) the mean of the conditional distribution over the visible p(s i |𝜏 ) is stretched by the variance σ v . It might be useful to remove this effect by renormalizing the weight matrix and the visible biases as in Ref. [14]: 𝑤 → 𝑤/σ 2 v and θ i → θ i /σ 2 v . Using this parametrization, we obtain

p(s i |𝜏 ) ≈ 𝒩 (θ i + ∑ a w ia τ a , σ 2 v ),
where 𝒩 represents the normal distribution. Note that it is possible to include the learning of these parameters in the likelihood ascent as in Ref. [14]. It is however important to stress here that even if appealing, the possibility to tune the variance of each visible node doesn't solve the problem of learning individual variances of separated clusters in a dataset. Indeed, consider the problem where the dataset is formed of many well-separated clusters with distinct variances. For a given visible node i, its variance computed over the whole dataset or instead over a given cluster have no reason to coincide. And the the prior variance if properly learned will only account for the global variance of this node. This should involve a more complex setting of the RBM which we won't discuss here in order to account for individual variances of clusters in a complex dataset.

Bernoulli-Bernoulli RBM

The last model here is traditionally the one which is implied when speaking of RBM. In that case both the visible and hidden nodes are in {0, 1} with the following priors The activation functions are sigmoid functions, for both the hidden and visible nodes

q v (s i ) = 1 2 (δ s i ,0 + δ s i ,1 ) , q h (τ a ) = 1 2 (δ τ a ,0 + δ τ a ,1
p(s i = 1|𝜏 ) = 1 1 + exp(-∑ a w ia τ a + θ i ) , (27)
p(τ a = 1|𝑠) = 1 1 + exp(-∑ i w ia s i + η a ) . (28)
In that case, the prior distribution has the advantage of not having any free parameter to be determined. In practice this model is used when dealing with a discrete dataset while the Gaussian-Benoulli is for continuous ones.

Rectified Linear Units (RELU) -Let us briefly mention how the Bernoulli prior on the hidden nodes can be linked to the RELU activation function [53] for the RBM. In a work by Teh et al., [START_REF] Teh | Rate-coded restricted Boltzmann machines for face recognition[END_REF] was highlighted one important shortcoming with Bernoulli prior. With the hidden variable in {0, 1}, a given pattern can be expressed τ = 1, or not τ = 0. Therefore the influence of a feature is binary, either 0, either a fixed amount given by the value of 𝑤: it is not possible to tune this amount as a function of how strongly a hidden node responds to a visible configuration. Of course it is possible for the machine to learn many times the same pattern, but this doesn't seem very efficient. A simple idea to correct this problem is to duplicate many time a hidden node, keeping the same features and bias values. Then, if the probability of turning on this hidden node is p, the average number of activated hidden nodes for this feature will be N p giving the possibility to tune the intensity of the feature.

Generalizing this idea, it is possible to construct an infinite number of replica, [START_REF] Nair | Rectified linear units improve restricted Boltzmann machines[END_REF] adjusting the bias for each of them such that in order to activate more and more neurons it is necessary that the signal ∑ i w i s i is stronger and stronger. Let us focus for a moment on a single hidden nodes with a feature w i and a bias η along with its replicas a ′ = 1, . . . , N r . We denote r = ∑ i w i s i + η the potential associated for this neuron given the signal 𝑠. The number of activated replica will be given by

1 √ N r N r -1 ∑ a ′ =0 sig r(1 -a ′ / √ N r ) ≈ N r →∞ log(1 + exp(r)) (29)
where we defined the sigmoid function sig(x) = (1 + exp(-x)) -1 . The r.h.s. of Eq. ( 29) is very close to the RELU activation function RELU(x) = max(0, x), hence showing that having all these replicas gives a similar activation function as RELU. In practice, it is not very efficient to have a large number of sigmoids for the training algorithm. An approximation is found by using the truncated Gaussian distribution. The average number of activated replica is then given by τ a = max(0, r + 𝒩 (0, σ a ))

where now τ a is a RELU hidden node and σ a is the variance associated with the number of activated replicas for the hidden node a. Eq. ( 30) can now be seen as an approximation of the Truncated-Gaussian prior for the hidden nodes

q h (τ a ) ∝ δ τ a >0 exp(- τ 2 a 2σ h ) (31)
In the following section, we will focus mainly on the Bernoulli-Bernoulli setting, its equilibrium phase diagram and its learning dynamics in the mean-field regime.

Phase diagram of the Bernoulli-Bernoulli RBM

In this section, we discuss various aspects of the phase diagram of the Bernoulli-Bernoulli RBM. In the rest of the section we will use {±1} instead of the usual {0, 1} for commodity. There are (at least) two series of works dealing with the RBM in the thermodynamic limit, each of them making different hypothesis on the statistical ensemble from which the RBM is taken. In the first one [START_REF] Barra | [END_REF]57] the weight matrix is taken from a simple statistical ensemble with iid elements and possibly additional sparse constraints on the patterns as will be explained in Section 4.1. . In the second one [36,37] it is assumes that the weight matrix contains a structured part of rank K = 𝒪(1) in addition to a random matrix corresponding to noise; the main results of this approach will be exposed in Section 4.2. . Both approaches are based on the replica computation [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF] of the free energy. For systems with quenched disorder, this is a classical approach (the replicas or its equivalent formulation) to find the macroscopic behavior. [START_REF] Barra | [END_REF][58][START_REF] Tubiana | Restricted Boltzmann machines: from compositional representations to protein sequence analysis[END_REF][START_REF] Agliari | [END_REF][61] 4.1. Mean-field approach, the random-RBM This MF approach to the macroscopic behavior of the RBM is based on statistical ensembles with iid elements of the weight matrix. Here, a random ensemble for the weight matrix is defined as follows. The weight matrix will be constructed using binary

w ia = ξ ia √ N v
. Now, each pattern is selected to be

ξ ia =    0 p r ∼ 1 -p i +1 p r ∼ p i /2 -1 p r ∼ p i /2 (32)
Using this definition, the degree of sparsity of the system is p = ∑ i p i /N v . The term random-RBM was coined by Tubiana et al. [57] but Agliari et al. [62,63] worked on a similar model although with a different theoretical approach. In particular, they computed the phase diagram in Ref. [56]. We start by reproducing here the argument of Agliari that was developed for the RBM with a finite number of patterns before switching to the replica computation done in Tubiana's thesis. [START_REF] Tubiana | Restricted Boltzmann machines: from compositional representations to protein sequence analysis[END_REF] Parallel retrievingthe usual definition of the Hopfield model (which we recall here is analogous to a Binary-Gauss RBM, see Section 3.), consists in using extensive pattern ξ a i = ±1/ √ N v for all i = 1, . . . , N v , where a = 1, . . . , P, P being the number of patterns. The Hopfield model in the low storage regime, where the number of pattern is fixed, or scales logarithmically with the system size, is characterized by a low temperature regime made of configurations with an extensive overlap with one of the patterns. This model can be recovered from a binary-binary RBM where the number of hidden nodes has the same scaling. Hence, having N h ∼ log(N v ), we can write exactly the partition function of the binary-binary RBM in the limit of large system size

Z = ∑ {𝑠},{𝜏 } exp β ∑ i,a s i w ia τ a = ∑ {𝑠} ∏ a cosh β √ N v ∑ i w ia s i ≈ ∑ {𝑠} exp β 2 2N v ∑ i, j ∑ a w ia w ja s i s j = ∑ {𝑠} exp N v β 2 2 ∑ a m a (𝑠) 2
recovering the Hopfield model with a square inverse temperature, and where we define the magnetization along the pattern a as m a (𝑠). In Ref. [62], the authors considered a weight dilution as in Eq. ( 32) applied to the above binary-binary RBM, or equivalently to a Hopfield model with a rescaled temperature. It is important to mention that it is a different procedure than diluting the network itself, see Ref. [64] for more details on the other case. Having sparse patterns allows the network to retrieve more than one pattern at a time. In particular, global minima of the free energy can have an overlap with many patterns and locally stable states can be composed of a complex mixture of patterns. We reproduce below on Figure 3 the plot from the article [62] showing the overlap over three and six patterns in the (almost) zero temperature limit. We observe on the left panel that one pattern is fully retrieved when the dilution is low. Then, when increasing p i , more and more patterns are retrieved together until the system enters a paramagnetic phase at high dilution. Left: a case with 3 patterns where we can observe how at small dilution, only one pattern is fully retrieved while the second and third one appear for larger dilution. Right: a case with 6 patterns where the figure is zoomed in the high dilution region where the branching phenomena is occurring and all the overlaps converge toward the same value.

Replica approach of the random-RBM -We will now follow the approach of Tubiana et al. [START_REF] Tubiana | Restricted Boltzmann machines: from compositional representations to protein sequence analysis[END_REF] and give more details on the derivation. This approach is based on a Bernoulli-RELU architecture giving the possibility to have continuous positive value for the hidden variables. The characterization of the phase diagram is based on the determination of the free energy in thermodynamic limits. Given the weight ensemble (see Eq. ( 32)), the weight matrix is now made of independent and sparse elements. In this context, the replica analysis can be used to perform the quenched average. The replicated interaction term can be first easily computed and gives

E 𝑤 exp β ∑ p s p i w ia τ p a ≈ exp p i β 2 2N ∑ pq s p i s q i τ p a τ q a
for the interaction term (ia). The interaction between the visible and hidden nodes can be decoupled using the HS transformation exp p 2 i β 2 2N ∑ pq s p i s q i τ p a τ q a ∼ ∏ pq dq pq d qpq 2π exp -Nβ (q pq qpqq pq p i p s p i s q i -

1 2 β pτ p a τ q a )
introducing the spin-glass order parameter over the replicas (we denoted by p, q, . . . , the replica index):

qpq ∼E w [⟨τ q a τ p a ⟩] q pq ∼E w [⟨(p i /p)s q i s p i ⟩]
where we note that the parameters over the visible nodes are weighted by the sparsity of the network. Using the replica symmetric ansatz, the quenched free energy can be computed and from it a set of order parameters emerges. An new order parameter is introduced in their derivation: the number L of hidden nodes that have a macroscopic activation ∼ 𝒪(m √ N), while the other ones remain silents (of order 1). This parameter is reminiscent from the replica approach of the Hopfield model. [24][25][26] In this approach, the number of pattern that can be expressed is fixed, in order to investigate the stability of retrieving one or more patterns. The important difference here is that the sparsity p imposes that the diluted weights can let many patterns to be expressed at the same time. Hence, the phase diagram will be characterized by, the value of the weight sparsity p and L the number of activated hidden nodes. The phase diagram is computed numerically, by scanning the possible value for the order parameters (see Ref. [59] for more details). It is found that when • p = 1 and L = 1: no sparsity and only one hidden node is activated. At low temperature, it gives back the behavior of the well-known Hopfield model having a recall phase of the patterns. An interesting additional result when using ReLu activations is that the capacity of the network can be increased by playing with the bias on the hidden nodes, at the cost of reducing the basin of attraction of the patterns.

• p < 1 a ferromagnetic transition is found when imposing L = 1, where one pattern is recalled at a time.

• p < 1, when all the hidden nodes are all weakly activated, a SG phase is found.

• p < 1 and L is such that 1 ≪ L ≪ N h ; a compositional phase is numerically identified. It is characterized by an intermediate number of hidden nodes strongly activated.

In this analysis, it is demonstrated that in the possible equilibrium behaviors of the random-RBM, an interesting phase mixing many patterns is present that characterizes in some way the efficient working regime of a learned RBM. It is of course a simplified case where the patterns are {±1} with a certain dilution factor. Now, the fact that there exists a family of weights where this phase exists is quite different from showing that the learning dynamics converges toward such a phase and how. In Tubiana's thesis, a stability analysis of the different phases is done showing that for a range of parameters of the RBM, the compositional phase is indeed the dominant one. Then, a certain number of numerical results are provided on the MNIST dataset which tends to confirm that the behavior of the learned RBM looks similar to a "compositional phase". It would therefore be of great interest to characterize the learning curve theoretically in order to understand how this phase is reached. It is also interesting to mention a recent work investigating the role of the diluted weights [65] during the learning in a RBM with one hidden node. In this article, it is shown that the proportion of diluted weights tends to vanish during the learning procedure. This might be a signal that when the number of hidden features is very low, the RBM automatically adjusts itself in the ferromagnetic phase described above, learning a global pattern of the dataset. The difficulty with the RBM is to be able to study the phase diagram of the model without discarding the fact that during the learning, the weights w ia become correlated between each others: starting from independently distributed w ia , we can observe how the spectrum of the weight matrix is modified during the learning (see Figure 11 for instance). Classical approaches in statistical mechanics consider a set of independent weights, all identically distributed, before trying to compute the quenched free energy of the system by the replica trick. In the present case the hypothesis of independent weights cannot hold, as can be seen by looking at the spectrum of the weight matrix at the beginning of the learning and a few iterations later. The absorption of information by the machine prompts the development of strong correlations. This phenomena is illustrated in Section 5.3. on Figure 11. In order to understand how these eigenvalues affect the phase diagram of the system, it is reasonable to assume a particular statistical ensemble of the weight matrix of the form

w ia = K ∑ α=1 u α i w α v α a + r ia (33)
where K ≪ N v , assuming a low-rank decomposition of the weight matrix plus some random noise r ia , where 𝑟 is a random matrix with iid centered Gaussian elements with variance σ . With this decomposition we assume that the eigenvalues w α correspond to some intrinsic property of a learned dataset, while the matrix, 𝑟, 𝑢 and 𝑣 can be treated as quenched disorder and averaged over. The set of vectors 𝑢 α and 𝑣 α correspond approximately to the left and right eigenvectors of the matrix 𝑤. We can thus start to average Z n over all these variables. Starting with the average over the random matrix 𝑟, it introduces the following interaction term ∑ ia,p̸ =q s p i s q i τ p a τ q a , where p, q runs over the n replicas. In this term, it is possible to decouple the interaction between the visible and the hidden nodes by introducing the overlap parameters

Q pq ∼E r,v,u [⟨τ q a τ p a ⟩] (34) Qpq ∼E r,v,u [⟨s q i s p i ⟩] (35)
Then, the form of the weight matrix, Eq. ( 33) leads to the following change of variable

s α = 1 √ L ∑ i s i u α i τ α = 1 √ L ∑ a τ a v α a
where L = √ N v N h . It corresponds to the projection of the visible and hidden variables over the matrix 𝑢 and 𝑣 coming from the SVD of 𝑤. With this projection we will be able to define the order parameters of the system as the condensation of the visible and hidden nodes over the SVD modes of 𝑤. Using again the Hubbard-Stratanovitch (HS) transformation in order to define the replicated magnetization

exp ∑ ia s p i w ia τ p a = exp ∑ α w α s p α τ p α ∝ ∏ α dm p α d mp α 2π exp -L ∑ α w α (m p α mp α -m p α s p α -mp α τ p α )
we obtain two additional order parameters

m p α ∼E r,v,u [⟨τ p a ⟩] mp α ∼E r,v,u [⟨s p i ⟩]
namely the condensation of the visible (resp. hidden) nodes along the SVD modes of 𝑤. After some computation, the replicated free energy is obtained

E u,v,r [Z n ] = ∏ p,α dm p α d mp α 2π ∏ p̸ =q dQ pq d Qpq 2π × exp -L ∑ p,α w α m p α mp α + σ 2 2 ∑ p̸ =q Q pq Qpq - 1 √ κ A[m, Q] - √ κB[ m, Q] ( 36
)
in few words, considering the quantity Z n for a given (integer) n, where Z is the partition function, for small n, we can develop Z n ≈ 1 + n log(Z). The key point here is that it is generally possible to compute the quenched Z n and then making a small n expansion.

xxxx-16

A c c e p t e d M a n u s c r i p t

where E indicates an average over the variables that are in subscripts and κ = N h /N v . The quantities A and B are given by

A[m, Q] ≡ log ∑ S a ∈{-1,1} E u e √ κσ 2 2 ∑ p̸ =q Q pq S p S q +κ 1 4 ∑ p,α (w α m p α -η α )u α S p , (37)
B[ m, Q] ≡ log ∑ S p ∈{-1,1} E v e √ κσ 2 2 ∑ p̸ =q Qpq τ p τ q +κ -1 4 ∑ p,α (w α mp α -θ α )v α τ p . ( 38
)
In order to avoid more cumbersome computations we will skip the details, the interested reader being referred to Ref. [36]. The phase diagram of the model is based on the behaviour of the order parameters Q, Q, m and m. After taking the saddle point of the free energy in the limit L → ∞ keeping κ fixed, using the replica symmetric ansatz and letting the number of replica go to zero, it is possible to distinguish different phases according to the values of the order parameters solutions to the saddle point equations. The order parameters of the systems in the replica symmetric phase are the condensation over the SVD modes (both for the visible and hidden nodes) mα and m α and the overlaps q and q. The saddle point equations of the free energy leads to the following self-consistent equations for the order parameters:

m α = κ 1 4 E v,x v α tanh h(x, v) q = E v,x tanh 2 h(x, v) , ( 39
) mα = κ -1 4 E u,x u α tanh h(x, u) q = E u,x tanh 2 h(x, u) , (40)
where

h(x, u) = κ 1 4 σ √ qx + ∑ γ (w γ m γ -η γ )u γ , h(x, v) = κ -1 4 σ √ qx + ∑ γ (w γ mγ -θ γ )v γ .
A first look at the equations for the magnetization over the mode α tells us that they correspond to the usual mean-field equations of the Sherrington-Kirkpatrick model [66] projected on the SVD decomposition of the weight matrix. The same is true for the overlap, with the difference that we have an overlap parameter for each layer. Analyzing these equations, we can distinguish three phases.

• A Paramagnetic phase: it correponds to the case where q = 0, q = 0, m α = 0 and mα = 0. In the high temperature phase there exists only one minimum to the free energy.

• A Ferromagnetic phase: given by q, q, m α , mα ̸ = 0. In this phase, the magnetization of the system is polarized toward one or many modes α.

• A Spin Glass phase: where q, q ̸ = 0, but m α = mα = 0. In that phase, the system is trapped into one of the many minima of the free energy that are completely uncorrelated with the SVD modes of the weight matrix.

xxxx-17

A c c e p t e d M a n u s c r i p t

On left panel of Figure 4 is shown the phase diagram as a function of 1/σ and of w max /σ , the ratio of strongest mode of 𝑤 to the variance σ of the noise. Fig. 4. Left: the phase diagram of the model. The y-axis corresponds to the variance of the noise matrix, the x-axis to the value of the strongest mode of 𝑤. We see that the ferromagnetic phase is characterized by having strong mode eigenvalues. In this phase, the system can behave either by recalling one eigenmode of 𝑤 or by composing many modes together (compositional phase). For the sake of completeness, we indicate the AT region where the replica symmetric solution is unstable, but for practical purpose we are not interested in this phase. Right: An example of a learning trajectory on the MNIST dataset (in red) and on a synthetic dataset (in blue). It shows that starting from the paramagnetic phase, the learning dynamics brings the system toward the ferromagnetic phase by learning a few strong modes.

From the learning perspective, the interesting phase is the ferromagnetic one. It seems also important that the learning avoid entering into the spin glass (SG) phase. The SG phase, apart from being uncorrelated with the SVD of 𝑤, can affect very badly the MCMC that is used to compute the gradient. By inspecting the phase diagram on the left panel of Figure 4, we understand that at the beginning of the learning it is important to start with a weight matrix with a small variance σ in order to avoid starting from the SG phase. Then, we expect during learning that one or many eigenvalues w α will be expressed and that the trajectory will drift toward the ferromagnetic phase.

The nature of the ferromagnetic phaseit is instructive to look more in details at the ferromagnetic phase to understand the behavior of the RBM. We can distinguish two cases: in the first one only one eigenvalue w α is learned (w α > σ ) and the other ones are close to zero; in the second scenario, many eigenvalues are expressed. In fact the first case is quite simple. Since only one mode has been learned the system will condense along this mode and it will be very similar to a ferromagnet. In the second scenario, we may have many w α that have been learned, i.e. which are above noise threshold. The question then is whether the system will preferentially condense along one single mode taken out of the learned ones or whether it will be able to make compositions by condensing on several modes at the same time. In order to analyze this second scenario, it is important to recall that in order to derive the phase diagram one has to perform the quenched averaging over the matrices 𝑢 and 𝑣. The results will depend on the distribution that is used for the averaging. In Ref. [36], it is shown that depending on the kurtosis of the distribution taken over 𝑢 and 𝑣, the system can behave in different ways. Denoting with γ the relative kurtosis (w.r.t. the normal distribution) three different behaviors are identified:

• γ = 0, e.g. the Gaussian distribution. In that case, only the strongest mode is stable, and the weaker ones are unstable w.r.t.

to the strongest one. Here, the system will condense along the strongest mode only.

• γ < 0 e.g. the uniform or the Bernoulli distribution. Here the weaker modes can be metastable if they are not "too far away" from the strongest one. However the system will condense only toward one mode.

• γ > 0 e.g. a sparse Bernoulli, or the Laplace distribution. In that case, the strongest mode is unstable w.r.t. weaker ones, leaving the possibility to have condensation over many modes at the same time. This corresponds to a dual compositional phase, by reference to the terminology introduced in Ref. [57] which corresponds to combination of features instead of modes.

Hence depending on the form that will take the matrices 𝑢 and 𝑣 during the learning, different types of condensation may appear. This give us some insight on the way the statistical properties of the SVD of the weight matrix are reflected on the recall phase. In some cases the system might recall one macroscopic state, in another one an equilibrium state can be made of a mixture of modes. We illustrate on the right panel of Figure 4 the learning trajectory on the phase diagram obtained both on artificial and MNIST data.

Learning RBM

Let us now discuss possible mechanisms at work during the learning of a RBM, which as we expect should have something to do with pattern formation mechanisms. [START_REF] Amari | A c c e p M a n u s c r i p t Chin[END_REF] We start by summarizing what is understood in exactly solvable models such as the Gaussian-Gaussian and Gaussian-Spherical RBMs. Then we will review a recent work [68] showing how the learning dynamics on a simple dataset for the Bernoulli-Bernoulli RBM with one hidden node can be cast into a spatial diffusion equation. Then we will investigate numerically the behavior of the RBM on the MNIST dataset. In particular, how the learned features at short time can be interpreted using the SVD of the weight matrix and how, at later time, they seem to change completely. Then leaving aside the classical approach based on the Monte-Carlo computation of the gradient -contrastive divergence, [START_REF] Hinton | [END_REF] persistent contrastive divergence, [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF] parallel tempering [START_REF] Tubiana | Restricted Boltzmann machines: from compositional representations to protein sequence analysis[END_REF]69,[START_REF] Desjardins | Parallel tempering for training of restricted Boltzmann machines[END_REF] -we will show how to use the MF self-consistent equations in order to compute the negative term to perform the Finally, we will focus on the ensemble average equations for the learning, where we show how the MF theory developed in Section 4.2. can be integrated numerically and lead to the learning curve of the weight matrix 𝑤.

Learning dynamics for exactly solvable RBMs

Gaussian-Gaussian RBM -We have already seen, in Section 3.1. that the gradient of the Gaussian-Gaussian RBM can be computed exactly and how to characterize the growth of the eigenmodes of the weight matrix when freezing the rotation of the matrices 𝑢 α and 𝑣 α . We put additional results here, operated on an artificial dataset [37] containing 4 well-separated Gaussian clusters. Recall that the modes of the SVD of the dataset that are higher than the intrinsic variance of the visible modes σ 2 v will be expressed, and the vectors of rotation 𝑢 α will aligned themselves with the principal directions of the dataset owing to Eqs. (14,16,17). We can observe on Figure 5 the learning curve obtained for the first eigenmodes of the system coming out of the bulk. We can also see that the first eigenvectors 𝑢 α , associated to the expressed eigenvalues of 𝑤, are aligning with the first principal directions of the SVD of the dataset. In parallel, we see that the likelihood -that can be computed exactly here-of the system increases stepwise after each new mode is learned. v are learned and reach the threshold indicated by Eq. ( 15). In the inset, the alignment of the first four principal directions of the matrix 𝑢 α of the SVD of 𝑤 and of the dataset. In red, we observe that the likelihood function is increasing each time that a new mode emerges.

Gaussian-Spherical RBM -In the case of the Gaussian-Spherical case, it is again possible to obtain an exact analytical expression for the response function of the RBM ⟨s α τ β ⟩ [40] for both the positive and negative terms, where the average is performed respectively over the dataset and the model distribution. The qualitative pictures is very similar to the previous one. As for the linear model, linear correlations between different modes vanish and therefore the matrix 𝑢 has to rotate until it is properly aligned with the principal directions of the dataset. At the same time singular values get either amplified or damped. In contrary to the linear case they do not evolve independently. Instead, as seen on the left panel of Figure 6 lower modes willing to condense exert some pressure on higher modes and accumulate at the top of the spectrum, hence pushing the whole spectrum upward. On the right panel of Figure 6, to illustrate the result of mode condensation, we show a scatter plot containing data from the training dataset and data generated on the trained model when two modes condense.

xxxx-19

A c c e p t e d M a n u s c r i p t Fig. 6. Left: the learning curves for the modes w α using an RBM with (N v , N h ) = (100, 100) learned on a synthetic dataset distributed in the neighborhood of a 20d ellipsoid embedded into a 100d space. Here the modes interact together: the weaker modes push the stronger ones higher, and they all accumulate at the top of the spectrum, as explained in Section 3.2. . Right: a scatter plot projected on the first two SVD modes of the training (blue) and sampled data from the learned RBM (red) for a problem in dimension N v = 50 with two condensed modes. We can see that the learned matrix 𝑢 captures relevant directions and that the RBM data perfectly similar to the one of the training set.

Pattern formation in the 1D Ising chain

In a recent work, [68] the formation of features is studied analytically on a RBM with one hidden node. The training dataset is generated from a 1D Ising chain with a uniform coupling constant and periodic boundary conditions. The model used for generating the data has a translational symmetry which is exploited to solve the learning dynamics exactly. There is indeed available a closed form expression for the correlation function. Thanks to the translation invariance this depends only on the relative distance between the variables. Numerically, it is found that:

• the weights 𝑤 function of the visible node index have a peak value for one of the visible node and decay with distance to this node. Since the position of the center breaks the translation symmetry it tends to diffuse over the system during the learning.

• Using more hidden nodes (but still few), it is observed that each feature is peaked at different places an repel each other to encode the correlation patterns of the data. Again, the position of the peaks diffuse with time even though some repulsive interaction seems to forbid them to cross. See Figure 7 taken from [68] illustrating this phenomena. Right: figure from Ref. [68], the position of the three peaks as a function of the number of training epochs. We observe that the peaks diffuse while repelling each others. The diffusion aims at reproducing the correlation patterns of the translational symmetry, while the repelling interaction ensure that two peaks will not overlap.. Now in Ref. [68], the author compute the gradient of a system with one hidden node

∂ log ℒ ∂ w i = s i tanh β ∑ j s j w j data -tanh(w i )
This expression can be developed up to the fourth order in w (and β ), giving in the case of the 1D Ising chain It is easy to identify in the first two terms the 1D discrete spatial diffusion. This equation can be cast into a spatial diffusion equation with additional term in the continuous time limit (see Ref. [68] for more details). From this small coupling expansion it is also possible to study the stationary solution in the one hidden node case and show that it is consistent with experimental results: it describes a peaked function decreasing rapidly as the distance from the center increases. An approximated weak coupling equations can also be derived in the case of two hidden units. In that case, an effective coupling between the two features vectors w 1 and w 2 is present and responsible for a repulsive interaction between the two peaks. This illustrates nicely how the features learned by the RBM tend to describe local correlations between variables. In addition, these features diffuse over the whole system during the learning to restore the translational symmetry without crossing thanks to a repulsive interactions between them. In the next section, we will focus on the learning behavior on the MNIST dataset and see that in that case, the learned features similarly describe local correlations.

∂ log ℒ ∂ w i ≈ β (w i+1 + w i-1 ) -w i ∑ k w 2 k + w 3 i + 𝒪(w 4 , β w 3 ) xxxx-20

Pattern formation in MNIST: from SVD to ICA?

The pattern formation mechanism can be studied numerically on the MNIST dataset. MNIST [8] is one of the most used real dataset in Machine Learning, it contains 60000 images of black and white handwritten digits of 28 × 28 pixels, ranging from 0 to 9. The digits are about all the same size and are at the center of the image. They are illustrated on Figure 8. To investigate how the patterns emerge from the learning process, we inspect the features during the learning on the Bernoulli-Bernoulli RBM. The first phase of the learning can be understood thanks to a standard linear stability analysis. [36,37] For this let us recall the learning behavior of the Gaussian-Gaussian RBM analyzed in Section 3.1. . In this simple case, the learning was triggered by the SVD of the dataset, and the growth of the modes w α was controlled by how strong was the mode projected in the principal direction of the matrix 𝑢. Consider now the Bernoulli-Bernoulli RBM with {±1} visible and hidden variables (to simplify), and expand the log-likelihood gradient in the limit of small 𝑤 (putting the local biases to zero):

∂ ℒ ∂ w ia = 1 M ∑ d s (d) i tanh ∑ j s (d) j w ja -⟨s i τ a ⟩ ℋ ≈ 1 M ∑ d s (d) i ∑ j s (d) j w ja -w ia = ∑ j C i j w ja -w ia
If we project these equations on the SVD modes of 𝑤 as in Section 3.1. , we obtain the learning dynamics (14) in the Gaussian-Gaussian case, when σ v = σ h = 1. Hence, at the beginning of the learning, this RBM follows the same trajectory as the Gaussian-Gaussian one, where the mode w α are amplified by the principal modes of the dataset. Similarly, it can be shown that the matrix 𝑢 will start to align with the principal direction of the dataset. To see how the features evolve in the non-linear regime, we train an RBM with a very low learning rate and 500 hidden nodes on MNIST.

dw α dt = w α ⟨ ŝ2 α ⟩ -
On Figure 9 we observe as expected from the linear stability analysis, that at the beginning of the learning the first modes of the weight matrix are almost identical to the one of the SVD of the dataset. We see in particular that the features themselves correspond to modes of the dataset, meaning that the RBM starts by learning global features.

Fig. 9. Left: the first 10 modes of the MNIST dataset (top) and the RBM (bottom) at the beginning of the learning. The similarity between of them is clearly visible. Right: 100 random features of the RBM at the same moment of the learning. We can see that most features correspond to a mode of the dataset when comparing with the left-top panel.

Additionally, at this stage of the learning the MC samples obtained from the RBM are typically prototypes: each sample is almost identical (or have a large overlap) with a learned feature. In fact, during the training, if we monitor samples at each epoch (keeping a low learning rate), we can see that the samples have a high overlap with one mode at the beginning, then later on with combinations of modes. To be more precise, we can distinguish different stages of the learning by inspecting the features, the produced samples and the distance between the discretize features (taking the sign of each feature and computing the overlap).

We illustrate these different stages on Figure 10.

xxxx-22

A c c e p t e d M a n u s c r i p t Fig. 10. The column represents respectively (i) the first hundred learned features, (ii) the histogram of distances between the binarized features: W ±1 = sign(W) and (iii) 100 samples generated from the learned RBM. The first row corresponds to the beginning of the learning when only one feature is learned. Looking at the histogram, we see that most of the features have a high overlap. Also, the MC samples are all similar to the learned features. On the second row, the RBM has learned many features, and therefore the histogram is wider but still centered at zero. The MC sampling however is only capable of reproducing one of the learned features. On the last row the learning is much more advanced. The features tend to be very localized and the samples correspond now to digits.

Finally at the end of the learning we recover localized features as in the study-case of the previous section. It has been noticed many times that these localized features are very similar to the ones given by an ICA. To which extent this aspect of learning is affected by the dataset that is considered is an open and interesting question. If we push further the learning, we observe that the RBM keeps learning more and more modes. It is not clear if the system enters in another phase (spin-glass or something else) or if it just overfits the dataset. To end up with these numerical experiments, let us look at the spectrum of 𝑤, at the beginning, at an intermediate stage and at the end of the learning. On Figure 11, we see that starting from a Marchenko-Pastur law, coming from the spectrum of a Gaussian random matrix, quite quickly, many eigenvalues get out of the bulk as they are learned by the machine.

To summarize we have identified the following stages:

• Stage 1: at initialization, the features are completely random and therefore the histogram of distances is Gaussian and centered at zero. The spectrum of 𝑤 follows the Marchenko-Pastur distribution. The RBM starts from the paramagnetic phase.

• Stage 2: the RBM enters the ferromagnetic phase, the first strongest mode of the SVD is learned by all features, giving a high positive or negative overlap in the inter-features distances while the generated samples have a high overlap with the learned features.

xxxx-23

A c c e p t e d M a n u s c r i p t

• Stage 3: where many have emerged, but the learned features remain global and close to the modes of the dataset. The histogram of distances becomes much broader but the generated sample correspond basically to the learned features with few variety. The RBM is in a pure Mattis phase analogous to the recall phase of the Hopfield model.

• Stage 4: finally, after a much longer period, we observed that the learned features are much alike an ICA decomposition while the distances between features is still centered in zero but with a much smaller variance. Finally the generated samples look very similar to the provided dataset. The RBM is in a compositional phase, both regarding the features and the modes (the dual one).

• Stage 5: empirically, we observe that the learning of the modes of 𝑤 never stops. Hence, a macroscopic number of modes is expressed and it is not clear anymore what would be the behavior of the machine in this regime, whether this corresponds to a standard spin-glass phase [61] or another unknown disordered phase.

In future works, it could be interesting to understand the mechanism leading to the localization of the features, in particular whether this is related to some specific tail distribution of the weight matrix spectrum. An aspect of RBMs completely absent from the previous description of the learning process is the behaviour of the biases associated to hidden nodes. These are very important since they determine the threshold above which the features are activated and their learning dynamics is quite intertwined with the modes dynamics. This aspect of the learning could be worth studying especially to improve present learning algorithms.

Learning RBM using TAP equations

The difficulty of learning an RBM comes as already said from the negative term which requires to compute the thermal average of correlations between a visible and hidden nodes. In particular, when the machine starts to learn many modes, it becomes more and more difficult to estimate this term correctly using Monte-Carlo methods due to the eventually large relaxation time. In addition, to get a precise measurement, it is necessary to get many statistically independent samples in order to reduce the statistical error.

In this section we will derive the mean-field self-consistent equations that can be used to approximate the negative term by using a high-temperature expansion of the Boltzmann measure. We illustrate the method showing the result of Gabrié et. al [START_REF] Desjardins | Parallel tempering for training of restricted Boltzmann machines[END_REF] where a RBM has been trained by using the TAP equations.

High-Temperature (Plefka) expansion -We review here a famous approach using a high-temperature expansion of the system in order to compute the mean-field magnetization. This method is both very simple to implement and also provides a way to approximate the free energy of the system in the weak couplings regime. Recent successful approaches [START_REF] Gabrié | Training restricted Boltzmann machine via the Thouless-Anderson-Palmer free energy[END_REF][START_REF] Tramel | [END_REF] showed how it is possible to train a RBM using these mean-field equations. For this subsection, we will use {±1} binary variables for simplicity.

For the Ising model, it is well-known that the (na?ve) mean-field (nMF) approximation can be written as a set of selfconsistent equations on the magnetizations, and the associated approximation of the free energy can be computed as a function of these magnetizations:

m i = tanh ∑ j J i j m j + h i , ∀i xxxx-24
F[𝑚] = ∑ i 1 -m i 2 log 1 -m i 2 + 1 + m i 2 log 1 + m i 2 -∑ i< j J i j m i m j -∑ i m i h i
These equation can be translated directly to the case of the RBM, with the only need to specify clearly which variables are the visible and hidden ones. One gets the following:

m i = tanh ∑ a w ia m a + θ i m a = tanh ∑ i w ia m i + η a F[𝑚] = ∑ i 1 -m i 2 log 1 -m i 2 + 1 + m i 2 log 1 + m i 2 + ∑ a 1 -m a 2 log 1 -m a 2 + 1 + m a 2 log 1 + m a 2 -∑ i,a w ia m i m a -∑ i m i η i -∑ a m a θ a
Here, we remind the reader that we use the indices i, j, k for the visible nodes, and a, b, c for the hidden ones. We recognize on the first two lines of the free energy the entropy terms S(m i ) and S(m a ) of the model for respectively the visible and the hidden nodes. Note first that the nMF approximation corresponds to a first order development in β (or in small 𝑤), but it can be generalized to higher orders, recovering the so-called TAP [73] equations at the second order. Second, we can generalize this scheme to any order using the Pfleka expansion. [74,75] Let us demonstrate first how to obtain the first and second order approximation in the case of ±1 variables. To simplify the computation, we center all the terms around their mean value and make the computation for a case without local bias

ℋ = -∑ i,a s i w ia τ a = -∑ ia (s i -m i )w ia (τ a -m a ) -∑ i (s i -m i ) ∑ a w ia m a -∑ a (τ a -m a ) ∑ i w ia m i -∑ ia m i w ia m a
Using this expression, we can follow [75] and compute the magnetization in the infinite temperature limit of the following free energy

-β A = log ∑ {s,τ} exp(-β ℋ + ∑ i λ i (β )(s i -m i ) + ∑ a λ a (β )(τ a -m a )
The relation between the magnetization and the Lagrange multipliers λ are obtained by imposing m i = ⟨s i ⟩ β =0 = λ i (0) and similar constraints for the hidden nodes. Then, we expand the free energy in a high temperature series

-β A = -β A β =0 -β ∂ β A ∂ β β =0 - β 2 2 ∂ 2 β A ∂ β 2 β =0 + . . .
With our Hamiltonian, we can compute the first and second order easily

- ∂ β A ∂ β β =0 = ⟨ℋ⟩ = ∑ ia m i w ia m a - ∂ 2 β A ∂ β 2 β =0 = 1 2 ∑ ia w 2 ia (1 -m 2 i )(1 -m 2 a )
where we used the following identities for the second order computation

∂ 2 β A ∂ β ∂ m i β =0 = -∑ a w ia m a xxxx-25
∂ 2 β A ∂ β ∂ m a β =0 = -∑ i w m i
As show in Ref. [71], the expansion can be easily extended to the third order without a big computational cost due to the particular topology of the RBM. Deriving the free energy obtained at this order w.r.t. the magnetization, we obtain the self-consistent set of equations defining the TAP equations for the RBM

m i = tanh ∑ a w ia m a -∑ a w 2 ia m i (1 -m 2 a ) (41)
m a = tanh ∑ i w ia m a -∑ i w 2 ia m a (1 -m 2 i ) (42)
Hence, a solution of the TAP equations should satisfied Eqs. ( 41) and ( 42) and give us at the same time the approximated free energy associated to this solution:

F[𝑚] = ∑ i S(m i ) + ∑ a S(m a ) -∑ ia w ia m i m a + ∑ ia w 2 ia 2 (1 -m i ) 2 (1 -m a ) 2 (43)
We can now use these mean-field equations to learn the RBM. First, we need to take into account the fact that many solutions to Eqs. ( 41) and ( 42) exist, each one with a given value of the free energy. Hence, the partition function can be approximated by

Z = ∑ γ e -F(m (γ) i ,m (γ) a )
where the sum runs over all the possible solution to the mean-field equations ( 41)- (42), weighted by the free energy given in Eq. (43). Using this approximation in the computation of the likelihood we obtain the following gradient

∂ ℒ ∂ w ia = ⟨s i τ a ⟩ data -⟨m i m a + w 2 ia (1 -m 2 i )(1 -m 2 a )⟩ MF
where

⟨O⟩ MF = ∑ γ O γ e -F γ ∑ γ ′ e -F γ ′ (44)
correspond to the model average over all the solutions of the mean-field equations. We can see here a notable difference with the approach developed in Ref. [71]. In their work, Gabrié et al. runs the sums over all obtained fixed point from the mean-field equations divided by the number of fixed points only. The risk is that if the mean-field equations converge toward a fixed point that is suboptimal (have a high free energy) or even spurious (being a a maximum of the free energy) the estimation of the negative term will be polluted by such fixed points. More details on the Plefka expansion on bipartite Ising model can be found here. [76] As a final remark, let us insist on the fact that, even if the convergence of the TAP equations is not guaranteed, problems of convergence are practically not met in the ferromagnetic phase. On the contrary, such problems occur quite often in the spin glass phase which we wish to avoid in the context of learning the RBM.

Experiment with TAP learning

We show here some results obtained on MNIST using the same parameters as above but with the mean-field approximation taken from. [START_REF] Gabrié | Training restricted Boltzmann machine via the Thouless-Anderson-Palmer free energy[END_REF] Here, the comparison is done using the persistent chain algorithm, where a set of MC chains is maintained all along the learning whenever using CD, nMF or the TAP approximation (in the case of nMF or TAP, the chain is updated using the corresponding self-consistent equations), see Figure 12.

xxxx-26

A c c e p t e d M a n u s c r i p t Fig. 12. Top: -figure taken from [START_REF] Gabrié | Training restricted Boltzmann machine via the Thouless-Anderson-Palmer free energy[END_REF] -Samples taken from the permanent chain at the end of the training the RBM. The first two lines correspond to samples generated using PCD. The second two lines to samples obtained using the P-nMF approximation and the last two, using P-TAP. Bottom: A 100 features obtained after the training, we can see that they are qualitatively very similar to the ones obtain when training the RBM with P-TAP.

First we see that the samples generated by all the three methods are qualitatively similar. Second, the features learned are also qualitatively similar to the PCD case. Therefore, on the MNIST dataset the two machines are hardly distinguished by just looking at the generated samples and learned features indicating that the MF/TAP approximation is working very well. It is also important to point out here that, the advantage of the mean-field approximation in that case does not rely on any speed up with regard to the learning procedure. But, more importantly, it provides complementary tools such as the fixed points as local maxima of the free energy and their associated free energy. For instance, in Ref. [77] a RBM is used as a prior distribution in the context of compressed sensing where the mean-field equations are used to infer equilibrium values of the variables. In Ref. [78], the RBM is used to reconstruct images from partial observations, again using the mean-field formulation to infer the states of the missing information.

Mean-Field learning: ensemble average

The mean-field equations derived in Section 4.2.for the RBM, where the weight matrix is constructed as a low rank decomposition, can be integrated numerically in order to learn the parameters of the RBM. By contrast to the TAP equations described in Section 5.4. which are solved on single instances, they correspond to the ensemble average (over the parameters 𝑢, 𝑣 and the noise), i.e. are meant to represent an average case of learning.

In the approach developed in Ref. [36], using the statistical ensemble defined in Section 4.2.it is possible to have a mean-field estimate of the response functions involved in the gradient of the log-likelihood. For the response term on the data we get ⟨s α τ β ⟩ data = ⟨s α (s β w βθ β )(1q β [𝑠])⟩ data where the parameter q β [𝑠] is a variant attached to mode β of the spin-glass parameter taken as a function of m in equation (39), when the visible nodes are pinned to the dataset (see Ref. [36] for details). The negative term is more complicated to compute. It depends on the fixed points obtained through equations ( 39) and (40) for a given set of model parameters. Once the fixed point are obtained, the response terms of the RBM can be written

⟨s α τ β ⟩ ℋ = 1 Z MF ∑ γ e -L f (m γ , mγ ,q γ , qγ ) m γ α mγ β = ⟨m γ α mγ β ⟩ MF Z MF = ∑ γ e -L f (m γ , mγ ,q γ , qγ )
where γ runs over the set of fixed points; f is the mean-field free energy that can be derived from Eq. (36). These response terms allows one also to compute the skew-symmetric rotation generators of the visible and hidden singular vectors of 𝑤 through We see again the different eigenvalues emerging one by one, and that each newly learned eigenvalue is triggering a jump of the likelihood together with a jump in the number of fixed points. At the end of the learning, the obtained mean-field fixed points are located at the center of each cluster of the dataset, as can be seen on the scatter plots. In Ref. [36], it is also shown that the behavior is qualitatively similar to what is routinely obtained when performing a standard learning based on PCD.

Other mean-field approach

Other approaches using for instance, message-passing technique such as BP have been developed in order to infer the magnetization of RBM instances. These approaches usually are correct in the limit of weak couplings, and can be used on single instance by iteratively updating a set of messages, here 𝒪(N h N v ) until convergence (see for instance [49,[START_REF] Huang | [END_REF] ). In these works, it is shown how BP can be used to infer the magnetization or the free energy in few well-chosen cases. However, as far as practical learning tasks are concerned, it is not clear that this can be used in general when dealing with the ferromagnetic phase, as can be expected when considering structured data. In fact, it has been showed in many works that BP can have very bad convergence properties in a ferromagnetic phase when the underlying factor-graph is not a tree [80] (particularly if the couplings are strong). This would be most probably the case with RBMs. It is also worth mentioning that in the case of the inverse Ising model, BP approaches never manage to succeed because of the convergence problems [80] and the TAP solution was preferred when inferring the couplings. [81,82] However, some attempts [58,83] using BP and the replica theory on a RBM with one hidden unit were done. In that setting, it is possible to compute the marginal over the weight matrix using BP and therefore to compute its maximum likelihood given some observed datapoints. The results tend to show that, as the number of data increases, the learned features become more localized as is observed in many experiments. Managing to extend this result to the case of many hidden nodes would open the possibility to study the pattern formation using message-passing techniques. An even more recent study, [84] using a variational approach to approximate the posterior distribution of the patterns given the RBM and a dataset shows on artificial xxxx-28 data that the patterns are recovered during the learning. Once again, the missing convincing piece in that case is the applicability real dataset and the ability to samples complex distribution.

Conclusion

With this review, we strive at showing that not only is RBM part of a hectic field of study, but it is also an intriguing puzzle with pieces which are missing in order to be able to understand the way these models can/could assimilate complex information/more complex information. While the black box nature of the learning process starts to fade away very slowly, there are still many key aspects that we do not understand or master for such simple models. We try to list interesting leads for the future.

• Learning quality: Despite the fact that we are maximizing a likelihood function (which can not be computed) it is very hard to obtain a good indicator for comparing two learned RBM. Even if many methods exist to compute the likelihood approximately [START_REF] Salakhutdinov | On the quantitative analysis of deep belief networks[END_REF][START_REF] Krause | [END_REF] the obtained scores are in general not commented in regard to robust statistical analysis. If for very hard cases of image generation, it is easy to compare the results by eye inspection, there are no general method that manage to assess the quality of the samples in terms of how well the learned distribution reproduces the dataset distribution. Some recent work [START_REF] Yale | Generation and evaluation of privacy preserving synthetic health data[END_REF] introduced the notion of "ressemblance" and "privacy" that test the geometric repartition of the true data against the generated samples. This could be a first step defining scores according to different criteria (actually, this problem is not specific to the RBM but concerns actually most of the unsupervised learning models (GANs, VAEs, ...).

• The number of hidden nodes: It is striking that we are still unable to have a principled manner of deciding how many hidden nodes are necessary to learn datasets which are not too complex. For instance, on MNIST, it is possible to learn a machine with only 50 hidden nodes and it somehow manages to produce decent samples. The understanding on how much hidden nodes are necessary to reach a given sample quality is completely missing. In addition, the number of hidden nodes influences a lot the learning behavior of the machine, again in a way that is not fully understood.

• The landscape of free energy: When using statistical mechanics to understand RBMs, the natural question that comes in mind is about the landscape of free energy of the learned machine. It is easy to observe the mean-field fixed point obtained in the ferromagnetic phase and that they do correspond to prototypes of the dataset. Still, we do not know how these many fixed points are organized: are there low free energy paths relating them one from each others? do these paths define a network structure or instead separated clusters of low free energy?

• The landscape of learned RBMs: This is a generic question in Machine Learning : what is the landscape of "good" learned machines in parameter space (here the weight matrix). For supervised tasks, some consensus seems to describe a space which is globally flat where all the good model are next to one another. However this is true for deep models, in the case of RBM, apart from the permutation symmetry of the hidden nodes, we have no clue about what this landscape looks like.

• Link between the dataset and the learned features: We have seen that in the Gaussian-Gaussian case there is a direct link between the eigen-decomposition of the dataset and the learned features. However, for the non-linear model, we do not understand how the modes of the weight matrix are linked to the dataset, nor to the associated rotation matrices.

Fig. 1 .

 1 Fig. 1. Bipartite structure of the RBM.

  (d) i = 0, ∀i. Then, reducing slowly the variance of each component of the mixture, we can look for the moment at which point the degenerate solution corresponding to all the centers placed at the center of masses of the dataset becomes unstable. Linearizing the EM xxxx-10 A c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx equations (

xxxx- 11 A

 11 c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

xxxx- 13 A

 13 c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx pattern:

Fig. 3 .

 3 Fig. 3. From Ref. [62]. Overlap with different patterns when varying the dilution factor p (named d on the figure) at low temperature.Left: a case with 3 patterns where we can observe how at small dilution, only one pattern is fully retrieved while the second and third one appear for larger dilution. Right: a case with 6 patterns where the figure is zoomed in the high dilution region where the branching phenomena is occurring and all the overlaps converge toward the same value.

xxxx- 14 A

 14 c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

xxxx- 15 A

 15 c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx 4.2. Mean-field approach using rank K weight matrix

xxxx- 18 A

 18 c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

Fig. 5 .

 5 Fig.5. On this artificial dataset, we observe that eigenvalues that follows ⟨s α ⟩ 2 > σ 2 v are learned and reach the threshold indicated by Eq. (15). In the inset, the alignment of the first four principal directions of the matrix 𝑢 α of the SVD of 𝑤 and of the dataset. In red, we observe that the likelihood function is increasing each time that a new mode emerges.

Fig. 7 .

 7 Fig. 7. Left: figure from Ref.[68], the value of w i for each visible site of a RBM with 3 hidden nodes trained on the dataset of the 1D homogeneous Ising model with periodic boundary condition. We see three similarly peak shaped potentials with a decreasing magnitude of similar order for the three. Each peak intends to reproduce the correlation pattern around a central node, and therefore cannot reproduce the translational symmetry of the problem. Right: figure from Ref.[68], the position of the three peaks as a function of the number of training epochs. We observe that the peaks diffuse while repelling each others. The diffusion aims at reproducing the correlation patterns of the translational symmetry, while the repelling interaction ensure that two peaks will not overlap..

A c c e p

  t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

Fig. 8 .

 8 Fig. 8. A subset of the MNIST dataset.

Fig. 11 .

 11 Fig. 11. (a) Singular values distribution of the initial random matrix compared to the Marchenko-Pastur law. (b) As the training proceeds we observe singular values passing above the threshold set by the Marchenko-Pastur law. (c) Distribution of the singular values after a long training: the Marchenko-Pastur distribution has disappeared and been replaced by a fat tailed distribution of eigenvalues mainly spread above threshold and a peak of below-threshold singular values near zero. The distribution of eigenvalues do not get close to any standard random matrix ensemble spectrum.

  A c c e p t e d M a n u s c r i p tChin. Phys. B Vol.30, No. x (2021)

A c c e p

  t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

Ω u αβ = w β w 2 α -w 2 β 2 α -w 2 βw α w 2 α -w 2 β 2 α -w 2 βFig. 13 .

 2222222213 Fig. 13. Top panel: Results for a RBM of size (N v , N h ) = (1000, 500) learned on a synthetic dataset of 10 4 samples having 20 clusters randomly located on a sub-manifold of dimension d = 15. The learning curve for the eigemodes w α (left) and the associated likelihood function (right-red) together with the number of obtained fixed point at each epoch. We can see that, before the first eigenvalue is learned there is one single fixed point, then as modes are learned, the number of fixed points increases. Bottom panel: Results for a RBM of size (N v , N h ) = (100, 50) learned on a synthetic dataset of 10 4 samples having 11 clusters randomly defined a sub-manifold of dimension d = 5. On the left, the scatter plot of the training data together with the position of the fixed points projected on the first two directions of the SVD of 𝑤. On the right, the projection along the third and fourth axis. The results are shown after learning 5 modes and where 16 fixed points are found (in fact more than the number of hidden clusters.

  A c c e p t e d M a n u s c r i p t Chin. Phys. B Vol. 30, No. x (2021) xxxx

  Chin. Phys. B Vol.30, No. x (2021) xxxxThe gradient w.r.t. the different parameters will then take a simple form. Let us detail the computation of the gradient w.r.t. the weight matrix. By deriving the log-likelihood w.r.t. the weight matrix we get

e d M a n u s c r i p t

  Chin. Phys. B Vol. 30, No. x (2021) xxxx

	xxxx-12 A c c e p t e d M a n u s c r i p t

) .

  1 , Chin. Phys. B Vol. 30, No. x (2021) xxxx identical at first order in α to

	xxxx-21 A c c e p t e d M a n u s c r i p t

* AD was supported by the Comunidad de Madrid and the Complutense University of Madrid (Spain) through the Atracción de Talento program (Ref. 2019-T1/TIC-13298).