A Decelle

Giancarlo Fissore
email: giancarlo.fissore@inria.fr

C Furtlehner

Luigi Gresele
email: luigi.gresele@tuebingen.mpg.de

Adrián Javaloy

Bernhard Schölkopf

Aapo Hyvärinen

Spectral Dynamics of Learning Restricted Boltzmann Machines

Keywords: Modélisation générative : physique statistique des Machines de Boltzmann Restreintes, apprentissage avec informations manquantes Machines de Boltzmann Restreintes, Normalizing Flows, Modèles génératifs, Champ moyen, Informations manquantes, Réseaux de neurones Generative modeling : statistical physics of Restricted Boltzmann Machines, learning with missing information and scalable training of Linear Restricted Boltzmann Machines, Normalizing Flows, Generative models, Mean field, Missing information, Neural networks Relative backpropagation . . . . . . . . . . . . . PACS 02.70.Hm -Spectral methods PACS 02.30.Zz -Inverse problems PACS 89.75.-k -Complex systems

Neural network models able to approximate and sample high dimensional probability distributions are known as generative models. In recent years this class of models has

Résumé en français

 qui ont caractérisé le paysage scientifique de ces années-là.

Le premier écart par rapport à ce paradigme apparaît dès les années 1950, lorsqu'une nouvelle approche est proposée pour "procéder sur la base de la conjecture selon laquelle chaque aspect de l'apprentissage ou de toute autre caractéristique de l'intelligence peut en principe être décrit si précisément qu'une machine peut être faite pour le simuler" (McCarthy, 1955). Ce changement de paradigme est une conséquence directe de l'avènement de l'ordinateur moderne et il s'agit d'une approche plus abstraite de l'IA, car elle repose sur l'hypothèse que l'intelligence peut être construite et simulée sur la base de la logique et des systèmes formels. Les premières mises en oeuvre réussies de systèmes d'IA sont dues à cette approche, qui est devenue connue sous le nom d'IA symbolique et a représenté le paradigme dominant pendant des décennies. Cependant, malgré une série de succès initiaux, l'approche de l'IA symbolique a eu du mal à faire face à certains problèmes centraux de la cognition humaine, tels que la perception, l'apprentissage et la reconnaissance des formes. Ces problèmes fondamentaux semblaient être plus facilement abordés par les réseaux neuronaux artificiels, systèmes d'unités de calcul du type du neurone de McCulloch et Pitts, qui ont suscité un renouvellement de l'intérêt pour ce type de modèles dans les années 1980.

C'est à ce moment que l'intérêt pour les réseaux neuronaux a commencé à s'étendre à d'autres domaines, notamment en ce qui concerne la science nouvellement établie de la complexité. Dans ce contexte, des contributions substantielles aux fondements théoriques et au développement des réseaux neuronaux sont venues du domaine de la physique statistique. Parmi celles-ci, Hopfield (1982) a introduit le premier modèle de réseau neuronal fonctionnel pour une mémoire associative. Peu à peu, l'objectif initial de réaliser un cerveau artificiel s'est estompé, et toute une classe de modèles utilisant des réseaux neuronaux artificiels complétés par des outils mathématiques et statistiques est apparue comme un domaine riche en soi. Désigné sous le nom de apprentissage automatique statistique, ce domaine se concentre désormais sur l'ingénierie de systèmes capables d'exploiter les données et les algorithmes de calcul pour résoudre des problèmes complexes.

.

Organisation de la thèse Cette thèse est structurée en deux parties principales. Dans la première partie, nous adoptons le point de vue de la physique statistique sur l'analyse des réseaux de neurones. Cette approche est principalement basée sur la théorie des systèmes désordonnés (Mezard, 1987), un type prototypique de système complexe particulièrement adapté à la description des réseaux de neurones ; trouvant son origine dans les années 1980, cette ligne de recherche est très active de nos jours. Dans ce contexte, l'étude ici présentée concerne la Machine de Boltzmann Restreinte (RBM), Proceeding from geometrical considerations, we propose a stochastic gradient descent optimization algorithm that exploits the matrix structure of fully connected neural networks without imposing any constraints on their structure other than the fixed dimensionality required by invertibility. This algorithm is computationally efficient and can scale to very high dimensional datasets. We demonstrate its effectiveness in training a multilayer nonlinear architecture employing fully connected layers.

Résumé : Les modèles de réseaux neuronaux capables d'approximer et d'échantillonner des distributions de probabilité à haute dimension sont connus sous le nom de modèles génératifs. Ces dernières années, cette classe de modèles a fait l'objet d'une attention particulière en raison de son potentiel à apprendre automatiquement des représentations significatives de la grande quantité de données que nous produisons et consommons quotidiennement. Cette thèse présente des résultats théoriques et algorithmiques relatifs aux modèles génératifs et elle est divisée en deux parties.

Dans la première partie, nous concentrons notre attention sur la Machine de Boltzmann Restreinte (RBM) et sa formulation en physique statistique. Historiquement, la physique statistique a joué un rôle central dans l'étude des fondements théoriques et dans le développement de modèles de réseaux neuronaux. La première implémentation neuronale d'une mémoire associative (Hopfield, 1982) est un travail séminal dans ce contexte. La RBM peut être considérée comme un développement du modèle de Hopfield, et elle est particulièrement intéressante en raison de son rôle à l'avantgarde de la révolution de l'apprentissage profond [START_REF] Hinton | [END_REF]. En exploitant sa formulation de physique statistique, nous dérivons une théorie de champ moyen de la RBM qui nous permet de caractériser à la fois son fonctionnement en tant que modèle génératif et la dynamique de sa procédure d'apprentissage. Cette analyse s'avère utile pour dériver une stratégie d'imputation robuste de type champ moyen qui permet d'utiliser la RBM pour apprendre des distributions empiriques dans le cas difficile où l'ensemble de données à modéliser n'est que partiellement observé et présente des pourcentages élevés d'informations manquantes.

Dans la deuxième partie, nous considérons une classe de modèles génératifs connus sous le nom de Normalizing Flows (NF), dont la caractéristique distinctive est la capacité de modéliser des distributions complexes à haute dimension en employant des transformations inversibles d'une distribution simple et traitable. L'inversibilité de la transformation permet d'exprimer la densité de probabilité par un changement de variables dont l'optimisation par Maximum de Vraisemblance (ML) est assez simple mais coûteuse en calcul. La pratique courante est d'imposer des contraintes architecturales sur la classe de transformations utilisées pour les NF, afin de rendre l'optimisation par ML efficace. En partant de considérations géométriques, nous proposons un algorithme d'optimisation stochastique par descente de gradient qui exploite la structure matricielle des réseaux de neurones entièrement connectés sans imposer de contraintes sur leur structure autre que la dimensionnalité fixe requise par l'inversibilité. Cet algorithme est efficace en termes de calcul et peut s'adapter à des ensembles de données de très haute dimension. Nous démontrons son efficacité dans l'apprentissage d'une architecture non linéaire multicouche utilisant des couches entièrement connectées. received tremendous attention due to their potential in automatically learning meaningful representations of the vast amount of data that we produce and consume daily. This thesis presents theoretical and algorithmic results pertaining to generative models and it is divided in two parts.

In the first part, we focus our attention on the Restricted Boltzmann Machine (RBM) and its statistical physics formulation. Historically, statistical physics has played a central role in studying the theoretical foundations and providing inspiration for neural network models. The first neural implementation of an associative memory (Hopfield, 1982) is a seminal work in this context. The RBM can be regarded to as a development of the Hopfield model, and it is of particular interest due to its role at the forefront of the deep learning revolution [START_REF] Hinton | [END_REF]). Exploiting its statistical physics formulation, we derive a mean-field theory of the RBM that allows us to characterize both its functioning as a generative model and the dynamics of its training procedure. This analysis proves useful in deriving a robust mean-field imputation strategy that makes it possible to use the RBM to learn empirical distributions in the challenging case in which the dataset to model is only partially observed and presents high percentages of missing information.

In the second part we consider a class of generative models known as Normalizing Flows (NF), whose distinguishing feature is the ability to model complex high-dimensional distributions by employing invertible transformations of a simple tractable distribution. The invertibility of the transformation allows expressing the probability density through a change of variables, whose optimization by Maximum Likelihood (ML) is rather straightforward but computationally expensive. The common practice is to impose architectural constraints on the class of transformations used for NF, in order to make the ML optimization efficient. Proceeding from geometrical considerations, we propose a stochastic gradient descent optimization algorithm that exploits the matrix structure of fully connected neural networks without imposing any constraints on their structure other than the fixed dimensionality required by invertibility. This algorithm is computationally efficient and can scale to very high dimensional datasets. We demonstrate its effectiveness in training a multilayer nonlinear architecture employing fully connected layers.

un réseau neuronal simple qui a été parmi les protagonistes de la phase initiale de la révolution de l'apprentissage profond (Hinton, 2006). Ce réseau est particulièrement intéressant du point de vue de la physique statistique car il peut être interprété comme une généralisation du modèle de Hopfield [START_REF] Barra | On the equivalence of Hopfield networks and Boltzmann Machines[END_REF]. Dans le Chapitre 1, nous présentons le modèle de Hopfield et discutons des principes de base permettant de dériver une théorie de champ moyen pour les réseaux neuronaux. Dans le Chapitre 2, nous présentons la RBM. Nous commençons à discuter de nos travaux originaux dans le Chapitre 3, avec la présentation d'une description empirique des propriétés dynamiques du RBM. Cela pose la base de l'analyse théorique de champ moyen développée dans le Chapitre 4, en discutant des propriétés d'équilibre et dynamiques. Enfin, dans le Chapitre 5, nous montrons comment adapter la RBM au scénario difficile de l'apprentissage avec des informations manquantes, en tirant parti de l'image théorique présentée dans les sections précédentes.

Dans la deuxième partie, nous nous intéressons aux Flux Normalisateurs (NF) [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF], un type de modèle plus récent qui trouve des applications dans des domaines similaires à ceux de la RBM. Nous nous éloignons ici de la physique statistique : dans le Chapitre 6, nous proposons un modèle NF qui utilise des réseaux neuronaux entièrement connectés et nous montrons comment, à partir de simples considérations de géométrie riemannienne, nous pouvons dériver un algorithme efficace pour entraîner le modèle proposé.

Generative modeling: statistical physics of Restricted Boltzmann Machines, learning with missing information and scalable training of Linear Flows Giancarlo Fissore June 29, 2022 Abstract Neural network models able to approximate and sample high dimensional probability distributions are known as generative models. In recent years this class of models has received tremendous attention due to their potential in automatically learning meaningful representations of the vast amount of data that we produce and consume daily. This thesis presents theoretical and algorithmic results pertaining to generative models and it is divided in two parts.

In the first part, we focus our attention on the Restricted Boltzmann Machine (RBM) and its statistical physics formulation. Historically, statistical physics has played a central role in studying the theoretical foundations and providing inspiration for neural network models. The first neural implementation of an associative memory (Hopfield, 1982) is a seminal work in this context. The RBM can be regarded to as a development of the Hopfield model, and it is of particular interest due to its role at the forefront of the deep learning revolution [START_REF] Hinton | [END_REF]). Exploiting its statistical physics formulation, we derive a mean-field theory of the RBM that allows us to characterize both its functioning as a generative model and the dynamics of its training procedure. This analysis proves useful in deriving a robust mean-field imputation strategy that makes it possible to use the RBM to learn empirical distributions in the challenging case in which the dataset to model is only partially observed and presents high percentages of missing information.

In the second part we consider a class of generative models known as Normalizing Flows (NF), whose distinguishing feature is the ability to model complex high-dimensional distributions by employing invertible transformations of a simple tractable distribution. The invertibility of the transformation allows expressing the probability density through a change of variables, whose optimization by Maximum Likelihood (ML) is rather straightforward but computationally expensive. The common practice is to impose architectural constraints on the class of transformations used for NF, in order to make the ML optimization efficient.

Introduction A brief history of Artificial Intelligence

The field of artificial intelligence (AI) has changed and evolved considerably over the last decades. The same term "artificial intelligence" has taken on various different meanings following the advancements of the epoch, being more often then not only loosely defined. We identify the genesis of the field with the first attempt to model biological neurons mathematically, proposed by Mc-Culloch and Pitts [START_REF] Mcculloch | A Logical Calculus of Ideas Immanent in Nervous Activity[END_REF] in the 1940s. Their approach took inspiration from the biological brain and it is framed in a set of efforts to build an artificial brain capable of "human" intelligence, fueled by novel advances and the emergence of modern scientific disciplines such as the theory of computation, neuroscience and cybernetics [START_REF] Waldrop | The dream machine: J.C.R. Licklider and the revolution that made computing personal[END_REF][START_REF] Wiener | Cybernetics: or Control and Communication in the Animal and the Machine[END_REF] that have characterized the scientific landscape in those years. The first departure from this paradigm comes as early as in the 1950s, when a new approach is proposed to "proceed on the basis of the conjecture that every aspect of learning or any other feature of intelligence can in principle be so precisely described that a machine can be made to simulate it" ( [START_REF] Mccarthy | A PROPOSAL FOR THE DARTMOUTH SUMMER RESEARCH PROJECT ON ARTIFICIAL INTELLIGENCE[END_REF]). This paradigm shift is a direct consequence of the advent of the modern computer and it is a more abstract approach to AI, as it is based on the assumption that intelligence can be constructed and simulated on the basis of logic and formal systems. The first successful implementation of AI systems are due to this approach [START_REF] Newell | The logic theory machine-A complex information processing system[END_REF], which became known as Symbolic AI and represented the dominant paradigm for decades. However, in spite of a series of initial successes, the Symbolic AI approach struggled to cope with certain problems that are central to human cognition, such as perception, learning and pattern recognition. These fundamental problems appeared to be more easily approached by artificial neural networks [START_REF] Rumelhart | Learning Representations by Back-propagating Errors[END_REF], systems of computational units of the kind of the McCulloch and Pitts neuron, which sparked a renewed interest in this kind of models in the 1980s. It is at this point that interest in neural networks started to spread to other fields, especially with regard to the newly established science of complexity [START_REF] Waldrop | Complexity: the emerging science at the edge of order and chaos[END_REF]. In this context, substantial contributions to the theoretical foundations and the development of neural networks came from the field of statistical physics. Among these, [42] introduced the first functioning neural network model for an associative memory, subsequently analyzed in exquisite detail in [5,7,6], and [30,29] produced seminal work on the storage capacity of neural networks. Gradually, the original goal of realizing an artificial brain has v faded away, and a whole class of models employing artificial neural networks complemented with mathematical and statistical tools has emerged as a rich field in its own right. Referred to as statistical machine learning, its focus has become that of engineering systems capable of leveraging data and computations to solve complex problems. In the 2010s, the technological advancement has determined the deep learning revolution: the availability of more powerful hardware and network infrastructures for computation, as well as access to the massive amount of data constantly generated by myriads of interconnected technological devices, have made it possible to build increasingly larger and more powerful neural network models leveraging larger amounts of data to solve the most disparate problems with increasing success. The relation to AI is to be found in the fact that deep learning models have been employed to automatically solve problems that traditionally needed the employment of human intelligence to be solved, like object recognition [48], speech synthesis [START_REF] Van Den Oord | Parallel WaveNet: Fast High-Fidelity Speech Synthesis[END_REF] and natural language translation [23].

Thesis organization

This thesis is structured in two main parts. In the first part, we adopt the statistical physics point of view on the analysis of neural networks. This approach is mainly based on the theory of disordered systems [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF], a prototypical kind of complex system that is particularly suited to describe neural networks; finding its origin in the 1980s, this line of research is very active nowadays [17]. In this context, the investigation here presented concerns the Restricted Boltzmann Machine (RBM), a simple neural network that has been among the protagonists in the initial phase of the deep learning revolution [39]. This is particularly interesting from the statistical physics point of view as it can be interpreted as a generalization of the Hopfield model [9]. In Chapter 1 we introduce the Hopfield model and discuss the basic principles to derive a mean-field theory for neural networks. In Chapter 2 we introduce the RBM. We start discussing original work in Section 3, with the presentation of an empirical description of the dynamical properties of the RBM. This poses the basis for the theoretical mean-field analysis developed in Chapter 4, discussing equilibrium and dynamical properties. Finally, in Chapter 5 we show how to adapt the RBM to the challenging scenario of learning with missing information, taking advantage of the theoretical picture discussed in previous Sections.

In the second part, we turn our attention to Normalizing Flows (NF) [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF], a more recent kind of models that finds applications in similar domains as the RBM. Here we depart from the statistical physics focus: in Section 6.2 we propose a NF model that employs fully connected neural networks and in Sections 6.3 and 6.4 we show how proceeding from simple Riemannian geometry considerations we can derive an efficient algorithm to train the proposed model. vi

Contributions

Contributions

A reading guide Let's start with a stylistic note. Throughout the thesis, I use "we" as the subject. This is to signal that all the work discussed has been done with the help and guidance of my supervisors and collaborators. This is the only Section in which I use "I", and this is to point out in a precise way my personal contributions.

In Part I, Chapter 1 serves as an introduction to the tools and models that we use in the following Chapters. While the topics here discussed are wellknown, the presentation is rather original. The whole treatment is meant to be self-contained, ultimately introducing the RBM in a principled way. Chapter 3 reports the empirical study of the RBM that I performed in the initial phase of my PhD. In the development of the mean-field theory of Chapter 4 I took on a more marginal role. The main theoretical analysis has been developed by my supervisors while I contributed in general discussions, in performing numerical analysis and writing the paper. This Chapter is supposed to be reasonably selfcontained, and it represents my attempt to more gently introduce the related paper (Reprint B). For the work in Chapter 5 and Part II, I played a more central role both in the design of the proposed algorithms and their development. These Chapters are not self-contained; they represent a brief introduction to the associated papers (Reprints C and D), that I suggest reading in full.

Publications Parts I and II of this thesis are meant to introduce and complement the papers below, whose content is reported in full in the Reprints section. To investigate the mechanisms underlying intelligent behavior, we need a model of the brain that is simple enough to be studied yet rich enough to be functional. More precisely, we want to come up with models and algorithms that are able to realize complex tasks like pattern formation, object recognition and associative memories. A fundamental assumption in this context is that solutions to these complex tasks emerge from the interaction dynamics of a large number of simple elementary constituents, so we can abstract from the biological implementation details and build models that are useful and that we can analyze thoroughly.

The first successful attempt to model neurons mathematically has been proposed by McCulloch and Pitts [START_REF] Mcculloch | A Logical Calculus of Ideas Immanent in Nervous Activity[END_REF]. Their model of the brain consists of a network of "on/off" binary nodes n i with associated thresholds µ i and connected by weights w ij (connecting nodes i and j). The behavior of the network is determined by a simple time-dependent update rule:

n i (t + 1) = Θ   j w ij n j (t) -µ i   . (1.1)
This update rule models the firing mechanisms observed in networks of real neurons. Here the nodes are modeled as individual processing units that activate when the inputs from the other units is above a certain threshold; the weights w ij mimic the role of synaptic connections, with positive and negative values respectively modeling excitatory and inhibitory synapses, and the threshold values µ i account for the need of an activation potential. Many biological details of real neurons are not directly included in this model, and we can assume that their role is not fundamental in determining the emergence of complex behavior in the network. Nonetheless, some features of the model are not biologically plausible, most notably the implicit need for a global synchronization mechanism (note that activations at time t + 1 depend on the state of the full network at time t) which has not been observed in the animal brain.

Background

We will see in Section 1.2 how in the Hopfield model we relax this constraint to obtain a functioning associative memory; here let us stress the fact that ultimately our objective is not to build a model of the brain which is as close to reality as possible but to study the abstract mechanisms that give rise to complex behaviors, and it is in this sense that the McCulloch-Pitts network of neurons is a success.

In the following Section we will introduce the Hopfield network, a simple generalization of the McCulloch-Pitts model that gives rise to an associative memory. Subsequently we will introduce the tools that let us analyze the behavior and the limits of such networks, and we will see how to apply them to the Hopfield model. Finally, we will introduce the main subject of interest of this thesis, the Restricted Boltzmann Machine (RBM).

The Hopfield model

One of the simplest tasks that we can realize with a neural network is that of implementing an associative memory, i.e. a memory system that is contentaddressable: given a partial observation of a pattern as input, the memory retrieves a previously stored pattern that is similar to the input. A popular neural network system to implement an associative memory is the Hopfield model [42], consisting in a network of N binary nodes s i = ±1 with associated thresholds θ i and connected by weights w ij .

Storage of patterns.

A pattern is represented by a specific configuration of the nodes s i and storage of patterns in the network is achieved through Hebbian learning [36], which originates from the observation that synaptic connections are stronger among neurons that activate simultaneously. The Hebbian rule is commonly stated as "neurons that fire together, wire together" and in the context of the Hopfield model with p binary patterns x µ to store it takes the form

w ij = 1 N p µ=1
x µ i x µ j , i = j.

(1.2)

We note that the model can be described by a symmetric weight matrix W ; furthermore, connections of a node with itself are not allowed so that we have w ii = 0.

Patterns retrieval. In contrast to a network of McCulloch-Pitts neurons, the Hopfield network is updated asynchronously: at each time step we select a node s i at random and we update its value following the rule

s i = sign   j w ij s j + θ i   , sign(x) = -1 x ≤ 0 +1 x > 0 (1.3)
We assume without loss of generality that θ i = 0. Further assuming that the number of nodes in the network is large enough, stored patterns are stable 1.2. The Hopfield model under the application of the above update rule. To show this, let's call h ν i the input to a node i when the network is presented with a stored pattern x ν

h ν i = j w ij x ν j = 1 N j µ x µ i x µ j x ν j .
If the number of stored patterns is small w.r.t. the number of nodes and recalling that the patterns are binary, we have

h ν i = x ν i + 1 N j µ =ν x µ i x µ j x ν j
x ν i and thus the following stability condition holds sign(h ν i ) = x ν i . (1.4) This let us propose a possible strategy for retrieval of stored patterns. Starting from an initial configuration, we let the network evolve under the update rule (1.3); if we start from a configuration that is close to a stored pattern x ν , we might expect the network to stabilize on the configuration corresponding to x ν , thus realizing the objective of a content-addressable memory. We now go on by presenting a simple analysis of the convergence properties of the proposed strategy, to show that it is indeed a sensible alternative to realize an associative memory.

Convergence. A simple argument for the convergence of a Hopfield network to a stable state is provided in [13]. We start by representing the network as a bipartite graph; one partition includes the nodes whose value is +1 and the other partition includes the nodes whose value is -1. Each node will have "internal" connections linking it to the nodes in the same partition and "external" connections linking it to the nodes in the other partition. At each step we select a node and we compute the sum of internal and external connections; if the sum of the external connections is higher, we flip the node's value and we move it to the opposite partition. This process is equivalent to adopting the Hopfield update rule (1.3) but from this point of view the convergence of the network amounts to finding a minimum cut of the graph representing the network. Given that at every step the number of connections linking the two partitions either decreases or stays the same and that the number of possible states of the network is finite (so we would eventually explore all nodes with high probability) we are guaranteed to end up in a stable state.

Unfortunately, what is not guaranteed in the Hopfield model is that a stable state corresponds to a stored pattern. It is easy to see, for instance, that the reversed patterns -x ν are also stable. As another example of a stable state, we consider the special case in which the patterns are random and we take the mixed state x m defined as x m i = sign(x µ1 i + x µ2 i + x µ2 i ).

Background

As the patterns are random, 3 times out of 4 we expect to observe x m i = x µ1 i (and the same is valid for x µ2 i and x µ3 i ), so that for N large the following identity holds

j x µ1 j x m j = 3 4 N - N 4 = N 2 .
Assuming again that the number of patterns is small (w.r.t. the number of nodes) we get

h m i = j w ij x m j = 1 2 x µ1 i + 1 2 x µ2 i + 1 2 x µ3 i + 1 N µ k ∈(µ1,µ2,µ3) j µ =µ k x µ i x µ j x m i 1 2 x µ1 i + 1 2 x µ2 i + 1 2 x µ3 i
which is equivalent to the stability condition (1.4) and thus x m is a stable configuration of the network even though it is not a stored pattern. In general, the argument can be extended to show that linear combinations of an odd number of patterns are stable, and these are called mixture states.

As a final remark, let us note that patterns stability has been assessed under the assumption that the number of patterns is small compared to the number of nodes in the network. In principle this is not very problematic as we can always assume the network to be large enough, but from a practical perspective it would be useful to know what is the capacity of the model, i.e. how many patterns we can store as a function of the size N of the network. We will deal with this matter and other details about the conditions to guarantee the functioning of the Hopfield model as an associative network in Section 1.7.

The energy function. To each configuration s of a Hopfield network we can associate a scalar function E defined as follows

E(s) = - 1 2 N i=1 N j=1 w ij s i s j + N i=1 θ i s i . (1.5)
This is called the energy function of the model, for reasons that will be clarified in Section 1.3. For now we are interested in showing that the dynamics induced by the update rule (1.3) minimizes this quantity. Starting from a configuration s, at each step we select a node k and either we leave its value unchanged or we transition to a new configuration s with s k = -s k . In the first case the energy doesn't change, while in the second case we have

∆E = E(s ) -E(s) = (s k -s k )   j w kj s j -θ k   1.3. Statistical description of a neural network
which is necessarily negative to be consistent with the update rule (1.3). Similarly, when the network attains a stable configuration we have ∆E > 0 if we flip the value of any node.

We can now paint an interesting picture of the Hopfield model, that will be the starting point for the more advanced analysis presented later on. The configurations of the network live in an energy landscape with many minima, represented by stable configurations of the system. The network evolves under rule (1.3) to minimize the energy of the system, attaining the stable configurations. As the energy decreases at each step of the dynamical evolution, we can think of the stable states as attractors of the system; initializing the system close to a stable configuration, the network will be attracted to it. In this sense it seems reasonable that under the right conditions the system can function as an associative memory; Section 1.7 is devoted to delineating these conditions.

Statistical description of a neural network

We turn to a stochastic description of a network of computational neurons, which amounts to introduce some uncertainty in the model by associating probabilities to the possible states of the system. This let us reason directly on the expected probabilistic behavior of the model and not on specific instances of it, allowing us to reach general conclusions. Moreover, by analogy to neuronal networks it seems reasonable to introduce some noise as a way to take into account minor effects due to the complexities of the actual biological systems.

As a neural network of choice, we take the Hopfield model introduced in Section 1.2 and define a probability distribution P (s) over the nodes of the network. On this basis, we derive the exact form of this probability distribution and we introduce the necessary tools to analyze the expected behavior of the system in great details. We remark that the Hopfield model as described in Section 1.2 can already be considered weakly stochastic; indeed, its dynamics depend on the random choice of nodes, but the update rule (1.3) is deterministic.

Defining a probability distribution over the nodes s we allow the equilibrium configurations to fluctuate around their expected value. In this context the description of the model in terms of the energy function is useful, as allowing the configurations to fluctuate is equivalent to allowing the energy to fluctuate. But what is a good functional form for the probability distribution over the nodes? Our strategy is to look for a distribution that is as general as possible, i.e. we want to avoid making unnecessary assumptions. This can be formulated as an instance of the Occam's razor, a principle stating that the most economical theories are to be preferred. The logical and quantitative reasons motivating the success of this principle lie in the fact that from a Bayesian standpoint the most economical explanation consistent with the observations at hand is more probable ( [55], Chapter 28). Intuitively, in our case the only observation that the appropriate probability distribution should encode is that the average value of the energy is fixed and constant; as for the rest, avoiding assumptions 1. Background means that we are left in a state of maximum uncertainty about the fluctuations of the system around the average. The information-theoretic expression that quantifies uncertainty is the entropy of the distribution [START_REF] Shannon | A mathematical theory of communication[END_REF]18] H[P ] =s P (s) log(P (s)).

A sensible strategy to derive the functional form of the probability P (s) is thus to apply the principle of maximum entropy [45]: we look for the normalized distribution whose entropy is maximal under the only constraint of constant average energy maximize

H[P ] subject to s P (s) = 1, s P (s)E(s) = E .
We can solve this constrained optimization problem by using Lagrange multipliers. Introducing the multipliers λ and β we can write the Lagrangian

L = H(s) -λ s P (s) -1 -β s P (s)E(s) -E
that we extremize (eliminating λ) to get

P (s) = e -βE(s) Z , Z = s e -βE(s) . (1.6)
Here β is a parameter that regulates the noise. To see this, we discuss the stochastic dynamics of the network. An equilibrium configuration s of the network is determined by a certain energy E(s) with associated probability P (s). Being the system at equilibrium, the total probability of the network to transition to a different state s must be counterbalanced by the probabilities of transitioning to state s when found in a state s = s. Calling T s→s the transition probability from state s to state s , a sufficient condition to guarantee equilibrium is given by the detailed balance condition

T s→s P (s) = T s →s P (s )
which we can write as

T s→s T s →s = P (s ) P (s) = e -β(E(s)-E(s )) = e -β∆E . (1.7)
A stochastic update rule for the Hopfield network must then obey condition (1.7). A particularly fitting choice is the update rule defined by Glauber dynamics [32]

s i = +1 with probability f β (h i ) -1 with probability 1 -f β (h i ) (1.8)
1.4. Free energy where

f β (h i ) = 1 1 + e -2βhi h i = j w ij s j + θ i .
Given rule (1.8) we see that the evolution of the system is completely random for β → 0 (maximum noise) while it falls back to the deterministic rule (1.3) for β → +∞ (no noise).

Equation (1.6) defines the Boltzmann-Gibbs distribution, with Z being called the partition function. Such a distribution is used to model a wide variety of physical systems at thermal equilibrium, where the energy function effectively represents the energy of the system. In this context, the temperature is a proxy for the noise and it is related to β as

β = 1 T .
We will refer to β as the inverse temperature and to T as the temperature of our model. While it might seem strange to define a temperature for a neural network, we note that it is in truth pretty natural. In physical systems the temperature is generally interpreted as a quantity that specifies the level of noise, abstracting from the actual microscopic mechanisms generating the noise; here, we introduced the inverse temperature in the exact same way.

Free energy

We have seen in Section 1.2 that in the deterministic case the energy of the Hopfield network is minimized; the model evolves towards the states of minimum energy and such states are stable. In the stochastic case, due to the noise, the system is allowed to move away from the states of lowest energy. The actual value of the energy fluctuates, while its average value is constant. The quantity that is now extremized is the entropy, which is required to be maximal. Rewriting the entropy to explicitate the average energy

H[P ] = - s P (s) log P (s) = β E + log Z
we see that the quantity that gets minimized is now log Z, which is defined to be the free energy

F def = -T log Z = E(s) -T H[P ]. (1.9)
This let us refine our picture of the Hopfield model. The system still evolves towards states of minimal energy, but this tendency is counteracted by the entropic term that allows other energy configurations at equilibrium. The entropic term is weighted by the temperature: when T → 0 we get back to the 1. Background deterministic case and the energy of the system is minimized, while for T → +∞ the entropic term dominates the dynamics and all energy configurations are allowed.

The free energy fully determines the behavior of the system. From its expression we can derive a variety of quantities of interest to describe the system, like the average values of the nodes and their correlations

s i = - ∂F ∂θ i s i s j = - ∂F ∂w ij .
A tractable form of the free energy is thus a very powerful tool to analyze the system in depth. Unfortunately, this is not easy to compute as the partition function Z involves a sum over an exponential number of terms. A standard procedure is to express the free energy in a more general form and compute its value in an approximated way, as we'll see in Section 1.6.

Order parameters and phase transitions

The Hopfield network can operate in different regimes, depending on the aggregated characteristic of the actual configurations of the nodes. These regimes are called phases, and the terminology for the regimes comes from the theory of the Ising model. The Ising model is an abstract model of a ferromagnet, in which the nodes are called spins and the spins model the magnetic moments of atoms in magnetic material. The specificity of the Ising model is that the interactions among spins (the weights, in the Hopfield model) are all positive and constant. At high temperature the entropic term dominates and the spins are randomly oriented; the average value of the spins is thus zero, and this phase is called paramagnetic in analogy to the paramagnetic material. When the temperature is low, the interaction among the spins drive the system to the minimum of the energy, which is identified by the configuration in which all the spins are aligned. The system is thus found in a different phase in which the spins are ordered, called the ferromagnetic phase; the magnetization is now different from 0 and the transition from a null value to a finite value signals the transition to a different phase of the system. The magnetization is called order parameter as it signals the phase transition. For the Hopfield model, at high temperature we can observe the same behavior and we have again a paramagnetic phase. The behavior in the ferromagnetic phase is more subtle; given that the weights can have both positive and negative values, the system can be found in situations in which there are no ways to coordinate all the nodes to attain orientations that minimize the energy. These are frustrated configurations that evidence how the magnetization is not a good order parameter to describe the Hopfield model, and more refined alternatives need to be used. Edwards-Anderson order parameters are appropriate to describe the separation between paramagnetic and spin glass phase.

1.6. Effective free energy

Effective free energy

In this section we approach the problem of computing a good approximation to the free energy by exploiting the saddle point method, a general methodology that is applicable to a wide variety of models. We recall that the exact form of the free energy (1.9) contains an exponential number of terms, making the computation of its value a computationally hard problem. Adopting a statistical description of the Hopfield model, we know from Section 1.4 that the system evolves towards equilibrium states that minimize the free energy. Furthermore, the equilibrium configurations are found in specific phases of the system, characterized by the appropriate order parameters. Exploiting these observations, the strategy consists in defining an effective free energy as a function of the order parameters. By construction, the minimum of the effective free energy will be equivalent to the true free energy of the system, reducing the problem of computing the free energy to a minimization problem. Further assuming that the system is "large enough" we can compute the minimization with the saddle point approximation. To exemplify this methodology, we will consider the Ising ferromagnet and the associated magnetization. This choice is made to simplify the presentation of the method, as the Hopfield model necessitates of more elaborate order parameters (Section 1.5) and some more involved mathematical manipulations. The general approach is nonetheless the same; details about the Hopfield model are discussed in Section 1.7 and a more elaborate method is presented in Section 4.2 with the introduction of the Replica Method.

