A Theory of Solving TAP Equations for Ising Models with General Invariant Random Matrices

Introduction

TAP equations provide generalized mean field equations for statistical physics models with random, infinite range interactions which (under certain conditions) are assumed to be exact in the limit of an infinite system . In recent years there has been an increasing interest in such equations within and also outside of the statistical physics community. This is partly due to the fact that the TAP approach can be applied to statistical inference in probabilistic models in information theory , , statistics and machine learning , . Originally developed by D.J. Thouless, P. Anderson and R. Palmer for the Sherrington Kirkpatrick (SK) model of an Ising spin–glass , the TAP approach has been generalised to a variety of other problems . This includes models with continuous variables rather than Ising spins but also cases where the independent random interactions are replaced by other, more structured statistical ensembles that allow for certain dependencies.

While methods for deriving the TAP approach for different models are now well established it is not necessarily clear how the resulting system of nonlinear equations can be solved efficiently. A naive algorithm based on a simple iteration of the equations usually fails to achieve convergence. This problem has been addressed by a paper of Bolthausen for the case of the SK model . He has analyzed the dynamics of iterations rigorously and shown how the iterations can be altered in order to achieve exponential convergence (above the so–called AT line of stability). Other ideas to arrive at a convergent method are based on taking the limit of dense couplings in belief propagation algorithms, see , , and , (for rigorous analyses and ). Unfortunately, for this approach it is necessary to augment the original variables by auxiliary ones, such that the interactions in the new model are independent. For example, the Hopfield model can be represented by a bipartite graph of Ising spins and continuous variables. It is not clear how such a method should be set up for a matrix of interactions with more general statistical dependencies. Also, taking the dense coupling limit of an approximate message passing (AMP) algorithm valid for sparse coupling problems will not always lead to the correct dense coupling algorithm. Here we will construct a theory for dynamics using dense couplings as the starting point.

In this paper we will address the problem of solving TAP equations for Ising models with dense random coupling matrices with a general invariant probability distribution. Our analysis is based on dynamical mean field theory which allows to study the dynamics of iterative algorithms in the thermodynamic limit by a suitable average over the ensemble of couplings. It turns out that by including certain memory terms in the iteration, the effective field in the dynamics becomes a simple Gaussian random variable suggesting that the dynamics might converge. The explicit form of the memory terms depends explicitly on the statistical ensemble of couplings. We show that our method reproduces previous convergent algorithms for SK and Hopfield models. We also work out the details of our theory for the spin model with orthogonal random couplings . Simulations of the resulting algorithms show exponential convergence above a line of stability which can be identified with the so–called AT line.

The paper is organized as follows: in Section II we introduce the general random matrix formulation for the TAP equations. In Section III we present the results of dynamical functional theory. In Section IV, we introduce “the single-step memory construction” iterative algorithm for solving TAP equations. Section V is devoted to the derivation of the AT stability condition. Discussions and outlooks are presented in Section VI. Lengthy technical derivations are deferred to the Appendix.

General Invariant Random Matrix Ensembles

We will consider Ising models with pairwise interactions given by the Gibbs distribution for the spins {\mathchoice{\mbox{\boldmath\displaystyle S}}{\mbox{\boldmath\textstyle S}}{\mbox{\boldmath\scriptstyle S}}{\mbox{\boldmath\scriptscriptstyle S}}}=(S_{1},\ldots,S_{N})

We will later need the generating function of asymptotically invariant random matrices given by

with symmetric matrix $\textstyle Q$ having a finite rank. Since we believe that our paper might be of interest to researchers with an information theory background, we will briefly mention how ${\rm G}$ is related to quantities which are well known in the theory of free probability , which is a powerful approach to random matrix theory. Setting

one can show that $\rm R$ equals the so–called R-transform (in the theory of free probability) of the limiting spectrum. Its formal definition can be given in terms of the Cauchy transform: let ${\rm P}$ denote the limiting spectrum of $\textstyle J$. Moreover let

