Acknowledgments↩︎

I would like to thank a few people without whom this thesis would not have been possible. My PhD supervisor, Daniele, for guiding me throughout the thesis with impeccable patience and rigor while also giving me a lot of freedom to choose the topics. My PhD coordinator, Michele, for kindly advising me and overseeing my progress. My colleague Francesco for our fruitful collaboration. My partner, Daniela, for supporting me unconditionally and encouraging me to go out of my comfort zone. My parents, Lyne and Charles, for fostering my curiosity and encouraging me to pursue my goals. My brother, Lucas, for keeping in touch and playing games with me despite the physical distance and the time constraints. My McGill supervisor, Paul, for having introduced me to scientific research.

1 Introduction↩︎

=[thick, draw = black, circle, minimum size = 22]

1.1 Background↩︎

Artificial Neural Networks (NNs) and Machine Learning (ML) algorithms have, with stunning speed, gone from niche tools to being the state-of-the-art in solving numerous fundamental and practical problems. To name a few of their achievements, they have long surpassed humans in image recognition tasks dan2011committee?, they are able to emulate written language so well that it is possible to have an open-ended conversation with them openai2024gpt4?, and they can predict the structure of complex molecules such as proteins jumper2021highly?, Abramson2024accurate?, which has promising applications in drug discovery and in medicine as a whole yang2023alpha?. In the process, they have also become much larger and more sophisticated, to the point where even models with hundreds of millions of parameters such as latent diffusion Rombach2022high? are considered compact by today’s standards. This level of complexity makes their behavior very challenging to explain theoretically. However, the general idea behind them is relatively simple.

In essence, Machine Learning (ML) is a set of computer algorithms for fitting high-dimensional data \(\mathbf{x}\) with a function \(F_{\mathbf{W}}\) by adjusting its parameters \(\mathbf{W}\), which are also known as the function’s weights. This process is commonly refer to as training or learning the weights. Here, the term high-dimensional data means data with many components \(x_1, x_2, ..., x_N\), which are usually concrete features of the data. For example, digital grayscale images are a type of high-dimensional data whose components are the intensities of their pixels. An artificial Neural Network (NN) is a function \(F_{\mathbf{W}}\) whose weights \(\mathbf{W}\) can be interpreted as the strengths of connections in a graph or network. For example, the linear model \(y = \sum_{i = 1}^N w_i x_i\) of \(\mathbf{x} = \left( x_1, x_2, ..., x_N \right)\) is an NN whose weights \(w_1, w_2, ..., w_N\) can be “machine learned” from a set of high-dimensional data points \(\left( \mathbf{x}^1, y^1 \right), \left( \mathbf{x}^2, y^2 \right), ..., \left( \mathbf{x}^M, y^M \right)\) using linear regression (see Fig. 1). Intuitively, each connection \(i\) is like a pipe with a width of \(w_i\) that controls how much of the input component \(x_i\) flows to the output component \(y\) (see Fig. 2, left panel). In general, NNs have a more sophisticated network of connections, some of which perform non-linear operations (see Fig. 2, right panel), which is what makes them much more complicated and powerful than the linear model HORNIK1989feedforward?. The nodes at the junctions of these NN connections are typically called neurons or units.

ML algorithms are particularly useful for learning probability distribution functions that describe how data points are placed relative to each other. The goal of such models is usually to generate new data points that are close enough to the existing ones to appear realistic, such as coherent sentences in the case of (large) language models openai2024gpt4?. Simpler NNs are also used to predict a noisy output from a deterministic input, such as for categorizing an image into one of many classes lecun1998gradient?. Learning a generic function and a probability distribution are two sides of the same coin. For example, linear regression is arguably equivalent to fitting data \(\left( \mathbf{x}, y_{\text{data}} \right)\) with the “deterministic” distribution \(p \left( y \mid \mathbf{x} \right) = \delta \left( y - \sum_{i = 1}^N w_i x_i \right)\) such that \(y\) is equal to \(\sum_{i = 1}^N w_i x_i\) with probability \(1\) and cannot take any other values (see Fig. 1, right panel). On the flip side, we can also fit the data with a proper probability distribution such as the Gaussian distribution \(p \left( y \mid \mathbf{x} \right) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{1}{2 \sigma^2} \left[ y - \sum_{i = 1}^N w_i x_i \right]^2 \right)\) of \(y\) given \(\mathbf{x}\) (illustrated in Fig. 3), then use the resulting \(w_i\) as the weights of the linear fit. Generating data points amounts to sampling the probability distribution fit to the data. In the case of the linear model, it consists of calculating \(y_{\text{model}} = \sum_{i = 1}^N w_i x_i\) for a given \(\mathbf{x}\) and adding some random noise \(\varepsilon\) to the result, for example according to the Gaussian distribution mentioned above (see Figs. 3 and 1). In this context, the variance \(\sigma^2\) represents the level of noise injected into \(y\). When it is sufficiently small, the generated values of \(y\) are practically indistinguishable from the \(y_{\text{model}}\) obtained without adding noise.

Due to the difficulty in studying them theoretically and how quickly they have evolved, ML and NNs are still poorly understood, making it difficult to fully realize their potential and overcome their weaknesses. On the one hand, it was observed that ML algorithms train NN weights by moving them along a low-dimensional subspace, or manifold, of their allowed values frankle2018lottery?, mao2024training?, which could potentially be used to design more efficient algorithms if it were better understood. To use an analogy, it is more efficient to locate objects on the surface of the Earth with the two-dimensional latitude-longitude coordinate system than with the three-dimensional Cartesian coordinate system. For simplicity, we call this phenomenon implicitly low-dimensional learning for the rest of the Introduction.

On the other hand, NNs that classify data with high accuracy can often be fooled by modifying the data with carefully crafted perturbations that are either meaningless or completely invisible to humans biggio2013evasion?, szegedy2013intriguing? (see Fig. 4 for an example). This data corruption process, called an adversarial attack, can be dangerous in some cases. For example, it can confuse the traffic sign recognition algorithms of self-driving vehicles pavlitska2023adversarial?. Adversarial attacks are not well understood, which makes them hard to circumvent.

In the same way that a bridge model made of ice cream sticks can give a rough idea of the physics of a real bridge, many aspects of ML and NNs can be understood by studying simpler high-dimensional models. For example, goldt2020modelling?, gerace2021generalization? modeled high-dimensional data as lying on a low-dimensional manifold and tanner2024high?, vilucchio2024geometry?, vilucchio2025existence? studied adversarial attacks in linear models. The calculations and theoretical modeling in both of these two lines of work are based on statistical mechanics, and this is not a coincidence. In fact, the field of statistical physics (also known as statistical mechanics) was originally developed to study phase transitions in materials with a large number of particles, but it can also be used to describe data points with many components, datasets with many data points, NNs with many neurons, etc. In particular, when we change their properties and those of their training data, NNs can undergo phase transitions between different regimes of data representation (uninformative, universal, specialized, etc.) barbier2025optimal? analogous to those between the states of matter (gaseous, liquid, solid, etc.). The scientific literature on the topic is rich and flourishing charbonneau2023spin?, recent highlights being the 2024 Nobel Prize awarded to John Hopfield and Geoffrey Hinton for their highly influential work on two models of NNs inspired by statistical physics: the Hopfield networks hopfield1982neural? and restricted Boltzmann machines hinton2002training?, respectively.

As suggested in the previous paragraph, Restricted Boltzmann Machines (RBMs) and Hopfield Networks (HNs) are simple enough to be studied with analytical calculations, yet still exhibit some of the most intriguing characteristics of modern ML and NNs. For example, RBMs and generalizations of HNs called dense Hopfield networks chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, horn1988capacities?, krotov2016dense? can fit data with a rich underlying structure by undergoing a multistage learning process whose first stages were found to be low dimensional decelle2017spectral?, decelle2018thermodynamics?, bachtis2025cascade?, boukacem2024waddington?, a property that was then exploited to accelerate RBM training bereux2025fast?. Moreover, dense HNs were observed to be either vulnerable or robust to adversarial attacks in different controlled scenarios krotov2018dense?, suggesting that it could be possible to improve our understanding of adversarial attacks by studying them.

RBMs and HNs fit data \(\mathbf{x}\) using an energy function \(H\), which is also known as a Hamiltonian. In ML terminology, they are energy-based models. The energy function has weights analogous to those of other NNs, which we will denote as \(\mathbf{J}\), and it is designed so that more realistic data points typically have lower energy than less realistic ones once the weights \(\mathbf{J}\) have been determined. Intuitively, a physical object leaking energy with the environment, such as a ball rolling down a mountain, eventually converges to a configuration where its energy cannot decrease anymore, which in the case of the ball corresponds to standing still at the bottom of a valley. Such a state is called an equilibrium configuration or a local minimum of the energy. The local minima \(\mathbf{x}\) of an energy-based model with given weights \(\mathbf{J}\) can be found analogously, and these states generally correspond to realistic data provided that the model and the strategy used to find \(\mathbf{J}\) are appropriate to represent its underlying structure hopfield1982neural?, hinton2002training?, song2021score?. For example, the minimum energy configurations of HNs are typically preexisting data points stored in the weights by a “Machine Memorization” (MM) procedure (described in Section 1.2), and RBMs can learn to generate new data on which they were not trained.

This thesis studies the regimes of data representation and the phase transitions of RBMs and dense HNs chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, horn1988capacities?, krotov2016dense?, ramsauer2020hopfield? to advance our understanding of MM, ML and NNs, with a focus on implicitly low-dimensional learning and adversarial attacks. Our methods are based on statistical mechanics calculations that are enabled by the relative simplicity of these energy-based models. On the way towards this objective, we investigate deep relationships between different versions of these models and their phase transitions barra2017phase?, decelle2021inverse?, albanese2022replica?, hou2019minimal? that considerably simplify some of the calculations alemanno2023hopfield?, manzan2025effect?.

The rest of this Chapter aims to introduce readers with a background in probability theory and physics to the models that we study and the tools that we use to do so. In Section 1.2, we use HNs as an example to present the general idea behind energy-based models. In Section 1.3, we build upon this foundation to describe RBMs and dense HNs. In Section 1.4, we present the teacher-student setting that we use in our analytical calculations and explain what kind of results can be derived with it. Finally, Section 1.5 uses the notions presented in Sections 1.2, 1.3 and 1.4 to put together a comprehensive outline of the other Chapters of the thesis, which present the results of our research.

Here is a high-level preview of these results. In Chapter 2, we study the adversarially vulnerable and adversarially robust regimes of dense HNs using a formal relationship between MM and ML that we investigate in detail theriault2024dense?. In Chapter 3, we study the various learning strategies of RBMs as a function of the properties of the training data theriault2025modeling?. In particular, we show how RBMs that are larger than necessary to learn the data adopt an implicitly low-dimensional learning strategy. In Chapter 4, we develop a theory of implicitly low-dimensional learning in a kind of RBM that belongs to the class of Dense Associative Memory (DAM), which is a generalization of dense HNs theriault2026saddle?. We then use our theory to greatly accelerate training.

1.2 The Hopfield network: a simple example of energy-based model↩︎

In this Section, we present HNs as simple examples of energy-based models. This exercise is useful for building intuition before tackling RBMs and dense HNs, which are more complicated.

HNs use their units \(\boldsymbol{\sigma} = \left\{ \sigma_i \right\}_{i = 1}^N = \left( \sigma_1, \sigma_2, ..., \sigma_N \right)\) to represent data whose components can take only two values, such as binary black-and-white images (see Fig. 6, left panel). Without loss of generality, we take these two values to be \(+1\) and \(-1\). In this case, the HN Hamiltonian takes the form \[\label{eq:Hopfield95Hamiltonian} H_{\text{HN}} \left[ \boldsymbol{\sigma} ; \mathbf{J} \right] = -\sum_{i = 1}^N \sum_{j = i + 1}^N J_{i j} \sigma_i \sigma_j = -\sum_{i < j = 1}^N J_{i j} \sigma_i \sigma_j,\tag{1}\] where the weights \(J_{i j}\) are the strengths of the connections between different components \(\sigma_i\) and \(\sigma_j\). Contrary to those of the NNs shown in Fig. (2), these connections are bidirectional: they are meant to represent how \(\sigma_i\) and \(\sigma_j\) are correlated rather than how one influences the other (see Fig. 6, right panel). The sums in the Hamiltonian (Eq. 1 ) are restricted to \(j > i\) so that the contribution of every connection to the total energy is counted exactly once.

There are many ways to fit the weights \(J_{i j}\) to data kanter1987associative?, storkey1997increasing?, barra2012equivalence?, serricchio2025daydreaming?, the simplest one being Hebb’s rule hopfield1982neural?, hebb2005organization? \[\begin{align} \label{eq:Hebb} J_{i j} &= R_{i j} \left( \boldsymbol{\xi} \right) = \frac{1}{N} \sum_{\mu = 1}^P \xi^\mu_i \xi^\mu_j, \addtocounter{equation}{1}\end{align}\tag{2}\] where \(i < j\) and each pattern \(\boldsymbol{\xi}^\mu = \left\{ \xi^\mu_i \right\}_{i = 1}^N\) in the set \(\boldsymbol{\xi} = \left\{ \boldsymbol{\xi}^\mu \right\}_{\mu = 1}^P\) is a data point given to the model. This rule establishes strong positive connections between neurons \(\sigma_i\) and \(\sigma_j\) that consistently take the same value (\(\sigma_i = \sigma_j\)), strong negative connections between neurons that consistently take opposite values (\(\sigma_i = -\sigma_j\)) and weak connections between neurons that do not synchronize in either of these ways. It is often summarized by the mantra “neurons that fire together, wire together. Neurons that do not sync, fail to link”.

HNs can retrieve the patterns stored in their weights, thus acting as models of associative memory hopfield1982neural? and making Hebb’s rule a form of “Machine Memorization” (MM). To be more precise, once its weights have been determined, an HN can be driven to local minima of the energy as described in hopfield1982neural?. When the weights are given by Hebb’s rule with a sufficiently small number \(P\) of patterns hopfield1982neural?, amit1985storing?, the units \(\boldsymbol{\sigma}\) usually converge to one of the patterns \(\boldsymbol{\xi}^\mu\) or to a mixture of them amit1987statistical?, amit1987information?. Conversely, when the number \(P\) of patterns exceeds a specific threshold proportional to the number \(N\) of components per data point, the HN typically no longer converges to meaningful results. In other words, it has limited storage capacity hopfield1982neural?, amit1985storing?.

As we argued in Section 1.1, learning a generic function and a probability distribution are two sides of the same coin. This observation suggests that we could learn the patterns \(\boldsymbol{\xi}^\mu\) by fitting a probability distribution of the energy \(H_{\text{HN}} \left[ \boldsymbol{\sigma} ; R \left( \boldsymbol{\xi} \right) \right] = -\frac{1}{N} \sum_{i < j = 1}^N \sum_{\mu = 1}^P \xi^\mu_i \xi^\mu_j \sigma_i \sigma_j\) (see Eqs. 1 and 2 ) to the data, as studied in decelle2021inverse?, alemanno2023hopfield?, barra2012equivalence?, rather than setting the patterns directly equal to the data. For this purpose, we turn to statistical mechanics.

In statistical mechanics, a physical system exchanging energy with the environment often does not converge to a single well-defined equilibrium configuration, but rather to an equilibrium distribution of configurations such that more probable configurations have lower energy. The ambient temperature \(T\) controls the level of noise in the distribution in the same way as the variance \(\sigma^2\) in the Gaussian distribution (see Section 1.1). In particular, as \(T\) decreases, minimum-energy configurations become more probable. Let us illustrate this point using water. At high temperature, water molecules are randomly distributed in space as vapor. With decreasing temperature, they eventually condense into ice, which is their minimum-energy configuration. Formally, the equilibrium distribution of configurations \(c\) with energy \(H \left[ c \right]\) is proportional to \(\exp \left( -\frac{1}{T} H \left[ c \right] \right)\) at temperature \(T\). We say that they follow the (Boltzmann)-Gibbs distribution, a property that is assumed to hold analogously for energy-based models. In HNs, this choice leads to the distribution \[\label{eq:prelim95Hopfield95distribution} \mathop{\mathrm{P}}\left( \boldsymbol{\sigma} \mid \boldsymbol{\xi} \right) = \frac{\exp \left( -\frac{1}{T} H_{\text{HN}} \left[ \boldsymbol{\sigma} ; R \left( \boldsymbol{\xi} \right) \right] \right)}{\sum_{\boldsymbol{\sigma}^\prime \in \left\{ -1, 1 \right\}^N} \exp \left( -\frac{1}{T} H_{\text{HN}} \left[ \boldsymbol{\sigma}^\prime ; R \left( \boldsymbol{\xi} \right) \right] \right)}\tag{3}\] for the state \(\boldsymbol{\sigma}\) given the patterns \(\boldsymbol{\xi}\), where \(\left\{ -1, 1 \right\}^N\) is the set of all possible values of \(\boldsymbol{\sigma}\), which we recall contain only \(+1\)’s and \(-1\)’s. The denominator of Eq. (3 ) normalizes the distribution so that the total probability \(\sum_{\boldsymbol{\sigma} \in \left\{ -1, 1 \right\}^N} \mathop{\mathrm{P}}\left( \boldsymbol{\sigma} \mid \boldsymbol{\xi} \right)\) is equal to \(1\), and the temperature \(T\) is an abstract quantity representing the level of noise injected into \(\boldsymbol{\sigma}\). Following statistical mechanics conventions, we define the partition function \(Z \left( \boldsymbol{\xi} \right) = \sum_{\boldsymbol{\sigma}^\prime \in \left\{ -1, 1 \right\}^N} \exp \left( -\frac{1}{T} H_{\text{HN}} \left[ \boldsymbol{\sigma}^\prime ; R \left( \boldsymbol{\xi} \right) \right] \right)\) and the inverse temperature \(\beta = 1/T\). We then write the Gibbs distribution as \[\label{eq:Hopfield95distribution} \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma} \mid \boldsymbol{\xi} \right) = Z_\beta \left( \boldsymbol{\xi} \right)^{-1} \exp \left( -\beta H_{\text{HN}} \left[ \boldsymbol{\sigma} ; R \left( \boldsymbol{\xi} \right) \right] \right),\tag{4}\] where we use the subscript \(\beta\) to explicitly mark the dependence of \(\mathop{\mathrm{P}}\left( \boldsymbol{\sigma} \mid \boldsymbol{\xi} \right)\) and \(Z \left( \boldsymbol{\xi} \right)\) on \(\beta\). Now that we have derived this probability distribution, we explain how to learn the patterns. Note that the inverse temperature \(\beta\) and the number of patterns \(P\) are not learned. We call such parameters hyperparameters.

A relatively intuitive way to fit Eq. (4 ) to a dataset \(\mathcal{D} = \left\{ \boldsymbol{\sigma}^a \right\}_{a = 1}^M = \left( \boldsymbol{\sigma}^1, \boldsymbol{\sigma}^2, ..., \boldsymbol{\sigma}^M \right)\) is to find the patterns that make the dataset the most probable, a procedure known as maximum likelihood estimation. The probability that the Gibbs distribution (Eq. 4 ) produces a specific data point \(\boldsymbol{\sigma}^a\) from such a dataset is \(\mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^a \mid \boldsymbol{\xi} \right)\). Similarly, the probability that it generates the whole dataset is \(\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right) = \prod_{a = 1}^M \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^a \mid \boldsymbol{\xi} \right) = \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^1 \mid \boldsymbol{\xi} \right) \times \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^2 \mid \boldsymbol{\xi} \right) \times ... \times \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^M \mid \boldsymbol{\xi} \right)\). Maximum likelihood estimation amounts to finding the patterns \(\boldsymbol{\xi}\) that maximize \(\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right)\). It is usually formulated as maximizing the log likelihood \[\label{eq:log-likelihood} L \left( \boldsymbol{\xi} \right) = \frac{1}{M} \log \mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right) = \frac{1}{M} \sum_{a = 1}^M \log \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\sigma}^a \mid \boldsymbol{\xi} \right),\tag{5}\] which is mathematically equivalent and easier to handle numerically. Maximum likelihood estimation is popular for training energy-based models and NNs because there are various efficient ways to perform it numerically hinton2002training?, rosenblatt1958perceptron?, dempster1977maximum?, ruder2017overview?, which we will not explain in detail here. In particular, it was used to train HNs in decelle2021inverse?. For historical reasons, maximum likelihood training of NNs is generally phrased as minimizing a loss function equal to minus the log likelihood, and the shape of the loss as a function of the learnable parameters is called the loss landscape.

Alternatively, we can also train the patterns \(\boldsymbol{\xi}\) by sampling the posterior distribution \(\mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \mid \mathcal{D} \right)\), which represents the probability that an HN with patterns \(\boldsymbol{\xi}\) generates the data \(\mathcal{D}\). Although this algorithm is generally less efficient from a numerical point of view, it can be studied analytically to characterize the fundamental limits of NN training zdeborova2016statistical?, which has been applied to HNs alemanno2023hopfield?, RBMs hou2019minimal?, manzan2025effect?, linear models gardner1989unfinished?, gyorgyi1990first-order? and far beyond charbonneau2023spin?. The analytical calculations presented in this thesis follow this approach, which also ends up providing information about the loss landscape in Chapter 4 theriault2026saddle?. To understand how to calculate \(\mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \mid \mathcal{D} \right)\) from \(\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right)\), we first study a similar problem in a more familiar setting. Suppose that we know the probability \(\mathop{\mathrm{P}}\left( \text{the ground is wet} \mid \text{it rains} \right)\) that the ground is wet when it rains. We want to know the probability \(\mathop{\mathrm{P}}\left( \text{it rains} \mid \text{the ground is wet} \right)\) that it rains when the ground is wet to predict the weather. The probability that it rains and that the ground is wet is equal to the probability that it rains when the ground is wet multiplied by the probability that the ground is wet. Formally, this relationship is written as \[\begin{align} \mathop{\mathrm{P}}\left( \text{it rains}, \text{the ground is wet} \right) &= \mathop{\mathrm{P}}\left( \text{it rains} \mid \text{the ground is wet} \right) \mathop{\mathrm{P}}\left( \text{the ground is wet} \right). \end{align}\] Crucially, it also works in the other direction, that is \[\begin{align} \mathop{\mathrm{P}}\left( \text{it rains}, \text{the ground is wet} \right) &= \mathop{\mathrm{P}}\left( \text{the ground is wet} \mid \text{it rains} \right) \mathop{\mathrm{P}}\left( \text{it rains} \right). \end{align}\] Combining both of these equalities, we find \[\begin{align} \mathop{\mathrm{P}}\left( \text{it rains} \mid \text{the ground is wet} \right) &= \frac{\mathop{\mathrm{P}}\left( \text{the ground is wet} \mid \text{it rains} \right) \mathop{\mathrm{P}}\left( \text{it rains} \right)}{\mathop{\mathrm{P}}\left( \text{the ground is wet} \right)}, \end{align}\] which is known as Bayes’ theorem. In the case of HNs, it gives \[\label{eq:Hopfield95posterior} \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \mid \mathcal{D} \right) = \frac{\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right) \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \right)}{\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \right)},\tag{6}\] where the prior \(\mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \right)\) represents the knowledge that we have about the distribution of the patterns before training and \(\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \right) = \sum_{\boldsymbol{\xi} \in \left\{ -1, 1 \right\}^{N \times P}} \mathop{\mathrm{P}}_\beta \left( \mathcal{D} \mid \boldsymbol{\xi} \right) \mathop{\mathrm{P}}_\beta \left( \boldsymbol{\xi} \right)\) normalizes the distribution in the same way as \(Z_\beta \left( \boldsymbol{\xi} \right)\) in the Gibbs distribution (Eq. 4 ). By analogy, we call the normalization constant \(\mathop{\mathrm{P}}_\beta \left( \mathcal{D} \right)\) the posterior partition function and write it as \(\mathcal{Z}_\beta \left( \mathcal{D} \right)\).

1.3 Restricted Boltzmann machines and dense Hopfield networks↩︎

As for HNs, RBMs and dense HNs are energy-based models whose corresponding Hamiltonians and Gibbs distributions can be fit on data using maximum likelihood estimation or posterior sampling. Here, we summarize the differences between them.

The generalized HNs that are nowadays called dense were developed shortly after the original to improve upon its storage capacity chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, horn1988capacities?. These models have been receiving renewed attention under the name dense associative memory (DAM) since they were used for pattern recognition in krotov2016dense?. In the process, the meaning of the term dense associative was also broadened to encompass more general models with large storage capacities krotov2021large?, krotov2025modern?, which were in turn related to various other ML paradigms ramsauer2020hopfield?, hoover2023memory?, ota2023attention?, ambrogioni2024search?, karakida2024hierarchical?. This evolution in terminology is the reason why we adopt the convention of calling the generalized HNs first considered in chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, horn1988capacities? dense HNs rather than DAMs. The dense HN Hamiltonian is \[\begin{gather} \label{eq:DHN95Hamiltonian} H_{\text{DHN}} \left[ \boldsymbol{\sigma} ; \mathbf{J} \right] = -\sum_{i_1 < ... < i_p = 1}^N J_{i_1 ... i_p} \sigma_{i_1} ... \sigma_{i_p}, \addtocounter{equation}{1}\\ \text{where} \quad J_{i_1 ... i_p} = R_{i_1 ... i_p} \left( \boldsymbol{\xi} \right) = \frac{p!}{N^{p-1}} \sum_{\mu = 1}^M \xi_{i_1}^\mu ... \xi_{i_p}^\mu. \end{gather}\tag{7}\] The weights of Eq. (7 ) connect the units \(\boldsymbol{\sigma} = \left\{ \sigma_i \right\}_{i = 1}^N = \left( \sigma_1, \sigma_2, ..., \sigma_N \right)\) within groups of size \(p\). These connections, which are called \(p\)-body interactions, reduce to those of the original HN when \(p = 2\).¹

Dense HNs bridge the gap between ML and MM in the sense that they learn prototypes of the data krotov2016dense? when trained with an ML algorithm. For example, their patterns become prototypes of digits when trained on the MNIST dataset of handwritten digits (see Fig. 7). These dense HNs were observed to have a variable level of adversarial robustness as a function of \(p\) krotov2018dense?, offering an interesting opportunity to study adversarial attacks and how to mitigate them.

Although they have \(2\)-body interactions like the original HN, RBMs differ from it in that they represent data on two levels of abstraction. Intuitively, these two levels are roughly analogous to how we see objects and how we describe them mentally (for example, thinking of a fire as “hot”, “red” and “bright”). In practice, however, the abstract concepts that RBMs associate with data are different from our own mental representations. Formally, the RBM Hamiltonian is \[\label{eq:RBM95Hamiltonian} H_{\text{RBM}} \left[ \mathbf{v}, \mathbf{h} ; \mathbf{J} \right] = -\sum_{i = 1}^N \sum_{\mu = 1}^P J^{\mu}_i v_i h_\mu,\tag{8}\] where the visible layer \(\mathbf{v} = \left\{ v_i \right\}_{i = 1}^N = \left( v_1, v_2, ..., v_N \right)\) is the data expressed in terms of its components, the hidden layer \(\mathbf{h} = \left\{ h_\mu \right\}_{\mu = 1}^P\) is a high-level representation of the data in terms of abstract concepts, and the learnable weights \(\mathbf{J} = \left\{ J^\mu_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\) are connections that the RBM makes between \(\mathbf{v}\) and \(\mathbf{h}\). As in HNs, these connections are bidirectional, which means that they represent how \(\mathbf{v}\) and \(\mathbf{h}\) are correlated rather than how one influences the other (see Fig. 8). The units \(v_i\) and \(h_\mu\) of the two layers can be bound to \(\pm1\) as those of HNs, but various other choices are also possible. For example, they can be integers within an interval Salakhutdinov2007restricted? or arbitrary real numbers, thus defining multiple types of RBMs barra2012equivalence?. Every such RBM has a Gibbs distribution of the form \[\label{eq:RBM95distribution} \mathop{\mathrm{P}}_\beta \left( \mathbf{v}, \mathbf{h} \mid \mathbf{J} \right) = Z_\beta \left( \mathbf{J} \right)^{-1} \mathop{\mathrm{P}}_0 \left( \mathbf{v} \right) \mathop{\mathrm{P}}_0 \left( \mathbf{h} \right) \exp \left( -\beta H \left[ \mathbf{v}, \mathbf{h} ; \mathbf{J} \right] \right),\tag{9}\] where the inverse temperature \(\beta\) regulates the strength of the interactions that \(\mathbf{v}\) and \(\mathbf{h}\) exchange through \(H\), the priors \(\mathop{\mathrm{P}}_0 \left( \mathbf{v} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} \right)\) are the distributions of \(\mathbf{v}\) and \(\mathbf{h}\) when they are decoupled, i.e. \(\beta = 0\), and the partition function \(Z_\beta \left( \mathbf{J} \right)\) is again defined so that the total probability is equal to \(1\). The priors on both layers, which represent the knowledge that we have on the data before training, usually have a simple form such that their units are (a priori) independent and identically distributed (i.i.d.) with respect to each other. For example, neurons restricted to \(\pm1\) usually have a Rademacher prior \[\begin{align} \mathop{\mathrm{P}}_0 \left( v_i \right) = \begin{cases} 1/2 &\quad \text{if v_i = \pm1} \\ 0 &\quad \text{otherwise}, \end{cases} \end{align}\] and units that can be arbitrary real numbers usually have a Gaussian prior \(\mathop{\mathrm{P}}_0 \left( v_i \right) = \frac{\exp \left( -\frac{1}{2 \sigma^2} v_i^2 \right)}{\sqrt{2\pi \sigma^2}}\). Simple priors with non-i.i.d. units are also possible. For example, RBM layers whose units are restricted to always be at a distance of \(\sqrt{\sum_{i = 1}^N v_i^2} = 1\) from the origin can be described using the hyperspherical prior \[\begin{align} \mathop{\mathrm{P}}_0 \left( \mathbf{v} \right) = \begin{cases} \frac{1}{\Omega_N \left( \beta \right)} &\quad \text{if \sqrt{\sum_{i = 1}^N v_i^2} = 1} \\ 0 &\quad \text{otherwise}, \end{cases} \end{align}\] where \(\Omega_N \left( 0 \right)\) is the size of the region of space where \(\sqrt{\sum_{i = 1}^N v_i^2} = 1\), which is called the unit hypersphere. Additionally, the categorical (or Potts potts1952some?, wu1982potts?) prior \[\mathop{\mathrm{P}}_0 \left( \mathbf{v} \right) = \begin{cases} \frac{1}{N + 1} &\quad \text{if \mathbf{v} \in \left\{ 0, 1 \right\}^N and \sum_{i = 1}^N v_i \in \left\{ 0, 1 \right\}} \\ 0 &\quad \text{otherwise}, \end{cases}\] is used in RBM layers where at most one unit \(v_i\) can be activated at once (equal to \(1\)), while all the others are turned off (equal to \(0\)).

