Inverse problems for structured datasets using parallel TAP equations and RBM

TL;DR

This paper introduces a parallel TAP-based algorithm for solving inverse problems in binary clustered datasets, demonstrating its effectiveness in the Hopfield model and exploring its relation to RBMs, with advantages over traditional methods.

Contribution

It presents a novel parallel TAP equations approach for inverse problems, capable of retrieving exact patterns and handling large systems, and compares it with RBM methods.

Findings

Parallel TAP equations accurately retrieve teacher patterns.

AMP equations do not match expected behavior in the direct problem.

The approach scales well to large system sizes.

Abstract

We propose an efficient algorithm to solve inverse problems in the presence of binary clustered datasets. We consider the paradigmatic Hopfield model in a teacher student scenario, where this situation is found in the retrieval phase. This problem has been widely analyzed through various methods such as mean-field approaches or the pseudo-likelihood optimization. Our approach is based on the estimation of the posterior using the Thouless-Anderson-Palmer (TAP) equations in a parallel updating scheme. At the difference with other methods, it allows to retrieve the exact patterns of the teacher and the parallel update makes it possible to apply it for large system sizes. We also observe that the Approximate Message Passing (AMP) equations do not reproduce the expected behavior in the direct problem, questioning the standard practice used to obtain time indexes coming from Belief Propagation (BP). We tackle the same problem using a Restricted Boltzmann Machine (RBM) and discuss the analogies between the two algorithms.

Equations

$$\mathcal{L}(\{{\underline{\xi }}^{\mu }{\}}_{\mu =1,\dots ,P}|\mathcal{D})=\frac{\beta }{2N}\sum _{ij}\sum _{\mu =1}^P\sum _{a=1}^M{s}_i^a{s}_j^a{\xi }_i^{\mu }{\xi }_j^{\mu }-M\log Z(\{{\xi }_i^{\mu }{\}}_{\mu =1,\dots ,P}).$$

$$m_i^{t+1}=\tanh \left(\beta \sum _{j=1}^N{J}_{ij}{m}_j^t-\frac{\alpha \beta }{1-\beta (1-q^t)}{m}_i^t\right),$$

$$\hat{\mathcal{L}}=\frac{1}{2N}\sum _{ij}\sum _{\mu =1}^{P^{\prime }}\sum _{a=1}^M{s}_i^a{s}_j^a{m}_i^{\mu }{m}_j^{\mu },$$

$$P(\{W_i^{\mu }{\}}_{\mu =1,\dots ,P}|\mathcal{D})\propto$$

$$Z^{-1}(\underline{W})\prod _{a=1}^M\int d\underline{\lambda }{e}^{-{\sum }_{\mu }{\lambda }_{\mu }^2/2+{\sum }_{i,\mu }{W}_i^{\mu }{s}_i^a{\lambda }_{\mu }},$$