where $m_{n}$ is $n$th order moments of the limiting spectrum $\rm P$, i.e. $m_{n}=\int{\rm dP}(\lambda)\;\lambda^{n}$. Then the R-transform of $\rm P$ is given by

with ${\rm M}^{-1}$ denoting the composition inverse of $\rm M$. It admits the power series expansion

where $c_{n}$ are known as the free cumulants of $\rm P$. For example, the first two free cumulants $c_{1}$ and $c_{2}$ are the mean and variance of the distribution $\rm P$, respectively, i.e. $c_{1}=m_{1}$ and $c_{2}=m_{2}-m_{1}^{2}$. For details we refer the reader to .

In the sequel we give the R-transforms for some random matrix ensembles that we will discuss later. We first point out the simple identity

2 TAP Equations for General Invariant Couplings

TAP equations are a set of self-consistent equations for the vector of magnetisations {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}=\langle{\mathchoice{\mbox{\boldmath\displaystyle S}}{\mbox{\boldmath\textstyle S}}{\mbox{\boldmath\scriptstyle S}}{\mbox{\boldmath\scriptscriptstyle S}}}\rangle where the brackets denote expectation w.r.t. the Gibbs distribution (1). For a general invariant ensemble they have been derived first in using the large $N$ scaling of a perturbation expansion. A second derivation using the cavity method and the large $N$ limit of the ‘adaptive’ TAP equations can be found . We have provided a more rigorous derivation of the transition from “adaptive TAP” to the self-averaging limit using random matrix theory in Appendix B. The resulting TAP equations read

where q\triangleq\frac{1}{N}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}} and $\textstyle h$ is the vector of nonrandom external fields. Note, that the only dependency on the random matrix ensemble is via the R–transform ${\rm R}(1-q)$ in the so–called Onsager term which is a correction to the naive mean field term $\textstyle J$$\textstyle m$. One can show that $\Psi_{i}$ is the mean of the cavity field. Furthermore, following the calculations of one finds that $\Psi_{i}$ is Gaussian distributed (in the large $N$ limit) with respect to the random couplings $\bm{J}$ with mean $h_{i}$ and variance

Hence, the subtraction of the Onsager term {\rm R}(1-q){\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}} from the mean field $\textstyle J$$\textstyle m$ makes the remainder Gaussian. We will next transfer the idea of a Gaussian field from the static solutions to the dynamics of an algorithm.

The Results of Dynamical Functional Theory

Dynamical properties of disordered systems can be computed by the method of dynamical functionals . In the limit $N\to\infty$ this method provides us with exact results for the marginal distribution of a trajectory of a single variable (in our case a magnetization $m_{i}(t)$), when we define a dynamics of an algorithm for the solution of the TAP equations. As a typical result of such a calculation one finds that the ‘field’ {\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t) becomes a sum of a Gaussian term and a memory term which includes the magnetizations at all previous times. This memory often makes the dynamics of disordered systems highly complex allowing e.g. for a persistent dependency on the initial conditions and thus a failure to converge to a unique fixed point. Hence, we propose to introduce explicit memory terms which are chosen to cancel the implicit memory terms derived from the dynamical functional theory. In such a way, at each time step, the update of the magnetization for the algorithm involve a Gaussian distributed random field only and we expect that we might obtain good convergence results. This Gaussian property of the effective dynamical field was already shown for a Hopfield model in and (and proved in ) and reappeared in Bolthausen’s iterative construction of solutions to the TAP equations for the SK model in .

We start with defining a set of dynamical equations which could serve as a candidate algorithm for solving the TAP equations for invariant random coupling matrices