Given that the priors have a simple form, using a higher temperature \(T = 1/\beta\) encourages the RBM to learn simpler weights \(J^\mu_i\) during training, which is known as regularization in ML. Intuitively, regularization encourages NNs not to overthink problems with relatively simple solutions, which can help them perform better.

RBMs are trained by fitting to the data their marginal distribution \(\mathop{\mathrm{P}}_\beta \left( \mathbf{v} \mid \mathbf{J} \right) = \int_{\mathbf{h}} d \mathbf{h} \mathop{\mathrm{P}}_\beta \left( \mathbf{v}, \mathbf{h} \mid \mathbf{J} \right)\), which embodies how they represent the data with all their hidden units simultaneously. The marginal distribution often does not have a simple closed form and must be approximated, such as by repeatedly sampling the conditional distributions \(\mathop{\mathrm{P}}_\beta \left( \mathbf{h} \mid \mathbf{v}, \mathbf{J} \right)\) and \(\mathop{\mathrm{P}}_\beta \left( \mathbf{v} \mid \mathbf{h}, \mathbf{J} \right)\) following the contrastive divergence algorithm introduced in hinton2002training?. In some cases, the marginal distributions of RBMs with specific priors reduce to the Gibbs distributions of other energy models Decelle2021restricted?. For example, the marginal distribution of an RBM with a Rademacher prior on its visible units and a Gaussian prior on its hidden units simplifies to \[\mathop{\mathrm{P}}_\beta \left( \mathbf{v} \mid \mathbf{J} \right) = Z_\beta \left( \mathbf{J} \right)^{-1} \exp \left( \beta \sum_{i < j = 1}^N \sum_{\mu = 1}^P J^\mu_i J^\mu_j v_i v_j \right),\] which is the Gibbs distribution of a Hopfield network with inverse temperature \(\beta N\), memorized patterns \(\mathbf{J}^\mu = \left\{ J^\mu_i \right\}_{i = 1}^N\) and units \(\mathbf{v}\). In other words, RBMs with these priors are formally equivalent to HNs barra2012equivalence?. By analogy, the rows \(\mathbf{J}^\mu\) of the weights are occasionally called patterns and written as \(\boldsymbol{\xi}^\mu\) even for other choices of priors.

RBMs with a hyperspherical prior on their visible units and a categorical prior on their hidden units become Dense Associative Memory (DAM) models that learn prototypes of the data like dense HNs. Formally, the marginal distribution of such an RBM simplifies to \[\label{eq:DAM95distribution} \mathop{\mathrm{P}}_\beta \left( \mathbf{v} \mid \mathbf{J} \right) = \frac{1}{P + 1} \sum_{\mu = 1}^P \frac{\exp \left( \beta \sum_{i = 1}^N J^\mu_i v_i \right)}{\Omega_N \left( \beta \right)} + \frac{1}{P + 1} \frac{1}{\Omega_N \left( 0 \right)},\tag{10}\] where the normalization constant \(\Omega_N \left( \beta \right)\) is defined in Chapter 4. In that Chapter, we discuss such DAMs and their affinity for prototype learning in more detail.

1.4 Statistical mechanics of networks in the teacher-student setting↩︎

We now present the main theoretical tool that we use to study the statistical mechanics of RBMs and dense HNs: the teacher-student setting gardner1989unfinished?. This tool has been used numerous times to solve various fundamental problems in the field of ML and NNs charbonneau2023spin?.

In the teacher-student setting, a teacher NN with weights \(\mathbf{J}^*\) generates a sample dataset \(\mathcal{D}\) on which a student NN then trains its own weights \(\mathbf{J}\). In other words, the teacher shows the dataset \(\mathcal{D}\) to the student to teach it how to mimic its contents. This framework provides a controlled environment for studying NN behavior with analytical calculations. Although NNs with suitable weights can represent arbitrarily complex data HORNIK1989feedforward?, the teacher weights \(\mathbf{J}^*\) are usually determined according to a simple criterion to simplify calculations. For example, they can be fixed samples from a simple probability distribution. Moreover, in the case of dense HNs, we assume that the weights \(\mathbf{J}^*\) and \(\mathbf{J}\) decompose into patterns \(\boldsymbol{\xi}^{* \mu}\) and \(\boldsymbol{\xi}^{\nu}\) following the second line of Eq. (7 ).

In our work, the student has limited knowledge of the architecture of the teacher. To be more precise, it knows whether the teacher is a dense HN or an RBM and, in the latter case, the priors on the RBM’s neurons. However, unless explicitly stated otherwise, it does not know the hyperparameters of the teacher, such as its inverse temperature \(\beta^*\), its number of patterns \(P^*\), its interaction order \(p^*\), etc. Therefore, it adopts the same form as the teacher, but with hyperpameters that are generally different, which we write without an asterisk (*) to distinguish them from those of the teacher. Following statistical mechanics conventions, we call the special case where the student has the same hyperparameters as the teacher the Nishimori line nishimori1980exact?, nishimori2001statistical?, contucci2009spin?, iba1999nishimori?.

Given that the teacher and the student are two instances of the same type of model, the degree of similarity between their weights is a good indication of the performance of the student in learning the data. We measure this degree of similarity using the overlaps \(Q \left( \boldsymbol{\xi}^{* \mu}, \boldsymbol{\xi}^{\nu} \right) = \frac{1}{N} \sum_{i = 1}^N \xi^{* \mu}_i \xi^{\nu}_i\) between the patterns of the student and those of the teacher, which we recall are the rows \(\mathbf{J}^\nu\) and \(\mathbf{J}^{* \mu}\) of the weights for RBMs. In the case where all components \(\xi^{* \mu}_i\) and \(\xi^{\nu}_i\) of the patterns \(\boldsymbol{\xi}^{* \mu}\) and \(\boldsymbol{\xi}^{\nu}\) are equal to \(\pm1\), \(Q \left( \boldsymbol{\xi}^{* \mu}, \boldsymbol{\xi}^{\nu} \right) = 1\) means that the student learns the teacher pattern \(\boldsymbol{\xi}^{* \mu}\) perfectly. We use other types of overlaps to probe various other aspects of learning. For example, the overlaps \(Q \left( \boldsymbol{\xi}^{1 \mu}, \boldsymbol{\xi}^{2 \nu} \right) = \frac{1}{N} \sum_{i = 1}^N \xi^{1 \mu}_i \xi^{2 \nu}_i\) between two learning attempts \(\boldsymbol{\xi}^{1} = \left\{ \xi^{1 \mu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\) and \(\boldsymbol{\xi}^{2} = \left\{ \xi^{2 \mu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\) of the student represent its tendency to stay frozen in specific pattern configurations rather than visiting all possible values of \(\boldsymbol{\xi}\). \(Q \left( \boldsymbol{\xi}^{1 \mu}, \boldsymbol{\xi}^{2 \nu} \right)\) can sometimes be relatively large even when \(Q \left( \boldsymbol{\xi}^{* \mu}, \boldsymbol{\xi}^{\nu} \right)\) is close to \(0\), which represents the situation where the student mistakenly learns a very different pattern from that of the teacher, or in other words misunderstands what the teacher is trying to teach it. In general, the overlaps change with the hyperparameters of the teacher-student setting and the number of samples \(M\) in the training dataset \(\mathcal{D}\), defining different regimes of data representation.

Our theoretical analyses are based on computing mean overlaps, which are also known as order parameters, using the replica method (explained at length in nishimori2001statistical?). This method is well-established in the statistical mechanics community, where it is used to study a type of physical system called spin glass charbonneau2023spin?, nishimori2001statistical?. Spin glasses are characterized by disordered components that are perpetually frozen, or quenched, as the system converges to equilibrium. In the case of NNs in the teacher-student setting, the quenched disorder is the patterns of the teacher and the dataset that it generates, which are kept fixed during learning. Conversely, the student patterns, which are allowed to evolve throughout training, are not quenched. Although they can sometimes end up freezing in disordered spin-glass states characterized by \(Q \left( \boldsymbol{\xi}^{* \mu}, \boldsymbol{\xi}^{\nu} \right) \approx 0\) and \(Q \left( \boldsymbol{\xi}^{1 \mu}, \boldsymbol{\xi}^{2 \nu} \right) > 0\), it is a consequence of the learning task rather than a hard constraint.

We use these calculations to characterize the learning phases of RBMs and dense HNs as a function of the number of samples \(M\) in the training dataset \(\mathcal{D}\) and the hyperparameters of the teacher-student setting. In particular, we calculate phase diagrams representing the different phases and phase transitions of RBMs and dense HNs as a function of the hyperparameters and the size of the training dataset. For example, Fig. (9), which is taken from alemanno2023hopfield?, is the phase diagram of HNs in the teacher-student setting when the student and the teacher have a single pattern and the same inverse temperature \(\beta = \beta^*\). We see three phases and two phase transitions (black lines) as a function of the amount of data \(\gamma = M/N\) and the temperature \(T = 1/\beta\), which we will now explain in terms of the mean teacher-student overlap (subsequently \(m\)) and the mean overlap between the teacher pattern and the samples that it generates (subsequently \(r\)). Below the first phase transition, at the bottom of the plot, the temperature \(T = 1/\beta\) is low enough for the teacher to generate highly informative training samples (\(r > 0\)) and for the student HN to perfectly learn the teacher pattern through them (\(m = 1\)). Between the two phase transitions, the training dataset is much noisier (\(r = 0\)), but the student can still learn the teacher pattern imperfectly (\(0 < m < 1\)). Finally, above the second phase transition, the temperature \(T = 1/\beta\) is too high and the dataset too small for the student HN to learn anything about the teacher pattern (\(m = 0\)).

1.5 Outline of the thesis↩︎

In Chapter 2 (based on a work with Daniele Tantari, published as theriault2024dense?), we study dense HNs with a single pattern in the teacher-student setting. In particular, we investigate the relationship between a dense HN with a single pattern trained using \(M\) data points (the inverse model) and a dense HN with \(M\) random memorized patterns (the direct model), the latter of which is well understood thanks to previous studies chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, horn1988capacities?, krotov2016dense?, albanese2022replica?. To be more precise, we show analytically that the phase where the inverse model with \(\beta = \beta^*\) and \(p = p^*\) is capable of learning coincides with the phase where the direct model freezes in random states (called the spin-glass phase), refining the well-established quiet planting bound that the learning phase of the inverse model extends at least as far as the spin-glass phase of the direct model zdeborova2016statistical?, achlioptas2008algorithmic?, krzakala2009hiding?, zdeborova2011quiet?. Pushing our analysis of this relationship further, we study noise tolerance and adversarial attacks in dense HNs memorizing as patterns the \(M\) examples generated by the teacher. In particular, we derive a mathematical formula quantifying the adversarial robustness of dense HNs at zero temperature, and we clarify why the adversarial robustness of the dense HNs studied in krotov2018dense? changes as a function of the interaction order \(p\). Our results lay the groundwork for a broader statistical mechanics theory of adversarial attacks, contemporaneously with complementary studies of adversarial attacks in linear models tanner2024high?, vilucchio2024geometry?, vilucchio2025existence?.

In Chapter 3 (based on a work with Francesco Tosello and Daniele Tantari, published as theriault2025modeling?), we study RBMs in the teacher-student setting where the numbers of hidden units \(P^*\) and \(P\) are free to be any natural numbers much smaller than the number of samples in the training dataset. This work generalizes the results of hou2019minimal?, which studied the teacher-student setting with \(P^* = P = 2\) hidden units. First, we investigate the relationship between our framework and RBMs in the teacher-student setting with \(P^* = P = 1\) hidden units. To be more precise, we show that an RBM with \(P\) hidden units learns data generated by a teacher with \(P^* = P\) random patterns with the same overlaps as \(P\) RBMs with one hidden unit each, thus validating a conjecture formulated in barra2017phase?. Going one step further, we study the scenario where the student has more hidden units than necessary to learn the training data, i.e. \(P > P^*\). We show that, in this context, only \(P^*\) of its patterns end up learning the data, while the others freeze in disordered states. Empirically, we investigate the connection between this low-dimensional learning strategy and the lottery ticket hypothesis according to which a generic randomly initialized overparameterized NN contains subnetworks that fit data with similar accuracy as the entire trained network when they are extracted from it and trained independently frankle2018lottery?. Finally, we show that, when the teacher patterns are correlated together rather than purely random, the student undergoes successive phase transitions where it learns the data in increasingly detailed ways as the size of the training dataset increases.

In Chapter 4 (based on a work with Daniele Tantari, published as theriault2026saddle?), we study the teacher-student setting of an RBM that fits into the class of Dense Associative Memory (DAM) models. This DAM generalizes Eq. (10 ), which is itself a close relative of the DAM studied in ramsauer2020hopfield?, ota2023attention?, lucibello2024exponential?, to a form capable of both supervised and unsupervised pattern classification. As in Chapter 3, we study a model with a generic number of hidden units. Based on our results, we propose a novel regularization scheme that makes training significantly more stable. Inspired by the observation that dense HNs cross multiple saddle points of the loss landscape during training boukacem2024waddington?, we study the saddle points of our DAM. We demonstrate that the weights learned by DAMs with a relatively small number of hidden units are saddle points of the loss landscape of larger DAMs (see Eq. 5 ), thus extending studies of an analogous property in related models Decelle2021restricted?, kappen1993using?, KAPPEN1995deterministic?, nijman1997symmetry?, rose1990statistical?, kloppenburg1997deterministic?, akaho2000nonmonotonic?. We exploit this hierarchy by letting our DAM grow during learning wu2019splitting?, wang2019energy?, which makes it significantly faster to train as a function of its final size.

2 Dense Hopfield networks in the teacher-student setting↩︎

2.1 Introduction↩︎

Hopfield networks are artificial neural networks that model associative memory hopfield1982neural?. In the Hopfield model, examples \(\sigma \in \left\{ -1, 1 \right\}^N\) of memories \(\xi^{\mu} \in \left\{ -1, 1 \right\}^N\), \(\mu=1,\ldots,M\), are retrieved by sampling the Gibbs distribution of a \(2\)-body Hamiltonian \(H \left[ \sigma \mid \xi \right]\) at a given temperature \(T\) amit1985spin?. Hopfield networks can be trained in a biologically plausible way using Hebb’s rule hopfield1982neural?, hebb2005organization?, which leads to \(H \left[ \sigma \mid \xi \right] = -\frac{1}{N} \sum_{\mu = 1}^M \left( \sum_{i = 1}^{N} \xi^{\mu}_i \sigma_i \right)^2\). However, they can only store up to \(M \sim \mathcal{O} \left( N \right)\) i.i.d. random memories in the limit of large \(N\) hopfield1982neural?, cover1965geometrical?, amit1985storing?. One way to find this scaling is to study the phase diagram of \(H \left[ \sigma \mid \xi \right]\) as a function of the temperature \(T\) and load \(\alpha = \frac{M}{N}\) amit1985storing?, where the so-called ferromagnetic phase, which extends up to \(\alpha \approx 0.14\), corresponds to accurate retrieval.

Since Hopfield’s seminal work, several generalizations have been investigated in relation to their critical storage capacity and retrieval capabilities. For example, parallel retrieval has been studied in relation to pattern sparsity agliari2012multitasking?, agliari2013immunemedium?, agliari2013immune?, sollich2014extensive?, agliari2017retrieving? or hierarchical interactions agliari2015retrieval?, agliari2018non?, agliari2015hierarchical?, agliari2014metastable?, agliari2015topological?, and non-universality has been shown with respect to more general pattern entries and unit priors barra2017phase?, barra2018phase?, barra2015multi?, agliari2017neural?, barra2012glassy?, genovese2020legendre?, rocchi2017high?. Efforts to overcome the \(\mathcal{O} \left( N \right)\) limitation of the capacity led to the development of a novel class of modern Hopfield networks ramsauer2020hopfield?, widrich2020modern?, krotov2021large?, which are sometimes called dense due to their faculty to store much more memories than the original Hopfield model krotov2016dense?. These neural networks surpass \(\mathcal{O} \left( N \right)\) storage capacity by using higher-order interactions instead of the original \(2\)-body couplings chen1986high?, psaltis1986nonlinear?, baldi1987number?, gardner1987multiconnected?, abbott1987storage?, horn1988capacities?. In particular, Gardner gardner1987multiconnected? calculated the Replica-Symmetric (RS) phase diagram of the Hamiltonian \(H \left[ \sigma \mid \xi \right] = -\sum_{i_1 < ... < i_p = 1}^{N} J_{i_1 ... i_p} \sigma_{i_1} ... \sigma_{i_p}\) with \(p\)-body interactions \(J_{i_1 ... i_p} = \frac{p!}{N^{p - 1}} \sum_{\mu=1}^{M} \xi^{\mu}_{i_1} ... \xi^{\mu}_{i_p}\) conditioned on i.i.d. random memories \(\xi^{\mu} \in \left\{ -1, 1 \right\}^N\), finding a \(M=\mathcal{O} \left( N^{p-1} \right)\) storage capacity. These calculations were later extended to include the effects of \(1\)-step Replica Symmetry Breaking (1RSB) albanese2022replica?.

Although they draw a rather detailed picture of the retrieval of individual i.i.d. random memories, these results are not the end of the story. First of all, 1RSB calculations allegedly struggle to find the paramagnetic to spin-glass phase transition accurately at large \(p\) because of numerical instability issues albanese2022replica?. Second of all, dense Hopfield networks have been rapidly gaining a renewed attention for reasons other than their storage capacity since a recent paper krotov2016dense? by Krotov and Hopfield (K & H), where they were used as a trainable machine learning architecture. For instance, they have been related to transformers ramsauer2020hopfield?, hoover2023energy? and diffusion models hoover2023memory?, ambrogioni2024search?, and they were found to be significantly more explainable and adversarially robust than feedforward neural networks with ReLU activation functions krotov2016dense?, krotov2018dense?.

One such aspect of dense Hopfield networks that is still poorly understood is their performance as generative models for unsupervised learning, where they are trained over some given dataset to reproduce its probability distribution. As far as we are aware, this problem has not yet been studied theoretically for \(p\)-body models with \(p \geq 3\). However, it was studied for the original \(2\)-body Hopfield network by using the teacher-student setting alemanno2023hopfield? first described in barra2017phase?, barra2018phase?, decelle2021inverse?. In the teacher-student setting, which is also called inverse problem in opposition to the direct problem of random pattern retrieval, a student model \(H \left[ \xi \mid \sigma \right]\) is trained with \(M\) teacher examples \(\sigma^a \sim H \left[ \sigma^a \mid \xi^* \right]\) conditioned on the planted pattern \(\xi^*\). In other words, the student tries to infer the pattern \(\xi^*\) of the teacher using a structured set of examples \(\sigma^a\).

At finite load \(\alpha = \frac{M}{N}\), two regimes of pattern retrieval were found: example Retrieval (\(eR\)) and signal Retrieval (\(sR\)). In the \(eR\) phase, the student tries to reconstruct \(\xi^*\) by directly retrieving the examples \(\sigma^a\), which is a good strategy provided that they are strongly correlated with \(\xi^*\). In the \(sR\) phase, on the other hand, retrieval is done by extracting subtle cues from weakly correlated examples. The two types of examples used in these two retrieval strategies are respectively called prototypes and features of \(\xi^*\) krotov2016dense?. Interestingly, a prototype regime and a feature regime were also observed by K & H in dense Hopfield networks trained to classify real data krotov2016dense?, where it was found that the prototype regime is significantly more adversarially robust than the feature regime. In other words, the prototype regime is more resistant than the feature regime to small data perturbations that are specifically designed to cause incorrect classification biggio2013evasion?, szegedy2013intriguing?. This prototype approach is arguably a big step towards designing adversarially robust neural networks, a long-standing problem that still lacks a fully satisfying solution goodfellow2015explaining?, madry2018towards?, muhammad2022survey?.

In this work, we study the performance of \(p\)-body Hopfield networks in the teacher-student setting, revealing a prototype regime and a feature regime as in the \(2\)-body model. In Section 2.2, we review Gardner’s main results in studying \(p\)-body Hopfield models and summarize what the rest of the literature on spin-glass models with \(p\)-body interactions tell us about the paramagnetic to spin-glass phase transition in \(p\)-body Hopfield models. In Section 2.3, we compute the phase diagram of these \(p\)-body models in the teacher-student setting. In Section 2.4.1, we discuss the transition to the retrieval phase in the inverse problem. In Section 2.4.2, we compare this retrieval transition against the transition to the spin-glass phase in the direct problem. Despite their different nature, we show that these two transitions are equivalent on the Nishimori line where the teacher and the student have the same \(p\) and \(T\) nishimori1980exact?, nishimori2001statistical?, contucci2009spin?, iba1999nishimori?. In Section 2.4.3, we discuss the phase diagram on the Nishimori line in more details. In Section 2.4.4 and Section 2.4.5, we discuss the phase diagram outside of the Nishimori line. First of all, we investigate the effect of using an inference temperature different from the dataset noise. Second of all, we reveal that using a larger \(p\) for the student than the teacher gives the student an extensive tolerance against both teacher noise and pattern interference. Finally, in Section 2.4.6, we derive a closed-form expression that measures the adversarial robustness of the student at zero temperature and explain what our results reveal about the nature of adversarial attacks.

2.2 Overview of Gardner’s results↩︎

Consider the \(p\)-body Hamiltonian \[\label{eq:H} H \left[ \sigma \mid \xi \right] = -\sum_{i_1 < ... < i_p = 1}^{N} J_{i_1 ... i_p} \sigma_{i_1} ... \sigma_{i_p}= -\frac{p!}{N^{p - 1}}\sum_{i_1 < ... < i_p = 1}^{N} \sum_{\mu=1}^{M} \xi^{\mu}_{i_1} ... \xi^{\mu}_{i_p} \sigma_{i_1} ... \sigma_{i_p} \,,\tag{11}\] conditioned on a set of \(M=\frac{\alpha N^{p-1}}{p!}\) quenched memories \(\xi^\mu\in\{-1,1\}^N\), \(\mu=1,\ldots,M\), sampled i.i.d. from the Rademacher distribution \(\frac{1}{2} \left[ \delta \left( \xi^\mu_i-1 \right) + \delta \left(\xi^\mu_i +1 \right) \right]\). In the direct model, patterns \(\sigma\) are in turn sampled from the equilibrium Gibbs distribution \(P(\sigma \mid \xi)=Z^{-1}e^{-\beta H[\sigma \mid \xi]}\), where \(\beta\geq 0\) is the inverse temperature and \(Z=\sum_\sigma e^{-\beta H[\sigma \mid \xi]}\) is the system’s partition function. The so-called direct problem studied by Gardner gardner1987multiconnected? consists of quantifying the performance of this model as a method of memory retrieval. In that context, the overlap \(\frac{1}{N} \sum_i \xi^{\mu}_i \sigma_i\) is a good measure of retrieval accuracy, and its expected value can be derived from the quenched free entropy \(f = \frac{1}{N} \left\langle \log Z \right\rangle_{\xi}\) in the thermodynamic limit \(N\to\infty\). At finite \(p\), Gardner used the (non-rigorous) replica trick charbonneau2022replica? to evaluate the RS approximation of \(f\) (see also Appendix 2.7) in terms of a variational principle of the form \[\begin{align} f = \lim_{N\to\infty}\frac{1}{N} \left\langle \log Z \right\rangle_\xi & = \lim_{N\to\infty, L \rightarrow 0} \left( \frac{\partial}{\partial L} \left[ \frac{1}{N} \log \left\langle Z^L \right\rangle_\xi \right] \right)=\mathop{\mathrm{Extr}}_{m,k,q,k,r} \;f(m,k,q,r)\,, \end{align}\] whose solution is \[\begin{align} \label{eq:saddlepointdirect} q & = \int_{\mathbb{R}} dx \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) \tanh^2 \left( \beta \left[ \sqrt{\alpha r} x + k \right] \right) \nonumber \,, \\ m & = \int_{\mathbb{R}} dx \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) \tanh \left( \beta \left[ \sqrt{\alpha r} x + k \right] \right) \nonumber \,, \\ r & = p q^{p-1} \,,\nonumber \\ k & = p m^{p-1}\,, \end{align}\tag{12}\] and the order parameters \(m\) and \(q\) are to be interpreted as expected overlaps. To be more precise, \(m\) can be shown to be the expected overlap of a retrieval attempt \(\sigma\) against one memory in the thermodynamic limit, i.e. \(m=\lim_{N\to\infty} \left\langle \frac{1}{N} \sum_{i} \xi^\mu_i \sigma_i \right\rangle_{\xi,\sigma}\). Similarly, \(q\) is the expected overlap between two retrieval attempts \(\sigma^1\) and \(\sigma^2\), i.e. \(q =\lim_{N\to\infty} \left\langle \frac{1}{N} \sum_i \sigma^{1}_i \sigma^{2}_i \right\rangle_{\xi,\sigma}\) or equivalently \(q=\lim_{N\to\infty}\left\langle\frac{1}{N} \sum_i \left\langle\sigma_i\right\rangle^2_\sigma \right\rangle_\xi\). Intuitively, \(q\) measures the tendency of the system to stay frozen in specific configurations rather than visiting all possible values of \(\sigma\).

The resulting RS phase diagram (see Fig. 10) are derived from the value of the order parameters as a function of three hyperparameters: the interaction order \(p\), temperature \(T=1/\beta\) and load \(\alpha = \frac{M p!}{N^{p-1}}\). There are four different phases:

• In the Paramagnetic phase (\(P\)), the overlaps \(m\) and \(q\) both vanish. The network does not retrieve any specific pattern: sampled configurations are completely random.

• In the Spin-Glass phase (\(SG\)), \(m\) vanishes but \(q > 0\). In other terms, the network does not retrieve individual stored memories but rather converges to spurious patterns depending on all the memories in a nontrivial way.

• In the signal Retrieval phases (\(lR\) and \(gR\)), \(m \neq 0\) and \(q > 0\), which means that the network is able to retrieve the stored memories. \(lR\) and \(gR\) are respectively locally stable and globally stable. In other words, local retrieval \(lR\) is only attainable from initial conditions in a limited neighborhood of a memory \(\xi^\mu\), while global retrieval \(gR\) is accessible from any initial conditions given enough time. These two phases are said to be ferromagnetic.

Gardner also calculated the exact \(p \rightarrow \infty\) phase diagram without making any assumptions about replica symmetry gardner1987multiconnected?. In this limit, the resulting Paramagnetic to Spin-Glass (\(P\)-\(SG\)) phase transition occurs at a temperature \(T_E(\alpha)\) that coincides with the boundary of the region where the total entropy of the paramagnetic phase becomes negative, given by \(\beta^2 \alpha = 2 \log 2\) (white dashed line in Fig. 10).

At finite \(p\), Gardner’s results only tell us that the model cannot be in the paramagnetic phase below \(T_E(\alpha)\). Therefore, a spin-glass transition should occur at a temperature \(T_s\left(\alpha,p\right) \geq T_E(\alpha)\). Since the RS spin glass solution of Eqs. (12 ) exists only below \(T_E(\alpha)\) (violet region in Fig. 10), the spin-glass transition must be towards a RSB spin-glass phase.

Outside of the signal retrieval phases, the free entropy of the direct model is the same as for the spin-glass model with \(p\)-body Gaussian interactions where the temperature is rescaled by a factor of \(\sqrt{2 \alpha}\) gross1984simplest?, gardner1985spin?. Therefore, the spin-glass and paramagnetic solutions are the same in the direct model as in this Gaussian spin-glass model, and we expect the exact phase diagrams of both models to be identical when the direct model is not in its signal retrieval phases. According to previous work on the Gaussian model with finite \(p\) gardner1985spin?, a 1RSB solution with \(m = k = 0\) exists and is globally stable throughout a whole phase below \(T_s(\alpha,p)\geq T_E(\alpha)\) (see Fig. 10). This solution becomes unstable at a lower transition temperature \(T_G(\alpha, p)\) (see Fig. 11), below which multiple steps of RSB are required. In the limit of \(p\to\infty\), it holds that \(T_s(\alpha,p)\to T_E(\alpha)\) and \(T_G(\alpha,p)\to 0\). In other terms, the direct model becomes 1RSB, which is consistent with the fact that it is converging to a random energy model with temperature rescaled by \(\sqrt{2 \alpha}\) derrida1980random?, gross1984simplest?, gardner1987multiconnected?. Finally, we mention that this type of models exhibits a random first order transition phenomenology monasson1995structural?, montanari2003nature?, crisanti2005complexity?, franz2017universality?: there is in fact a range of temperatures \(T_s(\alpha,p)\leq T\leq T_d(\alpha,p)\) where the dynamics get trapped in an exponential number of metastable clusters, with an emerging RSB structure that does not affect the free energy (see Fig. 11). This range of temperatures thus defines a so-called dynamical 1RSB (d1RSB) phase. Below \(T_s(\alpha,p)\), the number of clusters is no longer exponential, and the system undergoes the thermodynamic 1RSB phase transition that we mentioned previously. The critical temperatures \(T_G(\alpha,p)\), \(T_s(\alpha,p)\) and \(T_d(\alpha,p)\) can all be obtained by standard RSB methods, but the resulting saddle-point equations can be prone to numerical instability at large \(p\) albanese2022replica?. In Sections 2.4.2 and 2.4.3, we discuss an alternative way to obtain \(T_s(\alpha,p)\) and \(T_d(\alpha,p)\).

2.3 Teacher-student setting↩︎

On our end, we study a dense Hopfield network with Hamiltonian (11 ) as a generative model for unsupervised learning. In that context, the memories \(\xi\) are model parameters that have to be trained in such a way that the examples of a given dataset \(\{\sigma^a\}_{a=1}^M\) result as typical network configurations.

In particular, we study a controlled teacher-student setting in which the examples are sampled from the probability distribution \(P \left( \sigma^a \mid \xi^* \right)\) of a so-called teacher dense Hopfield network conditioned on a single planted pattern \(\xi^* \in \left\{ -1, 1 \right\}^N\) whose entries are quenched Rademacher random variables. A student dense Hopfield network, also known as the inverse model, then samples its own student pattern \(\xi\) from the posterior distribution \[\begin{align} \label{eq:post} P \left( \xi \mid \sigma \right) & = \frac{P(\xi)\prod_{a=1}^{M} P \left( \sigma^a \mid \xi \right)}{P(\sigma)}= \frac{P(\xi)}{P(\sigma)} \prod_{a=1}^M Z^{-1}\exp\left(-\beta H[\sigma^a \mid \xi]\right)\,, \end{align}\tag{13}\] where \(P \left( \sigma^a \mid \xi \right)\) is the Gibbs distribution of the direct model with a single memory \(\xi\), and \(P \left( \xi \right)\) is the prior on \(\xi\) that is chosen to be uniform. Since the direct model has only a single pattern, \(Z\) does not depend on \(\xi\) (see Appendix 2.8), and the posterior simplifies to \[\begin{align} P \left( \xi \mid \sigma \right) & = \mathcal{Z}^{-1}(\sigma) \exp\left(-\beta H[\xi \mid \sigma]\right)\,. \end{align}\] In sum, the student posterior distribution is that of a dense Hopfield network where \(\xi\) plays the role of the sampled pattern and the examples \(\sigma\) act like the \(M\) quenched memories. Our task, called the inverse problem, consists of quantifying the student’s capability to infer the teacher pattern, which we will also call the signal. Like Gardner, we calculate a free entropy of the form \(f = \frac{1}{N} \left\langle \log \mathcal{Z} \right\rangle_{\sigma}\) in the thermodynamic limit \(N\to\infty\). This time, however, the average \(\left\langle \cdot\right\rangle_{\sigma}\) is over a structured set of examples \(\sigma\). In fact, we recall that, unlike the i.i.d. memories studied by Gardner, the examples \(\sigma^a\) are sampled from the teacher distribution \(P \left( \sigma^a \mid \xi^* \right)\).

In general, the student does not have access to the teacher generative model. In our controlled teacher-student setting, the student knows that the correct model for \(P \left( \sigma^a \mid \xi \right)\) is a dense Hopfield network. Nevertheless, it does not necessarily have access to the interaction order \(p^*\) and inverse temperature \(\beta^*\) used by the teacher. Therefore, we denote the student hyperparameters by \(p\) and \(\beta\) and emphasize that they are not necessarily equal to \(p^*\) and \(\beta^*\). As previously stated, we calculate the free entropy \[\begin{align} f = \frac{1}{N} \left\langle \log \mathcal{Z} \right\rangle_{\sigma} & = 2^{-N}\sum_{\xi^*}\sum_{\sigma} \left[ Z^* \right]^{-M} \exp\left( \beta^* \frac{p^*!}{N^{p^*-1}} \sum_{a=1}^M \sum_{i_1 < ... < i_{p^*}} \xi^*_{i_1} ... \xi^*_{i_{p^*}} \sigma^{a}_{i_1} ... \sigma^{a}_{i_p} \right)\nonumber \\ & \quad\times \log \sum_{\xi} \exp\left( \beta \frac{p!}{N^{p-1}} \sum_{a =1}^M \sum_{i_1 < ... < i_p} \xi^b_{i_1} ... \xi^b_{i_p} \sigma^{a}_{i_1} ... \sigma^{a}_{i_p} \right)\,, \end{align}\] in the thermodynamic limit \(N\to\infty\). We then draw phase diagrams of the inverse problem as a function of \(p^*\), \(T^* = 1/\beta^*\), \(p\), \(T = 1/\beta\) and \(\alpha\), where \(\alpha\) is \(M\) normalized to \(\mathcal{O} \left( 1 \right)\). Unless explicitly specified otherwise, we use \(\alpha = \frac{M p!}{N^{p - 1}}\).

2.3.1 Matched interaction orders↩︎

We first consider the case where \(p^*=p\) and the only possible mismatch between the teacher and student networks is in the inverse temperature, i.e. \(\beta^*\neq\beta\). At low \(T^*\), the student’s task is easy. In fact, below the critical temperature \(T_ {\text{crit}}\) of the direct problem with one pattern (see Fig. 10, \(\alpha = 0\) axis), the teacher produces examples \(\sigma^a\) that cluster around \(\xi^*\). Therefore, the student can infer \(\xi^*\) by aligning its pattern \(\xi\) with the examples \(\sigma^a\). This retrieval strategy works even when using a very small amount of examples (see alemanno2023hopfield?). Since the size of our dataset is extensive, the retrieval accuracy is maximum in the thermodynamic limit. We call this region the (accurate) example Retrieval phase (\(eR\)).

Conversely, when \(T^*\) is above \(T_{\text{crit}}\), the examples in the training set are very noisy and we do not observe a finite overlap between \(\sigma^a\) and \(\xi^*\) (see Fig. 10, \(\alpha = 0\) axis). In this regime, we find that the RS approximation of the \(p^* = p\) free entropy can be computed (see Appendix 2.9) in terms of the variational principle \[\begin{align} \label{eq:tsfree} \nonumber f & = \mathop{\mathrm{Extr}}_{m,k,q,r,q^*,r^*} \Bigg\{ \beta^* \beta \alpha \left[ q^{*} \right]^p - \frac{1}{2} \beta^2 \alpha q^p + \beta m^p - \beta^* \beta \alpha r^{*} q^{*} \\ & \quad+ \frac{1}{2} \beta^2 \alpha r q - \frac{1}{2} \beta^2 \alpha r - \beta m k + \frac{1}{2} \beta^2 \alpha + \log 2 \nonumber \\ & \quad+ \int dx \frac{1}{\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} x^2 \right\} \Big\langle \log \left[ \cosh \left( \beta \left[ \sqrt{\alpha r} x + \beta^* \alpha r^* + k z \right] \right) \right] \Big\rangle_z \Bigg\}\,, \end{align}\tag{14}\] whose solution is the saddle-point equations \[\begin{align} \label{eqn:saddle-point95p95s61p} q^* & = \int_{\mathbb{R}} dx \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) \Big\langle \tanh \left( \beta \left[ \sqrt{\alpha r} x + \beta^* \alpha r^* + k z \right] \right) \Big\rangle_z \,, \\ q & = \int_{\mathbb{R}} dx \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) \Big\langle \tanh^2 \left( \beta \left[ \sqrt{\alpha r} x + \beta^* \alpha r^* + k z \right] \right) \Big\rangle_z \,, \\ m & = \int_{\mathbb{R}} dx \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) \Big\langle z \tanh \left( \beta \left[ \sqrt{\alpha r} x + \beta^* \alpha r^* + k z \right] \right) \Big\rangle_z \addtocounter{equation}{1}\,, \\ r^* & = p \left[ q^{*} \right]^{p-1} \,, \\ r & = p q^{p-1} \,, \\ k & = p m^{p-1}\,, \end{align}\tag{15}\] where \(z\) is a Rademacher random variable and \(\alpha = \frac{M p!}{N^{p - 1}}\). As in the direct model described in Section 2.2, the order parameters \(m\) and \(q\) have a clear interpretation in terms of expected overlaps. \(m=\lim_{N\to\infty}\left\langle \frac{1}{N} \sum_{i} \xi_i \sigma^a_i \right\rangle_{\xi^*,\sigma,\xi}\) is the expected overlap of a retrieval attempt with an example \(\sigma^a\), and \(q=\lim_{N\to\infty}\left\langle \frac{1}{N} \sum_{i} \left\langle\xi_i\right\rangle^2_{\xi} \right\rangle_{\xi^*,\sigma}\) is the expected overlap between two retrieval attempts. Similarly, \(q^*\) is the expected overlap between the teacher and student patterns, i.e. \(q^*=\lim_{N\to\infty}\left\langle \frac{1}{N} \sum_{i} \xi^*_i\xi_i \right\rangle_{\xi^*,\sigma,\xi}\). Therefore, it is a good measure of inference performance. The free entropy (Eq. 14 ) is expected to be exact in absence of mismatch between the teacher and the student, i.e. \(\beta^*=\beta\). This condition is known as the Nishimori line nishimori1980exact?, nishimori2001statistical?, contucci2009spin?, iba1999nishimori?. Outside of the Nishimori region, RSB corrections are expected. Like the direct problem, the inverse problem with \(T^* > T_{\text{crit}}\) has different phases characterized by the values of the order parameters:

• In the Paramagnetic phase (\(P\)), the overlaps \(m\), \(q^*\) and \(q\) all vanish.

• In the signal Retrieval phases (\(lR\) and \(gR\)), \(m = 0\) but \(q^*\neq 0\) and \(q>0\). \(lR\) and \(gR\) are respectively locally stable and globally stable. In other words, local retrieval \(lR\) is only attainable from initial conditions in a limited neighborhood of \(\xi^*\), while global retrieval \(gR\) is accessible from any initial conditions given enough time. These two phases are also said to be ferromagnetic.

• In the (inaccurate) example Retrieval phase (\(eR\)), \(m\neq 0\) and \(q>0\) but \(q^*=0\).

• In the Spin-Glass phase (SG), \(q>0\) but \(q^*\) and \(m\) vanish.

In sum, when \(T^*\) is above \(T_{\text{crit}}\), the student can only learn the teacher pattern in the signal retrieval phases. In all the other phases, the student pattern is uncorrelated with the signal, being either a random guess (\(P\) phase), aligned with a noisy example (inaccurate \(eR\) phase), or aligned with a spurious low energy state (\(SG\) phase). We stress that we cannot have \(m \neq 0\) and \(q^* \neq 0\) at the same time (accurate \(eR\) phase) when \(T^* > T_{\text{crit}}\) because \(\lim_{N\to\infty}\left\langle \frac{1}{N} \sum_{i} \xi^*_i \sigma^a_i \right\rangle_{\xi^*,\sigma} = 0\) in that regime (see Fig. 10, \(\alpha = 0\) axis).

2.3.2 Mismatched interaction orders↩︎

We also investigate the \(T^* > T_{\text{crit}}\) regime in the presence of a mismatch between the interaction orders of the teacher and student networks, i.e. \(p^*\neq p\). We focus on the case of \(p^* = 2\) and even \(p \geq 3\) to study the consequences of fitting the teacher of alemanno2023hopfield? using a student with higher order interactions. We find two different scaling regimes of the training set size \(M\) and inverse temperature \(\beta^*\) that make retrieval possible (see Appendix D):

• a large-noise scaling where \(\beta^* \sim \mathcal{O} \left( N^{2/p - 1} \right)\) and \(M \sim \mathcal{O} \left( N^{p - 1} \right)\), such that \(\alpha = \frac{M p!}{N^{p-1}}\) and \(\lambda = \frac{\left[ \beta^* \right]^{p/2}}{(p/2)!} N^{p/2 - 1}\) are finite;

• a finite-noise scaling where \(\beta^* \sim \mathcal{O} \left( 1 \right)\) and \(M \sim \mathcal{O} \left( N^{p/2} \right)\), such that \(\alpha = \frac{M (p/2 + 1)!}{N^{p/2}}\) is finite.

In the large-noise scaling, we obtain saddle point equations similar to Eqs. (15 ) but with \(\beta^*\) replaced by \(\lambda\) (see Appendix 2.9). Conversely, the finite noise scaling leads to \[\begin{align} \label{eqn:saddle-point95p95s612} q^* & = \Big\langle \tanh \left( \beta \left[ \eta \alpha r^* + k z \right] \right) \Big\rangle_z \,, \\ m & = \Big\langle z \tanh \left( \beta \left[ \eta \alpha r^* + k z \right] \right) \Big\rangle_z \addtocounter{equation}{1}\,, \\ r^* & = p \left[ q^{*} \right]^{p-1} \,, \\ k & = p m^{p-1}\,, \end{align}\tag{16}\] where \(\eta\) generally depends on \(\beta^*\) and \(p\) in a nontrivial way, but we find that \(\eta = \frac{2 \left[ \beta^* \right]^2}{\left( 1 - 2\beta^* \right)^2}\) when \(p = 4\) (see Appendix 2.9). These equations can also be derived by extrapolating the large-noise equations to \(\alpha_{\text{large noise}} \rightarrow 0\) and \(\lambda \rightarrow \infty\) with fixed \(\lambda \alpha_{\text{large noise}} = \eta \alpha_{\text{finite noise}}\).

2.4 Results and Discussion↩︎

2.4.1 Retrieval transition at large interaction order↩︎

The paramagnetic solution of Eqs. (15 ) always exists and is globally stable in the part of the phase diagram where the temperature \(T\) is relatively large and \(\alpha = \frac{M p!}{N^{p - 1}}\) is relatively small. On the other hand, the \(gR\) phase exists when \(\beta^2 \alpha p\) and \(\beta^* \beta \alpha p\) are both large. In fact, in that limit, \(q^* = q = 1\) is a fixed point of Eqs. (15 ). The critical line where \(gR\) becomes globally stable instead of \(P\) is not clear from this analysis alone, but we can at least find it analytically in the limit of infinite \(p\). As for the direct model, the free entropy and the total entropy of the paramagnetic phase are respectively \(\frac{1}{2} \beta^2 \alpha + \log 2\) and \(-\frac{1}{2} \beta^2 \alpha + \log 2\) gardner1987multiconnected?. At the same time, the \(p \rightarrow \infty\) free entropy takes the form \[\begin{align} f & = \mathop{\mathrm{Extr}}\Bigg\{ \beta^* \beta \alpha \;\theta\left( q^{*} - 1 \right) - \frac{1}{2} \beta^2 \alpha \;\theta\left( q - 1 \right) - \beta^* \beta \alpha r^{*} q^{*} + \frac{1}{2} \beta^2 \alpha r q - \frac{1}{2} \beta^2 \alpha r + \frac{1}{2} \beta^2 \alpha \\ & \quad+ \log 2 + \int dx \frac{1}{\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} x^2 \right\} \log \left[ \cosh \left( \sqrt{\beta^2 \alpha r} x + \beta^* \beta \alpha r^* \right) \right] \Bigg\}\,, \end{align}\] where \(\theta \left( q - 1 \right) := \lim_{p \rightarrow \infty} q^p\), \(q \in \left[ 0, 1 \right]\), is the Heaviside step function jumping at \(q=1\), i.e. \(\theta(1)=1\) and \(\theta(q)=0\) \(\forall q\in[0,1)\). In this limit, the ferromagnetic phase is characterized by \(q=q^*=1\), and its free entropy is then \[\begin{align} f & = \beta^* \beta \alpha - \beta^* \beta \alpha p + \int dx \frac{1}{\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} x^2 \right\} \log \left[ 2 \cosh \left( \sqrt{\beta^2 \alpha p} x + \beta^* \beta \alpha p \right) \right] \\ & \approx \beta^* \beta \alpha - \beta^* \beta \alpha p + \int dx \frac{1}{\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} x^2 \right\} \left( \sqrt{\beta^2 \alpha p} x + \beta^* \beta \alpha p \right) \\ & = \beta^* \beta \alpha\,. \end{align}\] The corresponding total entropy is \(s = f - \beta \frac{\partial f}{\partial \beta} = 0\), as expected from a ferromagnetic phase with \(q^* = q = 1\). On the Nishimori line, \(f = \beta^* \beta \alpha\) becomes larger than the free entropy of the paramagnetic phase, which triggers a phase transition, if and only if \[\begin{align} \label{eq:crlinelargep} T & < \sqrt{\frac{\alpha}{2 \log 2}}\,, \end{align}\tag{17}\] where \(T_E= \sqrt{\frac{\alpha}{2 \log 2}}\) is also the temperature below which the total entropy of the paramagnetic phase becomes negative. Outside of the Nishimori line, this inequality generalizes to \(\beta^* \beta \alpha > \frac{1}{2} \beta^2 \alpha + \log 2\), leading to \[\begin{align} \beta^* - \sqrt{\left[ \beta^* \right]^2 - \frac{2 \log 2}{\alpha}} < \beta & < \beta^* + \sqrt{\left[ \beta^* \right]^2 - \frac{2 \log 2}{\alpha}}\,, \end{align}\] while the temperature where the paramagnetic total entropy becomes negative stays the same.

2.4.2 Transition to the ordered phases: Universality↩︎

In the \(p \rightarrow \infty\) limit, the transition towards \(gR\) of the inverse model on the Nishimori line is identical to the exact \(P\)-\(SG\) transition of the direct model gardner1987multiconnected?. We claim that these two critical lines are actually closely related for any \(p\). In the Hopfield model with \(p=2\), they were already shown to be identical alemanno2023hopfield?. We will now argue that they overlap for any \(p\) and \(\beta\) such that \(T > T_{\text{crit}}\) (see Figs. 11 and 10). In the case of \(p=2\), both lines can be obtained exactly from the RS approximation of either the direct model or the inverse model, so there is no obvious advantage to using this equivalence in calculations. In general, while the inverse problem on the Nishimori line is replica symmetric, the direct problem is not, and the \(p\geq 3\) replica symmetric \(P\)-\(SG\) transition is not exact. Moreover, even the critical line calculated using 1RSB may be inaccurate due to numerical instability albanese2022replica?. In this situation, the knowledge of the \(gR\) transition in the replica-symmetric inverse problem can be used to locate the exact \(P\)-\(SG\) transition of the direct problem, where symmetry breaking occurs.

For that purpose, we will argue that, given \(T > T_{\text{crit}}\), the direct model is in the paramagnetic phase if and only if the inverse model is in the paramagnetic phase.

The converse implication comes from the fact that since (see Appendix 2.8) \[P \left( \sigma \right) = \frac{1}{2^{M N}} \frac{\mathcal{Z}(\sigma)}{\left\langle \mathcal{Z} \right\rangle}\,,\] the example distribution \(P \left( \sigma \right)\) of the inverse problem is contiguous roussas1972contiguity? to the uniform distribution, i.e. the memory distribution of the direct problem, when \[\lim_{N \rightarrow \infty} \left\{ \frac{\log \mathcal{Z} - \log \left\langle \mathcal{Z} \right\rangle}{N} \right\} = 0\,.\] As determined in Appendix 2.8 and 2.9, the annealed expression \(\frac{1}{N} \log \left\langle \mathcal{Z} \right\rangle\) is equal to the free entropy of the paramagnetic phase. Therefore, when the inverse model is in the paramagnetic phase, \(P \left( \sigma \right)\) is contiguous to the uniform distribution. This property is called quiet planting and is known to occur more generally in mean-field paramagnets achlioptas2008algorithmic?, krzakala2009hiding?, zdeborova2011quiet?, zdeborova2016statistical?. In our problem setting, it means that if the inverse model is in the paramagnetic phase, then it is equivalent to the direct model. In particular, if the inverse model is in the paramagnetic phase, then so is the direct model. In more intuitive terms, the \(gR\) transition temperature of the inverse model must be greater than or equal to the \(P\)-\(SG\) transition temperature of the direct model because the ensemble of examples \(\sigma^a\) generated by the teacher model is on average at least as structured as the set of i.i.d. random memories stored in the direct model.

For the direct implication, notice that the average replicated partition function of the direct model in the paramagnetic phase can be approximated as (see Appendix 2.10) \[\begin{align} \left\langle Z^L \right\rangle & \approx \frac{1}{\left\langle Z \right\rangle} \Bigg\langle \sum_{\sigma} \exp \Bigg( \beta N \sum_{\gamma} \sum_{\mu \in \Gamma_{\gamma}} \left[ \frac{1}{N} \sum_i \xi^{\mu}_i \sigma^{\gamma}_i \right]^p \\ & \quad+ \beta \sum_{\gamma} \sum_{\mu \in \Bar{\Gamma}} \frac{p!}{N^{p-1}} \sum_{i_1 < ... < i_p} \xi^{\mu}_{i_1} ... \xi^{\mu}_{i_p} \sigma^{\gamma}_{i_1} ... \sigma^{\gamma}_{i_p} \Bigg) \\ & \quad \sum_{\sigma_0} \exp \left( \beta \sum_{\mu \in \Bar{\Gamma}} \frac{p!}{N^{p-1}} \sum_{i_1 < ... < i_p} \xi^{\mu}_{i_1} ... \xi^{\mu}_{i_p} \sigma^{0}_{i_1} ... \sigma^{0}_{i_p} \right) \Bigg\rangle\,. \end{align}\] This expression is identical to the replicated partition function of the inverse model with \(T > T_{\text{crit}}\), which therefore must also be in the paramagnetic phase.

As a consequence, when \(T > T_{\text{crit}}\), the \(P\)-\(SG\) transition line of the direct model must be identical to the \(gR\) transition line of the inverse model on the Nishimori line.

2.4.3 Phase diagram on the Nishimori line↩︎

On the Nishimori line, the student is fully informed about the teacher generative model and uses \(\beta=\beta^*\) and \(p=p^*\). In this scenario, thanks to the Nishimori identities nishimori2001statistical?, it is well known that \(\xi^*\) and \(\xi\) play symmetric roles and that \(q^* = q\). For the same reason, the overlaps \(\frac{1}{N} \sum_i \xi^*_i\xi_i\) and \(\frac{1}{N} \sum_i \xi^1_i\xi^2_i\) have the same distribution. From the self-averaging of \(\frac{1}{N} \sum_i \xi^*_i\xi_i\), it follows that the system is expected to be replica symmetric, and Eqs. (14 ) and (15 ) are expected to hold. Fig. (11) shows the phase diagrams obtained by solving the saddle-point equations numerically on the Nishimori line. Both \(q^*=q\) and the replica symmetry condition are verified. In particular, numerical solutions of a few values of \(p \geq 3\) show that the \(gR\) transition occurs at a higher \(T\) than the line \(\beta^2 \alpha = 2 \log 2\) where the total entropy of the paramagnetic phase becomes negative. In other terms, the phase transition towards \(gR\) prevents the total entropy from becoming negative when \(T\) decreases below \(\sqrt{\frac{\alpha}{2\log2}}\), which is consistent with the RS solution being exact on the Nishimori line.

At low \(T\), the student can learn efficiently within the accurate \(eR\) regime. In this phase, learning is possible (\(q^*\neq0\)) because the examples are correlated with the signal and the student can retrieve it by simply being aligned with them (\(m\neq 0\)).

At high \(T\), learning is possible only if the amount of examples, i.e. the size of the dataset, is sufficiently large. When \(\alpha\) is too small, Eqs. (15 ) have only a paramagnetic fixed point because the amount of information carried by the dataset is not large enough. Numerical solutions suggest that the paramagnetic fixed point always exist and it is actually locally stable in the whole high-temperature regime. When \(\alpha\) is sufficiently large, the signal retrieval fixed point appears as a locally stable attractor (\(lR\) phase). It becomes globally stable (\(gR\) phase) as the size of the dataset is increased further or the student temperature decreases.

As per the previous Section, the critical boundary of the \(gR\) phase obtained by solving Eqs. 15 is identical to the 1RSB \(P\)-\(SG\) transition temperature \(T_s(\alpha,p)\) of the direct model. Similarly, we observe that the metastable \(lR\) phase coincides with the d1RSB phase of the direct model (see Fig. 11). Our results are also consistent with the fact that \(T_s(\alpha,p)\to T_E(\alpha)\) in the \(p \rightarrow \infty\) limit. In fact, we find that the analytical limit boundary closely agrees with the numerical solution of the saddle-point equations with \(p^* = p = 10\) and remains a good approximation even down to \(p^* = p = 4\).

In the student model, \(\sigma\) plays a similar role as the weights of the trainable dense Hopfield network model that K & H designed for classification of data krotov2016dense?. In that context, \(\xi\) is analogous to the test data whose labels are being predicted (see Fig. 12). In fact, the computation performed by K & H’s model to recover labels is similar to the update rule used by the student to infer the teacher pattern (see Appendix 2.6). Moreover, the \(eR\) and \(gR\) phases are respectively reminiscent of the prototype and feature regimes of K & H’s networks. Therefore, we believe that the student can act as a toy model of label prediction in these two regimes.

Comparing instead the phase diagrams of our inverse model with that of the inverse \(2\)-body Hopfield model, we see that the \(eR\) and \(gR\) phases of the inverse \(p\)-body model with \(p \geq 3\) are respectively analogous to the \(eR\) and \(sR\) (signal Retrieval) phases presented in alemanno2023hopfield?. One of the key differences between \(p = 2\) and \(p \geq 3\) is that the paramagnetic to signal retrieval phase transition of the \(p\)-body model is second order for \(p = 2\) but first order for \(p \geq 3\). On the one hand, the second order phase transition of \(p = 2\) indicates that its paramagnetic fixed point is never locally stable and sets an unambiguous boundary between the \(sR\) phase where \(\xi^*\) can be recovered starting from any initial conditions and the paramagnetic phase where pattern retrieval is impossible zdeborova2016statistical?. On the other hand, the first order phase transition of \(p \geq 3\) allows the retrieval and paramagnetic regimes to coexist. The \(lR\) phase is locally stable precisely because it coexists with the paramagnetic phase and has a lower free entropy. Meanwhile, the \(gR\) phase also coexists with the paramagnetic phase, but has a larger free entropy. In the presence of phase coexistence, an algorithm trying to retrieve \(\xi^*\) starting from random initial conditions can get stuck in the paramagnetic phase instead. In fact, it has been conjectured that there is no algorithm with random initial conditions that can find such a ferromagnetic fixed point in a tractable amount of time zdeborova2016statistical?, antenucci2019glassy?. That kind of metastable region was thus given the name hard phase zdeborova2016statistical?, zdeborova2007phase?. In summary, we expect that \(p \geq 3\) models in the \(gR\) phase can only recover partially corrupted patterns whereas \(p = 2\) can recover them entirely.

Fig. (13) shows results from Monte Carlo simulations with \(p = 3\), where \(L\) replicas of the student pattern \(\{\xi^b\}_{b=1}^L\) are initialized to the teacher pattern \(\xi^*\) corrupted by some Rademacher noise \(\varepsilon\). In other words, the initial values of \(\xi^b_i\) are sampled from the distribution \(\left( 1 - \varepsilon \right) \delta \left(\xi_i - \xi^*_i \right) + \frac{\varepsilon}{2} \left[ \delta \left( \xi_i+1 \right) + \delta \left( \xi_i - 1 \right) \right]\) with \(\varepsilon \in [0, 1]\). The value of \(\varepsilon\) is tuned so that the simulations start relatively close to the saddle-point solutions. As explained previously, \(gR\) is a hard phase, so this initialization is necessary to make \(\xi^b\) converge to \(gR\) in a reasonable amount of time. Additionally, it is also used to make \(\xi^b\) converge to the \(lR\) phase rather than the \(P\) phase when desired. Once the simulations are over, the overlaps are averaged over all \(L\) replicas. If we fix \(\varepsilon = 0\), then the simulations generally converge to the \(lR\) phase when it is a fixed point. If instead we initialize them to the saddle-point solutions by handpicking \(\varepsilon\), then they stay near the initial overlaps. In either case, the simulations converge to \(eR\) when it is globally stable. Some simulation data points might be systematically shifted up with respect to the saddle-point solutions. However, this difference decreases with the system size \(N\), so finite size effects seem sufficient to explain it (see Fig. 18 in Appendix 2.11). Overall, the Monte Carlo simulations are in very good agreement with the \(p = 3\) overlap landscape obtained by solving the saddle-point equations numerically.

2.4.4 Inference temperature vs dataset noise↩︎

In the two next Sections, we will discuss the phase diagram when the student is only partially informed about the teacher generative model, i.e. when the Nishimori conditions do not hold. We start with the case where \(p=p^*\) but \(\beta\neq\beta^*\), i.e. the inference temperature \(T\) is different from the dataset noise \(T^*\). As we argued in Section 2.3.1, the student accurately retrieves \(\xi^*\) when \(T^* < T_ {\text{crit}}\). On the other hand, we must solve the saddle-points equations (see Eqs. 15 ) to study \(T^* > T_ {\text{crit}}\).

We show the phase diagram of this region on Fig. (14). At high inference temperature \(T\), the situation is similar to Fig. (11): retrieval is possible if the data load \(\alpha\) is sufficiently large, but the paramagnetic phase is always locally stable. The situation is different when the inference temperature is low. In that case, there are two phases that we did not see for \(\beta = \beta^*\): the inaccurate \(eR\) phase and the \(SG\) phase. When \(\alpha\) is relatively small, the student falls in the inaccurate \(eR\) phase. In this regime, it has finite overlap with one of the noisy examples and cannot retrieve the signal \(\xi^*\). When \(\alpha\) is larger, the interference among the noisy examples prevents the student to be aligned with them. In this regime, the \(SG\) phase, the student locally converge to spurious patterns that are uncorrelated with the signal.

Accurate pattern retrieval is only possible in the \(lR\) and \(gR\) phases where \(\alpha\) is so large that the student can gather enough information from the dataset to become very close to \(\xi^*\). The phase diagrams indicate that pattern retrieval is optimal on the Nishimori line in the sense that \(\beta = \beta^*\) is the inverse temperature where the student needs the least examples to recover \(\xi^*\). In other words, the student’s performance is nonmonotonic in \(T\) and peaks at \(T = T^*\). These properties were also observed in the teacher-student setting of the \(p = 2\) Hopfield network alemanno2023hopfield?.

Contrary to what one would expect to see on the exact phase diagram nishimori1980exact?, nishimori2001statistical?, the Nishimori line \(T=T^*\) does not to cross a triple point on the RS phase diagram. The issue is that the RS phase diagram is not exact outside of the Nishimori line. In particular, the \(SG\) phase boundary is not exact. Outside of the retrieval regime, the free entropy of the inverse model is the same as the direct model. Since the transition towards \(gR\) of the inverse model on the Nishimori line overlaps the exact \(P\)-\(SG\) transition of the direct model (see Section 2.4.3), we deduce that it must also overlap the exact \(P\)-\(SG\) transition of the inverse model outside of the \(gR\) phase. Plotting it on the RS phase diagrams, we see that it indeed crosses the Nishimori line and the \(gR\) phase boundary at the same point, which therefore becomes a triple point, as expected.

2.4.5 Interaction order and noise tolerance↩︎

So far, we assumed that the student is informed about the interaction order used by the teacher, i.e. \(p=p^*\). In this Section, we investigate the role of the student’s choice of \(p\) when the task is to learn from a dataset sampled by a \(2\)-body Hopfield network, i.e. \(p^*=2\). We study two different non trivial scalings regimes of \(M\) and \(\beta^*\) that make pattern inference possible (see Appendix 2.9).