For the Ising ferromagnet, the magnetization m is the only order parameter that we need to characterize the system

m = 1 N i s i . (1.10)
At equilibrium, the system will have a specific magnetization m eq for which the free energy is defined. To extend the notion of free energy, we want to consider all possible values of magnetization. We start by rewriting the partition function in terms of the free energy

Z = s e -βE(s) = e -βF .
Now we can fix the value of the magnetization and consider all configurations s m of the system that are consistent with such value, i.e. for which we have

1 N i s i = m.
Summing over all possible values of the magnetization, we can rewrite the partition function as (sm) from which a natural generalization of the free energy for a specific magneti-1. Background zation is defined .11) We call this magnetization-dependent quantity the effective free energy and we will also be interested in its per-node value

Z = m sm e -βE
F (m) def = -T log Sm e -βE(Sm) , Z = m e -βF (m) . ( 1
f (m) = F (m) N .
For a large system (N → +∞) the magnetization can take an almost-continuous set of values (flipping one node determines a change of 2 N in the magnetization) so we can rewrite the partition function as an integral

Z = N 2 +1 -1 dm e -βN f (m) .
As it turns out, for large N the value of the above integral is dominated by the minimum value of f . We can thus substitute the integral to obtain

Z e -βN f (mmin) .
This is called the saddle point approximation and it lets us link the effective free energy to the equilibrium free energy

F F (m min ), m min ≡ m eq .
The problem of computing the partition function can be solved by deriving a tractable form for f (m) and determining m eq .

Mean-field theory of the Hopfield model

A complete mean-field theory of the Hopfield model has been derived in [5,7] under the assumption that the variables describing the patterns are independent. In this section we summarize the derivation and clarify the role of the patterns in defining a good order parameter for the Hopfield model. To begin with, we will follow the strategy presented in Section 1.6 and derive a tractable effective free energy. Subsequently, we discuss the need to average over the patterns to describe a general statistical behavior of the system. We will only present the results of this averaging; a more detailed account can be found in [37,[START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF].

Effective free energy. To better understand why the partition function is hard to compute, we rewrite the energy function of the Hopfield model as

E(s) = - 1 2N p µ=1 i s i ξ µ i 2 + const.
where we made explicit the dependence of the weights w ij on the patterns and we neglected the thresholds θ i . The constant term comes from the diagonal weights w ii and we can drop it as it doesn't impact the energy minimization. We observe that the energy function can be interpreted as a sum of contributions from each stored pattern, where each contribution is measured by the square of the overlap of the system configuration with the considered pattern.

The partition function takes the form

Z = S e βE(S) = S exp   - β 2N p µ=1 i s i ξ µ i 2  
which is hard to compute because the sum in the exponential includes quadratic terms that cannot be factored. A simple yet powerful strategy to deal with this problem is to introduce a Gaussian integral to verify the identity exp

  - β 2N i s i ξ µ i 2   = βN 2π p +∞ -∞ dm µ × exp - βN 2 (m µ ) 2 + i s i ξ µ i m µ
where we have introduced an auxiliary variable m µ for each pattern. This is a Hubbard-Stratonovich transformation, useful to introduce the order parameters m = {m µ , µ = 1, . . . , p}. We will give the interpretation for these order parameters in the next paragraph; for now we proceed to rewrite the partition function by exploiting the above identity

Z = s βN 2π p +∞ -∞ µ dm µ exp - βN 2 (m µ ) 2 + β i s i ξ µ i m µ .
(1.12) The overlap terms in the exponential are now linear, so we can introduce a tractable effective free energy f (β, m)

Z = βN 2π p +∞ -∞ dm e -βN f (β,m) f (β, m) = 1 2 m 2 - 1 βN i log 2 cosh β µ m µ ξ µ i . (1.13)
Employing the saddle point approximation we calculate the equilibrium free energy

F = N f (β, m eq )
and we can derive an implicit expression for m µ eq ∂f ∂m µ

m µ eq = 0 =⇒ m µ eq = 1 N i ξ µ i tanh β ν m ν eq ξ ν i .

Background

Order parameters. In the previous paragraph we have omitted the thresholds θ i to simplify the notation. Reintroducing them in a Hebbian fashion, we get a set of external fields θ i = µ h µ ξ µ i correlated to the patterns, where h µ is a scalar value determining the strength of the threshold. The partition function and the effective free energy (1.13) become

log Z = log S exp   β 2N p µ=1 i s i ξ µ i 2 + h µ i ξ µ i s i   f (β, m) = 1 2 m 2 - 1 βN i log 2 cosh β µ (m µ + h µ ) ξ µ i .
Deriving both expressions by the thresholds' strenghts h µ we get

∂F ∂h µ = - 1 β ∂ log Z ∂h µ = - i s i ξ µ i ∂F ∂h µ = N ∂f ∂h µ = - i ξ µ i tanh β ν (m ν + h ν ) ξ ν i
and similarly to the previous paragraph we can get the equilibrium value of the magnetization

m µ eq = 1 N i ξ µ i tanh β ν m ν eq + h ν ξ ν i .
Putting together all of the above equations we finally get the equivalence

m µ = 1 N i s i ξ µ i (1.14)
where we see that m µ represents the overlap between the network's configuration and pattern ξ µ . Its value depends on the correlation of the nodes s i to the pattern: we have m µ = 1 when the configuration of the network reproduces pattern ξ µ exactly and m µ = 0 when there is no correlation. The order parameters m µ can thus be thought of as the magnetizations of the network correlated to the µth pattern, providing a concrete physical interpretation.

Patterns averaging and phase diagram. The expression that we have derived for the effective free energy depends explicitly on the specific patterns to store. To describe our model in full generality, we want to get rid of the patterns by averaging out their effect on the free energy. This averaging process requires fixing the statistical properties of the patterns that we can store, to come up with a description of a statistical ensemble rather than a specific model instance. The way this ensemble is chosen is fundamental to describe functional systems, and the assumptions made in this context define the range of systems that are well described by the mean-field picture. We stress this point here as it is fundamental in our treatment of a realistic model of RBM.

For the Hopfield model, the assumption is that the variables representing the patterns are i.i.d. with ξ µ i = ±1 at random. Here, the difficulty lies in the fact that it is not clear how to perform the average over the patterns in Equation (1.13). In Section 4.2 we will introduce the Replica Method to address this problem; for now, we skip the details on the averaging to list the results for the Hopfield model.

In Section 1.2 we have seen how the stored patterns are stable states of the network. We have also seen how certain kinds of mixture states are stable, and we couldn't roll out the possibility of stable spurious states, i.e. states that are not correlated to the patterns. Moreover, we worked with the assumption that the number of stored patterns is small and we didn't provide any indications about the capacity of the network. By analyzing the stability of the minima of the averaged mean-field effective free energy, we can address all of the above problems. A detailed discussion is provided in [37], and the results are summarized in the phase diagram reported in Figure 1.1. We define the capacity α of the network as the number of patterns per unit node: α = p N . The patterns are stable only in region A and B of the phase diagram; in both regions we also find stable mixture and spurious states, with the mixture states always presenting higher free energy than the stored patterns and thus the mixture 1. Background states are local minima in regions A and B. In region A the patterns represent global minima and the spurious states are local minima, while in region B the situation is reversed. Patterns are thus stable in the A region and metastable in the B region. We call the AB region the retrieval phase. In region C only the spurious states are stable, while in region D the noise dominates and only states with m = 0 are possible. We note that there is a critical value α c for the capacity above which the patterns are not stable, and that determines a bound on the capacity of a functional Hopfield network.

The picture that emerges from the mean-field analysis is rather complete and detailed. We pass from a not very well-defined energy landscape to a more refined free energy landscape in which the functioning of the network is described in detail. Finally, for T → 0 we recover the behavior of the deterministic Hopfield model.

As a final note, let us remind that all of the above results are obtained for random patterns. To describe a Hopfield model able to memorize real-world data, we need more realistic assumptions. We will deal with this fundamental problem in the chapters about the RBM.

Beyond memory: Boltzmann Machines

In the previous section we have seen how a stochastic description of the Hopfield network is useful to analyze the model in great detail. Nonetheless, we have seen that the deterministic rule (1.3) is sufficient to robustly retrieve stored patterns, successfully realizing an associative memory. In this section we will show how the introduction of stochastic units is not only useful for the purpose of analyzing the model, but it can serve to approach the task of pattern formation and data generation. We will then discuss the need to introduce hidden units to obtain more expressive networks, and define the Boltzmann machine model in general terms.

Stochastic units. Let's consider binary nodes v i = 0, 1 as opposed to s i = ±1; we will keep this notation consistent throughout the chapters. The units are stochastic, such that similarly to rule (1.8) we have

v i = 1 with probability p(x) 0 with probability 1 -p(x) with p(x) = 1 1 + exp -β j w ij x j + θ i = sigm   β j w ij x j + θ i   1.8.
Beyond memory: Boltzmann Machines where x represents the data (i.e. the patterns) and the sigmoid function sigm(x) = 1/(1 + exp(-x)) appears naturally due to the use of v i = 0, 1 binary units. Similarly to (1.6) the probability over the nodes is given by v) .

P (v) = e -βE(v) Z , Z = v e -βE(
(1.15)

The above expression defines a parametrized probability distribution over the data, with parameters W , θ. Assuming that the dataset is described by a ground truth distribution, our goal is to determine the set of parameters that best approximates such a distribution. Adding the requirement that we must be able to efficiently sample from the learned distribution, we obtain a generative model of the data.

Hidden units. The probabilistic model that we introduced above, in the current formulation, presents a major limitation. The energy function only contains first-and second-order terms in the nodes, meaning that the network can enforce only up to pairwise correlations; higher-order correlations are not contemplated. As an example, consider the configurations of even parity for a network with 3 nodes: (000), ( 011), ( 101), (110). There doesn't exist a set of values for the weights of the network such that all of those configurations can be represented by equilibrium states of the network. While in this simple case the problem could be solved by adding an extra node with the appropriate value for each configuration, in practice the extra information is not present in the dataset to model. We need to add an extra set of hidden nodes to the network, which are not used to represent the data but to introduce higher-order dependencies among the original nodes. The original nodes are called visible by virtue of the fact that they are used to directly represent the data. The kind of stochastic network with visible and hidden nodes that we described above is known as Boltzmann machine; its energy function includes three different weight matrices L, J , W to model visible-to-visible, hidden-tohidden and visible-to-hidden correlations. Denoting by v i and h j respectively the visible and hidden nodes, the energy function reads

E(v, h) = - 1 2 i,j v i l ij v j - 1 2 ij h i j ij h j - ij v i w ij h j - i θ i v i - j η j h j .
(1.16) One final problem that we need to address to use the Boltzmann Machine to model the empirical distribution of a dataset is to define a procedure to learn the parameters of the model. A training algorithm for the Boltzmann Machine has been proposed in [1]; in next chapters, we will see how dropping some parameters let us define the Restricted Boltzmann Machine, for which a more efficient training algorithm is readily obtained.

1. Background

Statistical Physics approach to neural networks

This Chapter serves as an example of how statistical physics can be useful to study neural network models in great detail. Scientific work adopting this approach has been flourishing in the last decade [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF]17], tackling the challenge of studying the theoretical foundations of deep learning. The work on the RBM that is the object of next Chapters finds its spot in this context. In Section 1.8 we have discussed the similarity between the RBM and the Hopfield model; a detailed analysis of this connection is found in [9,8,[START_REF] Mézard | Mean-field message-passing equations in the Hopfield model and its generalizations[END_REF]. What makes the RBM a reference model from the statistical physics point of view is its similarity to the Sherrington-Kirkpatrick (SK) model of a spin glass, that has been studied in detail in [START_REF] Sherrington | Solvable Model of a Spin-Glass[END_REF][START_REF] Thouless | Solution of 'Solvable model of a spin glass[END_REF]3]. Indeed, the RBM can be seen as a bipartite variant of the SK model, and its theoretical analysis is based on this observation.

Chapter 2

Restricted Boltzmann Machines (RBM)

Definition

In section 1.8 we have seen how Boltzmann Machines (BM) provide a model to go beyond an associative memory (Hopfield model) to build a generative model of the data. While this is a success from the theoretical point of view, BMs are not very practical as their training is very expensive for high dimensional data. To improve on this situation, the RBM forbids connections among visible nodes and among hidden nodes, making it possible to efficiently train the model.

The actual RBM model then consists in a bipartite graph with a layer of hidden units h j and a layer of visible units v i . The units in one layer are not connected among them but are connected to all the units in the other layer.We restrict our treatment to the case of binary units h i , v i = 0, 1. The energy function of the BM is readily adapted to the bipartite case

E(v, h) = - i,j v i w ij h j - i θ i v i - j η j h j (2.1)
where θ i and η j are external fields, or biases, acting respectively on the visible and hidden units. The probability of a configuration (1.15) becomes

P (v, h) = e -βE(v,h) Z , Z = v,h e -βE(v,h) . (2.2)
To simplify the equations, we set β = 1 without loss of generality. We recall that in Section 1.3 we discussed how β is a parameter that regulates the noise in the network; for the RBM, a proxy for β is given by the variance of the weights W , as we'll see in Section 2.3. The crucial simplification in the RBM lies in the fact that units in the same layers are not connected: they are thus conditionally independent given the nodes in the other layer, and we can compute the individual probabilities in 2. Restricted Boltzmann Machines (RBM) closed form

P (v i = 1|h) = 1 1 + e -θi-j wij hj = sigm   θ i + j w ij h j   (2.
3)

P (h j = 1|v) = 1 1 + e -ηj -i wij vi = sigm η j + i w ij v i .
(2.4)

What we are interested in is the probability for the visible units, which is the layer we use to represent data. It can be easily defined as

P (v) = h P (v, h) = e -Fc(v) Z , Z = v e -Fc(v) (2.5)
where we have introduced the clamped free energy

F c (v) = -log h e -E(h,v) = - i a i v i - j log 1 + e (bj+ i wij vi) . (2.6)
To use the RBM as a generative model, we want to maximize P (v) for the samples belonging to the training set. This is done by performing gradient ascent over the log-likelihood log P (v), whose derivative with respect to the weights can be computed to be

∂ log P (v) ∂w ij = v i P (h j = 1|v) data -v i h j model (2.7)
where • data denotes an average over the empirical distribution of the training set (P

data (v) = 1 N vn∈data δ(v -v n )
) and • model denotes the average over the distribution (2.5). Introducing the learning rate α (as a parameter for gradient ascent) we obtain an update rule for the weight matrix ∆W = α vh T datavh T model .

(2.8)

In the same way we can get the update rules for the external fields

∆θ = α (v data -v model ) (2.9) 2.2. Contrastive divergence training ∆η = α (h data -h model ) . (2.10)
Given the update rules, the training consists in actually performing the gradient ascent. Once this is done, it is possible to sample the equilibrium configurations of the RBM to obtain samples which are generated according to the approximated probability distribution of the training data. Unfortunately, the average over the model distribution • model is intractable as such term is exponential in the number of visible units. Approximations are then necessary to train and sample from an RBM; in particular, Monte Carlo based algorithms are generally employed, such as k-steps contrastive divergence (CDk) [38] and persistent contrastive divergence (PCD) [START_REF] Tieleman | Training Restricted Boltzmann Machines Using Approximations to the Likelihood Gradient[END_REF]. Approximate algorithms based on mean-field methods have also been employed, originally in [START_REF] Welling | A New Learning Algorithm for Mean Field Boltzmann Machines[END_REF] and more recently in [28].

Contrastive divergence training

The standard strategy to compute an approximation of the gradient (2.8) is to use a Montecarlo sampling procedure. The first term vh T data , also called the positive term, is easily computed: here v is represented by the actual data while h is computed with (2.4). The term vh T model , also called the negative term, is more problematic: to average over the actual parameters of the model, we can initialize v with the data and perform Gibbs sampling (Algorithm 1) for a high number of steps k, i.e. until convergence is reached.

Algorithm 1 Gibbs sampling 1: Init: take a random configuration v 2: for i = 0 to k do 3:

h ∼ P (h|v)

4:

v ∼ P (v|h) 5: end for 6: Set v i = 1 with probability P (v i = 1|h) This would give a good unbiased estimate of the negative term, at the cost of a lengthy iterative procedure. In practice, instead of running the Montecarlo procedure until equilibrium is reached, the negative term is usually estimated by performing only a small number of Gibbs steps k, and this training procedure is called k-steps Contrastive Divergence (CDk) [38]. CDk uses a rather crude approximation to the gradient; a slightly refined training strategy consists in using the data to initialize the Gibbs sampling procedure only at the first gradient update, and keep iterating over the same chain for subsequent updates. The rationale is that every gradient update perturbs the weight matrix only slightly, meaning that after many minibatches and epochs of training we can expect the sampled configurations to better approximate the equilibrium state of the model. The sampled configurations are called the persistent chain and the procedure k-steps Persistent Contrastive Divergence (PCDk) [START_REF] Tieleman | Training Restricted Boltzmann Machines Using Approximations to the Likelihood Gradient[END_REF]. Both

Restricted Boltzmann Machines (RBM)

CDk and PCDk have been analyzed in detail in [22], which shows that the RBM can operate in two different regimes depending on the interplay between the number of steps k and the mixing time of the Gibbs sampling procedure. For small values of k (relative to the mixing time) the RBM operates in an out-of-equilibrium regime in which it is possible to produce good qualitative samples but whose data generation abilities don't properly match the training dataset. To obtain a properly equilibrated model, high values of k are necessary and in this case the RBM is able to produce good qualitative samples that also respect the expected statistics of the training set. In Algorithm 2 we detail the CDk training strategy in vectorized (matrix) form. Here, we choose to randomly initialize the negative term's sampling chain as proposed in [21] while traditionally the negative chain is initialized with the dataset samples as for the positive term [40,26]. Moreover, we note that in principle the number of samples in the negative chain can be freely chosen and we expect to better approximate the negative term by using a high number of samples. However, the number of samples in the negative chain is generally chosen to be equal to the batch size, which can be motivated by the need to keep the statistical errors in between positive and negative terms comparable. for all minibatches m do Positive term:

6: V m = [v 1 v 2 . . . v n ] 7:
sample

H p ∼ P (H p | V m ) (Equation (2.4)) 8: sample V p ∼ P (V p | H p ) (Equation (2.3)) 9:
compute the positive term vh T data = end for 15: end for

Mean-field training

We have seen in Section 1.9 how mean-field theories are a powerful tool to study the kind of neural network models we are dealing with. Here we follow [28] to introduce a mean-field strategy to compute the gradient of the log-2.3. Mean-field training likelihood of the RBM. The advantage of this method is that we can substitute the Montecarlo sampling with deterministic iterative mean-field equations.

To derive an effective free energy we explicitly reintroduce the inverse temperature β

P (v, h) = e -βE(v,h) Z . (2.11)
Exploiting the similarity of the RBM to the SK model discussed in Section 1.9, we follow the derivation in [31] and adapt it to the bipartite case. This consists in expanding the free energy at high temperature by setting β → 0 to obtain the Thouless-Anderson-Palmer (TAP) expression for the effective free energy [START_REF] Thouless | Solution of 'Solvable model of a spin glass[END_REF] that, truncated at second order, is given by

F T AP (m v , m h ) = + S(m v ) + S(m h ) - i a i m v i - j b j m h j - i,j w ij m v i m h j + i,j w 2 ij 2 m v i -m v i 2 m h j -m h j 2 (2.12) with S(m) = -i [m i log m i + (1 -m i ) log(1 -m i )]
and with m v , m h being the equilibrium magnetizations of visible and hidden nodes. The minimization of the above Equation (2.12) gives a valid approximation to the free energy F F F T AP ( mv , mh ), dF T AP dm mv , mh = 0.

(2.13)

To obtain mv , mh it is then necessary to extremize (2.12) to obtain the following coupled equations

m v i sigm    θ i + j w ij m h j -w ij m v i - 1 2 m h j -m h j 2    (2.14) m h j sigm η j + i w ij m v i -w ij m h j - 1 2 m v i -m v i 2 (2.15)
that can be solved by iteration [11]. Finally, given the approximation to the free energy (2.13), the optimization problem over the log-likelihood (2.7) is greatly simplified: the average over the training samples • data is unchanged while the intractable average over the model distribution • model is substituted by the maximization of (2.13), giving

∆W = α vh T data - ∂F T AP ( mv , mh ) ∂w ij (2.16) with ∂F T AP ∂w ij = - i,j m v i m h j + w ij (m v i -m v i 2 )(m h i -m h i 2 ).

Restricted Boltzmann Machines (RBM)

Summarizing, the training procedure based on TAP approximation is reported in Algorithm 3. We note that similarly to CD/PCD, instead of running the TAP equations (2.14)(2.15) to convergence we choose a fixed number of iterations k to approximate the equilibrium results. for all minibatches m do Positive term:

6: V m = [v 1 v 2 . . . v n ] 7:
sample

H p ∼ P (H p | V m ) (Equation (2.4)) 8: sample V p ∼ P (V p | H p ) (Equation (2.3)) 9:
compute the positive term vh

T data = 1 n V p H T p
Negative term: end for 15: end for We conclude by noting how the reintroduction of the inverse temperature β is just a formal passage; in the context of the RBM the high-temperature expansion β → 0 is substituted by a weak-couplings expansion, under the assumption that the weights w ij are small enough. The variance of the weight matrix can then serve as an effective inverse temperature [START_REF] Monasson | Emergence of Compositional Representations in Restricted Boltzmann Machines[END_REF] T ef f = 1 V ar(W ) .

10: initialize the batched magnetizations M v = V m , M h = P (H | V m )
(2.17)

Generalized RBM

In Section 2.1 we introduced the RBM using binary units for both the visible and the hidden layers. While this is the original and the most commonly used formulation of the model, we can generalize the definition to use different kinds of units. This is done by introducing a prior probability distribution over visible and hidden variables, that we denote respectively by q v and q h :

p(v, h) = 1 Z q v (v)q h (h)e -E(v,h) .
(2.18)

Generalized RBM

The conditional probability for visible and hidden layers read (Bayes formula)

p(h|v) = p(v, h) h p(v, h) = q h (h)e -E(v,h) h q h (h)e -E(v,h) (2.19) p(v|h) = p(v, h) v p(v, h) = q v (v)e -E(v,h) v q v (v)e -E(v,h) .
(2.20)

For the standard binary-binary RBM, a Bernoulli prior is used for both visible and hidden units

q v (v i ) = 1 2 (δ vi,0 + δ vi,1 ) (2.21) q h (h j ) = 1 2 (δ hj ,0 + δ hj ,1 ) (2.22)
which leads to Equations (2.3),(2.4) for the conditional activations.

When dealing with data whose features are better represented by real numbers, we can substitute the binary variables in the visible layer with Gaussian units. The prior distribution in this case is

q v (v i ) = 1 2πσ 2 v e - v 2 i 2σ 2 v (2.23)
and the conditional activation of the visible layer becomes

p(v i , h) = N   v i ; j w ij h j + θ i , σ 2 v   (2.24)
where N (x; µ, σ 2 ) is the Gaussian distribution of variable x with mean µ and variance σ 2 . In practice, when using Gaussian visible units, the data are centered and normalized as a preprocessing step and the variance σ 2 v is set to 1 and kept fixed.

Other choices of priors for the units of the RBM are possible, as well as different combinations of visible and hidden units types; a detailed account can be found in [21]. In this document we are only concerned with the Bernoulli-Bernoulli and Gaussian-Bernoulli models described above.

Chapter 3 Spectral Learning Dynamics of the RBM

Concerning the learning procedure of neural networks, many recent statistical physics based analyses have been proposed, most of them within teacherstudent setting [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF]. This imposes a rather strong assumption on the data: it is assumed that these are generated from a model belonging to the parametric family of interest, hiding as a consequence the role played by the data themselves in the procedure. From the analysis of related linear models [START_REF] Tipping | Mixtures of Probabilistic Principal Component Analyzers[END_REF]12], it is already a well established fact that a selection of the most important modes of the Singular Values Decomposition (SVD) of the data is performed in the linear case. In fact in the simpler context of linear feed-forward models the learning dynamics can be fully characterized by means of the SVD of the data matrix [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural network[END_REF], showing in particular the emergence of each mode by order of importance with respect to the corresponding singular values.

In this section we present some empirical results to qualify the role of the SVD of the weight matrix of the RBM during learning. We show how the learning procedure follows the dynamics of the strongest SVD modes, in a data-driven fashion. This will let us single out the information content of the RBM, leading us to formulate the assumption that the SVD spectrum is split in a continuous bulk of singular vectors corresponding to noise and a set of outliers that represent the information content. This assumption forms the basis for the mean-field analysis presented in Chapter 4.

The results presented here have been originally detailed in [19] (Reprint A). We will briefly introduce the SVD technique and its relation to Principal Components Analysis. Subsequently, we show how the linearized mean-field equations for the RBM naturally suggest the use of the SVD. We then go on to observe the SVD dynamics in a real-world scenario, characterizing the SVD modes and the statistical assumptions over the weight matrix.

Spectral Learning Dynamics of the RBM

Principal Component Analysis (PCA)

The PCA technique can be introduced by considering the covariance matrix of a dataset. Given a data matrix X of dimension n × d with n > d, where n is the number of samples and d is the dimension of each sample, and further assuming that samples are centered (i.e. column means have been subtracted, as data are arranged by rows), we can define an unbiased estimator for the related covariance matrix (square and symmetric)

C = X T X n -1 (3.1)
that can be diagonalized

C = V LV T (3.2)
where the columns of V are eigenvectors of C and L is the diagonal matrix of the eigenvalues λ α . Projecting the samples over the eigenvectors of the covariance matrix (also called principal directions in the context of PCA) we obtain the principal components: new, independent variables that account for the maximum possible variability in the data. More precisely, the first principal component maximizes the variance of the projections of the data (i.e. it has the highest possible variance) and the succeeding components maximize the variance while satisfying the constraints of being orthogonal to the preceding components. A rigorous demonstration of the properties of principal components is given in [10], Section 12.1.

Singular Value Decomposition (SVD)

The SVD is the generalization of eigenmodes decomposition to rectangular matrices, and it is given by

X = U ΣV T (3.3)
where U is an orthogonal n × d matrix whose columns are the left singular vectors u α , V is an orthogonal d × d matrix whose columns are the right singular vectors v α and Σ is a diagonal d × d matrix whose elements are the singular values w α . The separation into left and right singular vectors is due to the rectangular nature of the decomposed matrix, and the similarity with eigenmodes decomposition is revealed by the following SVD equations

Xv α = w α u α (3.4) X T u α = w α v α . (3.5)
Plugging the SVD of X into the definition of covariance matrix

C = X T X n -1 = V ΣU T U ΣV T n -1 (3.6) = V Σ 2 n -1 V T 3.
3. Linearized mean-field equations for a RBM we see how the right singular vectors can be identified as the principal directions and a relation between the singular values w α and the eigenvalues of the covariance matrix λ α is easily found

λ α = w 2 α n -1 (3.7)
Finally, the principal components are given by U Σ (XV = U ΣV T V = U Σ).

Linearized mean-field equations for a RBM

To see how the SVD enters the picture, we linearize the mean-field equations [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF] of our model. As there are no connections among variables in the same layer, we can initialize the magnetizations m v i , m h j with Equations (2.3),(2.4) and write the iterative mean-field equations

m v i = sigm   θ i + j w ij m h j - j w ij   (3.8) m h j = sigm η j + i w ij m v i - i w ij . (3.9)
At initialization the weights w ij are small, and we can get rid of the external visible field by centering the training data. The external hidden field is instead initialized to zero and it varies slowly, so it doesn't have any effects at the beginning of the training. Thus, neglecting both the external fields we can linearize the mean-field equations to obtain (defining mv i = m v i -1/2, mh j = m h j -1/2 for convenience)

mv i 1 4 j w ij mh j (3.10) mh j 1 4 i w ij mv i . (3.11)
We can now express the weights w ij in terms of the SVD as (u i,α identifies the i th component of the α th columns of U , and analogous notation is used for V )

w ij = α w α u i,α v j,α (3.12)
and expand the magnetizations over the singular vectors Combining Equations (3.10),(3.12),(3.13) and recalling that the columns of U form an orthonormal basis, we get

mv α = i u i,α mv i (3.13) mh α = j v j,α mh j . ( 3
mv α = 1 4 α w α δ α,α mh α = 1 4 w α mh α .
We can proceed in an analogous way with mh α to finally obtain the coupled iterative equations

mv α = 1 4 w α mh α (3.15) mh α = 1 4 w α mv α . (3.16)
These equations show that at the linear level the magnetizations aligned to the singular vectors with a strong w α are amplified, while magnetizations related to small w α are penalized. We then expect the samples generated by a trained RBM to be affine to the strongest singular vectors, and we can try to understand how. To this end, we can better specify the role of the SVD matrices in the context of a RBM:

• U encodes the singular vectors related to the visible layer; these can be visualized in the pixel space and basically consist in the principal components of W .

• V is related to the hidden layer; it is a square orthogonal matrix that can be interpreted as a rotation and its columns are the principal directions of W .

• The singular values w j contained in Σ can be thought of as scaling factors whose action is to weigh the singular vectors composing W .

Given the above characteristics we focused our attention on Σ and U , tracking the distribution of the singular values and looking at the corresponding left singular vectors during the training.

Distribution of the singular values

Distribution of the singular values

We now turn to the analysis of a training run: we use PCD to train the RBM over the MNIST dataset [50] and we monitor the evolution of the SVD components. The MNIST dataset is composed by 70000 images of handwritten digits (60000 for training and validation, 10000 for testing) of size 28×28 pixels. Some samples from the dataset are shown in Figure 3.1. The weight matrix W is initialized as a Gaussian random matrix with variance σ w (and zero mean).

The eigenvalues distribution of the corresponding symmetric square matrix W T W is known to be given by the Marchenko-Pastur law [56] in its canonical form. The singular values w α are related to the eigenvalues λ j by (3.7), and by defining the parameter r = N h /N v , with N h the size of the hidden layer and N v the size of the visible layer (we recall that W is a N v × N h matrix), the expression of the Marchenko-Pastur law is given by (in the limit N h , N v → ∞ with r finite)

ρ(λ) = 1 2πσ 2 w (λ -r -)(r + -λ) rλ (3.17)
where the higher and lower bounds r ± are

r ± = σ 2 w 1 ± √ r 2 . (3.18)
Figure 3.2a shows the agreement between the empirical distribution and the theoretical distribution. In particular, we note how all w α have values below the threshold set by the Marchenko-Pastur law, forming a bulk of singular values. Starting with the training, we see that many singular values increase in magnitude and overcome the threshold for a Gaussian random matrix; these are outliers leaving the bulk, shown in Figure 3.2b-3.2d. During the first epochs of training this process is very fast and many w α values are easily extracted from the bulk, growing by many orders of magnitude. The bulk is instead shrunk to low values, meaning that the w α values that do not overcome the threshold decrease in magnitude. Going on with the training this process slows down, but it does not stop: outliers keep growing slowly, and the bulk keeps shrinking to approach a spike around zero. It is important to note that a kind of hierarchy is maintained in the process: the first outliers are never overcome by the newly extracted w α , and this is made clear by looking at the corresponding left singular vectors (see next Section). In the final epochs of training, the singular values w α are separated into two categories: a concentrated set of almost-null singular values and a set of outliers spread above the threshold, as shown in Figure 3.2e. The evolution of the w α distribution described above suggests that the training process is able to discern between the most important singular vectors, that are brought above threshold first and heavily strengthened, and a bulk of less important singular vectors, that end up above threshold but whose w α reach values order of magnitude smaller than those of the strongest singular vectors. Moreover, the below-threshold singular vectors are practically eliminated by cutting down the corresponding w α . These observations give a good indication about what are the dynamics of the learning process, but a couple of matters need to be addressed: (i) it is not clear what singular vectors actually represent, (ii) the meaning of more and less important singular vectors has to be specified. We will deal with these problems in the next Section.

Dynamics of the singular vectors

To understand the role of the left singular vectors of an RBM we must keep in mind the interpretation for the SVD decomposition of W given previously.

We have seen how the matrix of singular values Σ is shaped during learning, and we recall that V is interpreted as a rotation in the space of the hidden units. We then expect to recover the structure of the training data into the U matrix; using the MNIST dataset proves useful in this context as we can visualize the u α vectors as images in the pixel space. Before focusing on the left singular vectors, we note that also the external visible field can be visualized in the pixel space. Following [40] we can initialize the field to where p i is the proportion of training samples in which unit i is active. Using the above rule (3.19) we are able to encode into the field the mean activations of the visible layer, which is clearly shown in Figure 3.3b in the pixel space.

θ i = log[p i /(1 -p i )] (3.19)
If we instead initialize the visible field with a null vector, the mean activation pattern is learned very effectively as the strongest left singular vector u α . The striking resemblance between the mean activation pattern computed from the training data and the one learned by the RBM is shown in Figure 3.3a-3.3b and it serves as a first example of what the u α vectors represent. It seems then equivalent to either encode the mean activation pattern into the visible field since the beginning or letting the RBM learn such a pattern as a left singular vector. In practice the second case is not desirable as the RBM associates to the mean activation pattern a very strong singular value, many orders of magnitude higher than the strongest outliers. This results in a bias in the sampling from the trained machine, such that the samples whose activation pattern is nearest to the mean are sampled with a higher frequency (in the worst case, those are the only configurations sampled at equilibrium). The first 10 left singular vectors of a trained RBM are shown in Figure 3.3e-3.3f. They are all composed by a homogeneous background on the borders and a set of alternating dark and light traits in the center, highlighting the fact that each singular vector acts globally on the visible layer. Even if the pictures seen in Figures 3.3e-3.3f are quite different one from another, an interesting trend is found: a higher number of alternating traits is present in the successive vectors.

Spectral Learning Dynamics of the RBM

These observations suggest that the RBM is able to learn the modes that compose the activation patterns of the data, starting with the low frequency modes and proceeding with the high frequency ones. Moreover, with reference to the w α distribution (Figure 3.2e), we note that the low frequency modes are given a higher weight. Recalling the connection between SVD and PCA (Section 3.2) these modes are understood as the principal modes of variation of the data. In Figures 3.3d-3.3e we compare the SVD modes extracted from the data and those of the W matrix, which prove to be very similar.

The dynamics described here present some similarities to the learning dynamics of deep linear neural networks [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural network[END_REF]. Going on with the training we expect non-linear effects to kick-in and this is seen in Figure 3.3f where the SVD of the data is not comparable to the SVD of W anymore.

Analyzing the learning dynamics more in detail, we observed that the modes take shape one by one as the corresponding singular value w α is brought above threshold. The subsequent strengthening of the w α values corresponds to refinements and rotations, with little effects on the characterization of the modes as high or low frequency modes. For what concerns the modes below threshold, they present a dark border and a random configuration in the center; in this case the only effect of the training is to discern what are the units which are never activated and no information about the actual structure of the data is found.

Summarizing, some insights on the behavior of an RBM in the linear regime were given by looking at the SVD-like equations (3.15)- (3.16), where we have seen how the magnetizations aligned to the strongest SVD modes are amplified. These magnetizations are thus unstable and they drive the formation of new mean-field fixed points during learning, that correspond to the magnetizations affine to the samples in the training set. These observations highlight a connection between the SVD of the data in the training set and the SVD of the weight matrix W , at least in the linear regime.

Characterization of the modes

By looking at the singular values distribution of W we have seen that there seem to be more and less important singular vectors. In the previous Section, we have then refined this observation by highlighting how the lowest-frequency modes are given the highest weights. We can then identify the more (less) important modes as the low (high) frequency ones. To gain some intuition about the meaning of this separation we can think about the Fourier decomposition of a square wave; in such a case the superposition of the low frequency harmonics is sufficient to build a good approximation of a square wave, while the role of the high frequency harmonics is that of sharpening the waveform at the discontinuity points. In the context of a trained RBM, we then expect that good approximations to the training data are obtained by exploiting only the low frequency modes, while the high frequency modes should represent minor corrections. To discern between high and low frequency modes, we look at the w α values in decreasing order on a log-log plot (Figure 3.4): the strongest w α values are located far above threshold and are of comparable magnitude, followed by a tail of exponentially damped values. This picture is consistent across the training, the only difference being the damping cutoff that is increasing with the epochs. After a relatively long training, however, the increase in the cutoff is very slow and this could serve as a signal to stop the learning, as such situation amounts to slowly strengthening singular values which are exponentially less important than the already learned modes. We are then driven to define the more important low frequency modes as the modes before cutoff, and the exponentially less important high frequency modes as those after cutoff. A consistency check is shown in Figure 3.5, where just the 100 strongest modes are retained to construct the samples and the remaining modes are shown to encode boundary corrections. Choosing the first 100 modes is arbitrary; in Figure 3.4 we can see how the cutoff is well below 100, so with this choice we are sure that we included in the reconstruction of the samples all the strong modes before cutoff plus a small number of modes giving boundary corrections.

As a conclusion, the above observations suggest that the weight matrix of a trained RBM is composed by two classes of modes; recalling the expansion (3.12) we can express the components of W as

w ij = α∈bulk w α u α i v α j + α∈outliers w α u α i v α j (3.20)
where the modes of the bulk are those that remain random after the training (below threshold or after cutoff) while the outliers correspond to the modes that actually encode the structure of training data. This separation sets the basis for the mean-field analysis that we develop in the next section.

Chapter 4

Mean-field theory

In Section 1.9 we have discussed how statistical physics tools have been fruitfully applied to model neural network models. In particular, various works [2,44,43,8,[START_REF] Monasson | Emergence of Compositional Representations in Restricted Boltzmann Machines[END_REF] have studied the equilibrium properties of the RBM in order to understand collective phenomena in the latent representation. The common assumption is that the components of the weight matrix are i.i.d., and some insights about the equilibrium phase of a trained RBM are given. However, this approximation is problematic since the learning procedure introduces correlations among the weights making the assumption of independent weights unsuitable to describe a realistic model of the RBM.

In this Chapter we introduce a statistical ensemble for the weight matrix W in which we drop the assumption of independent components and we apply the Replica Method to derive a mean-field theory of the RBM and characterize the equilibrium states of the model. These results have been presented originally in [20] (Reprint B).

Statistical ensemble

In Chapter 3 we have seen how the evolution of the singular values of the trained model departs strongly from a Marchenko-Pastur distribution, which would be expected if the weights were i.i.d. Gaussian. In particular, we have seen how the information content of the trained RBM is encoded into a finite number of spectral modes, while a residual part of the spectrum constitutes a noise source. 1 We can then assume the weight matrix W to be composed by K spectral modes and a portion of random noise n

W ij = K α=1 w α u α i v α j + r ij . (4.1)

Mean-field theory

The vectors u α , v α represent the K strongest components of the SVD of the weight matrix, while the components r ij represent Gaussian noise. The above formula describes a class of models that we can analyze statistically, i.e. the statistical ensemble that we choose to model the RBM. In the following, we derive the static and dynamical properties of this ensemble.

Order parameters and effective free energy

To obtain an effective free energy (Section 1.6) for the RBM we apply the Replica method [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF], which is based on the following identity

E u,v,r [log Z] = lim p→0 d dp E u,,v,r [Z p ] (4.2)
where u, v are the singular vectors components and r the Gaussian noise defined in Section 4.1. The v α and u α vectors are orthonormal,

giving v ∼ O 1/ √ N h and u ∼ O 1/ √ N v
, and the noise terms are i.i.d. r ij = N (0, σ 2 /L) where N v , N h are the numbers of visible and hidden units and L = √ N v N h . We see from this equation that our objective is to average the free energy over the noise and the components of the singular vectors; the second average is the crucial one, as it represents the average over the structured part of the weight matrix, i.e. the part that encodes the correlations among the weights. Nonetheless, the components of the vectors u α and v α are assumed to be i.i.d., and we will see that the choice of their distribution determines the properties of the equilibrium states.

An important assumption in the replica method is that the exponent p assumes integer values, hence the replicated partition function Z p can be thought of as the partition function of a system comprising p copies, i.e. the replicas, of the original system. We will see that exploring the correlations among the different replicas of the system will be crucial to fully characterize the equilibrium behavior of the system. For now the only concern is that Equation (4.2) requires taking the limit p → 0 while here we are assuming p to be integer; we will discuss this incongruence at the end of the Section.

We can start by writing down the replicated partition function

Z p Z p = p a=1 s a ,σ a exp    - i,j,α s a i u α i w α v α j σ a j - i,j s a i σ a j r ij    × exp    - i η i s a i - j θ j σ a j    . (4.3)
By isolating the part of Z p that contains the noise terms r ij we can compute the noise average to obtain exp

  σ 2 2L a s a i σ a j 2   = exp   σ 2 2L   p + a =b s a i s b i σ a j σ b j     .

Order parameters and effective free energy

Considering now the expansion of the units over the singular modes

s a α = 1 √ L i s a i u α i , σ a α = 1 √ L j σ a j v α j (4.4)
and neglecting the fields η i , θ j for a moment we can rewrite the replicated partition function as

Z p ∝ p a=1 s a ,σ a exp    -L α w α s α σ α + σ 2 2L i,j,a =b s a i s b i σ a j σ b j    . (4.5)
We observe that both sums in the exponential include higher order terms that we cannot factorize to simply compute the equilibrium averages over u and v.

Following the strategy presented in Section 1.7 for the Hopfield model, we can introduce an appropriate set of order parameters to linearize the exponents.

For the first term we can consider the magnetization of the system m a α , ma α correlated with the spectral modes, expressed as

m a α ∼ E u,v,r ( σ a α ) , ma α ∼ E u,v,r ( s a α ) (4.6)
and introduced with the following integral identity (Hubbard-Stratonovich transformation):

exp L α w α s a α σ a α ∝ α dm a α d ma α 2π × exp -L α w α (m a α ma α -m a α s a α -ma α σ a α ) .
We note that the magnetizations m a α , ma α play the same role as the magnetization over the stored patterns in the Hopfield model; for the RBM we effectively substitute the patterns with the spectral modes. The second sum in (4.5) stems from the average over the noise r ij , whose effect has been to couple the different replicas of the system. To treat the now interacting replicas, we take into consideration their correlations by introducing the following order parameters:

Q ab ∼ E u,v,r σ a i σ b i , Qab ∼ E u,v,r s a j s b j . (4.7)
These are the Edwards-Anderson order parameters and are fundamental to paint a detailed description of the equilibrium states of the system. They are introduced with the following Hubbard-Stratonovich transformation exp

  σ 2 2L   i,j,a =b s a i s b i σ a j σ b j     = a =b dQ ab d Qab 2π × exp   - Lσ 2 2 a =b   Q ab Qab - Q ab N v i s a i s b i - Qab N h j σ a j σ b j    

Mean-field theory

whose utility is again to simplify the expression by linearizing the sum in the exponent w.r.t. the nodes of single replicas. Putting together the above manipulations and including the fields, we rewrite the full average as

E u,v,r [Z p ] = a,α dm a α d ma α 2π a =b dQ ab d Qab 2π × exp   -L   a,α w α m a α ma α + σ 2 2 a =b Q ab Qab - 1 √ κ A(m, Q) - √ κB( m, Q)     (4.8) with κ = N h Nv and A(m, Q) def = log   s a ∈-1,1 E u   exp   √ κσ 2 2 a =b Q ab s a s b + κ 1 4 a,α (w α m a α -η α )u α s a       B( m, Q) def = log   σ a ∈-1,1 E v   exp   √ κσ 2 2 a =b Qab σ a σ b + κ 1 4 a,α (w α m a α -θ α )v α σ a       .
where the fields are also expanded over the spectral modes and the projections are assumed to be O(1):

η α def = 1 √ L i η i u α i = O(1) (4.9) θ α def = 1 √ L j θ j v α j = O(1). (4.10)
The expression in Equation (4.8) makes it possible to compute the average over u, v with the saddle point method by letting L → ∞ (thermodynamic limit), which amounts to considering the average over a large system (we recall that

L = √ N v N h ).
In order to proceed, we need to address a couple of issues: i) we need to specify an explicit dependence of the order parameters on the replica indices a, b;

ii) given the solution to the saddle point equations, we need to consider the analytic continuation p → 0 as in (4.2) to obtain the actual average free energy. At this point, it is unclear how to perform this limit.

As for point i) the simplest assumption is that the order parameters should not depend on the replica index; this seems to make sense given that replicas have been introduced artificially in Equation (4.2) for mathematical convenience and therefore we expect that they do not have any effects on the physics of the system. Point ii) is more subtle, and a possible approach is to consider a 4.2. Order parameters and effective free energy specific structure for the order parameters such that the saddle point equations are analytic in p [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF].

The above assumptions are realized by the replica symmetric ansatz

Q ab = δ ab + (1 -δ ab )q (4.11) Qab = δ ab + (1 -δ ab )q (4.12)
in which the dependence on the replica index is dropped and the set of Edwards-Anderson order parameters is reduced to the couple of q, q parameters. The resulting Q matrices are symmetric, filled with unity on the diagonal and with q, q off-diagonal, and they constitute a simple example of an ultrametric matrix. This simple structure makes it easy to plug (4.11), (4.12) into (4.8) to transform the sums involving the Q ab , Qab terms into analytic expressions involving p and q, q. We can now take the limit p → 0 to obtain the effective free energy

f (m, m, q, q) = α w α m α mα - σ 2 2 q q + σ 2 2 (q + q) - 1 √ κ E u,x [log 2 cosh(h(x, u))] - √ κ E v,x log 2 cosh h(x, v) . (4.13)
where x ∼ N (x; 0, 1). Finally, in the limit L → ∞ we extremize the effective free energy to obtain the saddle point equations

m α = κ 1 4 E v,x [v α tanh h(x, v) ] q = E v,x [tanh 2 ( h(x, v))] (4.14) mα = κ -1 4 E u,x [u α tanh(h(x, u))] q = E u,x [tanh 2 (h(x, u))] (4.15) with h(x, u) def = κ 1 4 σ √ qx + γ (w γ m γ -η γ )u γ h(x, v) def = κ -1 4 σ √ qx + γ (w γ mγ -θ γ )v γ .
The numerical solution of (4.14), (4.15) characterizes the equilibrium states of the system and in the next Section we will see what are the different phases of operation of the RBM and what is the impact of the choice of distribution for the u, v components.

To conclude this Section, let us remark that to obtain Equation (4.13) we need to take the limit p → 0 coming from Equation (4.2) and the thermodynamic limit L → ∞ by carefully applying the Replica method's recipe. Indeed, there is no mathematical justification for taking the analytic continuation p → 0 after having treated p as an integer; moreover, we freely exchanged 4. Mean-field theory the order of the p → 0 and L → ∞ limits. While lacking a rigorous mathematical justification, the Replica method has been shown to provide the correct equilibrium description for the SK model and other similar models to which it has been applied [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF]; for the RBM, our experimental results (Section 4.8) are compatible with the theoretical description.

Phase diagram

The state of the system described by the solutions to Equation (4.14)(4.15) is directly characterized by the value of the order parameters. The stable solutions are determined by looking at the values of the order parameters for which the Hessian of the effective free energy is positive definite. In turn, this let us identify the unstable modes that define the lines of separation among the different phases. In the basic case with no biases we identify three different phases:

• a paramagnetic phase (q = q = m α = mα = 0) (P),

• a ferromagnetic phase (q, q, m α , mα = 0) (F),

• a spin glass phase (q, q = 0; m α = mα = 0) (SG).

The full phase diagram is found in Figure 4.1.

Learning phase

For an RBM learned on real-world data we expect that at equilibrium the nodes condense over the spectral modes, determining a certain magnetization of the system. We are thus interested in fully characterizing this learning phase, which corresponds to the ferromagnetic phase introduced in Section 4.3 with m α , mα = 0. Our objective is describing the magnetizations of the system in terms of a combination of spectral modes, in order to obtain compositional states that are suited to represent realistic data. This is in analogy to [START_REF] Monasson | Emergence of Compositional Representations in Restricted Boltzmann Machines[END_REF] in which the equilibrium states are described in terms of a composition of hidden nodes.

To our end, we are going to show how the distribution of the singular vector components u, v determines the structure of the learning phase and how the choice of an appropriate distribution is fundamental to describe a realistic system. As a starting point, we consider Gaussian i.i.d. components. Once again we discard the biases, so that the Gaussian averaging of the magnetization part of the saddle point equations (4.14)(4.15) reads Reproduced from [20].

m α = w α mα (1 -q) (4.16) mα = w α m α (1 -q) (4.
from which we see that the singular value of a mode contributing a non-zero magnetization respects (singular values being non-negative)

w α = 1 (1 -q)(1 -q) . (4.18)
With non-degenerate singular values the above eq. ( 4.18) can be satisfied by a single spectral mode, revealing how the i.i.d. Gaussian averaging falls short of realizing the compositional phase we are looking for.

To see why the Gaussian averaging fails and how to fix it, we can approach the averaging in a more general fashion. Assuming the distribution of u and v to be even and introducing the auxiliary distributions

p * (u) def = - u -∞ xp(x)dx = ∞ |u| xp(x)dx, p(x) = p(-x) p α (u) def = p * (u α ) β =α p(u β )

Mean-field theory

we can define the new variables q α , qα

q α = dx e -x 2 /2 √ 2π dvp α (v) tanh 2 κ -1/4 σ √ qx + γ w γ mγ v γ (4.19) qα = dx e -x 2 /2 √ 2π dup α (u) tanh 2 κ 1/4 σ √ qx + γ w γ m γ u γ (4.20)
and rewrite the averaging in the saddle point equations (4.14)(4.15) as

m α = w α mα (1 -q α ) (4.21) mα = w α m α (1 -qα ). (4.22)
from which we get a generalized version of Equation (4.18) for the singular value of a mode contributing to the magnetization

w α = 1 (1 -q α )(1 -qα ) def = w(q α , qα ). (4.23)
To explicitate the conditions under which multiple modes are able to condense and contribute to the magnetization, we can assume that a single mode α has condensed and consider the onset of a mode β(β = α). The stability condition for mode β coexisting with mode α is given by setting w γ = 0, ∀γ = α, β in the expression of the effective free energy (4.13) and looking for the positive definiteness of its Hessian. This gives a stability gap ∆w α def = w(q, q)w(q α , qα ) (4.24)

for which an arbitrary β mode is unstable if

w β < w α + ∆w α . (4.25)
In the Gaussian case the distribution p * that we introduced is still Gaussian (p * (u) = p(u)), meaning that the q α , qα parameters don't really depend on α and resulting in a null gap ∆w α = 0. This discards the possibility of a compositional phase as all modes β with w β < w α are unstable and thus only the strongest spectral mode contributes to the magnetization of the system. Considering other distributions for the u, v components results in obtaining a set of different q α , qα parameters for each α mode, impacting the value of the stability gap. With a negative gap, the possibility of multiple spectral modes being stable is allowed and thus a compositional phase is possible. In the inset of Figure 4.1 the stability gap is traced for various simple distributions. Of particular interest is the Laplace distribution for the spectral components, as it presents a negative gap, making it a good candidate for a realistic theoretical description of the RBM. We will exploit this fact in the following.

Learning equations

Learning equations

Building on the insights given in Chapter 3 and the above Sections, we will exploit the decomposition of the learning equations (2.8), (2.9), (2.10) over the SVD components. Introducing a time variable t we rewrite (4.1) as

w ij (t) = α w α (t)u i,α (t)v j,α (t) (4.26)
where we have discarded the noise term r ij in (4.1) as at the onset of the learning the weight matrix W is random and thus we don't need to consider a separate set of noisy modes. Moreover, in order to define deterministic dynamics, we also discard the noise due to the stochastic ascent optimization by taking the continuous limit of eqs. (2.8), (2.9), (2.10) to obtain

dw ij dt = s i σ j data -s i σ j model dη i dt = s i data -s i model dθ j dt = σ j data -σ j model .
Projecting the above equations on the SVD basis we obtain

1 L dW dt αβ = s α σ β data -s α σ β model (4.27) 1 √ L dη dt α = s α data -s α model (4.28) 1 √ L dθ dt α = σ α data -σ α model (4.29)
with

s α = 1 √ L i s i u i,α , σ α = 1 √ L j σ j v j,α . (4.30)

Mean-field theory

The left-hand side of (4.27),(4.28),(4.29) can be expressed in a simple form by expanding it over the basis defined by the SVD

dW dt αβ = ij u i,α dw ij dt v j,β = δ α,β dw α dt + (1 -δ αβ ) w α dv α,T dt v β + w β u α,T duβ dt = δ α,β dw α dt + (1 -δ αβ ) w α Ω u αβ + w β Ω v βα (4.31) 1 √ L dη dt α = dη dt - β Ω u αβ η β (4.32) 1 √ L dθ dt α = dθ dt - β Ω v αβ η β (4.33)
where we have defined the generators of rotations in both u α and v α bases

Ω u αβ (t) = du α,T dt u β (4.34) Ω v αβ (t) = dv α,T dt v β . (4.35)
Finally, we can express the rotation generators in terms of the other quantities

Ω u αβ (t) = - 1 w α + w β dW dt A αβ + 1 w α -w β dW dt S αβ (4.36) Ω v αβ (t) = 1 w α + w β dW dt A αβ + 1 w α -w β dW dt S αβ (4.37) with dW dt A,S αβ def = 1 2 ( s α σ β Data ± s β σ α Data ∓ s β σ α RBM -s α σ β RBM ) . (4.38)

Linear regime

The random initialization of W generally entails that the magnetization of the system is small in the beginning, so that we can describe the onset of learning by considering an expansion of the mean-field free energy up to second order in the magnetizations. For the RBM, we can adapt the expression of the effective free energy (1.11) to contemplate a bipartite structure and non-constant w ij 4.6. Linear regime couplings to obtain

F (m v , m h ) 1 2 N i=1 (1 + m v i ) log(1 + m v i ) + (1 -m v i ) log(1 -m v i ) (4.39) + 1 2 M j=1 (1 + m h j ) log 1 + m h j + (1 -m h j ) log 1 -m h j (4.40) - i,j w ij m v i m h j + N i=1 a i m v i + M j=1 b j m h j (4.41) 1 2 N i=1 (m v i ) 2 + 1 2 M j=1 m h j 2 - ij w ij m v i m h j + N i=1 a i m v i + M j=1 b j m h j .
(4.42)

The corresponding probability measure at fixed magnetization is Gaussian (1.11), for which we can write the covariance matrix in the following form (σ v , σ h being the variances of visible and hidden magnetizations)

cov(m v , m h ) =   σ -2 h σ -2 v σ -2 h -WW T W 1 σ -2 v σ -2 h -W T W W T 1 σ -2 v σ -2 h -WW T σ -2 h σ -2 v σ -2 h -WW T   . (4.43)
In this mean-field picture, the values of visible and hidden nodes are identified with the magnetizations. This consists in treating a RBM with Gaussian variables, and further assuming the external fields and the mean of the data to be null (normalization and rescaling of the training set can determine such conditions) we can rewrite the empirical expectation in (4.27) as

s α σ β data = σ 2 h w β s α σ β data = σ 2 h w β cov(s α , s β ) (4.44)
where we see that the covariance matrix of the data comes out. The deterministic learning equation of the weight matrix (4.27) can thus be written explicitly: the empirical average is given by Equation (4.44), while the model average is given by the covariance matrix (4.43). The diagonal part of the equation then reads

dw α dt = σ 2 h w α s 2 α data - σ 2 v 1 -σ 2 v σ 2 h w 2 α . (4.45)
while for the off-diagonal part we retain the rotation generators Ω u,v αβ rewriting them as

Ω u,v αβ = (1 -δ αβ )σ 2 h w β -w α w α + w β ∓ w β + w α w α -w β s α s β Data . (4.46)
In the equation above we see that the rotations are null only when s α s β is diagonal, and this happens when the spectral modes are aligned to the principal

w 2 α = s 2 α data -σ 2 v σ 2 v σ 2 h s 2 α data s 2 α data > σ 2 v 0 s 2 α data < σ 2 v (4.47)
where we see how the evolution of the singular values is driven by the SVD modes of the training data. The strongest modes, those above the threshold σ 2 v , are selected and learned while the modes below threshold are damped. Equations (4.45)(4.46)(4.47) let us simply summarize the behavior of the RBM in the linear regime: first, the learning equations drive the spectral modes of the weight matrix to align to the principal components of the data; subsequently, the singular values of the aligned modes are selected and amplified to match the dataset. Figure 4.2 shows the behavior of a linear RBM, in agreement with our theoretical analysis.

Nonlinear regime

To derive the dynamics in the non-linear regime, we will make use of the meanfield results of Sections 4.

σ α Data = w α s α (1 -q α [s]) Data (4.48) s α σ α Data = s α w β s β (1 -q β [s]) Data (4.49)
where q[s] is the empirical counterpart of the q α variable (4.19), here defined as

q[s] def = dx e x 2 2 √ 2π dvp α (v) tanh 2 κ -1/4 σx + γ w γ s γ v γ . (4.50)
The above averages are relatively easy to compute for a given dataset, but it is important to follow the prescriptions of Section 4.4 and select for the u α and v α components an appropriate distribution that leads to a compositional phase. The model average . . . RBM is instead computed by averaging the stable mean-field solutions weighted by the corresponding free energies. The mean-field partition function is thus defined as

Z T herm def = ω e -Lf (m ω , mω ,q ω ,q ω ) (4.51)
where ω indexes all the solutions of the saddle point equations (4.14),(4.15). The model averages then read

s α RBM = 1 Z T herm ω e -Lf (m ω , mω ,q ω ,q ω ) mω α def = mα T herm (4.52) s α σ β RBM = 1 Z T herm ω e -Lf (m ω , mω ,q ω ,q ω ) mω α m ω β def = mα m β T herm (4.53)
and we can put together all of the above to write a set of nonlinear dynamical equations that we are able to solve numerically

1 L dw α dt = s α w α s α (1 -q α [s]) Data -mα w α mα (1 -q α ) RBM (4.54) dη dt = mα T herm -s α Data + β Ω u αβ η β (4.55) dθ α dt = w α mα (1 -q α ) T herm -w α s α (1 -q α [s]) Data + β Ω v αβ θ β .
(4.56)

Empirical dynamics

To experimentally analyze the mean-field dynamics of the RBM given by Equations (4.54)-(4.56) we introduce a simple synthetic dataset, composed by C clusters of data with low dimensionality d embedded in a high dimensional space represented by N binary variables s i = ±1, with N d. To generate the data, we select a set of d orthonormal random vectors b α as a basis and we fix the magnetizations

m c i = d α=1 m c α b α i (4.57)
where the m c α factors are drawn at random between [-1, 1] and normalized. The multimodal distribution associated to the dataset is

P (s) = C c=1 p c N i=1 e h c i si 2 cosh(h c i ) (4.58)
where p c is a probability weighting cluster c and h c i is determined from the magnetizations m c i = tanh h c i . Samples are generated by choosing a cluster according to p c and setting the visible variables to ±1 according to

p(s i = 1) = 1 1 + e -2h c i .
The generated dataset is used to compute the empirical averages and the Equations (4.54)-(4.56) are integrated numerically, choosing the Laplace distribution for u α and v α components to compute the q, q, q α terms. The resulting datadriven analytical dynamics are shown in observe how the evolution is analogous to the linear case in which the modes emerge from the bulk one by one, with the difference that while in the linear regime the modes were evolving independently, here we notice how the expressed modes are interacting, with lower modes exerting a repulsive pressure on the singular values above threshold. The number of fixed point solutions used to compute the model averages increases in steps, roughly doubling each time a new mode is expressed, and the learning trajectory on the phase diagram (Figure 4.4) clearly shows how the RBM is found in the paramagnetic phase at initialization and the learning dynamics drive the model towards the learning phase, supporting our expectations that the model works in such a phase at equilibrium.

Using the same synthetic dataset, we also performed a standard training of the RBM (using CDk, see Section 2.2) and observed the spectral dynamics shown in Figure 4.5. We see how the evolution of the singular values and the emergence of new fixed point solutions follow closely the analytical evolution described above, with the difference that the RBM training conserves a lower number of fixed point solutions. This is probably due to the fact that the biases are expanded over an incomplete basis in Equations (4.9)(4.10), meaning that we are neglecting a residual part perpendicular to the K-modes spectral basis; this represents a limitation of our theoretical analysis but the overall behavior is well described, meaning that our mean-field model is sound.

Finally, in Figure 4.6 we show the empirical dynamics over the MNIST dataset. We see that the qualitative behavior follows the analytical prescriptions, showing that the mean-field model is well adapted to real-world scenarios. The picture of the RBM model that emerges is now rather complete. We have seen how the learning equations drive the selection of a certain number of Chapter 5

Missing information

A key advantage of the RBM is that it makes it possible to efficiently compute the conditional probability of a set of variables given the other variables. This property is fundamental to define its learning strategy and it adapts naturally to the case of missing information, i.e. datasets in which the visible variables are not fully observed. In this context, of particular interest is the case of learning with missing information and missing labels. We will show that the RBM can easily be adapted to model the joint probability of features and labels of a dataset, and this modeling results beneficial in imputing missing values for both features and labels collaboratively.

In the next Section, we will introduce the problem more in detail and underline its relevance. We will then introduce our approach and discuss our results.

Multi-output learning with incomplete data

Modern machine learning techniques usually require large sets of fully observed and well labelled data for training, which are seldom available in real-world applications. Sometimes a random subset of features is absent (e.g. failed sensors of a monitoring system), sometimes the data are insufficiently annotated (limitation due to the difficulties of human annotation) meaning that a lot of labels are missing. In fact the most common situation is that we have to train a model with both missing features and labels. Machine learning models that are able to deal with missing observations and missing labels are therefore highly desired.

In our work, we will consider multi-class and multi-label learning; the former associates an input instance to one class of a finitely defined class set, while the latter allows one instance to be associated to multiple labels simultaneously. Though both the problems of data imputation and multi-output learning with semi-supervised labels have been discussed in previous works [START_REF] Sun | Multi-Label Learning with Weak Label[END_REF]14,52], the proposed models are generally not designed to handle both challenges concurrently. In our work we shall consider that both features and labels are missing 5. Missing information completely at random, meaning that the mask associated to a sample data is assumed to be stochastically independent of the data.

Recently, methods based on Deep Latent Variables Models (DLVM) have been proposed to deal with missing data. In [START_REF] Mattei | MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets[END_REF], the Variational Autoencoder [46] has been adapted to be trained with missing data and a sampling algorithm for data imputation is proposed. Other approaches based on Generative Adversarial Networks (GAN) [33] are proposed in [START_REF] Yoon | GAIN: Missing Data Imputation using Generative Adversarial Nets[END_REF] and [51]. Impressive results on image datasets are displayed for these models, at the price of a rather high model complexity and the need for a large training set. In addition these works are focused on features reconstruction, and additional specifications and fine-tuning would be necessary to be able to take partially observed labels into account. The models specifications are quite involved, and any new specificity of the dataset may increase both the cost and the difficulty in training.

We chose to address the missing data problem in a more economical and robust manner. We consider the simple architecture of the Restricted Boltzmann Machine and adapt it to the multi-output learning context (RBM-MO) with missing data. The RBM-MO method serves as a generative model which collaboratively learns the marginal distribution of features and label assignments of input data instances, despite the incomplete observations. Building on the ideas expressed in [START_REF] Nijman | Using Boltzmann Machines to Fill in Missing Values[END_REF] we adapt the approach to the PCD training procedure (Section 2.2) and we propose a mean-field imputation method. Convincing results are shown on various real-world datasets. The advantage of the RBM-MO model is that of providing a robust and flexible method to deal with missing data, with little additional complexity with respect to the classic RBM. Moreover, it works seamlessly with multi-class and multi-label tasks, providing a unified framework for multi-output learning.

Lossy-CDk

In this section, we consider an RBM model with a layer of binary hidden nodes h j = 0, 1. The visible nodes are indicated by v and they will be either binary (v i = 0, 1) or Gaussian, according to a prior distribution p prior as described in Section 2.4, determining the following probability distribution

P (v, h) = e -E(v,h) Z p prior (v). (5.1)
To model the presence of missing values in the dataset, we separate the visible nodes in a set O of observed values and a set M of missing values. The visible layer then separates into the observed nodes v o = {v i , i ∈ O} and the missing nodes v m = {v i , i ∈ M}. The probability distribution over the observed variables can be written as

P (v o ) = Z O Z (5.2)

Lossy-CDk

where Z O is given by a marginalization over the missing variables

Z O = i∈M p prior (v i )dv i ×e k∈V a k v k j 1 + exp k∈V w kj v k + b j , V = O ∪ M. (5.
3)

The log-likelihood gradient w.r.t. the weights takes the form

∂L(v) ∂w ij = I o (i)v i hj h j p(h j |v O ) data (5.4) + (1 -I o (i))
hj

dv i v i h j p(h j |v o ) data (5.5) -v i h j RBM (5.6)
where I o is the indicator function of the set of observed nodes O. We note that w.r.t. the standard case the empirical average is now split in two terms, and we need to integrate over the missing variables v i , i ∈ M to compute the conditional probability p(h j |v o ) that we need to perform the Gibbs sampling. This integration is akin to the computation of the negative term, meaning that we can estimate the full positive term with a similar strategy: we pin the observed visible nodes to their value and we iterate the Gibbs sampling equations (Algorithm 1) for a prescribed number of steps k. for all minibatches m do Positive term:

6: pin variables v i , i ∈ O to their observed value 7: initialize v i , i ∈ M randomly 8: sample h, v using p(v | h) and p(h | v) for k steps 9:
compute the positive terms in (5.4) and (5.5) Negative term:

Mean-field imputation

Given a trained RBM modeling an empirical data distribution, this can be used to impute missing values for samples belonging to the same data distribution. Similarly to the Lossy-CDk training, the idea is to fix the values of the observed features and sample the missing features; fixing the known variables has the effect of biasing the sampling procedure towards the correct equilibrium configuration and helps in speeding up convergence, given that the fraction of missing features is small enough. For high percentages of missing features, we might expect this biased sampling to provide degraded solutions, as the low information given by the observed features can be supplanted by the sampling noise. To mitigate this effect, we can average multiple imputations weighted by the learned distribution. Recalling the mean-field description of the RBM in Chapter 4, we expect the equilibrium configurations to be represented by a limited number of fixed point magnetizations. Averaging values imputed from the fixed point magnetizations correlated with the observed variables provides imputations that are weighted by the model distribution (assuming the free energy differences of the fixed point solutions to be negligible), giving more reliable solutions. A set of self-consistent equations to compute the mean-field imputations is obtained by writing the mean-field equations at lowest order for the marginal probabilities of visible variables m i and the marginal probabilities of hidden variables q j

m i =   j w ij q j + a i   σ 2 v , i / ∈ O (5.7) q j = σ i w ij m i + b j (5.8)
that we can iterate to convergence while keeping fixed the observed values v i , i ∈ O. To obtain the weighted imputations, we simply run the above equations N f times (take N f ∼ 10) starting from random initial conditions and average the results

mi = 1 N f N f n=1 m (n) i .
(5.9)

Multi-output classification with missing information

A rather flexible approach to introduce labels in the RBM model is to simply model the joint probability distribution of both features and labels of an annotated dataset [49]. This boils down to separating the set of visible nodes V into a set of features nodes V f and a set of label nodes V l and let the model learn the correct label structure and the correlations between labels and features. The flexibility of this approach lies in the fact that multi-class and multi-label scenarios can be handled seamlessly with the same model, testing the capacity 5.4. Multi-output classification with missing information of the RBM to properly learn the constraints on labels structure in the different cases.

The learning algorithm is thus unmodified, while the mean-field imputation includes the p i magnetizations for label nodes giving the following set of iterative equations

m i =   j w ij q j + a i   σ 2 v , i / ∈ O (5. 10
)
p i = σ   j w ij q j + a i   , i / ∈ O (5.11) q j = σ   i∈V f w ij m i + i∈V l w ij p i + b j   .
(5.12)

Qualitative results of the RBM learning with missing information are shown in Figure 5.1 for the MNIST dataset; quantitative results for classification with missing labels and missing features in the multiclass setting are reported in Table 5.1, again for the MNIST dataset. The accuracy score in the base case in which only 30% of both features and labels are missing comes close to the literature result for RBM classification with fully observed samples [49], providing a good sanity check for our proposed mean-field imputation strategy. In the extreme case in which 80% of features and labels are missing the model still provides meaningful results, testifying the robustness of the approach. Other datasets and experimental settings (the multilabel learning case in particular) are discussed in [27] (Reprint C). Table 5.1: Classification results for the MNIST dataset with missing information. q ml % is the percentage of missing labels, and q f ea % the percentage of missing features. The AUC score is the Area Under the ROC curve.

Relative gradient optimization

In Part I about the RBM we have seen how a generative model of the data can be expressed as an energy-based model. Here we proceed from the observation that fundamentally we introduced generative models as probability distributions that we can learn from a dataset and subsequently sample from (Section 1.8). In this Part we present a simple invertible neural network architecture to model empirical distributions, together with an efficient training algorithm. This work has been originally presented in [35] (Reprint D).

Invertible transformations of probability densities and Normalizing Flows

We proceed from the assumption that there exists a "ground truth" probability distribution that assigns a normalized density to each sample x of the dataset of interest. The samples are considered to be high dimensional, with dimension d 1, and our goal is to learn a parametrized and tractable approximation of the distribution P (x) = P (x; θ). The strategy we adopt here is to consider a transformation T that maps a tractable base distribution P s (s) into the data distribution P (x). That is, T transforms a set of samples s drawn according to P s into the dataset samples x x = T (s), s ∼ P s (s).

Restricting the transformation T to be invertible and differentiable (with a differentiable inverse too), we can rewrite the data distribution in terms of the s samples by introducing a change of variable ([10], Section 1.2.1). This introduces the Jacobian matrix of the inverse transformation, that we denote J T -1 , and let us write the data log-likelihood as

log P (x) = log P s (s) + log | det J T -1 (x)|, s = T -1 (x). (6.1)

Relative gradient optimization

which we can optimize by Maximum Likelihood

P (x) = P (x; θ * ) θ * = arg max θ E x∈Data [log P (x; θ)] .
Choosing a tractable distribution P s , the complexity of computing and optimizing the data distribution P (x) is encapsulated into the determinant of the Jacobian matrix, to which we will refer as the Jacobian term in the following. We remark that restricting the transformation T to be invertible and differentiable not only let us employ the change of variable (6.1), but it also gives us the freedom to choose a tractable distribution P s (s) without restricting the class of distributions P (x) that we can model. Indeed, under mild conditions over P s and P , Equation (6.1) can represent any distribution P (x) if T is a diffeomorphism, i.e. it is continuous and both T and T -1 are differentiable [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF]. On the other hand, requiring T to be a diffeomorphism entails that the transformation must preserve the topological properties of the input space, i.e. the spaces in which x and s are embedded will have the same topology.

Concretely, this comes with two major constraints: i) x and s have the same dimension d;