for $\tau=0,\ldots,t$ which depend on the field {\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t) and the previous local magnetizations {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(\tau). Here $f_{t}$ is an appropriate sequence of non-linear scalar functions. Our goal is to get the statistics of a single trajectory of (14)-(15), when $\textstyle J$ is a random matrix with generating function (3). To do so we make use of the dynamical functional theory (DFT) analysis as described in which is a discrete time version of the method of . We also refer the reader to where DFT was used to analyze the AMP algorithm in the context of the CDMA communication algorithm.

We introduce the generating functional corresponding to the dynamics (14)–(15) as

Notice that Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)={\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}}}\})=1. The statistics of the variables can be computed from the averaged generating functional \left<Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)\})\right>_{{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}}. In the large $N$ limit we obtain that (see C)

with $\mathcal{N}(\cdot;\mu,\Sigma)$ denoting the multivariate normal distribution with mean $\mu$ and covariance $\Sigma$. This result shows that in the large $N$ limit single trajectories can be treated as independent following the effective stochastic dynamical process given by

Here {\mathchoice{\mbox{\boldmath\displaystyle\phi}}{\mbox{\boldmath\textstyle\phi}}{\mbox{\boldmath\scriptstyle\phi}}{\mbox{\boldmath\scriptscriptstyle\phi}}}(t) is a vector of independent Gaussian random variables with covariance matrix $\mathcal{C}_{\phi}$ given by

where $\mathcal{G}$ and $\mathcal{C}$ are $T\times T$ the response and the correlation matrices, respectively. With slight abuse of notation, their $(t+1,\tau+1)$ indexed entries are given by

Moreover the specific random matrix ensemble enters the result through the coefficients $c_{n}$, see (7), and the memory matrix $\mathcal{\hat{G}}$ given by

So far we have not yet referred to the TAP equations in the DFT analysis. Instead we have considered a somewhat general dynamical system with disorder and memory. Such a formulation gives us enough freedom to construct a convenient dynamics which asymptotically converges to the solution of the TAP equations. We will define the dynamics to be of the form

where the variables $\psi_{i}(t)$ must be chosen to become independent Gaussian fields in the resulting effective single variable dynamics (19). In fact, there are actually various methods for doing so. In the sequel we will limit our attention to a method that we call the single step memory construction.

The Single Step Memory Construction

In the single step memory algorithm we will construct the update in such a way that the resulting memory term (23) satisfies the equation

Hence, if (25) holds, then using (19) we find that the variable

becomes a a Gaussian field. We will choose the field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) in (24) as a linear combination of the Gaussian fields {\mathchoice{\mbox{\boldmath\displaystyle\phi}}{\mbox{\boldmath\textstyle\phi}}{\mbox{\boldmath\scriptstyle\phi}}{\mbox{\boldmath\scriptscriptstyle\phi}}}(\tau), $\tau=1,\ldots,t$ of the form

where we have to construct the non-random terms ${\mathcal{A}}(t+1,\tau)$ to make the dynamical order parameters consistent with the single step memory condition (25). This condition leads to a very simple result for the response function (21) because there is no complicated propagation in time of a response to an external field. In fact, from (24) we obtain for the response function (21)

with q(t)\triangleq\frac{1}{N}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t). Thus we have the explicit result

Finally, using (23) we get an explicit result for the response function in terms of the memory terms $\mathcal{\hat{G}}(t,t-1)$. Note that by construction of the single step memory matrix $\mathcal{\hat{G}}$ (25) we can write (23) as

where the coefficients $a_{n}$ are obtained from the power series expansion of the composition inverse of the R-transform:

By definition the trace of $\textstyle J$ is zero, i.e. $\rm R(0)=0$. Hence, the power series expansion in (32) starts from the first order term.

To complete the specification of the single step memory construction we only need to specify $\mathcal{\hat{G}}(t,t-1)$. This will be chosen such that the method is asymptotically consistent with the static TAP equations. Specifically, from (12) we should have

which assuming convergence $q(t)\to q$ as $t\to\infty$ leads to (33). This form has also the advantage, that in (30), for the update an unwanted factor $1-q(t+1)$, which would make {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) depending on the future state {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t+1), cancels.

Putting everything together the single-step memory algorithm for $t\geq 0$ is defined as

