Non-linear PDEs approach to statistical mechanics of Dense Associative Memories

Artificial intelligence (AI) is rapidly changing the face of our society thanks to its impressive abilities in accomplishing complex tasks and extracting information from non-trivially structured, high-dimensional datasets. The successful applications of modern AI range from hard sciences and technology to more practical scenarios (e.g. medical sciences, economics and finance, and many daily tasks). Its important success is primarily due to the ascent of deep learning LeCun2015 ; Schmidhuber2015 , a set of semi-heuristic techniques consisting in stacking together several minimal building blocks in complex architectures with extremely-high learning performances. Despite its success, a rigorous theoretical framework guiding the development of such machine-learning architectures is still lacking. In this context, statistical mechanics of complex systems offers ideal tools to study neural network models from a more theoretical (and rigorous) point of view, thus drawing a feasible path which makes AI less empirical and more explainable.

In statistical mechanics, a primary classification of physical systems is the following. On the one side, we have simple systems, which are essentially characterized by the fact that the number of equilibrium configurations does not depend on the system size $N$. A paradigmatic (mean-field) realization of this situation is the Curie-Weiss (CW) model, in which all the spins $\sigma_{i}$, $i=1,\dots,N$, making up the system interact pairwisely by a constant, positive (i.e. ferromagnetic) coupling $J$. Below the critical temperature, in fact, the system exhibits ordered collective behaviors, and the equilibrium configurations of the system are characterized by only two possible values of the global magnetization $m(\boldsymbol{\sigma}):=\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}$ (which are further related by a spin-flip symmetry $\sigma_{i}\to-\sigma_{i}$ for each $i=1,\dots,N$). On the opposite side, we have complex systems, in which the number of equilibrium configurations increases according to an appropriate function of the system size $N$ due to the presence of frustrated interactions mezard1988spin . The prototypical example of mean-field spin-glass is the Sherrington-Kirkpatrick (SK) model sherrington1975 , in which the interaction strengths between the spin pairs are i.i.d. Gaussian variables. With respect to simple systems, spin-glass models exhibit a richer physical and mathematical structure, as shown by the presence of the spontaneous replica-symmetry breaking and an infinite number of phase transitions (e.g. see Parisi1979rsb1 ; Parisi1979rsb2 ; Parisi1980rsb3 ; Parisi1980sequence ; mezard1984replica ; guerra2003broken ; ghirlanda1998general ; talagrand2000rsb ) as well as the ultrametric organization of pure states (e.g. see panchenko2010ultra ; panchenko2011ultra ; panchenko2013parisi ). Statistical mechanics of spin glasses has acquired a prominent role during the last decades due to its ability to describe the equilibrium dynamics of several paradigmatic models for AI, in particular thanks to the work by Amit, Gutfreund, and Sompolinsky amit1985 for associative neural networks. For our concerns, the relevant ones are the Hopfield model hopfield1982hopfield ; pastur1977exactly and its $p$-spin extensions, the Dense Associative Memories (DAMs) Krotov2016DenseAM ; Krotov2018DAMS ; AD-EPJ2020 , exhibiting features which are peculiar both of ferromagnetic (simple) and spin-glass (complex) systems. In these models, the interactions between the spins are designed in order to store $K$ “information patterns”, each encoded by $K$ binary vectors of length $N$ and denoted by $\{\boldsymbol{\xi}^{\mu}\}_{\mu=1,\dots,K}$, with $\boldsymbol{\xi}^{\mu}=(\xi_{1}^{\mu},\xi_{2}^{\mu},...,\xi_{N}^{\mu})\in\{-1,+1\}^{N}$ and $\xi^{\mu}_{i}$ is a Rademacher random variable for any $i=1,\dots,N$ and $\mu=1,\dots,K$; the $\mu$-th pattern is said to be stored if the configuration $\boldsymbol{\sigma}=\boldsymbol{\xi}^{\mu}$ is an equilibrium state and the relaxation to this configuration, starting from a relatively close one (i.e. a corrupted version of $\boldsymbol{\xi}^{\mu}$), is interpreted as the retrieval of that pattern. The Hamilton function (or the energy in physics jargon) of these systems can be expressed as