2.4.5.1 Large noise scaling

We first consider a large noise scaling where \(\beta^* \sim \mathcal{O} \left( N^{2/p - 1} \right)\) and \(M \sim \mathcal{O} \left( N^{p - 1} \right)\), such that \[\alpha = \frac{M p!}{N^{p-1}}\,, \quad and \quad \lambda = \frac{\left[ \beta^* \right]^{p/2}}{(p/2)!} N^{p/2 - 1} \,,\] are finite. In this scaling, a \(p \geq 3\) network requires \(\mathcal{O} \left( N^{p - 2} \right)\) more training examples than a \(p = 2\) network with finite load \(\gamma = \frac{M}{N}\), but also has a higher tolerance to teacher noise. For instance, a student with \(p = 4\) interactions is able to retrieve the pattern of a teacher with \(T^* \sim \mathcal{O} \left( N^{1/2} \right)\) noise when it is shown enough examples \(M \sim \mathcal{O} \left( N^3 \right)\) to be in the \(gR\) phase (see Fig. 15).

\(\mathcal{O} \left( N^{1/2} \right)\) noise tolerance was also observed in the \(p = 4\) direct model, where it is a consequence of the redundancy stemming from storing \(\mathcal{O} \left( N \right)\) memories rather than the \(\mathcal{O} \left( N^3 \right)\) needed to saturate the storage capacity agliari2020neural?.

Our \(p = 4\) inverse model exploits a different kind of redundancy by learning from \(\mathcal{O} \left( N^3 \right)\) examples whereas \(p = 2\) only needs \(\mathcal{O} \left( N \right)\). In other terms, both storing extensively less memories than the maximum allowed amount and generating extensively more examples than the minimum required amount provide enough redundancy to recover a pattern muddled in an extensive amount of noise. In both cases, there is an \(\mathcal{O} \left( N^2 \right)\) gap between the number of patterns used in the noise-tolerant and noise-susceptible regimes. Going beyond \(p = 4\), the inverse model has \(\mathcal{O} \left( N^{1 - 2/p} \right)\) noise tolerance as a function of \(p\). In particular, our theory predicts that the tolerance saturates at \(T^* \sim \mathcal{O} \left( N \right)\) as \(p \rightarrow \infty\), but at the cost of using an intractable number of examples. This behavior is different from the \(\mathcal{O} \left( N^{1/2 - p/4} \right)\) tolerance of the direct \(p\)-body model in the noisy-learning regime studied in agliari2020tolerance?. In other terms, the dataset noise that we are facing is of a different nature than the learning noise of agliari2020tolerance?. In any case, it is interesting that both the direct and inverse models are able to tolerate an extensive amount of noise. Overall, our results suggest that it could be advantageous to use a student network with a relatively large \(p\) to learn from a large but noisy dataset when the \(p^*\) of the teacher generative model is unknown.

An unavoidable drawback of large teacher noise is that it always lead to uncorrelated examples, which makes accurate example retrieval impossible. Instead, it is replaced by the inaccurate example retrieval phase where the student has finite overlap \(m\) with a noisy example generated by the teacher but no overlap with the signal (see Fig. 15). Depending on \(T\) and \(\alpha\), this phase can be either globally stable or locally stable. For the sake of clarity, we plot only the globally stable phase on our phase diagram in Fig. (15). The locally stable phase is arguably less important to plot because it is identical to the locally stable ferromagnetic phase previously reported in the direct model when assuming replica symmetry (see albanese2022replica? and Fig. 10).

Given \(m = 0\), the free entropy of the inverse model with \(p \geq 3\), \(p^* = 2\) and \(\beta = \lambda\) is the same as on the Nishimori line (see Eq. 15 and Appendix 2.9). As a direct consequence, the total entropy is positive outside of the \(eR\) phase (see Fig. 15). Additionally, the \(p^* = 2\), \(p \geq 3\) phase diagrams with \(\beta \neq \lambda\) are identical to the \(p = p^*\) phase diagrams with \(\beta \neq \beta^*\), which suggests that \(\beta = \lambda\) is optimal for \(p^* = 2\), \(p \geq 3\) in the same sense as \(\beta = \beta^*\) is optimal for \(p = p^*\) (see Fig. 14). Monte Carlo simulations confirm that a student with \(p \geq 3\) is able to retrieve the pattern of a teacher with \(p = 2\) and \(T^* \sim \mathcal{O} \left( N^{1/2} \right)\) (see Fig. 16). However, the \(lR\) phase transition is at a higher \(T\) in the simulations than on the \(\beta = \lambda\) RS phase diagram (see Fig. 14), which means that RSB is necessary to describe it accurately. One could check where replica symmetry holds by evaluating the stability of the RS saddle point throughout the phase diagram.

2.4.5.2 Finite noise scaling

We also consider a different scaling regime where \(\beta^* \sim \mathcal{O} \left( 1 \right)\) and \(M \sim \mathcal{O} \left( N^{p/2} \right)\), such that \[\alpha = \frac{M (p/2 + 1)!}{N^{p/2}}\,,\] is finite. In this finite-noise scaling, \(p \geq 3\) requires \(\mathcal{O} \left( N^{p/2 - 1} \right)\) more training examples than \(p = 2\), which is a lot less than the first scaling. For instance, a student with \(p = 4\) needs \(\mathcal{O} \left( N^2 \right)\) examples to retrieve \(\xi^*\). As before, the phase transitions are all first order, the overlap \(q^*\) stays high throughout the \(gR\) and \(lR\) phase of \(p = 4\) and \(gR\) is a hard phase. The saddle-point equations (see Eqs. 16 ) are free from the pattern interference term \(\sqrt{\alpha r} x\) present in their \(p^* = p\) counterparts (see Eqs. 15 ) until \(\beta^*\) becomes so small that is approaches \(\mathcal{O} \left( N^{2/p - 1} \right)\). Therefore, contrary to \(p^* = p = 2\), the network is never in the \(SG\) phase. Practically, it means that \(p \geq 3\) gives more freedom than \(p = 2\) for tuning \(\beta\) and \(\alpha\). The only remaining restriction is that choosing \(\alpha\) and \(T\) too small puts the network into the inaccurate \(eR\) phase resulting from the \(k z\) term (see Fig. 15). The saddle point equations can be derived without the RS ansatz because they do not involve \(q\) and \(r\). Consequently, we expect them to yield an exact solution. Like on the Nishimori line, the total entropy of the paramagnetic phase is always positive, which is consistent with the solution being exact.

2.4.6 Robustness against adversarial attacks↩︎

Inverse models with \(p^* = 2\) and \(p \geq 3\) offer an opportunity to study adversarial attacks in a simple setting because their phase diagrams have regions where the signal retrieval phases (\(gR\) and \(lR\)) overlap with the inaccurate \(eR\) phase. Recall that, in the \(lR\) phase, a noisy student pattern \(\xi\) either converges to \(\xi^*\) or falls in the paramagnetic phase, depending on the amount of noise that \(\xi\) contains initially. The quantity of noise needed to prevent pattern retrieval becomes smaller as one approaches the \(lR\) to \(P\) phase transition and the basin of attraction of \(lR\) shrinks. Similarly, in the region of inaccurate \(eR\) where signal retrieval is metastable, patterns \(\xi\) that are corrupted by replacing some of their entries \(\xi_i\) by the components \(\sigma^a_i\) of an example \(\sigma^a\) may converge to \(\sigma^a\) when enough entries are replaced. The fraction \(\varepsilon\) of entries that need to be replaced becomes smaller as the basin of attraction of inaccurate \(eR\) expands and overtakes that of signal retrieval. In practice, an adversary can use this strategy to trick the student into converging to a pattern other than \(\xi^*\). This scenario is similar to an adversarial attack targeting the input of K & H’s dense Hopfield network model because the student pattern \(\xi\) plays a similar role in the inverse model as the test data in K & H’s dense Hopfield networks (see Fig. 12, Section 2.4.3 and Appendix 2.6). In that analogy, the examples \(\sigma\) are acting like the neural network weights rather than taking the role of the training data.

We will now investigate what values of the perturbation size \(\varepsilon\) are a threat by deriving a formula for the largest \(\varepsilon\) such that the student converges to the signal at zero temperature. This largest \(\varepsilon\) will be denoted \(\varepsilon^*\), and we expect it to be a good measure of adversarial robustness. The saddle-point equations with \(T = 0\) indicate that the student converges to one of the signal retrieval phases if and only if \(k < \eta \alpha r^*\) (see Eqs. 16 ). Sampling the initial conditions of \(\xi_i\) from \(\left( 1 - \varepsilon \right) \delta \left( \xi_i-\xi^*_i \right) + \varepsilon \;\delta \left( \xi_i-\sigma^a_i \right)\) with \(\varepsilon \in \left[ 0, 1 \right]\), we get \[\begin{align} r^* & = p \left[ \frac{1}{N} \sum_{i = 1}^{\left( 1 - \varepsilon \right) N} \xi^*_i \xi^*_i + \frac{1}{N} \sum_{i = 1}^{\varepsilon N} \xi^*_i \sigma^a_i \right]^{p - 1}\,, \\ \quad k & = p \left[ \frac{1}{N} \sum_{i = 1}^{\left( 1 - \varepsilon \right) N} \xi^*_i \sigma^a_i + \frac{1}{N} \sum_{i = 1}^{\varepsilon N} \sigma^a_i \sigma^a_i \right]^{p - 1}\,. \end{align}\] By the law of large numbers, \(\frac{1}{\varepsilon N} \sum_{i = 1}^{\varepsilon N} \xi^*_i \sigma^a_i\) and \(\frac{1}{\left( 1 - \varepsilon \right) N} \sum_{i = 1}^{\left( 1 - \varepsilon \right) N} \xi^*_i \sigma^a_i\) are both typically close to \(m^* = \frac{1}{N} \sum_i^N \xi^*_i \sigma^a_i \approx 0\) as \(N \rightarrow \infty\). If we take \(\sigma^a\) to be a typical example, then \(r^*\) and \(k\) reduce to \[\begin{align} r^* & \approx p \left( 1 - \varepsilon \right)^{p - 1}\,, \\ k & \approx p \varepsilon^{p - 1}\,. \end{align}\] Substituting these expressions back in \(k < \eta \alpha r^*\) yields \[\begin{align} \varepsilon^{p - 1} & < \eta \alpha \left( 1 - \varepsilon \right)^{p - 1} \,, \\ \varepsilon & < \frac{\left[ \eta \alpha \right]^{\frac{1}{p - 1}}}{\left[ \eta \alpha \right]^{\frac{1}{p - 1}} + 1}\,. \end{align}\] In other terms, the inverse model with \(p^* = 2\) and even \(p \geq 3\) is resistant to adversarial attacks of size \(\varepsilon^* = \frac{\left[ \eta \alpha \right]^{\frac{1}{p - 1}}}{\left[ \eta \alpha \right]^{\frac{1}{p - 1}} + 1}\) and smaller. For \(p = 4\), \(\varepsilon^*\) is in good agreement with Monte Carlo simulations of the inverse model corrupted by a typical example (see Fig. 17). This comparison is good evidence that our solution of the finite-noise scaling is indeed exact. Additionally, \(\varepsilon^*\) is a decent approximation of empirical robustness even when the inverse model is corrupted by the example that has the largest overlap with \(\xi^*\). A similar construction with the perturbation sampled uniformly at random gives \(k \sim \mathcal{O} \left( N^{1/2 - p/2} \right) \approx 0\), so adversarial attacks are much more efficient at fooling the model than random noise. Just like adversarial attacks targeting more complicated neural networks biggio2013evasion?, szegedy2013intriguing?, our example-based attack can be hard to detect at low \(\varepsilon\) because a few adversarially perturbed entries \(\xi_i\) do not look very different from a low amount of meaningless noise.

Moreover, \(\varepsilon^*\) grows monotonically with \(\alpha\), which is consistent with the common observation that larger neural networks are also more adversarially robust madry2018towards?, gowal2020uncovering?, huang2021exploring?, bubeck2021law?, bubeck2021universal?, puigcerver2022adversarial?, ribeiro2023overparameterized?. At first glance, this effect can be counter-intuitive because adversarial vulnerability looks like a form of overfitting goodfellow2015explaining?. In our model, however, all examples work together to stabilize the \(lR\) phase, and the best way to push the student into the \(eR\) phase is to perturb it with a single example. Therefore, it is not surprising that increasing \(\alpha\) makes the student more robust. We recall that the examples \(\sigma\) are a feature-based representation of \(\xi^*\). Interestingly, it means that the underlying mechanism of our example-based attack is conceptually similar to gradient-based attacks targeting many common types of neural networks goodfellow2015explaining?. In fact, gradient-based attacks find features stored in neural network weights and add them to the data in order to fool the network goodfellow2015explaining?, jetley2018friends?, ilyas2019adversarial?, tsipras2019robustness?. It would be interesting to investigate, both empirically and theoretically, if only a small number of weights are involved in constructing these adversarial attacks. If it is the case, it could explain why larger neural networks are often more robust. In general, we expect this kind of one-example attack to be possible in any region of signal retrieval that overlaps with the inaccurate \(eR\) phase. Using \(p \neq p^*\) may not be a necessary ingredient of adversarial vulnerability in more general models with other sources of mismatch, but in our case it ensures that the signal retrieval phases intersect the inaccurate \(eR\) phase. Conversely, the accurate \(eR\) phase is by definition robust to adversarial attacks since retrieving an example \(\sigma^a\) is the same as recovering \(\xi^*\). This distinction clarifies why the dense Hopfield networks designed by K & H are adversarially robust in the prototype phase despite being adversarially vulnerable in the feature phase. In fact, K & H observed that adversarial attacks are unsuccessful in the prototype phase specifically because they retrieve stored examples that are semantically meaningful krotov2018dense?. In summary, our model yields two main results concerning adversarial examples. First of all, it suggests a reason why large feature-based neural networks are more adversarially robust than smaller ones. Second of all, it clarifies why dense Hopfield networks are much more robust in the prototype phase than in the feature phase.

2.5 Conclusion↩︎

In this work, we derive the exact phase diagram of the \(p\)-dense networks in the teacher-student setting barra2017phase?, barra2018phase?, gardner1987multiconnected?, alemanno2023hopfield?. On the Nishimori line, we find an example Retrieval phase (\(eR\)) and a global Retrieval phase (\(gR\)) reminiscent of the prototype and feature regimes observed empirically in dense Hopfield networks krotov2016dense?. We show that the phase transition towards \(gR\) of the inverse model overlaps the Paramagnetic to Spin-Glass (\(P\)-\(SG\)) transition of the direct model, which allows us to locate the \(P\)-\(SG\) transition much more precisely than before gardner1987multiconnected?, albanese2022replica?. On the other hand, we discover that inverse models outside of the Nishimori line are able to resist an extensive amount of noise. In fact, a student with \(p \geq 3\) is able to learn from a teacher with \(p^* = 2\) even when the teacher’s inverse temperature \(\beta^*\) is as low as \(\mathcal{O} \left( N^{2/p - 1} \right)\). Moreover, such a student is immune to pattern interference until \(\beta^*\) reaches \(\mathcal{O} \left( N^{2/p - 1} \right)\). In this setting, we derive a formula measuring the adversarial robustness of the student with \(p \geq 3\) and \(T = 0\). We then use this formula to describe how making a neural network larger can potentially increase its robustness to adversarial attacks constructed with only a few learned weights madry2018towards?, gowal2020uncovering?, huang2021exploring?, bubeck2021law?, bubeck2021universal?, puigcerver2022adversarial?, ribeiro2023overparameterized?. Our model also clarifies why the prototype phase of dense Hopfield networks is adversarially robust krotov2018dense?. We compare our key results against Monte Carlo simulations.

Dense networks with exponential interactions have been argued to be the \(p \rightarrow \infty\) limit of the \(p\)-body models demircigil2017model?. It would be interesting to see if they can achieve \(\mathcal{O} \left( N \right)\) noise tolerance at the cost of an exponential number of training examples. More generally, studying exponential models in the teacher-student setting would be an interesting extension of this work and could be used to complement existing studies of the direct model demircigil2017model?, lucibello2023exponential?. A caveat of our model is that the teacher has only one pattern. In fact, we would need to use a teacher with at least two patterns to describe more completely the kind of adversarial attack aiming to misclassify data. It should be possible to study this kind of teacher by using an approach similar to hou2019minimal?. In particular, barra2017phase? and hou2019minimal? argue that the performance of restricted Boltzmann machines with a finite number \(P\) of i.i.d. planted patterns is independent of \(P\) in the teacher-student setting. It would be interesting to investigate whether this characteristic also holds for \(p\)-body dense networks. On the practical side, we highlight the untapped benefits of using \(p\)-body models to either resist an extensive amount of noise in the feature phase or improve adversarial robustness in the prototype phase. Overall, we stress that further investigations of dense Hopfield networks could unlock their true potential.

Data availability↩︎

The figures can be reproduced using the code available on this public Github repository.

3 Modeling structured data learning with Restricted Boltzmann machines in the teacher–student setting↩︎

3.1 Introduction↩︎

Restricted Boltzmann Machines (RBMs) ackley1985learning?, smolensky1986information?, freund1991unsupervised?, hinton2002training? are empirically known to fit complicated data and then sample new instances that are faithful to the underlying distribution hinton2002training?, leroux2011learning?, tubiana2017emergence?. For example, Salakhutdinov2007restricted? trained RBMs to predict the interests of Netflix users by fitting a database of movie ratings, kiviken2012multiple? used them to generate realistic textures and srivastava2013modelling? employed them for topic modelling. There are even universal approximation theorems stating that an RBM with an arbitrary number of hidden units can approximate any distribution with binary support arbitrarily well leroux2008representational?, montufar2011refinements?. The current theoretical understanding of RBMs is largely based on statistical mechanics. For example, decelle2017spectral?, decelle2018thermodynamics? used statistical mechanics to model RBM training on real data. The statistical mechanics community also studied RBMs trained on simpler, synthetic datasets decelle2017spectral?, decelle2018thermodynamics?. In particular, many works investigated the teacher-student setting where a student RBM is trained with data produced by a teacher RBM huang2016unsupervised?, huang2017statistical?, barra2017phase?, barra2018phase?, huang2018role?, hou2019minimal?, decelle2021inverse?, manzan2025effect?. Such studies are crucial to isolate individual characteristics of structured datasets and Neural Network (NN) design choices in a controlled environment and explain their effects on NN training. For example, barra2017phase?, barra2018phase?, manzan2025effect? investigated the effects of the prior chosen for the data on RBM learning. In this paper, we study the effects of data correlations and number of hidden units on RBMs in the teacher-student setting. In particular, we show that RBMs with a few hidden units in the teacher-student setting can serve as a toy model of the lottery ticket hypothesis frankle2018lottery?, zhou2019deconstructing?, ramanujan2020hidden?, eran2020proving?, frankle2020linear?.