such that $Q(-1)=1$. The memory term $\hat{\mathcal{G}}(t,t-1)$ is given by (34). Moreover the algorithm initializes with {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)={\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}}} for $t\in\{-1,0\}$.

2 Asymptotic Consistency with TAP Equations

In the sequel we show that if the single step memory algorithm (35)-(37) converges, it solves the TAP equations (11)-(12). Let us assume that {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)\to{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}} as $t$ tends to infinity. To have the convergence to the TAP equations we solely need to show that the sum in (36) converges to the proper limit. From (27) and (30) we must have

We make the so-called weak long-term response assumption that

Hence, for sufficiently large $t$ and $\tau^{\prime}<t$ such that $t/\tau^{\prime}$ being finite as $t\to\infty$, we can write

Next we will provide the details of the single-step memory algorithm for the SK, Hopfield and random orthogonal models.

3 Example 1 The SK-Model

Recall that, for the standard SK model we have ${\rm R}(x)=\beta^{2}x$; so that ${\rm R}^{-1}(x)=x/\beta^{2}$. Hence, $a_{1}=1/\beta^{2}$ and $a_{n}=0$ for $n>1$. Thus, the single-step memory algorithm may written as

At first glance, these dynamical equations are similar but not exactly equal to those proposed by Bolthausen . The difference is that instead of the dynamical order parameter $q(t)$ the fixed point solution of $q$ appears. Using the explicit form of the covariance of the field $\psi_{i}(t)$ given by (60) in the next section, one finds for the field variance $\langle\left(\psi_{i}(t)-\langle\psi_{i}(t)\rangle\right)^{2}\rangle=\beta^{2}q(t)$. Hence, if we start the iteration (as in ) with $m_{i}(1)=\sqrt{q}$ such that $q(1)=q$, then we find that in the large $N$ limit, we also have $q(t)=\left\langle\tanh^{2}(\psi_{i}(t-1)\right\rangle=q$ for all times $t$ and we get agreement with .

4 Example 2 The Hopfield Model

Thus the memory coefficients are given as $a_{n}=1/(\alpha\beta^{n+1})$ for $n\geq 1$. In the sequel we show that the single step memory algorithm for the Hopfield model coincides with AMP algorithm which was introduced in the context of the CDMA problem in and compressed sensing in . From (36) we first write

Notice that from (37) we may write (49) in the form of

Then, defining {\mathchoice{\mbox{\boldmath\displaystyle z}}{\mbox{\boldmath\textstyle z}}{\mbox{\boldmath\scriptstyle z}}{\mbox{\boldmath\scriptscriptstyle z}}}(t)\triangleq{\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t)-A(t){\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t), we write the single step memory algorithm as

where $\textstyle\bf I$ is the identity matrix of appropriate dimension. Note, that ({\mathchoice{\mbox{\boldmath\displaystyle\bf I}}{\mbox{\boldmath\textstyle\bf I}}{\mbox{\boldmath\scriptstyle\bf I}}{\mbox{\boldmath\scriptscriptstyle\bf I}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}/\beta) asymptotically coincides with the corresponding central Wishart matrix (see Section 2). Thereby we exactly obtain the AMP iteration steps as introduced in . We also refer the reader to the related works and , where the dynamics of the AMP algorithm is analyzed by means of DFT and Bolthausen’s conditioning technique , respectively.

Bolthausen’s conditioning technique for SK model and Hopfield model are based on the assumption that the entries of the underlying coupling matrix are defined via zero-mean iid and Gaussian distributed random variables, see Section 2.1. Recently, it has been shown in that the same analyses can be obtained without the need of Gaussian distribution assumption but a sub-Gaussian tail condition of the distribution is required. Indeed, thanks to the central limit theorem, one can show that the generating function (3) also yields the same result regardless of whether the Gaussian distribution assumption is considered or not, see [23, Section 5].

5 Example 3 The Random Orthogonal Model

For the random orthogonal model from (10) we have ${\rm R}^{-1}(x)=x/(\beta^{2}-x^{2})$. This yields the memory coefficients as

we illustrate the convergence of the single-step memory algorithm for the random orthogonal model obtained by running simulations. Notice that after few iteration steps the convergence becomes exponentially fast. The flat lines around $(-300)$dB, i.e. $10^{-30}$, are the consequence of the machine precision of the computer which was used.

The convergence improves with increasing temperature parameter $1/\beta$. On the other hand, for large enough $\beta$, the algorithm fails to converge. To estimate the critical parameter, we study the inverse decay time measured by the angle $\theta$ as illustrated in Figure 1 and extrapolate the simulational data to $\theta=0$ using a convenient range of $\beta$. One might expect that the critical $\beta$ would coincide with the one obtained from a de Almeida–Thouless (AT) stability condition which can also be derived from the TAP approach . The AT line is given by the equation

where the random variable $\psi_{i}$ is a Gaussian with mean $h_{i}$ and variance $q{\rm R}^{\prime}(1-q)$. Note, that contains a typo in the corresponding expression. In Figure 3,

we present a comparison between the simulations and (55). This coincidence can be understood from a dynamical point of view by analyzing the stability of the dynamics close to the fixed point. The details will be postponed to Section V. Finally, it also worth noting that the trajectories of the algorithm show a self-averaging behaviour above the AT line for large $N$. On the other hand, we find that below the AT line, there are strong sample to sample fluctuations. However, by averaging order parameters over many samples, we get a good agreement with the theory. Since the main goal of this paper is to present a convergent algorithm we will leave a more careful investigation of this point to future publications.

6 The Field Covariance Matrix

In order to compare simulations of systems with analytical results obtained from the dynamical functional approach in the limit $N\to\infty$ and to study the stability of TAP fixed–points we have to perform expectations over the Gaussian random variables $\psi_{i}(t)$ (see (27)). Specifically we write

where the order parameter $\mathcal{C}$ is given by

Here ${\rm P}$ denotes a two-dimensional Gaussian distribution with zero mean. The mean of the field $\psi_{i}(t)$ follows from (27) as

Hence, we need to compute the corresponding covariance matrix which is defined as

In this expression, for a power series $f(x,y)=\sum_{n,k\geq 0}a_{n}b_{k}x^{n}y^{k}$, we have introduced the symbol ${\rm Co}_{x^{n}y^{k}}[f(x,y)]\triangleq a_{n}b_{k}$ for its coefficients. Moreover the function $\rm A$ is defined as

The function ${\rm A}$ has a relatively simple form for the three random matrix ensembles considered in this paper. For the SK and Hopfield models we have ${\rm A}(x,y)=xy/\beta^{2}$ and ${\rm A}(x,y)=xy/(\beta^{2}\alpha)$, respectively. Moreover, for the random orthogonal model, from ${\rm R}^{-1}(x)=x/(\beta^{2}-x^{2})$ we have ${\rm Co}_{x^{n}y^{k}}[{\rm A}(x,y)]=\delta_{nk}(-1)^{n+1}/\beta^{2n}$. We next compare our simulations with theoretical results. We used the initializations {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)={\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}}} for $t\in\{0,1\}$, hence we assign $\mathcal{C}(1,0)=0$. In Figure 4 and 5

