Moshir Harsh

Jérôme Tubiana

Simona Cocco

Rémi Monasson

'Place-cell' emergence and learning of invariant data with restricted Boltzmann machines: breaking and dynamical restoration of continuous symmetries in the weight space

Distributions of data or sensory stimuli often enjoy underlying invariances. How and to what extent those symmetries are captured by unsupervised learning methods is a relevant question in machine learning and in computational neuroscience. We study here, through a combination of numerical and analytical tools, the learning dynamics of Restricted Boltzmann Machines (RBM), a neural network paradigm for representation learning. As learning proceeds from a random configuration of the network weights, we show the existence of, and characterize a symmetry-breaking phenomenon, in which the latent variables acquire receptive fields focusing on limited parts of the invariant manifold supporting the data. The symmetry is restored at large learning times through the diffusion of the receptive field over the invariant manifold; hence, the RBM effectively spans a continuous attractor in the space of network weights. This symmetry-breaking phenomenon takes place only if the amount of data available for training exceeds some critical value, depending on the network size and the intensity of symmetry-induced correlations in the data; below this 'retarded-learning' threshold, the network weights are essentially noisy and overfit the data.

.

Our paper is organized as follows. RBM and their learning algorithms are introduced in Section II. We consider the case of a data distrbution with a single invariance in Section III, and with two symmetries in Section IV. A detailed theoretical analysis of the learning dynamics and of the receptive field emerging through the symmetry-breaking transition can be found in Section V. Conclusions and speculative connections with experiments in neuroscience are proposed in Section VI.

II. RESTRICTED BOLTZMANN MACHINES

A. Definition and log-likelihood

I. INTRODUCTION

Many high-dimensional inputs or data enjoy various kinds of low-dimensional invariances, which are at the basis of the socalled manifold hypothesis [START_REF] Bengio | Representation Learning: A Review and New Perspectives[END_REF]. For instance, the pictures of somebody's face are related to each other through a set of continuous symmetries corresponding to the degrees of freedom characterizing the relative position of the camera (rotations, translations, changes of scales) as well as the internal deformations of the face (controlled by muscles). While well-understood symmetries can be explicitely taken care of through adequate procedures, e.g. convolutional networks, not all invariances may be known a priori. An interesting question is therefore if and how these residual symmetries affect the representations of the data achieved by learning models. This question does not arise solely in the context of machine learning, but is also of interest in computational neuroscience, where it is of crucial importance to understand how the statistical structure of input stimuli, be they visual, olfactive, auditory, tactile, ... shapes their encoding by sensory brain areas and their processing by higher cortical regions. Information theory provides a mathematical framework to answer this question [START_REF] Laughlin | A simple coding procedure enhances a neuron's information capacity[END_REF], and was applied, in the case of linear models of neurons, to a variety of situations, including the prediction of the receptive fields of retinal ganglion cells [START_REF] Atick | Could information theory provide an ecological theory of sensory processing?[END_REF], the determination of cone fractions in the human retina [START_REF] Garrigan | Design of a trichromatic cone array[END_REF] or the efficient representation of odor-variable environments [START_REF] Tesileanu | Adaptation of olfactory receptor abundances for efficient coding[END_REF]. In the case of natural images, which enjoy approximate translational and rotational invariances, non-linear learning rules resulting from adequate modification of Oja's dynamics [START_REF] Brito | Nonlinear Hebbian learning as a unifying principle in receptive field formation[END_REF] or sparse-representation learning procedures [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF] produce local edge detectors, such as do independent component analysis [START_REF] Hyvärinen | Independent component analysis: algorithms and applications[END_REF]. These detectors bear strong similarities with the neural receptive fields measured in the visual cortex (V1 area) in mammals.

It is therefore natural to wonder whether the existence of localized receptive fields is a general feature to be expected from representations of invariant distributions of inputs. Gardner's theory of optimal learning for single-layer neural network (perceptron) predicts that spatially correlated patterns, e.g. drawn from a translationally-invariant distribution, lead to a localized profile of weights [START_REF] Monasson | Storage of spatially correlated patterns in auto-associative memories[END_REF]. Further supporting evidence was recently brought by several works, focusing on the production of such receptive fields in the context of unsupervised learning. Learning of symmetric data with similarity-preserving representations [START_REF] Sengupta | Manifold-tiling Localized Receptive Fields are Optimal in Similarity-preserving Neural Networks[END_REF] or with auto-encoders [START_REF] Benna | Are place cells just memory cells? Memory compression leads to spatial tuning and history dependence[END_REF] both led to localized receptive fields tiling the underlying manifold, in striking analogy with place cells and spatial maps in the hippocampus. In turn, such high-dimensional place-cell-like representations have putative functional advantages: they can be efficiently and accurately learned by recurrent neural networks, and thus allow for the storage and retrieval of multiple cognitive low-dimensional maps [START_REF] Battista | Capacity-resolution trade-off in the optimal learning of multiple low-dimensional manifolds by attractor neural networks[END_REF].

The present work is an additional effort to investigate this issue in a highly simplified and idealized framework of unsupervised learning, where both the data distribution and the machine are under full control. Similarly to previous studies [START_REF] Mehta | An exact mapping between the Variational Renormalization Group and Deep Learning[END_REF][START_REF] Koch-Janusz | Mutual information, neural networks and the renormalization group[END_REF], we consider synthetic data with controlled invariances generated by standard models instatistical physics, such as the Ising and XY models. These data are then used to train Restricted Boltzmann Machines (RBM), a simple albeit powerful framework for representation learning, where a layer of hidden (latent) units account for the correlation structure in the data configurations. We show how the receptive fields of the hidden units undergo a symmetry-breaking transition in the space of couplings: units individually cover localized regions of the input space, but concur to tile the space as best as possible, in much the same way as hippocampal place cells do. This symmetry breaking is dynamically restored if we let the training algorithm run for very long times (well beyond the training time needed to saturate the log-likelihood of the test set): while keeping their localized A Restricted Boltzmann Machine (RBM) is a bipartite, undirected stochastic neural network with two layers, see Fig. 1: • the visible layer includes N units v i , i = 1, ..., N, which carry the configurations of data. For simplicity, we assume here that visible units take binary values, v i = ±1.

• the hidden layer includes M units h µ , µ = 1, ..., M, on which are expressed the representations of the data configurations. Hidden, or latent variables h µ can take real or binary values.

FIG. 1. The two-layer structure of RBM, with weigts w iµ connecting N visible units v i to M hidden units h µ . These binary-valued units are subjected to local fields, called b i and c µ for, respectively, the visible and hidden layers.

The model is formally defined by a Gibbs probability distribution over the sets of visible (v) and hidden (h) variable configuration:

p(v, h) = 1 Z e -E(v,h) , where Z = v dh e -E(v,h) (1)
is the partition function, such that p is normalized to unity, and the energy function E(v,h) is given by

E(v, h) = - N i=1 M µ=1 w iµ v i h µ - N i=1 b i v i + M µ=1 U µ ( h µ ) . (2)
In the formula above, w iµ is the real-valued weight (coupling) connecting the hidden unit h µ and the visible unit v i , b i are real-valued bias terms, also called fields and U µ are the hidden unit potentials. We consider two possible choices for U µ :

• For binary (±1) valued hidden units, a regular field term U µ (h µ ) = -c µ h µ similar to the visible units. In that case, Eqn. 2 is a special case of Ising distribution, with only couplings between units belonging to different layers.

• For real valued hidden units, the symmetric double well potential

U µ (h µ ) = 1 2 h 2 µ + θ µ |h µ |.
For θ µ = 0, the potential is quadratic and the corresponding variable is Gaussian and for θ µ < 0 the potential has two minimas at ±θ µ ; this choice of potential effectively interpolates between Gaussian (θ µ = 0) and binary (θ µ → -∞) hidden units [START_REF] Barra | Phase transitions in Restricted Boltzmann Machines with generic priors[END_REF].