where $p$ is the interaction order (for the Hopfield model $p=2$, while $p>2$ for the DAMs) and $m_{\mu}(\boldsymbol{\sigma}):=\frac{1}{N}\sum_{i=1}^{N}\xi^{\mu}_{i}\sigma_{i}$ is the so-called Mattis magnetizations measuring the retrieval of the $\mu$-th pattern. It has been shown that the number of storable patterns scales, at most, as a function of the system size, more precisely, $K<\alpha_{c}(p)N^{p-1}$, where $\alpha_{c}(p)\in\mathbb{R}$ depends on the interaction order $p$ and is referred to as critical storage capacity baldi1987number ; bovier2001spin . By a statistical-mechanics investigation of these models one can highlight the macroscopic observables (order parameters) useful to describe the overall behavior of the system, namely to assess whether it exhibits retrieval capabilities, and the natural control parameters whose tuning can qualitatively change the system behavior; such knowledge can then be summarized in phase diagram.

Regarding the methods, statistical mechanics offers a wide set of techniques for analyzing the equilibrium dynamics of complex systems, and in particular to solve for their free-energy. Historically, the first method (which was applied to the SK model and the Hopfield model amit1985 ; steffan1994replica ) is the replica trick, which – despite being straightforward and effective – is semi-heuristic and suffers from delicate points, see for example tanaka2007moment . Alternative, rigorous approaches were developed during the years and, among these, the relevant one for our concerns is Guerra’s interpolating framework. In this case, we can take advantage of rigorous mathematical methods by applying sum rules guerra2001sum or by mapping the relevant quantities (the free energy or the model order parameters) of the statistical setting to the solutions of PDE systems. Indeed, differential equations involving the partition functions (or related quantities) of thermodynamic models have been extensively investigated in the literature, see for example agliari2012notes ; barra2013mean ; barra2008mean ; barra2010replica ; barra2014proc ; Moro2018annals ; AMT-JSP2018 ; Moro2019PRE ; ABN-JMP2019 ; fachechi2021pde ; AAAF-JPA2021 . In particular, they allow us to express the equation of state (or the self-consistency equations) governing the equilibrium dynamics of the system in terms of solutions of non-linear differential equations, and to describe phase transition phenomena as the development of shock waves, thus linking critical behaviours to gradient catastrophe theory Barra2015Annals ; DeNittis2012PrsA ; Giglio2016Physica ; Moro2014Annals . In a recent work fachechi2021pde , a direct connection between the thermodynamics of ferromagnetic models with interactions of order $p$ and the equations of the Burgers hierarchy was established by linking the solution of the latter as the equilibrium solution of the order parameter of the former (i.e. the global magnetization $m$). In the present paper, we extend these results to complex models, in particular to the Hopfield model and the DAMs.

The paper is organised as follows. In Section 2, we introduce the relevant tools for our investigations, in particular Guerra’s interpolating scheme for the PDE duality. In Section 3, as a warm up, we review some basic results about the $p$-spin ferromagnetic models. In Section 4, we extend our results to the Derrida models (constituting the $p-$spin extension of the SK spin glass) derrida1981rem . In Section 5, we merge our results in a unified methodology for dealing with the DAMs, especially in the so-called high storage limit, and re-derive the self-consistency equations for the order parameters by means of PDE technology.

In this Section, we present the thermodynamic objects we aim to study. We start with a system made up of $N$ spins whose configurations $\boldsymbol{\sigma}\in\Sigma_{N}\equiv\{-1,+1\}^{N}$ are the nodes of a hypercube and interacting via a suitable tensor $\boldsymbol{J}$ of order $p$. The Hamilton functions of the system we will consider in this paper are of the form

where $D_{p,N,J}$ is a normalization factor ensuring the linear extensivity of the energy with the system size. Once the Hamiltonian is fixed, we introduce the partition function in the usual Boltzmann-Gibbs form. Thus, given $\beta\in\mathbb{R}_{+}$ the level of thermal noise of the system, the partition function is defined as

For simple systems, the partition function can be computed exactly for any $\boldsymbol{J}$ coupling matrix. As is standard in statistical mechanics, it is convenient to compute intensive quantities which are well-defined in the thermodynamic limit $N\to\infty$. Since the partition function is a sum of $2^{N}$ contribution, it is sufficient to take the intensive logarithm of the partition function, i.e.