$${\mathcal{S}}_r=\frac{1}{M}\sum _{a=1}^M\log ({⟨p(s_r^a|\underline{\lambda }⟩}_{P(\underline{\lambda }|{\underline{s}}^a,\{W\})}).$$

$$m_i^{t+1}=\tanh \left(\beta \sum _{j=1}^N{J}_{ij}{m}_j^t-\frac{\alpha \beta }{1-\beta (1-q^t)}{m}_i^t\right)$$

$$H_i^{t+1}=\frac{1}{1-u^t}\left[\sum _{j\ne i}{J}_{ij}{m}_j^t-u^t{H}_i^t-\frac{\alpha u^t}{1-u^{t-1}}{m}_i^{t-1}\right],$$

$$m_i^{t+1}\simeq \beta \sum _{j=1}^N{J}_{ij}{m}_j^t+O(\alpha )$$

$${\tilde{m}}_{\lambda }^{t+1}\simeq \beta \lambda {\tilde{m}}_{\lambda }^t,$$

$${\tilde{m}}_{\lambda }^{t+1}=\frac{\beta }{1-\beta }\left[(\lambda -1){\tilde{m}}_{\lambda }^t-O(\alpha )\right].$$

$$m_i^{t+1}$$

$$z_i^t$$

$$Z(\{{\xi }_i^1,{\xi }_i^2\})=\sum _s{e}^{\frac{\beta }{N}{\sum }_{ij}({\xi }_i^1{\xi }_j^1+{\xi }_i^2{\xi }_j^2)s_i{s}_j},$$

$$H_{\xi }(s)=\frac{2\beta }{N}\sum _{i\in S,j\in S}{\xi }_i^1{\xi }_j^1{s}_i{s}_j+\frac{2\beta }{N}\sum _{i\in \overline{S},j\in \overline{S}}{\xi }_i^1{\xi }_j^1{s}_i{s}_j,$$

$$Z(\{{\xi }_i^1,{\xi }_i^2\})=Z_{cw,\beta (1+q)}^{N(1+q)/2}{Z}_{cw,\beta (1-q)}^{N(1-q)/2},$$

$$\varphi (q)=\frac{1+q}{2}{f}_{cw}(\beta (1+q))+\frac{1-q}{2}{f}_{cw}(\beta (1-q))$$

$$P(\underline{s},\underline{\lambda }|\{W\})=Z^{-1}(\{W\})e^{-E_W(\underline{s},\underline{\lambda })}$$

$$Z(\{W\})=\int \prod _{\mu =1}^{N_h}dP({\lambda }_{\mu })\sum _{\underline{v}}{e}^{-E_W(\underline{s},\underline{\lambda })}.$$

$$P(\underline{s}|\{W\})=\int \prod _{\mu =1}^{N_h}dP({\lambda }_{\mu })P(\underline{s},\underline{\lambda }|\{W\})$$

$$\mathcal{L}=\frac{1}{M}\sum _{a=1}^M\left(-\log Z(\underline{W})+\log \int dP(\underline{\lambda })e^{-E_W({\underline{s}}^a,\underline{\lambda })}\right).$$

$$W_i^{\mu }=W_i^{\mu }+\left({⟨{\lambda }_{\mu }{s}_i⟩}_D-{⟨{\lambda }_{\mu }{s}_i⟩}_{RBM}\right)$$

$${⟨{\lambda }_{\mu }{s}_i⟩}_D=M^{-1}\sum _a^M\int dP(\underline{\lambda })P(\underline{\lambda }|{\underline{s}}^a,\{W\}){\lambda }_{\mu }{s}_i^a$$

$${⟨{\lambda }_{\mu }{s}_i⟩}_{RBM}=\frac{\partial }{\partial W_i^{\mu }}\log Z(\underline{W}).$$

$${\mathcal{S}}_r=\frac{1}{M}\sum _{a=1}^M\log ({⟨p(s_r^a|\underline{\lambda }⟩}_{P(\underline{\lambda }|{\underline{s}}^a,\{W\})}),$$

$${⟨p(s_r^a|\underline{\lambda }⟩}_{P(\underline{\lambda }|{\underline{s}}^a,\{W\})}=\int dP(\underline{\lambda })p(s_r^a|\underline{\lambda })P(\underline{\lambda }|{\underline{s}}^a,\{W\})$$

$${⟨p(s_r^a|\underline{\lambda }⟩}_{P(\underline{\lambda }|{\underline{s}}^a,\{W\})}={\mathcal{N}}_a^{-1}\int dP(\underline{\lambda })e^{{\sum }_{\mu k}{W}_k^{\mu }{s}_k^a{\lambda }_{\mu }},$$

$${\mathcal{N}}_a=\sum _{s_r^a}\int dP(\underline{\lambda })e^{{\sum }_{k\mu }{W}_k^{\mu }{s}_k^a{\lambda }_{\mu }}.$$

$$\frac{\partial {\mathcal{S}}_r}{\partial W_r^{\nu }}=\frac{1}{M}\sum _{a=1}^M\left(\frac{\int {\prod }_{\mu }dP({\lambda }_{\mu })s_r^a{\lambda }_{\nu }{e}^{{\sum }_{k\mu }{W}_k^{\mu }{s}_k^a{\lambda }_{\mu }}}{\int dP(\underline{\lambda })e^{{\sum }_{k\mu }{W}_k^{\mu }{s}_k^a{\lambda }_{\mu }}}-\frac{\partial }{\partial W_r^{\nu }}\log \sum _{s_r^a}\int dP(\underline{\lambda })e^{{\sum }_{k\mu }{W}_k^{\mu }{s}_k^a{\lambda }_{\mu }}\right).$$

$$P(\underline{\lambda }|{\underline{s}}^a,\{W\})=\frac{e^{{\sum }_{k,\mu }{W}_k^{\mu }{\lambda }_{\mu }{s}_k^a}}{\int dP(\underline{\lambda })e^{{\sum }_{k,\mu }{W}_k^{\mu }{\lambda }_{\mu }{s}_k^a}}$$

$${⟨{\lambda }_{\nu }{s}_r⟩}_D=\frac{1}{M}\sum _{a=1}^M\int dP(\underline{\lambda }){\lambda }_{\nu }{s}_r^{a}P(\underline{\lambda }|{\underline{s}}^a,\{W\}).$$

Full text

Inverse problems consist in inferring information about the structure of a system from the observation data of configurations. Cases where the system’s variables $s_{i}$ are binary can be studied in the framework of the inverse Ising model, whose parameters $\{J_{ij},h_{i}\}$ are tuned in order to describe the observed configurations according to the Boltzmann weight $P(s)\sim\exp(\sum_{i<j}J_{ij}s_{i}s_{j}+\sum_{i}h_{i}s_{i})$. This is the simplest distribution emerging when using the maximum entropy approach in order to reproduce exactly the one and two points statistics of the data. One of the most successfull application of this method is the 3D reconstruction of protein structures morcos2011direct . Inverse problems arise also in collective behavior bialek2014social , immunology wood2012mechanism , neural activity schneidman2006weak ; cocco2009neuronal and financial time series bury2013market . In general, inferring the parameters of the model is a challenging problem because maximizing the likelihood involves the computation of the partition function $Z=\sum_{s}P(s)$, which is impossible in most of the realistic cases. On the other hand, when dealing with time-series, it is possible to use a simples approach modeling the sampling process by a stochastic parallel dynamics analysed in roudi2011mean , optimized in decelle2015inference and generalized in dunn2013learning ; campajola2018inference . A recent review on this subject can be found in nguyen2017inverse .

The original attempt to solve the problem was based on an Expectation-Maximization algorithm known as Boltzmann learning ackley1985learning . This method is unpractical on large systems unless heuristic methods, like Monte Carlo (MC) sampling, are used to estimate correlations huang2010reconstructing . Nevertheless MC is slow and thus many sophisticated techniques coming from statistical mechanics and machine learning have been proposed as alternative approaches kappen1998efficient ; tanaka2000information ; sohl2011new ; cocco2011adaptive ; aurell2012inverse ; ricci2012bethe ; nguyen2012mean ; cocco2012adaptive ; raymond2013mean ; decelle2014pseudolikelihood ; lokhov2018optimal ; franz2019fast . These methods, however, share one or both of the following shortcomings: i) they require a large number of observations and ii) the overall performance drops significantly when the dataset is structured. This is often the case when data is produced from a (sub)set of many attractive states, or data is collected in different regimes, e.g. quiescent and spiking regimes in neural networks. This problem becomes particularly relevant at low temperatures and it has already been studied both in the sparse braunstein2011inference and in the dense case cocco2011high ; decelle2016solving . Pseudo-likelihood besag1977efficiency based methods decelle2016solving were shown to be the best options in a wide range of temperatures. Here, we present two algorithms to compete with the existing state-of-the-art by posing the problem in a Bayesian framework using the Thouless-Anderson-Palmer (TAP) equations thouless1977solution and the Restricted Boltzmann Machine (RBM) hinton2002training . Our TAP-based algorithm will be shown to achieve a better quality of the results by observing far fewer configurations in the clusterized phase. Moreover, it allows to consider larger system size with respect to those studied in cocco2011high ; decelle2016solving . As a side aspect of this work, we discuss the behavior of the Approximate Message Passing (AMP) equations for the Hopfield model mezard2017mean , and show that they do not do not reproduce the thermodynamic expected behavior.

We consider a dataset with many clusters by drawing configurations from the Hopfield model hopfield1982neural ; amit1985storing . Given a set of $N$-dimensional binary independent patterns $\{\underline{\zeta}^{\mu}\}_{\mu=1,\dots,P}$, teacher’s patterns, the coupling matrix of the associated Hopfield model is defined as $J_{ij}=N^{-1}\sum_{\mu}\zeta^{\mu}_{i}\zeta^{\mu}_{j}$ and its Hamiltonian is $H_{\zeta}(\underline{s})=-1/2\sum_{ij}J_{ij}s_{i}s_{j}$. Next, we construct a set of equilibrium configurations $\mathcal{D}=\{\underline{s}^{a}\}_{a=1,\dots,M}$ sampled from the Boltzmann distribution $P(\underline{s})=Z^{-1}\exp[-\beta H_{\zeta}(\underline{s})]$, $\beta$ being the inverse temperature. The set $\mathcal{D}$ is given to a student whose task is to infer the teacher’s patterns. This task cannot be accomplished with the methods developed in cocco2011high ; decelle2016solving , where the focus was the inference of the coupling matrix $J$. Using a uniform prior on teacher’s patterns, the posterior distribution is proportional to the likelihood.

For $P=1$, the Hopfield model is nothing but a Curie-Weiss model. In this case, the problem is called the dual Hopfield model barra2017phase ; barra2018phase , and the log-likelihood is $\mathcal{L}(\underline{\xi}|\mathcal{D})=\frac{\beta}{2N}\sum_{ij}\sum_{a=1}^{M}s^{a}_{i}s^{a}_{j}\xi_{i}\xi_{j}\>$, where $\underline{\xi}$ denote the student’s pattern. This is readily established by absorbing the $\xi$-dependence of the partition function into a redefined set of variables $\underline{s}$ via $s^{\prime}_{i}=\xi_{i}s_{i}$. On the other hand, for $P>1$ this transformation is not feasible. In this case, the log-likelihood reads

[TABLE]

TAP equations mezard1990spin ; nakanishi1997mean ; shamir2000thouless ; kabashima2001tap ; mezard2017mean describe the stationary points of the free energy and their use as an inference method has been pioneered in opper1996mean ; kappen1998boltzmann ; tanaka1998mean , as well as in kabashima1998belief ; kabashima2001tap , where their relationship between the Bayes theorem and the Belief Propagation (BP) equations was discussed. These works paved the way to their applications in a number of other problems saad2009line such as error correcting codes, compressed sensing and learning in neural networks, as discussed in the recent review zdeborova2016statistical .

For the direct problem, where the coupling matrix is $J_{ij}=N^{-1}\sum_{\mu}\zeta^{\mu}_{i}\zeta^{\mu}_{j}$, TAP equations estimate the local magnetizations. For the inverse problem, where we are interested in retrieving the model’s parameters given the dataset $\mathcal{D}$, we propose an algorithm based on these equations to estimate the posterior in Eq. 1. Even if TAP and mean field methods have already been used to solve this inverse problem braunstein2011inference ; cocco2011high ; decelle2016solving ; nguyen2012mean ; gabrie2015training , this approach is completely different from the previous ones: dealing with the dual model is the key to improve the quality of the reconstructed network. On the dual model, the role of spins and patterns is exchanged: $m_{i}=\left<\xi_{i}\right>$ and the $M$ sampled configurations play the role of the student’s patterns. We notice that a similar approach has been indipendently proposed in hou2019minimal for an RBM with $2$ hidden binary units, using BP.

We start by considering the simplest case $P=1$. We introduce a naive time indexing for an iterative scheme of the TAP equations,

[TABLE]

where $J_{ij}=N^{-1}\sum_{a}s^{a}_{i}s^{a}_{i}$, $\alpha=P/N$ and $Nq^{t}=\sum_{i}(m_{i}^{t})^{2}$. The entire set of magnetizations $\underline{m}^{t}$ are updated to achieve $\underline{m}^{t+1}$ in a parallel way. In principle, any sophisticated time indexing schemes can be employed as long as it achieves the convergence to a physical state. Particularly, the so-called Approximate Message Passing (AMP) equations has been the focus of many studies in inference problems zdeborova2016statistical . This scheme is inspired by the convergence issues of naive indexing in SK model even in the replica symmetric phase kabashima2003propagating . An explanation to this behavior can be found in bolthausen2014iterative , where a less trivial time index setting is shown to improve convergence properties outside the glassy phase. It was later realized that this time indexing emerges naturally when deriving AMP equations from BP equations, keeping track of BP time indexes. This approach requires considering the fully connected limit of BP equations, that are usually written on sparse graphs kabashima2001tap ; kabashima2003cdma ; donoho2009message ; krzakala2012probabilistic . To our surprise, the AMP equations exhibit convergence issues for the case of Hopfield model in the direct problem, when the initial condition is chosen at random. A different approach to derive AMP equations is proposed in opper2016theory . This new approach, in our case, does not improve the performances obtained by the naive time indexing Eq. 2. These issues are discussed in detail in appendix A. Thus, in the following, we mainly adopt Eq. 2, unless stated otherwise.

In Fig. 1, we present the teacher-student overlap $q=N^{-1}\sum_{i}m_{i}\zeta_{i}$ for $N=1000$. We observe that in the ferromagnetic-retrieval phase $\beta>1$, a perfect reconstruction may be realized already with a small number of samples. This is due to the large signal contained in the correlation matrix of the data $\overline{c}$. In particular, we notice that in the ferromagnetic phase the student is able to find a pattern correlated with the teacher’s one even at $M=1$. On the other hand in the paramagnetic phase the signal in $\overline{c}$ is weaker and reconstruction is possible only exploiting finite size effects, at the price of observing an extensive number of samples. As discussed in barra2018phase , the critical fraction $M/N$ of samples necessary for reconstruction corresponds to the paramagnetic-spin glass transition of the direct problem, as long as we restrict the analysis to the Bayes optimal scenario.

The $P>1$ case is more difficult because of the presence of the denominator in Eq. 1. However, we argue in appendix B that this term is effectively a (soft) orthogonality constraint over inferred patterns. This observation allows us to design an inference algorithm accordingly. First, let us construct a time evolution of the coupling matrix $J^{\tau}_{ij}$ with its initial condition given by $J^{\tau=0}_{ij}=N^{-1}\sum_{a}s_{i}^{a}s_{j}^{a}$. At each time step $\tau$, we consider $P^{\prime}$ TAP students trying to learn the teacher’s patterns independently. Namely, the magnetizations $m_{i}^{\mu}=\left<\xi_{i}^{\mu}\right>$ for each student evolve according to Eq. 2 from a randomly initialized configuration 111To escape the (unstable) fixed point at $m=0$, the absolute value of the local magnetization is chosen to be in the interval $[0.1,1]$. . Upon convergence, we evaluate the $P^{\prime}$ solutions with the score $S_{\mu}=\sum_{ij}\sum_{a=1}^{M}s^{a}_{i}s^{a}_{j}m^{\mu}_{i}m^{\mu}_{j}$. These scores characterize the quality of the TAP fixed points and we pick the one with the largest score 222This trick is closely reminiscent of the algorithm presented in franz2019fast , where the iterative steps are performed by evaluating the likelihood gain obtained moving in different directions and choosing the one with the largest payoff. The corresponding magnetization selected by this criterion at time $\tau$ are denoted by $\underline{m}^{\tau}$. Finally, in order to learn the remaining contributions, we remove the rank-1 part associated to the retrieved state $\underline{m}^{\tau}$ from the coupling matrix. When the student knows the actual number of patterns, this correspond to the rule $J^{\tau+1}_{ij}=J^{\tau}_{ij}-\gamma N^{-1}m^{\tau}_{i}m^{\tau}_{j}$, where $\gamma=M/P$ (assuming that different states are uniformly sampled in the dataset). We repeat these steps until no further patterns are found.

We stress that our algorithm finds solutions correlated with the patterns without any prior information, i.e. we start iterating the TAP equations from a random initial configuration. This is a rather remarkable property in comparison to the method used in braunstein2011inference , where BP equations were guided to converge to the fixed points associated with the patterns using a reinforcement term aligned with the magnetizations of the states. In Fig. 2 we compare the learned $P^{\prime}$ TAP fix points and the $P$ teacher’s patterns in a system of $N=1000$ with $P=20$ at $\beta=2$. In this regime data is generated in the retrieval phase mezard2017mean . The students observe $M=200$ samples, i.e. $10$ samples per state. We clearly see that all the $20$ patterns are successfully retrieved from the first 20 students. In addition, let us define two performance measures, the simplified likelihood and the reconstruction error

[TABLE]

and $\epsilon=[N(N-1)/2]^{-1}\sqrt{\sum_{i<j}(J^{r}_{ij}-J^{*}_{ij})^{2}/\sum_{i<j}(J^{*}_{ij})^{2}}$, where $J^{*}$ denotes the teacher’s coupling matrix, and $J^{r}$ is the inferred matrix at time $\tau$, $J^{r}_{ij}=N^{-1}\sum_{t=1}^{\tau}m_{i}^{t}m_{j}^{t}$. The simplified likelihood is defined by neglecting the partition function in Eq. 1. In Fig. 3 their behaviors are reported as a function of iteration time. As expected, $\epsilon$ decreases as the students learn the patterns but then increases when the students start to learn the remaining noise. Similarly, the simplified likelihood $\hat{\mathcal{L}}$ develops a kink at the point where students learn all the patterns, that can be used as a stopping condition of the algorithm. In Fig. 4 we study the behavior of $\epsilon$ for different values of the temperature. As a function of $\beta$, the system sweeps through different regions of the phase diagram. Data is generated with a sequential Glauber dynamics and states are equally sampled. In Fig. 4 we show the behavior of the error computed using the criterion discussed above with the number of observations in different regions of the phase diagram. As expected, perfect reconstruction is obtained only in the retrieval phase.

Another method to perform inference is the RBM, which is closely related to the Hopfield model barra2012equivalence ; barra2018phase ; mezard2017mean ; agliari2018non ; decelle2018thermodynamics . In fact, the posterior distribution of Eq. 1 can be recast to

[TABLE]

with $W_{i}^{\mu}=\sqrt{\frac{\beta}{N}}\xi^{\mu}_{i}$. At the same time, Eq. 4 defines an RBM with $N_{v}=N$ binary visible units, and $N_{h}=P$ Gaussian hidden units. Compared to existing methods sohl2011new ; aurell2012inverse ; nguyen2012mean ; cocco2012adaptive ; decelle2014pseudolikelihood ; lokhov2018optimal ; franz2019fast , the RBM is both time and space efficient as the number of parameters to be optimized scales as $N_{v}N_{h}$ rather than $N^{2}$. The number of hidden units $N_{h}$ plays the role of $P^{\prime}$. We consider the general setting $N_{h}\geq P$ in the following. Following the standard practice hinton2012practical , the weights $W_{i}^{\mu}$ are learned maximizing the log-likelihood using the PCD-$T$ algorithm, where PCD stands for Persistent Constrastive Divergence tieleman2008training and $T$ for the number of Monte Carlo steps used to estimate the part of the log-likelihood derivative involving the partition function. We used $T=10$. RBM learns a set of weights $J^{r}_{ij}=\beta^{-1}\sum_{\mu}W_{i}^{\mu}W_{j}^{\mu}$ that we can compare with the teacher coupling matrix. The error between $J^{r}_{ij}$ and $J^{*}_{ij}$ decreases during learning but it never achieves the values found with TAP. In order to monitor learning, we study the pseudo-likelihood $\mathcal{S}$ besag1977efficiency , i.e a proxy for the likelihood that can be easily computed. In our case, it is defined by $\mathcal{S}=\sum_{r=1}^{N_{v}}\mathcal{S}_{r}$, where

[TABLE]

More details on these points as well as on the implementations are given in appendix C. In Fig. 5 we show the behavior of these quantities for the dataset generated by the teacher at $\beta=2$, $N=1000$, $P=10$, when observing a different number $M$ of samples. An RBM with $N_{v}=1000$ visible units and $N_{h}=15$ hidden units is used. The minimum of $\epsilon$ is achieved when the pseudo-likelihood flattens. This happen when all of the relevant ($P=10$) modes of the data have been learned, as can be seen in Fig. 6. Unless learning starts in the vicinity of the teacher’s patterns, final RBM weights do not reproduce them, contrarily to the TAP-based algorithm discussed above. In fact, the Hopfield model is invariant under a rotation in the pattern space cocco2011high : the student RBM can learn, at most, the subspace spanned by teacher’s patterns. To prove it, we consider the Singular Value Decomposition (SVD) of the dataset, and the SVD of the weights. We denote by $\{\sigma^{\alpha}\}$ the singular values of the matrix $W^{\mu}_{i}$ and by $t^{\alpha}$ the error in reconstructing the singular vector of the data, indexed by $\alpha$, using only the singular vectors of the weight matrix. In Fig. 6, we show the emergence of different modes during learning. When the singular values $\sigma^{\alpha}$ of the coupling matrix emerge, the error $t^{\alpha}$ decreases. The first $P=10$ principal modes of the dataset are well represented by the subspace spanned by the singular vectors of the weight matrix $W$.

In summary, we discussed a new method to solve inverse problems with a clusterized dataset. We analyzed the fully connected Hopfield model in a teacher-student scenario and proposed an inference method based on the TAP equations working directly on the posterior distribution, i.e. the dual problem. We discussed a retrieval algorithm based on the parallel updating of the TAP equations with a naive indexing, showing that in our case it gives good results. Contrarily to previous methods, our algorithm is able at retrieving patterns, besides couplings, because TAP equations allows to reduces the continuous symmetry under rotation to a simple symmetry under permutation over the pattern labels. As a side result, we provide the analysis of the failure of AMP equations in our case, when iterated in a parallel manner starting from a random initial condition. Finally we compare these results with those obtained with RBM, exploiting their analogies with the Hopfield model. RBM is a good candidate model to perform inference with many variables, a task that would require a much longer execution time to methods based on the optimization of the pseudo-likelihood of an associate pairwise Ising model. Their ability to perform inference tasks systematically, as well as their performance on inferring sparse models, will be addressed elsewhere.

We thank S. Franz, F. Ricci Tersenghi and D. Saad for interesting discussions. J.R. and S.H. acknowledge the support of a grant from the Simons Foundation (No. 454941, Silvio Franz). D.T. acknowledges GNFM-Indam for financial support of the project ”Unsupervised Learning with Multi-Layer Boltzmann Machines”.

Appendix A linear stability analysis of TAP and AMP in the paramagnetic state

Here we present the linear stability analysis of TAP and AMP equations in the paramagnetic state. We will focus on the direct problem where $J_{ij}$ is constructed from $M$ random patterns. While the complete analysis is possible for arbitrary $\alpha$, we find it more instructive to focus on the limit $\alpha\to 0$ as it greatly simplifies the discussion. As will be shown below, our results are valid if $|1-\beta|\sim O(1)$, which is larger than $O(\alpha)$.

From now on, we denote by TAP the simple iterative updating scheme discussed in the text, reported here for convenience

[TABLE]

and by AMP the iterative scheme derived in [32],

[TABLE]

where $m^{t}_{i}=\tanh(\beta H_{i}^{t})$, $u^{t}=\beta(1-q^{t})$ and $Nq^{t}=\sum_{i}(m^{t}_{i})^{2}$. This time index setting naturally emerges from the expansion of the BP equations in the large connectivity limit [40].

Let us first consider the linear stability of Eq. 6. Near the paramagnetic state $M_{i}\sim 0$, this equation may be expanded into

[TABLE]

where the second term is neglected as it is of $O(\alpha)$. Performing the coordinate change with the eigenvectors of $J_{ij}$ as its basis, one obtains

[TABLE]

where $\lambda$ is an eigenvalue of $J_{ij}$. This implies that the paramagnetic solution becomes unstable when $\beta\lambda_{\text{max}}>1$. The spectrum of coupling matrix follows the Marchenko-Pastur law [32]. Namely, $P$ eigenvalues are $1+O(\sqrt{\alpha})$ and their eigenvectors span the same space spanned from the set of patterns. The remaining $N-P$ eigenvalues are zero. Thus we find that the critical temperature is $T_{c}=1+O(\sqrt{\alpha})$ (the true value is $T_{c}=1/(1+\sqrt{\alpha})$, found expanding TAP equation beyond the $\alpha\rightarrow 0$ limit, which is identical to the result of replica theory [32].

Similarly, AMP Eq. 7 can be expanded as follows:

[TABLE]

Because of the $\lambda-1$ term, in the limit $\alpha\rightarrow 0$, the $N-P$ eigenvalues equal to zero give the largest $O(1)$ contribution to the instability of the paramagnetic fixed point. In particular, the modes associated with patterns, with eigenvalues $1+O(\sqrt{\alpha})$, give a vanishing contribution. From the infinite temperature limit, the first $T$ where this equation becomes unstable is given by $-\frac{\beta}{1-\beta}=-1$, i.e. $T_{c}=1/2$. Nevertheless, this unstable direction is orthogonal to the patterns and the magnetization either converges to an unphysical state or never converge (see Fig. 7). The negative value of the leading eigenvalue for $\beta\in[1/2,1]$ leads to an oscillating behavior starting from the paramagnetic solution, as can be seen in the second plot in Fig. 7. Similar issues with parallel updating of the AMP equations were discussed in [46], and they can be alleviated by updating spins sequentially and introducing a strong dumping. Nevertheless their sequential updating leads to a much slower algorithm, without showing any improvement in the quality of inference in comparison to the parallel updating scheme of Eq. 6.

A different updating scheme of the TAP equations has been recently proposed by [54]. This approach does not require to consider the fully connected limit of the BP equations and it is suitable to be applied in systems with dense random coupling matrices. It is based on a dynamical mean field theory which allows to study the dynamics of iterative algorithms in the thermodynamic limit by averaging over the noise contained in the couplings. For the Hopfield model, the updating scheme turns out to be

[TABLE]

where $A_{t}=\beta/(1+\alpha u_{t})$ and $u_{t}$ is the same quantity introduced in the AMP Eq. 7. It is possible to see that this updating scheme does not present the issues of the AMP algorithm by repeating the same $\alpha\rightarrow 0$ analysis presented above.

In Fig. 8 we compare the performances of these three algorithms for different system sizes. We define $P_{c}$ as the probability to converge to one of the patterns of the system with overlap greater than $0.7$ when the initial condition is chosen at random. Sequential AMP were iterated with a damping term $d$, i.e. $m_{i}^{t+1}=(1-d)\tanh\beta H_{i}^{t+1}+dm_{i}^{t}$, and $d=0.95$. For the two parallel TAP equations, (Eq. 6-Eq. 12), the iteration is stopped when the average difference between $m_{i}^{t+1}$ and $m_{i}^{t}$ is smaller than $0.001$. For the AMP sequential algorithm the iteration is stopped when the average between $\tanh\beta H_{i}^{t+1}$ and $\tanh\beta H_{i}^{t}$ is smaller than $0.001$. In all the cases we observe that convergence to patterns is achieved in the retrieval phase. For small values of $N$, due to finite size effect, convergence regime extends in the metastable retrieval phase too.

The instability issue of the AMP equations presented above holds for the direct problem, but it can be extended also to the inverse, dual problem. In this last case, where $J_{ij}=N^{-1}\sum_{a=1}^{M}s^{a}_{i}s^{a}_{j}$, if there is enough signal in the data and $\lambda_{\text{max}}>2$, inference is possible also with parallel AMP equations. Nevertheless the analysis shows that obtaining time indexes from BP does not necessarily lead to good algorithms. TAP, as well as BP, equations describe only fixed points of the associated free energy and, in principle, any updating scheme could be used to solve these equations in an iterative manner, as shown in [54]. The relevance of this observation for other problems requires further analysis and, given that the AMP convergence issues are usually mitigated by considering a sequential updating with a strong dumping, it would be interesting to study whether a similar improvements is achieved when iterating TAP equations with the naive time indexing sequentially and with a strong damping, in problems where their parallel updating was failing.

Appendix B Posterior for $P>1$

We discuss the role of the difficult term arising in the posterior distribution when $P>1$. We show that for the simple case $P=2$, it has a clear interpretation in terms of a constraint on the orthogonality of the inferred patterns. In fact, let us consider

[TABLE]

and let us define $S=\{i:\xi^{1}_{i}=\xi^{2}_{i}\}$, such that $|S|=N(1+q)/2$, where $q$ is the mutual overlap between the two patterns, $Nq=\sum_{i}\xi^{1}_{i}\xi^{2}_{i}$. The exponent in Eq. 13 reads

[TABLE]

where we indicate with $\overline{S}$ the complement of set $S$. Using again the gauge transformation $s^{\prime}_{i}=\xi^{1}_{i}s_{i}$, Eq. 14 leads to

[TABLE]

where we indicate with $Z_{cw,\beta}^{N}$ the partition function of a ferromagnetic Curie-Weiss model at inverse temperature $\beta$. We observe that the interaction depends only on their mutual overlap. If we define $\phi=-N\log Z$, we obtain

[TABLE]

where $f_{cw}(\beta)$ is the free energy of the Curie-Weiss model at inverse temperature $\beta$. It is easy to check that $\phi(q)$ is a convex function with a minimum in $q=0$. Thus the term $-M\log Z(\underline{\xi})$ in the posterior can be interpreted as a soft regularizer for patterns orthogonality.

Appendix C Restricted Boltzmann Machine

A Restricted Boltzmann Machine (RBM) is a particular kind of Boltzmann machines in which units are divided in two layers, formed by visible $\{s_{i}\}$ and hidden $\{\lambda_{\mu}\}$ units, and only interactions $W^{\mu}_{i}$ between units of different layers are allowed, such that the proxy probability distribution reads

[TABLE]

where $E_{W}(\underline{s},\underline{\lambda})=-\sum_{i,\mu}W^{\mu}_{i}s_{i}\lambda_{\mu}$, and $Z(\{W\})$ is the partition function,

[TABLE]

For our purposes, $P(\underline{\lambda})$ denotes a generic distribution over hidden units, while visible units are $\pm 1$ binary variables. We indicate with $N_{v}$ the number of visible units and with $N_{h}$ the number of hidden units. RBM has the property that the two conditional probabilities, $P(\underline{s}|\underline{\lambda},\{W\})$ and $P(\underline{\lambda}|\underline{s},\{W\})$, factorize over the visible (resp. hidden) units. These machines are used to learn weights such that the distribution over the visible units reproduce the distribution of the data. In other words

[TABLE]

should reproduce as close as possible $P_{D}(\underline{s})=M^{-1}\sum_{a=1}^{M}\delta_{\underline{s},\underline{s}^{a}}\>.$ Weights can be found minimizing the KL distance between the two distribution, which is equivalent to maximizing the likelihood $\prod_{a=1}^{M}P(\underline{s}^{a}|\{W\})$ or the log-likelihood

[TABLE]

Optimal weights can be learned by gradient ascent:

[TABLE]

where the first average, usually referred to as positive phase, is

[TABLE]

and the second average, usually referred to as negative phase, is

[TABLE]

The second one is known to be difficult and it can be computed with approximate methods. One way to estimate it is to use a Monte Carlo (MC). Depending on the number of steps $T$ of the MC Markov chain, this method is referred to as CD-$T$, where CD stands for Contrastive Divergence. In the text, we discuss results obtained with $T=10$. When the positive term is computed over a sub set (mini-batch) of the dataset, the direction indicated by the gradient does not correspond to the correct one obtained considering the whole dataset. This trick introduces a source of randomness in the path to the solution, and the associate learning algorithm is called Stochastic Gradient ascent. In our experiments we use a mini-batch size equal to $100$. Since mini-batch samples are independent, different parallel MC can be used. In our experiments we used $100$ MC chains, one per mini-batch sample. Their initial conditions can be chosen to be the considered samples, but this quickly results in over fitting the parameters, since the MC dynamics spend all the time in the phase space regions close to the samples. When the initial condition of the MC dynamics is chosen at random and we keep track of their positions through different batches and epochs, this method is called Persistent CD (PCD). Our results are obtained using PCD.

As stated above, the likelihood function cannot be easily computed. Thus, we introduce the pseudo-likelihood that, for a model with hidden units, is defined by $\mathcal{S}=\sum_{r=1}^{N_{v}}\mathcal{S}_{r}$, where

[TABLE]

where the term inside the $\log$ is defined by

[TABLE]

and it is equal to

[TABLE]

where $\mathcal{N}_{a}$ is a sample dependent normalization factor,

[TABLE]

The Pseudo-likelihood is optimized by the same set of parameters $\{W\}$ that optimize the likelihood. In order to show this property, we can take derivatives of $\mathcal{S}$:

[TABLE]

The definition

[TABLE]

allows to write the first term of Eq. 28 as

[TABLE]

The second term is given by the average over samples of

[TABLE]

and similar manipulations on the second term lead to

[TABLE]

In the infinite sampling limit,

[TABLE]

In fact, it is easy to show that

[TABLE]

and, on the other hand, that

[TABLE]

Thus the gradient of the pseudo-likelihood $\mathcal{S}$ vanishes on the same set of parameters $\{W\}$ that solve $0=\partial_{W}\mathcal{L}$, and this is the reason the pseudo-likelihood can be used to control the learning state. In practice, the probability in Eq. 25, can be estimated after one step of Monte Carlo:

[TABLE]

where $\underline{\lambda}$ is sampled from the distribution $P(\underline{\lambda}|\underline{s}^{a},\{W\})$.

Finally, we discuss the learning of the RBM compared to our TAP based algorithm. As mentioned in the text, unless learning starts in the vicinity of the teacher’s patterns, final RBM weights do not reproduce them. In fact, the Hopfield model is invariant under a rotation in the pattern space. The dataset $\mathcal{D}$ analyzed by the student could have been produced by another set of patterns $\{\underline{\hat{\zeta}}^{\mu}\}_{\mu=1,\dots,P}$ given by $\hat{\zeta}^{\mu}_{i}=\sum_{\nu=1}^{P}O_{\mu,\nu}\zeta^{\nu}_{i}$ where $O$ is an orthogonal matrix. This symmetry implies that the student RBM cannot learn exactly the teacher’s patterns. One could think that the singular vectors of $W$ should learn at least the principal vectors of the data (that, given the spherical symmetry, are not necessarily aligned along the teacher’s patterns), as discussed in [59]. Nevertheless this is true only during the initial steps of learning, when couplings are small. This is reminiscent of the results discussed in [27], where the posterior of the problem is analyzed in a perturbative expansion. At the first order, corresponding to the small couplings regime, the student’s patterns are aligned along the singular vectors of the data at zero order. Anyway, computing higher order corrections, this relation breaks down.

In the following we show that RBM is learning the subspace spanned by the singular vectors of the data. To prove it, we consider the Singular Value Decomposition SVD of the dataset, $D=U_{D}\Sigma_{D}V_{D}^{T}$, where, considering $N<M$, $D$ is a $N\times M$ matrix, $U_{D}$ is a orthogonal $N\times N$ matrix, $\Sigma_{D}$ is a $N\times M$ matrix, with only $N$ diagonal elements different from zero, and $V_{D}$ is a orthogonal $M\times M$ matrix. $D$ represent the matrix of the dataset $\mathcal{D}$, where each column is a sample. Similarly, we consider the SVD of the weight matrix, $W^{T}=U_{W}\Sigma_{W}V_{W}^{T}$, where $W^{T}$ is $N_{v}\times N_{h}$, $U_{W}$ is a $N_{v}\times N_{v}$ orthogonal matrix, $\Sigma_{W}$ is a $N_{v}\times N_{h}$ matrix, with $N_{h}$ diagonal elements different from zero, and $V_{W}$ is a $N_{h}\times N_{h}$ orthogonal matrix. We consider $N_{v}=N$ and we decompose all of the data modes $u^{(\alpha),D}_{i}=[U_{D}]_{i\alpha}$ onto the subspace spanned by the first $N_{h}$ singular vectors of the weights, $u^{(\mu),W}_{i}=[U_{W}]_{i\mu}$, $\mu=1,\ldots,N_{h}$:

[TABLE]

where $\{\vec{u}^{(\mu),W}\}_{\mu=1,\ldots,N_{h}}$ are orthogonal vectors normalized to one. We measure the behavior of

[TABLE]

at the beginning and at the end of learning. This quantity measures the difference between the original vector and its projection onto the subspace spanned by the basis $\{\vec{u}^{(\mu),W}\}_{\mu=1,\dots N_{h}}$. The results of this analysis are found in the insets of Fig. 6, where we plot these quantities at the initial stage of learning and at the end.