ii) the base distribution P s and the data distribution P (x) model the same number of disconnected modes.

Constraint ii) tells us that, for instance, to choose an appropriate base distribution P s we need to know the number of separate clusters that the dataset is separated into, which is generally unknown and it is not clear how to determine it. This problem might seam daunting, but in practice we can ignore it and let the data distribution P (x) connect the separate clusters with a negligible amount of density mass; in practice, the employment of a unimodal base distribution has been shown to be able to approximate as well as possible some complex input topologies ( [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF]25], Figure 6.3). Constraint i) poses instead some strong limitations on the approximation capacity of the model; for instance, it is known that for transformations realized through neural networks, the ability to arbitrarily choose the model's width is necessary to achieve universal approximation capacity for arbitrary functions [53]. Notwithstanding this limitation, we recall that models employing a diffeomorphic transformation T are universal approximators for distributions, which seams to suggest that such models are powerful enough.

Another crucial property of diffeomorphic transformations T is that they are naturally composable and they present a simple expression for the Jacobian term in Equation (6.1)

T = T 1 • T 2 =⇒ det J T -1 1 •T -1 2 = det J T -1 1 • det J T -1 2 .
This let us build complex transformations of the base samples s by chaining together simpler transformations in successive steps, as exemplified in Figure 6.1. This idea is at the base of the class of models known as Normalizing Flows (NF), where the name comes from the fact that we can consider the inverse transformation taking the data samples and "normalizing" them into samples of the base distribution through a "flow" of transformations.