which is the intensive statistical pressure (which is, apart for a factor $-\beta$, the usual free energy) of the system. However, when dealing with spin-glass systems, the coupling tensor $\boldsymbol{J}$ is a multidimensional random variable, thus the partition function defines a random measure on the configuration space. For good enough probability distributions of the coupling matrix, the intensive logarithm of the partition function is expected to converge to its expectation value in the thermodynamic limit $N\to\infty$ by virtue of self-averaging theorems ShcherbinaPastur-JSP1991 ; Bovier-JPA1994 , so it is natural to consider the quenched intensive pressure associated to the partition function (2), which is defined as

where $\mathbb{E}_{\boldsymbol{J}}$ denotes the average over the quenched disorder $\boldsymbol{J}$ (we stress that, in this case, the free energy does not depend any longer on the coupling matrix because of the average operation).

Rather than working with the quantity (2), we will use as a fundamental object Guerra’s interpolated partition function and its associated interpolating intensive pressure. For instance, for spin-glass systems we would have

where $H_{N,p,\boldsymbol{J}}(t,\boldsymbol{x})$ denotes the interpolating Hamiltonian satisfying the properties that, at $\boldsymbol{x}=0$ and $t\neq 0$, it recovers the Hamiltonian (1) times $\beta$, and at $t=0$ and $\boldsymbol{x}\neq\boldsymbol{0}$ it corresponds to an exactly-solvable $1$-body system, namely a system where spins interact only with an external field that has to be set a posteriori. The interpolating parameters $t$ and $\boldsymbol{x}$ are interpreted, in a mechanical analogy, as spacetime coordinates with suitable dimensionality.

The interpolating structure (4) implies a generalized measure, whose related Boltzmann factor is

Thus, for an arbitrary observable $O(\boldsymbol{\sigma})$ in the configuration space $\Sigma_{N}$, we can introduce the Boltzmann average induced by the partition function (4) as

Usually, in spin-glass systems, the quenched average is performed after taking the Boltzmann expectation values on the $s$-replicated space $\Sigma_{N}^{(s)}=(\Sigma_{N})^{\otimes s}\equiv\{-1,+1\}^{sN}$, which is naturally endowed with a random Gibbs measure corresponding to the partition function $Z_{N,p,\boldsymbol{J}}^{(s)}(t,\boldsymbol{x})=Z_{N,p,\boldsymbol{J}}(t,\boldsymbol{x})^{s}$. Given a function $O:\Sigma_{N}^{(s)}\to\mathbb{R}$, the Boltzmann average in the $s$-replicated space are straightforwardly defined as

where $\underline{\boldsymbol{\sigma}}\in\Sigma_{N}^{(s)}$ is the global configuration of the replicated system, and $B_{N,p,\boldsymbol{J}}^{(s)}(t,\boldsymbol{x})$ is the Boltzmann factor associated to the $s$-replicated partition function. Of course, in spin-glass theory, the relevant quantities are the quenched expectation values, which are defined as

For the sake of simplicity, we dropped the index $s$ from the quenched averages, as it would be clear from the context.

With all these definitions in mind, we are then able to find the link between the resolution of the statistical mechanics of a given spin-like model and a specific PDE problem in the fictitious space $(t,\boldsymbol{x})$. Before concluding this section it is worth recalling that here we will work under the replica-symmetry (RS) assumption, meaning that we assume the self-averaging property for any order parameter $X$, i.e. the fluctuations around their expectation values vanish in the thermodynamic limit. In distributional sense, this corresponds to

where $\bar{X}=\langle X\rangle_{t,\boldsymbol{x}}$ is the expectation value w.r.t. the interpolating measure $\mathcal{P}_{t,\boldsymbol{x}}(X)$. Typically, for simple systems this assumption is correct, conversely, for complex systems this is not always the case, for instance, in spin-glasses the RS is broken at low temperature mezard1988spin . When dealing with neural-network models, RS constitutes a standard working assumption as it usually applies (at least) in a limited region of the parameter space, while elsewhere it yields only small quantitative discrepancies with respect to the exact solution amit1989 ; SpecialIssue-JPA . The latter, accounting for RSB phenomena, can be obtained by iteratively perturbing the RS interpolation scheme (e.g. see AABO-JPA2020 ; AAAF-JPA2021 ; albanese2021 ), thus, our results find direct application on the practical side and provide the starting point for further refinements on the theoretical side.