we show such a comparison above the AT line. Note, that no averaging over coupling matrices was used for the simulations. The integration over two-dimensional correlated Gaussian distribution used to calculate $\mathcal{C}(t,t-1)$ was performed numerically. However the accuracy of the numerical method limits us for providing very precise results as $t$ grows. In Figure 4 we illustrate the theoretical prediction of \frac{1}{N}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(m-1) by the order parameter $\mathcal{C}(t,t-1)$ for a large range of $t$. Below the AT line,

the single-step memory algorithm diverges and simulated trajectories show strong sample fluctuations. However, by taking an average over a large number of trajectories we obtain a good agreement with the theory (see Figure 6).

7 Asymptotic consistency with Cavity Variance

In Section 4.2 we have demonstrated the convergence of the single step memory algorithm to the TAP equations. In a similar way one can show that (60) converges to the variance of the static field variance in (13) as $t$ and $t^{\prime}$ tend to infinity. Specifically, by invoking the weak long-term response assumption (40) in (20) and following similar steps as in (D), for sufficiently large $t$, $\tau<t$, $t^{\prime}$ and $\tau^{\prime}<t^{\prime}$ such that $\tau/t$ and $\tau^{\prime}/t^{\prime}$ being finite as $t$ and $t^{\prime}$ tend infinity, one can show that