RBMs have a visible layer \(\sigma = \left\{ \sigma_i \right\}_{i = 1}^N\), a hidden layer \(\tau = \left\{ \tau_{\mu} \right\}_{\mu = 1}^P\) and a set of internal connections \(\xi = \left\{ \xi^{\mu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\), which are commonly referred to as weights. The visible layer is a set of concrete features found directly in the data, the hidden layer is an internal representation of the data in terms of abstract concepts, and the weights represent how the input features and abstract concepts are correlated with one another. To give a caricatural example, the presence of pointy ears as an input feature could be correlated with the abstract concept of a cat in the weights of a particular RBM. Given \(\xi\), the visible and hidden layer follow the joint distribution \[\begin{align} \label{eq:Gibbs} \mathop{\mathrm{P}}_\beta \left( \sigma, \tau \mid \xi \right) &= Z_\beta \left( \xi \right)^{-1} \mathop{\mathrm{P}}\left( \sigma \right) \mathop{\mathrm{P}}\left( \tau \right) \exp \left(-\beta H \left[ \sigma, \tau ; \xi \right] \right), \end{align}\tag{20}\] where the Hamiltonian \(H \left[ \sigma, \tau ; \xi \right] = -\frac{1}{\sqrt{N}} \sum_{\mu = 1}^P \tau_{\mu} \sum_{i = 1}^N \xi^{\mu}_i \sigma_i\) weighs the cost of every RBM configuration, \(\mathop{\mathrm{P}}\left( \sigma \right)\) and \(\mathop{\mathrm{P}}\left( \tau \right)\) are priors on the visible and hidden layers, respectively, and \(Z_\beta \left( \xi \right) = \mathbb{E}_{\sigma, \tau} \left[ \exp \left( -\beta H \left[ \sigma, \tau ; \xi \right] \right) \right]\) is the partition function normalizing the distribution, with \(\mathbb{E}_{\sigma, \tau}\) the joint expectation over the visible and hidden unit priors. Intuitively, the priors are the default distributions of \(\sigma\) and \(\tau\) when the Hamiltonian does not contribute to \(\mathop{\mathrm{P}}_\beta \left( \sigma, \tau \mid \xi \right)\), i.e. when \(\beta\) is zero. Gibbs ditributions of the form (20 ) have been deeply investigated in the mathematical physics community for their link with the Parisi theory of spin glasses barra2015multi?, panchenko2015free?, barra2012glassy?, genovese2020legendre?, genovese2017overlap?, genovese2016non?. An RBM with a given \(\xi\) can generate data \(\sigma\) by sampling the marginal distribution \[\label{eq:RBM95direct95distribution} \mathop{\mathrm{P}}_\beta \left( \sigma \mid \xi \right) = Z_\beta \left( \xi \right)^{-1} \mathop{\mathrm{P}}\left( \sigma \right) \psi_\beta \left( \sigma ; \xi \right)= Z_\beta \left( \xi \right)^{-1} \mathop{\mathrm{P}}(\sigma) \prod_{\mu=1}^P \phi_\beta(\xi^\mu\cdot \sigma),\tag{21}\] where \(\phi_\beta \left(\xi^\mu\cdot \sigma \right) = \mathbb{E}_{\tau_\mu} \left[ \exp \left( \frac{\beta}{\sqrt{N}} \tau_{\mu} \sum_{i = 1}^N \xi^{\mu}_i \sigma_i \right) \right]\), \(\psi_\beta \left( \sigma ; \xi \right)\) factorizes as \(\psi_\beta \left( \sigma ; \xi \right) = \prod_{\mu=1}^P \phi_\beta(\xi^\mu\cdot \sigma)\) and \(Z_\beta \left( \xi \right) = \mathbb{E}_{\sigma} \left[ \psi_\beta \left( \sigma ; \xi \right) \right]\). Marginal Gibbs distributions of the form (21 ) are also known as generalized Hopfield networks barra2017phase?, barra2018phase?, hopfield1982neural?, sollich2014extensive?, agliari2018non?, agliari2015retrieval?, agliari2015hierarchical?. Conversely, following Bayes’ theorem, the weights of an RBM can be trained on a dataset \(\boldsymbol{\sigma} = \left\{ \sigma^a_i \right\}_{1 \leq i \leq N}^{1 \leq a \leq M}\) of \(M\) examples by sampling the posterior distribution \[\label{eq:inverse95distribution} \mathop{\mathrm{P}}_{\beta} \left( \xi | \boldsymbol{\sigma} \right) = \mathcal{Z}_{\beta} \left( \boldsymbol{\sigma} \right)^{-1} \mathop{\mathrm{P}}\left( \xi \right) \prod_{a = 1}^M \mathop{\mathrm{P}}_\beta \left( \sigma^a | \xi \right),\tag{22}\] where \(\mathop{\mathrm{P}}\left( \xi \right)\) is a prior on the weights and \(\mathcal{Z}_{\beta} \left( \boldsymbol{\sigma} \right) = \mathbb{E}_\xi \left[ \prod_{a = 1}^M \mathop{\mathrm{P}}_\beta \left( \sigma^a | \xi \right) \right]\) is the posterior partition function normalizing the distribution. In the teacher-student setting, the data \(\boldsymbol{\sigma}\) used to train the RBM is produced by another RBM huang2016unsupervised?, huang2017statistical?, barra2017phase?, barra2018phase?, huang2018role?, hou2019minimal?, decelle2021inverse?, manzan2025effect?. In other words, a student RBM is trained using a dataset supplied by a teacher RBM. The student’s ability to fit the teacher’s data can then be evaluated in terms of the so-called overlaps \(Q(\xi^{*\mu},\xi^\nu)=\frac{1}{N} \sum_{i = 1}^N \xi^{* \mu}_i \xi^{\nu}_i\) between the teacher’s weights \(\xi^* = \left\{ \xi^{* \mu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P^*}\) and the student’s weights \(\xi = \left\{ \xi^{\nu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\), whose rows \(\xi^{* \mu} = \left\{ \xi^{* \mu}_i \right\}_{i = 1}^N\) and \(\xi^{\mu} = \left\{ \xi^{\nu}_i \right\}_{i = 1}^N\) are also called patterns. As such, we refer to the expected value of \(Q(\xi^{*\mu},\xi^\nu)\) as the student’s performance. This framework was used to predict the critical amount of data needed to train multiple variants of RBMs huang2016unsupervised?, huang2017statistical?, barra2017phase?, barra2018phase?, huang2018role?, hou2019minimal?, manzan2025effect?.

In barra2017phase?, it was observed empirically and conjectured that the performance of a student RBM is independent of the number of hidden units when the teacher patterns are uncorrelated. It was later shown analytically that a student RBM taught by a teacher with uncorrelated patterns achieves the same performance whether the teacher and the student have \(1\) or \(2\) hidden units each hou2019minimal?. In a nutshell, this simplification occurs because an RBM with \(2\) uncorrelated hidden units effectively factorizes into \(2\) RBMs with \(1\) hidden unit each. The conjecture of barra2017phase? was explicitly rephrased in terms of this factorization property in hou2019minimal?². Beyond this idealized scenario, it has long been known that machine learning models benefit from the correlations found in structured data to build an efficient internal representation gardner1988space?. This phenomenon has received considerable attention in recent theoretical studies decelle2017spectral?, decelle2018thermodynamics?, dobriban2018high?, goldt2020modelling?, wu2020optimal?, ghorbani2021neural?, liao2021random?, chen2021multiple?, nakkiran2021optimal?, richards2021asymptotics?, ascoli2021interplay?, loureiro2021learning?, refinetti2021classifying?, ichikawa20222statistical?, notably for RBMs with \(2\) hidden units in the teacher-student setting hou2019minimal?. Extending the latter study to an arbitrary number of hidden units remains an intriguing open problem. In fact, although increasing the number of hidden units in the teacher model may bring us closer to the complexity of real data—complexity that RBMs are known to capture due to universal approximation theorems leroux2008representational?, montufar2011refinements?—there is still no theoretical framework that fully explains the learning performance of RBMs for an arbitrary number of hidden units.

In this work, we evaluate the student’s learning performance in the teacher-student setting where both the teacher and the student have an arbitrary finite number of hidden units and the teacher patterns are allowed to be correlated with one another. In particular, we evaluate the critical data load \(\alpha_{\text{crit}} = \frac{M}{N}\) above which learning becomes possible. In Section 3.2, we introduce the teacher-student setting and the replica method used in our calculations. In Section 3.3.1, we present the so-called saddle-point equations governing the performance and the critical load that we obtain from it. In Section 3.3.2, we discuss the case where there are no correlations between the teacher patterns. This Section is divided into three Subsections. In Subsection 3.3.2.1 and 3.3.2.2, we show that the student’s performance is independent of the number of hidden units when the teacher patterns are uncorrelated. In Subsection 3.3.2.3, we argue that our teacher-student setting with uncorrelated teacher patterns can serve as a toy model of the lottery ticket hypothesis. Next, in Section 3.3.3, we discuss the effects of uniform teacher pattern correlations on the student’s performance. In particular, in Subsection 3.3.3.1, we discuss the different learning phases of the teacher-student problem as a function of the correlations and the number of hidden units. Finally, in Section 3.3.4, we discuss random correlations and compare their effect on the performance to that of uniform correlations. Throughout the paper, we compare key results against Monte Carlo simulations. The code and hyperparameter values of the training algorithms used to make the figures are available at the following public Github repository theriault2025modellingsoftware?.

3.2 Model↩︎

In the teacher-student setting, a student RBM with marginal likelihood \(\mathop{\mathrm{P}}_\beta \left( \sigma | \xi \right)\) (Eq. 21 ) is trained using a dataset \(\boldsymbol{\sigma} = \left\{ \sigma^a_i \right\}_{1 \leq i \leq N}^{1 \leq a \leq M}\) of \(M\) examples \(\sigma^a = \left\{ \sigma^a_i \right\}_{i = 1}^{N}\) of dimension \(N\), generated by a teacher RBM with a prescribed marginal likelihood \(\mathop{\mathrm{P}}_{\beta^*} \left( \sigma | \xi^* \right)\). We call \(\alpha=\frac{M}{N}\) the ratio between the size of the training set and the input dimension. In this scenario, the student knows that the correct model for the data is an RBM. However, it does not necessarily know the number of hidden units \(P^*\) and the inverse temperature \(\beta^*\) used by the teacher. Therefore, unless explicitly stated otherwise, we will assume that the number of hidden units \(P\) and the inverse inference temperature \(\beta\) of the student are not necessarily the same as the teacher’s. For convenience, we will frequently state our results in terms of the temperatures \(T^* = 1/\beta^*\) and \(T = 1/\beta\) rather than in terms of \(\beta^*\) and \(\beta\).

For simplicity, we assume that the visible and hidden units of both the student and the teacher take values in \(\left\{ -1, +1 \right\}\) with a uniform prior, i.e. they are binary random variables with no prior bias towards \(-1\) or \(+1\). We impose structure in the data by taking the teacher patterns to be random variables with a fixed covariance matrix \(\mathcal{Q}\). To be more precise, we assume that the columns \(\xi^*_i = \left\{ \xi^{* \mu}_i \right\}_{\mu = 1}^{P^*}\) of \(\xi^*\) are i.i.d. random variables, with mean \(0\) and a well-defined \(P^* \times P^*\) covariance matrix \(\mathcal{Q}\). Uncorrelated teacher patterns are obtained by setting \(\mathcal{Q} = \mathbf{I}\).

As previously mentioned, the student learns its patterns \(\xi\) by sampling them from the posterior distribution \[\begin{align} \label{eq:RBM95posterior} \mathop{\mathrm{P}}_{\beta} \left( \xi | \boldsymbol{\sigma} \right) &=& \mathcal{Z}_{ \beta} \left( \boldsymbol{\sigma} \right)^{-1} \mathop{\mathrm{P}}\left( \xi \right) \prod_{a = 1}^M \mathop{\mathrm{P}}_\beta \left( \sigma^a | \xi \right) = \mathcal{Z}_{ \beta} \left( \boldsymbol{\sigma} \right)^{-1} \mathop{\mathrm{P}}\left( \xi \right) \prod_{a = 1}^M \left[ Z_{\beta} \left( \xi \right)^{-1} \psi_\beta \left( \sigma^a ; \xi \right) \right]\nonumber \\ &=& \mathcal{Z}_{ \beta} \left( \boldsymbol{\sigma} \right)^{-1}\mathop{\mathrm{P}}\left( \xi \right) \prod_{a = 1}^M \left[ Z_{\beta} \left( \xi \right)^{-1} \mathop{\mathrm{P}}(\sigma^a)\prod_{\mu=1}^P\phi_\beta \left( \sigma^a \cdot \xi^\mu \right) \right]\nonumber \\ &=& \mathcal{Z}_{ \beta} \left( \boldsymbol{\sigma} \right)^{-1}\mathop{\mathrm{P}}\left( \sigma \right)Z_{\beta} \left( \xi \right)^{-M}\prod_{\mu=1}^P \psi_\beta(\xi^\mu;\sigma). \end{align}\tag{23}\] Compared with Eq. (21 ), the posterior distribution is composed of \(P\) generalized Hopfield distributions, one for each hidden unit. The data now plays the role of dual patterns, which interact only through the term \(Z_{\beta} \left( \xi \right)^{-M}\). The latter is therefore responsible for encoding mutual correlation between student patterns, otherwise independent barra2017phase?, decelle2021inverse?, alemanno2023hopfield?. The student has no access to the structure of the teacher patterns, so we assume an identity covariance matrix for the student pattern prior \(\mathop{\mathrm{P}}(\xi)\).

In general, the partition functions \(Z_{\beta} \left( \xi \right)\) and \(\mathcal{Z}_{ \beta} \left( \boldsymbol{\sigma} \right)\) are intractable, which makes the properties of Eq. (23 ) difficult to study. However, our assumptions about the priors on the visible units, hidden units and patterns make analytical computations possible in the limit of \(N, M \rightarrow \infty\).

The overlaps \(Q(\xi^{*\mu},\xi^\nu)=\frac{1}{N} \sum_{i = 1}^N \xi^{* \mu}_i \xi^{\nu}_i\) are a good measure of the student’s learning performance because they quantify how close the student patterns are to the teacher patterns. We exploit techniques from statistical mechanics to compute the expected value \(\mathbb{E}_{\xi^*,\boldsymbol{\sigma},\xi} \left[Q(\xi^{*\mu},\xi^\nu)\right]\) of the overlaps with respect to the distribution of the teacher patterns, the generated dataset and the inferred student patterns. Specifically, in the limit of large dataset and data dimension \(M,N\to \infty\), we obtain the expected value of the overlaps as a byproduct of the limiting quenched free entropy \[\label{eq:def95free95entropy} f (\alpha, \beta^*, \beta, P^*, P, \mathcal{Q}) = \lim_{M,N\to\infty} \frac{1}{N} \mathbb{E}_{\xi^*,\boldsymbol{\sigma}} \log \left[ \mathcal{Z}_{ \beta}(\boldsymbol{\sigma}) \right],\tag{24}\] where \(\mathbb{E}_{\xi^*,\boldsymbol{\sigma}}\) is the expected value w.r.t. the joint distribution of the teacher patterns and generated dataset. In fact, we show in Section 3.3.1 that Eq. (24 ) can be expressed as the result of a variational principle w.r.t. a set of order parameters, whose optimum gives the expected value of the overlaps. We then use this result to investigate the effects of \(P\) and \(\beta\) on the student’s ability to learn a dataset characterized by \(P^*, \beta^*\) and \(\mathcal{Q}\), as well as the impact of these hyperparameters on the critical threshold \(\alpha_{\text{crit}}\) beyond which learning becomes possible. We focus our quantitative discussion on Gaussian and binary priors \(\mathop{\mathrm{P}}(\xi), \mathop{\mathrm{P}}(\xi^*)\) (see 3.5.1), but our results can be generalized to other pattern distributions that meet the other assumptions stated earlier in this Section. The Gaussian case is particularly interesting because it is closely related to RBMs used in practical applications kiviken2012multiple?, Salakhutdinov2007restricted?, srivastava2013modelling?, which usually have continuous weights.

3.3 Results and discussion↩︎

3.3.1 Free entropy and saddle point equations↩︎

The free entropy \(f\) can be computed by exploiting a well-established statistical mechanics technique called the replica trick, which is based on the identity \[\begin{align} f = \lim_{N\to\infty}\frac{1}{N} \mathbb{E}\log \left[ \mathcal{Z} \right] &= \lim_{N\to\infty, L \rightarrow 0} \left( \frac{1}{LN} \log \mathbb{E} \left[\mathcal{Z}^L \right] \right). \end{align}\] Calculations are shown in the Appendices. In 3.6, we calculate the average replicated partition function \(\mathbb{E}_{\xi^*,\boldsymbol{\sigma}}\left[ \mathcal{Z}^L \right]\) in the limit of \(N \rightarrow \infty\). In 3.7, we use \(\mathbb{E}_{\xi^*,\boldsymbol{\sigma}}\left[ \mathcal{Z}^L \right]\) to evaluate the quenched free entropy \(f\) under the so-called Replica-Symmetric (RS) approximation. We find that the RS free entropy can be expressed in terms of the variational principle \[\begin{align} \label{eq:RBM95free95entropy} f &= \mathop{\mathrm{Extr}}_{m, \hat{m}, q, \hat{q}, s, \hat{s}} \Bigg\{ -\sum_{\mu = 1}^{P^*} \sum_{\nu = 1}^{P} m^{\mu \nu} \hat{m}^{\mu \nu} - \frac{1}{2} \sum_{\mu \neq \nu}^{P} s^{\mu \nu} \hat{s}^{\mu \nu} + \frac{1}{2} \sum_{\mu, \nu = 1}^P q^{\mu \nu} \hat{q}^{\mu \nu} \\ &\quad+ \mathbb{E}_{\xi^*} \mathbb{E}_{z} \log \left[ \mathcal{Z} \left( \mathcal{L}^C \right) \right] + \alpha \left\langle \mathbb{E}_z \log \left[ \mathcal{Z} \left( \mathcal{L}^O \right) \right]\right\rangle_{\mathcal{M}_*} - \alpha \log \left[ \mathcal{Z} \left( \mathcal{M} \right) \right] \Bigg\}, \addtocounter{equation}{1}\end{align}\tag{25}\] where \(z = \left[ z_{\mu \nu} \right]_{\mu,\nu = 1}^{P}\) is a set of i.i.d. standard Gaussian random variables, and where the effective energies \(\mathcal{M}_*\), \(\mathcal{M}\), \(\mathcal{L}^O\) and \(\mathcal{L}^C\) are defined as \[\begin{align} \mathcal{M}_* \left( \tau_* \right) &=& \frac{1}{2} \left[ \beta^* \right]^2 \sum_{\mu, \nu = 1}^{P^*} \mathcal{Q}_{\mu \nu} \tau^*_{\mu} \tau^*_{\nu} \tag{26} \\ \mathcal{M} \left( \tau \right) &=& \frac{1}{2} \beta^2 \sum_{\mu, \nu = 1}^{P} s^{\mu \nu} \tau_{\mu} \tau_{\nu} \tag{27} \\ \mathcal{L}^O \left( \tau; \tau^*, z \right) &=& \mathcal{L}_{\beta^*, \beta} \left( \tau, \tau^*, z ; m, s, q \right)\tag{28} \\ \mathcal{L}^C \left( \xi; \xi^*, z \right) &=& \mathcal{L}_{1, 1} \left( \xi, \xi^*, z ; \hat{m}, \hat{s}, \hat{q} \right) \tag{29} \\ \text{with} \quad \mathcal{L}_{\lambda_1, \lambda_2} \left( \xi, \xi^*, z ; m, s, q \right) &=& \frac{1}{2} \left[ \lambda_2 \right]^2 \sum_{\mu, \nu = 1}^P \left( s^{\mu \nu} - q^{\mu \nu} \right) \xi^{\mu} \xi^{\nu} \tag{30} \\ &\quad&+ \lambda_1 \lambda_2 \sum_{\mu = 1}^{P^*} \sum_{\nu = 1}^{P} m^{\mu \nu} \xi^{* \mu} \xi^{\nu} \\ &\quad&+ \lambda_2 \sum_{\mu, \nu = 1}^P A_{\mu \nu} \left( q \right) z_{\mu \nu} \frac{\xi^{\mu} + \xi^{\nu}}{2} \nonumber \\ \text{and} \quad A_{\mu \nu} \left( q \right) &=& \sqrt{2 q^{\mu \nu} - \delta_{\mu \nu} \sum_{\eta = 1}^P q^{\mu \eta}} \tag{31}. \end{align}\] In particular, we use the notation \(\mathcal{Z}(f) = \mathbb{E}_x \left[ \exp \left\{ f \left( x ; \cdot \right) \right\} \right]\) for the partition function and \(\left\langle g \right\rangle_f=\mathbb{E}_x \left\{ g(x) \mathcal{P} \left[ f \right] \left( x ; \cdot \right) \right\}\) for the expectation value of an observable \(g\) w.r.t. the Gibbs distribution \(\mathcal{P} \left[ f \right] \left( x ; \cdot \right) = \mathcal{Z}(f)^{-1} \exp \left[ f \left( x ; \cdot \right) \right]\) of an effective energy \(f\), where \(\mathbb{E}_x\) is the expectation over the prior of \(x\). As explained in 3.5.3, \(q^{\mu \nu}\) is assumed symmetric in Eq. (30 ). The solution of this variational principle (see 3.8 for a detailed derivation) must obey the saddle-point equations \[\begin{align} \label{eq:RBM95saddle-point} m^{\mu \nu} &= \mathbb{E}_{\xi^*} \mathbb{E}_{z} \left[ \xi^{* \mu} \left\langle \xi^{\nu} \right\rangle_{\mathcal{L}^C} \right] \\ s^{\mu \nu} &= \mathbb{E}_{\xi^*} \mathbb{E}_{z} \left[ \left\langle \xi^{\mu} \xi^{\nu} \right\rangle_{\mathcal{L}^C} \right] \\ q^{\mu \nu} &= \mathbb{E}_{\xi^*} \mathbb{E}_{z} \left[ \left\langle \xi^{\mu} \right\rangle_{\mathcal{L}^C} \left\langle \xi^{\nu} \right\rangle_{\mathcal{L}^C} \right],\\ \hat{m}^{\mu \nu} &= \beta^* \beta \alpha \;\left\langle \mathbb{E}_{z} \left[ \tau^*_{ \mu} \left\langle \tau_{\nu} \right\rangle_{\mathcal{L}^O} \right] \right\rangle_{\mathcal{M}_*} \addtocounter{equation}{1}\\ \hat{s}^{\mu \nu} &= \beta^2 \alpha \left( \left\langle \mathbb{E}_{z} \left[ \left\langle \tau_{\mu} \tau_{\nu} \right\rangle_{\mathcal{L}^O} \right]\right\rangle_{\mathcal{M}_*} - \left\langle \tau_{\mu} \tau_{\nu} \right\rangle_{\mathcal{M}} \right) \\ \hat{q}^{\mu \nu} &= \beta^2 \alpha \left\langle \mathbb{E}_{z} \left[ \left\langle \tau_{\mu} \right\rangle_{\mathcal{L}^O} \left\langle \tau_{\nu} \right\rangle_{\mathcal{L}^O} \right] \right\rangle_{\mathcal{M}_*}. \end{align}\tag{32}\] The optimal order parameters \(m\in \mathbb{R}^{P^*\times P}\), \(s\in \mathbb{R}^{P\times P}\) and \(q\in\mathbb{R}^{P\times P}\) solving Eqs. (32 ) are to be interpreted as expected overlaps. First of all, \(m^{\mu \nu}\), which is commonly known as the Mattis magnetization, is the limiting expected overlap between the teacher pattern \(\xi^{* \mu}\) and the student pattern \(\xi^{\nu}\), i.e. \[m^{\mu \nu} = \lim_{N,M \rightarrow \infty} \mathbb{E}_{\xi^*,\boldsymbol{\sigma},\xi} \left[Q(\xi^{*\mu},\xi^\nu)\right]\;\;\;\;\mu=1,\ldots, P^* \; \nu=1,\ldots, P,\] where \(\mathbb{E}_{\xi^*,\boldsymbol{\sigma},\xi}\) is the expectation w.r.t. the joint distribution \(\mathop{\mathrm{P}}(\xi^*)\mathop{\mathrm{P}}_{\beta^*}(\boldsymbol{\sigma}|\xi^*)\mathop{\mathrm{P}}_\beta(\xi|\boldsymbol{\sigma})\) of teacher patterns, dataset and student patterns. Second of all, \(s^{\mu \nu}\) with \(\mu \neq \nu\) is the limiting expected overlap between any two student patterns \(\xi^{\mu}\) and \(\xi^{\nu}\) from the same sample \(\xi \sim \mathop{\mathrm{P}}_\beta \left( \xi | \boldsymbol{\sigma} \right)\), i.e. \[s^{\mu \nu} = \lim_{N,M \rightarrow \infty} \mathbb{E}_{\xi^*,\boldsymbol{\sigma},\xi} \left[Q(\xi^\mu,\xi^\nu)\right]\;\;\;\mu,\nu=1,\ldots,P.\] Finally, \(q^{\mu \nu}\) is the limiting expected overlap between any two student patterns \(\xi^{1 \mu}\) and \(\xi^{2 \nu}\) from two independent posterior samples \(\xi^1\) and \(\xi^2\), i.e. \[\begin{align} q^{\mu \nu} &= \lim_{N,M \rightarrow \infty} \mathbb{E}_{\xi^*,\boldsymbol{\sigma}, \xi^1\times\xi^2} \left[Q(\xi^{1 \mu},\xi^{2 \nu})\right] \\ &= \lim_{N,M \rightarrow \infty} \mathbb{E}_{\xi^*,\boldsymbol{\sigma}} \left[Q(\mathbb{E}_{\xi|\boldsymbol{\sigma}}\left[\xi^{\mu}\right],\mathbb{E}_{\xi|\boldsymbol{\sigma}}\left[\xi^{\nu}\right])\right]\;\;\;\mu,\nu=1,\ldots,P, \addtocounter{equation}{1}\end{align}\] where \(\mathbb{E}_{\xi|\boldsymbol{\sigma}}\) indicates the expectation w.r.t. the posterior distribution (Eq. 23 ). While \(s\) measures the effective correlation between student patterns, the so-called spin-glass order parameter \(q\) quantifies the tendency of the student to stay frozen in specific configurations rather than sampling all possible patterns uniformly. For simplicity’s sake, we will usually call \(m\) the magnetization and \(q\) the Spin-Glass (SG) overlap. As in the Introduction, we will also occasionally refer to \(m\) as the student’s performance.

The RS saddle point equations (Eqs. 32 ) can be solved by numerical iteration for any values of the hyperparameters \(\beta^*\), \(\beta\), \(\alpha\), \(P^*\) and \(P\) (see 3.13). We expect the RS solution to be exact when the student is fully informed about the teacher’s prior and hyperparameters and matches them with its own, i.e. \(\beta = \beta^*\), \(P = P^*\) and \(\mathop{\mathrm{P}}\left( \xi \right) = \mathop{\mathrm{P}}\left( \xi^* \right)\). This regime of complete information is commonly referred to as the Nishimori line nishimori1980exact?, iba1999nishimori?, nishimori2001statistical?, contucci2009spin?. When \(\beta = \beta^*\), we find two different phases (see Figs. 20, 21, 28 and 34) :

• the Paramagnetic (P) phase, where the order parameters \(m\), \(s\) and \(q\) all vanish;

• the Ferromagnetic (F) phase, where they are all larger than zero.

Intuitively, the student RBM can partially learn the teacher patterns in the F phase but becomes unable to do so in the P phase. The corresponding P-F phase transition is thus also the onset of learning. When \(\beta \neq \beta^*\), we also find

• the Spin-Glass (SG) phase, where \(m= 0\) but \(q>0\).

In this phase, the student converges to spurious low-energy states (\(q > 0\)) unrelated to the teacher patterns (\(m = 0\)) (see Figs. 29 and 30). Looking at Figs. (20), (21), (28) and (34), the P-F phase transition appears to be second order. In general, such second-order phase transitions coincide with the onset of instability of the paramagnetic solution in the saddle-point equations (Eqs. 32 ). The instability condition (see 3.10) reads as \[\begin{align} \label{eq:critical95load} \alpha \geq \alpha_{\text{crit}} = \frac{1}{\left[ \beta^* \beta \right]^2 \lambda^{\mathcal{S}}_{\text{max}}}, \end{align}\tag{33}\] where \(\lambda^{\mathcal{S}}_{\text{max}}\) is the largest eigenvalue of the matrix \(\mathcal{S} = \mathcal{Q} \mathcal{R}\), with the covariance matrix \(\mathcal{R}\) defined as \(\mathcal{R}_{\mu\nu}=\left\langle \tau^{*\mu}\tau^{*\nu}\right\rangle_{\mathcal{M}_*}\). As expected, the \(\alpha_{\text{crit}}\) of Eq. (33 ) coincides with the P-F phase transition of Figs. (28) and (34). Interestingly, \(\alpha_{\text{crit}}\) does not depend on the number of hidden units \(P\) of the student nor on its prior \(\mathop{\mathrm{P}}\left( \xi \right)\). Outside the Nishimori line, i.e. when \((\beta,P,P(\xi)) \neq (\beta^*,P^*,P(\xi^*))\), the RS approximation is not expected to be exact. Therefore, in principle, one would need to calculate the so-called Replica Symmetry Breaking (RSB) corrections of \(\alpha_{\text{crit}}\). However, we find that Eq. (33 ) is consistent with Monte Carlo simulations even when \(P\neq P^*\) (see Fig. 23) and \(\mathop{\mathrm{P}}(\xi)\neq \mathop{\mathrm{P}}(\xi^*)\) (see Fig. 32). As such, the RS approximation of \(\alpha_{\text{crit}}\) seems robust to the priors and number of hidden units in practice.

3.3.2 Learning uncorrelated patterns↩︎

In this Section, including Subsections 3.3.2.1, 3.3.2.2 and 3.3.2.3, we take \(\mathcal{Q} = \mathbf{I}\), i.e. there are no correlations between the teacher patterns. In that regime, we find that the teacher-student setting exhibits an effective factorization property on the student patterns, which makes the behavior of the student RBM explainable in terms of that with a single hidden unit. A first indication of this result is given by the critical load \(\alpha_{\text{crit}}\) (see Eq. 33 ). As we explained in the previous Section, no matter \(\mathcal{Q}\), Eq. (33 ) does not depend on the number \(P\) of student patterns. When \(\mathcal{Q} = \mathbf{I}\), and thus \(\mathcal{S} = \mathbf{I}\) (see Eqs. 33 ), we find that \[\begin{align} \label{eq:critical95load95without95correlations} \alpha_{\text{crit}} = \frac{1}{\left[ \beta^* \beta \right]^2}, \end{align}\tag{34}\] which does not depend on the number \(P^*\) of teacher patterns either. Eq. (34 ) generalizes to arbitrary finite \(P^*\) and \(P\) the critical load \(\alpha_{\text{crit}} = \beta^{-4}\) previously found for \(P^* = P = 1\) hidden units huang2016unsupervised?, huang2017statistical? and \(P^* = P = 2\) uncorrelated hidden units hou2019minimal? when \(\beta^* = \beta\). As we show in the next Subsections, the previous universality result also holds for the order parameters, which we recall are a good measure of the student’s learning performance.

3.3.2.1 Independence of the number of hidden units: binary patterns

In the case of binary patterns, Eqs. (32 ) with \(P=P^*=1\) simplify to \[\begin{align} \label{eq:reduced95saddle-point95P1} m &= \mathbb{E}_{z} \left[ \tanh \left( \hat{m} + \sqrt{\hat{q}} z \right) \right] \\ q &= \mathbb{E}_{z} \left[ \tanh^2 \left( \hat{m} + \sqrt{\hat{q}} z \right) \right] \addtocounter{equation}{1}\\ \hat{m} &= \beta^* \beta \alpha \;\mathbb{E}_{z} \left[ \tanh \left( \beta^* \beta m + \beta \sqrt{q} z \right) \right] \\ \hat{q} &= \beta^2 \alpha \;\mathbb{E}_{z} \left[ \tanh^2 \left( \beta^2 m + \beta \sqrt{q} z \right) \right], \end{align}\tag{35}\] whose solution is summarized in the phase diagrams of Fig. (19). To be more precise, the left and right panels of Fig. (19) show the fixed-\(\beta^*\) regime and the Nishimori regime \(\beta=\beta^*\), respectively. For clarity’s sake, we will discuss Fig. (19) and the solution of Eq. (35 ) in terms of \(T^* = 1/\beta^*\) and \(T = 1/\beta\) for the remainder of this paragraph. We recall that the student can learn the teacher patterns in the Ferromagnetic (F) phase (\(m \neq 0\)), but not in the Paramagnetic (P) phase (\(m = q = 0\)) nor in the Spin-Glass (SG) phase (\(m = 0\) and \(q > 0\)). For \(T < T^*\), we find the SG phase between the P phase and the F phase. The corresponding P-SG transition line \(\alpha_{\text{P-SG}} = T^4\) extends until the Nishimori line \(T = T^*\), where it meets the F phase. Therefore, the Nishimori line crosses the triple point of the P, SG and F phases, as expected from spin glass theory nishimori1980exact?, nishimori2001statistical?. On the Nishimori line, it is straightforward to verify that \(m = q\), which is also expected from spin glass theory nishimori1980exact?, nishimori2001statistical?. In particular, \(m\) and \(q\) simultaneously become nonzero above the critical load \(\alpha_{\text{crit}} = \left[ T^* T \right]^2 = T^4\) of the P-F phase transition on the Nishimori line, which prevents the SG phase from forming (see Fig. 19, right panel). As in some related inference problems in the teacher-student setting, \(\alpha_{\text{P-SG}} = T^4\) is identical to the \(\alpha_{\text{crit}}\) found on the Nishimori line manzan2025effect?, alemanno2023hopfield?, theriault2024dense?. For \(T > T^*\), the SG phase does not exist either, and the P phase transitions directly to the F phase at \(\alpha_{\text{crit}} = \left[ T^* T \right]^2\).

Looking at the SG-F transition, we see that decreasing the inference temperature \(T\) too much can make it harder for the student to learn the teacher patterns. In particular, a student with enough data to learn the teacher patterns at a given inference temperature \(T_1\) may fail to do so at a lower inference temperature \(T_2 < T_1\). The RS approximation is probably not completely accurate in this regime, so one would need to calculate RSB corrections to study this behavior quantitatively.

When \(P, P^*\) are arbitrary finite numbers and \(\mathcal{Q} = \mathbf{I}\), we search for solutions of Eqs. (32 ) by making the ansatz \[\label{eq:no95correlations95ansatz} \begin{align} &m^{\mu \nu} = \delta_{\mu \nu} m, &\quad &\hat{m}^{\mu \nu} = \delta_{\mu \nu} \hat{m}, \\ &s^{\mu \nu} = \delta_{\mu \nu}, &\quad &\hat{s}^{\mu \nu} = 0, \\ &q^{\mu \nu} = \begin{cases} \delta_{\mu \nu} q &\quad \mu, \nu \leq P^* \\ \delta_{\mu \nu} g &\quad \text{otherwise}, \end{cases} &\quad &\hat{q}^{\mu \nu} = \begin{cases} \delta_{\mu \nu} \hat{q} &\quad \mu, \nu \leq P^* \\ \delta_{\mu \nu} \hat{g} &\quad \text{otherwise}, \end{cases} \end{align}\tag{36}\] which describes the situation where the student learns the teacher patterns one-to-one. Following the nomenclature introduced in hou2019minimal?, we will call Eq. (36 ) the Permutation Symmetry Breaking (PSB) ansatz. We have assumed without loss of generality that the \(\min \left\{ P, P^* \right\}\) nonzero magnetizations are on the main diagonal, as any hidden unit permutation would give an equivalent solution. According to this ansatz, the first \(\min \left\{ P, P^* \right\}\) student patterns converge one-to-one to the first \(\min \left\{ P, P^* \right\}\) teacher patterns with magnetization \(m\) and SG overlap \(q\). When \(P>P^*\), the remaining \(P - P^*\) student patterns (i.e. \(P\geq\nu>P^*\)) are aligned in spurious directions (i.e. \(m^{\mu\nu}=0\; \forall \mu=1,...,P^*\)) with a different SG overlap \(g \neq q\). Conversely, when \(P<P^*\), \(P^*-P\) of the teacher patterns are not learned for lack of student hidden units. The latter case is described in terms of only \(m\) and \(q\), and so is \(P = P^*\). Under the PSB ansatz, Eqs. (32 ) decouple into \(P\) independent systems of equations for the \(P\) hidden units \(\nu = 1, ..., P\) (see 3.11). In other words, the student factorizes into \(P\) students with one hidden unit each. The systems of equations for the first \(\min \left\{ P, P^* \right\}\) patterns are all identical to each other and equivalent to Eqs. (35 ). When \(P>P^*\), the \(P - P^*\) remaining systems of equations all take the form \[\begin{align} \label{eq:reduced95saddle-point} g &= \mathbb{E}_{z} \left[ \tanh^2 \left( \sqrt{\hat{g}} z \right) \right]\\ \hat{g} &= \beta^2 \alpha \;\mathbb{E}_{z} \left[ \tanh^2 \left( \beta \sqrt{g} z \right) \right].\nonumber \end{align}\tag{37}\] In summary, Eqs. (32 ) with binary patterns and \(\mathcal{Q} = \mathbf{I}\) reduce to Eqs. (35 , 37 ) under the PSB ansatz. In particular, the solution of Eqs. (35 , 37 ), which we will call the PSB solution, has \(m = q\) even in the region outside the Nishimori line where \(\beta = \beta^*\) and \(P \neq P^*\). In Fig. (20), we verify that the PSB solution can be found by iterating the original saddle-point equations (Eqs. 32 ) with binary student patterns (see 3.5.1), \(\mathcal{Q} = \mathbf{I}\) and near-diagonal initial conditions \(m^{\mu \nu} = \delta_{\mu \nu} m_0 + \left( 1 - \delta_{\mu \nu} \right) \varepsilon\), where \(\varepsilon \ll m_0\). In other words, we verify that the solution of Eqs. (35 , 37 ) is stable for the full saddle-point iteration. We plot only the case of \(P \geq P^*\) because \(P < P^*\) yields similar results.

The original saddle-point equations (Eq. 32 ) with binary patterns, \(\mathcal{Q} = \mathbf{I}\) and \(P > P^*\) also have other solutions that cannot be expressed as cleanly by a set of simplified saddle-point equations (see Fig. 21). These solutions can be found by initializing \(m^{ \mu\nu}\) with \(m_0 \gg \varepsilon\) on the diagonal and at least one off-diagonal entry close to \(m_0\), a process we will refer to as off-diagonal initialization. They correspond to a distributed representation where some of the student patterns learn the same teacher pattern, so we dub these solutions partial PSB. Throughout this paper and its figures, we focus on the case where at most two student patterns learn the same teacher pattern, but larger numbers are also possible. Within a partial PSB solution, some student patterns may still learn teacher patterns one-to-one. These student patterns \(\xi^{\mu}_{\text{PSB}}\) have the same \(m\) and \(q\) as Eq. (35 ), and in particular satisfy \(m = q\) (see Fig. 21). Comparison of the RS free entropy (Eq. 25 ) of the PSB and partial PSB solutions of Eqs. (32 ) suggests that the former is always favored in the limit \(N,M\to\infty\). The free entropy difference between them decreases with increasing \(T = T^*\) and vanishes at the onset of the paramagnetic phase (see Fig. 35). These results are confirmed by numerical simulations. In low \(T = T^*\) simulations, students with binary patterns always converge to the PSB solution, even when we use random initial conditions (see Figs. 23), which is solid evidence that PSB is favored at low \(T = T^*\). When \(T = T^*\) is relatively close to the paramagnetic phase, the simulations have relatively large error bars and repeatedly jump between the different solutions rather than converging to a single mode. This kind of instability can be attributed to finite-size fluctuations.

The PSB solution is independent of \(P^*\) and \(P\), which means that, in terms of the student’s performance \(m\), learning \(P^*\) i.i.d. patterns stored in a single teacher RBM is as difficult as learning \(P^*\) i.i.d. patterns from \(P^*\) separate RBMs with one hidden unit each. Such independence of the student’s performance from the number of hidden units was first conjectured in barra2017phase? based on empirical observations. The performance on the Nishimori line was then shown to be the same for \(P = P^* = 2\) hou2019minimal? as for \(P = P^* = 1\) huang2016unsupervised?, huang2017statistical?. Our PSB solution extends this result to \(\beta = \beta^*\) and any finite \(P, P^*\). On the Nishimori line, i.e. when \(\beta = \beta^*\) and \(P = P^*\), the teacher-student setting is replica symmetric, so we expect the PSB solution to be exact, which is confirmed by simulations (see Fig. 22). When \(\beta = \beta^*\) and \(P \neq P^*\), replica symmetry is not guaranteed, but the PSB solution is still in good agreement with simulations (see Fig. 23). Finally, we expect RSB corrections when \(\beta \neq \beta^*\).

Interestingly, we observe a weaker form of independence from the number of hidden units in the suboptimal partial PSB solutions, in the sense that the partial PSB solutions that we find at given \(P^*\) and \(P\) are also present at larger \(P^*\) and \(P\). These results seem related to the embedding principle stating that a neural network contains all minima of narrower neural networks with the same architecture zhang2021embedding?.

3.3.2.2 Independence of the number of hidden units: Gaussian patterns

When the student patterns \(\xi\) are real-valued variables with an i.i.d. standard Gaussian prior, Eqs. (32 ) under the PSB ansatz (Eqs. 36 ) simplify to (see 3.11) \[\begin{align} \label{eq:reduced95normal95saddle-point} m &= \frac{\hat{m}}{1 + \hat{q}} \\ q &= \frac{\hat{m}^2}{\left( 1 + \hat{q} \right)^2} + \frac{\hat{q}}{\left( 1 + \hat{q} \right)^2} \\ g &= \frac{\hat{g}}{\left( 1 + \hat{g} \right)^2} \addtocounter{equation}{1}\\ \hat{m} &= \beta^* \beta \alpha \;\mathbb{E}_{z} \left[ \tanh \left( \beta^* \beta m + \beta \sqrt{q} z \right) \right] \\ \hat{q} &= \beta^2 \alpha \;\mathbb{E}_{z} \left[ \tanh^2 \left( \beta^2 m + \beta \sqrt{q} z \right) \right] \\ \hat{g} &= \beta^2 \alpha \;\mathbb{E}_{z} \left[ \tanh^2 \left( \beta \sqrt{g} z \right) \right]. \end{align}\tag{38}\] This regime is interesting because RBMs used in practical applications usually have continous weights kiviken2012multiple?, Salakhutdinov2007restricted?, srivastava2013modelling?. Similarly as before, we compare Eqs. (38 ) against the PSB and partial PSB solutions of Eqs. (32 ) for real-valued student patterns with an i.i.d. standard Gaussian prior and teacher pattern covariance \(\mathcal{Q} = \mathbf{I}\). The resulting plots are qualitatively similar to Figs. (20) and (21), so we show them in 3.14 rather than in the main text (see Figs. 36 and 37). As before, we observe that the solutions of the Gaussian equations are independent of \(P^*\) and \(P\). One important difference between the student model with binary patterns and the student model with real-valued patterns is that they are simulated differently. In fact, binary patterns are obtained using the standard random walk Metropolis-Hastings algorithm, while real-valued patterns are learned via underdamped stochastic Langevin dynamics welling2011bayesian?, zhang2023improved? of the RBM marginal likelihood (Eq. 21 ) gradient. In a nutshell, the latter algorithm is a variant of stochastic gradient ascent where noise is added to the momentum vector to sample the RBM posterior (Eq. 22 ) rather than to find one of its modes. We use contrastive divergence, i.e. alternate Gibbs sampling of the visible and hidden units, to evaluate the marginal likelihood hinton2002training?, hinton2012practical?. At low temperature, simulations of the student model with real-valued patterns usually converge to the solution of Eqs. (38 ) (see Fig. 24). They can also sometimes stay stuck near the partial PSB solution without converging to it when the learning rate is reduced too aggressively during training. This behavior is shown in Fig. (24) by the simulation data points that are scattered around the partial PSB solution but do not agree with it within their error bars. In fact, although simulation error is useful for measuring the convergence of the Langevin learning algorithm, its sensitivity to the learning rate schedule makes it poorly representative of the equilibrium distribution of the magnetization \(m\).

3.3.2.3 A simple model of the lottery ticket hypothesis

Many works studying feedforward neural networks found a training regime where some of the hidden units learn the underlying data distribution while the others take a back seat chizat2018global?, ma2020quenching?, luo2021phase?, zhang2021embedding?. Interestingly, this behavior is similar to our PSB solution with \(P > P^*\), where only a subset of the student patterns learn the teacher patterns (see Sections 3.3.2.1 and 3.3.2.2). As proposed in zhang2021embedding?, we investigate the relationship between this kind of training regime and the lottery ticket hypothesis frankle2018lottery?.

The lottery ticket hypothesis states that a generic randomly initialized overparameterized neural network contains subnetworks that fit data with similar accuracy as the entire trained network when they are extracted from it and trained independently frankle2018lottery?, zhou2019deconstructing?, ramanujan2020hidden?, eran2020proving?, frankle2020linear?. These special subnetworks, which are commonly referred to as winning lottery tickets, are thought to have fortuitous initial conditions that facilitate training frankle2018lottery?. The PSB solution studied in Sections 3.3.2.1 and 3.3.2.2 makes it clear that any student network of size \(P > P^*\) contains subnetworks of size \(\tilde{P} \in \left\{ P^*, ..., P - 1 \right\}\) that learn the teacher patterns to the same extent, and thus fit the data at least as well. It is not obvious, however, whether any of these subnetworks can have lucky initial conditions such that they train more easily than with i.i.d. random initial conditions. As such, we apply to our teacher-student setting a variant of the magnitude pruning algorithm traditionally used to identify winning tickets frankle2018lottery?, frankle2020linear? and check if it finds a subnetwork that converges especially quickly. Consider one teacher and three students with real-valued patterns, hereby referred to as the teacher, student 0, student A and student B, respectively, where A is a control network and B is a candidate winning ticket obtained from 0. We perform the following numerical experiment:

• Initialize student 0 with \(P = 8\) i.i.d. Gaussian patterns and the teacher with \(P^* = 4\) i.i.d. Gaussian patterns. Save the initial value of the patterns of student \(0\) as \(\xi_0\).

• Train student 0 on data generated by the teacher. Compute the Euclidean norms of the learned patterns.

• Initialize student A with \(P = 4\) i.i.d. Gaussian patterns and student B with the \(P = 4\) patterns of \(\xi_0\) that evolved to have the largest Euclidean norm when trained.

• Train student A and student B on data generated by the teacher. Record their respective magnetizations \(m_A\) and \(m_B\) as a function of the training epochs in order to compare their convergence speeds.

This procedure is slightly different from the usual magnitude pruning algorithm in that it prunes the patterns with the smallest norms rather than the entries of the weight matrix with the smallest absolute values frankle2018lottery?, frankle2020linear?. As in Section 3.3.2.2, we train students 0, A and B with underdamped stochastic Langevin dynamics welling2011bayesian?, zhang2023improved? of the RBM marginal likelihood (Eq. 21 ). Our results are shown in Fig. (25). In the left panel, we verify that student B converges to the PSB solution (i.e. the solution of Eqs. 38 ) in the range of \(\alpha\) that we are studying. As in Section 3.3.2.2, the error in \(m\) is useful for measuring the convergence of the learning algorithm, but poorly representative of the equilibrium distribution. There seems to be a small systematic shift in the simulations at low \(\alpha\), which could be due to finite-size effects or order \(\frac{P}{M}\) corrections barra2017phase?. In the right panel, we plot the median of \(m_B - m_A\) over \(\alpha\) as a function of the epochs. B initially converges more quickly than A, and the lead of B eventually shrinks to zero because of the ergodicity of the training algorithm. In other words, magnitude pruning is able to identify a winning ticket that converges more quickly than networks with i.i.d. random initial conditions. The patterns of student 0 that converge to the teacher patterns the fastest have a small initial magnetization (see Fig. 25, epoch 0). Therefore, the shape of the loss function at the initial conditions and the basin of attraction in which they are located may be more important than the initial conditions themselves, as suggested by frankle2018lottery?.

3.3.3 Learning uniformly correlated patterns↩︎

In this Section, including Subsection 3.3.3.1, we introduce uniform correlations in the teacher patterns by fixing the covariance matrix \(\mathcal{Q}\) of \(\xi^*\) to \(\mathcal{Q}_{\mu \nu} = \delta_{\mu \nu} + \left( 1 - \delta_{\mu \nu} \right) c\), where \(c \in \left[ 0, 1 \right)\) controls the correlation strength. In the presence of uniform correlation, the Hamiltonian \(\mathcal{M}_*\) (Eq. 26 ) is that of the Curie-Weiss model with coupling constant \(\frac{1}{2} \left[ \beta^* \right]^2 c\). Therefore, its correlation matrix \(\mathcal{R}\) has the same form as \(\mathcal{Q}\) but with a different off-diagonal element \(d\). In 3.12, we show that the maximum eigenvalue \(\lambda^{\mathcal{S}}_{\text{max}}\) of Eq. (33 ) is \[\begin{align} \label{eq:max95eigval} \lambda^{\mathcal{S}}_{\text{max}} = \left( P^* - 1 \right)^2 c d + \left( P^* - 1 \right) \left( c + d \right) + 1. \end{align}\tag{39}\] As expected, the corresponding critical load is also the onset of nonzero Mattis magnetization in the regime of \(\beta=\beta^*\) where only the paramagnetic and ferromagnetic phases exist, which we show explicitly in Fig. (28) for binary patterns. We show in 3.14 that the same holds for real-valued student patterns with a standard Gaussian prior (see Fig. 38). Eq. (39 ) extends to arbitrary finite \(P^*\) the critical load obtained for \(P^* = 2\) in hou2019minimal?. In fact, when \(P^* = 2\), we find \(d = \tanh \left( \left[ \beta^* \right]^2 c \right)\) and Eq. (33 ) reduces to the critical load of hou2019minimal?.

Plotting Eq. (39 ) (see Fig. 27) and the corresponding \(\alpha_{\text{crit}}\) (Fig. 26) as a function of \(T\), \(c\) and \(P^*\), we see that \(\alpha_{\text{crit}}\) decreases with \(P^*\) and \(c\). For \(c \gg \left[ T^* \right]^2\), the spins of Curie-Weiss Hamiltonian \(\mathcal{M}_*\) all align, so their correlation is \(d = 1\), and we obtain \[\begin{align} \lambda^{\mathcal{S}}_{\text{max}} = \left( \left( P^* - 1 \right) c + 1 \right) P^*. \end{align}\] In this limit, \(\alpha_{\text{crit}}\) is roughly inversely proportional to both \(c\) and \(\left[ P^* \right]^2\). Conversely, for \(c \ll \left[ T^* \right]^2\), the correlation of \(\mathcal{M}_*\) is \(d = 0\), and we obtain \[\begin{align} \lambda^{\mathcal{S}}_{\text{max}} = \left( P^* - 1 \right) c + 1. \end{align}\] In this limit, \(\alpha_{\text{crit}}\) decreases less quickly by a factor of \(P^*\) than for \(c \gg \left[ T^* \right]^2\). However, relatively small correlations can still significantly decrease \(\alpha_{\text{crit}}\) when \(P^*\) is large. In other words, the critical load benefits from structured data with many correlated underlying abstract concepts even when the correlation between the different concepts is rather small. Overall, these results shed light upon the way \(\alpha_{\text{crit}}\) depends on \(P^*\), which was previously unclear given that previous work hou2019minimal? focused on \(P = P^* = 2\).

Neural networks are often regularized by feature decorrelation techniques that reduce or eliminate correlations in their inputs and hidden layers lecun1998efficient?, ioffe2015batch?, huangi2018decorrelated?, zhang2021stochastic?. Decorrelation is known to increase training speed ioffe2015batch?, huangi2018decorrelated?, zhang2021stochastic?, wiesler2011convergence?. However, our work suggests that it can also increase the critical load as a drawback, which may hinder performance when the data is noisy (i.e. \(\beta^*\) finite). galloway2019batch?, fagbohungbe2022effect? observed that batch normalization ioffe2015batch?, which was motivated by decorrelation huangi2018decorrelated?, lecun1998gradient?, makes neural networks less robust to noise. This effect could potentially be related to our findings. Arguably, the main caveat to our analysis is that we assumed that \(P^*\) was finite when deriving Eq. (39 ), so the critical load is probably different when \(P^*\) is \(\mathcal{O} \left( N \right)\) barra2017phase?.

Figs. (29) and (30) display the magnetization \(m\) and the SG overlap \(q\), respectively, found by solving Eqs. (32 ) at fixed \(\beta^*\). As usual, we find the paramagnetic (P) phase (\(m = q = 0\)) at small \(\alpha\) and low \(T\). At low \(c\) and \(P\), the spin glass (SG) phase (\(m = 0\) and \(q > 0\)) occupies the medium \(\alpha\) and low \(T\) region of the \(\alpha, T\) plane. It transitions to the ferromagnetic (F) phase (\(m \neq 0\)) as \(\alpha\), \(c\) and \(P\) increase. At high \(c\), \(P\) and \(\alpha\), the SG overlap \(q\) seems to be nonmonotonic in \(T\). As expected, the critical load of Eq. (33 ) with \(\lambda^{\mathcal{S}}_{\text{max}}\) given by Eq. (39 ) follows the onset of nonzero magnetization corresponding to the P-F phase transition, but not the SG-F phase transition. The critical load of \(\beta = \beta^*\) approximately follows the P-SG phase transition when \(c\) is small, which is consistent with Section 3.3.2.1 and previous works alemanno2023hopfield?, theriault2024dense?. As in Section 3.3.2.1, decreasing the inference temperature \(T\) too much can make it harder for the student to learn the teacher patterns (see Fig. 29, top right panel).

3.3.3.1 Permutation symmetry breaking transitions

As in Section 3.3.2, the critical load marking the onset of learning (see Eq. 33 ) again does not depend on the number of hidden units \(P\) of the student RBM. Despite this, a single wide RBM does not learn teacher patterns as would multiple separate RBMs with one hidden unit each. Moreover, distinct phases emerge based on the level of correlation \(c\) and data load \(\alpha\), each characterized by a different learning strategy.

RBMs with \(P=P^*=2\) hidden units learning correlated patterns were previously found to have three distinct ferromagnetic phases hou2019minimal?:

• in the Spontaneous Symmetry Breaking (SSB) phase, both student patterns converge to a single configuration that has the same overlap with the two teacher patterns. In other words, the student learns only the features that the teacher patterns have in common;

• in the student Permutation Symmetry Breaking (\(\text{PSB}_{\text{s}}\)) phase, the student patterns converge to two distinct configurations that have the same overlap with the teacher patterns. In other words, the student is able to learn distinct features of the teacher patterns;

• finally, in the teacher Permutation Symmetry Breaking (\(\text{PSB}_{\text{t}}\)) phase, the student can tell the two teacher patterns apart and learns them one-to-one as in the PSB phase of Section (3.3.2.1).

These three phases can be identified by comparing the values of the order parameters \(m^{\mu\nu}\) and \(q^{\mu\nu}\) on and off the diagonal. As shown in the left panel of Fig. (31), bifurcations occur as \(\alpha\) increases, marking second-order phase transitions at the critical loads \(\alpha^{\text{SSB}}_{\text{crit}} = \alpha_{\text{crit}}\), \(\alpha^{\text{PSB}_{\text{s}}}_{\text{crit}} > \alpha^{\text{SSB}}_{\text{crit}}\) and \(\alpha^{\text{PSB}_{\text{t}}}_{\text{crit}} > \alpha^{\text{PSB}_{\text{s}}}_{\text{crit}}\).

In this context, the student does not know the level of correlation \(c\) of the teacher patterns. As such, it uses a uniform prior \(\mathop{\mathrm{P}}\left( \xi \right) \neq \mathop{\mathrm{P}}\left( \xi^* \right)\), and we are always outside the Nishimori regime. Despite this, the RS critical thresholds and magnetizations are in good agreement with Monte Carlo simulations (see Fig. 31, left panel).

As shown in Figs. (31) and (32), the same three phases also appear at larger \(P\) and \(P^*\). \(\alpha^{\text{SSB}}_{\text{crit}}\) once again coincides with the onset of the ferromagnetic phase, which occurs at the same load \(\alpha_{\text{crit}}\) regardless of the number of hidden units \(P\) of the student. On the other hand, \(\alpha^{\text{PSB}_{\text{s}}}_{\text{crit}}\) is smaller for \(P = 3\) than for \(P = 2\) (see Fig. 31). We think that the term \(A_{\mu \mu} \left( q \right) = \sqrt{2 q^{\mu \mu} - \sum_{\eta = 1}^P q^{\mu \eta}}\) in \(\mathcal{L}_{\lambda_1 \lambda_2}\) (see Eqs. 30 and 31 ) penalizes values of \(q\) with a large off-diagonal sum, i.e. \(\sum_{\eta \neq \mu}^P q^{\mu \eta}\), compared to the diagonal \(q^{\mu \mu}\). As the number of terms in the sum grows with \(P\), the individual off-diagonal coefficients \(q^{\mu\nu}\) may be encouraged to become smaller, pushing permutation symmetry breaking of the student patterns (\(\text{PSB}_{\text{s}}\)) to occur at a lower critical threshold \(\alpha^{\text{PSB}_{\text{s}}}_{\text{crit}}\).

In hou2019minimal?, the critical load \(\alpha^{\text{PSB}_{\text{t}}}_{\text{crit}}\) increases with the correlation strength \(c\). As such, increasing \(c\) can trigger a phase transition from \(\text{PSB}_{\text{t}}\) to \(\text{PSB}_{\text{s}}\) and thus decrease the magnetization. In other words, large correlations in the teacher patterns can undermine the student’s ability to learn them accurately. Based on these findings, Hou et al. hou2019minimal? formulated the hypothesis that the best \(c\) for learning a relatively low load of data is nonzero but still relatively small. In the case of \(P, P^* \geq 3\) (see Fig. 32), we also find that \(\alpha^{\text{PSB}_{\text{t}}}_{\text{crit}}\) grows with \(c\). However, unlike \(P = P^* = 2\), Monte Carlo simulations of \(P, P^* \geq 3\) do not completely agree with the RS approximation at large \(c\). In particular, increasing \(c\) does not seem to decrease the magnetization significantly in the simulations, even when it does so according to the RS approximation (see Fig. 32, black crosses).

Repeating the lottery ticket experiment of Section 3.3.2.3 for \(c = 0.05\), we find that the lead of the winning ticket B over its randomly-initialized counterpart A has a smaller maximum (see Fig. 33). In other words, small correlations in the data appear to make the student less sensitive to initial conditions, which also reduces the benefits of magnitude pruning (see Section 3.3.2.3). As in Sections 3.3.2.2 and 3.3.2.3, the error in \(m\) (see Fig. 33, left panel) is useful for measuring the convergence of the learning algorithm, but poorly representative of the equilibrium distribution.

Outside the teacher-student setting, RBMs typically undergo a sequence of second-order phase transitions at the beginning of training decelle2017spectral?, decelle2018thermodynamics?, bachtis2025cascade?, bereux2025fast?. After the first transition, they learn a rank-one matrix, which is reminiscent of our SSB phase. The second transition then breaks the rank-one symmetry like our \(\text{PSB}_{\text{s}}\) transition. The subsequent phase transitions are harder to interpret from the point of view of permutation symmetry breaking, and it is not clear whether a transition analogous to \(\text{PSB}_{\text{t}}\) can occur without an explicit teacher. \(\text{PSB}_{\text{s}}\) phase transitions were also observed in Gaussian mixture models rose1990statistical?, kloppenburg1997deterministic?, akaho2000nonmonotonic?, Decelle2021restricted? and modern Hopfield networks boukacem2024waddington?, which suggests that they could be common in machine learning models.

3.3.4 Random correlations↩︎

In this section, we take the covariance matrix \(\mathcal{Q}\) to be random. By definition, \(\mathcal{Q}\) must be positive semi-definite. Moreover, when the teacher weights are binary, \(\mathcal{Q}\) must have ones on the diagonal. We sample \(\mathcal{Q}\) from the projected Wishart distribution \(\mathcal{W} \left( c, D \right)\) defined in 3.5.2 because it satisfies these two requirements. By the law of large numbers, \(\mathcal{Q} \sim \mathcal{W} \left( c, D \right)\) approaches \(\mathcal{C} = \delta_{\mu \nu} + \left( 1 - \delta_{\mu \nu} \right) c\) in probability as \(D \rightarrow \infty\) (see 3.5.2). Therefore, \(\lambda^{\mathcal{S}}_{\text{max}}\) and \(\alpha_{\text{crit}}\) are the same as in Section 3.3.3 when \(D\) is large. At finite \(D = P\), the dependence of \(\alpha_{\text{crit}}\) on \(T\) is still qualitatively similar to that of uniform correlations. However, the arithmetic mean \(\overline{\alpha_{\text{crit}}}\) and the harmonic mean \(\left[ \overline{1 / \lambda^{\mathcal{S}}_{\text{max}}} \right]^{-1}\) over many independent samples \(\mathcal{Q} \sim \mathcal{W} \left( c, P = D \right)\) are very different from the \(\alpha_{\text{crit}}\) and \(\lambda^{\mathcal{S}}_{\text{max}}\) of uniform correlations as a function of \(c\) and \(P\) (see Figs. 26 and 27). For instance, \(\overline{\alpha_{\text{crit}}}\) decreases much more slowly with \(c\) than the critical load of uniform correlations. At small \(c\) and high \(T\), \(\overline{\alpha_{\text{crit}}}\) tends to be smaller than the critical load of uniform correlations. Conversely, at large \(c\) and low \(T\), the critical load of uniform correlations is usually smaller. At any given \(c\), the entries of a random correlation matrix can sometimes be larger than \(c\), which can make learning possible even when \(T\) is too high for correlations of size \(c\) or smaller to be picked up on by the student. However, this advantage disappears at larger \(c\) and lower \(T\) where correlations of size \(c\) and smaller are no longer muddled in noise. In summary, the behavior of the critical load as a function of \(P\) and \(c\) is in general very different for random correlations and for uniform correlations.

3.4 Conclusion↩︎

In this paper, we theoretically study the learning performance of Restricted Boltzmann Machines (RBM) ackley1985learning?, smolensky1986information?, freund1991unsupervised?, hinton2002training? with a finite number of hidden units in the teacher-student setting huang2016unsupervised?, huang2017statistical?, barra2017phase?, barra2018phase?, huang2018role?, hou2019minimal?, decelle2021inverse? using the replica method under the Replica-Symmetric (RS) approximation charbonneau2022replica?. Given \(M\) data points and \(N\) visible units, we compute the critical data load \(\alpha_{\text{crit}} = \frac{M}{N}\) above which learning becomes possible. Our findings validate the conjecture that the student’s performance is independent of the number of hidden units when the patterns of the teacher are uncorrelated barra2017phase?, generalizing the results of hou2019minimal? to any number of hidden units \(P\) and teacher patterns \(P^*\) much smaller than \(M\) and \(N\). In particular, we confirm that an RBM with \(P\) uncorrelated hidden units factorizes into \(P\) RBMs with one hidden unit each hou2019minimal?. Additionally, we show that the teacher-student setting without correlations has both a Permutation Symmetry Breaking (PSB) solution in which the student learns the teacher patterns one-to-one and a metastable partial PSB solutions in which multiple student patterns can converge to the same teacher pattern. We argue that the teacher-student setting complies with the lottery ticket hypothesis frankle2018lottery?, zhou2019deconstructing?, ramanujan2020hidden?, eran2020proving?, frankle2020linear?, and demonstrate that the student can be trained efficiently using a variant of the magnitude pruning algorithm frankle2018lottery?, frankle2020linear?. Given teacher patterns with uniform correlations \(c \in [0, 1)\), we find a closed-form expression for \(\alpha_{\text{crit}}\), which generalizes the one found in hou2019minimal? for \(P = P^* = 2\), and show that it decreases with both \(c\) and \(P^*\). Still in the case of uniform correlations, we find that the teacher-student setting undergoes a sequence of Spontaneous Symmetry Breaking (SSB) and permutation symmetry breaking transitions that generalizes the one found in hou2019minimal? to \(P^*, P \geq 3\). Both with and without correlations, we find that decreasing the inference temperature \(T\) too much can prevent the student from learning the teacher’s solution even when the dataset is relatively large. Throughout the paper, we compare key results against Monte Carlo simulations and observe that the RS ansatz is often a solid approximation even in the region outside the Nishimori line nishimori1980exact?, iba1999nishimori?, nishimori2001statistical?, contucci2009spin? where the teacher and the student have the same temperature, but a different number of hidden units and a different prior on their patterns.

One could study the learning dynamics of the magnitude pruning experiment described in Section 3.3.2.3 or characterize the distribution of winning ticket initial conditions to gain additional insight into the lottery ticket hypothesis. More generally, it would be interesting to extend our study to the case where \(P\) and \(P^*\) are of order \(N\). In fact, RBMs trained on real data typically have a large number of hidden units hinton2002training?, leroux2011learning?, tubiana2017emergence?, kiviken2012multiple?, fischer2014training?. One could also introduce correlations between the columns \(\xi^*_i = \left\{ \xi^{* \mu}_i \right\}_{\mu = 1}^P\) of the teacher’s weight matrix \(\xi^* = \left\{ \xi^{* \mu}_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\) in addition to the correlations between patterns \(\xi^{* \mu} = \left\{ \xi^{* \mu}_i \right\}_{i = 1}^N\) that we study in this work. Such a modification could be achieved by using a weight matrix with a low-rank SVD structure like in decelle2017spectral?, decelle2018thermodynamics? and would allow the teacher to generate data with even more structure than in this paper. One could use this framework to investigate whether the distinction between \(\text{PSB}_{\text{s}}\) and \(\text{PSB}_{\text{t}}\) persists in real data and has a measurable effect on the critical slowing down patterns observed in RBM training decelle2017spectral?, decelle2018thermodynamics?, bachtis2025cascade?, bereux2025fast?. Another promising research avenue would be to study other generative models with a finite number of hidden units in the teacher-student setting, such as modern (a.k.a. dense) Hopfield networks hopfield1982neural?, theriault2024dense?, chen1986high?, psaltis1986nonlinear?, krotov2016dense?. In sum, the teacher-student setting still has a lot of untapped potential for studying generative models with many hidden units.

Data availability↩︎

The code and hyperparameter values of the training algorithms used to make the figures are available at the following public repository theriault2025modellingsoftware?.

Acknowledgements↩︎

I am grateful to Francesco Tosello, who made a significant contribution to this work. He wrote a first draft of the calculations presented in the Appendices, which Daniele Tantari and I then generalized and expanded upon to write this Chapter.

4 Saddle hierarchy in dense associative memory↩︎

4.1 Introduction↩︎

Studying the stationary points of machine learning algorithms is crucial to understand how they work. For example, auffinger2013random?, choromanska2015loss? demonstrated that local minima in the loss landscape of large artificial Neural Networks (NNs) are relatively close to the global minimum, explaining why they generalize well in practice. Moreover, razvan2014saddle?, dauphin2014identifying? showed that saddle points are much more numerous than local minima in large NNs. These breakthroughs were made by establishing deep connections between the loss landscape of NNs and the energy landscape of disordered systems studied in statistical mechanics. Beyond the broad insights provided by these studies, and despite the progress made by FUKUMIZU2000local?, fukuzumi2019semi-flat?, simsek2021geometry?, zhang2021embedding?, the classification of stationary points in machine learning algorithms remains an open problem.

Energy-based models have been playing a central role in studies of NNs and their theoretical properties. The Hopfield network, one of the most historically important energy-based models, was originally introduced as a paradigmatic model of biological associative memory hopfield1982neural?. Generalized Hopfield networks chen1986high?, psaltis1986nonlinear? were then developed to improve upon the limited storage capacity of the original hopfield1982neural?, amit1985storing?. The scale of the improvement was determined rigorously in following studies baldi1987number?, gardner1987multiconnected?, abbott1987storage?. A few years ago, these generalized networks, commonly referred to as Dense Associative Memory (DAM) or modern Hopfield networks, were made into trainable machine learning models capable of accurate pattern classification by Krotov and Hopfield (K & H) krotov2016dense?. In a nutshell, K & H’s DAM learns prototypes of patterns in a trainable weight matrix. Each prototype casts a vote for a class, and the patterns awaiting classification are assigned based on the votes of the prototypes that most closely resemble them. The resulting classification scheme is considerably more adversarially robust krotov2018dense?, theriault2024dense? and interpretable krotov2016dense? than that of feedforward NNs with ReLU activation functions. Since their debut as trainable machine learning architectures, deep connections have been made between modern Hopfield networks and other well-known and established machine learning paradigms, such as attention ramsauer2020hopfield? and generative diffusion hoover2023memory?, ambrogioni2024search?. In particular, modern Hopfield networks were used to implement the attention mechanism of transformers ramsauer2020hopfield?, which has been attracting a lot of interest in fundamental and applied research. See krotov2023new?, krotov2025modern? for a review of recent advances on the topic. Recently, it was observed that the trainable weights of K & H’s DAM are channeled toward minima by a low-dimensional network of valleys in the loss landscape boukacem2024waddington?. Moreover, the points where valleys branch out from one another were identified as saddles in the simple case where the DAM has two patterns to learn. In general, it is not straightforward to classify the stationary points of machine learning algorithms zhang2021embedding?. However, the results of boukacem2024waddington? and the interpretability of DAMs suggest that their stationary points are both fundamental to their learning dynamics and easier to characterize than that of generic NNs. With this goal in mind, we revisit dense associative memory for pattern classification krotov2016dense? using the framework of Boltzmann Machines (BMs) ackley1985learning?, smolensky1986information?, freund1991unsupervised?, hinton2002training? and statistical mechanics.

The body of this paper is structured as follows. In Section 4.2, we introduce the DAM model that we study and the analytical tools that we use throughout our work. In particular, in Section 4.2.1, we derive the DAM in question from a BM template, and, in Section 4.2.2, we explain the setting that we use to analyze DAM stationary points.

Section 4.3 presents our theoretical contributions. In Section 4.3.1, we formulate saddle-point equations that characterize DAM stationary points, and, in Section 4.3.2, we leverage these equations to establish a saddle-point hierarchy principle stating that weights learned by DAMs of a given width are embedded as saddle points in wider DAMs.

Section 4.4 explains how the theoretical insights of Section 4.3 can be used to enhance training and explores how DAMs represent data once trained. In Section 4.4.1, we introduce a regularization method that facilitates supervised learning. Next, in Section 4.4.2, we show that our DAM, although designed for supervised tasks, can discover interpretable solutions in both supervised and unsupervised classification settings. Finally, in Section 4.4.3, we link our findings to the learning dynamics shaped by valleys and saddles studied in boukacem2024waddington?, and we implement a network-growing algorithm wu2019splitting? that exploits the saddle-point hierarchy to significantly reduce DAM training costs.

The code and hyperparameter configurations used in our experiments are available in the following public repository theriault2025saddlesoftware?.

4.2 Model↩︎

A Boltzmann Machine is a canonical graphical model of correlations in discrete data ackley1985learning?. It is customary to partition BMs into a visible layer \(\mathbf{v} = \left\{ v_i \right\}_{i = 1}^N \in \mathbb{R}^N\) and a hidden layer \(\mathbf{h} = \left\{ h_\mu \right\}_{\mu = 1}^P \in \mathbb{R}^P\) such that connections between the two layers are allowed, but connections within them are prohibited smolensky1986information?. In this case, the visible layer represents concrete features of the data, whose mutual correlations are encoded in connections with the hidden layer. The Restricted Boltzmann Machine (RBM) obtained using this partition is much easier to train than a generic BM freund1991unsupervised?, hinton2002training? and still has considerable generating power freund1991unsupervised?, leroux2008representational?, making it more practical in machine learning applications Salakhutdinov2007restricted?, kiviken2012multiple?, srivastava2013modelling?. The visible and hidden units of an RBM follow the Gibbs distribution \[\begin{align} \mathop{\mathrm{P}}_\beta \left( \mathbf{v}, \mathbf{h} \big| \mathbf{J} \right) &= Z_\beta \left( \mathbf{J} \right)^{-1} \mathop{\mathrm{P}}_0 \left( \mathbf{v} \right) \mathop{\mathrm{P}}_0 \left( \mathbf{h} \right) \exp \left( -\beta H \left[ \mathbf{v}, \mathbf{h} ; \mathbf{J} \right] \right), \end{align}\] where \(\beta \geq 0\) is known as the inverse temperature, \(\mathop{\mathrm{P}}_0 \left( \mathbf{v} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} \right)\) are priors on \(\mathbf{v}\) and \(\mathbf{h}\), \(H \left[ \mathbf{v}, \mathbf{h} ; \mathbf{J} \right] = -\sum_{i = 1}^N \sum_{\mu = 1}^P J^{\mu}_i v_i h_\mu\) is called the energy function or Hamiltonian, \(\mathbf{J} = \left\{ J^\mu_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\) are trainable weights, and \(Z_\beta \left( \mathbf{J} \right)\) is a normalization constant called the partition function. The inverse temperature \(\beta\) represents the absolute strength of the RBM connections, or equivalently controls the amount of noise \(T = 1/\beta\) in the RBM. In this regard, \(\mathop{\mathrm{P}}_0 \left( \mathbf{v} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} \right)\), which restrict the form of the Gibbs distribution to help the RBM represent the data, are the marginal laws of \(\mathbf{v}\) and \(\mathbf{h}\) when there are no connections, i.e. \(\beta = 0\). Their contribution to the Gibbs distribution can be tuned with \(\beta\), which can therefore be interpreted as a regularization parameter.

4.2.1 A Dense Associative Memory (DAM) model↩︎

As mentioned in the Introduction, we will now derive a DAM model for classification from a BM template. We will explain why our model is a DAM at the end of this Section, once we have clearly defined it. We make three basic assumptions on the distribution of data to be classified:

1. the data is scale invariant, i.e. for all positive scalars \(c\), data points \(\mathbf{x}\) and \(c \mathbf{x}\) are equivalent;

2. the data can be partitioned in disjoint clusters;

3. the clusters can be grouped into mutually exclusive classes.

In order to exploit these three assumptions, we study a BM partitioned into three layers with different roles: the data layer \(\mathbf{x}\), the hidden layer \(\mathbf{h}\) and the class layer \(\mathbf{q}\), which represent data, cluster membership and class membership, respectively. The corresponding energy is \[\begin{align} \label{eq:RBM95energy} - H \left[ \mathbf{x}, \mathbf{q}, \mathbf{h} ; \mathbf{J} \right] =\sum_{i = 1}^N \sum_{\mu = 1}^P w^\mu_i x_i h_\mu + \sum_{y = 1}^C \sum_{\mu = 1}^P u^\mu_y q_y h_\mu + \sum_{\mu = 1}^P h_\mu b^\mu, \addtocounter{equation}{1}\end{align}\tag{47}\] where \(\mathbf{J} = \left\{ \mathbf{w}, \mathbf{u}, \mathbf{b} \right\}\) is the set of the trainable weights \(\mathbf{w} = \left\{ w^\mu_i \right\}_{1 \leq i \leq N}^{1 \leq \mu \leq P}\), \(\mathbf{u} = \left\{ u^\mu_i \right\}_{1 \leq i \leq C}^{1 \leq \mu \leq P}\) and \(\mathbf{b} = \left\{ b^\mu \right\}_{\mu = 1}^P\). There are no direct interactions between the visible layer and the class layer. In other words, conditional on the cluster layer, the visible layer and the class layer are independent. Therefore, this BM is a Deep Boltzmann machine (DBM) with 3 layers salakhutdinov2009deep?, which can also be thought of as an RBM whose visible layer \(\mathbf{v}\) is further divided into \(\mathbf{x}\) and \(\mathbf{q}\).

Since the data is scale invariant (Assumption [assum:invariance]), the scale of individual data points contains no information about the classification task, so we normalize them by their (Euclidean) norm in the data layer. In other terms, we take the data units \(x_i\) to be continuous variables with unit norm \(\sqrt{\sum_{i = 1}^N \left( x_i \right)^2} = 1\). We assume no further knowledge about \(\mathbf{x}\), so we take the prior \(\mathop{\mathrm{P}}_0 \left( \mathbf{x} \right)\) to be the uniform distribution on the \(N-1\) dimensional unit hypersphere \(S^{N-1}\), i.e. is the set of all \(\mathbf{x}\) with \(\sqrt{\sum_{i = 1}^N \left( x_i \right)^2} = 1\). Data normalization is a very common practice in machine learning. For example, normalization by the Euclidean norm has been popular in text document clustering even since its introduction in the 1980s salton1983introduction?. Various types of normalization also occur in the brain and retina Carandini2012normalization?.

Since the hidden layer and the class layer aim to represent disjoint clusters and classes, respectively (Assumptions [assum:clustering] and [assum:labeling]), we take their respective units to be mutually exclusive binary variables, i.e. \(\mathbf{h}\in \{0,1\}^P\) with \(\sum_{\mu=1}^P h_\mu\in\{0,1\}\) and \(\mathbf{q}\in \{0,1\}^C\) with \(\sum_{y=1}^C q_y\in\{0,1\}\). At any given time, at most one unit per layer can take the value \(1\), representing the fact that the clusters and classes are disjoint. In other words, we take each of the two layers to be the vector representation of a single categorical (or Potts potts1952some?, wu1982potts?) variable with \(P+1\) and \(C+1\) categories, respectively. As such, \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{q} \right)\) simplify to probability mass functions \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} = \mathbf{e}_{\gamma} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{q} = \mathbf{e}_{y} \right)\), where we introduce \(\mathbf{e}_\gamma = \left\{ \delta_{\gamma \mu} \right\}_{\mu = 1}^P\) for \(\gamma \in \left\{ 0,\ldots, P \right\}\) and define \(\mathbf{e}_y\in\{0,1\}^C\) analogously for \(y\in\{0,\ldots , C\}\). In particular \(\mathbf{e}_0=\mathbf{0}\) represents a state outside the \(P\) clusters or the \(C\) classes.