The present section is a compendium of the results reported in fachechi2021pde , so we refer to that work for a detailed derivation. In $p$-spin ferromagnets, the interaction between spins is fixed by the requirement $J_{i_{1},i_{2},\dots,i_{p}}=J$ for each $i_{1},\dots,i_{p}=1,\dots,N$ and $J>0$; without loss of generality, one can set $J=1$, since it corresponds to a rescaling of the thermal noise. Thus, the Hamilton function of the model simply reads as

with

being the global magnetization of the system. By following the same lines of fachechi2021pde , Guerra’s interpolating partition function reads as

where $(t,x)\in\mathbb{R}^{2}$. The starting point is to notice that the interpolating statistical pressure associated to the partition function (10) has spacetime derivatives

The expectation value of monomials of the global magnetization satisfies the following relation fachechi2021pde :

This means that we can act on the expectation value $\omega_{t,x}(m(\boldsymbol{\sigma}))$ to generate higher momenta. In particular, calling $u(t,x)=\omega_{t,x}(m(\boldsymbol{\sigma}))$ and setting $s=p-1$, we directly get the Burgers hierarchy

This duality also allows us to analyse the thermodynamic limit, corresponding to the inviscid scenario for the Burgers hierarchy. Indeed, posing $\bar{u}(t,x)=\lim_{N\to\infty}u(t,x)=\lim_{N\to\infty}\omega_{t,x}(m(\boldsymbol{\sigma}))$, we have the initial value problem

where the initial profile is easily computed by straightforward calculations (since it is a 1-body problem). This system describes the propagation of non-linear waves, and can be solved by assuming a solution in implicit form $\bar{u}(t,x)=\tanh(x-v(t,x)t)$, where $v(t,x)=p\bar{u}(t,x)^{p-1}$ is the effective velocity. Recalling that the thermodynamics of the original $p$-spin model associated to the Hamilton function (9) is recovered by setting $t=-\beta$ and $x=0$, we directly obtain

where $\bar{m}=\bar{u}(-\beta,0)$. This is precisely the self-consistency equation for the global magnetization for the $p$-spin ferromagnetic model fachechi2021pde . The phase transition of the system is expected to take place where the gradient of the solution explodes, which, on the Burgers side, corresponds to the development of a shock wave at $x=0$. Since the temporal coordinate $t$ is directly related to the thermal noise at which the phase transition occurs, with standard PDE methods we can analytically determine the critical temperature according to the simple system

where $F(\xi)=p\tanh(\xi)^{p-1}$. This prediction is in perfect agreement with the numerical solutions of the self-consistency equation (17).

In this section, we adapt the previous methodologies to treat complex systems with $p$-spin interactions. The paradigmatic case is given by the $p$-spin SK model, also referred to as Derrida model, defined as follows

As stated in Sec. 1, when dealing with spin glasses we need to enlarge our analysis to the $s$-replicated version of the configuration space. To this aim, we use the following

In the following, in order to lighten the notation, the replica index $s$ of the Boltzmann average $\Omega^{(s)}_{t,x}$ can be dropped, since it is understood directly from the function to be averaged.

We can now focus on the PDE approach to the statistical mechanics of the $p$-spin SK model. To this aim, we compute the spacetime derivative of the quenched intensive pressure, as given in the following

The proof of this lemma can be found in Appendix A.

In order to proceed, we have now to make some physical assumptions on the model. As standard in spin-glass theory, the simplest requirement is the RS in the thermodynamic limit. In fact, as we are going to show, this makes the PDE approach feasible, due to the fact that we can express non-trivial expectation values of function of the replicas in a very simple form.

At this point, we have at our disposal all the ingredients needed for making explicit our approach. Using (41) in (28), we get

Deriving (42) with respect to the spatial coordinate $x$, we have

where $V_{N}(t,x):=-\partial_{x}(R_{N}(t,x)-\frac{1}{2}Q_{N}^{(p-1)}(t,x))$, vanishing in the $N\to\infty$ limit. We can then write the following equation

On the l.h.s. we recognize a Burgers hierarchy structure, while on the r.h.s. we have a source term (which further vanishes in the thermodynamic limit). The equilibrium dynamics of the Derrida models is then realized taking the limit $N\to\infty$, as summarized in the following

As a comparison, in Fig. 1, we reported the solutions of the self-consistency equations (52) for $p=2,\dots,8$ (solid curves), and the critical temperatures as predicted by the system (53) (dashed lines).

Going beyond the pure spin-glass case, in this Section we will approach the DAMs that are the main focus of this work.