Due to the absence of connections between the units within a layer, the conditional probability of hidden units given the visible units factorizes as follows:

p(h|v) = M µ=1 p h µ |I µ (v) , (3)
where I µ (v) = i w iµ v i is the total input received from the visible layer by hidden unit µ in the absence of fields on visible units, and p(h µ |I) ∝ e U µ (h µ )+h µ I . Therefore, sampling from the conditional distribution is simply done by first computing the hidden layer inputs I µ , then sampling independently each hidden unit given its input according to its hidden unit potential. Similarly, the average activity of a hidden unit given the visible units, h µ |v , is a non-linear function of the input I µ (v) ; for binary hidden units, we have h µ |v = tanh( i w iµ v i + c µ ). Therefore, RBM can be viewed as linear-nonlinear model similar to other feature extraction methods such as Independent Component Analysis. Symmetric formulas can be written for the conditional probability of visible units given the hidden units.

In addition, the marginal distribution over the visible units p(v) can be written in closed form:

p(v) = dh p(v, h) = 1 Z e N i=1 b i v i N µ=1 dh µ e -U µ (h µ )+h µ I µ (v) = 1 Z exp N i=1 b i v i + M µ=1 Γ µ I µ (v) -E eff (v)
, where Γ µ (I) = log dh e -U µ (h)+h I is the cumulant generative function, or log Laplace transform, associated to the potential U µ ; for binary hidden units, Γ µ (I) = log 2 cosh(I + c µ ). Note that by construction, Γ µ (I µ ) is the average value of the hidden unit given its input I µ ; therefore the hidden unit potential determines the transfer function of the hidden unit. Importantly, although the joint distribution is pairwise, the marginal distribution is not in general as Γ µ functions are not quadratic. Therefore, RBM generate effective high-order interactions between the units v i , and are capable of expressing complex measures over the visible configurations [START_REF] Le Roux | Representational power of restricted Boltzmann machines and deep belief networks[END_REF][START_REF] Tubiana | Emergence of compositional representations in restricted Boltzmann machines[END_REF].

B. Training algorithm

Training the RBM is the process of fitting the parameters Θ = {w iµ , b i , c µ /θ µ } to maximize the average log-likelihood of the S data items v data assumed to be independently drawn from p(v). While this may be done with the gradient ascent method, calculating the likelihood is computationally intensive as it requires evaluating the partition function, and sampling methods like Markov Chain Monte Carlo (MCMC) in the form of Gibbs sampling are used.

Gradient of log-likelihood

For the model with parameters Θ, the log-likelihood of a single training example

v data is log L(v data |Θ) = log p v data = -E eff (v data ) -log Z = -E eff (v data ) -log        v e -E eff (v)        . (4)
Taking the partial derivative with respect to Θ gives

∂ log L(v data |Θ) ∂Θ = - ∂E eff (v data ) ∂Θ + 1 Z ∂Z ∂Θ = - ∂E eff (v data ) ∂Θ + ∂E eff (v) ∂Θ RBM , (5)
where (. (v) (.) denotes the average according to the marginal distribution over the visible units with parameter values Θ.

) RBM = 1 Z v e -E eff
In particular, for the weights w iµ , we have according to [START_REF] Garrigan | Design of a trichromatic cone array[END_REF],

∂E e f f (v) ∂w iµ = -v i Γ µ I µ (v) ≡ -v i h µ |v .
The gradient of the total log-likelihood is then

∂ log L(v data |Θ) data ∂w iµ = v data i h µ |v data data -v i h µ |v RBM . (6)
Equation ( 6) is an example of moment-matching condition, as it imposes that the correlation between the variables v i and h µ computed from the data coincides with its counterpart defined by the RBM model distribution p(v, h). The gradients of L over b i and c µ lead to similar moment-matching conditions for, respectively, the average values of v i and of h µ .

Approximating the log-likelihood gradient

In the gradient of the log-likelihood of Eqn. ( 6), the model-distribution moment is not computationally tractable, as it requires to sum over all values of the visible and the hidden variables. In practice, an approximate value for this term is obtained by Markov Chain Monte Carlo (MCMC) methods. The Markov Chain is defined by repeated iterations of Gibbs sampling, which consists in sampling h from v and v from h using Eqn. [START_REF] Atick | Could information theory provide an ecological theory of sensory processing?[END_REF]. In principle, one should run a full MCMC simulation at each gradient step, but this is computationally prohibitive. For our RBM training we use the Persistent Contrastive Divergence (PCD) algorithm [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]: Markov Chains are initialized at the beginning of the training and updated with only a few Gibbs Monte Carlo steps between each evaluation of the gradient, see [START_REF] Fischer | An introduction to restricted Boltzmann machines[END_REF] for a more detailed review. This approximation works very well for the data distribution studied here because they are in a paramagnetic phase (= monomodal), hence the Markov Chains mix very rapidly.

Stochastic Optimization

The RBM is trained using Stochastic Gradient Ascent (SGA), the golden standard for neural network optimization. SGA is a variant of ordinary gradient ascent where at each step, only a small subset of the data set (the minibatch), of size B ∼ 10 -100 examples, is used to evaluate the average log-likelihood, see Eqn 7 where ν is the learning rate which dictates how much to change the parameter in the direction of the steepest gradient. The dataset is divided into S /B mini-batches Batch(t), and for each epoch of training t, we perform one SGA update for each mini-batch. An epoch consists of using all the subsets for the update such that each data sample is used once. After every epoch the subsets are again drawn randomly. Several dozens of epochs are usually required to reach convergence.

Θ t+1 = Θ t + ν         1 B b ∈ Batch(t) ∇ Θ log L(v data,b |Θ)         (7)
Compared to ordinary gradient ascent, SGA serves several purposes. First and foremost, its computational cost is significantly reduced as only a small batch is used per update; yet the update is usually reliable thanks to data redundancy. Second, the stochastic evaluation of the gradient introduces noise in the learning process. This prevents the dynamic from getting trapped in local maxima, which is crucial for non-convex optimization landscapes, and it also directs the dynamics toward minima with wider basins of attraction [START_REF] Jastrzkebski | Three factors influencing minima in SGD[END_REF]. It has been argued that the later effect contributes in improving generalization performance [START_REF] Hochreiter | Flat minima[END_REF][START_REF] Keskar | On large-batch training for deep learning: Generalization gap and sharp minima[END_REF][START_REF] Chaudhari | Entropy-SGD: Biasing gradient descent into wide valleys[END_REF]. Though the convergence rate of SGA has a slower asymptotic rate than ordinary gradient descent, it often does not matter in practice for finite data sets, as the performance on the test set usually does not improve anymore once the asymptotic regime is reached [START_REF] Bottou | The tradeoffs of large scale learning[END_REF].

The noise level of the SGD is directly related to the batch size and learning rates parameters, see for instance [START_REF] Smith | A bayesian perspective on generalization and stochastic gradient descent[END_REF]. Briefly speaking, assuming i.i.d. and infinite number of samples, the SGA parameter increment has mean value ν∇ Θ log L(v data |Θ) data , and variance proportional to ν 2 /B; in the large B limit it is also Gaussian distributed according to the central limit theorem. In comparison, the increments of a continuous time Langevin equation with energy landscape E =log L(v data |Θ) data and noise covariance matrix ∝ σ 2 S GA , integrated over a time step ν has the same mean value and a covariance proportional to σ 2 S GA ν. Identifying both noises gives the following scaling law for the SGA noise, σ S GA ∝ ν B . Reducing the learning rate and increasing the batch size therefore decrease the noise level, and vice-versa. In all our experiments, both learning rates and batch sizes are kept fixed throughout a training session.