A comprehensive account of NF models is found in [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF]. The unifying aspect is that all NF models propose some kind of neural network architecture to parametrize the transformation T and optimize the log-likelihood (6.1) through backpropagation. The differences lie in the strategy that each model proposes to deal with the Jacobian term, which is generally expensive to compute. Indeed, being the Jacobian a d × d matrix, its determinant can be computed in O(d 3 ) time and the optimization with backpropagation also takes O(d 3 ) times. For high-dimensional datasets, it becomes quickly unfeasible to directly optimize the Jacobian term; NF models solve this problem by imposing constraints on the structure of the Jacobian matrix or by computing the Jacobian term in an approximate way, to make the optimization of the Jacobian term feasible.

Linear flows as Fully Connected neural network layers

One of the simplest transformations that we can consider in the context of Equation (6.1) is an invertible linear transformation

s = W x (6.2)
where W is an invertible d × d weight matrix parametrizing the inverse transformation T -1 . In the literature, the above transformation is known as a linear flow [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF] and its Jacobian matrix is simply the weight matrix W itself, implying that the Jacobian term is given by the determinant of W . Computation of the determinant is an O(d 3 ) operation, as well as the computation of its derivative

∂ log |det W | ∂W = W T -1 (6.3)
where the O(d 3 ) time complexity comes from the need of inverting the matrix W T . Maximum likelihood optimization of a linear flow with backpropagation

Relative gradient optimization

and Stochastic Gradient Descent (SGD) is thus a O(d 3 ) process, as backpropagation essentially computes the gradient (6.3) with the chain rule.

A more general nonlinear transformation can be obtained by identifying a linear flow as a Fully Connected (FC) layer of a multilayer neural network. Denoting by σ the activation function of the network, the k-th layer g k is expressed as

g k (x) = σ (W k x) . (6.4)
Choosing σ to be invertible and differentiable, we obtain the diffeomorphic transformation

T -1 k = g k . The composition of multiple FC blocks let us define a multilayer FC neural network g = g 1 • g 2 • • • • • g L with z k = g k (z k-1 ), k = 1, . . . , L (6.5)
where L is the number of layers. The input to the network will be the data samples x and the output represents the base samples s, so that we have z 0 = x and z L = s. Finally, for the base distribution we choose a standard multivariate normal with identity covariance, meaning that P s factorizes over the components of s and we can write the loglikelihood of the multilayer transformation as

log P (x) = i log N ([g(x)] i ; 0, 1) + log |det J g (x)| (6.6) = i log N ([g(x)] i ; 0, 1) + L k=1 log |det J g k (z k-1 )|
where we stress again that the computational complexity in estimating the probability density of the data lies in the computation of the Jacobian term; indeed, a forward pass through the network g has O(d 2 ) complexity, which becomes negligible w.r.t. the O(d 3 ) complexity of computing the Jacobian for large d. We note however that the computation of the Jacobian term has to be performed only once, and its result can be reused to compute the sample densities; in practice this mitigates the problem, as we can perform a single O(d 3 ) operation in advance and cache it to subsequently compute sample densities in O(d 2 ) time. Unfortunately this strategy doesn't work in the Maximum Likelihood optimization phase, in which we have to compute the gradients (6.3) at every layer for each step of the SGD update.

To obtain a better scaling behavior and make the training of linear flows feasible, we can parametrize the invertible matrix W in specific ways. [24,47] proposed to compute the P LU decomposition of W , where P is a permutation matrix (with unitary determinant) and L, U are lower-and upper-triangular matrices, whose determinant is the product of the diagonal elements. Computation of the Jacobian term is thus O(d), and the overall optimization is O(d 2 ). While this parametrization is in principle able to represent any invertible matrix W [START_REF] Papamakarios | Normalizing Flows for Probabilistic Modeling and Inference[END_REF], in practice we cannot optimize the permutation matrix P , thus limiting the class of matrices that we can represent. A more flexible 6.3. Relative gradients alternative is to consider the QR decomposition of W , where Q is an orthogonal matrix and R is upper triangular. Computing Q in full generality requires O(d 3 ) operations, but [START_REF] Jakub | Improving Variational Auto-Encoders using Householder Flow[END_REF] showed that we can apply the Q transformation as a sequence of at most d symmetry transformations each taking linear time. This makes it possible to compute and optimize the QR decomposition of W in O(d 2 ) time; note however that the sequential nature of the computation makes the method unsuitable for parallel optimization, making it rather slow in practice. An experimental comparison of the performance of the P LU and QR decompositions against the direct optimization of W is found in [41].

Relative gradients

The gradient of a function is generally introduced with the implicit assumption that the objects to optimize are defined w.r.t. a Euclidean space, i.e. a n-dimensional vector space in which the notions of distance and "orientation" (angle) are well defined. Computations are thus performed in the usual Cartesian coordinate system over the real space1 R n . In this context, invertible d × d matrices can be represented in the R d 2 vector space, i.e. a real space of dimension d 2 . This is possible due to the isomorphism between the space of invertible matrices R d×d = R d ⊗ R d , where ⊗ is the Kronecker product, and the real space R d 2 . However, while the isomorphism preserves the algebraic structure thus letting us perform valid computations, it doesn't guarantee that the geometric notions (such as distance) are preserved in the transformation. In the general case, invertible matrices live on a space with its own specific metric, i.e. a specification of the geometric notions of distance, angle, curvature and others. The Euclidean space has a flat metric; in dimension 3, this corresponds to the intuitive and evident properties of the physical world that we perceive, i.e. "naturally flat" surfaces and the interior angles of a triangle summing up to 180 • . The space of invertible matrices is a Riemannian manifold, a smooth space that is locally equivalent to a Euclidean space but with a different metric. Riemannian manifolds are generally thought of as curved spaces; the prototypical example is the surface of a sphere, which is naturally a 2-dimensional curved space. It is well known that it is not possible to project this surface on a 2-dimensional plane without distorting distances, which motivated the use of different maps and projections to represent the Earth surface over the centuries.

In the context of optimization, the use of the ordinary gradient defined in the Euclidean space becomes problematic when it is used with objects defined on a Riemannian manifold, as it doesn't respect the intrinsic curvature properties of such spaces. For invertible matrices, [4] showed how the steepest direction is not represented by the ordinary gradient but by the Natural Gradient (NG), which is introduced as the gradient computed directly on the 6. Relative gradient optimization manifold and not on its projection in Euclidean space. In the context of blind source separation, [16] derived a Relative Gradient (RG) for invertible matrices that stems from the properties of the multiplicative form of linear transformations (of the same kind as Equation (6.2)) and enjoys desirable convergence properties. As it turns out, the NG is identical to the RG as far as invertible matrices are concerned, with the NG being applicable in more general terms to any smooth manifold [15,[START_REF] Martens | New Insights and Perspectives on the Natural Gradient Method[END_REF].

For the layers of the FC network defined in Section 6.2, instead, the RG is not equivalent to the NG; consequently, in this context the RG doesn't represent the steepest direction and it doesn't enjoy the theoretical guarantees of learning efficiency that the NG does [4]. Nevertheless, the same properties that motivated the introduction of the RG are satisfied for the layers of the FC network we introduced, making it an attractive alternative to the ordinary gradient. We follow [16,35] to introduce the RG in a simple way and show that it provides a valid descent direction for a linear transformation.

A general way to compute the gradient of a function f is to perturb its input and consider the first-order term of the resulting Taylor expansion. When the input to the function is a matrix W , the perturbation will be represented by a matrix E with infinitesimal entries. The ordinary gradient ∇f (W ) is computed by adding the perturbation to the input

W → W + E (6.7)
and computing the Taylor expansion

f (W + E) -f (W ) = ∇f (W ), E + o(E) (6.8)
where A, B = ij a ij b ij is the Euclidean scalar product of matrices. The corresponding SGD update rule is obtained by aligning the update opposite to the gradient, i.e. setting E = -λ∇f (W ) to get

W → W -λ∇f (W ) (6.9)
where λ is the step size (corresponding to the learning rate in the context of neural networks training). The relative gradient ∇f rel (W ) is obtained by considering a multiplicative perturbation instead

W → W + EW = (I + E) W (6. 10
)
where I is the identity matrix and the Taylor expansion reads

f ((I + E)W ) -f (W ) = ∇f (W ), EW + o(E) = ∇f (W )W T , E + o(E) (6.
11) from which we obtain the equation defining the relative gradient2 ∇f rel (W ) = ∇f (W )W T . (6.12)

Relative backpropagation

Considering again a perturbation opposite to the gradient E = -λ∇f rel (W ) and applying it in a multiplicative fashion, we obtain the relative update rule

W → W -λ∇f (W )W T W . (6.13)
This update rule is valid, in the sense that it determines a decrease in the value of f . This is easily seen by taking λ small and computing the value of f for the updated value of W

f ((I -λ∇f rel (W ))W ) = f (W -λ∇f rel (W )W ) = f (W ) -λ ∇f (W ), ∇f rel (W )W + o(λ) = f (W ) -λ ∇f (W ), ∇f (W )W T W + o(λ) f (W ) -λ ∇f (W )W T , ∇f (W )W T = f (W ) -λ ∇f rel (W ) 2 < f (W )
where • 2 is the Frobenius norm.

Relative backpropagation

The crucial advantage of the RG is that it let us simplify the expression of the derivative of the Jacobian term for linear flows. We will show how this makes it possible to optimize linear flows efficiently and in full generality, i.e. without imposing constraints on the structure of the weight matrix W as it is usually done in NF literature. Subsequently, we will consider a more complex nonlinear model (the multilayer FC network introduced in Section 6.2) and show how the RG can play nicely with backpropagation, letting us introduce a Relative Backpropagation algorithm that makes it possible to efficiently optimize models employing unconstrained linear flows, with a cost that is equivalent to ordinary backpropagation. Finally, we will show some experimental results to empirically test the computational and statistical efficiency of the proposed algorithm.

Optimizing Linear Flows. Taking a linear flow f (x) = W x, we recall from Section 6.2 that the Jacobian term D f is given by

D f (x) = log |det J f (x)| = log |det W |
and its gradient is

∇D f (W ) = W T -1 . (6.14)
Exploiting Equation (6.12) we see that the RG is given by the identity matrix

∇D f rel = ∇D f W T = I (6.15)
which is a remarkable result as it let us avoid the computation of the gradient altogether and directly use the very simple update rule

W → (I -λ) W (6. 16
)
where λ is the step size. This makes the optimization of the Jacobian term of linear flows computationally efficient, as we avoid the matrix inversion. Moreover, we stress the fact that this optimization strategy can be applied in full generality without imposing any constraints on the structure of the weight matrix W , in stark contrast to the NF literature.

Backpropagating relative gradients. We consider the model g introduced in Section 6.2 and defined by Equations (6.4)(6.5). Defining the two objectives

L p (x) = i log N ([g(x)] i ; 0, 1) L J (x) = log |det J (x)|
where L J is the Jacobian term of the full network, the loglikelihood (6.6) can be written as L(x) = L p (x) + L J (x). (6.17)

We further split the Jacobian term over the layers and separate the linear transformations from the nonlinearities to get

L J (x) = L k=1 L k (z k ) L k (z k ) = L k,1 (z k ) + L k,2 (z k ) L k,1 (z k ) = d i=1 log |[σ(y k )] i | L k,2 (z k ) = log |det W k |
where y k = W k z k-1 and L k,2 is the Jacobian determinant of a linear flow.

Optimizing the full objective L(x) in an efficient manner is not straightforward. Indeed, with backpropagation we can calculate the ordinary gradients of each term of the loglikelihood, with the problem that the L k,2 term would be expensive to compute (Section 6.2). Substituting the ordinary gradient with the relative gradient using Equation (6.12) we can optimize L k,2 efficiently with the update rule (6.16), but the update rules for the other terms are

W k → (I -λ∇f (W k )W T k )W k ∼ O(d 3 ), f = {L p , L k,1 }
which are in turn expensive to compute due to the O(d 3 ) matrix multiplications, as both ∇f (W k ) and W k are d × d matrices. However, this problematic is only apparent and using some care it is possible to efficiently compute the relative gradient updates for all the terms of the loglikelihood. To see this, we use the chain rule to compute the gradient of the terms f = {L p , L k,1 } w.r.t. the parameters of layer k, i.e. the weight matrix W k :

∂f (x) ∂W k = ∂y k ∂W k ∂f (x) ∂y k = z k-1 δ T k δ k = ∂f (x) ∂y k T
where δ k is the error vector3 that is computed with automatic differentiation in the backpropagation algorithm. The relative gradient update rule then reads

W k → W k -λz k-1 ((δ T k W T k )W k ) ∼ O(d 2 ) (6.18)
which can be computed in O(d 2 ) time by carefully avoiding the matrix-matix multiplications (parentheses are important in Equation (6.18)). The full objective L(x) can thus be efficiently optimized in O(d 2 ) time with a hybrid approach: the L k,2 term is directly optimized with the update rule (6.16), while for L p and L k,1 we first compute the backpropagated errors δ k with automatic differentiation and then apply the update rule (6.18).

Experimental results. To test the computational efficiency and the empirical capacity of the model proposed in Section 6.2, we performed a couple of experiments:

i) a benchmark comparing optimization through the ordinary and relative gradients;

ii) probability density estimation over some simple 2-dimensional toy datasets.

In both cases, we used a model

g = g 1 • g 2 • • • • • g L with L layers comprising
the composition of a linear flow (6.2) and an invertible nonlinearity σ k defined as a smooth alternative to the Leaky ReLU function [54] commonly used in neural network architectures

σ k (x) = αx + (1 -α) log(1 + e x )
where α is a leakage coefficient that we treat as a hyperparameter. We note that the above expression is not invertible in closed form; in practice, this doesn't represent a problem as we can invert it numerically with the Newton method in a fixed number of iterations. For the last network layer we set σ L (x) = x, in an attempt to help the network to more easily map the outputs to a centered Normal distribution (the motivation comes from noticing that using the smooth Leaky ReLU defined above, we have σ k (0) > 0). The results of experiment i) are summarized in Figure 6.2. Even though the theoretical asymptotic difference between O(d 3 ) complexity for the ordinary gradient and O(d 2 ) complexity of the relative gradient is not clearly identified, the relative gradient optimization shows to be two order of magnitudes faster at both low and high dimensionality. We report benchmarking results for dimensionality up to d ∼ 10 4 ; for higher dimensionality, the real bottleneck becomes the number of parameters required by a FC layer, and the associated memory requirements.

Experiments ii) consists in training our model over data sampled from some simple 2-dimensional distributions and visualizing the learned probability densities as a heatmap. In Figure 6.3 we can see that the model achieves very of the gradient of the log-likelihood over a batch of 100 samples; the standard error of the mean is not reported as it is orders of magnitude smaller than the scale of the plot. We set the number of layers to L = 2 and performed the experiment using a Tesla P100 Nvidia GPU. Reproduced from [35].

good results in this simple setting. We remark that the original distributions to model are composed by disconnected clusters, and we know from Section 6.2 that these cannot be mapped to a unimodal Normal distribution by smooth invertible transformations. In spite of this, we see how the model is able to seemingly connect the cluster with a negligible amount of density mass.

Further experimental results over real-world datasets are reported in the original work ([35], Reprint D).

Conclusions

In Part I we have presented an extensive empirical analysis of the RBM during training, with particular focus on the dynamics of the SVD components of the weight matrix. We observed how the learning procedure selects a finite number of spectral modes to amplify, and a linear analysis of this process elucidates the role of the data in driving this selection. These observations motivated the choice of a statistical ensemble for the weight matrix of the RBM composed of a structured part and a random part. Mean-field analysis of the chosen ensemble let us derive a set of deterministic equations to describe both the dynamical evolution of the RBM during training and its equilibrium properties. The main outcome of this mean-field analysis is a clustering interpretation of the learning process, showing how the RBM is able to cluster the data around the solutions of its mean-field equations. While the mean-field equations showed some shortcomings in properly reproducing this clustering, a trained RBM fit this interpretation effectively.

When dealing with missing information, the simple formulation of the RBM naturally adapts to imputation problems. The clustering interpretation of the learning procedure motivated an imputation strategy in which multiple meanfield fixed-point solutions are averaged to obtain more robust results. This imputation algorithm is economical and successful in practice. A distinctive feature is its applicability to a wide range of different datasets including moderately low-dimensional data with only few hundred samples in the training set, as exemplified by an Internet-of-Things dataset discussed in [27] (Reprint C). Data-hungry deep learning models would most probably struggle in such a scenario. These considerations suggest that the proposed RBM-based imputation algorithm is a convenient alternative in real-world scenarios.

A couple of remarks, which are not discussed in this thesis, provide interesting avenues for further research. First, we note that a characterization of the free energy landscape of the RBM could in principle provide more refined imputations. The rationale is that the free energy of the mean-field fixed-point solutions can be used to weight the imputations over which we average; for simplicity, we disregarded this aspect in our work. Second, splitting the data features in observed and missing variables sets is akin to introducing a second layer of hidden units in the RBM architecture. From this point of view, with fully observed data and a RBM with 2 layers of hidden units the mean field imputation strategy can be analyzed as a training algorithm in its own right.

Relative gradient optimization

In Part II we have introduced the relative gradient and shown how it can be used to train Linear Flows without trading off expressive capacity for computational efficiency. We further proposed the use of Linear Flows as building blocks of a nonlinear multilayer neural network, and shown how the relative gradient plays nicely with standard backpropagation to provide an efficient training algorithm for the full architecture. The novelty of the proposed model is that it can scale to very high dimensionality without the need to impose structural constraints in the architecture or employing approximations during training, providing a new flexible model for density estimation. This feature makes it interesting from the theoretical point of view as it suggests a broader field of applications. In particular, the common "Gaussianization" operated by Normalizing Flows is conceptually close to Independent Component Analysis (ICA). In our work we took inspiration from the ICA literature in which the relative gradient is widely applied; inverting point of view, an interesting research proposition is to study how and if our model can be applied to perform ICA.

To conclude we note how the main subjects of this thesis have been combined before [34], suggesting that the investigation of the interplay between Restricted Boltzmann Machines and Normalizing Flows can be a fruitful one.

Introduction. -A Restricted Boltzmann machine (RBM) [1] constitutes nowadays a common tool on the shelf of machine learning practitioners. It is a generative model, in the sense that it defines a probability distribution, which can be learned to approximate any distribution of data points living in some N -dimensional space, with N potentially large. It also often constitutes a building block of more complex neural network models [2,3]. The standard learning procedure called contrastive divergence [4] is well documented [5] although being still a not so well understood fine empirical art, with many hyperparameters to tune without much guidelines. At the same time an RBM can be regarded as a statistical physics model, being defined as a Boltzmann distribution with pairwise interactions on a bipartite graph. Similar models have been already the subject of many studies in the 80's [6][7][8][9] which mainly concentrated on the learning capacity, i.e. the number of independent patterns that could be stored in such a model. The second life of neural networks has renewed the interest of statistical physicists for such models. Recent works actually propose to exploit its statistical physics formulation to define mean-field based learning methods using TAP equations [10][11][12]. Meanwhile some analysis of its static properties, assuming a given learned weight matrix W , have been proposed [13,14] in order to understand collective phenomena in the latent representa-(a) E-mail: aurelien.decelle@lri.fr tion [15], i.e. the way latent variables organize themselves to represent actual data. One common assumption made in these works is that the weights of W are i.i.d. which as we shall see is unrealistic. Concerning the learning procedure of neural networks, many recent statistical physics based analysis have been proposed, most of them within teacher-student setting [16] which imposes a strong assumption on the data, namely that these are generated from a model belonging to the parametric family of interest, hiding as a consequence the role played by the data themselves in the procedure. From the analysis of related models [17,18], it is already a well established fact that a selection of the most important modes of the singular value decomposition (SVD) of the data is performed in the linear case. In fact in the simpler context of linear feed-forward models the learning dynamics can be fully characterized by means of the SVD of the data matrix [19], showing in particular the emergence of each mode by order of importance regarding singular values.

In this work we follow this guideline in the context of a general RBM. We propose to characterize both the learned RBM and the learning process itself by the SVD spectrum of the weight matrix in order to isolate the information content of an RBM. This allows us then to write a deterministic learning equation leaving aside the fluctuations. This equation is subsequently analyzed first in the linear regime to identify the unstable deformation modes of W ; secondly at equilibrium assuming the learning is convergp-1 arXiv:1708.02917v2 [cond-mat.dis-nn] 30 Nov 2017 A. Spectral dynamics of learning in restricted Boltzmann machines ing, in order to understand the nature of the non-linear interactions between these modes and how these are determined from the input data. In the first section we recall the RBM model and associated learning algorithm.

In the second section we show how this algorithm can be described by a generic learning equation. Then we first analyze the linear regime and thereafter we describe what happens with the binary RBM. A set of dynamical parameters is shown to emerge naturally from the SVD decomposition of the weight matrix. The convergence toward equilibrium is analyzed and illustrated later with actual tests on the MNIST dataset.

The RBM and associated learning procedure. -An RBM is a Markov random field with pairwise interactions defined on a bipartite graph formed by two layers of non-interacting variables: the visible nodes and the hidden nodes representing respectively data configurations and latent representations. The former noted s = {s i , i = 1 . . . N v } correspond to explicit representations of the data while the latter noted σ = {σ j , j = 1 . . . N h } are there to build arbitrary dependencies among the visible units. They play the role of an interacting field among visible nodes. Usually the nodes are binary-valued (of boolean type or Bernoulli distributed) but gaussian distributions or more broadly arbitrary distributions on real-valued bounded support are also used [20], ultimately making RBMs adapt for more heterogeneous data sets.

Here to simplify we assume that visible and hidden nodes will be taken as binary variables s i , σ j ∈ {-1, 1} (using ±1 values has the advantage of symmetrizing the equations hence avoiding to deal with "hidden" biases on the variables when considering binary {0, 1} variables). Like the Hopfield model [6] which can actually be cast into an RBM [21] an energy function is defined for a configuration of nodes

E(s, σ) = - i,j s i w ij σ j - Nv i=1 η i s i - Nh j=1 θ j σ j (1)
and this is exploited to define a joint distribution between visible and hidden units, namely the Boltzmann distribution

p(s, σ) = e -E(s,σ) Z ( 2
)
where W is the weight matrix and η and θ are biases, or external fields on the variables. Z = s,σ e -E(s,σ) is the partition function of the system. The joint distribution between visible variables is then obtained by summing over hidden ones. In this context, learning the parameters of the RBM means that, given a dataset of M samples composed of N v variables, we ought to infer values to W , η and θ such that new generated data obtained by sampling this distribution should be similar to the input data. The general method to infer the parameters is to maximize the likelihood of the model, where the pdf (2) has first been summed over the hidden variables

L = j log(2 cosh( i w ij s i + θ j )) -log(Z). (3)
Different methods of learning have been set up and proven to work efficiently, in particular the contrastive divergence (CD) algorithm from Hinton [4] and more recently TAP based learning [10]. They all correspond to expressing the gradient ascent on the likelihood as

∆w ij = γ ( s i σ j p(σ j |s) Data -s i σ j pRBM ) ( 4
)
where γ is the learning rate. Similar equations can be derived for the biases. The main problem is the second term on the rhs of (4) which is not tractable, and various methods basically differ in their way of estimating this term (Monte-Carlo chains, mean field, TAP . . . ). For an efficient learning the first term also has to be approximated by making use of random mini batches of data at each step.

Deterministic dynamics of the learning. -In order to understand the dynamics of the learning we first project the CD equation ( 4) onto the basis defined by the SVD of W . As a generalization of eigenmodes decomposition to rectangular matrices, the SVD for a RBM is given by

W = UΣV T ( 5
)
where U is an orthogonal N v ×N h matrix whose columns are the left singular vectors u α , V is an orthogonal N h ×N h matrix whose columns are the right singular vectors v α and Σ is a diagonal matrix whose elements are the singular values w α . The separation into left and right singular vectors is due to the rectangular nature of the decomposed matrix, and the similarity with eigenmodes decomposition is revealed by the following SVD equations

Wv α = w α u α W T u α = w α v α
We consider the usual situation where N h < N v , which means that the rank of W is at most N h . W (t) represents the learned weight matrix at time t.