In this expression, by abuse of notation, we denote $\lim_{y\to x}{\rm A}(x,y)$ by ${\rm A}(x,x)$. Its explicit form is given by

The stability of TAP fixed points

In order to analyze the stability of the fixed points of the single step algorithm we resort to a linear stability analysis. We add Gaussian white noise to the dynamics, i.e. we set $\psi_{i}(t)\rightarrow\psi_{i}(t)+\epsilon_{i}(t)$ with $\langle\epsilon_{i}(t)^{2}\rangle=\epsilon$ and discuss the limit $\epsilon\to 0$. If the static TAP fixed point is stable, then the system should asymptotically show only small stationary fluctuations around this and we can work in the Fourier domain. Hence, we assume

Inserting these Fourier representations into (60), for large $t$ and $t^{\prime}$ we may write (see (62)-(64))

For small noise $\epsilon\to 0$, the assumption of stability translates into small fluctuations around the static solution and we can write

where we have separated fluctuations into static and dynamical parts. We will analyse the dynamical part next, but note, that also the static part $q$ will have contributions from $\epsilon$. Thus for $\omega\neq 0$ we have

where now the value of $q$ is computed for $\epsilon=0$. We next express $\hat{c}(\omega)$ in terms of $\hat{c}_{\psi}(\omega)$ for small $\epsilon$. The calculation in E is based on expanding

up to first order in $\epsilon$. The brackets denote expectations over the two dimensional Gaussian field $(u(t),u(t^{\prime}))$ with $\langle u(t)u(t^{\prime})\rangle\simeq s_{0}+\epsilon(s(t-t^{\prime})+\delta_{t,t^{\prime}})$ for $\epsilon\to 0$, where $s_{0}=q{\rm R}^{\prime}(1-q)$ and $s(t-t^{\prime})=c_{\psi}(t,t^{\prime})$. For $t=t^{\prime}$ the integral is over a single Gaussian only. The calculation shows that

with $\alpha$ is defined as in (55). Combining this relationship with (74) we have

In fact, for the SK, Hopfield and random orthogonal models, we have ${\rm Co}_{x^{n}y^{k}}[{\rm A}(x,y)]=0$ $\forall n\neq k$. Therefore (77) is actually independent of $\omega$ and from (65) we explicitly have that

In general, the right hand side of (77) must be non–negative to have a valid representation as a Fourier-transform of a time dependent correlation function. While the term ${\rm A}(e^{-i\omega}{\rm R}(1-q),e^{i\omega}{\rm R}(1-q))$ is always positive, see (71), the second term is small and positive for sufficiently small $\beta$. But it changes sign and diverges. One expects that the divergence will occur first for the long range fluctuations, i.e. for the limit of low frequencies. Taking the limit yields

for the onset of instability agrees with the well–known AT stability criterion .

Discussion and Outlook

Acknowledgment

The authors would like to thank Florent Krzakala for inspiring us to do this study. This work was partially supported by the European Commission in the framework of the FP7 Network of Excellence in Wireless COMmunications NEWCOM$\sharp$ (Grant agreement no. 318306).

where for the random variables $X$ and $Y$, ${\rm Var}[X]$ and ${\rm Cov}[X,Y]$ denoting the variance $X$ and the covariance of $X$ and $Y$, respectively. For the proof, we basically need to show that