These priors on the hidden layer and the class layer can also be obtained by introducing fixed inhibitory connections within the hidden layer and the class layer, respectively kappen1993using?, KAPPEN1995deterministic?. Since at most one hidden unit \(h_\mu\) can be activated at once, the hidden layer is a very sparse representation of the visible layer. In machine learning, sparsity can improve interpretability mozer1988skeletonization?, generalization, computational efficiency hoefler2021sparsity?, and adversarial robustness guo2018sparse?. The sparsity of the brain suggests that it is also beneficial for biological neural networks friston2008hierarchical?.

Given these priors \(\mathop{\mathrm{P}}_0 \left( \mathbf{x} \right)\), \(\mathop{\mathrm{P}}_0 \left( \mathbf{h} \right)\) and \(\mathop{\mathrm{P}}_0 \left( \mathbf{q} \right)\), we derive the marginal distribution of the visible layer \(\left( \mathbf{x}, \mathbf{q} \right)\) (see Appendix 4.6 for details). We start our derivation by showing that the conditional distribution of the data layer given the hidden layer has the form \[\begin{align} \label{eq:vmf95distribution} \mathop{\mathrm{P}}_{\beta} \left( \mathbf{x} | \mu, \mathbf{J} \right) &:= \mathop{\mathrm{P}}_{\beta} \left( \mathbf{x} | \mathbf{h} = \mathbf{e}_\mu, \mathbf{J} \right) \addtocounter{equation}{1}\\ &\propto \exp \left( \beta \sum_{i = 1}^N w^\mu_i x_i \right) \quad \forall \;\mathbf{x}\in S^{N-1} \;\text{and} \;\mu \in \left\{ 1, ..., P \right\}. \end{align}\tag{48}\] In other words, the probability density \(\mathop{\mathrm{P}}_{\beta} \left( \mathbf{x} | \mu, \mathbf{J} \right)\) corresponding to each cluster \(\mu > 0\) is a von Mises-Fisher (vMF) distribution centered on the direction \(\mathbf{w}^\mu = \left\{ w^\mu_i \right\}_{i = 1}^N\) (see Appendix 4.7). In order to interpret the \(\mathbf{w}^\mu\) as centroids for their respective clusters, we assume that they belong to \(S^{N - 1}\) like the data layer. Under this assumption, we find the normalization constant of \(\mathop{\mathrm{P}}_{\beta} \left( \mathbf{x} | \mu, \mathbf{J} \right)\) to be \(\Omega_N \left( \beta \right) = \frac{\left( 2\pi \right)^{N/2} I_{N/2 - 1} \left( \beta \right)}{\beta^{N/2 - 1}}\), where \(I_n \left( x \right)\) is the modified Bessel function of the first kind of order \(n\) (see Appendix 4.7).