Let {w α (t) ∈ [0, +∞[}, {u α (t) ∈ R Nv } and {v α (t) ∈ R Nh } such that the following decomposition w ij (t) = α u α i (t)w α (t)v α j (t) holds. Discarding stochastic fluctuations usually inherent to the learning procedure and letting the learning rate γ → 0, the continuous version of (4) can be recast as follows:

p-2 dw dt αβ = δ α,β dw α dt (t) + (1 -δ α,β ) w β (t)Ω v βα (t) + w α (t)Ω h αβ = s α σ β Data -s α σ β RBM (6) Ω v αβ (t) = -Ω v βα def = du α,T dt u β = -1 w α + w β dw dt A αβ + 1 w α -w β dw dt S αβ (7) Ω h αβ (t) = -Ω h βα def = dv α,T dt v β = 1 w α + w β dw dt A αβ + 1 w α -w β dw dt S αβ (8)
Here everything is expressed in the reference frame defined by singular vectors of W . s α = i u α i s i and σ α = j v α j σ j represent spin configurations in this frame. Note that one has to keep track of the original reference frame to be able to evaluate the data and RBM average in particular when the basic variables are discrete. We have introduced the skew-symmetric rotation generators Ω v,h αβ (t) of the basis vectors induced by the dynamics. These tell us how the data rotate relatively to this frame. The superscript S,A indicate the symmetric (resp. antisymmetric) part of the matrix. Note that these equations become singular when some degeneracy occurs in W because then the SVD is not uniquely defined. This is not really a problem since we are interested in rotations among non-degenerate modes, the rest corresponding to gauge degrees of freedom. Similar equations can be derived for the fields η α (t)

def = i η i (t)u α i (t) and θ α (t) def = j v α j (
t)θ j (t) projected onto the SVD modes. At this point we make the assumption that the learning dynamics is represented by a trajectory of ({w α (t), η α (t), θ α (t), Ω v,h αβ (t)}, while the specific realization of the u α i and v α j is considered to be irrelevant, and can be averaged out with respect to some simple distributions, as long as this average is correlated with the data. This means that the decomposition ŝα = i u α i ŝi of any given sample configuration is assumed also to be kept fixed while averaging. What matters mainly is the strength given by w α (t) and the rotation given by Ω v,h αβ (t) of these SVD modes. Assuming for example i.i.d centered normal distribution with respective variance 1/N v and 1/N h for u α i and v α j , the empirical term takes the simple form:

s α σ β Data = 1 N h s α (s β w β -θ β )V 1 N h γ (w γ s γ -θ γ ) 2 Data where V (x) = dy e -y 2 /2 √ 2π sech 2 ( √ xy), (9)
which actually depends on the activation function (an hyperbolic tangent in this case). The main point here is that the empirical term defines an operator whose decomposition onto the SVD modes of W functionally depends solely on w α , θ α and on the projection of the data on the SVD modes of W . This term is precisely driving the dynamics. The adaptation of the RBM to this driving force is given by the second term which can be as well estimated in the thermodynamic limit, as a function of w α , θ α and η α alone.

Linear instabilities. -First let us consider the linear regime which can be analyzed thoroughly. It can be obtained by rescaling all the weights and fields by a common "inverse temperature" β factor and let this go to zero in equations ( 6). This limit can be understood by keeping up to quadratic terms in the mean field free energy and should correspond to the first stages of the learning. In this limit, magnetizations (µ v , µ h ) of visible and hidden variables have Gaussian fluctuations with covariance matrix

C(µ v , µ h ) def = σ -2 v -W -W T σ -2 h -1
with σ 2 v = σ 2 h = 1 introduced for sake of generality when considering general linear RBM. To simplify the exposition, we discard the biases of the data and related fields (θ α , η α ) of the RBM. In that case the empirical term in (6) involves directly the covariance matrix of the data expressed in the frame defined by the SVD modes of W s α σ β Data = σ 2 h w β s α s β Data .

From C(µ v , µ h ) we get the other terms yielding the following equations:

dw α dt = w α σ 2 h s 2 α Data - σ 2 v 1 -σ 2 v σ 2 h w 2 α Ω v,h αβ = (1 -δ αβ )σ 2 h w β -w α w α + w β ∓ w β + w α w α -w β s α s β Data
Note that these equations are exact for a linear RBM, since they can be derived without any reference to the coordinates of u α and v α over which we average in the non-linear regime. These equations tell us that, during the learning the vectors u α (and also v α ) will rotate until being aligned to the the principal components of the data, i.e. until s α s β Data becomes diagonal. Then calling ŵ2 α the corresponding empirical variance given by the data, 1: Time evolution of the eigenvalues in the linear model and of the likelihood. We observe very clearly how the different modes emerge from the bulk and how the likelihood increases at each eigenvalue learned. In the inset, the scalar product of the vectors u obtained from the SVD of the data and of w. The us of w are aligned with the SVD of the data at the end of the learning.

the system reach the following equilibrium values:

w 2 α =    ŵ2 α -σ 2 v σ 2 v σ 2 h ŵ2 α if ŵ2 α > σ 2 v , 0 if ŵ2 α ≤ σ 2 v .
From this we see that the RBM selects the strongest SVD modes in the data. The linear instabilities correspond to directions for which the variance of the data is above the threshold σ 2 v . This determines the deformations of the weight matrix which can develop during the learning and will eventually interact, following the usual mechanism of non-linear pattern formation like e.g. in reaction-diffusion processes [22]. Other possible deformations are damped to zero. The linear RBM will therefore learn all (up to N h ) principal components that passed the threshold but it is important to remember that the resulting distribution will still be unimodal. Note that this selection mechanism is already known to occur for linear auto-encoders [18] or some other similar linear Boltzmann machines [17]. On Fig. 1 we can see the eigenvalues being learned one by one in a linear RBM. For non-linear RBM when the system escapes the linear regime, a well suited mean-field theory is required to understand the dynamics and the steadystate regime.

Non-linear regime. -During the linear regime some specific modes are selected and at some point these modes start to interact in a non-trivial manner. The empirical term in (6) involves higher order statistics of the data as exemplified by (9) and the Gaussian estimation with σ 2 v = σ 2 h = 1 of the RBM response term s α σ β RBM is no longer valid. In order to estimate this term in the thermodynamic limit, some assumptions on the form of the weight matrix are needed. A common assumption consists in considering i.i.d. random variables for the weights w ij and this, like for example in [13][14][15], generally leads to a Marchenko-Pastur distribution of the singular values of W , which as we shall see in the next section is unrealistic. Instead, based on our experiments such distribution corresponds to the noise of the weight matrix, while its information content is better expressed by the presence of SVD modes outside of the bulk. This leads us to write the weight matrix as

w ij = K α=1 w α u α i v α j + r ij (10)
where the w α = O(1) are isolated singular values (describing a rank K matrix), the u α and v α are the eigenvectors of the SVD decomposition and the r ij = N (0, σ 2 /L) where L = √ N h N v are i.i.d. corresponding to noise. To be consistent with the linear analysis, these modes are assumed to span the (left) subspace corresponding to the part of the empirical SVD above threshold while r spans the complementary space of empirical modes below threshold. We limit the analysis here to the case where K is finite. This then allows us to assume simple distributions p u and p v for the components of u α and v α considered i.i.d. for instance. This altogether defines our statistical ensemble of RBM to which we restrict ourselves to study the learning procedure. For K extensive we should instead average over the orthogonal group which would lead to a slightly different mean-field theory [23,24]. In the present form our model of RBM is similar to the Hopfield model and recent generalizations [25], the patterns being represented by the SVD modes outside the bulk. The main difference, in addition to the bipartite structure of the graph, is the non-degeneracy of the singular values w α . Still the analysis in the thermodynamic limit follows classical treatments like [7,26] for the Hopfield model or [14] for bipartite models. The starting point is to express the average over u, v and weights r ij of the log partition function Z in (2) with the help of the replica trick:

E u,v,r [log(Z)] = lim p→0 d dp E u,v,r [Z p ].
After averaging over the iid weights, 4 sets of order parameters {(m a α , ma α ), a = 1, . . . p, α = 1, . . . K} and {(Q ab , Qab ), a, b = 1, . . . p, a = b} are introduced with help of two distinct Hubbard-Stratonovich transformations. These variables represent the following quantities:

m a α ∼ 1 √ L E u,v,r σ a α ma α ∼ 1 √ L E u,v,r s a α Q ab ∼ E u,v,r σ a i σ b i Qab ∼ E u,v,r s a j s b j ,
namely the correlations of the hidden [resp. visible] states with the left [resp. right] singular vectors and the Edward-Anderson order parameters measuring the correlation between replicas of hidden or visible states. E u and E v denote an average wrt to the rescaled components u

√ N v u α i p-4
and v √ N h v α j of the SVD modes. The transformations involve pairs of complex integration variables because of the asymmetry introduced by the two-layers structure by contrast to fully connected models. They lead to the following representation:

E u,v,r [Z p ] = a,α dm a α d ma α 2π a =b dQ ab d Qab 2π exp -L a,α w α m α mα + σ 2 2 a =b Q ab Qab - 1 √ κ A[m, Q] - √ κB[ m, Q] with A[m, Q] def = log S a ∈{-1,1} E u e √ κσ 2 2 a =b QabS a S b +κ 1 4 a,α (m a α wα-ηα)u α S a , κ = N h /N v and B[ m, Q] obtained from A[m, Q
] by replacing u by v, η by θ and κ by 1/κ. The thermodynamic properties are obtained by first letting L → ∞ allowing for a saddle point approximation and then the limit p → 0 is taken. We restrict here the discussion to replica symmetric (RS) saddle points [27]. The breakdown of RS can actually be determined by computing the so-called AT line [28] and will be detailed somewhere else [29]. In the RS case the set {(Q ab , Qab } reduces to a pair (q, q) of spin glass parameters, while quenched magnetization towards the SVD directions are now represented by {(m α , mα ), α = 1, . . . K}. Letting x = N (0, 1) and skipping some details, the saddlepoint equations are given by

(m α , mα ) = E κ 1 4 v α tanh h(x, v) , κ -1 4 u α tanh h(x, u) (11)
(q, q) = E tanh 2 h(x, v) , tanh 2 h(x, u) , (12) with E denoting the average over (u, v, x) and

h(x, u) def = κ 1 4 σ √ qx + γ (w γ m γ -η γ )u γ h(x, v) def = κ -1 4 σ √ qx + γ (w γ mγ -θ γ )v γ .
These fixed point equations can be solved numerically to tell us how the variables condensate on the SVD modes within each equilibrium state of the distribution and whether a spin glass phase is present or not. The important point here is that with K finite and a nondegenerate spectrum the mode with highest singular value dominates the ferromagnetic phase. The phase diagram looks in fact similar to the one of the SK model with ferromagnetic coupling, when 1/σ is interpreted as a temperature and w max /σ the ferromagnetic coupling. Some subtleties arise when considering various ways of averaging over singular vectors components [29]. In [15,30] it is underlined the importance of the capability of networks to produce compositional states structured by combination of hidden variables. In our representation, we don't have direct access to this property, but to the dual one in some sense, namely states corresponding to combination of modes. Their presence and their structure, are rather sensitive to the way the average over u and v is performed. In this respect the case where u α and v α are Gaussian i.i.d distributed is very special: all other fixed points associated to lower modes can be shown to be unstable as well as fixed points associated to combinations of modes. Instead, for other distributions with smaller kurtosis, like uniform or Bernoulli, stable fixed points associated to many different single modes or combinations of modes can exist and contribute to the thermodynamics.

Coming back to the learning dynamics, the first thing which is expected, already from the linear analysis, is that the noise term in (10) vanishes by condensing into a delta function of zero modes. Then the term corresponding to the response of the RBM in ( 6) is estimated (in absence of bias) in the thermodynamic by means of the order parameters defined previously:

s α σ β RBM = L Z MF C q=1 e -Fq m(q) α m (q) β , Z M F def = q e -Fq
where the index q run over all stable fixed point solutions of (11,12) weighted accordingly to their free energy. These are the dominant contributions as long as free energy differences are O(1), internal fluctuations given by each fixed point are comparatively of order O(1/L). Note that this is the reason why the RBM needs to reach a ferromagnetic phase with many states to be able to match the empirical term in (6) in order to converge. For instance, in the case of a multimodal data distribution with many well separated clusters, the SVD modes of W which will develop are the one pointing in the direction of the magnetizations defined by these clusters. In this simple case the RBM will evolve as in the linear case to a state such that the empirical term becomes diagonal, while the singular values adjust themselves until matching the proper magnetization in each fixed point. More precise statements about the phase diagram of the RBM and the behaviour of our dynamical equations including the dynamics of the external fields η α and θ α will be given in [29].

Tests on the MNIST dataset. -We illustrate our results on the MNIST dataset. The MNIST dataset is composed of 60000 images of handwritten digits of 28 × 28 p-5 The hidden field is put to zero. The learning rate is chosen to be ≈ 0.01. With these parameters we verified that our machine was able to sample digits in a satisfactory way after 20 epochs. Now we can investigate the value of some observables introduced previously. First, we look at the SVD modes of the matrix w during the learning on Fig. 2. We see that, after seeing only few updates the system has already learned many SVD modes from the data. On Figure 2a-2c, we observe what is expected from the linear regime. Some modes escape from the Marchenko-Pastur bulk of the eigenvalues while other condense down to zero. In particular, we can see that the modes at the beginning of the learning correspond exactly to the SVD modes of the data, see Fig. 3. On this figure, we notice that the modes of the W matrix are the same as the ones of the data at the beginning of the learning as predicted by the linear theory.

A. Spectral dynamics of learning in restricted Boltzmann machines

After many epochs, we observe on Fig. 3-f that nonlinear effects have deformed the SVD modes of W with respect to the beginning of the learning. We can also look at the evolution of the eigenvalues of W . On Fig. 4 we observe their evolution and when they start to be amplified (or dumped). On the inset, we see how the strongest mode get out of the bulk and increase while the lowest ones are dumped after many epochs. We also observe that the top part of the spectrum of W appear flattened as compared to empirical SVD spectrum. This presumably favors the expression of many states of similar free energy related to various digit configurations, able to contribute to RBM response term in (6).

Discussion. -The equations obtained for the dynamics and the MF theory that allows us to compute them constitute a phenomenological description of the learning of an RBM. This is assumed to represent a typical learning trajectory in the limit of infinite batch size. These equations have been obtained by averaging over the components of left and right SVD vectors of the weight matrix, keeping fixed a certain number of quantities considered to be the relevant ones, fully characterizing a typical RBM during the learning process. This averaging corresponds actually to a standard self-averaging assumption in a RS phase. The singular values spectrum {w α } is playing the main role. The projections (η α , θ α ) of the bias onto the eigenmodes of W are also considered as intrinsic quantities. Finally the rotation vectors {Ω v,h α,β } give us the relative motion of the data w.r.t the time dependent frame given by the singular vectors of W . In our phenomenological description the learning dynamics is represented by a trajectory of {w α (t), η α (t), θ α (t), Ω v,h αβ (t)} which is uniquely determined by our equations once an initial condition specified by the decomposition of the data on the singular vectors of W is given. By contrast to usual approaches which rely on the teacher-student scenario, we may obtain generic learning curves of non-linear neural networks, which are driven by intrinsic properties of the data. The point is to give insights into the relationship between model and data. This allows us to give some elements of understanding on which properties of the data drive the learning and how they are represented in the model. Eventually this will lead us to identify and cure p-6 Fig. 4: Log-log plot of the singular values represented as discrete abscissas (in decreasing order) with their magnitude reported on the ordinates. The RBM contained 400 hidden variables. A cutoff is highlighted by the onset of the linear behaviour and the SVD modes of the data in black. We qualitatively observe that beyond some α tresh the modes are dumped while before they are amplified. In the inset, the time evolution of the modes 1, 2, 10, 100, 350, 400 during the learning as a function of the number of epochs, we see that for large value of α, the modes are decreasing. We observe that the linear cutoff (around α ≈ 50 seems different from the one observed when going deep into the non-linear regime (α ≈ 250). some flaws of present learning methods.

1 Introduction 2 connections between RBMs and well known disordered systems in statistical physics.

As an example, when data are composed by vectors with binary components the discrete RBM takes the form of an heterogeneous Ising model composed of one layer of visible units (the observable variables) connected to one layer of hidden units (the latent or hidden variables building up the dependencies between the visible ones), in which couplings and fields are obtained from the training data through a learning procedure. In order to build more powerful models, RBMs can be stacked to form "deep" architectures. In such a case, they can form a multi-layer generative model known as a Deep Boltzmann Machine (DBM) [2] or they can be stacked and trained layerwise as a pre-training procedure for neural networks [3]. The standard learning algorithms in use are the contrastive divergence [4] (CD) and the refined Persistence CD [5] (PCD), which are based on a quick Monte Carlo estimation of the response function of the RBM and are efficient and well documented [6]. Nevertheless, despite some interesting interpretations of CD in terms of non-equilibrium statistical physics [7], the learning of RBMs remains a set of obscure recipes from the statistical physics point of view: hyperparameters (like the size of the hidden layer) are supposed to be set empirically without any theoretical guidelines. Historically, statistical physics played a central role in studying the theoretical foundations of neural networks. In particular, during the 1980s many works on the Hopfield model [8,9,10,11] managed to define its learning capacity and to compute the number of independent patterns that it could store. It is worth noticing that, as RBMs are ultimately defined as a Boltzmann distribution with pairwise interactions on a bipartite graph, they can be studied in a way similar to that used for the Hopfield model. The analogy is even stronger since connections between the Hopfield model and RBMs have been made explicit when using Gaussian hidden variables [12], here the number of patterns of the Hopfield model corresponding to the number of hidden units. Motivated by a renewed excitement for neural networks, recent works actually propose to exploit the statistical physics formulation of the RBM to understand what is its learning capacity and how mean-field methods can be exploited to improve the model. In [13,14,15], mean-field based learning methods using TAP equations are developed. TAP solutions are usually expected to define a decomposition of the measure in terms of pure thermodynamical states and are useful both as an algorithm to compute the marginals of the variables of the model and to identify the pure states when they are yet unknown. For instance, in a sparse explicit Boltzmann machine (i.e. without latent variables) this implicit clustering can be done by means of belief propagation fixed points1 with simple empirical learning rules [16]. In [17,18], an analysis of the static properties of RBMs is done assuming a given weight matrix W , in order to understand collective phenomena in the latent representation, i.e. the way latent variables organize themselves in a compositional phase [19,20] to represent actual data. These analysis make use of the replica trick (or equivalent) making the common assumption that the components of the weight matrix W are i.i.d.; despite the fact that this approach may give some insights into the retrieval phase, this approximation is problematic since, as far as a realistic RBM is concerned (an RBM learned on data), the learning mechanism introduces correlations within the weights of W and then it seems rather crude to continue to assume the independence and hope to understand the realistic statistical properties of the model.

Concerning the learning procedure of neural networks, many recent statistical physics based analyses have been proposed, most of them within teacher-student set-ting [21]. This imposes a rather strong assumption on the data in the sense that it is assumed that these are generated from a model belonging to the parametric family of interest, hiding as a consequence the role played by the data themselves in the procedure. From the analysis of related linear models [22,23], it is already a well established fact that a selection of the most important modes of the singular values decomposition (SVD) of the data is performed in the linear case. In fact in the simpler context of linear feed-forward models the learning dynamics can be fully characterized by means of the SVD of the data matrix [24], showing in particular the emergence of each mode by order of importance with respect to the corresponding singular values.

First steps to follow this guideline have been done in [25], in the context of a general RBM and to address the shortcomings of previous analyses, in particular concerning the assumptions over the weights distribution. To this end it has been proposed to characterize both the learned RBM and the learning process itself by means of the SVD spectrum of the weight matrix in order to single out the information content of the RBM. It is assumed that the SVD spectrum is split in a continuous bulk of singular vectors corresponding to noise and a set of outliers that represent the information content. By doing this it is possible to go beyond the usual unrealistic assumption of i.i.d. weights made for analyzing RBMs. Proceeding along this direction, in the present work we first present a thermodynamic analysis of RBMs under the more realistic assumptions over the weight matrix that we propose. Then, on the same basis, the learning dynamics of RBMs is studied by direct analysis of the dynamics of the SVD modes, both in the linear and non-linear regimes.

s 1 s n v σ 1 σ n h s i σ j

Hidden layer

Visible layer

W ij

Fig. 1: bipartite structure of the RBM.

The paper is organized as follows: in Section 2 we introduce the RBM model and its associated learning procedures. Section 3 presents the static thermodynamical properties of the RBM with realistic hypothesis on its weights: a statistical ensemble of weight matrices is discussed in Section 3.1; mean-field equations in the replica-symmetric (RS) framework are given in Section 3.2 and the corresponding phase diagram is studied in Section 3.3 with a proper delimitation of the RS domain where the learning procedure is supposed to take place. The ferromagnetic phase is studied in great details in 3.4 by looking in particular at the conditions leading to a compositional phase. Section 4 is devoted to the learning dynamics. In Section 4.1, a deterministic learning equation is derived in the thermodynamic limit and a set of dynamical parameters is shown to emerge naturally from the SVD of the weight matrix. This equation is analyzed for linear RBMs in Section 4.2 in order to identify the unstable deformation modes of W that result in the first emerging patterns at the beginning of the learning process; the non-linear regime is described in Section 4.3, on the basis of the thermodynamic analysis, by numerically solving the effective learning equations in simple cases. Our analysis is finally illustrated and validated in Section 5 by actual tests on the MNIST 2 The RBM and its associated learning procedure An RBM is a Markov random field with pairwise interactions defined on a bipartite graph formed by two layers of non-interacting variables: the visible nodes and the hidden nodes representing respectively data configurations and latent representations (see Figure 1). The former noted s = {si, i = 1 . . . Nv} correspond to explicit representations of the data while the latter noted σ = {σj, j = 1 . . . Nh} are there to build arbitrary dependencies among the visible units. They play the role of an interacting field among visible nodes. Usually the nodes are binary-valued (of Boolean type or Bernoulli distributed) but Gaussian distributions or more broadly arbitrary distributions on real-valued bounded support are also used [26], ultimately making RBMs adapted to more heterogeneous data sets. Here to simplify we assume that visible and hidden nodes will be taken as binary variables si, σj ∈ {-1, 1} (using ±1 values gives the advantage of working with symmetric equations hence avoiding to deal with the "hidden" biases on the variables that appear when considering binary {0, 1} variables). Like in the Hopfield model [8], which can actually be cast into an RBM [12], an energy function is defined for a configuration of nodes

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

E(s, σ) = - i,j siWijσj + Nv i=1 ηisi + N h j=1 θjσj ( 1
)
and this is exploited to define a joint distribution between visible and hidden units, namely the Boltzmann distribution

p(s, σ) = e -E(s,σ) Z ( 2
)
where W is the weight matrix and η and θ are biases, or external fields on the variables. Z = s,σ e -E(s,σ) is the partition function of the system. The joint distribution between visible variables is then obtained by summing over hidden ones. In this context, learning the parameters of the RBM means that, given a dataset of M samples composed of Nv variables, we ought to infer values to W , η and θ such that new generated data obtained by sampling this distribution should be similar to the input data. The general method to infer the parameters is to maximize the log likelihood of the model, where the pdf (2) has first been summed over the hidden variables

L = j log(2 cosh( i Wijsi -θj)) Data - i ηi si Data -log(Z). (3)
Different learning methods have been set up and proven to work efficiently, in particular the contrastive divergence (CD) algorithm from Hinton [4] and more recently TAP based learning [13]. They all correspond to expressing the gradient ascent on the likelihood as

∆Wij = γ ( siσjp(σj|s) Data -siσj p RBM ) (4) ∆ηi = γ ( si p RBM -si Data) (5) ∆θj = γ ( σj p RBM -σjp(σj|s) Data) ( 6
)
where γ is the learning rate. The main problem are the • • • p RBM terms on the right hand side of (4-6). These are not tractable and the various methods basically differ in their way of estimating those terms (Monte-Carlo Markov chains, naive mean-field, TAP. . . ). For an efficient learning the • • • Data terms must also be approximated by making use of random mini-batches of data at each step.

3 Static thermodynamical properties of an RBM

Statistical ensemble of RBMs

When analyzing the thermodynamical properties of RBMs, it is common to assume that the weights Wij are i.i.d. random variables, like for example in [20,17,18]. This generally leads to a Marchenko-Pastur (MP) distribution [27] of the singular values of W , which is unrealistic.

In order to clarify our notation, let us recall the definition of the singular value decomposition (SVD). As a generalization of eigenmodes decomposition to rectangular matrices, the SVD for a RBM is given by

W = UΣV T ( 7
)
where U is an orthogonal Nv × Nh matrix whose columns are the left singular vectors u α , V is an orthogonal Nh×Nh matrix whose columns are the right singular vectors v α and Σ is a diagonal matrix whose elements are the singular values wα. The separation into left and right singular vectors is due to the rectangular nature of the decomposed matrix, and the similarity with eigenmodes decomposition is revealed by the following SVD equations

Wv α = wαu α W T u α = wαv α
In [25] it is argued that the MP distribution of SVD modes actually corresponds to the noise of the weight matrix, while the information content of the RBM is better expressed by the presence of SVD modes outside of this bulk. This leads us to write the weight matrix as

Wij = K α=1 wαu α i v α j + rij ( 8
)
where the wα = O(1) are isolated singular values (describing a rank K matrix), the u α and v α are the dominant eigenvectors of the SVD decomposition and the rij = N (0, σ 2 /L) are i.i.d. terms corresponding to noise, with L = √ NhNv. The {u α } and {v α } are two sets of respectively Nv and Nh-dimensional orthonormal vectors, which means that their components are respectively O(1/ √ Nv) and O(1/ √ Nh), and K ≤ Nv, Nh. We assume Nh < Nv to be the rank of W and wα > 0 and O(1) for all α. Note that in the limit Nv → ∞ and Nh → ∞ with κ def = Nh/Nv fixed and K/L → 0, W W T has a spectrum density ρ(λ) composed of a Marchenko-Pastur bulk of eigenvalues and of set of discrete modes:

ρ(λ) = L 2πσ 2 (λ + -λ)(λ -λ -) κλ 1 1 {λ∈[λ -,λ + ]} + K α=1 δ(λ -w 2 α ),

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

with λ ± def = σ 2 κ 1 4 ± κ -1 4 2 .
The interpretation for the noise term rij is given by the presence of an extensive number of modes at the bottom of the spectrum, along which the variables won't be able to condense but that still contribute to the fluctuations. In the present form our model of RBM is similar to the Hopfield model and recent generalizations [28], the patterns being represented by the SVD modes outside of the bulk. The main difference, in addition to the bipartite structure of the graph, is the non-degeneracy of the singular values wα. The choice made here is to consider K finite, giving Wij = O(1/N ) which means that the thresholds θj (having the meaning of feature detectors) should be O(1) because feature j is detected when an extensive number of spins Si is aligned with Wij.

In addition, this allows us to assume simple distributions for the components of u α and v α (for instance, considering them i.i.d.). Altogether, this defines the statistical ensemble of RBM to which we restrict our analysis of the learning procedure. Another approach would be to consider K = Nh extensive, thereby assuming that all modes can potentially condense even though they are associated to dominated singular values. In that case, the separation between the condensed modes and the rest should be made when order parameters are introduced and the noise would then correspond to uncondensed modes. If the number of condensed modes is assumed to be extensive, then we should instead consider an average over the orthogonal group which would lead to a slightly different mean-field theory [29,30].

Replica symmetric Mean-field equation

Our analysis in the thermodynamic limit follows classical treatments using replicas, like [31,9] for the Hopfield model or [17] for bipartite models. The starting point is to express the average over u, v and rij of the log partition function Z in (2) with the help of the replica trick:

Eu,v,r[log(Z)] = lim p→0 d dp Eu,v,r[Z p ].
î σ 2 2L Ä i,j,a =b s a i s b i σ a j σ b j äó = a =b dQabd Qab 2π × exp î - Lσ 2 2 a =b Ä Qab Qab - Qab Nv i s a i s b i - Qab Nh j σ a j σ b j äó .
The second one is aimed at extracting magnetization's contributions correlated with the modes:

exp Ä L α wαs a α σ a α ä ∝ α dm a α d ma α 2π × exp Ä -L α wα m a α ma α -m a α s a α -ma α σ a α ä , with s a α def = 1 √ L i siu α i and σ a α def = 1 √ L j σ a j v α j , (9)
These variables represent the following quantities: The transformations involve pairs of complex integration variables because of the asymmetry introduced by the two-layers structure in contrast to fully connected models.

m a α ∼ Eu,v,r σ a α ma α ∼ Eu,v,r s a α Qab ∼ Eu,v,r σ a i σ b i Qab ∼ Eu,v,
We obtain the following representation:

Eu,v,r[Z p ] = a,α dm a α d ma α 2π a =b dQabd Qab 2π × exp ¶ -L Ä a,α wαmα mα + σ 2 2 a =b Qab Qab - 1 √ κ A[m, Q] - √ κB[ m, Q] ä© with κ = Nh/Nv and A[m, Q] def = log î S a ∈{-1,1} Eu Ä e √ κσ 2 2 a =b Q ab S a S b +κ 1 4 a,α (wαm a α -ηα)u α S a äó , (10)
B[ m, Q] def = log î S a ∈{-1,1} Ev Ä e √ κσ 2 2 a =b Qab σ a σ b +κ -1 4 a,α (wα ma α -θα)v α σ a äó , (11) (12)
with

θα def = 1 √ L j θjv α j = O(1).
Since {v α } is an incomplete basis we also need to take care of the potential residual

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

transverse parts η ⊥ and θ ⊥ , such that the following decompositions hold:

ηi = η ⊥ i + √ L α ηαu α i , (13)
θj = θ ⊥ j + √ L α θαv α j . ( 14
)
To keep things tractable, both η ⊥ and θ ⊥ will be considered negligible in the sequel. Taking into account these components would lead to the addition of a random field to the effective RS field of the variables and eventually to a richer set of saddle point solutions. Note that the order of magnitude of ηα and θα is at this stage an assumption. If ηi and u α i (or θj and v α j ) were uncorrelated they would scale as 1/ √ L. Moreover, regarding the ensemble average, we will consider ηα and θα fixed in the sequel.

The thermodynamic properties are obtained by first making a saddle point approximation possible by letting first L → ∞ and taking the limit p → 0 afterwards. We restrict here the discussion to RS saddle points [32]. The breakdown of RS can actually be determined by computing the so-called AT line [33] (see Appendix A). At this point we assume a non-broken replica symmetry. The set {Qab, Qab} reduces then to a pair (q, q) of spin glass parameters, i.e. Qab = q and Qab = q for all a = b, while quenched magnetizations on the SVD directions are now represented by {(mα, mα), α = 1, . . . K}.

Taking the limit p → 0 yields the following limit for the free energy:

f [m, m, q, q] = α wαmα mα - σ 2 2 q q + σ 2 2 (q + q) - 1 √ κ Eu,x î log 2 cosh h(x, u) ó - √ κEv,x î log 2 cosh h(x, v) ó . ( 15
)
Assuming a replica-symmetric phase, the saddle-point equations are given by

mα = κ 1 4 Ev,x î v α tanh h(x, v) ó , q = Ev,x î tanh 2 h(x, v) ó ( 16
) mα = κ -1 4 Eu,x î u α tanh h(x, u) ó , q = Eu,x î tanh 2 h(x, u) ó ( 17
)
where

h(x, u) def = κ 1 4 σ √ qx + γ (wγmγ -ηγ)u γ h(x, v) def = κ -1 4 σ √ qx + γ (wγ mγ -θγ)v γ ,
and κ = Nh/Nv, with Eu,x and Ev,x denoting an average over the Gaussian variable x = N (0, 1) and the rescaled components u ∼ √ Nvu α i and v ∼ √ Nhv α j of the SVD modes. We note that the equations are symmetric under the exchange κ → κ -1 , simultaneously with m ↔ m, q ↔ q and η ↔ θ, given that u and v have the same distribution. In addition, for independently distributed u α i and v α j and vanishing fields (η = θ = 0), solutions corresponding to non-degenerate magnetizations have symmetric counterparts: each pair of non-vanishing magnetizations can be negated independently as (mα, mα) → (-mα,mα), generating new solutions. So to one solution presenting n condensed modes, there correspond 2 n distinct solutions.

Phase Diagram

The fixed point equations (16,17) can be solved numerically to tell us how the variables condensate on the SVD modes within each equilibrium state of the distribution and whether a spin-glass or a ferromagnetic phase is present. The important point here is that with K finite and a non-degenerate spectrum the mode with highest singular value dominates the ferromagnetic phase.

In absence of bias (η = θ = 0) and once 1/σ is interpreted as temperature and wα/σ as ferromagnetic couplings, we get a phase diagram similar to that of the Sherrington-Kirkpatrick (SK) model with three distinct phases (see Figure 2)

• a paramagnetic phase (q = q = mα = mα = 0) (P),

• a ferromagnetic phase (q, q, mα, mα = 0) (F),

• a spin glass phase (q, q = 0; mα = mα = 0) (SG).

In general, the lines separating the different phases correspond to second order phase transitions and can be obtained by a stability analysis of the Hessian of the free energy. They are related to unstable modes of the linearized mean-field equations and correspond to an eigenvalue of the Hessian becoming negative.

The (SG-P) line is obtained by looking at the Hessian in the (q, q) sector:

Hqq = m=0 q=0 - 1 2 ï σ 2 σ 4 √ κ √ κσ 4 σ 2 ò
from what results that the spin glass phase develops when σ ≥ 1 2 . This transition line is understood tacking directly into account the spectral properties of the weight matrix. Classically, this is done with the help of the linearized TAP equations and exploiting the Marchenko-Pastur distribution [32]. In our context, the linearized TAP equations read

ï µ ν ò = ï - √ κσ 2 W T W -σ 2 √ κ ò ï µ ν ò
given the variance σ 2 /L of the weights in absence of dominant modes. Then we can show that the paramagnetic phase becomes unstable when the highest eigenvalue of the matrix on the rhs is equal to 1: if λ is a singular value of W , the corresponding eigenvalues Λ ± verify the relation

Λ ± √ κ ± σ 2 √ κΛ ± ± σ 2 = λ 2 .
from which it is clear that the largest eigenvalue Λmax corresponds to the largest singular value λmax. Owing to the Marchenko-Pastur

distribution λmax = σ 2 ( √ κ + 1)(1 + 1/ √ κ) so Λmax verifies Λmax √ κ + σ 2 √ κΛmax + σ 2 = σ 2 ( √ κ + 1) 1 √ κ + 1 .
Λmax = 1 is readily obtained for σ 2 = 1.

2 Note that in [17]

a dependence κ(1 -κ) Ä α(1 -α)in their notation ä is found. This
dependence is hidden in our definition of σ 2 giving L = √ NvN h times the variance of r ij instead of Nv + N h as in their case.

For the (F-SG) frontier we can look at the sector (mα, mα) corresponding to the emergence of a single mode α (written in the spin-glass phase):

Hαα = wα w 2 α Ev,x î (v α ) 2 sech 2 h(x, v)) ó w 2 α Eu,x î (u α ) 2 sech 2 h(x, u) äó wα = mα=0 ï wα w 2 α (1 -q) w 2 α (1 -q) wα ò
From this it is clear that the first mode to become unstable is the mode α with highest singular value wα and this occurs when q and q, solutions of (16,17), verify

(1 -q)(1 -q)w 2 α = 1.
As for the SK model, this line appears to be well below the de Almeida-Thouless (AT) line, which is the line above which the RS solution is stable (see Figure 2, and Appendix A for the computation of the AT line). This means that in principle a replica symmetry breaking treatment would be necessary to properly separate the two phases. However, we will leave aside this point as we are mainly interested in the practical aspects, namely the ability of the RBM to learn arbitrary data, and so we are mostly concerned with the ferromagnetic phase above the AT line.

For the (P-F) line we consider the same sector of the Hessian but now written in the paramagnetic phase, i.e. setting q = 0 in the above equation, and this simply yields the emergence of the single mode α for wα = 1.

Note that all of this is independent on how the statistical average over u and v is performed. Instead, as we shall see later on, the way of averaging influences the nature of the ferromagnetic phase.

Regarding the stability of the RS solution, the computation of the AT line reported in Appendix A is similar to the classical one made for the SK model, though slightly more involved. In fact we were not able to fully characterize, in replica space, all the possible instabilities of the Hessian which would potentially lead to a breakdown of the replica symmetry. At least the one responsible for the ordinary SK model RS breakdown has a counterpart in the bipartite case that gives a necessary condition for the stability of the RS solution:

1 σ 2 > Ex,u Ä sech 4 h(x, u) ä Ex,v Ä sech 4 h(x, v) ä ,
For κ = 1 the terms below the radical become identical and the condition reduces to the one of the SK model, except for the u averages which are not present in the SK model. In Figure 2, is shown the influence on the phase diagram of the value of κ and of the type of average made on u and v.

Nature of the Ferromagnetic phase

Some subtleties arise when considering various ways of averaging over the components of the singular vectors. In [19,20] is emphasized the importance for networks to be able to reproduce compositional states structured by combination of hidden variables. In our representation, we don't have direct access to this property but, in some sense, to the dual one, which is given by states corresponding to combinations of modes. Their presence and their structure are rather sensitive to the way the average over u Fig. 2: Phase diagram in absence of bias and with a finite number of modes, with Gaussian and Laplace distributions for u and v. The dotted line separates the spin glass phase from the ferromagnetic phase under the RS hypothesis. The RS phase is unstable below the AT line. The influence of κ on the AT and SG-F lines is shown. In all cases, the hypothetical SG-F line lies well inside the broken RS phase. Inset: high temperature (σ = 0) stability gap ∆w α corresponding to a fixed point associated to a mode β, expressed as a function of w α and considering various distributions.

and v is performed. In this respect the case in which u α and v α have i.i.d. Gaussian components is very special: all fixed points associated to dominated modes can be shown to be unstable and fixed points associated to combinations of modes are not allowed. To see this, first notice that in such a case the magnetization's part of the saddle point equations (16,17) read

mα = (wα mα -θα)(1 -q) (18) mα = (wαmα -ηα)(1 -q). (19)
Since the role of the bias is mainly to introduce some asymmetry between otherwise degenerated fixed points obtained by sign reversal of at least one pair (mα, mα), let us analyze the situation without fields, i.e. by setting η = θ = 0. We immediately see that as long as the singular values are non degenerate, only one single mode may condense at a time. Indeed if mode α condenses we necessarily have

w 2 α (1 -q)(1 -q) = 1,

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

and this can be verified only by one mode at a time. Looking at the stability of the fixed points, we see that only the fixed point associated to the largest singular value is actually stable (details reported after the introduction of lemma 3.1).

For other distributions like uniform Bernoulli or Laplace, instead, stable fixed points associated to many different single modes or combinations of modes can exist and contribute to the thermodynamics. In order to analyze this question in more general terms we first rewrite the mean-field equations in a convenient way which require some preliminary remarks. We restrict the discussion to i.i.d. variables so that we can consider single variable distributions. Joint distributions will be distinguished from single variable distributions by the use of bold: u = {u α , α = 1, . . . , K}, K being the (finite) number of modes susceptible of condensing.

Given the distribution p and assuming it to be even, we define a related distribution p attached to mode α:

p (u) def = - u -∞ xp(x)dx = ∞ |u| xp(x)dx, (20)
This distribution has some useful properties.

Lemma 3.1. Given that p is centered with unit variance and kurtosis κu, p is a centered probability distribution with variance

∞ -∞ u 2 p (u)du = κu 3 .
Proof. Consider the moments of p . For n odd they vanish while for n even they read:

+∞ -∞ u n p (u)du = 2 ∞ 0 u n p (u)du = 2 ∞ 0 duu n ∞ u xp(x)dx = 2 ∞ 0 xp(x)dx x 0 u n du = 1 n + 1 ∞ -∞ x n+2 p(x)dx,
i.e. the nth even moments of p relate to moments of order n + 2 of p. The lemma then follows from the fact that p has unit variance.

In this respect, the Gaussian averaging is special because we have κu = 3 and p = p. Then the mean-field equations (16,17) corresponding to the magnetizations can be rewritten in a form similar to (18,19) by introducing the variables qα and qα:

mα = (wα mα -θα)(1 -qα), (21) mα
= (wαmα -ηα)(1 -qα), (22)
with

qα = dx e -x 2 /2 √ 2π dvpα(v) tanh 2 Ä κ -1 4 σ √ qx + γ (wγ mγ -θγ)v γ ä , (23) qα
= dx e -x 2 /2 √ 2π dupα(u) tanh 2 Ä κ 1 4 σ √ qx + γ (wγmγ -ηγ)u γ ä , (24)
where pα(u)

def = p (u α ) β =α p(u β ).
This rewriting will prove very useful also in the next section when analyzing the learning dynamics.

Let us now assume, in absence of bias, a non-degenerate fixed point associated to some given mode β with finite (mβ, mβ) and mα = mα = 0, ∀α = β. The fixed point equation imposes the relation

wβ = 1 (1 -qβ)(1 -qβ) def = w(qβ, qβ). (25)
The stability of such a fixed point with respect to any other mode α is related to the positive definiteness of the following block of the Hessian

Hαα =   wα w 2 α Ev,x î (v α ) 2 sech 2 h(x, v) ó w 2 α Eu,x î (u α ) 2 sech 2 h(x, u) ó wα   with, in the present case h(x, u) = κ 1 4 σ √ qx + wβ mβu β and h(x, v) = κ -1 4 σ √ qx + wβ mβv β ,
This reduces to

Hαα = ñ wα w 2 α (1 -q) w 2 α (1 -q) wα ô .
Therefore for the Gaussian averaging case, since qβ = q, qβ = q and given (25), we necessarily have

1 -(1 -q)(1 -q)w 2 α = 1 - w 2 α w 2 β < 0 for wα > wβ,
i.e. the Hessian has negative eigenvalues. This means that if the mode β is dominated by another mode α, the magnetization (mα, mα) will develop until (1-q)(1-q)w 2 α = 1, while mβ will vanish.

For the general case of i.i.d. variables, assuming u α and v α obey the same distribution p, let F and Fα be the cumulative distributions associated respectively to p and pα

F (u) def = u -∞ p(x)dx Fα(u) def = du θ(u -u α )pα(u)dx = - u ∞ du α u α -∞ xp(x)dx.

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

Given the values of (q, q) obtained from the fixed point associated to mode β, we have the following property:

Proposition 3.2. If (i) Fβ(u) < F (u), ∀u ∈ R + then qβ > q and qβ > q, (ii) Fβ(u) > F (u), ∀u ∈ R + then qβ < q and qβ < q,
which in turn implies w(q, q) < wβ (i) and w(q, q) > wβ (ii) with w(q, q)

def = 1 (1 -q)(1 -q) .
Proof. This is obtained by straightforward by parts integration respectively over u and v in equations (16,17), relative to magnetizations.

In other words if Fβ dominates F on R + then there is a positive stability gap defined as

∆wβ def = w(q, q) -wβ (26)
such that there is a non-empty range for higher values of wα ∈ [wβ, w(q, q)[ for which the fixed point associated to mode β corresponds to a local minimum of the free energy. Note that property (i) [resp. (ii)] is analogous (in the sense that it implies it) to pβ having a larger [resp. smaller] variance than p, i.e. κu > 3 [resp. κu < 3]. Therefore distributions p with negative relative kurtosis (κu -3) will tend to favor the presence of metastable states, while the situation will tend to be more complex for probabilities with positive relative kurtosis. Indeed, in the latter case the fixed point associated to the highest mode αmax might not correspond to a stable state if lower modes in the range [w(q, q), wα max [ are present, and fixed points associated to combinations of modes have to be considered. Note that in contrary with the Gaussian case, this can happen because qα is different for each mode and therefore more flexibility is offered by equations (21,22) than from equations (18,19). Let us give some examples. Denote by γu def = κu -3 the relative kurtosis. As already said the Gaussian distribution is a special case with γu = 0. In addition, for instance for p corresponding to Bernoulli, Uniform or Laplace, we have the following properties illustrated in the inset of Figure 2:

• Bernoulli (γu = -2): p(u) = 1 2 δ(u + 1) + δ(u -1) , F (u) = 1 2 θ(u + 1) + θ(u -1) pα(u) = 1 2 θ(1 -u 2 ), Fα(u) = 1 2 θ(1 -u 2 )(u + 1) + θ(u -1)
then Fα(u) > F (u) for u > 0, yielding a positive stability gap.

• Uniform (γu = -6/5):

p(u) = 1 2 √ 3 θ(3 -u 2 ), F (u) = 1 2 √ 3 θ(3 -u 2 )(u + √ 3) + θ(u - √ 3) pα(u) = 1 4 √ 3 θ(3 -u 2 )(3 -u 2 ), Fα(u) = 1 4 √ 3 θ(3 -u 2 )(3u - u 3 3 + 2 √ 3) + θ(u - √ 3).
It can be verified that Fα(u) > F (u) for u > 0, yielding again a positive stability gap.

• Laplace (γu = 3):

p(u) = 1 √ 2 e - √ 2|u| , F (u) = 1 2 + u 2|u| 1 -e - √ 2|u| pα(u) = 1 2 |u| + 1 √ 2 e - √ 2|u| , Fα(u) = F (u) - u 2 √ 2 e - √ 2|u|
.

Here we have Fα(u) < F (u) for u > 0, yielding a negative stability gap. These three examples fall either in condition (i) or (ii), with a stability gap ∆wβ that is either always positive or always negative, independently of wβ. We can also provide examples for which the stability condition may vary with wβ. Consider for instance a sparse Bernoulli distribution, with r ∈ [0, 1] a sparsity parameter:

p(u) = r 2 δ(u + 1 √ r ) + δ(u - 1 √ r ) + (1 -r)δ(u).
The relative kurtosis is in this case

γu(r) = 1 r -3.
Looking at F (u) and Fα(u) it is seen that both conditions (i) and (ii) are not fulfilled, except for r = 1 which corresponds to the plain Bernoulli case. As we see in the inset of Figure 2, for r < 1/3 the stability gap is always negative, meaning that a unimodal ferromagnetic phase is not stable, and it is replaced by a compositional ferromagnetic phase at all temperatures. Instead, for r > 1/3 and at sufficiently high temperature (low wα) the single mode fixed point dominate the ferromagnetic phase.

Laplace distribution: let us look at the properties of the phase diagram in the case of singular vectors' components being Laplace i.i.d., case in which a negative stability gap is expected and it may lead to a compositional phase. For this we need the expression for a sum of Laplace variables to compute the averages involved in (16,17). For this purpose, we define the following distributions:

f (s) = γ du γ λγ 2 e -λγ |u γ | δ(s - γ u γ ), gα(s) = du α λα 4 (λα|u α | + 1)e -λα|u α | γ =α du γ λγ 2 e -λγ |u γ | δ(s - γ u γ ).
Their Laplace transform upon decomposing into partial fractions reads:

f (ω) = γ λ 2 γ λ 2 γ -ω 2 = γ Cγ λ 2 γ λ 2 γ -ω 2 and gα(ω) = λ 2 α λ 2 α -ω 2 γ λ 2 γ λ 2 γ -ω 2 = Cα λ 4 α (λ 2 α -ω 2 ) 2 + γ =α Cγ λ 2 γ λ 2 α λ 2 α -λ 2 γ Ä 1 λ 2 γ -ω 2 - 1 λ 2 α -ω 2 ä .

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

where

Cγ def = δ =γ λ 2 δ λ 2 δ -λ 2 γ
.

From these decompositions we immediately identify

f (s) = 1 2 γ Cγλγe -λγ |s| , gα(s) = λαCα 4 (λα|s| + 1)e -λα|s| + 1 2 γ =α Cγ λγλα λ 2 α -λ 2 γ λαe -λγ |s| -λγe -λα|s| .
This results in the following decomposition of the EA parameters:

q = dxds e - √ 2|s|-x 2 /2 2 √ π γ Cγ[ m] tanh 2 hγ(x, s) (27)
qα = dxds e - √ 2|s|-x 2 /2 2 √ π î 1 √ 2 (|s| + 1 √ 2 )Cα[ m] tanh 2 hα(x, s) (28)
+ γ =α Cγ[ m] (wγ mγ -θγ) 2 tanh 2 hγ(x, s) -(wα mα -θα) 2 tanh 2 hα(x, s) (wγ mγ -θγ) 2 -(wα mα -θα) 2 ó (29) with hγ(x, s) def
= κ -1 4 σ √ qx + (wγ mγ -θγ)s
and

Cγ[ m] def = δ =γ (wγ mγ -θγ) 2 (wγ mγ -θγ) 2 -(wδ mδ -θδ) 2 .
This allows for an efficient resolution of the mean-field equations (16,17,21,22), which let us observe the appearance of a purely compositional phase in the ferromagnetic domain when the modes at the top of the spectrum get close enough. In order to characterize this phase, we consider the stability gap ∆ (n) (wα) for which the range [wa -∆ (n) (wα), wa] lies below the highest mode wa, such that the ferromagnetic states correspond to the condensation of n distinct modes present in this interval, including the highest.

In addition, this will prove useful when analyzing the learning dynamics described in the next section.

4 Learning dynamics of the RBM 4.1 Learning dynamics in the thermodynamic limit A mean field analysis of the learning dynamics has been proposed in [25], in the form of phenomenological equations obtained after averaging over some parameters of the RBM, i.e. by choosing a well defined statistical ensemble of RBMs and using selfaveraging properties in the thermodynamic limit. Here we rederive these equations, we add some details and then explore their properties in the light of the preceding section.

First we project the gradient ascent equations (4-6) onto the bases {uα(t) ∈ R Nv } and {vα(t) ∈ R N h } defined by the SVD of W . Discarding stochastic fluctuations usually inherent to the learning procedure and letting the learning rate γ → 0, the continuous version of (4-6) can be recast as follows:

1 L Ä dW dt ä αβ = sασβ Data -sασβ RBM, (30)
1 √ L Ä dη dt ä α = sα RBM -sα Data, (31)
1 √ L Ä dθ dt ä α = σα RBM -σα Data, (32)
with sα and σα given in (9). We also have

dW dt αβ = δα,β dwα dt + (1 -δα,β) Ä wβ(t)Ω v βα (t) + wα(t)Ω h αβ ä 1 √ L dη dt α = dηα dt - β Ω v αβ ηβ 1 √ L dθ dt α = dθα dt - β Ω h αβ θβ
where

Ω v αβ (t) = -Ω v βα def = du α,T dt u β Ω h αβ (t) = -Ω h βα def = dv α,T dt v β
1 L dwα dt = sασα Data -sασα RBM (33)
dηα dt = sα RBM -sα Data + β Ω v αβ ηβ (34)
dθα dt = σα RBM -σα Data + β Ω h αβ θβ (35)
along with the infinitesimal rotation generators of the left and right singular vectors Ä

Ω v αβ (t) = - 1 wα + wβ dW dt A αβ + 1 wα -wβ dW dt S αβ (36)
Ω h αβ (t) = 1 wα + wβ dW dt A αβ + 1 wα -wβ dW dt S αβ (37)

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

sασβ Data ± sβσα Data ∓ sβσα RBM -sασβ RBM ä .
The dynamics of learning is now expressed in the reference frame defined by the singular vectors of W . The skew-symmetric rotation generators Ω v,h αβ (t) of the basis vectors (induced by the dynamics) tell us how data rotate relatively to this frame. Given the initial conditions, these help us keeping track of the representation of data in this frame. Note that these equations become singular when some degeneracy occurs in W because then the SVD is not uniquely defined. Except from the numerical point of view, where some regularizations might be needed, this does not constitute an issue. In fact only rotations among non-degenerate modes are meaningful, while the rest corresponds to gauge degrees of freedom.

At this point our set of dynamical equations (33)(34)(35)(36)(37) is written in a general form. Our goal is to find the typical trajectory of the RBM within a certain statistical ensemble. For this reason, we make the hypothesis that the learning dynamics is represented by a trajectory in the space {wα(t), ηα(t), θα(t), Ω v,h αβ (t)}, while the specific realization of u α i , v α j and rij in ( 8) can be considered irrelevant and only the way they are distributed is important. We are then allowed to perform an average over u α i , v α j and rij with respect to some simple distributions, as long as this average is correlated with the data. By this we mean that the components sα of any given sample are kept fixed while averaging. In the end, what really matters are the strength and the rotation of the SVD modes, respectively determined by wα(t) and Ω v,h αβ (t). As a simplification and also by lack of understanding of what intrinsically drives their evolution, the distributions of u α i and v α j will be considered stationary in the sequel. Concerning rij, we allow its variance σ 2 /L to vary with time in order to give a minimal description of how the MP bulk evolves during the learning. The detailed dynamics of σ will be derived later in Section 4.3. Using the same notation of Section 3.4 and in particular using the rescaling v ∼ √ Nhv α i , the empirical terms take the form:

σα Data = (sαwα -θα) 1 -qα[s] Data (38)
sασβ Data = sα(sβwβ -θβ) 1 -qβ[s] Data (39)
where

qα[s] def = dx e -x 2 2 √ 2π dvpα(v) tanh 2 Ä κ -1 4 σx + γ (wγsγ -θγ)v γ ä ,
Note that the last equation actually depends on the activation function (hyperbolic tangent in this case), and the term σx corresponds to k rkjsk and is obtained by central limit theorem from the independence of the rkj. qα[s] is the empirical counterpart of the EA parameters q and qα already encountered in Section 3.4, and for simple i.i.d. distributions like Gaussian or Laplace it can be estimated easily. The main point here is that the empirical terms (38,39) define operators whose decomposition over the SVD modes of W functionally depends only on wα, θα and on the projection of the data over the SVD modes of W . These terms are driving the dynamics in a precise way. The adaptation of the RBM to this driving force is given by the . . . RBM terms in (33,34,35), which can be estimated in the thermodynamic limit (see Section 4.3) as a function of wα, θα and ηα alone, by means of the order parameters (mα, mα) given In the inset, the scalar product of the vectors u obtained from the SVD of the data and from the weights w. The us of w are aligned with the SVD of the data at the end of the learning.

in Section 3.2 and once the mean-field equations (16,17) have been solved. Of course, all of this is based on the hypothesis that the RBM stays in the RS domain during learning. Experimental evidence supports this hypothesis (see Section 5).

Linear instabilities

At the beginning of the learning, the elements of the weight matrix W are usually small; therefore, we can analyze the linear behavior of the RBM in order to understand what happens. In particular, we will see that the dynamics of a non-linear RBM at the beginning of the learning can be understood by looking at the stability analysis of the learning process. The purpose of this analysis is to identify which "deformation modes" of the weight matrix are the most unstable, and how they are related to the input data. Additionally, a good feature of the linear case is that no averaging is needed, the dynamics being actually independent on the particular realization of the components u α i and v β j . Also, always relative to the linear case, no distinction has to be made between dominant modes and other modes to be treated as the noise component of equation ( 8), we can simply put all of the modes on the same footing.

Let us analyze the linear regime for an RBM with binary units. The derivation is done by rescaling all the weights and fields by a common "inverse temperature" β and letting this go to zero in equation ( 4). In principle, the stability analysis would lead to assume both the weights and the magnetizations to be small. However, we can assume only the magnetizations to be small and consider a slightly more general case with no approximations. Such a case is analogous to a linear RBM whose magnetizations undergo Gaussian fluctuations, and it is derived by keeping up to quadratic terms of the magnetizations in the mean field free energy:

FMF (µ, ν) 1 2 N i=1 (1 + µi) log(1 + µi) + (1 -µi) log(1 -µi) + 1 2 M j=1 (1 + νj) log(1 + νj) + (1 -νj) log(1 -νj) - i,j Wijµiνj - 1 2 W 2 ij (µ 2 i + ν 2 j ) + N i=1 ηiµi + M j=1 θjνj = 1 2σ 2 v N i=1 µ 2 i + 1 2σ 2 h M j=1 ν 2 i - ij Wijµiνj + N i=1 ηiµi + M j=1 θjνj.
where the variances (σ 2 v , σ 2 h ) of respectively visible and hidden variables read (Nh < Nv):

σ -2 v = 1 + j W 2 ij 1 + α w 2 α (40) σ -2 h = 1 + i W 2 ij = 1 + α w 2 α . (41)
We omitted the quadratic term in Wij coming from the TAP contribution to the free energy, which is optional for our stability analysis. In absence of this term the modes evolve strictly independently, while taking it into account leads to a correction to individual variances which couples the modes. Magnetizations (µ, ν) of visible and hidden variables have now Gaussian fluctuations with covariance matrix

C(µv, µh) def = ñ σ -2 v -W -W T σ -2 h ô -1
We can discard the biases of the data and the related fields (θα, ηα) with a proper centering of the variables, and we consider equation (33) directly involving the covariance matrix of the data expressed in the frame defined by the SVD modes of W sασβ Data = σ 2 h wβ sαsβ Data.

From C(µv, µh) we get the other terms yielding the following equations:

dwα dt = wασ 2 h Ä s 2 α Data - σ 2 v 1 -σ 2 v σ 2 h w 2 α ä Ω v,h αβ = (1 -δαβ)σ 2 h Ä wβ -wα wα + wβ ∓ wβ + wα wα -wβ ä sαsβ Data
Note that these equations are exact for a linear RBM, since they can be derived without any reference to the coordinates of uα and vα over which we average in the non-linear regime. These equations tell us that the learning dynamics drives the rotation of the vectors u α (and v α ) until they are aligned to the principal components of the data, i.e. until sαsβ Data becomes diagonal. Calling ŵ2 α the empirical variance of the data, the system reaches the following equilibrium values:

w 2 α =    ŵ2 α -σ 2 v σ 2 v σ 2 h ŵ2 α if ŵ2 α > σ 2 v , 0 if ŵ2 α ≤ σ 2 v .
assuming (σv, σh) fixed. From this we see that the RBM selects the strongest SVD modes of the data. The linear instabilities correspond to directions along which the variance of the data is above the threshold σ 2 v , and they determine the development of the unstable deformation modes of the weight matrix; during the learning process, these modes will eventually interact following the usual mechanism of non-linear pattern formation encountered for instance in reaction-diffusion processes [34]. Other possible deformations are damped to zero. The linear RBM will therefore learn all the principal components that passed the threshold (up to Nh). Note that this selection mechanism is already known to occur for linear auto-encoders [23] or other similar linear Boltzmann machines [22]. On Fig. 3 we can see the eigenvalues being learned one by one in a linear RBM.

If we take into account the expressions (40,41) for (σv, σh), we see that the system cannot reach a stable solution except for the case in which all the modes are below the threshold at the beginning. Otherwise the modes that are excited first will eventually grow like √ t for a large time, and the excitation threshold will tend to zero for all modes.

In any case, by the definition of a multivariate Gaussian, this simple non-linear analysis describes a unimodal distribution. In order to properly understand the dynamics and the steady-state regime of a non-linear RBM, a well suited mean-field theory is required.

Non-linear regime

In the linear regime, some specific modes are selected and at some point they start to interact in a non-trivial manner. As seen explicitly in (39), the empirical terms in (4-6) involve higher order statistics of the data and then the Gaussian estimation with σ 2 v = σ 2 h = 1 of the RBM response terms sα RBM and sασβ RBM is no longer valid when the interactions kick in. Schematically, the linear regime is valid as long as the RBM is found in the paramagnetic phase. But as soon as one mode passes the linear threshold, the system enters the ferromagnetic phase. Then the proper estimation of the response terms follows from the thermodynamic analysis performed in Section 3, and depends on the assumptions made on the statistical properties of the components of the singular vectors of the weight matrix. In the case of Gaussian i.i.d. components, given the analysis proposed in Section 3.4, we know that the mode with the highest singular value completely dominates the ferromagnetic phase: we expect one single ferromagnetic state characterized by magnetizations aligned to this mode only, while magnetizations correlated to other modes vanish. To be precise, this is the correct picture without fields (η = θ = 0) but we don't expect this picture to drastically change in the case of non-vanishing fields. In fact, solving the mean-field equations in presence of the fields show the appearance of meta-stable states correlated with single dominated modes; however, the free energy difference with respect to the ground state, i.e. the state correlated with the mode with the highest singular value, is of order O L(wαwmax) , which means that the contribution of those meta-stable states become rapidly negligible with large system size.

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

To draw a realistic picture of the learning process we now consider Laplace i.i.d. components for the SVD modes that, as seen in Section 3.4, allow the ferromagnetic phase to be of compositional type. The reason for this is that the Laplace distribution leads to less interference among modes than the Gaussian distribution, so that the modes will weakly interact in the mean-field equations. Solving equations (21,22,27,29) in absence of fields yields the following picture: one fixed point solution will typically have non-vanishing magnetizations {mα, mα} for all α such that wα ∈ [wmax -∆w, wmax], where ∆w is approximately the gap ∆w(q, q) defined in (26). This solution is a degenerate ground state, all other solutions being obtained by independently reversing the signs of the condensed magnetizations (mα, mα). Hence for K condensed modes we get a degeneracy of 2 K . When the fields are included, all these fixed points are displaced in the direction of the fields, and some of them may disappear. In the end we are left with a potentially large amount of nearly degenerate states able to cover the empirical distribution of the data, at least in some simple cases.

Coming back to the learning dynamics the terms corresponding to the response of

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

the RBM in (4,6) are estimated in the thermodynamic limit by means of the previously defined order parameters:

sα RBM = 1 ZTherm ω e -Lf (m ω , mω ,q ω ,q ω ) mω α def = mα Therm, sασβ RBM = 1 ZTherm ω e -Lf (m ω , mω ,q ω ,q ω ) mω α m ω β def = mαmβ Therm.
Here . . . Therm denotes the thermodynamical average and the partition function is expressed, in the thermodynamic limit, as

ZTherm def = ω e -Lf (m ω , mω ,q ω ,q ω )
The index ω runs over all the stable fixed point solutions of (16,17) weighted accordingly to the free energy given by (15). These are the dominant contributions as long as free energy differences are O(1), and the internal fluctuations given by each fixed point are comparatively of order O(1/L). In addition, the dynamics of the bulk can be characterized by empirically defining σ 2 :

σ 2 = 1 L ij r 2 ij ,
whose evolution is:

dσ 2 dt = 1 L ij rij dWij dt , = 1 L ij rij ñ si tanh Ä k rkjsk + κ -1 4 α (wαsα -θα)v α j √ L ä Data -siσj RBM ô
given the independence of ri * (resp. r * j) and u α i (resp. v α i ). Exploiting the self-averaging properties of both the empirical and the response terms with respect to rij, u α i and v α j yields

1 L 2 ij rij siσj Data = σ 2 L 1 -q[s] Data 1 L 2 ij rij siσj RBM = σ 2 L 1 -q Therm , with q[s] def = dx e -x 2 2 √ 2π dvp(v) tanh 2 Ä κ -1 4 σx + γ (wγsγ -θγ)v γ ä .
Summarizing, our equations take the suggestive form

1 L dwα dt = sα(wαsα -θα)(1 -qα[s]) Data -mα(wα mα -θα)(1 -qα) Therm, (42)
dηα dt = mα Therm -sα Data + β Ω v αβ ηβ, (43)
dθα dt = (wα mα -θα)(1 -qα) Therm -(wαsα -θα)(1 -qα[s]) Data + β Ω h αβ θβ, (44)
dσ 2 dt = σ 2 Ä q Therm -q[s] Data ä , (45)
with Ω v,h taking the form of a difference between a data averaging . . . Data and a thermodynamical averaging . . . Therm involving only order parameters. Note here that the wα variables, with respect to the other variables, evolve on a faster time scale. This is our final and main result, which might possibly help improving current learning algorithms of RBMs. From this, it is clear what the learning of an RBM is aimed at: the equations will converge once the dataset is clustered in such a way that each cluster is represented by a solution of the mean-field equations with magnetizations mα and EA parameters qα corresponding respectively to their empirical counterparts sα and qα[s] representing cluster magnetization and variance. In particular, these clusters can somehow be regarded as the attractors in the context of feed-forward networks, defining a partition of the data. This can be seen by starting from random configurations and letting the system evolve using the TAP equations or a MCMC method. At the end the system will end up in one of those clusters (characterized by a fixed point of the mean-field equations). Note that this is the reason why the RBM needs to reach a ferromagnetic phase with many states to be able to match the empirical term in (4) and reach convergence. Additionally, the log likelihood (3) can be estimated in the thermodynamic limit (after normalization by L).

L = ¨√κEx,v î log cosh Ä κ -1 4 σx + α (wαsα -θα)v α äó∂ Data - α ηαsα Data - 1 L log ZTherm ,
As an example, for a multimodal data distribution with a finite number of clusters embedded in a high dimensional configuration space, the SVD modes of W that will develop are the one pointing to the directions of the magnetizations defined by these clusters (which will be almost surely orthogonal, given the high dimensionality of the embedding space). In this simple case the RBM will evolve, as in the linear case, to a state in which the empirical term becomes diagonal, while the singular values will adjust to match the proper magnetization in each fixed point.

We have integrated equations (42,43,44,45,36,37) in simple cases by using the Laplace averaging of the components of the SVD modes and using for the EA parameters the expressions given in (27,29). Basically, the hidden distribution to be B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics modeled is defined by

P (s) = C c=1 pc N i=1 e h c i s i 2 cosh(h c i ) , (46)
i.e. a multimodal distribution composed of C clusters of independent variables, where the magnetization of each variable i in cluster c is given by m c i = tanh(h c i ). Each cluster is weighted by some probability pc. In addition we assume these magnetization vectors m c to be embedded in a low dimensional space of dimension d << N . d defines the rank of W . The initial conditions for W are such that the left singular vectors {uα, α = 1, . . . d} span this low dimensional space. An example of the typical dynamics obtained in the case at hand is shown in Figure 4. In contrast to the linear problem where singular values evolve independently, here we distinctively witness the interaction between singular values: a kind of pressure is exerted by lower modes on higher ones resulting in successive bumps in the dynamics of the top modes. The number of states is roughly multiplied by two each time a mode condenses and get close enough to the top modes. Concerning the dynamics of the fields, we don't really observe convergence towards stable directions. Some (possibly numerical) instability is observed when many modes condense, with both the fields and the number of fixed point solutions becoming very noisy. It is also interesting to see how the magnetizations related to the states are distributed with respect to the dataset. On Figure 5 we see that the fixed points tend (as expected) to settle within dense regions of sample points. However, our coarse description shows some limitations for more complex situations, the number of adjustable parameters being too limited to be able to match arbitrary distributions of clusters. It is then appropriate to think about this behaviour in a mean sense; at least, it is able to reproduce a realistic learning dynamics of the singular values of the weight matrix.

Numerical Experiments

Given the comprehensive theoretical analysis of the RBM model given in the previous sections, we are now able to provide a meaningful description of the learning dynamics for a RBM trained with k-steps contrastive divergence (CDk) [4]. The observations presented in this section will serve as a validation for the theoretical analysis. First, to provide a more direct comparison to section 4.3, we will look at the learning dynamics of an RBM trained on a set of simple synthetic data. Subsequently, we will test the model against real world data by training on the MNIST dataset.

Synthetic dataset

As a simple case, we trained the RBM over the same dataset defined in fig. 4, derived from the simple multimodal distribution in eq. 46 (see Appendix B for details). Thus we set Nv = 1000, Nh = 500 and we trained using 10 4 samples with an effective dimension d = 15 organized in 20 separate clusters. The weights are initialized from a Gaussian distribution with standard deviation σ = 10 -3 , while the hidden bias is initialized to 0 and the visible bias is initialized with the empirical mean of the data

ηi = 1 2 log Å pi 1 -pi ã
where pi is the empirical probability of activation for the ith hidden node. The results of the analysis are shown in fig. 6. We see that the dynamics of the singular values obtained by direct integration of the mean-field equations (Fig. 4) are very well reproduced, the only difference being a slightly higher pressure on the strongest modes. The number of fixed point solutions also seems to follow the same trend but more noise is present, an indication of the fact that the RBM has a tendency to learn spurious fixed points during the training. The learning trajectory on the phase diagram is also of interest; we see that the RBM is initialized in the paramagnetic state as expected and the effect of the learning is to drive the model to the ferromagnetic phase. Once in the ferromagnetic phase, the trajectory slows down and the model is assessed near the critical line between paramagnetic and ferromagnetic states, where the estimate of the weights is most stable (according to [35]). Finally, in Fig. 8 we see how the RBM is able to generate a proper clustering of the data over the spectral modes. In particular, the TAP fixed points of the trained model are well distributed and able to cover the full data distribution, improving over the typical behaviour for Laplace distributed weights that emerged with our theoretical analysis (Fig. 5).

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

MNIST dataset

The MNIST dataset is composed by 70000 handwritten digits (60000 for training, 10000 for testing) of size 28 × 28 pixels. Being highly multimodal, we expect this dataset to push the limits of our spectral analysis. For the training, the initialization of the model is the same one used for the synthetic data, 10000 training samples are used (taken at random from the dataset) and the values of the other hyperparameters . With respect to the linear regime (described in section 4.2) we see in Fig. 9 how the RBM is able to learn the SVD of the dataset quite precisely at the beginning of the training, then the learning dynamics quickly enter the non-linear regime. Even in this highly multimodal scenario, our findings over simple synthetic data seem to be confirmed, as seen in Fig. 6. The high number of modes, however, determines an increase in the magnitude of the singular values of condensed modes and seems to destabilize a bit the learning, making the computation of fixed points less reliable. In fact, as a high number of modes are condensing, the model is not able to get rid of all the spurious fixed points. This problem can be mitigated by using an even smaller learning rate, at the cost of slowing down the training. Probably, using a variable learning rate could be a more practical solution (descreasing the learning rate from time to time to let the model eliminate unneeded fixed points). Concerning the (relative) kurtosis of the mode components distributions, we did not observe a very stable and systematic behavior.

Either we see small fluctuations around zero, either some excursions occur and a finite value in the range [0, 3] is building up either for the u or the v components, coherently to the compositional phase interpretation given previously. The latter is the case for MNIST, as shown in Fig. 7. Additionally the transverse part of the fields, meaning orthogonal to the condensed modes, is usually not completely negligible, in contrary to what we assume in (13,14). This clearly constitutes a limitation of our analysis. These transverse components offer more flexibility for generating and selecting fixed points and interfere in some non-trivial way with the kurtosis property, which possibly explains why we don't get a systematic behavior.

Discussion

Before drawing some perspectives, let us summarize the main outcomes of the present work:

• (i) thermodynamic properties of realistic RBMs: our analysis focused on a non-i.i.d. ensemble of weight matrices, whose derivation has been inspired by empirical observations obtained by training RBMs on real data.

• (ii) RS equations and compositional phase: we found a way of writing the RS equations for the RBM (in particular with equations (21,22,23,24)) which leads to a simple characterization of the ferromagnetic phase where the RBM is assumed to operate. Schematically, a negative relative kurtosis for the distribution of the singular vectors' components favors the proliferation of metastable states, while a positive one tends to favor a compositional phase. In particular, we were able to precisely address a concrete case presenting the compositional phase by considering a Laplace distribution for the singular vectors' components.

• (iii) a set of equations representing a typical learning dynamics that defines a trajectory in {wα(t), ηα(t), θα(t), Ω v,h αβ (t), σ 2 (t)}. The spectrum of the dominant singular values, represented by {wα(t)} and expressing the information content of the RBM, is playing the main role. The bulk of dominated modes corresponding to noise sees its dynamics summarized by the evolution of σ 2 (t). Rotations of dominant singular vectors during the learning process are given by Ω v,h while the projections of the biases along the main modes are given by η and θ. These equations have been obtained by averaging over the components of left and right SVD vectors of the weight matrix, while keeping fixed the quantities considered to be relevant. This averaging actually corresponds to a standard self-averaging assumption in a RS phase.

• (iv) a clustering interpretation of the training process is obtained through equations (42,43,44,45) where it is explicitly shown the kind of matching that the RBM is trying to perform between the order parameters obtained from the fixed point solutions and their empirical counterparts in the non-linear regime.

A natural clustering of the data can actually be defined by assigning to each sample the fixed point obtained after initializing the fixed point equations with a visible configuration corresponding to that same sample.

The main picture emerging from the present analysis is that of a set of clusters corresponding to the fixed points of the RBM, which try to uniformly cover the support of the dataset. A full understanding of the mechanism by which the RBM manages to properly cover the dataset is still lacking, even though the case of Laplace distributed singular vectors' components gives some insights. By comparison, real RBMs have more flexibility than the simple "mean Laplace RBM" considered in Section 3.4 and they can produce a good covering of the data manifold. We were not yet able to precisely pinpoint the main ingredients for that mechanism, even though we suspect the transverse biases (orthogonal to the modes) of the hidden units to be the missing ingredient in our analysis.

From the theoretical point of view we would like to see how these results can be adapted to more complex models like DBM or generative models based on convolutional networks. In particular we would like to understand whether adding more layers can facilitate the covering of the dataset by fixed points. From the practical point of view these results might help to orientate the choice of the hyper-parameters used for training an RBM and to refine the criteria for assessing the quality of a learned RBM. For instance, the choice of the number of hidden variables is dictated by two considerations: the effective rank of W , i.e. the number of relevant modes to be considered, and the level of interaction between these modes. Using less hidden variables gives more compact RBMs and reduces the rank of W to its needed value, but it also leads to modes with stronger interactions, which means less flexibility for generating a good covering of fixed points.

A AT line

The stability of the RS solution to the mean-field equations is studied along the lines of [33] by looking at the Hessian of the replicated version of the free energy and identifying eigenmodes from symmetry arguments. Before taking the limit p → 0 the free energy reads

f [m, m, Q, Q] = a,α wαm a α ma α + σ 2 2 a =b Qab Qab - 1 √ κ Ap[m, Q] - √ κBp[ m, Q],
with Ap and Bp given in (10,11). Assuming the small perturbations

m a α = mα + a α ma α = mα + ¯ a α Qab = q + ηab Qab = q + ηab,
around the saddle point (mα, mα, q, q), the perturbed free energy reads

Aαβ def = (δαβ -mαmβ)wαwβ Bαβ def = Ä Ex,v v α v β tanh 2 ( h(x, v)) -mαmβ ä wαwβ Cα def = κ 1/4 σ 2 2 mα(1 -q)wα Dα def = κ 1/4 σ 2 2 Ä Ex,v v α tanh 3 ( h(x, v)) -mαq ä wα E0 def = √ κσ 4 4 (1 -q 2 ) E1 def = √ κσ 4 4 q(1 -q) E2 def = √ κσ 4 4 Ä Ex,v tanh 4 ( h(x, v)) -q 2 ä with h(x, u) def = κ 1/4 √ qσx + α (mαwα -ηα)u α ,
Conjugate quantities are obtained by replacing mα by mα, q by q, u α by v α , ηα by θα and κ by 1/κ. As for the SK model, the 2Kp × 2Kp Hessian thereby defined can be diagonalized with the help of three similar sets of eigenmodes corresponding to different permutation symmetries in replica space. The first set corresponds to 2K + 2 replica symmetric modes defined by η a α = ηα

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

and ηab = η solving the linear system

wα 2 -λ ¯ α - 1 2 Āαα α + β Āαβ + (p -1) Bαβ β + (p -1) Cα + (p -1)(p -2) 2 Dα η = 0 wα 2 -λ α - 1 2 Aαα¯ α + β Aαβ + (p -1)Bαβ ¯ β + (p -1)Cα + (p -1)(p -2) 2 Dα η = 0 σ 2 2 -λ η + α Cα + p -2 2 Dα α + 2 Ē0 + 2(p -2) Ē1 + (p -2)(p -3) 2 Ē2 η = 0 σ 2 2 -λ η + α Cα + p -2 2 Dα ¯ α + 2 E0 + 2(p -2)E1 + (p -2)(p -3) 2 E2 η = 0
with eigenvalue λ solving a polynomial equation of degree 2K + 2 corresponding to a vanishing determinant in the above system. The second set corresponds to a broken replica symmetry where one replica a0 is different from the others

( a α , ¯ a α ) = ( α, ¯ α) for a = a0 (1 -p)( α, ¯ α) for a = a0 (ηab, ηab) = (η, η) for a, b = a0 (1 -p 2 )(η, η) for a = a0 or b = a0
This set has dimension (2K + 2)(p -1). Its parameterization is obtained by imposing orthogonality with the previous one. The corresponding system reads

wα 2 -λ ¯ α - 1 2 Āαα α + β ( Āαβ -Bαβ) β + p -2 2 Cα -Dα η = 0 wα 2 -λ α - 1 2 Aαα¯ α + β (Aαβ -Bαβ)¯ β + p -2 2 Cα -Dα η = 0 σ 2 2 -λ η + α ( Cα -Dα) α + 2 Ē0 + (p -4) Ē1 -(p -3) Ē2 η = 0 σ 2 2 -λ η + α (Cα -Dα)¯ α + 2 E0 + (p -4)E1 -(p -3)E2 η = 0
Finally the eigenmodes of the Hessian are made complete by considering a broken symmetry where two replicas a0 and a1 are different from the others, with the following parameterization dictated again by orthogonality constraints with the previous sets:

( a α , ¯ a α ) = 0, (ηab, ηab) =        (η, η) for a, b = a0 3-p 2 (η, η) for a ∈ a0, a1 or b ∈ a0, a1 (p-2)(p-3) 2 (η, η) for (a, b) = (a0, a1).
The dimension of this set is now p(p-3), and it represents eigenvectors iff the following system of equations is satisfied

σ 2 2 -λ η + 2( Ē0 -2 Ē1 + Ē2)η = 0 σ 2 2 -λ η + 2(E0 -2E1 + E2)η = 0
The corresponding eigenvalues read

λ = σ 2 2 ± 2 ( Ē0 -2 Ē1 + Ē2)(E0 -2E1 + E2),
with degeneracy p(p -3)/2. Finally the RS stability condition reads

1 σ 2 > Ex,u Ä sech 4 h(x, u) ä Ex,v Ä sech 4 h(x, v)
ä , which reduces to the same form of the AT line for the SK model when κ = 1, except for the u and v averages that are specific to our model. As seen in Figure 2 the influence of κ is very limited.

B Synthetic dataset

The multimodal distribution modeling the N-dimensional synthetic data is

P (s) = C c=1 pc N i=1 e h c i s i 2 cosh(h c i ) , ( 47
)
where C is the number of clusters, pc is a weight and h c is a hidden field for cluster c.

The values for pc are taken at random and normalized, while to compute h c i we take into account the magnetizations m c i = tanh(h c i ). Expanding over the spectral modes, we can set an effective dimension d by constraining the sum to the range α = 1, . . . , d

m c i = d α=1 m c α u α i (48)
Clusters' magnetizations m c α are drawn at random between [-1, 1] and normalized with the factor

Z = … α m 2 α d • r , r = tanh(η) ( 49
)
where r is introduced to decrease the clusters' polarizations (in our simulations, we used η = 0.3). The spectral basis u α i is obtained by drawing at random d N-dimensional vectors and applying the Gram-Schmidt process (which can be safely employed as N is supposedly big and thus the initial vectors are nearly orthogonal). The hidden fields are then obtained from the magnetizations

h c i = tanh -1 (m c i ) (50)
and the samples are generated by choosing a cluster according to pc and setting the visible variables to ±1 according to

p(si = 1) = 1 1 + e -2h c i (51)

B. Thermodynamics of Restricted Boltzmann Machines and Related Learning Dynamics

four methods above can handle all the challenges co-currently raised in our work. First, all of these assume the testing instances to have fully observed feature profiles. They don't consider coping with incomplete testing instances by design. Second, all of them are designed specifically for multi-label learning and adapting them to multi-class classification is not straightforward.

More recently, methods based on Deep Latent Variable Models (DLVM) have been proposed to deal with missing data. In [START_REF] Mattei | MIWAE: Deep generative modelling and imputation of incomplete data sets[END_REF], the Variational Autoencoder [START_REF] Diederik | Auto-encoding variational bayes[END_REF] has been adapted to be trained with missing data and a sampling algorithm for data imputation is proposed. Other approaches based on Generative Adversarial Networks (GAN) by [GPAM + 14] are proposed in [START_REF] Yoon | GAIN: Missing data imputation using generative adversarial nets[END_REF] and [START_REF] Cheng-Xian | Learning from incomplete data with generative adversarial networks[END_REF]. Impressive results on image datasets are displayed for these models, at the price of a rather high model complexity and the need for a large training set. In addition these works are focused on features reconstruction, and additional specifications and fine-tuning are required to be able to take partially observed labels into account. The models specifications are quite involved and any new specificity of the dataset may increase both the cost and the difficulty in training (especially for the approaches based on GANs).

In this paper we choose to address this problem in a more economical and robust manner. We consider the old and simple architecture of the Restricted Boltzmann Machine and adapt it to the multi-output learning context (RBM-MO) with missing data. The RBM-MO method serves as a generative model which collaboratively learns the marginal distribution of features and label assignments of input data instances, despite the incomplete observations. Building on the ideas expressed in [START_REF] Nijman | Using boltzmann machines to fill in missing values[END_REF][START_REF] Ghahramani | Supervised learning from incomplete data via an em approach[END_REF] we adapt the approach to the more effective contrastive divergence training procedure [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF] and provide results on various real-world datasets. The advantage of the RBM-MO model is that of providing a robust and flexible method to deal with missing data, with little additional complexity with respect to the classic RBM. Indeed, the trained model can be naturally applied to both transductive and inductive scenarios, achieving superior multi-output classification performance then state-of-the-art baselines. Moreover, it works seamlessly with multi-class and multi-label tasks, providing a unified framework for multi-output learning.

Overview of Restricted Boltzmann Machines

An RBM is a Markov random field with pairwise interactions defined on a bipartite graph formed by two layers of non-interacting variables: the visible nodes represent instances of the input data while the hidden nodes provide a latent representation of the data instances. V and H will denote respectively the sets of visible and hidden variables. In our setting, the visible variables will further split into two subsets V f and V corresponding respectively to features and labels, such that V = V f +V . The visible variables form an explicit representation of the data and are noted v = {v i , i ∈ V}. The hidden nodes h = {h j , j ∈ H} serve to approximate the underlying dependencies among the visible units.

In this paper, we will work with binary hidden nodes h j ∈ {0, 1}. The variables corresponding to the visible features will be either real with a Gaussian prior or binary, depending on the data to model, and labels variables will always be binary (v i ∈ {0, 1}). The joint probability distribution over the nodes is defined through an energy function

P (v, h) = e -E(v,h) Z p prior (v), E(v, h) = - i∈V,j∈H v i w ij h j - i∈V a i v i - j∈H b j h j (1)
where a i and b j are biases acting respectively on the visible and hidden units and w ij is the weight matrix that couples visible and hidden nodes. p prior is in product form and encodes the nature of each visible variable, either with a Gaussian prior p prior = N (0, σ 2 v ) or a binary prior p prior (v) = δ(s 2s). Z = v,h p prior (v)e -E(v,h) is the partition function. The classical training method consists in maximizing the marginal likelihood over the visible nodes P (v) = h P (v, h) by tuning the RBM parameters θ = {w ij , a i , b j } via gradient ascent of the log likelihood L(v; θ).

The tractability of the method relies heavily on the fact that the conditional probabilities P (v|h) and P (h|v) are given in closed forms. In our case these read:

P (v|h) = i∈Vf e j∈H viwij hj +aivi p prior (v i ) vi e j∈H viwij hj +aivi p prior (v i ) i∈V σ j∈H v i w ij h j + a i v i , (2)
P (h|v) = j∈H σ i∈V v i w ij + b j , (3)
where σ(x) = 1/(1 + e -x ) is the logistic function. The gradient of the likelihood w.r.t. the weights (and similarly w.r.t. the fields a i and b j ) is given by

∂L(v; θ) ∂w ij = v i h j p(h j |v) data -v i h j RBM (4)
where the brackets data and RBM respectively indicate the average over the data and over the distribution (1). The positive term is directly linked to the data and can be estimated exactly with (3), while the negative term is intractable. Many strategies are used to compute this last term: the contrastive divergence (CD) approach [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF] consists in estimating the term over a finite number of Gibbs sampling steps, starting from a data point and making alternate use of ( 2) and ( 3); in its persistent version (PCD) the chain is maintained over subsequent mini-batches; using mean-field approximation [START_REF] Marylou | Training restricted Boltzmann machines via the Thouless-Anderson-Palmer free energy[END_REF] the term is computed by means of a low-couplings expansion.

Learning RBM with incomplete data

The RBM is a generative model able to learn the joint distribution of some empirical data given as input. As such, it is intrinsically able to encode the relevant statistical properties found in the training data instances that relate features and labels, and this makes the RBM particularly suitable to be used in the multi-output setting in the presence of incomplete observations. In this sense, the most natural way to deal with incomplete observations is to marginalize over the missing variables; in this section we show how the contrastive divergence algorithm can be adapted to compute such marginals. Given a partially-observed instance v, we have a new partition of the visible space V = O + M, where O is a subset of observed values of v that can correspond both to features and labels. v o = {v i , i ∈ O} and v m = {v i , i ∈ M} denote respectively the observed and missing values of v. The probability over the observed variables v o is given by (θ representing the parameters of the model)

P (v o ) = Z O [θ] Z ∅ [θ] , Z O [θ] = i∈M p prior (v i )dv i × e k∈V akvk j∈H 1 + exp k∈V w kj v k + b j
Taking the log-likelihood and then computing the gradient with respect to the weight matrix element w ij (also similarly for the fields a i and b j ), we obtain two different expressions for i ∈ O and i ∈ M.

∂ log Z O [θ] ∂w ij = v i hj h j p(h j |v o ) i ∈ O, ∂ log Z O [θ] ∂w ij = hj dv i v i h j p(v i , h j |v o ) i ∈ M (5)
The gradient of the LL over the weights (4) now reads

∂L(v; θ) ∂w ij = I o (i)v i hj h j p(h j |v o ) data + 1 -I o (i) hj dv i v i h j p(h j |v o ) data -v i h j RBM (6)
where I o is the indicator function of the samples dependent set O. The observed variables v i , i ∈ O are pinned to the values given by the training samples. In terms of our model, the pinned variables play the role of an additional bias over the hidden variables of a RBM where the ensemble of visible variables is reduced to the missing ones.

With respect to the non-lossy case where p(h j |v) is given in closed form, here we need to sum over the missing variables in order to estimate p(h j |v o ). This means that also the positive term of the gradient ( 6) is now intractable and we need to approximate it. For CD training, we can simply perform Gibbs sampling over the missing variables (keeping fixed the observed variables). Details are reported in Alg. 1. We note that the extra computational burden of Lossy-CD with respect to standard CD is due only to the extra Gibbs sampling steps in the positive term. Given that the observed variables strongly bias the sampling procedure speeding up convergence, only few sampling steps are needed to compute this term. Indeed, in our experiments we observed that a single sampling step (Lossy-CD1) is enough, making the additional complexity minimal. Finally, we note that the same method can be applied to PCD and mean-field training procedures. In the first case, it is sufficient to keep track of an additional persistent chain, which requires little extra memory

C. Robust Multi-Output Learning with Highly Incomplete Data via Restricted Boltzmann Machines

and no extra computational complexity. In the second case, we only need to substitute Gibbs sampling with iterative mean-field equations. iterate eq. ( 2), (3) (k steps) to compute v i h j model Full update: end for 14: end for 4 Mean-field based imputation with RBM As a generative model, the trained RBM can be used to sample new data. For imputation of missing features and labels we just need to use the observed portions of our data to bias the sampling procedure in the same way as for the computation of the positive term in Alg. 1. Namely, we estimate p(v m |v o ) by pinning the observed variables and iterating CD/PCD or mean-field to approximate the equilibrium values of the missing variables. In case of a high percentage of missing observations, however, we might expect the observed variables to be correlated to many different equilibrium configurations, such that the sampling could be biased towards the wrong sample. To overcome this problem, we simply average over multiple mean-field imputations for each incomplete data instance.

More in details, let {p i , i ∈ V } and {q j , j ∈ H} be the marginal probabilities respectively of visible labels and hidden variables to be activated and {m i , i ∈ V f } the marginal expectation of the visible features variables. Mean-field equations at lowest order (O(1/N ), N being the size of the system) express self-consistent relations among these quantities

m i = j∈H w ij q j + a i σ 2 v ∀i ∈ V f \O p i = σ j∈H w ij q j + a i ∀i ∈ V \O (7)
q j = σ i∈Vf w ij m i + i∈V w ij p i + b j (8)
Higher order terms corresponding to TAP equations are discarded [M 17]. These equations can be efficiently solved by iteration starting from random configurations until a fixed point is reached. Observed variables are simply introduced by pinning their corresponding probabilities (0 or 1 for label variables) or their marginal expectation (for feature variables) to the observed values. In practice we run these fixed-point equations N f ∼ 10 times and the imputations are obtained by simple average mi

= 1 N f Nf n=1 m (n) i p i = 1 N f Nf n=1 p (n) i .
In the multi-label setting, the predictor is the indicator function pi = (p i > t) (t is learned, it is chosen to maximize the accuracy for known labels), while for class labels we have pi = 5 Experimental Study

Experimental configuration

To evaluate the efficiency of RBM-MO we compare its performance against CLE, NoisyIMC and MC-1, which provide state-of-the-art baselines.

For the transductive experiments we randomly hide features and labels of the whole dataset to generate incomplete data for training, and we compute appropriate scores for the reconstruction of missing features and labels. In the inductive test, instead, we split the whole dataset into non-overlapping training and testing sets. Concerning the training set the same protocol is used as in the transductive test. For the test set the difference is that now all labels are hidden. Once the classifier is trained, it is applied on the test set to predict the labels. We still randomly hide the entries of test features vectors, so as to form an incomplete testing set. Finally, in the splitting we use 70% of the data instances for training and the remaining 30% for testing.

We denote by q f ea , q ml and q mc the percentage of masked features, labels and classes labels respectively. Note that a masked class label means that all binary variables attached to the classes of a given label are masked together. These rates of masking are kept identical in the learning and test sets.

In the transductive test, we compute the Root Mean Squared Error (RMSE) to measure the reconstruction accuracy with respect to the missing feature values. Furthermore, for the reconstructed labels we calculate Micro-AUC scores and Hamming-accuracy [START_REF] Madjarov | An extensive experimental comparison of methods for multi-label learning[END_REF] in the multi-label scenario, and Averaged AUC plus Accuracy [START_REF] Yue Wu | Experimental comparisons of multi-class classifiers[END_REF] in the multi-class case. In the tables, we define Hamming-accuracy as 1-Hamming loss to keep a consistent variation tendency with the AUC scores. In the inductive test we only compute the scores on the reconstructed labels, since reconstructing missing features is not the goal of inductive classification.

We run the test as described 10 times with different realizations of the missing features and labels. Average and standard deviation of the computed scores are recorded to compare the overall performances. In the tables, we use red fonts to denote the best reconstruction and classification performances among all the algorithms involved in the empirical study. The bold black font is used to highlight the performance of the proposed RBM-MO method.

For the baselines, we used grid search to choose the optimal parameter combination following the suggested ranges of parameters as in [START_REF] Han | Multi-label learning with highly incomplete data via collaborative embedding[END_REF].

The RBM-MO is trained following the guidelines in [START_REF] Hinton | A practical guide to training restricted Boltzmann machines[END_REF]. We always use binary variables for the hidden layer, while in the visible layer we use binary variables for MNIST and Gaussian variables for the other datasets. In all the simulations, we fix the number of hidden nodes to 100. The learning rate η is fixed to 0.001 and the size of the mini-batches to 10. During training the number of Gibbs steps is set to k = 1 while for imputation we iterate the mean-field equations 10 times. As a stopping condition, we considered the degradation of the transductive AUC scores with a look-ahead of 500 epochs

Summary of datasets

We consider 3 publicly available datasets related to image processing. These datasets cover both multi-label and multi-class learning tasks, and they are popularly used as benchmark datasets in multi-output learning research.

In addition, we consider the challenging scenario of abnormality detection on IoT devices. The relevant dataset, that we call EventCat, consists in security telemetry data collected from various network appliances (e.g. smart watches, smartphones, driving assistance systems...), each reporting a features vector whose entries indicate the occurring frequency of a specific type of alert (e.g. downloading suspicious files, login failures, unfixed vulnerabilities...). Multiple labels are assigned to each device in the collected dataset, corresponding to a variety of categories of security threats. Some details about the datasets are reported in Table .5.

Qualitative results on MNIST

A qualitative evaluation of the performance of the RBM-MO model is given by looking at features reconstruction for the MNIST dataset, as reported in Fig. 1. The model at hand has been trained over a dataset in which 50% of the features were missing. To assess the robustness of the method, we computed the reconstructions in the highly challenging case in which 80% of the features were missing. Apart from some smoothing due to the employment of mean-field imputations, the reconstructed samples look reasonably realistic. In general, from the qualitative point of you the results are comparable to those obtained with more complex and expensive DLVMs like MIWAE and MisGAN [MF19, LJM19].

Empirical results

The transductive results for MNIST (multi-class) and Scene (multi-label) datasets are reported in tables 1 and 2. Going into the details, we first observe that RBM-MO is by a large margin more efficient than all of the Figure 1: Features reconstruction by RBM-MO trained over an incomplete dataset with 50% missing-at-random features, whose classification accuracy has been measured to be around 91%. The first block shows some complete testing instances. The second and third block show the same testing instances after hiding respectively 50% and 80% of the pixels. The last two columns show the results of the mean-field imputations over the incomplete testing instances.

baselines for the inference of class labels (table 1), probably because it is able to encode more complex statistical properties.

On the multi-label problems, the situation is still in favour of RBM-MO but with less margin (table 2), in particular at a larger percentage of missing features. Now if we look at the reconstruction error on these datasets we observe that RBM-MO generally achieves a higher reconstruction accuracy than the other opponents, especially on the MNIST dataset. The results verify empirically the basic motivation of using a generative model such as the RBM: incomplete features and labels can provide complementary information to each other, so as to better recover the missing elements. The variance of the results is omitted in the tables by lack of space. For RBM-MO the standard deviation of the derived RMSE, AUC and accuracy scores is not larger than 0.01 over the different datasets. Although the RMSE scores reported by the baseline methods look comparable to the RBM-MO ones, and in certain cases they are better, they also come with a slightly higher variance, such that the RBM-MO seems to be more efficient and robust for features reconstruction.

Except MC-1, all the baseline methods are used for inductive learning. As in the transductive test, we show only the mean of the derived metrics in the tables. Nevertheless, we have similar variance ranges for the computed scores as reported in the transductive test. Clearly RBM-MO is much better adapted to this setting than the baseline methods both for multi-class (table 3) and multi-label learning. The baseline inductive methods CLE and NoisyIMC are specifically designed for multi-label learning and their performance deteriorates significatively in the multi-class scenario. By comparison, RBM-MO can be adapted seamlessly to multi-class and multi-label learning, producing consistently good performances.

For the EventCat dataset, inductive results are reported in table 4. Even with highly incomplete training data, RBM-MO produces the best predictions over partially observed testing data instances.

Conclusion

Machine learning is witnessing a race to high complexity models eager for large data and computational power. In the context of multi-output classification in a challenging scenario -(i) learning with highly incomplete features and partially observed labels; ii) applying the learnt classifier with incomplete testing instances)we advocate instead for simple probabilistic and interpretable models. After refining the learning of the RBM model, we give empirical evidences that it can be efficiently adapted to this context on a great variety of datasets. Experiments are conducted on both public databases and a real-world IoT security dataset, showing various sizes of training sets as well as features and labels vectors. Our approach consistently outperforms the state-of-the-art robust multi-class and multi-label learning approaches with imperfect training data, indicating good usability for practical applications.

C. Robust Multi-Output Learning with Highly Incomplete Data via Restricted Boltzmann Machines

Unfortunately, however, typical strategies employed in neural networks training do not scale well for objective functions like the aforementioned ones; in fact, through the change of variable formula, the logarithm of the absolute value of the determinant of the Jacobian appears in the objective. Its exact computation, let alone its optimization, quickly gets prohibitively computationally demanding as the data dimensionality grows.

A large part of the research on deep density estimation, generally referred to under the term autoregressive normalizing flows, has therefore been dedicated to considering a restricted class of transformations such that the computation of the Jacobian term is trivial [14,44,15,34,25,12], thus imposing a tradeoff between computation and expressive power. While such models can approximate arbitrary probability distributions, the extracted features are strongly restricted based on the imposed triangular structure, which prevents the system from learning a properly disentangled representation. Other strategies involve the optimization of an approximation of the exact objective [5], and continuous-time analogs of normalizing flows for which the likelihood (or some approximation thereof) can be computed using relatively cheap operations [13,19].

In this work, we provide an efficient way to optimize the exact maximum likelihood objective for deep density estimation as well as for learning disentangled representations by latent variable models. We consider a nonlinear, invertible transformation from the observed to the latent space which is parameterized through fully connected neural networks. The weight matrices are merely constrained to be invertible. The starting point is that the parameters of the linear transformations are matrices; this allows us to exploit properties of the Riemannian geometry of matrix spaces to derive parameter updates in terms of the relative gradient, which was originally introduced as the natural gradient in the context of linear ICA [11,2], and which can be feasibly computed. We show how this can be integrated with the usual backpropagation employed to compute gradients in neural network training, yielding an overall efficient way to optimize the Jacobian term in neural networks. This is a general optimization approach which is potentially useful for any objective involving such a Jacobian term, and is likely to find many applications in diverse areas of probabilistic modelling, for example in the context of Bayesian active learning for the computation of the information gain score [48], or for fitting the reverse Kullback-Leibler divergence in variational inference [54,7]. The computational cost of our proposed optimization procedure is quadratic in the input sizeessentially the same as ordinary backpropagation-which is in stark contrast with the cubic scaling of the naive way of optimizing via automatic differentiation. The joint asymptotic scaling of forward and backward pass as a function of the input size is therefore the same that aforementioned alternative methods achieve by imposing strong restrictions on the neural network structure [44] and thus on the class of functions they can represent. In contrast, our approach allows to efficiently optimize the exact objective for neural networks with arbitrary Jacobians.

In sections 2 and 3 we review maximum likelihood estimation for latent variable models, backpropagation and the Jacobian term for neural networks, and discuss the complexity of the naive approaches for optimizing the Jacobian term. Then in section 4 we discuss the relative gradient, and show how it can be integrated with backpropagation resulting in an efficient procedure. We verify empirically the computational speedup our method provides in section 5.

Background 2.1 Maximum likelihood for latent variable models

Consider a generative model of the form

x = f (s) (1)
where s ∈ R D is the latent variable, x ∈ R D represents the observed variable and f : R D → R D is a deterministic and invertible function, which we refer to as forward transformation. Under the model specified above, the log-likelihood of a single datapoint x can be written as

log p θ (x) = log p s (g θ (x)) + log | det Jg θ (x)| , (2)
where g θ is some representation with parameters θ of the inverse transformation 2 of f ; Jg θ (x) ∈ R D×D its Jacobian computed at the point x, whose elements are the partial derivatives [Jg θ (x)] ij = ∂g i θ (x)/∂x j ; and p θ and p s denote, respectively, the probability density functions of x and of the latent variable s under the specified model. In many cases, it is additionally assumed that the distribution of the latent variable is sufficiently simple; for example, that it factorizes in its components, log p θ (x) =

In this case, the problem can be interpreted as nonlinear independent component analysis (nonlinear ICA), and the components of g θ (x) are estimates of the original sources s.

Another variant of this framework can be developed to solve the problem that nonlinear ICA is, in general, not identifiable without additional assumptions [29]; that means, even if the data is generated according to the assumed model, there is no guarantee that the recovered sources bear any simple relationship to the true ones. In order to obtain identifiability, it is possible to consider models [27,28,30,20] in which the latent variables are not unconditionally independent, but rather conditionally independent given an additional, observed variable u ∈ R d ,

log p θ (x|u) = i log p i (g i θ (x)|u) + log | det Jg θ (x)| , (4)
where d can be equal to or different from D depending on the model.

Maximum likelihood estimation for the model parameters amounts to finding, through optimization, the parameters θ * such that the expectation of the likelihood given by the expression in equation ( 3) is maximized. For all practical purposes, the expectation will be substituted with the sample average. Specifically, for optimization purposes, we will be interested in the computation of a gradient of such term on mini-batches of one or few datapoints, such that stochastic gradient descent can be employed.

Neural networks and backpropagation

Neural networks provide a flexible parametric function class for representing g θ through a sequential composition of transformations, g θ = g L • . . . • g 2 • g 1 , where L defines the number of layers of the network. When an input pattern x is presented to the network, it produces a final output z L and a series of intermediate outputs. By defining z 0 = x and z L = g θ (x), we can write the forward evaluation as z k = g k (z k-1 ) for k = 1, . . . , L .

Each module g k of the network involves two transformations, (a) a coupling layer C Wk , that couples the inputs to the layer with the parameters W k to optimize; (b) other arbitrary manipulations σ of inputs/outputs. Typically, these are element-wise nonlinear activation functions with fixed parameters; we can for simplicity think of them as operations of the form σ(x) = (σ (x 1 ) , . . . , σ (x n )) applied to vector variables.

The resulting transformation can thus be written as g k (z k-1 ) = σ(C Wk (z k-1 )).

We will focus on fully connected modules, where the coupling C W is simply a matrix-vector multiplication between the weights W k and the input to the k-th layer; overall, the transformation operated by such a module can be expressed as σ(W k z k-1 ). Another kind of coupling layer is given by convolutional layers, typically used in convolutional neural networks [36].

The parameters of the network are randomly initialized and then learned by gradient based optimization with an objective function L, which is a scalar function of the final output of the network. At each learning step, updates for the weights are proportional to the partial derivative of the loss with respect to each weight.

The computation of these derivatives is typically performed by backpropagation [47], a specialized instance of automatic differentiation. Backpropagation involves a two-phase process. Firstly, during a forward pass, the intermediate and final outputs of the network z 1 , . . . , z L are evaluated and a value for the loss is returned. Then, in a second phase termed backward pass, derivatives of the loss with respect to each individual parameter of the network are computed by application of the chain rule.

The gradients are computed one layer at a time, from the last layer to the first one; in the process, 3 D. Relative gradient optimization of the Jacobian term in unsupervised deep learning the intermediate outputs of the forward pass are reused, employing dynamic programming to avoid redundant calculations of intermediate, repeated terms. 3In matrix notation, the updates for the weights of the k-th fully connected layer W k can then be written as

∆W k ∝ z k-1 δ k , (6)
where δ k is the cumulative result of the backward computation in the backpropagation step up to the k-th layer, also called backpropagated error. We report the full derivation in appendix A. We adopt the convention of defining x, z k and δ k as column vectors.

Difficulty of optimizing the Jacobian term of neural networks

In the case of the objective function specified in Eq. ( 3), we have L(x) = log p θ (x). By defining

L p (x) = i log p i (g i θ (x)); L J (x) = log |det Jg θ (x)| , (7)
the objective can be rewritten as L(x) = L p (x) + L J (x). The evaluation of the gradient of the first term L p can be performed easily if a simple form for the latent density is chosen, as it only requires simple operations on top of a single forward pass of the neural network. Given that the loss is a scalar, as backpropagation is an instance of reverse mode differentiation [4], backpropagating the error relative to it in order to evaluate the gradients does not increase the overall complexity with respect to the forward pass alone.

In contrast, the evaluation of the gradient of the second term, L J , is very problematic, and our main concern in this paper. The key computational bottleneck is in fact given by the evaluation of the Jacobian during the forward pass. Since the Jacobian involves derivatives of the function g θ with respect to its inputs x, this evaluation can again be performed through automatic differentiation.

Overall, it can be shown [4] that both forward and backward mode automatic differentiation for a L-layer, fully connected neural network scale as O(LD 3 ), with L the number of layers. This is prohibitive in many practical applications with a large data dimension D.

Normalizing flows with simple Jacobians An approach to alleviate the computational cost of this operation is to deploy special neural network architectures for which the evaluation of L J is trivial. For example, in autoregressive normalizing flows [14,15,34,25] the Jacobian of the transformation is constrained to be lower triangular. In this case, its determinant can be trivially computed with a linear cost in D. Notice however that the computational cost of the forward pass still scales quadratically in D; the overall complexity of forward plus backward pass is therefore still quadratic in the input size [44].

Most critically, such architectures imply a strong restriction on the class of transformations that can be learned. While it can be shown, based on [29], that under certain conditions this class of functions has universal approximation capacity for densities [25], that is less general than other notions of universal approximation [23,24]. In fact it is obvious that functions with such triangular Jacobians cannot be universal approximators of functions, since, for example, the first variable can only depend on the first variable. This is a severe problem in learning features for disentanglement, for example by nonlinear ICA [27,30], which would usually require unconstrained Jacobians. In other words, such restrictions might imply that the deployed networks are not general purpose: [5] showed that constrained designs typically used for density estimation can severely hurt discriminative performance. We further elaborate on this point in appendix E. Note that fully connected modules have elsewhere been termed linear flows [42], and are a strict generalization of autoregressive flows. 43 Log-determinant of the Jacobian for fully connected neural networks

As a first step toward efficient optimization of the L J term, we next provide the explicit form of the Jacobian for fully connected neural networks. As a starting point, notice that invertible and differentiable transformations are composable; given any two such transformations, their composition is also invertible and differentiable. Furthermore, the determinant of the Jacobian of a composition of functions is given by the product of the determinants of the Jacobians of each function,

det J[g 2 • g 1 ](x) = det Jg 2 (g 1 (x)) • det Jg 1 (x) . (8)
The log-determinant of the full Jacobian for a neural network therefore simply decomposes in a sum of the log-determinants of the Jacobians of each module, L J (x) = L k=1 log | det Jg k (z k-1 )|. We will focus on the Jacobian term relative to a single submodule k with respect to its input z k-1 ; with a slight abuse of notation, we will call it L J (z k-1 ). As we remarked, fully connected g k are themselves compositions of a linear operation and an element-wise invertible nonlinearity; applying the same reasoning, we then have

L J (z k-1 ) = D i=1 log σ (y i k ) + log |det W k | =: L 1 J (y k ) + L 2 J (z k-1 ) . ( 9
)
where y k = W k z k-1 . The first term L 1 J is a sum of univariate functions of single components of the output of the module, and it can be evaluated easily with few additional operations on top of intermediate outputs of a forward pass; gradients with respect to it can be simply computed via backpropagation, not unlike the L p term introduced in section 2.3.

The second term L 2 J however involves a nonlinear function of the determinant of the weight matrix. From matrix calculus, we know that the derivative is equal to

∂ log | det W k | ∂W k = W k -1 . (10)
Therefore, the computation of the gradient relative to such term involves a matrix inversion, with cubic scaling in the input size. 5 For a fully connected neural network of L layers, given that we have one such operation to perform for each of the layers, the gradient computation for these terms alone would have a complexity of O(LD 3 ), thus matching the one which would be obtained if the Jacobian were to be computed via automatic differentiation as discussed in section 2.

It can therefore be seen that these inverses of the weight matrices are the problematic element in the gradient computation. In the next section, we show how this problem can be solved using relative gradients.

Relative gradient descent for neural networks

We now derive the basic form of the relative gradient, following the approach in [11]. 6 The starting point is that the parameters in a neural networks are matrices, in particular invertible in our case. Thus, we can make use of the geometric properties of invertible matrices, while they are usually completely neglected in gradient optimization in neural networks.

Relative gradient based on multiplicative perturbation In a classical gradient approach for optimization, we add a small vector to a point x in a Euclidean space. However, with matrices, we are actually perturbing a matrix with another, and this can be done in different ways. In the relative gradient approach, we make a multiplicative perturbation of the form

W k → (I + )W k (11)
where is an infinitesimal matrix. If we consider the effect of such a perturbation on a scalar-valued function f (W k ), we have

f ((I + )W k ) -f (W) = ∇f (W k ), W k + o(W k ) = ∇f (W k )W k , + o(W k ) (12
) which shows that the direction of steepest descent in this case is given by making = µ∇f (W k )W k where µ is an infinitesimal step size. Furthermore, when we combine this with the definition of a multiplicative update, we find that the best perturbation to W is actually given as

W k → W k + µ∇f (W k )W k W k (13)
That is, the classical Euclidean gradient is replaced by ∇f (W k )W k W k , i.e. it is multiplied by W k W k from the right. This is the relative gradient.

A further alternative can be obtained by perturbing the weight matrices from the right, as W k → W k (I + ). A similar derivation shows that in this case, the optimal is given by W k W k ∇f (W k ); we refer to this as transposed relative gradient. In the context of linear ICA, the properties of the relative and transposed relative gradient were discussed in [49]. This version of the relative gradient might be useful in some cases; for example, the transposed relative gradient can be implemented more straightforwardly in neural network packages where the convention is that vectors are represented as rows.

The relative gradient belongs to the more general class of gradient descent algorithms on Riemannian manifolds [1]. Specifically, relative gradient descent is a first order optimization algorithm on the manifold of invertible D × D matrices. Almost sure convergence of the parameters to a critical point of the gradient of the cost function can be derived even for its stochastic counterpart, with decreasing step size and under suitable assumptions (see e.g. [8]).

Jacobian term optimization through the relative gradient In section 3, we showed that the difficulty in computing the gradient of the log-determinant is in the terms L 2 J , whose gradient involves a matrix inversion. Now we show that by exploiting the relative gradient, this matrix inversion vanishes. In fact, when multiplying the right hand side of equation ( 10) by W k W k from the right we get

W k -1 W k W k = W k , (14)
and similarly when multiplying by W k W k from the left. Most notably, we therefore have to perform no additional operation to get the relative gradient with respect to this term of the loss; it is, so to say, implicitly computed -as we know that the update for the parameters in W k with respect to the error term L 2 J is proportional to W k matrix itself. As for the remaining terms of the loss, L p and L 1 J , simple backpropagation allows us to compute the weight updates given by the ordinary gradient in equation ( 6), which still need to be multiplied by W k W k to turn it into a relative gradient. We will next see that we can do this avoiding matrix-matrix multiplications, which would be computationally expensive. Note that backpropagation necessarily computes the δ k vector in equation ( 6) and for our model, by applying the relative gradient carefully, we can avoid matrix-matrix multiplication altogether by computing

(∆W k ) W k W k ∝ z k-1 δ k W k W k . (15)
Thus, we have a cheap method for computing the gradient of the log-determinant of the Jacobian, and of our original objective function. In appendix D we provide an explanation of how our procedure can be implemented with relative ease on top of existing deep learning packages.

While we so far only discussed update rules for the weight matrices of the neural network, our approach can be extended to include biases. Including bias terms in our multilayer network endows it with stronger approximation capacity. We detail how to do this in appendix F.

Complexity Note that the parentheses in equation ( 15) stress the point that the relative gradient updates only require matrix-vector or vector-vector multiplications, each of which scales as O(D 2 ), in a fixed number at each layer; that is, overall O(LD 2 ) operations. They therefore do not increase the complexity of a normal forward pass. Furthermore, the overall complexity with respect to the input size is quadratic, resulting in an overall quadratic scaling with the input size as in normalizing flow methods [44], but without imposing strong restrictions on the Jacobian of the transformation.

Extension to convolutional layers As we remarked in section 2.2, the formalism we introduced includes convolutional neural networks (CNNs) [36]. A natural question is therefore whether our approach can be extended to that case. The first natural question pertains the invertibility of convolutional neural networks; the convolution operation was shown [39] to be invertible under mild conditions (see appendix G), and the standard pooling operation can be by replaced an invertible operation [31]. We therefore believe that the general formalism can be applied to CNNs; this would require the derivation of the relative gradient for tensors. We believe that this should be possible but leave it for future work.

Invertibility and generation Given that invertible and differentiable transformations are composable, as discussed in section 3, invertibility of our learned transformation is guaranteed as long as the weight matrices and the element-wise nonlinearities are invertible. Square and randomly initialized (e.g. with uniform or normally distributed entries) weight matrices are known to be invertible with probability one; invertibility of the weight matrices throughout the training is guaranteed by the fact that the L 2 J terms would go to minus infinity for singular matrices (though high learning rates and numerical instabilities might compromise it in practice), as in estimation methods for linear ICA [6,11,26]. We additionally employ nonlinearities which are invertible by construction; we include more details about this in appendix H. If we are interested in data generation, we also need to invert the learned function. In practice, the cost of inverting each of the matrices is O(D 3 ), but the operation needs to be performed only once. As for the nonlinear transformation, the inversion is cheap since we only need to numerically invert a scalar function, for which often a closed form is available.

Experiments

In the following we experimentally verify the computational advantage of the relative gradient. The code used for our experiments can be found at https://github.com/fissoreg/ relative-gradient-jacobian.

Computation of relative vs. ordinary gradient As a first step, we empirically verify that our proposed procedure using the formulas in section 4 leads to a significant speed-up in computation of the gradient of the Jacobian term. We compare the relative gradient against an explicit computation of the ordinary gradient, as described in section 3, and with a computation based on automatic differentiation, as discussed in section 2.3, where the Jacobian is computed with the JAX package [10]. While the output and asymptotic computational complexity of the ordinary gradient and automatic differentiation methods should be the same, a discrepancy is to be expected at finite dimensionality due to differences in how the computation is implemented. In the experiment, we generate 100 random normally distributed datapoints and vary the dimensionality of the data from 10 to beyond 20,000. We then define a two-layer neural network and evaluate the gradient of the Jacobian. The main comparison is run on a Tesla P100 Nvidia GPU. For the main plots, we deactivated garbage collection. Plots with CPU and further details on garbage collection can be found in appendix H.1. For each dimension we computed 10 iterations with a batch size of 100. Results are shown in figure 1, left. On the y-axis we report the average of the execution times of 100 successive gradient evaluations (forward plus backward pass in the automatic differentiation case). It can be clearly seen that the relative gradient is much faster, typically by two orders of magnitude. Autodiff computations could actually only be performed for the smallest dimension due to a memory problem. We report additional details on memory consumption in appendix H.1.

Optimization by relative vs. ordinary gradient Since our paper is, to the best of our knowledge, the first one proposing relative gradient optimization for neural networks (though other kinds of natural gradients have been studied [2]), we want to verify that the learning dynamics induced by the 7 D. Relative gradient optimization of the Jacobian term in unsupervised deep learning relative as opposed to the ordinary gradient do not bias the training procedure towards less optimal solutions or create other problems. We therefore perform a deterministic (full batch) gradient descent for both the relative and the ordinary gradient. 7 We employ 1,000 datapoints of dimensionality 2 and a two-layer neural network. We take 10 initial points and initialize both kinds of gradient descent at those same points. On the x-axis we plot the training epoch, while on the y-axis we plot the value of the loss. Figure 1, right shows the results: there is no big difference between the two gradient methods. There may actually be a slight advantage for the relative gradient, but that is immaterial since our main point here is merely to show that the relative gradient does not need more iterations to give the same performance.

Combining these two results, we see that the proposed relative gradient approach leads to a much faster optimization than the ordinary gradient. Perhaps surprisingly, the results exhibit a rather constant speed-up factor of the order of 100 although the theory says it should be changing with the dimension D; in any case, the difference is very significant in practice.

Density estimation Although our main contribution is the computational speed-up of the gradient computation demonstrated above, we further show some simple results on density estimation to highlight the potential of the relative gradient used in conjuction with the unconstrained factorial approximation in section 2.1. We use a fairly simple feedforward neural network with a smooth version of leaky-ReLU as activation function. Our empirical results show that this system, despite having quite minimal fine-tuning (details in appendix H.3), achieves competitive results on all the considered datasets compared with existing models-which are all tailored and fine-tuned for density estimation. First, we show in Figure 2 different toy examples that showcase the ability of our method to convincingly model arbitrarily complex densities. Second, in order to show the viability of our method in comparison with well-established methods we perform, as in [43], unconditional density estimation on four different UCI datasets [16] and a dataset of natural image patches (BSDS300) [41], as well as on MNIST [37]. The results are shown in Table 1. To achieve a fair comparison across models, the number of parameters was tuned so that the number of trainable parameters are as similar as possible. Note that, as we can perform every computation efficiently, all the experiments are suitable to run on usual hardware, thus avoiding the need of hardware accelerators such as GPUs. As a final remark, the reported results make no use of batch normalization, dropout, or learning-rate scheduling. Therefore, it is sensible to expect even better results by including them in future work.

Table 1: Test log-likelihoods (higher is better) on unconditional density estimation for different datasets and models (same as in Table 1 of [43]). Models use a similar number of parameters; results show mean and two standard deviations. Best performing models are in bold. More details in appendix H. 3

Conclusions

Using relative gradients, we proposed a new method for exact optimization of objective functions involving the log-determinant of the Jacobian of a neural network, as typically found in density estimation, nonlinear ICA, and related tasks. This allows for employing models which, unlike typical alternatives in the normalizing flows literature, have no strong limitation on the structure of the Jacobian. We use modules with fully connected layers, thus strictly generalizing normalizing flows with triangular Jacobians, while still supporting efficient combination of forward and backward pass. These neural networks can represent a larger function class than autoregressive flows, which, despite being universal approximators for density functions, can only represent transformations with triangular Jacobians. Our method can therefore provide an alternative in settings where more expressiveness is needed to learn a proper inverse transformation, such as in identifiable nonlinear ICA models.

The relative gradient approach proposed here is quite simple, yet rather powerful. The importance of the optimization of the log-determinant of the Jacobian is well-known, but it has not been previously shown that there is a way around its difficulty without restricting expressivity. Now that we have shown that the optimization of this term can be done quite cheaply, a substantial fraction of the research in the field can be reformulated in stronger terms and with more generality.

Broader impact

As this paper presents novel theoretical results in unsupervised learning, the authors do not see any immediate ethical or societal concern. An important aspect of our paper is the improvement in computational efficiency with respect to naive methods. This can hopefully lead to reduced energy consumption to achieve comparable model performance.

i 1 1 . 3 1 . 4 1 . 5 18 2

 113141518 Introduction v A brief history of Artificial Intelligence . . . . . . . . . . . . . . . . . v Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii The mathematical neuron . . . . . . . . . . . . . . . . . . . . . 3 1.2 The Hopfield model . . . . . . . . . . . . . . . . . . . . . . . . 4 Statistical description of a neural network . . . . . . . . . . . . 7 Free energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Order parameters and phase transitions . . . . . . . . . . . . . 10 1.6 Effective free energy . . . . . . . . . . . . . . . . . . . . . . . . 11 1.7 Mean-field theory of the Hopfield model . . . . . . . . . . . . . 12 1.8 Beyond memory: Boltzmann Machines . . . . . . . . . . . . . . 16 1.9 Statistical Physics approach to neural networks . . . . . . . . . Restricted Boltzmann Machines (RBM) 19 2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Contrastive divergence training . . . . . . . . . . . . . . . . . . 21 2.3 Mean-field training . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4 Generalized RBM . . . . . . . . . . . . . . . . . . . . . . . . . . 24

1. 7 .Figure 1 . 1 :

 711 Figure 1.1: Phase diagram of the Hopfield model.

Algorithm 2 k 1 : 2 : 4 :

 2124 -steps Contrastive Divergence (CDk) Data: a training set of N data vectors v i Randomly initialize the weight matrix W 3: for t = 0 to T (# of epochs) do Divide the training set in m minibatches of n data vectors v i 5:

Algorithm 3 k 1 : 2 : 4 :

 3124 -steps Extended mean-field training (EMFk) Data: a training set of N data vectors v i Randomly initialize the weight matrix W 3: for t = 0 to T (# of epochs) do Divide the training set in m minibatches of n data vectors v i 5:

. 14 )Figure 3 . 1 :

 1431 Figure 3.1: Samples of handwritten digits from the MNIST dataset.

Figure 3 . 2 :

 32 Figure 3.2: (a) Singular values distribution of the initial random matrix compared to Marchenko-Pastur law. (b)-(d) With the training we can see some singular values strengthening and overcoming the threshold set by the Marchenko-Pastur law. (e) Distribution of the singular values after a long training: we can see many outliers spread above threshold and a spike of belowthreshold singular values near zero.

3. 5 .Figure 3 . 3 :

 533 Figure 3.3: (a) First mode learnt by the RBM with the external visible field initialized as a null vector. (b) External visible field initialized with rule (3.19). (c) First principal components extracted from the training set. (d) Principal components extracted from the training set (starting from the second). (e) The first 10 modes of a RBM trained for 1 epoch, using initialization rule (3.19). (f ) Same as (e) but after a 10 epochs training.

Figure 3 . 4 :

 34 Figure 3.4: Log-log plot of the singular values represented as discrete abscissas (in decreasing order) with their magnitude reported on the ordinates. A cutoff is highlighted by the onset of the linear behavior.

Figure 3 . 5 :

 35 Figure 3.5: (a) The image shows some samples obtained with the trained RBM (after a 40 epochs training) and then "filtered" by eliminating the 400 weakest modes (just the 100 strongest modes are retained). (b) The images are composed by eliminating the 100 strongest modes to see what the weakest modes actually encode (20 epochs training). (c) As in (b) but after a 40 epochs training.

Figure 4 . 1 :

 41 Figure 4.1: Phase diagram in absence of bias and with a finite number of modes, with Gaussian and Laplace distributions for u and v. Inset: high temperature (σ = 0) stability gap ∆w α corresponding to a fixed point associated to a mode β, expressed as a function of w α and considering various distributions.Reproduced from[20].

  2 and 4.4 to compute the . . . Data (empirical) and . . . RBM (model) averages in Equations (4.27)-(4.29). In the empirical average the value of the nodes of the network is driven by the data, so that we can directly substitute the nodes values in Equations (4.21)(4.22) to obtain

  Figure 4.3. For the singular values we

Figure 4 . 3 :

 43 Figure 4.3: Predicted mean evolution of an RBM for a synthetic dataset, with (N v , N h ) = (1000, 500) and intrinsic dimension d = 15. The singular values evolution is obtained by integration of Equation (4.54).

Figure 4 . 4 :

 44 Figure 4.4: Empirical trajectories of the RBM weight matrix in the phase diagram.

Figure 4 . 5 :

 45 Figure 4.5: Empirical evolution of the RBM for a synthetic dataset with (N v , N h ) = (1000, 500) and intrinsic dimension d = 15. The behavior is comparable to the predicted evolution shown in Figure 4.3.

Figure 4 . 6 :

 46 Figure 4.6: Empirical evolution of the RBM for the MNIST dataset.

3 Figure 4 . 7 :

 347 Figure 4.7: Projections of the mean-field magnetizations (in red) and samples from a synthetic dataset over the planes defined by the strongest spectral directions x 1 , x 2 , . . . x 5 .

Figure 4 . 8 :

 48 Figure 4.8: Empirical projections of the TAP fixed points of a trained RBM (in red) and samples from a synthetic dataset over the planes defined by the strongest spectral directions x 1 , x 2 , . . . x 5 .

4 :

 4 The full training algorithm becomes: Algorithm 4 Lossy-CDk (RBM training with incomplete data) 1: Data: a training set of N data vectors 2: Randomly initialize the weight matrix W 3: for t = 0 to T (# of epochs) do Divide the training set in m minibatches 5:

10 : initialize v randomly 11 :

 1011 sample h, v using p(v | h) and p(h | v) for k steps 12: compute vh T model Full update: 13:update W with equations (2.8) and (5.4)-(5.6)

Figure 5 . 1 :

 51 Figure 5.1: Features reconstruction by RBM-MO trained over an incomplete dataset with 50% missing-at-random features, whose classification accuracy has been measured to be around 91%. The first block shows some complete testing instances. The second and third block show the same testing instances after hiding respectively 50% and 80% of the pixels. The last two columns show the results of the mean-field imputations over the incomplete testing instances.

64

 64

6. 2 .Figure 6 . 1 :

 261 Figure 6.1: Example of a 4-step flow transforming samples from a standardnormal base density to a cross-shaped target density. Reproduced from [68].

Figure 6 . 2 :

 62 Figure 6.2: Comparison of the average computation times of a single evaluation of the gradient of the log-likelihood over a batch of 100 samples; the standard error of the mean is not reported as it is orders of magnitude smaller than the scale of the plot. We set the number of layers to L = 2 and performed the experiment using a Tesla P100 Nvidia GPU. Reproduced from[35].

Figure 6 . 3 :

 63 Figure 6.3: Illustrative examples of 2-dimensional density estimation. Samples from the true distribution and predicted densities are shown, in this order, side by side.

p- 3 A

 3 Fig.1: Time evolution of the eigenvalues in the linear model and of the likelihood. We observe very clearly how the different modes emerge from the bulk and how the likelihood increases at each eigenvalue learned. In the inset, the scalar product of the vectors u obtained from the SVD of the data and of w. The us of w are aligned with the SVD of the data at the end of the learning.

Fig. 2 :

 2 Fig. 2: (a) Singular values distribution of the initial random matrix compared to Marchenko-Pastur law. (b) With the training we can see some singular values strengthening and overcoming the threshold set by the Marchenko-Pastur law. (c) Distribution of the singular values after a long training: we can see many outliers spread above threshold and a spike of below-threshold singular values near zero.

Fig. 3 :

 3 Fig. 3: (a) First mode learnt by the RBM with the external visible field initialized as a null vector. (b) External visible field initialized on the empirical mean. (c) First principal components extracted from the training set. (d) Principal components extracted from the training set (starting from the second). (e) The first 10 modes of a RBM trained for 1 epoch. (f ) Same as (e) but after a 10 epochs training.

  10, 100, 350 and 400 (top to bottom)

2

  The RBM and its associated learning procedure 4 dataset.

First the average over rij yields exp î σ 2 2L.

 2 After this averaging, 4 sets of order parameters {(m a α , ma α ), a = 1, . . . p, α = 1, . . . K} and {(Qab, Qab), a, b = 1, . . . p, a = b} are introduced with the help of two distinct Hubbard-Stratonovich transformations. The first one corresponds to exp

  following set of dynamical equations:

Fig. 3 :

 3 Fig.3: Time evolution of the eigenvalues and of the likelihood in the linear model. We observe very clearly how the different modes emerge from the bulk and how the likelihood increases with each learned eigenvalue.In the inset, the scalar product of the vectors u obtained from the SVD of the data and from the weights w. The us of w are aligned with the SVD of the data at the end of the learning.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Predicted mean evolution of an RBM of size (N v , N h ) = (1000, 500) learned on a synthetic dataset of 10 4 samples of size N v = 1000 obtained from a multimodal distribution with 20 clusters randomly defined on a submanifold of dimension d = 15. The dynamics follows the projected magnetizations in this reduced space with help of 15 modes. We observe a kind of pressure on top singular values from lower ones.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Experimental evolution of an RBM during training for a synthetic dataset (top plots, to compare to Fig. 4) and for MNIST (central plots). The bottom left plot shows the learning trajectories in the phase diagram, while the bottom right image shows some examples of fixed point solutions for MNIST (we note the presence of some spurious fixed points).

Fig. 8 :

 8 Fig. 8: Scatter plots of samples (blue) and fixed points (red) in various plan projections defined by pairs of left eigenvectors of W. The dataset is the same as in Fig. 5 and in this case 5 modes have condensed and 7 fixed point solutions have been found.

Finally, the training

  set is divided into batches of size 20, 5 Gibbs sampling steps are used (CD5) and the learning rate γ is kept low in order to reduce noise, γ = 5 × 10 -8 .

Fig. 9 :

 9 Fig. 9: (a) Principal components extracted from the training set (starting from the second, as the first one is encoded into the visible bias). (b) The first 10 modes of a RBM trained for 1 epoch (with γ 0.1). (c) Same as (b) but after a 10 epochs training.

a 1 -

 1 =b,c,α(δab + δac) Cα + (1δacδbc) Dα c α ηab + CT + a =b,c =d δ (ab)(cd) Ē0 + 1 1 {a∈(cd)⊕b∈(cd)} Ē1 + 1 1 {(ab)∩(cd)=∅} Ē2 ηabηcd + CT ,where CT means "conjugate term" in the sense ↔ ¯ , Aαβ ↔ Āαβ. . . , where δab def = δab and the operators are given by

Algorithm 1 : 1 : 4 : 6 : 7 : 8 : 9 :

 1146789 Lossy-CDk (RBM training with Incomplete data) Data: a training set of N data vectors 2: Randomly initialize the weight matrix W 3: for t = 0 to T (# of epochs) do Divide the training set in m minibatches 5: for all minibatches m do Positive term: pin variables v i , i ∈ O to their correct value initialize v i , i ∈ M randomly sample h, v m using p(v m | h) and p(h | v) for k steps compute the positive terms in (5) Negative term: 10: initialize v randomly 11:

12 :

 12 update W with equation(6) 13:

i

  log p i (g i θ (x)) + log | det Jg θ (x)| .

Figure 1 :

 1 Figure1: Left: Comparison of the average computation times of a single evaluation of the gradient of the log-likelihood; the standard error of the mean is not reported as it is orders of magnitude smaller then the scale of the plot. Right: Time-evolution of the negative log-likelihood for deterministic full-batch optimization for the two methods with the same initial points.

Figure 2 :

 2 Figure 2: Illustrative examples of 2D density estimation. Samples from the true distribution and predicted densities are shown, in this order, side by side.

  Figure 4.2: Linear RBM. Reproduced from Reprint A. components of the data. In turn, the diagonal evolution expressed in Equation (4.45) involves only the singular values. The stable solutions of this diagonal dynamics are

	4. Mean-field theory
	5 10 15 20 EVs	1 0	0	1st EV 2nd EV 3rd EV 4th EV Likelihood u X .u w 3000 6000 t	Likelihood -A.U.
	0	0 1000 2000 3000 4000 5000 6000 7000 8000
				t

Table 1 :

 1 1 if i = argmax k (p k ) Transductive test on MNIST multi-class data set (our method in bold, best result in red)

	Model		RMSE		Averaged AUC		Accuracy
	qmc%	30%	50%	80%	30%	50%	80%	30%	50%	80%
	RBM-MO(qfea% =50%) 0.183 0.182 0.185 0.969 0.971 0.929 0.950 0.912 0.822
	CLE(qfea% =50%)	0.195	0.195	0.195	0.686	0.718	0.742	0.256	0.232	0.282
	NoisyIMC(qfea% =50%)	0.209	0.210	0.210	0.621	0.578	0.552	0.225	0.232	0.192
	MC-1(qfea% =50%)	0.334	0.335	0.337	0.495	0.493	0.500	0.110	0.111	0.112
	RBM-MO(qfea% =80%) 0.209 0.213 0.211 0.938 0.932 0.906 0.920 0.852 0.733
	CLE(qfea% =80%)	0.206	0.208	0.206	0.673	0.678	0.625	0.230	0.215	0.220
	NoisyIMC(qfea% =80%)	0.212	0.211	0.213	0.652	0.577	0.537	0.230	0.217	0.210
	MC-1(qfea% =80%)	0.334	0.334	0.335	0.500	0.501	0.500	0.112	0.110	0.110
	Model		RMSE			Micro-AUC		Hamming-Accuracy
	qml%	30%	50%	80%	30%	50%	80%	30%	50%	80%
	RBM-MO(qfea% =50%) 0.131 0.137 0.123 0.943 0.934 0.888 0.919 0.907 0.873
	CLE(qfea% =50%)	0.130	0.130	0.131	0.905	0.893	0.898	0.885	0.871	0.878
	NoisyIMC(qfea% =50%)	0.132	0.133	0.133	0.865	0.863	0.858	0.845	0.841	0.848
	MC-1(qfea% =50%)	0.258	0.255	0.267	0.522	0.528	0.527	0.826	0.817	0.824
	RBM-MO(qfea% =80%) 0.160 0.158 0.158 0.875 0.867 0.826 0.856 0.858 0.832
	CLE(qfea% =80%)	0.129	0.129	0.128	0.913	0.897	0.899	0.889	0.875	0.876
	NoisyIMC(qfea% =80%)	0.133	0.134	0.134	0.853	0.857	0.849	0.839	0.835	0.826

Table 2 :

 2 Transductive test on Scene multi-label data set (our method in bold, best result in red)

Table 3 :

 3 Inductive test on Pendigits multi-class dataset (our method in bold, best result in red) C. Robust Multi-Output Learning with Highly Incomplete Data via Restricted Boltzmann Machines

	Model		Micro-AUC		Hamming-Accuracy
	qml%	30%	50%	80%	30%	50%	80%
	RBM-MO(qfea% =50%) 0.839 0.826 0.793 0.970 0.970 0.965
	CLE(qfea% =50%)	0.705	0.707	0.706	0.700	0.724	0.719
	NoisyIMC(qfea% =50%)	0.704	0.702	0.700	0.710	0.717	0.718
	RBM-MO(qfea% =80%) 0.759 0.791 0.766 0.964 0.964 0.967
	CLE(qfea% =80%)	0.693	0.688	0.694	0.718	0.706	0.718
	NoisyIMC(qfea% =80%)	0.689	0.688	0.685	0.705	0.704	0.704

Table 4 :

 4 Inductive test on EventCat multi-label data set (our method in bold, best result in red)

	Dataset No. of Instances No. of Features No.of Labels No. of Classes
	Scene	2,407	294	6	-
	Pendigits	10992	16	-	10
	MNIST	70,000	784	-	10
	EventCat	5,93	72	6	-

Table 5 :

 5 Summary of 4 public multi-label and multi-class data sets.

The separation of the weight matrix spectrum into a structured part and a random part has been observed over different real-world datasets and for different configurations of the trained RBM model. While this separation might not always hold, it is a relevant case in practice.

Here we refer to the n-dimensional Euclidean space and the real space R n interchangeably; it is understood that we work in the real space R n equipped with Euclidean structure and a Cartesian coordinate system.

([35], Reprint D) calls the update term in (6.13) the relative gradient. Here we align to the original nomenclature by[16].

The derivative to compute δ k is written in numerator layout notation; the derivatives w.r.t. the matrix W k are naturally arranged in the column vector z k-1 .

a somewhat different form of the TAP equations

The forward transformation could also be parameterized, but here we only explicitly parameterize its inverse.2

Note that invertible neural networks provide the possibility to not save, but rather recompute the intermediate activations during the backward pass, thus providing a memory efficient approach to backpropagation[18].

Comprehensive reviews on normalizing flows can be found in[42, 35]. Other related methods are reviewed in appendix B.4

Though slightly more favorable exponents can in principle be obtained, see appendix C.

For linear blind source separation, this approach also corresponds to the natural gradient, which can be justified with an information-geometric approach[2].5D. Relative gradient optimization of the Jacobian term in unsupervised deep learning

Notice that there's no need to compare to autodiff in this case because the computed gradient should be exactly the same as the ordinary gradient with the formulas in section 3.

We thank B. Seoane helping us improving the quality of the letter.

Acknowledgments

A.H. was supported by a Fellowship from CIFAR, and by the DATAIA convergence institute as part of the "Programme d'Investissement d'Avenir", (ANR-17-CONV-0003) operated by Inria. L.G. started working on this project while on an ELLIS exchange, hosted by the Parietal team at Inria, Saclay. We thank Vincent Stimper, Patrick Putzky, Cyril Furtlehner and Ilyes Khemakhem for valuable comments on an earlier draft of this paper, and Isabel Valera and Guillaume Charpiat for helpful comments and tips. L.G. additionally thanks Roma Beaufret and Alexis Bozio for helpful support.

Reprints

Abstract

We investigate the thermodynamic properties of a Restricted Boltzmann Machine (RBM), a simple energy-based generative model used in the context of unsupervised learning. Assuming the information content of this model to be mainly reflected by the spectral properties of its weight matrix W , we try to make a realistic analysis by averaging over an appropriate statistical ensemble of RBMs. First, a phase diagram is derived. Otherwise similar to that of the Sherrington-Kirkpatrick (SK) model with ferromagnetic couplings, the RBM's phase diagram presents a ferromagnetic phase which may or may not be of compositional type depending on the kurtosis of the distribution of the components of the singular vectors of W .

Subsequently, the learning dynamics of the RBM is studied in the thermodynamic limit. A "typical" learning trajectory is shown to solve an effective dynamical equation, based on the aforementioned ensemble average and explicitly involving order parameters obtained from the thermodynamic analysis. In particular, this let us show how the evolution of the dominant singular values of W , and thus of the unstable modes, is driven by the input data. At the beginning of the training, in which the RBM is found to operate in the linear regime, the unstable modes reflect the dominant covariance modes of the data. In the non-linear regime, instead, the selected modes interact and eventually impose a matching of the order parameters to their empirical counterparts estimated from the data.

Finally, we illustrate our considerations by performing experiments on both artificial and real data, showing in particular how the RBM operates in the ferromagnetic compositional phase.

Introduction

The Restricted Boltzmann Machine (RBM) [1] is an important machine learning tool used in many applications, by virtue of its ability to model complex probability distributions. It is a neural network which serves as a generative model, in the sense that it is able to approximate the probability distribution corresponding to the empirical distribution of any set of high-dimensional data points living in a discrete or real space of dimension N 1. From the theoretical point of view, the RBM is of high interest as it is one of the simplest neural network generative models and the probability distribution that it defines presents a simple analytic form. Moreover, there are clear

Abstract

Learning expressive probabilistic models correctly describing the data is a ubiquitous problem in machine learning. A popular approach for solving it is mapping the observations into a representation space with a simple joint distribution, which can typically be written as a product of its marginals -thus drawing a connection with the field of nonlinear independent component analysis. Deep density models have been widely used for this task, but their maximum likelihood based training requires estimating the log-determinant of the Jacobian and is computationally expensive, thus imposing a trade-off between computation and expressive power. In this work, we propose a new approach for exact training of such neural networks. Based on relative gradients, we exploit the matrix structure of neural network parameters to compute updates efficiently even in high-dimensional spaces; the computational cost of the training is quadratic in the input size, in contrast with the cubic scaling of naive approaches. This allows fast training with objective functions involving the log-determinant of the Jacobian, without imposing constraints on its structure, in stark contrast to autoregressive normalizing flows.

Introduction

Many problems of machine learning and statistics involve learning invertible transformations of complex, multimodal probability distributions into simple ones. One example is density estimation through latent variable models under a specified base distribution [51], which can also have applications in data generation [14,33,19] and variational inference [44]. Another example is nonlinear independent component analysis (nonlinear ICA), where we want to extract simple, disentangled features out of the observed data [27,30].

One approach to learn such transformations, introduced in [50] in the context of density estimation, is to represent them as a composition of simple maps, the sequential application of which enables high expressivity and a large class of representable transformations. Deep neural networks parameterize functions of multivariate variables as modular sequences of linear transformations and componentwise activation functions, thus providing a natural framework for implementing that idea, as already proposed in [45]. *