To do this we make use of the so-called (orthogonal) Weingarten calculus that allows to for calculate joint moments of Haar entries (analogeous to Wick calculus for Gaussian matrices). For details we refer the reader to . From [27, Theorem 2.1, Example 2.1] we have

Furthermore, we have $\left<O_{ij}^{2}\right>=1/N$, $\forall i,j$. Thus, the variance (83) reads

Appendix B The self-averaging limit of adaptive TAP Equations

We will provide a derivation of the TAP equations (11)-(12) from the ‘adaptive TAP’ approach of . Under the assumption of Gaussian distributed cavity fields and an approximate linear response argument one finds

with the positive definite matrix {\mathchoice{\mbox{\boldmath\displaystyle\chi}}{\mbox{\boldmath\textstyle\chi}}{\mbox{\boldmath\scriptstyle\chi}}{\mbox{\boldmath\scriptscriptstyle\chi}}}\triangleq({\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}})^{-1}.

To obtain the TAP equations in (11)-(12) we basically need to show that $V_{i}\simeq{\rm R}(1-q)$. To that end we write (88) in the form of

where we denote the R-transform of the spectrum of an $N\times N$ symmetric matrix $\textstyle X$ by {\rm R}_{{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}}^{N}. Note that we have 1-q=\frac{1}{N}{\rm tr}({\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}})^{-1}. Then, from Lemma 1 below we have

where $V$ is defined through the implicit equation $(1-q)=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\Lambda_{ii}-V}$. This implies that {\rm tr}({\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}})^{-1}\simeq{\rm tr}({\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-V{\mathchoice{\mbox{\boldmath\displaystyle\bf I}}{\mbox{\boldmath\textstyle\bf I}}{\mbox{\boldmath\scriptstyle\bf I}}{\mbox{\boldmath\scriptscriptstyle\bf I}}})^{-1}. Here, from (94) to (95) we make use of the identity (106) below. Note that for a non-negative $N\times N$ matrix $\textstyle X$ we can write the Stieltjes transform of the eigenvalue distribution of $\textstyle X$, say {\rm P}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}}, as {\rm M}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}}(\omega)=\int{\rm dP}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}}(x)/(\omega-x) with $\omega\in(-\infty,0)$. Then, from (91) we use the subordination property [28, Chapter 22] as {\rm M}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}}(\omega)\simeq{\rm M}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}}(\omega+{\rm R}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}}(-{\rm M}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle\Lambda}}{\mbox{\boldmath\textstyle\Lambda}}{\mbox{\boldmath\scriptstyle\Lambda}}{\mbox{\boldmath\scriptscriptstyle\Lambda}}}-{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}}(\omega))) and take the limit $\omega\to 0$. Doing so yields $V\simeq{\rm R}(1-q)$ where we note that {\rm R}(x)=\lim_{N\to\infty}{\rm R}^{N}_{{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}}(x). Notice also that

which yields $V_{i}\simeq V={\rm R}(1-q)$.

Let an $N\times N$ matrix $\textstyle X$ be positive definite. Let Q=\frac{1}{N}{\rm tr}({\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}^{-1}). Then,

For convenience let \eta(\epsilon)\triangleq\frac{1}{N}{\rm tr}(({\mathchoice{\mbox{\boldmath\displaystyle\bf I}}{\mbox{\boldmath\textstyle\bf I}}{\mbox{\boldmath\scriptstyle\bf I}}{\mbox{\boldmath\scriptscriptstyle\bf I}}}+\epsilon{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}})^{-1}) for $\epsilon>0$. Since $\textstyle X$ is positive definite we can write

Applying the substitution $\omega\triangleq t\eta(t)$ to the following integral we have

with $Q(\epsilon)\triangleq\epsilon\eta(\epsilon)$. In other words we have

Taking the limit $\epsilon\to\infty$ we complete the proof.

Appendix C Derivation of DFT Results