After slightly more work, we find the marginal distribution of the visible layer to be \[\begin{align} \label{eq:direct95distribution} \mathop{\mathrm{P}}_\beta \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right) &:= \mathop{\mathrm{P}}_\beta \left( \mathbf{x}, \mathbf{q}=\mathbf{e}_y \big| \mathbf{J} \right) \\ &= \sum_{\mu = 1}^P p^\mu_y \frac{\exp \left( \beta \sum_{i = 1}^N w^\mu_i x_i \right)}{ \Omega_N \left( \beta \right)} + p^0_y \frac{1}{\Omega_N \left( 0 \right)} \quad \forall \;\mathbf{x}\in S^{N-1} \;\text{and} \;y \in \left\{ 0, ..., C \right\},\addtocounter{equation}{1}\end{align}\tag{49}\] where \(\mathbf{p} = \left\{ p^\gamma_y \right\}_{0 \leq y \leq C}^{0 \leq \gamma \leq P} = \left\{ \mathop{\mathrm{P}}_\beta\left(\mathbf{q}=\mathbf{e}_y,\mathbf{h}=\mathbf{e}_\gamma|\mathbf{J}\right) \right\}_{0 \leq y \leq C}^{0 \leq \gamma \leq P}\). \(\Omega_N \left( 0 \right) = \frac{2 \pi^{N/2}}{\Gamma \left( N/2 \right)}\) is the surface area of \(S^{N-1}\), where \(\Gamma \left( x \right)\) is the Gamma function, so \(\frac{1}{\Omega_N \left( 0 \right)}\) is the uniform distribution on \(S^{N - 1}\). The uniform distribution term of Eq. (49 ) encourages the model to ignore very noisy data during training, which may or may not be desirable depending on the application.

The detailed derivation of Eq. (49 ) (in Appendix 4.6) is inspired by the derivation of Gaussian mixtures from RBMs presented in Decelle2021restricted? based on the framework of kappen1993using?, KAPPEN1995deterministic?. In our case, the marginal distribution \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | \mathbf{w}, \mathbf{p} \right) = \sum_{y = 0}^C \mathop{\mathrm{P}}_\beta \left( \mathbf{x}, y | \mathbf{w}, \mathbf{p} \right)\) and conditional distributions \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | y ; \mathbf{w}, \mathbf{p} \right) = \frac{\mathop{\mathrm{P}}_\beta \left( \mathbf{x}, y | \mathbf{w}, \mathbf{p} \right)}{\sum_{\gamma = 0}^P \mathop{\mathrm{P}}_\beta \left( y, \gamma | \mathbf{w}, \mathbf{p} \right)}\) are convex mixtures between the uniform distribution on \(S^{N - 1}\) and the vMF distributions \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | \mu, \mathbf{J} \right)\) (see Eq. 48 ). Probabilistic models similar to \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | \mathbf{w}, \mathbf{p} \right)\) are notably used in text document clustering banerjee2005clustering?. Mixture distributions that have class-dependent weights like \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | y ; \mathbf{w}, \mathbf{p} \right)\) are also used in Gaussian mixture discriminant analysis hastie1996discriminant?.

The class weights \(\mathbf{p}\) depend on the trainable parameters \(\mathbf{u}\) and \(\mathbf{b}\) of Eq. (47 ) (see Appendix 4.6). Without loss of generality, we choose to directly study (and train) \(\mathbf{p}\) instead of \(\mathbf{u}\) and \(\mathbf{b}\). Recall that \(p^\gamma_y\) is a probability distribution, and in particular \(p^\gamma_y \geq 0\) (see Eq. 49 ). We constrain the marginal \(\sum_{y = 0}^C p^\gamma_y\), i.e. the fraction of data in each cluster, to a fixed distribution \(p_{\mathbf{h}} \left( \gamma \right)\) and the marginal \(\sum_{\gamma = 0}^P p^\gamma_y\), i.e. the proportion of data in each class \(y\), to another fixed distribution \(p_{\mathbf{q}} \left( y \right)\). In sum, we end up with the constraints \(p^\gamma_y \geq 0\), \(\sum_{\gamma = 0}^P p^\gamma_y = p_{\mathbf{q}} \left( y \right)\) and \(\sum_{y = 0}^C p^\gamma_y = p_{\mathbf{h}} \left( \gamma \right)\). Since each cluster belongs to a single class (Assumption [assum:labeling]), we expect the \(p^\gamma_y\) of a trained model to be close to \(\sum_{y^\prime = 0}^C p^\gamma_{y^\prime}\) for a given \(y\) and close to \(0\) otherwise.

Given a dataset of \(P^*\) patterns \(\left\{ \mathbf{x}^{* \mu} \right\}_{\mu = 1}^{P^*}\) with soft labels \(q^{* \mu}_y\) Szegedy2016rethinking?, we can train the weights \(\mathbf{w}\) and \(\mathbf{p}\) by minimizing the negative log-likelihood loss \[\label{eq:loss} L \left( \mathbf{w}, \mathbf{p} \right) = -\frac{1}{P^*} \sum_{\mu = 1}^{P^*} \sum_{y = 0}^C q^{* \mu}_y \log \mathop{\mathrm{P}}_\beta \left( \mathbf{x}^{* \mu}, y | \mathbf{w}, \mathbf{p} \right),\tag{50}\] which is a form of maximum likelihood estimation. We do so using constrained Stochastic Gradient Descent with momentum (simply called SGD in this paper). We explain how the constraints on \(\mathbf{w}\) and \(\mathbf{p}\) are enforced in Appendix 4.8, and we briefly discuss the initial conditions and the learning rate in Appendix 4.9.

The weights \(\mathbf{w}^\mu\) learned by SGD of the loss (Eq. 50 ) lecun1998gradient? are interpretable prototypes, or memories, of the data (see Fig. 40 in the case of the MNIST dataset of handwritten digits lecun1998gradient?), which is consistent with their role as cluster centroids. Once the model is trained, we can use its conditionals \(\mathop{\mathrm{P}}_\beta \left( y | \mathbf{x}; \mathbf{w}, \mathbf{p} \right)\) and \(\mathop{\mathrm{P}}_\beta \left( \mathbf{x} | y; \mathbf{w}, \mathbf{p} \right)\) to efficiently reconstruct \(y\) from \(\mathbf{x}\) or \(\mathbf{x}\) from \(y\), respectively. In particular, we can classify unseen patterns using the Bayes classification rule \(y = \mathop{\mathrm{argmax}}_{y^\prime} \left\{ \log \mathop{\mathrm{P}}_\beta \left( y^\prime | \mathbf{x} ; \mathbf{w}, \mathbf{p} \right) \right\}\) and reconstruct patterns of a given class by finding the local minima of the effective energy \(-\log \mathop{\mathrm{P}}_\beta \left( \mathbf{x} | y ; \mathbf{w}, \mathbf{p} \right)\) as a function of \(\mathbf{x} \in S^{N - 1}\).

High-dimensional probabilistic models that store or learn prototypes, such as our model, can typically reconstruct a limited number of patterns with nontrivial accuracy. In other words, they have limited storage capacity. For instance, Hopfield’s model of associative memory, the Hopfield network hopfield1982neural?, has a capacity of up to \(\mathcal{O} \left( N \right)\) patterns hopfield1982neural?, amit1985storing?. Dense Associative Memory (DAM) models, which are inspired by the Hopfield network, are a class of models with an asymptotically much higher capacity gardner1987multiconnected?, krotov2016dense?. The effective energy \(-\log \mathop{\mathrm{P}}_\beta \left( \mathbf{x} | y ; \mathbf{w}, \mathbf{p} \right)\) that we can minimize to reconstruct patterns with our model (see Eq. 49 ) is very similar to that of the DAM with exponential capacity introduced in ramsauer2020hopfield?, which uses energy minimization to implement the attention mechanism of transformers. In fact, according to lucibello2024exponential?, our model also belongs to the class of DAMs with exponential capacity.

4.2.2 Teacher-student setting↩︎

Among the extensive research on the properties of artificial Neural Networks (NNs) from the perspective of statistical mechanics (see charbonneau2023spin? for a review), there have been many studies of simple RBMs trained by maximum likelihood estimation (a line of work pioneered in decelle2017spectral?, decelle2018thermodynamics?, Decelle2020gaussian?) or by averaging samples from the posterior distribution of the weights \(\mathbf{J}\) given some observed data \(\mathcal{D}\). The latter approach has notably been used to characterize the fundamental limits of learning for many types of RBMs huang2016unsupervised?, huang2017statistical?, barra2017phase?, with different priors barra2018phase?, huang2018role?, manzan2025effect? and architectures hou2019minimal?, theriault2025modeling?, in the teacher-student setting gardner1989unfinished?, decelle2021inverse?, alemanno2023hopfield? where the data used to train \(\mathbf{J}\) is sampled from another RBM with planted weights \(\mathbf{J}^*\). This teacher-student setting can also be used to study our DAM. In this scenario, a teacher DAM with weights \(\mathbf{w}^*\) and \(\mathbf{p}^*\) generates a large amount \(M = \alpha N\) of noisy data \(\mathcal{D} = \left\{ \mathbf{x}^c, y^c \right\}_{c = 1}^M\) and feeds them to a student DAM, which then trains its weights \(\mathbf{w}\) and \(\mathbf{p}\) by averaging samples from the posterior distribution \[\begin{align} \label{eq:posterior} \mathop{\mathrm{P}}_\beta \left( \mathbf{w}, \mathbf{p} | \mathcal{D} \right) &= \mathcal{Z}_\beta \left( \mathcal{D} \right)^{-1} \mathop{\mathrm{P}}_0 \left( \mathbf{w}, \mathbf{p} \right) \prod_{c = 1}^M \mathop{\mathrm{P}}_\beta \left( \mathbf{x}^c, y^c | \mathbf{w}, \mathbf{p} \right), \addtocounter{equation}{1}\end{align}\tag{51}\] where \(\mathcal{Z}_\beta \left(\mathcal{D} \right) = \mathbb{E}_{\mathbf{w}, \mathbf{p}} \left[ \prod_{c = 1}^M \mathop{\mathrm{P}}_\beta \left( \mathbf{x}^c, y^c | \mathbf{w}, \mathbf{p} \right) \right]\) is the posterior partition function and \(\mathop{\mathrm{P}}_0 \left( \mathbf{w}, \mathbf{p} \right)\) is the prior on \(\mathbf{w}\) and \(\mathbf{p}\), which for simplicity we choose as uniform over the sets in which \(\mathbf{w}\) and \(\mathbf{p}\) are constrained.

We use a statistical mechanics approach to derive a single set of saddle-point equations that simultaneously characterize the weights that are stationary points of maximum likelihood estimation (Eq. 50 ) for generic data and the typical weight configurations obtained by averaging samples from Eq. (51 ) in the teacher-student setting. We assume that the student does not know the number of hidden units \(P^*\) and the inverse temperature \(\beta^*\) of the teacher, so it cannot match them with its own. In particular, we consider the case where the noise injected by the teacher in the data is relatively small, i.e. \(\beta^* / N > 0\) as \(N \rightarrow \infty\), while the student chooses a conservative inverse temperature \(\beta \ll N\) to avoid overfitting. Moreover, we fix \(\sum_{y = 0}^C p^{* 0}_y = p_{\mathbf{h}}^* \left( 0 \right) = 0\) and \(\sum_{y = 0}^C p^{* \mu}_y = p_{\mathbf{h}}^* \left( \mu \right) = 1/P^*\) for all \(\mu > 0\) so that each \(\mathbf{g}^{* \mu} := \mathop{\mathrm{P}}_\beta \left( y | \mu ; \mathbf{w}, \mathbf{p} \right) = \mathbf{p}^{* \mu} / p_{\mathbf{h}}^* \left( \mu \right) = P^* \mathbf{p}^{* \mu}\) is a soft label for the corresponding \(\mathbf{w}^{* \mu}\). On the contrary, we do not give \(\sum_{\gamma = 0}^{P^*} p^{* \gamma}_y = p_{\mathbf{q}}^* \left( y \right)\) a restrictive form. In other words, \(p_{\mathbf{q}}^* \left( y \right)\) is free to be any given probability mass function.

Table (2) of Appendix 4.10 summarizes, for quick reference, the most important symbols introduced in this Section and their meaning.

4.3 Theoretical results↩︎

4.3.1 Saddle-point equations↩︎

In this Section, we introduce a set of equations for the stationary points of the loss (Eq. 50 ), which we then relate to the saddle-point equations emerging from the statistical mechanics analysis of posterior sampling (see Eq. 51 ) in the teacher-student setting.

Let us first establish a few definitions that we will use frequently throughout the Section. Given two matrices \(\mathbf{w}^*\in\mathbb{R}^{P^*\times N}\) and \(\mathbf{w}\in\mathbb{R}^{P\times N}\), we define the overlap matrix \(\mathbf{m}(\mathbf{w}^*,\mathbf{w})=\mathbf{w}^*\mathbf{w}^T \in\mathbb{R}^{P^*\times P}\). We write its entries as \(m^{\mu_* \mu} \left( \mathbf{w}^*, \mathbf{w} \right) = \sum_{i = 1}^N w^{* \mu_*}_i w^\mu_i\) and its row vectors as \(m^{\mu_*}\left( \mathbf{w}^*, \mathbf{w} \right)\), where \(1 \leq \mu_* \leq P^*\) and \(1 \leq \mu \leq P\). Moreover, for any matrix \(\mathbf{m}\in \mathbb{R}^{P^*\times (P+1)}\) with entries \(m^{\mu_* \gamma}\) and row vectors \(m^{\mu_*}\), we use \[\sigma_\gamma(m^{\mu_*})= \frac{\exp \left( m^{\mu_* \gamma} \right)}{\sum_{\nu = 0}^P \exp \left( m^{\mu_* \nu} \right)}\] to represent the entry number \(\gamma \in \left\{ 0, 1, ..., P \right\}\) of the softmax function applied to the row vector \(m^{\mu_*}\). In this context, \(m^{\mu_*}\) has a zeroth component \(m^{\mu_*0}\), and so does its softmax.

In Appendix 4.11, we show that the stationary points of the negative log-likelihood loss (Eq. 50 ) satisfy \[\begin{align} \label{eq:stationarity} w^\mu_i &= \frac{\bar{w}^\mu_i}{\sqrt{\sum_{j = 1}^N \left[\bar{w}^\mu_j \right]^2}} \\ p^\gamma_y &= \frac{\bar{p}^\gamma_y}{\zeta^\gamma_y \left( \mathbf{\bar{p}} ; p_{\mathbf{h}} \right)} \addtocounter{equation}{1}\\ \text{with} \quad \bar{w}^{ \mu}_i &= \sum_{\mu_* = 1}^{P^*} x^{* \mu_*}_i \sum_{y = 0}^C q^{* \mu_*}_y \sigma_\mu \left( \beta m^{\mu_*} \left( \mathbf{x}^*, \mathbf{w} \right) + \log \left[ \mathbf{p}_y \right] \right) \\ \bar{p}^\gamma_y &= \sum_{\mu_* = 1}^{P^*} q^{* \mu_*}_y \sigma_\gamma \left( \beta m^{\mu_*} \left( \mathbf{x}^*, \mathbf{w} \right) + \log \left[ \mathbf{p}_y \right] \right) \end{align}\tag{52}\] for all \(1 \leq \mu \leq P\) and \(0 \leq \gamma \leq P\), where \(m^{\mu_* 0} \left( \mathbf{x}^*, \mathbf{w} \right) = \frac{1}{\beta} \log \left[ \Omega_N \left( \beta \right) / \Omega_N \left( 0 \right) \right]\) and the normalization constant \(\zeta^\gamma_y \left( \mathbf{\bar{p}} ; p_{\mathbf{h}} \right)\) is defined in Appendix 4.8.

Similar equations arise naturally in our statistical mechanics analysis, which amounts to using the replica method nishimori2001statistical?, charbonneau2022replica? to compute the limiting free entropy \[\begin{align} \label{eq:free95entropy} f \left( \varrho, \upsilon, \beta, P^*, P \right) = \lim_{\beta^*, M,N\to\infty} \frac{1}{N} \mathbb{E}_{\mathbf{w}^*, \mathbf{p}^*, \mathcal{D}} \log \left[ \mathcal{Z}_\beta \left( \mathcal{D} \right) \right] , \end{align}\tag{53}\] where \(\upsilon = \beta^* / N\), \(\varrho = \frac{M}{P^* N}\) and \(\mathbb{E}_{\mathbf{w}^*, \mathbf{p}^*, \mathcal{D}}\) is the joint expectation over the distribution of examples \(\mathcal{D} = \left\{ \mathbf{x}^c, y^c \right\}_{c = 1}^M\) generated by the teacher and the priors on the teacher weights. To be more precise, we show that the free entropy can be computed with a variational principle of the form \[\begin{align} f \left( \varrho, \upsilon, \beta, P^*, P \right) = \operatorname{Extr}_{\mathbf{m},\hat{\mathbf{m}},\mathbf{p}} f(\mathbf{m},\hat{\mathbf{m}},\mathbf{p}), \end{align}\] whose extremizer \(\mathbf{m} = \left\{ m^{\mu_* \mu} \right\}_{1 \leq \mu \leq P}^{1 \leq \mu_* \leq P^*}\) can be interpreted as the \(N \rightarrow \infty\) limit of the expected value of the teacher-student overlaps \(\mathbf{m}(\mathbf{w}^*,\mathbf{w})\). We show the derivation of the variational principle in Appendix 4.12 together with an explicit expression for the trial function \(f(\mathbf{m},\hat{\mathbf{m}},\mathbf{p})\). If we assume that there are \(P^* \ll N\) teacher memories and that the priors on \(\mathbf{w}^{* \mu_*}\) and \(\mathbf{g}^*\) are uniform like those of the student, we find that the expected teacher-student overlaps must satisfy the saddle-point equations \[\begin{align} \label{eq:uniform95saddle-point} m^{\mu_* \mu} &= \varsigma \left( 2 \beta_{\text{eff}} \varrho \sqrt{\sum_{\nu_* = 1}^{P^*} \left[ \hat{m}^{\nu_* \mu} \right]^2} \right) \frac{\hat{m}^{\mu_* \mu}}{\sqrt{\sum_{\nu_* = 1}^{P^*} \left[ \hat{m}^{\nu_* \mu} \right]^2}} \\ p^\gamma_y &= \frac{\bar{p}^\gamma_y}{\zeta^\gamma_y \left( \mathbf{\bar{p}} ; p_{\mathbf{h}} \right)} \addtocounter{equation}{1}\\ \text{with} \quad \hat{m}^{\mu_* \mu} &= \sum_{y = 0}^C p_{\mathbf{q}}^* \left( y \right) \sigma_\mu \left( \beta_{\text{eff}} m^{\mu_*} + \log \left[ \mathbf{p}_y \right] \right) \\ \bar{p}^\gamma_y &= p_{\mathbf{q}}^* \left( y \right) \sum_{\mu_* = 1}^{P^*} \sigma_\gamma \left( \beta_{\text{eff}} m^{\mu_*} + \log \left[ \mathbf{p}_y \right] \right) \end{align}\tag{54}\] for all \(1 \leq \mu \leq P\) and \(0 \leq \gamma \leq P\), where \(\varsigma \left( x \right) = \frac{x}{\sqrt{x^2 + 1} + 1}\), \(\beta_{\text{eff}} = \varsigma \left( 2 \upsilon \right) \beta\) and \(m^{\mu_* 0} = \frac{1}{\beta_{\text{eff}}} \log \left[ \Omega_N \left( \beta \right) / \Omega_N \left( 0 \right) \right]\).

If we instead clamp the teacher weights \(\mathbf{w}^{* \mu_*}\) and \(\mathbf{g}^{* \mu_*}\) to fixed patterns \(\mathbf{x}^{* \mu_*}\) and their corresponding soft labels \(\mathbf{q}^{* \mu_*}\) Szegedy2016rethinking?, respectively, to mimic a more general distribution for the data, then Eqs. (54 ) become (see Appendix 4.13) \[\begin{align} \label{eq:saddle-point} m^{\mu_* \mu} &= \varsigma \left( 2 \beta_{\text{eff}} \varrho \sqrt{\sum_{i = 1}^N \left[ \bar{x}^\mu_i \right]^2} \right) \frac{\sum_{i = 1}^N x^{* \mu_*}_i \bar{x}^\mu_i}{\sqrt{\sum_{i = 1}^N \left[ \bar{x}^\mu_i \right]^2}} \\ p^\gamma_y &= \frac{\bar{p}^\gamma_y}{\zeta^\gamma_y \left( \mathbf{\bar{p}} ; p_{\mathbf{h}} \right)} \addtocounter{equation}{1}\\ \text{with} \quad \bar{x}^\mu_i &= \sum_{\mu_* = 1}^{P^*} x^{* \mu_*}_i \sum_{y = 0}^C q^{* \mu_*}_y \sigma_{\mu} \left( \beta_{\text{eff}} m^{\mu_*} + \log \left[ \mathbf{p}_y \right] \right) \\ \bar{p}^\gamma_y &= \sum_{\mu_* = 1}^{P^*} q^{* \mu_*}_y \sigma_\gamma \left( \beta_{\text{eff}} m^{\mu_*} + \log \left[ \mathbf{p}_y \right] \right). \end{align}\tag{55}\] for all \(1 \leq \mu \leq P\) and \(0 \leq \gamma \leq P\), where we recall that \(m^{\mu_* 0} = \frac{1}{\beta_{\text{eff}}} \log \left[ \Omega_N \left( \beta \right) / \Omega_N \left( 0 \right) \right]\).