III. LEARNING DATA WITH A SINGLE INVARIANCE A. Data distribution: Ising model

Our first toy distribution for data v data is the celebrated one-dimensional ising model from statistical physics. Here each v i is a spin which can either be up or down, that is can take only ±1 binary values. The corresponding joint probability distribution of the visible units reads

p data (v 1 , v 2 , ..., v N ) = 1 Z ising e β N i=1 v i v i+1 (8)
where the partition function normalizes this probability over the 2 N visible configurations, and β > 0 is referred to as the inverse temperature. We enforce periodic boundary conditions through v N+1 ≡ v 1 .

As is well known, under distribution (8), all visible units v i have average values equal to zero, and the correlation function decays exponentially with the distance separating the corresponding units on the ring,

v i v j = v p data (v 1 , v 2 , ..., v N ) v i v j = e -|i-j|/ξ , where ξ = 1 ln coth β (9)
is the correlation length. The above expression for the correlation holds when N ξ. Formula (8) defines a simple example of invariant distribution under the set of translations (or, better, rotations) on the N-site ring. More precisely, for any integer k, we have

p data (v 1 , v 2 , ..., v N ) = p data (v k+1 , v k+2 , ..., v k+N ) , (10)
where i + k is to be intended modulo N. Figure 2 shows a number of configurations, drawn independently and at random from this probability distribution using the Gibbs sampling algorithm.

FIG. 2. 10,000 one-dimensional Ising model configurations with 100 spins each, sampled from distribution (8) at inverse temperature β = 1 and with periodic boundary conditions. Black and white dots represent units equal to, respectively, +1 and -1. The correlation length ξ may be though of as the typical length of black or white contiguous regions along the horizontal direction.

Here, β = 1, which corresponds to ξ 3.7.

B. Initial learning and emergence of place cells

Case of a single hidden unit

First we train the RBM with only M = 1 hidden unit, and N = 100 visible units. Such a limited machine is, of course, not expected to reproduce accurately the Ising model distribution underlying the data. However, this is an interesting limit

y(i) = {w i -max(w i )/2}
, where i is the site index; The roots of this spline were then determined numerically, and the width was defined as the modulus of the difference between the roots. This procedure reliably finds the Full Width at Half Maximum (FWHM). A linear fit (red line) of the form y = a x + b shows that the width of the place (receptive) field of the only hidden unit is proportional to the relevant characteristic length in the data. Notice that the intercept (b) is non zero, in agreement with the theoretical findings of Section V B in the β → 0 limit.

case to study how the RBM can make the most of its single set of weight attached to the unit. We use a large number of data configurations for training, which makes our distribution approximately invariant under rotations on the ring. We initialize the weights w i1 with small amplitude Gaussian random values; since the data are symmetric, we further impose b i = c µ = 0 ∀i, µ. The results of the training phase after 100 epochs, i.e. the weights w i1 are shown in Fig. 3(a). We observe that the weights are not uniform as could have been naively expected from rotational invariance, but focus on a limited portion (place) of the N-site ring. The position of the peak depends on the initial conditions for the weights; it may also be influenced by the small irregularities in the data set coming from the finite number of training configurations.

To understand what determines the width of the weight peak, we train different RBMs with data at different inverse temperatures β, and calculate their average peak widths over multiple runs. We plot the peak width as a function of the correlation length ξ in Fig. 3(b). We observe that the peak width scales proportionally to ξ. Interestingly, despite its very limited expression power, our single-unit RBM has correctly learned to coarse grain the visible unit configurations on the relevant scale length in the data, ξ. Having wider receptive, or place fields would not be as much as informative. For instance, with a set of uniform weights w i1 = w, the hidden unit would simply estimate the average magnetization of (mean value of all visible units in) the data configurations, which are all equal to zero up to fluctuations of the order of ±N -1/2 , and would completely miss the correlated structure of the data. Conversely, more narrow place fields would have lower signal-to-noise ratios: the strong correlations of visible units over the length ξ allows one to reliably estimate the local magnetization and the correlation structure on this scale.

Case of multiple hidden units

We next show results obtained when training RBM with M = 3 hidden units on the same data. Figure 4 shows that each one of the three sets of weights have roughly the same peaked structure (same width) as in the M = 1 case, but the peaks are centered at different places along the ring. The roughly equal distance between successive peaks shows the existence of an effective repulsion between the weights of any two hidden units. This phenomenon is easy to understand on intuitive grounds: having very overlapping place fields produces highly redundant hidden units,and would not help capturing the spatial correlation in the data spreading over the entire ring.

Training of RBMs with a large number of hidden units shows the same pattern of production of place fields attached to different hidden units, covering in a approximately uniform way the visible space (ring), see Fig. 4(c) in the case of M = 100 hidden units. The only notable difference is that the width of the place fields shrinks as M gets very large. This happens when Mξ N, i.e. when the single-hidden-unit peaks would start to largely overlap.

C. Long-time learning and restoration of invariance through place-field diffusion

We now let the training dynamics evolve for a much larger number of epochs. In the case of a RBM with one hidden unit only, the weight vector shows the overall peak structure of Fig. 3(a) at all times (after a short initial transient during which the localized peak emerges). However, the location of the peak may change on very long training time scales. Figure 5(a) shows ten trajectories of the center of the peak corresponding to ten random initialization of the weights (equal to small values drawn independently and randomly). We observe that the centers of the peaks undergo a seemingly random motion. When the number of data items used for training is very large (to erase any tiny non-homogeneity in the empirical distribution), this random motion looks like pure diffusion. shows that the Mean Square Displacement (MSD) of the peak center grows roughly linearly with the training time (number of epochs), which defines the effective diffusion coefficient of the weight peak. For intermediate numbers of data items, diffusion is activated: due to the inhomogeneities in the empirical data distribution, some places along the ring are preferred, and have a tendency to trap the weight peak for some time.

Repeating the same analysis for a RBM with M = 3 hidden units allows us to observe the diffusion of the three peak centers, see Fig. 6. We see that the motions of these centers are coupled to maintain a constant distance between each other. This is a clear signature of the effective repulsion between the hidden-unit weight vectors already discussed in Section III B 2.

D. Case of few data: Retarded learning transition

The emergence of a pronounced peak in the weight vector attached to a hidden unit reported above takes place only if the number of data items are sufficiently large. For very few data, the RBM weights do not show any clear spatial structure and seem to overfit the data. Similarly, for a fixed number of data samples, a transition is observed between the overfitting and spatially-structured regimes as the correlation length ξ (or the inverse temperature β), that is, the spatial signal in the data is increased. To distinguish these two regimes, we introduce the empirical order parameter