For the sake of notational compactness let

By the Fourier representation of the Dirac function we write

The derivation is separated into two parts: i) disorder average; ii) the saddle point method.

For convenience let us introduce $N\times T$ matrices $\textstyle X$ and $\textstyle\hat{X}$ with $X_{nt}=\frac{m_{n}(t+1)}{\sqrt{N}}$ and $\hat{X}_{nt}=\frac{\hat{\gamma}_{n}(t+1)}{i\sqrt{N}}$. We need to evaluate

with {\mathchoice{\mbox{\boldmath\displaystyle Q}}{\mbox{\boldmath\textstyle Q}}{\mbox{\boldmath\scriptstyle Q}}{\mbox{\boldmath\scriptscriptstyle Q}}}={\mathchoice{\mbox{\boldmath\displaystyle\hat{X}}}{\mbox{\boldmath\textstyle\hat{X}}}{\mbox{\boldmath\scriptstyle\hat{X}}}{\mbox{\boldmath\scriptscriptstyle\hat{X}}}}{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}^{\dagger}+{\mathchoice{\mbox{\boldmath\displaystyle X}}{\mbox{\boldmath\textstyle X}}{\mbox{\boldmath\scriptstyle X}}{\mbox{\boldmath\scriptscriptstyle X}}}{\mathchoice{\mbox{\boldmath\displaystyle\hat{X}}}{\mbox{\boldmath\textstyle\hat{X}}}{\mbox{\boldmath\scriptstyle\hat{X}}}{\mbox{\boldmath\scriptscriptstyle\hat{X}}}}^{\dagger}. Here (109) follows directly from (3)-(7). We will evaluate

Then by using cyclic invariance of the trace we obtain the expression

C.2 The saddle point calculation

By the Fourier representation of Dirac function we write the last line of (116) as

Here $c$ is a constant irrelevant for the saddle point calculation. We define the auxiliary single-site partition function

Thus, for large $N$ we get the factorization of the generating function

Finally, notice that $-i\langle m_{n}(t)\hat{\gamma}_{n}(s)\rangle_{\phi_{n}}=\langle\frac{\partial m_{n}(t)}{\partial\phi_{n}(s)}\rangle_{\phi_{n}}$. Thus, $\mathcal{G}$ equals the response–function (21). This completes the derivation.

Appendix D Derivation of equation (60)

First note that from (27) and (30) we have

where $\mathcal{C}_{\phi}$ is defined as in (20). For sake of compactness we let $f(x)={\rm R}^{-1}(x)$. By elementary combinatorics and using $\mathcal{G}={f}(\mathcal{\hat{G}})$ we can show that any power of the matrix $\mathcal{G}$ can be written as

where for any power series $f(x)=\sum_{n}a_{n}x^{n}$ we represent the coefficient via the definition $a_{k}\triangleq{\rm Co}_{x^{k}}(f(x))$. This means that we have

where we have extended the definitions of coefficients ${\rm Co}_{x^{n}y^{k}}[f(x,y)]\triangleq a_{n}b_{k}$ to double power series $f(x,y)=\sum_{n,k\geq 0}a_{n}b_{k}x^{n}y^{k}$. Summing up the geometric series, we have

Putting everything together completes the proof.

Appendix E Derivation of equation (76)

For the sake of compactness, without loss of generality, we may set $h_{i}=h$. Using the representation of the Gaussian density in terms of the characteristic function we have the expansion for $t\neq t^{\prime}$

The last line is obtained by representing $k_{1}k_{2}$ etc by derivatives with respect to $u_{1}$ and $u_{2}$. Repeating a similar equation for $t=t^{\prime}$, we get

Both expansions can be represented in the single equation

Note, that only the second term contributes to the dynamic part of the fluctuations. Hence, by taking the Fourier transform and noting that $\frac{\partial\tanh(u+h)}{\partial u}=1-\tanh^{2}(u+h)$ the result is obtained.

References