In the limit of \(\varrho, \upsilon \rightarrow \infty\), which indicates a large number of examples and a low level of teacher noise, Eqs. (55 ) become equivalent to Eqs. (52 ) if we make the identification \(\bar{x}^\mu_i=\bar{w}^\mu_i\), i.e. \[\label{eq:identification} w^\mu_i = \frac{\bar{x}^\mu_i}{\sqrt{\sum_{j = 1}^N \left[ \bar{x}^\mu_j \right]^2}}.\tag{56}\] Let us explain this limit step by step. For any \(\upsilon\) and \(\varrho\), the \(M\) examples \(\{\mathbf{x}^c\}_{c=1}^M\) given to the student are corrupted versions of the teacher patterns \(\mathbf{x}^{*}\) with a noise level of \(1/\upsilon\). Therefore, in the limit of \(\upsilon \rightarrow \infty\) with finite \(\varrho\), each example \(\mathbf{x}^c\) is one of the original patterns \(\mathbf{x}^*\), as in Eqs. (52 ). However, even in this case, the empirical distribution of the examples deviates from that of \(\mathbf{x}^*\), of which it is merely a bootstrap sample. At fixed \(P^*\) and \(P\), this mismatch disturbs the accurate learning of \(\mathbf{x}^*\) when \(\varrho = \frac{M}{P^* N}\) is finite, or equivalently the expected number of repetitions \(M / P^*\) of each pattern is not sufficiently large compared to \(N\) (see Eq. 55 ). This influence progressively weakens as \(\varrho\) grows and the \(\varsigma\) function approaches \(1\), reflecting the convergence of the empirical distribution of the examples to that of the teacher patterns.

In the limit of \(\varrho \rightarrow \infty\) with finite \(\upsilon\), Eq. (55 ) is still very similar to Eq. (52 ). In that case, the only difference between them is that Eq. (55 ) has \(\beta_{\text{eff}}\) instead of \(\beta\) in the argument of \(\sigma\) and in the denominator of the definition of \(m^{\mu_* 0}\). Using \(\varsigma\) as a shorthand for \(\varsigma \left( 2 \upsilon \right)\), we find that the fixed points of Eqs. (55 ) with \(\varrho \rightarrow \infty\) are related, through the identification made in Eq. (56 ), to the stationary points of the effective loss \[\begin{gather} \label{eq:effective95loss} \mathcal{L} \left( \mathbf{w}, \mathbf{p} \right) = -\frac{1}{P^*} \sum_{\mu = 1}^{P^*} \sum_{y = 0}^C q^{* \mu}_y \log \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu}, y | \mathbf{w}, \mathbf{p} \right), \addtocounter{equation}{1}\\ \text{where} \quad \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right) = \sum_{\mu = 1}^P p^\mu_y \frac{\exp \left( \varsigma \beta \sum_{i = 1}^N w^\mu_i x_i \right)}{ \Omega_N \left( \beta \right)} + p^0_y \frac{1}{\Omega_N \left( 0 \right)}. \end{gather}\tag{57}\] What distinguishes this equation from the standard loss (Eq. 50 ) is that \(\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right)\) has \(\beta_{\text{eff}} = \varsigma \beta\) in the argument of the exponential function instead of \(\beta\). As a consequence, \(\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right)\) is not a probability distribution unless \(\varsigma = 1\), in the limit of \(\upsilon \rightarrow \infty\) (see Appendix 4.7). A value of \(\varsigma\) less than one in the effective loss (Eq. 57 ) is reminiscent of the presence of noise in the data generation process, so we propose to use it as a regularizer for the weights. We discuss this point in more detail in Section 4.4.1.

4.3.2 Saddle-point hierarchy↩︎

As shown in zhang2021embedding?, the loss landscape of any NN with unconstrained weights contains the stationary points of narrower NNs with the same architecture. In Appendix 4.14, we show that this result also applies to the teacher-student setting with \(\varrho \rightarrow \infty\) and any nonzero \(\upsilon\). To be more precise, we show that, if the parameters \(\bar{x}^{\text{fixed}, \mu}_i\), \(\bar{p}^{\text{fixed}, \gamma}_y\), \(m^{\text{fixed}, \mu_* \gamma}\), \(p^{\text{fixed}, \gamma}_y\) with hidden unit prior \(p_{\mathbf{h}}^{\text{given}} \left( \gamma \right)\) are a fixed point of Eqs. (55 ) with \(P\) hidden units, then the duplicated parameters \[\begin{align} \label{eq:fixed95point} \bar{x}^{\text{dupli}, \mu}_i &= \begin{cases} \bar{x}^{\text{fixed}, \mu}_i &\quad 0 < \mu \leq P \\ \bar{x}^{\text{fixed}, \mu - P}_i &\quad P < \mu \leq P + R \\ \end{cases} \\ \bar{p}^{\text{dupli}, \gamma}_y &= \begin{cases} \bar{p}^{ \text{fixed}, 0}_y &\quad \gamma = 0 \\ \frac{1}{2} \bar{p}^{ \text{fixed}, \gamma}_y &\quad 0 < \gamma \leq R \\ \bar{p}^{ \text{fixed}, \gamma}_y &\quad R < \gamma \leq P \\ \frac{1}{2} \bar{p}^{ \text{fixed}, \gamma - P}_y &\quad P < \gamma \leq P + R \\ \end{cases} \\ m^{\text{dupli}, \mu_* \gamma} &= \begin{cases} m^{\text{fixed}, \mu_* 0} &\quad \gamma = 0 \\ m^{\text{fixed}, \mu_* \gamma} &\quad 0 < \gamma \leq P \addtocounter{equation}{1}\\ m^{\text{fixed}, \mu_*, \gamma - P} &\quad P < \gamma \leq P + R \\ \end{cases} \\ p^{\text{dupli}, \gamma}_y &= \begin{cases} p^{\text{fixed}, 0}_y &\quad \gamma = 0 \\ \frac{1}{2} p^{\text{fixed}, \gamma}_y &\quad 0 < \gamma \leq R \\ p^{\text{fixed}, \gamma}_y &\quad R < \gamma \leq P \\ \frac{1}{2} p^{\text{fixed}, \gamma - P}_y &\quad P < \gamma \leq P + R \\ \end{cases} \end{align}\tag{58}\] \[\begin{align} \text{along with} \quad p_{\mathbf{h}} \left( \gamma \right) &= \begin{cases} p_{\mathbf{h}}^{\text{given}} \left( 0 \right) &\quad \gamma = 0 \\ \frac{1}{2} p_{\mathbf{h}}^{\text{given}} \left( \gamma \right) &\quad 0 < \gamma \leq R \\ p_{\mathbf{h}}^{\text{given}} \left( \gamma \right) &\quad R < \gamma \leq P \\ \frac{1}{2} p_{\mathbf{h}}^{\text{given}} \left( \gamma - P \right) &\quad P < \gamma \leq P + R, \end{cases} \end{align}\]are a fixed point of the same saddle-point equations with \(P + R \in \left\{ P, ..., 2 P \right\}\) hidden units (See Appendix 4.14 for a detailed derivation and Fig. 41 for a illustrative example). In other words, we duplicate some of the weights solving Eqs. (55 ) to construct a fixed point for a wider network. In that sense, wide DAMs contain the fixed points of narrower DAMs. In particular (see Section 4.3.1), this property also holds for the stationary points of both the standard loss (Eq. 50 ) and the effective loss (Eq. 57 ).

The saddle-point equations (Eq. 55 ) with duplicated order parameters (Eqs. 58 ) are invariant to the permutation of any hidden unit \(\gamma \in \left\{ 1, ..., R \right\}\) and its duplicate \(\gamma + P\). This kind of symmetry can be spontaneously broken if it leads to a higher free entropy (see Eqs. 53 and 71 ). However, symmetry-breaking transitions can be prohibitively slow when the symmetric state is stable to local perturbations zdeborova2016statistical?. Interestingly, the DAM introduced in krotov2016dense? quickly undergoes many successive permutation symmetry-breaking bifurcations during training boukacem2024waddington?. This observation and additional empirical evidence suggest that the symmetric states are unstable, or in other words that they are saddles, which was verified analytically for DAMs with only two data points to memorize boukacem2024waddington?. In Appendix 4.14, we prove that, if \(\beta\) is large enough, Eqs. (58 ) is a saddle, which is a major step toward explaining why permutation symmetry breaking is relatively fast in DAMs trained with a large number of data points. We call this result the saddle-point hierarchy principle.

4.4 Empirical results↩︎

In this Section, we use our theoretical results to improve training, and we show empirically that our DAM learns interpretable solutions to both supervised and unsupervised classification problems. Unless explicitly stated otherwise, we perform our numerical experiments on the MNIST dataset of handwritten digits lecun1998gradient?. The code and hyperparameter values used in our numerical experiments are available in the following public repository theriault2025saddlesoftware?.

4.4.1 Learning by minimizing the effective loss↩︎

As the number of hidden units \(P\) increases, standard maximum likelihood estimation with fixed inverse temperature \(\beta\) becomes progressively less apt to train our DAM to its full potential. At high \(P\) and \(\beta\), many memories stay stuck in noisy states that do not contribute to classification, which is wasteful (see Fig. 42). Using a lower \(\beta\) helps the memories converge, but also reduces their resolution and diversity, possibly because it also lowers the DAM’s capacity lucibello2024exponential?. This change is far from being only cosmetic. In fact, it comes with a gradual reduction in classification accuracy (see Table 1).

Optimization algorithms with a parameter analogous to \(\beta\) often converge to better solutions when this parameter is increased from a small value during optimization. For example, annealing schedules kirkpatrick1983optimization?, rose1990statistical? can be incorporated into minimization algorithms to help them find deeper local minima than they could otherwise. This observation suggests that our DAM would learn better memories if \(\beta\) increased during training. It is tempting to do so by SGD of Eq. (50 ) with respect to \(\beta\), but it makes \(\beta\) increase so quickly that many memories still stay stuck in noisy states.

We find that optimizing the effective loss (Eq. 57 ) with a suitable \(\varsigma\) slows down the evolution of \(\beta\) enough to let the DAM learn much cleaner memories (see Fig. 43, top panel) and achieve much better classification accuracy (see Table 1) than with standard training (Eq. 50 ). Looking back at Section 4.3.1, we propose to interpret \(\varsigma\) as a regularization parameter that helps the DAM take noise from the data into account during training. There is no obvious theoretically motivated way to find the best value of \(\varsigma\) for a generic dataset, so we choose it by hand. Despite this limitation, we believe that the simplicity and interpretability of our method still make it an interesting alternative to annealing schedules.

When we use this training method, we compute the DAM classification accuracy from the effective predictions \[\begin{gather} \label{eq:effective95preds} \mathcal{P}_{\beta, \varsigma} \left( y \big| \mathbf{x} ; \mathbf{w}, \mathbf{p} \right) = \frac{\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right)}{\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x} \big| \mathbf{w}, \mathbf{p} \right)} \addtocounter{equation}{1}\\ \text{where} \quad \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x} \big| \mathbf{w}, \mathbf{p} \right) = \sum_{y = 0}^C \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right), \end{gather}\tag{59}\] instead of the true predictions \(\mathop{\mathrm{P}}_{\beta} \left( y \big| \mathbf{x} ; \mathbf{w}, \mathbf{p} \right) = \mathop{\mathrm{P}}_{\beta} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right) / \mathop{\mathrm{P}}_{\beta} \left( \mathbf{x} \big| \mathbf{w}, \mathbf{p} \right)\). This approach allows us to calculate both the accuracy and the loss through a single evaluation of \(\mathcal{P}_{\beta, \varsigma} \left( y \big| \mathbf{x} ; \mathbf{w}, \mathbf{p} \right)\), which is more efficient than computing \(\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right)\) and \(\mathop{\mathrm{P}}_{\beta} \left( \mathbf{x}, y \big| \mathbf{w}, \mathbf{p} \right)\) separately when monitoring the progress of training.

4.4.2 Dense associative memory is interpretable, even in unsupervised classification↩︎

Now that we understand how to train our DAM reasonably well, we investigate one of the most interesting properties of the solutions that it learns: their interpretability. As advertised in the Introduction, we will explain how to further improve training in the following Section. We already mentioned that the regularization parameter \(\varsigma\) of the effective loss is interpretable (see Section 4.4.1 and Eq. 57 ). Here, we point out that the learned weights are as well. In fact, each learned \(\mathbf{p}^\mu / p_{\mathbf{h}} \left( \mu \right)\) can be interpreted as a soft label for the corresponding \(\mathbf{w}^\mu\) (see Fig. 43). This property was also observed by boukacem2024waddington? in K & H’s DAM for pattern classification krotov2016dense?.

At test time, we observe that our DAM classifies approximately 98% of the test data points into the class \(y = \mathop{\mathrm{argmax}}_{y^\prime} \left\{ p^\mu_{y^\prime} \right\}\) of the memory \(\mathbf{w}^\mu\) to which they are the most similar. In other words, a 1-nearest neighbor classifier fix1989discriminatory? conditioned on the memories \(\mathbf{w}^\mu\) and their soft labels \(\mathbf{p}^\mu / p_{\mathbf{h}} \left( \mu \right)\) approximates the classification of our model with 98% fidelity. This behavior is reminiscent of K & H’s DAM krotov2016dense?, where only a few memories participate in the classification of each data point. The 1-nearest neighbor classifier approximates DAM classification more faithfully for correctly classified data points (approximately 99% fidelity) than incorrectly classified data points (approximately 70% fidelity).

We now show that the notion of interpretability described in this Section extends beyond supervised learning. Given a dataset of \(P^*\) unlabeled patterns \(\{\mathbf{x}^{* \mu}\}_{\mu=1}^{P^*}\), we find that we can train our DAM for unsupervised classification by replacing the soft labels \(q^{* \mu}_y\) in the effective loss (Eq. 57 ) with the softened DAM predictions \(\left( 1 - \varepsilon \right) \mathcal{P}_{\beta, \varsigma} \left( y | \mathbf{x}^{* \mu} ; \mathbf{w}, \mathbf{p} \right) + \varepsilon \frac{1}{C + 1}\), where \(\varepsilon \in \left[ 0, 1 \right]\) and \(\mathcal{P}_{\beta, \varsigma} \left( y | \mathbf{x}^{* \mu} ; \mathbf{w}, \mathbf{p} \right)\) is defined in Eq. (59 ). In this scenario, we are minimizing \[\label{eq:unsupervised95loss} \mathcal{L}_{\text{unsup}} \left( \mathbf{w}, \mathbf{p} \right) = -\frac{1}{P^*} \sum_{\mu = 1}^{P^*} \sum_{y = 0}^C \left[ \left( 1 - \varepsilon \right) \mathcal{P}_{\beta, \varsigma} \left( y | \mathbf{x}^{* \mu} ; \mathbf{w}, \mathbf{p} \right) + \varepsilon \frac{1}{C + 1} \right] \log \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu}, y | \mathbf{w}, \mathbf{p} \right).\tag{60}\] As before, we do so using constrained SGD (see Appendix 4.8) with the initialization and learning rate described in Appendix 4.9. Equivalently, we can also view this algorithm as minimizing the combined loss \[\begin{gather} \label{eq:total95loss} \mathcal{L}_{\text{total}} \left( \mathbf{w}, \mathbf{p} \right) = \mathcal{L}_{\text{margin}} \left( \mathbf{w}, \mathbf{p} \right) + \lambda \mathcal{L}_{\text{cond}} \left( \mathbf{w}, \mathbf{p} \right), \addtocounter{equation}{1}\\ \text{where} \quad \mathcal{L}_{\text{margin}} \left( \mathbf{w}, \mathbf{p} \right) = -\frac{1}{P^*} \sum_{\mu = 1}^{P^*} \log \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu} | \mathbf{w}, \mathbf{p} \right) \\ \text{and} \quad \mathcal{L}_{\text{cond}} \left( \mathbf{w}, \mathbf{p} \right) = -\frac{1}{P^*} \sum_{\mu = 1}^{P^*} \sum_{y = 0}^C \mathcal{P}_{\beta, \varsigma} \left( y | \mathbf{x}^{* \mu} ; \mathbf{w}, \mathbf{p} \right) \log \mathcal{P}_{\beta, \varsigma} \left( y | \mathbf{x}^{* \mu} ; \mathbf{w}, \mathbf{p} \right), \end{gather}\tag{61}\] where \(\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu} | \mathbf{w}, \mathbf{p} \right) = \sum_{y = 0}^C \mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu}, y | \mathbf{w}, \mathbf{p} \right)\) (see Eq. 59 ). \(\mathcal{L}_{\text{cond}} \left( \mathbf{w}, \mathbf{p} \right)\) is called the minimum entropy regularization term in unsupervised machine learning grandvalet2004semi?. Intuitively, it encourages the DAM to learn different class weights \(p^\gamma_y\) for each class \(y\), which is not possible by minimizing only \(\mathcal{L}_{\text{margin}} \left( \mathbf{w}, \mathbf{p} \right)\) because \(\mathcal{P}_{\beta, \varsigma} \left( \mathbf{x}^{* \mu} | \mathbf{w}, \mathbf{p} \right)\) does not depend on \(y\). Exploiting this characteristic, we train our DAM on patches of MNIST digits lecun1998gradient? and find that it learns reasonable latent classes \(y = \mathop{\mathrm{argmax}}_{y^\prime} \left\{ p^\mu_{y^\prime} \right\}\) for the memories \(\mathbf{w}^\mu\) (see Fig. 44 and Figs. 49, 50, 51 and 52 of Appendix 4.15). This approach could potentially be useful for feature extraction khalid2014survey?.

4.4.3 Fast training with splitting steepest descent↩︎

We now explain how to further improve the training method of Section 4.4.1. More specifically, we use the saddle-point hierarchy derived in 4.3.2 to accelerate training.

Machine learning models that experience permutation symmetry breaking as described in Section 4.3.2—such as K & H’s DAM boukacem2024waddington?, RBMs with binary units hou2019minimal?, theriault2025modeling? and Gaussian mixtures kappen1993using?, KAPPEN1995deterministic?, rose1990statistical? —have characteristic tree-shaped learning curves. Intuitively, permutation symmetry-breaking transitions are points where model parameters differentiate from each other, so they correspond to the bifurcations of the tree. In the left panel of Fig. (45), we show that the learning dynamics of our DAM also follows a tree of permutation symmetry-breaking transitions.

The saddle-point hierarchy principle introduced in Section 4.3.2 suggests the following idea to accelerate learning: train a relatively narrow DAM for cheap, then repeatedly duplicate (or “split”) some of the hidden units, escape the corresponding saddle point and continue training. Intuitively, we expect to save a lot of computing resources by using relatively few hidden units at the start of training, when the learning dynamics follows a few branches close to the root of the learning dynamics tree. The splitting steepest descent algorithm introduced in wu2019splitting? formalizes this idea in an efficient and theoretically motivated way. We implement a variant of this algorithm that takes the constraint \(\mathbf{w}^\mu \in S^{N - 1}\) into account, but is otherwise very similar to the original version (see Alg. 46), with fast splitting implemented as in wang2019energy?. The constraint \(\mathbf{w}^\mu \in S^{N - 1}\) only matters in steps 5 and 11, which are also arguably the most conceptually difficult parts of Alg. (46), so we explain them in detail in Appendix 4.16. Fig. (45) shows that the learning dynamics tree of splitting steepest descent has a more sparsely populated trunk than that of SGD without splitting (see Fig. 45). In other words, by using a relatively small number \(P\) of memories at the beginning of training, we reduce the total number of values \(P \times T\) that they take over a period of time \(T\), which is consistent with our intuition that it can save computational resources. The situation could have been different if using fewer memories slowed down training.

We now establish a methodology that we will use to compare the training times of DAMs trained on MNIST lecun1998gradient? and Fashion-MNIST xiao2017fashion? using the effective loss (Eq. 57 ) with and without splitting steepest descent. It is more interesting and meaningful to compare the training times of NNs with similar performance. Therefore, we pick hyperparameters such that DAMs trained with and without splitting have similar classification accuracy. For that purpose, we find the general region of hyperparameter space where DAMs trained without splitting steepest descent have the best classification accuracy. The accuracy does not change much in this region, so we pick generic hyperparameters inside. Next, we manually tune the hyperparameters of splitting steepest descent so that the resulting DAMs have comparable accuracy. The accuracy obtained for Fashion-MNIST is generically slightly higher with splitting steepest descent than without it. In contrast, splitting steepest descent does not significantly affect the accuracy obtained for MNIST.

Once the hyperparameters are set, we collect statistics of the accuracy and training time. We run our experiment on a CPU and manually set the seed of pseudorandom number generation to make training deterministic and reproducible. Two DAMs with the same hyperparameters and the same seed are guaranteed to be trained using the same number of computer operations and to have the same classification accuracy after training. However, background processes and other external factors can change between two runs with the same seed, which adds noise to the bare training time that we want to measure. As such, we approximate the bare training time as the minimum over many runs with the same seed. We then calculate the average and standard deviation of the classification accuracy and bare training time over multiple seeds.

In Fig. (47), we compare DAM training time and accuracy with and without splitting steepest descent, using round markers and error bars to show their respective means and standard deviations as a function of the maximum number of hidden units \(P_{\text{max}}\). Splitting stops at \(P_{\text{max}}\) hidden units, and DAMs without splitting have \(P_{\text{max}}\) hidden units from the start. All hidden units are split in each while-loop iteration of Alg. (46), except for the last iteration, where exactly \(P_{\text{max}} - P_{\text{cur}}\) hidden units are split. As wanted, DAMs have similar accuracy with and without splitting (left panels). Moreover, the reasonably small error bars of the training time indicate that the residual noise of the background processes is controlled and that the seed has a limited impact on the training speed (right panels). Splitting steepest descent has a significant speed advantage that scales very advantageously with \(P_{\text{max}}\) (right panels). Without splitting, the training time is proportional to \(P_{\text{max}}\). With splitting, it appears to be piecewise constant. The clearest sign that it increases with \(P_{\text{max}}\) is the jump between \(P_{\text{max}} = 1500\) and \(P_{\text{max}} = 1625\), where the number of while-loop iterations in Alg. (46) increases by 1. In other words, the run time of Alg. (46) is better explained as proportional to \(\left\lceil \log P_{\text{max}} \right\rceil\) than \(P_{\text{max}}\) in the range of \(P_{\text{max}}\) of our numerical experiment, which is a dramatic improvement over the run time without splitting. The benefit of the algorithm becomes evident primarily when the accuracy improves consistently across a wide range of hidden unit counts, i.e. when training a substantially large network is beneficial. This pattern, for example, does not apply to Fashion MNIST (see Fig. 47, bottom plot), where the accuracy gain from training a larger network is marginal. As a result, the need for strategies that accelerate training becomes less critical, although such acceleration still occurs.

4.5 Conclusion↩︎

In this work, we study a Dense Associative Memory model based on the framework of Boltzmann machines that is capable of learning interpretable solutions to both supervised and unsupervised problems.

We derive two sets of equations that respectively characterize the stationary points of DAMs trained on real data using maximum likelihood estimation and the fixed points of DAMs trained on synthetic data in the teacher-student setting. Guided by their similarity, we then establish that the maximum likelihood equations are a special case of the teacher-student equations. Building on this equivalence, we introduce an effective loss function that makes training significantly more stable thanks to a regularization parameter that mimics the noise present in the data generation process.

We show that the stationary points of the effective loss of DAMs with a given number of hidden units are saddle points of larger DAMs, and we considerably accelerate training using a custom implementation of the splitting steepest descent network-growing algorithm wu2019splitting? inspired by this saddle-point hierarchy principle. In our numerical experiments, the training time of splitting steepest descent is approximately logarithmic in the number of hidden units, which is a huge improvement over the linear training time that we observe without splitting steepest descent. Further research is needed to clearly determine the asymptotic functional form of the training time.

In Section 4.4.1, we briefly discuss how our DAM is interpretable. The ingredients needed to make an interpretable neural network are not well understood, and we believe that our DAM offers a promising opportunity to explore that area. Studying the teacher-student setting with teacher inverse temperature \(\beta^*\) less than order \(N\) is another interesting research direction that we leave to future work. Beyond the framework of our model, rose1993constrained?, rose1998deterministic? have designed a technique that trains Gaussian mixture models by replicating mixture components. In contrast with our work, it uses a handpicked cooling schedule and breaks permutation symmetries with random noise. By design, splitting steepest descent finds the directions in the loss landscape that allow NNs to escape saddle points with permutation symmetries as quickly as possible, but finding these directions (as described in Appendix 4.16 and wu2019splitting?, wang2019energy?) is more expensive than generating random noise. It would be interesting to study this tradeoff and identify the regimes where each method is preferable over the other. More generally, one could investigate how permutation symmetry breaking takes place in energy models related to DAMs, such as attentional BMs ota2023attention?, and study its relationship with the permutation symmetry-breaking hou2019minimal?, theriault2025modeling? and dynamical transitions decelle2017spectral?, decelle2018thermodynamics?, bachtis2025cascade? found in other types of RBMs, whose careful analysis was also used to improve training bereux2025fast?. Recent works boukacem2024waddington?, grishechkin2025mathematical? have made connections between DAM learning dynamics and cellular differentiation. Similarly, it could be possible to investigate whether there are similarities between cellular division and splitting steepest descent in our model. Finally, it would be interesting to see if splitting steepest descent and our proposed regularization technique can be used to improve the training of transformers and generative diffusion.

Data availability↩︎

The code for the numerical experiments presented in this work is available at this public repository theriault2025saddlesoftware?.

Acknowledgements↩︎

I am grateful to Carlo Lucibello and Matteo Negri for insightful discussions.

Conclusion↩︎

In this thesis, we study various kinds of dense associative memory (DAM) models and restricted Boltzmann machines (RBMs) in the teacher-student setting.

In Chapter 2, we study a DAM called dense Hopfield network (dense HN). On the Nishimori line, we show that the phase transition where student dense HNs with a single pattern become capable of learning a dataset with \(M\) samples generated by a teacher dense HN with a single pattern coincides with the spin-glass transition of dense HNs with \(M\) random patterns. Outside the Nishimori line, we study the resistance of dense HNs to noise and adversarial attacks. In particular, we derive a formula quantifying the adversarial robustness of dense HNs at zero temperature, and we clarify why the adversarial robustness of dense HNs depends on the learning regime, as observed in krotov2018dense?.

In Chapter 3, we study RBMs with any finite number of hidden units. When the patterns incident to the \(P^*\) hidden units of the teacher RBM are uncorrelated, we show that a student RBM with \(P\) hidden units learns the data in the same way as \(P\) separate RBMs with one hidden unit each, thus validating a conjecture formulated in barra2017phase?. Moreover, we show that student RBMs with more hidden units than their teacher learn a representation of the data where exactly \(P^*\) of their patterns align themselves with those of the teacher while the \(P - P^*\) remaining ones freeze in spin-glass states. We then argue that such RBMs can be used as toy models to study the lottery ticket hypothesis frankle2018lottery?. When the teacher patterns are correlated, we show that the student can adopt different learning strategies depending on the hyperparameters of the teacher-student setting and the number of samples in the training dataset.

In Chapter 4, we study a kind of DAM that fits into the framework of RBMs and is capable of both supervised and unsupervised classification. Based on insights from our theoretical analysis of the teacher-student setting, we propose a novel regularized loss function that makes training significantly more stable. Moreover, we show that the fixed points of relatively small DAMs are saddle points of larger DAMs, and we leverage this saddle-point hierarchy to considerably accelerate training using the splitting steepest descent algorithm introduced in wu2019splitting?.

A natural continuation of this work would be to study adversarial attacks in DAMs with many hidden units, such as the model studied in Chapter 4 theriault2026saddle?. In particular, it would be interesting to see if the methods that we use to study adversarial attacks in Chapter 2 can be generalized and used for this purpose in combination with those used to study adversarial attacks in linear models tanner2024high?, vilucchio2024geometry?, vilucchio2025existence?.

One could also extend this work by studying DAM and RBM training dynamics with analytical calculations, the statistical mechanics approach of decelle2017spectral?, decelle2018thermodynamics? being a possible starting point. First, it would be interesting to build a theoretical model of the training dynamics of the lottery ticket experiment presented in Chapter 3 theriault2025modeling?, as it could provide crucial insight into the lottery ticket hypothesis. In particular, it could help us to understand if the initial conditions that make certain subnetworks of large NNs converge especially quickly have common properties that can be exploited to design novel initialization schemes leading to faster training. Second, it would also be interesting to study analytically the learning dynamics of the DAM studied in Chapter 4 theriault2026saddle?, particularly near saddle points, for a rigorous comparison of the standard training time with that of splitting steepest descent.

To further our study of implicitly low-dimensional learning, it would be interesting to investigate models with low-rank weights, such as in decelle2017spectral?, decelle2018thermodynamics?. For example, one could study the teacher-student setting in which a student RBM (or DAM) with weights constrained to have a low rank learns data generated by a teacher whose weights are structured in a similar way for representing data lying on a low-dimensional manifold goldt2020modelling?, gerace2021generalization?. This kind of study could be used to investigate the benefits of low-rank constraints in exploiting the structure of data to improve the speed and quality of learning.

Finally, the interpretable learning dynamics decelle2017spectral?, decelle2018thermodynamics?, boukacem2024waddington? of RBMs and DAMs make them promising surrogate models of biological systems. It was notably shown that DAM learning dynamics is similar to cellular differentiation boukacem2024waddington?, and it would be interesting to see if the splitting steepest descent algorithm used to train DAMs in Chapter 4 theriault2026saddle? can be related to cellular division in an analogous way.

Funding information↩︎

This work was partially supported by project SERICS (PE00000014) under the MUR National Recovery and Resilience Plan funded by the European Union - NextGenerationEU. The work was also supported by the project PRIN22TANTARI "Statistical Mechanics of Learning Machines: from algorithmic and information-theoretical limits to new biologically inspired paradigms" 20229T9EAT – CUP J53D23003640001.

---

1. Except for a factor of two that can be removed by rescaling the inverse temperature.↩︎

2. See the end of Appendix C of the cited paper.↩︎