W = N i=1 w i1 , (11)
which is expected to be large when place fields emerge and the weights are spatially structured, and much smaller (and vanishingly small in the large-N limit) in the overfitting regime. Figure 7(a) shows the value of the order parameter W as a function of the intensity of spatial correlations for a fixed number of data samples. For small values of β (and ξ) W vanishes: the very weak spatial structure in the available data is not learned by the RBM. At large β, a place field emerges, focusing on a finite portion of the ring, and W is non zero. The same transition is observed when β is fixed and the number of training samples, S , is varied, see Fig. 7(b). For few samples or, equivalently, large noise levels r = N/S , the RBM overfits the data and W vanishes. For small values of r, W becomes non zero, signalling the emergence of a place field focusing on a finite portion of the ring. This transition is an example of the very general mechanism of the so-called retarded learning phenomenon [START_REF] Watkin | Optimal unsupervised learning[END_REF], also encountered in the context of random correlation matrices and the spiked covariance model. The connection with random matrices will be made explicit in Section V A. The classical XY model is a popular model in statistical physics, used in particular to study topological phase transitions in two dimensions. We consider here the one-dimensional version of this model, which shows no such phase transition but is nonetheless very useful for our study due to the additional symmetry with respect to the Ising model. In the XY model each lattice site i carries an angle θ i ∈ [0, 2π[ with respect to some arbitrary, fixed direction. The energy function reads, up to a scale factor that can absorbed in the temperature definition,

E(θ 1 , θ 2 , ..., θ N ) = - N i=1 cos(θ i -θ i+1 ) (12)
with periodic boundary condition θ N+1 = θ 1 . We then discretize the set of angle values in multiples of 2π/P, where P is an integer. The resulting model is a Potts model over the N integer-valued variables v i = 0, 1, 2, ..., P -1, with probability distribution (with periodic boundary conditions)

p data (v 1 , v 2 , ...., v N ) = 1 Z e β n i=1 M(v i ,v i+1 ) ( 13
)
where the interaction kernel M mimics the XY energy function,

M(v, v ) = cos 2π P (v -v ) , (14)
and the partition function normalizes the distribution p. This distribution enjoys two symmetries, compare to the single symmetry of the Ising model in ( 10): for any integers K and L we have,

p data (v 1 , v 2 , ..., v N ) = p data (v k+1 + L, v k+2 + L, ..., v k+N + L) , (15)
where i + k and v + L are to be intended, respectively, modulo N and P. We first train a RBM with a single hidden unit h 1 , which takes real values and is submitted to a double-well potential. Figure 9(a) shows the weights obtained after training from a very large number of configurations, starting from small white noise initial conditions for the w i,1,v . We observe a strong modulation of the weights in the space and angle directions, achieving peak values around some site i and angle v. Similar results were found for a binary-value hidden unit, h 1 = ±1, with a slightly weaker localization of the weights and at a different location, see Fig. 9(b). In the following, we show results obtained for the RBM with the real-value hidden unit only.

Since the interaction matrix M in the Potts model takes the cosine function form, our RBM should learn the same functional dependence from the data samples. We show in Fig. 9(c) the quantity which measures the angular modulation of the weights on each site i. We see a strong space localization around i = 39, because the weights only take non-zero values near that location. This location is arbitrary and similar to the place-field formation accompanying the breaking of translation symmetry over space observed for the Ising model. In addition, at the location of the maxima, the weight vector is very well approximated by a cosine function, see Fig. 9(d). The RBM has learned the correct frequency equal to 2π/10, and the phase takes an arbitrary value. Indeed, the phase in a free parameter due to the invariance against choices of L in [START_REF] Watkin | Optimal unsupervised learning[END_REF].

W angular i = P-1 v=0 w i,1,v , (16)
To obtain a more precise picture of the receptive field, we then consider, for each site i, the P-dimensional vector of the weights w i,µ=1,v . We then fit this vector with a cosine function of adjustable frequency and phase, referred to as, respectively, ω i and ϕ i . We show, as functions of the site index i, the periods ω i and the phases ϕ i in Fig. 9(e). We observe that the period takes the expected value 2π/P over the receptive field (sites ranging approximately between i = 30 and 50). Similarly, the phase is constant (and takes an arbitrary value) over the same region of space. Informally speaking, the hidden unit is on, all the XY spins supported by the sites in the receptive field point to the same direction.

Case of multiple hidden units

We also train a RBM with M = 5 Real valued hidden unit with double well potential, with results shown in Fig. 10. We see that the receptive fields of the hidden units are mutually separated in space, and show the same phenomenon of repulsion between the units observed for the Ising data. In addition, the angular dependence of the five weight vectors exhibit the same frequency (equal to 2π/10), but the phases show also a nice equi-separation due to repulsion along the angular direction.

Though we expected to see a diffusion of the receptive fields both along the spatial and angular dimensions for very large learning times, we did not observe this phenomenon even with RBM trained with 1,000,000 samples. This is likely due to the fact that the landscape is still rough for this amount of data, and diffusion remains activated. We have not tried to increase the number of samples because of the high computational cost. (similar to ( 16), but the hidden-unit index µ runs from 1 to 5) vs. space location i. All the parameters are same as those of Fig. 9 but with M = 5 Real valued hidden unit with double well potential. There is a clear strong space localization and the receptive fields for the different units show mutual repulsion of their weight vectors. (b) Weight vectors at the maximum visible unit index i max (µ) for the respective hidden units µ. The different curves each lie on a cosine function with frequency = 2π/10 as expected but with phases showing equal separations when ranked in increasing order.

C. Differentiated retarded learning transitions

In this section, we show that RBM trained with data generated by the discretized XY model shows retarded learning phase transitions. However, as there are two potential symmetry breaking directions in this model, one corresponding to the angular space and the other to the positional space, the breaking of symmetry along these direction may take place at two different values of the noise ratio r = N/S , i.e. for different number of samples in the data set used for training. The reason is that the number of Potts states in the angular direction, P, may largely differ from the number of sites on the lattice, N. Consequently, the effective system sizes along the two directions are different.

This phenomenon of differentiated retarded learning phase transitions is reported in Fig. 11. We show in panel (a) of the figure the spatial modulation defined through,

W spatial v = N i=1 w i,1,v , (17)
as a function of the Potts angular state variable v. We observe that for large r, the spatial modulation vanishes all over the angular space: low amount of data are not sufficient for the RBM to capture the angular correlations in the configurations. For large enough data set (r < 0.033) the spatial modulation shows a clear dependence on v. We then show in panel (b) of Fig. 11 the angular modulation W angular i as a function of the lattice site index i for varied levels of sampling noise, r. Again, for large r, no modulation is seen. However, for very small noise levels r < 0.002, we do observe that W angular i is peaked around some well defined site i. Interestingly, in the range 0.002 < r < 0.008, the angular modulation does not significantly vary over space, while the spatial modulation varies over angles, compare panels (a) and (b). We conclude that, for intermediate ranges of values of r, the RBM has created a place-field along the angular direction, but not along the spatial direction.

To test the generality of the phenomenon of differentiated transitions, we also generated data samples from variants of the discretized XY model. We modified the XY model in terms of changing the interaction matrix M in ( 13) from the cosine function to short range couplings, and also we changed the Hamiltonian to include not only nearest neighbor couplings but also long range couplings in the positional space. The resulting models display a variety of phase transitions in the RBM weights after training, with positional symmetry breaking arising before (for smaller amount of training data) angular ordering in some cases (not shown).

V. THEORETICAL ANALYSIS

Hereafter, we study analytically the dynamics of learning of the weights of the RBM with binary hidden units when trained with data. Two limit cases will be considered:

• The case of few data, which allows us to establish the connection with random matrix theory and the so-called retarded learning transition;

• The case of a large amount of data, with weak correlations, which we analyze in detail to understand the formation and shape of the place field, as well as the interactions between different place fields arising through learning.

While we will focus on the learning dynamics of the weights, we assume that the RBM has correctly learned the local fields, so we will set b i = c µ = 0 from the beginning in the case of unbiased binary data v i = ±1. In addition, we assume that hidden units are of Bernoulli type, h µ = ±1. The log-likelihood therefore reads

log L = M µ=1 log cosh        N i=1 w iµ v i        data -log Z({w iµ }) , (18)
where the partition function is

Z({w iµ }) = {v 1 ,v 2 ,...,v N } M µ=1 cosh        N i=1 w iµ v i        . ( 19
)
Taking the partial derivative with respect to w µi we get the following expression for the gradient of the log-likelihood:

∂ log L ∂w iµ = v i tanh         N j=1 w jµ v j         data - 1 Z({w iµ }) {v 1 ,v 2 ,...,v N } v i sinh         N j=1 w jµ v j         λ( µ) cosh         N j=1 w iλ v j         . (20)
The continuous-time dynamical equations for the evolution of the weights during training, assuming that the batch size is maximal, i.e. that all the data are used for training, are

dw iµ dt = ν ∂ log L ∂w iµ , ( 21
)
where ν is the learning rate.

A. Few data: Small weight expansion and the retarded learning transition

Linearized equations of the dynamics

In this Section, we assume that the weights have initially very small (random) values. For small enough learning times, we may linearize the dynamical equations [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]. We obtain

dw iµ dt = ν N j=1 C i j w jµ -w iµ , (22)
where

C i j = v i v j data (23)
is the empirical covariance matrix estimated from the data. Let Λ be the largest eigenvalue of C, and e the associated eigenvector, with components e i . As the diagonal elements C ii are equal to unity (v 2 i = 1), we have that Λ > 1, unless C is the identity matrix and the data shows no correlation at all. Hence, according to [START_REF] Fischer | An introduction to restricted Boltzmann machines[END_REF], all weight vectors w µ = {w 1µ , w 2µ , ..., w Nµ } align along e; this result holds within the linear approximation, and is therefore expected to be valid at short times only.

Let us consider the noise ratio r = N/S , equal to the number of visible units (system size) over the number of training samples. For bad sampling (large r), the empirical covariance matrix can be approximated by the covariance matrix of a null model, in which all N visible units are independent and unbiased: v i is equal to ±1 with equal probabilities (= 1/2), independently of the other v j 's. The asymptotic distribution of the eigenvalues of such a random matrix has a special form, called the Marcenko Pastur (MP) spectrum, whose right edge (top eigenvalue) is given by

Λ noise = Λ MP = 1 + √ r 2 ( 24
)
and the corresponding top eigenvector e has random, Gaussian distributed components. Conversely, for good sampling (small r), we expect the empirical covariance to be similar to the covariance matrix computed from the model distribution p from which data were generated. Due to the translational invariance of p, its top eigenvector e Model has the same symmetry: e model = (1, 1, ..., 1), up to a normalization factor. Hence, we expect e to be similar to e Model and be roughly uniform. In the double, large N and large S limit, the two regimes may be separated by a sharp transition, taking place at a critical value of r. To locate this value, we compute below the top eigenvalue of the model covariance matrix, and compare it to its MP counterpart [START_REF] Hochreiter | Flat minima[END_REF]. The crossover between the bad and good sampling regimes takes place when both eigenvalues are equal.

Case of Ising data

Let us consider the case of the one-dimensional Ising model. When a large number of configurations is available, we have C i j = (tanh β) |i-j| , see [START_REF] Monasson | Storage of spatially correlated patterns in auto-associative memories[END_REF]. Due to the rotational invariance, the top-eigenvector, e I sing has all its components equal. Therefore the top eigenvalue of the covariance matrix is

Λ I sing (β) = N j=1 C i j = 1 + 2 tanh β + tanh 2 β + tanh 3 β + ... ≈ 1 + tanh β 1 -tanh β = e 2 β . (25)
When the inverse temperature β is small, this 'signal' eigenvalue is smaller than the 'noise' eigenvalue Λ MP [START_REF] Hochreiter | Flat minima[END_REF], locating the right edge of the MP spectrum. In this case, we expect the top eigenvector e of the empirical covariance matrix C to be noisy, and not to capture the correlation between the Ising variables v i . In this regime, the RBM overfits the data and no receptive field with a localized weight structure can emerge. As β increases above

β(r) = log 1 + √ r , (26)
the signal eigenvalue Λ I sing (β) becomes larger than the MP edge, and we expect the top eigenvector of C to have comparable component and be similar to e I sing . The above statement is corroborated by the results shown in Fig. 12(Top), which shows the top eigenvalue Λ of the correlation matrix C (23) as a function of the noise ratio, r = N/S , where S is the number of samples. For large r (few samples), Λ is very well approximated by Λ noise , while, for small r (many samples), Λ gets very close to Λ I sing (β) as expected. The crossover between these two regimes takes place at values of r such that β β(r) [START_REF] Chaudhari | Entropy-SGD: Biasing gradient descent into wide valleys[END_REF]. Figure 12(c) shows how the top eigenvector of the data correlation matrix changes as more and more samples are considered. One clearly sees a phase transition from a random vector to the uniform eigenvector e I sing .

Case of XY data

For the discrete XY model, the correlation matrix in the r → 0 limit can be computed as well using the transfer matrix formalism. We find

C i j (v, v ) = 1 P -1 P-1 p=1 λ p (β) λ 0 (β) | j-i| cos 2πp(v -v ) P , (27)
where

λ p (β) = P-1 v=0 exp β cos 2πv P cos 2πvp P . (28)
C enjoys translational invariance along both axis, hence its eigenvectors are discrete 2D Fourier modes; after computation, we find that the top eigenvalue is

Λ XY (β) = P P -1 λ 0 (β) + λ 1 (β) λ 0 (β) -λ 1 (β) . (29)
with a corresponding eigenspace of dimension 2, spanned by

e 1 i (v) = 2 NP cos 2πv P , e 2 i (v) = 2
NP sin 2πv P . The top eigenvector is uniform over space, as for the Ising model, but not over the angular variables, see Fig. 13(c).

The 'noise' eigenvalue is similarly given by the MP spectrum, although slightly modified: the dimension to sample size ratio is now PN S = Pr, and in the S → ∞ limit, the correlation matrix has top eigenvalue P P-1 owing to the anticorrelations between Potts variables on the same site, C i,i (a, b) = -1 P-1 , ∀ a b. We obtain:

Λ noise = P P -1 + 1 + √ rP 2 (30)
Similarly to the case of Ising data, when β is small, the signal 'eigenvalue' is small compared to the 'noise' eigenvalue, and the empirical top eigenvector has a small projection in the space spanned by e 1 ,e 2 , see Fig. 13(b,c,d). The crossover between the two regimes takes place at values of r such that Λ MP Λ XY (β). The first retarded learning transition of the RBM occurs in the same range of r, see Fig. 11.

B. Many data: Small β expansion

After some training time, linear equation ( 22) for the weights breaks down and non linearities must be taken into account [START_REF] Decelle | Spectral dynamics of learning in restricted Boltzmann machines[END_REF]. We derive below an approximation to the RBM dynamic learning equation (with M = 1 or 2 hidden units) for the onedimensional Ising models, which is exact for small (but non vanishing) inverse temperature β. We show that this equation is free of any external parameters after appropriate rescaling of the weights. We compare the numerical solutions to this equation with the result of the training with RBM to find a parameter independent agreement with the shape and the structure of the weights. We also cast the equation into a continuous form, and formulate the system in terms of a standard Reaction-Diffusion instability problem with the weights as an inducer and the sum of weights squared as the repressor.

One Hidden Unit System: formation of receptive field

For one hidden unit, equations [START_REF] Tubiana | Emergence of compositional representations in restricted Boltzmann machines[END_REF][START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF] become, after some elementary manipulation,

∂ log L ∂w j = v i tanh         N j=1 w j v j         data -tanh w i , (31)
where we have dropped the µ = 1 index for the sake of clarity. Expanding the hyperbolic tangents to the third powers of their arguments, we obtain

∂ log L ∂w j = j v i v j data w j - 1 3 j,k,l v i v j v k v l data w j w k w l -w i + 1 3 w 3 i + O(w 4 ) . (32)
Let us now assume that a large number of samples is available. At the lowest order in β, we have

v i v j =          1 if i = j , β if i = j ± 1 , 0 otherwise . (33)
and for, respectively, the 2-and 4-point correlations. we therefore obtain

v i v j v k v l =              1 if i = j, k = l or any permutation , β if i = j ± 1, k = l or any permutation , β if i = j, k = l ± 1 or any permutation , 0 otherwise . (34)
∂ log L ∂w j = β(w i+1 + w i-1 ) -w i k w 2 k + w 3 i + O(w 4 , β w 3 ) . ( 35
)
Upon appropriate rescaling of the weights, w i → w i / √ β, and of the learning rate, ν → ν/β, we obtain, in the small β regime, the non trivial, parameter-free dynamical equation 1 ν

dw i dt = w i+1 + w i-1 -w i k w 2 k + w 3 i . ( 36
)
The stationary solution of this equation is shown in Fig. 14(a). This equation can be cast in a continuous form over space, where we use the Laplacian to describe spatial diffusion. The corresponding continuous partial differential equation reads

1 ν ∂w ∂t (x, t) = ∂ 2 w ∂x 2 (x, t) + 2 -b(t) w(x, t) + w(x, t) 3 , (37)
where

b(t) = L 0 w(x, t) 2 dx . ( 38
)
These coupled dynamical equations lead to a non-trivial spatial formation through the so-called Turing Reaction Diffusion instability mechanism [START_REF] Murray | Mathematical Biology. II. Spatial Models and Biomedical Applications[END_REF]. The field w(x, t) diffuses over space and activates itself (self-promoting, through the cubic term), but is inhibited by another species, b. This repressor is diffusing with an infinite diffusion coefficient, i.e. is spatially uniform, and depends on w through (38). As w grows due to self-activation, so does the repressor b, until w reaches a stationary profile. We show in Appendix A that the above dynamical equation satisfy the general criteria for stable pattern formation.

Multiplying by dw dx on both sides and integrating over x we obtain that

E(x) ≡ 1 2 dw dx 2 + w(x) 2 + 1 4 w(x) 4 - b 2 w(x) 2 (40)
has a uniform value E 0 , independent of x. When x → ±∞, both w(x) and dw dx (x) tend to 0, which sets E 0 = 0. We deduce that

dw dx (x) = ±w(x) b -2 - w(x) 2 2 . (41)
We now explicitly break the symmetry by fixing the centre of the peak of the weights in x = 0, with w(0) > 0, dw dx > 0 for x < 0, and dw dx < 0 for x > 0. Imposing that the derivative of the weight with respect to x vanishes at its maximum, i.e. that w is twice differentiable in x = 0 gives

w(0) = 2(b -2) . (42)
Integrating (41) with condition (42), we find

w(x) = √ 2(b -2) cosh x √ b -2 . ( 43
)
Using definition (38) for b we then find

b = ∞ -∞ w(x) 2 dx = 4 √ b -2 , (44)
whose solutions are b ± = 8 ± 4 √ 2. The corresponding profiles of the weights are shown in Fig. 14(c). We now study the stability of the solution under the time-dependent perturbation w(x) → w(x) + (x, t), where w(x) is given by (43). According to equations [START_REF] Chen | How vision and movement combine in the hippocampal place code[END_REF] and (39), we have

1 ν ∂ ∂t (x, t) = ∂ 2 ∂x 2 (x, t) + 2 -b (x, t) -2 dy w(y) (y, t) (x, t) + 3 w(x) 2 (x, t) . (45)
Multiplying by w(x) and integrating over x, we get the following equation

1 ν d dt dx w(x) (x, t) = - √ b -2 dx w(x) (x, t)        8 - b cosh x √ b -2 2        . (46)
We deduce that the weight profile is stable if and only if b < 8. Therefore, the b = b + solution is unstable against small variations of the peak amplitude near x = 0, and the solution b = b -is the correct, stable one. Notice that the width of the peak of the weight, in the β → 0 limit is finite according to expression (43). This phenomenon was also observed by the RBM training results in Fig 3(b), where the peak width obtained by linear fit (coefficient b) was also positive and finite.

Two hidden units: interaction between receptive fields

For two hidden units, equations [START_REF] Tubiana | Emergence of compositional representations in restricted Boltzmann machines[END_REF][START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF] become, after some simple manipulation,

∂ log L ∂w j1 = v i tanh         N j=1 w j1 v j         data - tanh(w i1 + w i2 ) 1 + N k=1 cosh(w k1 -w k2 ) cosh(w k2 + w k2 ) - tanh(w i1 -w 2i ) 1 + N k=1 cosh(w k1 + w k2 ) cosh(w k1 -w k2 ) , (47)
together with a similar equation for the weight vector µ = 2 obtained by swapping the hidden-unit indices 1 and 2. Notice that this equation simplifies to [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF] when the weight vector µ = 2 is set to zero, i.e. when w i2 = 0 for all visible units i, and the number of hidden units is effectively M = 1.

Let us now expand (47) in powers of the weights. The first term on the right hand side of the equation (involving the average over the data distribution) has the same expansion as in the M = 1 case above, see [START_REF] Tsodyks | Associative Memory and Hippocampal Place Cells[END_REF]. For the second term, using

1 + N k=1 cosh(w k1 -w k2 ) cosh(w k1 + w k2 ) = 2 -2 N k=1 w k1 w k2 + O(w 4 ) , (48)
and rescaling the weights, w → w/ √ β, and the learning rate, ν → ν/β, as before, we obtain

1 ν dw i1 dt = w i+1,1 + w i-1,1 -w i1 k w 2 k1 + w 3 i1 + w i1 w 2 2i -w i2 k w k1 w k2 . (49)
Similarly, we find

1 ν dw i2 dt = w i+1,2 + w i-1,2 -w i2 k w 2 k2 + w 3 i2 + w i2 w 2 i1 -w i1 k w k1 w k1 . (50)
The last two terms in the two equations above encode the couplings between the weight vectors attached to the two units. The stationary solutions of these equations are shown in Fig. 14(b). In practice, we observe that the numerical solutions for the weight profiles attached to the two units can have any relative separation between their peaks as long as it is larger than approximately one peak width. The repulsion between the peaks is indeed short range, hence the convergence to the diametrically opposed configuration is really slow. These equations can be turned into two partial differential equation over the continuous space continuous where we resort to the Laplacian to describe spatial diffusion:

1 ν ∂w 1 ∂t (x, t) = ∂ 2 w 1 ∂x 2 (x, t) + 2 -b 1 (t) w 1 (x, t) -c w 2 (x, t) + w 1 (x, t) w 1 (x, t) 2 + w 2 (x, t) 2 , 1 ν ∂w 2 ∂t (x, t) = ∂ 2 w 2 ∂x 2 (x, t) + 2 -b 2 (t) w 2 (x, t) -c w 1 (x, t) + w 2 (x, t) w 1 (x, t) 2 + w 2 (x, t) 2 , (51)
where

b 1 (t) = L 0 w 1 (x, t) 2 dx , b 2 (t) = L 0 w 2 (x, t) 2 dx , c(t) = L 0 w 1 (x, t) w 2 (x, t) dx , (52)
This system describes two diffusing species w 1 and w 2 which are, respectively, self-inhibited by b 1 and b 2 , while coefficient c corresponds to cross-inhibition. The diffusion coefficient for b 1 , b 2 , c tend to infinity and their concentrations are uniform in space. As in the single species (single hidden-unit) case, this dynamical system leads to the stable production of a non trivial pattern over space, corresponding to the emergence of two place fields, see Appendix A.

VI. CONCLUSION

In this work, we have studied the unsupervised learning of simple data distributions, enjoying one or two continuous symmetries, with a RBM. Contrary to standard approaches in machine learning, e.g. convolutional networks, we have not tried to factor out, and hardwire these symmetries in the network. On the contrary, our objective was to see how the symmetries affected the representations of the data and were effectively learned by the machine. This approach is motivated by the fact that most invariances in complex data are actually unknown and it is important to understand how well they can be captured in practice.

In the case of a single hidden (latent) variable, our main observation is that learning is accompanied by a symmetry breaking in the weight space. The hidden unit concentrates only on a small portion of the data configurations; the size of this receptive field is the length over which the variables in the data configurations are correlated. The symmetry is dynamically restored at long times through the diffusion of the receptive field, allowing it to span the whole data manifold. This phenomenon is strongly reminiscent of the concept of continuous attractor (CA) in the context of recurrent neural networks in computational neuroscience [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF][START_REF] Tsodyks | Associative Memory and Hippocampal Place Cells[END_REF][START_REF] Ben-Yishai | Theory of orientation tuning in visual cortex[END_REF], with the major differences that (1) CA usually refer to low-dimensional attractors in the (high-dimensional) space of neural activities, while the CA emerging here defines a low-dimensional manifold in the weight space, and (2) accordingly, the dynamics considered is the learning dynamics acting on weights and not the usual neural dynamics modifying activities. In the case of multiple hidden units, the CAs attached to these units are locked in: weight bumps diffuse coherently along their CA's (Fig. 6), maintaining their relative phases due to mutual repulsive interactions. The resulting multi-unit CA has therefore the same (low) dimension as the underlying symmetry in the data. In practice, however, repulsion is short ranged and may effectively lead to partial decoupling, see Fig. 14(b), and to an increase in the CA dimension. If the number of hidden units is sufficiently large (of the order of the number of visible units over the correlation length) the RBM hidden configurations are effective, coarsegrained version of the data configurations.

An important condition for this CA in the weight space to emerge is that the number of available data exceeds some critical value depending on the configuration size and the intensity of their intrinsic correlations. This phenomenon is a manifestation of the general mechanism of the so-called retarded-learning phase transitions [START_REF] Watkin | Optimal unsupervised learning[END_REF], in which a symmetry-breaking direction (here, in the weight space) is inferred when the ratio of the number of data and of the system size is larger than some signal-dependent (here, the level of intrinsic correlation in the data) threshold. Interestingly, in the presence of multiple invariances, the thresholds associated to these symmetries need not coincide. In such situations, the receptive field will be localized along one dimension in the input space and extended along the other, as seen for the XY model in this work.

It is tempting to make an analogy with recent experimental results on three-dimensional encoding by place cells [START_REF] Treves | [END_REF]. When a rodent explores a set of horizontal (x direction) and (z direction) vertical planes, place cells emerge with place fields localized in either or both planes [START_REF] Casali | Altered neural odometry in the vertical dimension[END_REF]. Yet, if the motion along the x and z axis is not independent, localization can be lost along one of the two directions. For instance, when motion takes place along a helicoidal ramp (x being the angle in the plane perpendicular to the helix axis z), place cells seem to be localized in the angular space and much less so along the vertical axis [START_REF] Hayman | Anisotropic encoding of three-dimensional space by place cells and grid cells[END_REF]. Due to the geometry of the helix, it is reasonable to assume that inputs related to path integration as well as to visual flow are strongly correlated for similar angles (corresponding to a small displacement on the ramp) and much less correlated for small translation along the z-axis, which requires a large physical displacement. It would be interesting to see what happens if the ramp axis is tilted and not vertical any longer. Based on the analogy with the differential retarded-learning transitions, one would expect that place fields are columnar along the ramp axis, and become therefore localized (albeit with different areas) along both x and z directions.

While the analogy with place cells and symmetry-broken hidden units is tempting, establishing a solid connection between our results and neuroscience is far from obvious. Though place cells are known to rely, for their establishment, on various sources of sensory information (including visual and path-integration inputs [START_REF] Chen | How vision and movement combine in the hippocampal place code[END_REF]), the mechanisms underlying the corresponding unsupervised learning processes are far from being elucidated. It is, from this point of view, remarkable that various unsupervised learning rules [START_REF] Brito | Nonlinear Hebbian learning as a unifying principle in receptive field formation[END_REF][START_REF] Sengupta | Manifold-tiling Localized Receptive Fields are Optimal in Similarity-preserving Neural Networks[END_REF][START_REF] Benna | Are place cells just memory cells? Memory compression leads to spatial tuning and history dependence[END_REF] agree with the two main features emerging from the Maximum Likelihood (ML) procedure for RBM studied here, namely the existence of (1) localized receptive field focusing on a subsets of strongly correlated inputs, and of (2) cross-inhibition between hidden units during the learning phase, which makes their place/receptive fields repell each other and forces them to cover as much as possible the input space (Fig. 4). Achieving a more precise understanding of how general this scenario is, and how it extends to deeper architectures i.e with more neural layers would be very interesting.

FIG. 3 .

 3 FIG. 3. (a) Structure of weights learned by a RBM trained weights after 100 epochs from 10,000 data configurations of the one-dimensional Ising model of size N = 100. Training parameters: rate ν = 0.001, batch size S batch = 10, PCD-20 learning. One observes the emergence of a peaked structure in the weights, centered around site i 68. Note the small fluctuations in the tails (small w i1 ), due to the finite (but large) number of data. (b) Width of the peak in the weight space as a function of the correlation length of the Ising model, ξ. We trained our RBM with one hidden unit 25 times on data generated at different temperatures, β, and then calculated the average peak width and the standard deviation (error bars) over the different samples. The width was calculated by fitting a cubic spline with one knot to the profile:y(i) = {w imax(w i )/2}, where i is the site index; The roots of this spline were then determined numerically, and the width was defined as the modulus of the difference between the roots. This procedure reliably finds the Full Width at Half Maximum (FWHM). A linear fit (red line) of the form y = a x + b shows that the width of the place (receptive) field of the only hidden unit is proportional to the relevant characteristic length in the data. Notice that the intercept (b) is non zero, in agreement with the theoretical findings of Section V B in the β → 0 limit.

FIG. 4 .

 4 FIG. 4. (a) Same as Fig. 3(a), but with a RBM having M = 3 hidden units.The weights attached to the same hidden units are shown with the same color. The receptive fields (peaks) for the three hidden units are roughly equally separated from each other. Note that due to the invariance of the probability distributions of the RBM under h µ → -h µ , w iµ → -w iµ (when c µ = 0), the overall sign of the weights attached to the same hidden unit does not matter. (b) & (c) Same as Fig. 3(a), but with a RBM having, respectively M = 25 (a) and M = 100 (c) hidden units. The color codes show the intensity of the weights w iµ as a function of the hidden (µ, y-axis) and visible (i, x-axis) unit labels. The hidden units have been arranged according to the centre of their respective receptive field.

FIG. 5 .

 5 FIG. 5. (a) Trajectories of the peak of the receptive field for a RBM with one hidden unit vs. number of epochs of training. Each trajectory corresponds to a run of the training procedure starting from a random initial condition for the weights. Training parameters: rate ν = 0.1, PCD-20 training and batch size S batch = 100, 1 million configurations of the one-dimensional Ising model at β = 1. Trajectories are corrected for periodic boundary conditions -when the peak crosses the boundary, we add or subtract the size N of the system. (b) Mean Square Displacement (MSD) of the peak of the receptive field vs. number t of epochs of training. Time t = 0 corresponds to the beginning of training. The red line is a linear fit = 2 D t obtained after removing the transitory period of the first 20 epochs of training. Results obtained from the trajectories shown in (a).

Figure 5 (

 5 Figure 5(b) shows that the Mean Square Displacement (MSD) of the peak center grows roughly linearly with the training time (number of epochs), which defines the effective diffusion coefficient of the weight peak. For intermediate numbers of data items, diffusion is activated: due to the inhomogeneities in the empirical data distribution, some places along the ring are preferred, and have a tendency to trap the weight peak for some time.Repeating the same analysis for a RBM with M = 3 hidden units allows us to observe the diffusion of the three peak centers, see Fig.6. We see that the motions of these centers are coupled to maintain a constant distance between each other. This is a clear signature of the effective repulsion between the hidden-unit weight vectors already discussed in Section III B 2.

FIG. 7 .

 7 FIG. 7. (a) Sum of weights W after training a one-hidden-unit RBM on 10,000 configurations of the one-dimensional Ising model at different β. Here, there are N = 100 visible units, i.e. noise ratio r = N/S = 0.01, other RBM training parameters are: ν = 0.001, S batch = 10, trained for 250 epochs with PCD-20. The error bars are calculated over 10 runs of the same training. (b) Sum of weights W after training a one-hidden-unit RBM on configurations of the one-dimensional Ising model at β = 0.5 as a function of the noise ratio r = N/S . Here, there are N = 20 visible units. All other parameters of training are same as (a) . The error bars are calculated over 10 runs of the same training.

  FIG. 8. 100 configurations (index along the y-axis) of the discretized XY model with P = 10 Potts states over N = 100 visible units (x-axis, index i of the units) at inverse temperature β = 1.5. Each colour refers to one of the 10 Potts states v = 0, 1, ..., 9.

FIG. 9 .

 9 FIG. 9. (a). Contour plot of the weights of the RBM with a single real-valued hidden unit with double-well potential trained on the XY model discretized by the Potts model. The x-axis shows space (index i of the visible units), while the y-axis refers to angles (Potts state v). Parameters: P = 10 Potts states, N = 100 visible units, M = 1 hidden unit, trained on 100,000 configurations, learning rate ν = 0.01 and batch size S = 100 trained over 100 epochs. (b). Same as (a) for a binary hidden unit. (c). Angular modulation W angular i of the weight vector fas a function of the space location i, see (16. Same parameters as in Fig. 9(b) for the RBM with real-valued hidden unit. There is a clear strong space localization with a peak centered in unit i = 39. (d). Weight vector w i=39,µ=1,v as a function of the angular-Potts variable v. The line represents the cosine function with frequency = 2π/10 as expected, with the best fit of the phase. (e). Phases ϕ i vs. site index i. Gray dotted line is the phase of the above fitting cosine. (Right:) Frequency ω i vs. site index i. Gray dotted line is again the frequency of 2π/10 of the cosine fit above, which is what one would expect from system with 10 Potts states. See text for the definition of the fitted frequencies and phases. The phase and the frequency is constant across the size of the receptive field, that is all the spins look in the same direction.

FIG. 10 .

 10 FIG. 10. (a) Angular modulation W angular µ,i

FIG. 11 .

 11 FIG. 11. The phenomenon of differentiated retarded learning phase transitions, for various noise ratios r = N/S , where S is the number of training samples. (a) Angular direction: The y-axis is the sum of the components w 1,i,v of the weight vector over all spatial location i = 1, ..., N. The x-axis shows the discretized angular states v = 0, 1, ..., P -1. (b) Spatial direction: The y-axis is the sum of the components w 1,i,v of the weight vector over all angular variables v = 1, ..., P. The x-axis shows the lattice site i = 1, 2, ..., N. Parameters for data generation and sampling: N = 100, P = 10, β = 1.5, M = 1 (Bernoulli hidden unit), ν = 0.1, S batch = 100, trained for 20 epochs.

FIG. 12 .

 12 FIG. 12. (a) Top Eigenvalue of the correlation matrix C of data generated from the one-dimensional Ising model at inverse temperature β = 0.3, as a function of the noise ratio r. The orange straight line is the top eigenvalue Λ I sing , corresponding to a perfectly sampled (infinite S ) Ising model. The blue curve is the top eigenvalue Λ noise of the correlation matrix of the null model with independent variables. (b) Same as panel (a) but for β = 0.5. (c) Top eigenvector of the correlation matrix C (23) of the configurations of the one-dimensional Ising model at fix β but with different numbers S of samples. Ising model samples to calculate the correlation matrix were generated at β = 1 for N = 100 spins.

FIG. 13 .

 13 FIG. 13. (a),(b),(c) Empirical top eigenvector of the correlation matrix C of data generated from the XY model at inverse temperature β = 1.5, for sample size S = 5 10 2 (a), S = 10 4 (b), S = 5 10 6 (c). (d). Corresponding top eigenvalue (dots) as a function of the sample size r = N S . The orange straight line is the top eigenvalue Λ XY (β = 1.5), corresponding to a perfectly sampled (infinite S ) XY model. The blue curve is the top eigenvalue Λ noise of the correlation matrix of the null model with independent variables. The eigenvalues cross r 0.066. (e). Norm of the projection of the empirical top eigenvector ê in the space spanned by the perfect-sampling top eigenvectors e 1 ,e 2 , (ê.e 1 ) 2 + (ê.e 2 ) 2

FIG. 14 . 2 .

 142 FIG. 14. (a)Stationary solution of the small β equations[START_REF] Hayman | Anisotropic encoding of three-dimensional space by place cells and grid cells[END_REF], describing the evolution of the weights w i (t) of a RBM with M = 1 hidden unit trained over many configurations of the one-dimensional Ising model. There are N = 100 visible units. The results shown were obtained with 500 integration steps, starting from small amplitude white noise initial conditions for the weights. (b) Same as (a), but for a RBM with M = 2 hidden units. The two colors shows the weights corresponding to the two units. While the two peaks should be in principle diametrically opposed, i.e. at distance 50 from each other, their mutual repulsion is short ranged; in practice deviations from stationarity smaller than the numerical accuracy cannot be detected. (c) Profiles of the stationary weight vector for a RBM with a unique hidden unit trained on data extracted from the one-dimensional Ising model at small inverse temperature β, see text. The two curves corresponds to the two candidate values for b. The solutions b = b + and b = b -are, respectively, unstable and stable against small fluctuations of the weights.

Acknowledgments. This work was partly funded by the ANR project RBMPro CE30-0021-01 and the HFSP project RGP0057-2016. M.H. benefited from a fellowship from the ICFP Labex of Department of Physics at ENS. J.T. acknowledges partial support by a fellowship from the Edmond J Safra Center for Bioinformatics at Tel Aviv University.

Appendix A: Conditions for pattern formation

The continuous partial differential equation [START_REF] Chen | How vision and movement combine in the hippocampal place code[END_REF] along with (38) describes the evolution of the field w(x, t) in space and time. This equation will only lead to a non-trivial stable steady state pattern if certain conditions that we make explicit below are fulfilled.

We start from (38) and differentiate this equation with respect to time (setting ν = 1 to lighten notations) to get

The non-trivial, uniform fixed point this equation is

One can rewrite the above equations in the following simple notation:

where d is eventually sent on infinity, since there is no spatial time lag for reaching the equilibrium value of b, and b is always spatially uniform. After linearization around the fixed point for small |w|:

these two equations can be written in vector form as follows

are, respectively, the diffusion and stability matrix. We impose first that the uniform fixed point should be stable in the absence of any spatial variation, as we demande that the instability solely comes from spatial interactions. Keeping the non-spatial part of the equation:

We look for solutions of the form w = e λt w 0 . Stability requires that Re(λ) be < 0, that is,

It is easy to check that these general conditions, once applied to the derivatives of f and g listed in (A4), are satisfied as soon as L > 1.

We then ask for the existence of ann instability resulting from the spatial part of the equation. As the Laplacian operator is translation invariant we look for a solution to the reaction-diffusion system (A4) that can be decomposed onto Fourier wave planes of momentum k multiple of 2π/L due to periodic boundary conditions:

where the constants c k are determined by a Fourier expansion of the initial conditions in terms of W k (x) and λ is the eigenvalue that determines the temporal growth of the instability. Inserting (A7) into (A4), we get for each k,

Hence, λ is the root of the following characteristic polynomial:

For the uniform steady state (w * , b * ) to be unstable against spatial fluctuations, we require Re(λ) > 0 for some k 0. The conditions for this can be easily worked out, with the result

and

Condition (A11) is always satisfied since d → +∞. To check (A10) we have to evaluate the coefficient of d, which is f w :

Hence, this condition is satisfied as soon as L > 1.