A very fast inference algorithm for finite-dimensional spin glasses: Belief Propagation on the dual lattice.

A very fast inference algorithm for finite-dimensional spin glasses: Belief Propagation on the dual lattice.

Alejandro Lage-Castellanos, Roberto Mulet Department of Theoretical Physics, Physics Faculty, University of Havana, La Habana, CP 10400, Cuba.    Federico Ricci-Tersenghi Dipartimento di Fisica, INFN – Sezione di Roma 1 and CNR – IPCF, UOS di Roma,
Università La Sapienza, P.le A. Moro 5, 00185 Roma, Italy
   Tommaso Rizzo Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 5, 00185 Roma, Italy
July 3, 2019

Starting from a Cluster Variational Method, and inspired by the correctness of the paramagnetic Ansatz (at high temperatures in general, and at any temperature in the 2D Edwards-Anderson model) we propose a novel message passing algorithm — the Dual algorithm — to estimate the marginal probabilities of spin glasses on finite dimensional lattices. We show that in a wide range of temperatures our algorithm compares very well with Monte Carlo simulations, with the Double Loop algorithm and with exact calculation of the ground state of 2D systems with bimodal and Gaussian interactions. Moreover it is usually 100 times faster than other provably convergent methods, as the Double Loop algorithm.

I Introduction

Inference problems are common to almost any scientific discipline. Often an inference problem can be recast in that of computing some marginal probability on a small subset of variables, given a joint probability distribution over a large number of variables, . A typical example is the computation of the expectation value of a variable where is the single-variable marginal probability defined as

Clearly the computation of the sum in the r.h.s. is as difficult as the computation of the partition function , which is the main subject of statistical mechanics. In the general (and the most interesting) case, this problem can not be solved exactly in a time growing sub-exponentially with the size . We are deemed to use some kind of approximation in order to compute marginals in a time growing linearly in . The approximation schemes used so far are mainly adopted from the field of statistical mechanics libroMarcAndrea (), where mean-field-like approximations are standard and well controlled tools for approximating the free-energy, .

For non-disordered models, like the Ising ferromagnetic model, these approximations work quite well and provide good results pelizzola (). However, much less is known for models with disorder, and for this reason we will focus on spin glasses in the present paper. In disordered models one usually deals with an ensemble of problems (e.g. in spin glasses each sample has its own couplings, randomly chosen from a given distribution), and the results obtained by the statistical mechanics tools refer to average quantities, i.e. those of the typical samples. In other words one is not concerned with the behavior of a specific sample, but rather one looks at the whole ensemble. On the contrary when doing inference one is interested on the properties of a single specific problem and thus the above approximation schemes (based on statistical mechanics mean-field approaches) need to be converted into an algorithm that can be run on such a specific sample. Computing marginals on a given sample clearly gives more information than computing averages over the ensemble.

To our knowledge, effective (i.e. linear time in ) algorithms for computing marginals can be essentially divided in two broad classes: stochastic local search algorithms, that roughly sample the configurational space according to , and algorithms based on some kind of mean-field approximation. The former are exact on the long run, but the latter can be much more useful if an approximated answer is required in a short time. Unfortunately the latter also have some additional drawback due to the mean-field nature of the underlying approximation, e.g. the appearance of spurious phase transitions, that may prevent the proper convergence of the algorithm.

One more reason why this latter class of algorithms has a limited scope of application is that the convergence of the algorithm may strongly depend on the presence of short loops in the network of variables interactions. In this sense the successful application of one of these algorithms to models defined on regular lattices (which have many short loops) would be a major achievement.

In this paper we introduce a fast algorithm for computing marginals in 2D and 3D spin glass models libroYoung () defined by the Hamiltonian (further details are given below)


The first non trivial mean-field approximation for the above model corresponds to the Bethe-Peierls approximation scheme and the Belief Propagation (BP) algorithm kabaSaad (). Unfortunately when BP is run on a given spin glass sample defined on a dimensional lattice, it seems to provide exactly the same output as if it were run on a random regular graph with fixed degree : that is, for it converges to a solution with all null local marginals (), while for it does not converge to a fixed point.

The next step in the series of mean-field approximations (also known as Kikuchi approximations or ‘cluster variation method’) is to consider joint probability distributions of the four spins belonging to the same plaquette kikuchi (); pelizzola (). Under this approximation an algorithm has been derived which is called Generalized Belief Propagation (GBP) yedidia (). To our knowledge, this algorithm has been applied to 2D spin glasses, but only in presence of an external magnetic field, which is known to improve the convergence properties of GBP. Our experience says that, running plain GBP on a generic sample of a 2D spin glass without external field, a fixed point is reached only for high enough temperatures, well above the interesting region. The lack of convergence of GBP (and other similar message-passing algorithms, MPA) is a well known problem, whose solution is in general far from being understood. For this reason, a new class of algorithms has been recently introduced HAK03 (), which provably converge to a fixed point: these algorithms use a double loop iterative procedure (to be compared to the single loop in GBP) and for this reason they are usually quite slow.

The algorithm we are going to introduce is as fast as BP (which is the fastest algorithm presenty available), and converges in a wider range of temperatures than BP. Moreover the marginal probabilities provided by our algorithm are as accurate as those which can be obtained by Double Loop algorithms in a much larger running time.

Our algorithm works in the absence of external field, that is in the situation where MPA have more problems in converging to a fixed point. In this sense it is a very important extension to presently available inference algorithms.

The rest of the work is organized as follows. First, in section II we derive the GBP equations for the 2D Edward Anderson model. Then, in section II.1 we rewrite these equations in terms of fields, a notation that has a nicer physical interpretation and that we are going to use in the rest of the work. Section III presents the novel algorithm, inspired by the paramagnetic Ansatz to the GBP equations. In section IV we show the results of running this algorithm on the 2D Edward-Anderson model. There we compare its performance with Monte Carlo simulations, with the Double Loop algorithm and with exact calculation of the ground state of systems with bimodal and Gaussian interactions. Then, in section V we generalized our message passing equation to general dimensions and present some results for the 3D Edward-Anderson model. Finally, some conclusions are drawn in section VI.

Ii Generalized Belief Propagation on the 2D EA model

Here we present the GBP equations for the Edwards-Anderson (EA) model on a 2D square lattice, and we refer the reader to yedidia (); pelizzola () for a more general introduction. In our case (as well as in many other cases) GBP is equivalent to Kikuchi’s approximation, known as the Cluster Variational Method (CVM) kikuchi (). We will try a presentation as physical as possible.

Consider the 2D EA model consisting of a set of Ising spins located at the nodes of a 2D square lattice, interacting with a Hamiltonian


where the sum runs over all couples of neighboring spins (first neighbors on the lattice). The are the coupling constants between spins and are supposed to be fixed for any given instance of the model. If the interactions are not random variables, i.e. , then the 2D ferromagnet is recovered. We will focus on the two most common disorder distributions: bimodal interactions, , and Gaussian interactions .

The statistical mechanics of the EA model, at a given temperature , is given by the Gibbs-Boltzmann distribution

The direct computation of the partition function

or any marginal distribution is a time consuming task, unattainable in practice, since it involves the addition of terms, and therefore an approximation is required111The 2D case is special: indeed, thanks to the small genus topology, the partition function can be computed efficiently. However we are interested in developing an algorithm for the general case, and we will not make use of this peculiarity..

The idea of the Region Graph Approximation to Free Energy yedidia () is to replace the real distribution by a reduced set of its marginals. The hierarchy of approximations is given by the size of such marginals, starting with the set of all single spins marginals [mean-field], then following to all neighboring sites marginals [Bethe], then to all square plaquettes marginals , and so on. Since the only way of knowing such marginals exactly is the unattainable computation of , the method pretends to approximate them by a set of beliefs , , etc. obtained from the minimization of a region based free energy. In the Region Graph Approximation to Free Energy, a set of regions, i.e. sets of variables and their interactions, is defined, and a free energy is written in terms of the beliefs at each region. The Cluster Variational Method does a similar job, but instead of starting from an arbitrary choice of regions, it starts by defining the set of largest regions, and smaller regions are defined recursively by the intersections of bigger regions. In this sense, CVM is a specific choice of region graph approximation to the free energy.

For the 2D EA model, we will consider the expansion of the free energy in terms of the marginals at three levels of regions: single sites (or spins), links, and plaquettes. By plaquettes we mean the square basic cell of the 2D lattice. This choice of regions corresponds to the CVM having the square plaquettes as biggest regions. The free energy of the system is therefore written as

where runs over all regions considered, and the free energy in a particular region depends on the marginals at that level :

The symbol refers to the set of spins in region , while is the energy contribution in that region. The counting numbers (also Moebius coefficients) are needed to ensure that bigger regions do not over count the contribution in free energy of smaller regions, and follow the prescription


where is any region containing completely region , as e.g. a plaquette containing a link or a link containing a site. In the case of the 2D lattice, the biggest regions are the square plaquettes, and therefore , while the links regions have (as each of them is contained in 2 plaquettes regions), and finally the spins regions have (as each spin belongs to 4 links and 4 plaquettes). So the actual approximation for the EA model on a 2D square lattice is

where the sums run over all plaquettes, links and sites respectively. Please note that we are using the following notation for region indices: lower-case for sites, upper-case for links and upper-case calligraphic for plaquettes. The energy term in the sites contribution is only relevant when an external field acts over spins, and can be taken as zero in our case, since no external field is considered. Notice that whenever the interactions are included in more than one region (in our case are included in Link and Plaquette regions), the counting numbers guarantee that the exact thermodynamical energy is obtained when the beliefs are the exact marginals of the Boltzmann distribution. On the other hand, the entropy contribution is intrinsically approximated, since the cutoff in the regions sizes imposes a certain kind of factorization of in terms of its marginals (see pelizzola () for an explanation of the region graph approximation in terms of cumulants expansions of the entropy).

The next step in the method, is to compute the beliefs from the minimization condition of the free energy. However, an unrestricted minimization will generally produce inconsistent solutions, since the beliefs (marginals) are not independent, as they are related by the marginalization conditions


where and . In order to minimize under the constraints in Eq. (5) and under the normalization condition for each belief, a set of Lagrange multipliers should be added to the free energy in Eq. (II). There are different ways of choosing the Lagrange multipliers yedidia (), and each of them will produce a different set of self consistency equations. We choose the so called Parent to Child scheme (see section XX in yedidia ()), in which constraints in Eq. (5) are imposed by two sets of Lagrange multipliers: relating the belief at link to that at site , and relating the one at plaquette to the one at link .

With constraints (5) enforced by Lagrange multipliers, the free energy stationary conditions for the beliefs are the following:


where the notation refers to all links containing site and to all plaquettes containing link . The upper indices in the sums are written just to help understanding how many terms are in each sum for the 2D case. The precise meaning of the indices in each summation can be understood from the graphical representation in Figure 1. Lagrange multipliers are shown as arrows going from parent regions to children regions: simple arrows correspond to and triple arrows to . Let us consider, for instance, the belief in a link region , depicted in the central picture of Fig. 1: the sum of the two Lagrange multipliers corresponds to the triple arrows from plaquettes on the left and right of the link , while the double sum over the three and the three correspond to the three arrows acting over the two spins.

Figure 1: Schematic representation of belief equations (6). Lagrange multipliers are depicted as arrows, going from parent regions to children regions.

In Eq.(6), the are normalization constants, and the terms and are the corresponding energies in plaquettes and links, and are represented in Fig. 1 by bold lines (interactions) between circles (spins). In our case is zero since no external field is acting on the spins.

The Lagrange multipliers are fixed by the constraints they were supposed to enforce, Eq. (5), and they must satisfy the following set of self-consistency equations:

Again, to help understanding these equations, we provide in Fig. 2 their graphical representation. Note that there is one of these equations for every pair of Link-Site and every pair of Plaquette-Link in the graph. With we refer to interactions in plaquette that are not in link .

Figure 2: Message passing equations (II), shown schematically. Messages are depicted as arrows, going from parent regions to children regions. On any link , represented as bold lines between spins (circles), a Boltzmann factor exists. Dark circles represent spins to be traced over. Messages from plaquettes to links are represented by a triple arrow, because they can be written in terms of three parameters , and , defining the correlation and magnetizations and , respectively.

For each link in the 2D lattice, there are two link-to-site multipliers, and . For each plaquette there are four plaquette-to-link multipliers , corresponding to the four links contained inside the plaquette. Let be the number of spins in the lattice, there are links and plaquettes. So the originally intractable problem of computing marginals, has been replaced by the problem of solving a set of coupled equations for the Lagrange multipliers as those in Eq. (II). Once these equations are solved, the approximation for the marginals is obtained from Eq. (6) for the beliefs, and all thermodynamic quantities are derived from them as in Eq. (II).

Minimizing a region graph approximation to free energy, as that in Eq. (II) with constraints Eq. (5), or equivalently solving the set of self-consistent equations in Eq. (II), is still a non trivial task. Let us consider two ways of doing it. The first method is the “direct” minimization of the constrained free energy, using a Double Loop algorithm HAK03 (). This method is quite solid, since it guarantees convergence to an extremal point of the constrained free energy, but it may be very slow to converge. The second method, which is generally faster but is not guaranteed to converge, is the family of the so called Message-Passing algorithm (MPA), in which the Lagrange multipliers are interpreted as messages going from plaquettes to links, and messages from links to sites. Self consistency equations (II) can be viewed as the update rules for the messages in the left hand side, in terms of those in the right hand side. A random order updating of the messages in the graph by Eq. (II) (message passing) can reach a fixed point solution, and therefore, to an extremal point of the constrained free energy yedidia (). Next, we show explicitly how the message-passing equations looks like in terms of fields.

ii.1 From multipliers to fields

A particularly useful way of representing the multipliers (messages), with a nice physical interpretation, is the one used in kabashima (), which we adopt here. In full generality kabashima (); TLMF10 (), these multipliers can be written in terms of effective fields


In particular, the field corresponds to the cavity field in the Bethe approximation yedidia (). Using Lagrange multipliers, messages or fields, is essentially equivalent. We will often refer to fields as -messages to emphasize their role in a message-passing algorithm, and we will refer to self-consistency equations (II) as the message-passing equations.

This parametrization of the multipliers has proved useful to other endeavors, like the extension of the replica theory to general region graph approximations TLMF10 (). Here, all the relevant information in the Lagrange multipliers is translated to “effective fields” and . Notice that in this representation every single field corresponds to an arrow in the schematic messages-passing equations in Figure 2. In particular, the messages going from plaquettes to links are characterized by three fields , and the field acts as an effective interaction term, that adds directly to the energy terms in the Boltzmann factor. For instance, the first message-passing Eq. (II) is


where the indices refer to the notation used in Fig. 2 and is the interaction coupling constant between spins and . This equation naturally defines the updating rule for the message :



Note that the usual cavity equation for fields in the Bethe approximation basicBethe () is recovered if all contributions from plaquettes and are set to zero.

Working in a similar way for the second equation in (II) we end up with the updating rule for the message sent from any given plaquette region to one of its children links (see right picture in Fig. 2)



Equations (11) and (12) are equivalent to equations in (II), once multipliers (messages) are parametrized in terms of fields. For instance, note that the multipliers in the left hand side of second equation in (II) appear now subtracted in the right hand side of Eq. (12).

The field notation is more comprehensible and has a clear physical meaning. Each plaquette is telling its children links that they should add an effective interaction term to the real interaction , due to the fact that spins and are also interacting through the other three links in the plaquette . Fields act like magnetic fields upon spins, and the complete message is characterized by the triplet , and will be referred from now on as message. Furthermore, it is clear that some fields enter directly in the message-passing equations like and in Eq. (11) and and in Eq. (12). Also note that since our model has no external field, the fields break the symmetry of the original Hamiltonian whenever they are non zero. For instance, in the ferromagnetic model, when all , these fields are zero at high temperature and become non zero at Kikuchi’s critical temperature kikuchi (), implying a spontaneous magnetization in the ferromagnet.

Iii The Dual approximation for the paramagnetic phase

Unfortunately, the iterative message-passing algorithm for solving the GBP equations (11) and (12) often does not converge on finite dimensional lattices. While this is expected if long range correlations are present, it is rather disappointing that it happens also in the paramagnetic phase, where one would like to find easily the solution to the model. Here we are going to focus only on the paramagnetic phase, and propose an improved solving algorithm based on physical assumptions.

In the paramagnetic phase of any spin model defined by the Hamiltonian in Eq. (2), that is with no external field, variables have no polarization or magnetization: this in turn implies that in the solution all fields must be zero, and only fields should be fixed self-consistently to non zero values.

This paramagnetic solution has some interesting properties. First, it is always a solution of the GBP equations, since Eq. (11) and Eq. (12) are self-consistent with all . This means that starting from unbiased messages (all ) the iterative GBP algorithm keeps this property. Second, the paramagnetic Ansatz is correct, from the GBP perspective, at least at high enough temperatures, meaning that even if we start with biased messages (), the iterative algorithm converges to all at high temperatures. And last, but not least, the well studied physical behavior for the 2D EA model with zero-mean random interactions , is expected to remain always paramagnetic, i.e. to have no transition to a spin glass phase at any finite temperature JLMM (). Therefore, the Ansatz is both physically plausible and algorithmically desirable.

Figure 3: Message passing of correlation messages in the dual approximation. In the right hand side the trace is taken over the black spins.

Under the paramagnetic Ansatz, which we shall also call Dual approximation for a reason to be explained soon, the message-passing equation (11) is irrelevant, as it is always satisfied given that , while Eq. (12) now turns into (see Fig. 3)


The only relevant messages now are those associated to the multipliers , and they will be refereed to as messages. Eq. (13) can be interpreted as a correlation message-passing equation, giving the new interaction field that a certain link shall experience as a consequence of the correlations transmitted around the plaquette. The belief Eq. (6) also simplify. Obviously for every spin in the graph, and the link and plaquette beliefs are


The Dual algorithm we are proposing to study the paramagnetic phase of the EA model, is a standard message passing algorithm for the -messages, which works as follows.

1:  Start with all -messages null
2:  repeat
3:              Choose randomly one plaquette and one of its children links
4:              Update the field according to Eq. (13) as in Fig. 3
5:  until The last change for any -message is less than (we use typically )
6:  return  The beliefs defined in Eq. (14) for every pair of neighboring spins

Some damping factor can be added in the update step in order to help convergence.

iii.1 Mapping to the dual model

It is worth noticing that Eq. (13) is nothing but the BP equation for the corresponding dual model (hence the name of the algorithm).

The dual model has a binary variable on every link of the original model, and the original coupling constants play now the role of an external polarizing (eventually random) field

This Hamiltonian looks like the sum of independent variables, but this is not the case. The dual variables must satisfy a constraint for each cycle (or closed path) in the original graph, enforcing that their product along the cycle must be equal to 1. On a regular lattice any closed path can be expressed in terms of elementary cycles of 4 links (the plaquettes) and so it is enough to enforce the constraint on every plaquette: . The Gibbs-Boltzmann probability distribution for the dual model is then given by


where the product runs over all elementary plaquette.

The model described by the probability measure in Eq. (15) can be viewed as a constraint satisfaction problem with a non uniform prior (given by ). It is straightforward to derive the BP equations for such a problem. Indeed by defining the marginal for the variable on link in the presence of the solely neighboring plaquette as , the BP equations read


In the second summation the terms containing one or two variables sums to zero, while the other two terms are those written in the last expression. Equating the first and the last expressions, this equation is manifestly equal to Eq. (13).

iii.2 Average Case Solution

GBP in general, and the Dual approximation in particular, are methods for the study of the thermodynamical properties of a given problem. However, in the limit of large systems (, thermodynamical limit), we expect a typical behavior to arise. This is the so called self-averaging property of disordered systems. By typical we mean that almost every realization of the interactions will result in a system whose thermodynamical properties (free energy, energy, entropy) are very close to the average value.

Normally, in disordered systems, we cope with the limit and with the average over the random by the replica method. The application of the replica trick to regions graph approximations is a challenging task TLMF10 (). However, we can still grasp the average case behavior with a cavity average case solution of the dual message-passing equations, at the price of neglecting the local structure of the graph (beyond plaquettes).

The idea is to represent the set of -messages flowing in any given graph, by a population of messages . Then the message-passing Eq. (13) is used to obtain such population in a self-consistent way. More precisely, in every iteration three messages are randomly drawn from the population and a new message is computed by Eq. (13) using three couplings randomly selected from . The obtained message is put back into the population, and the iteration is repeated many times, until the population stabilizes.

Once we have the self consistent population of messages, we can compute the average energy


by a random sampling of the population and of the interactions. The average case solution is supposed to be very good whenever the network of interactions has no or few short loops. This is not the case in any finite dimensional lattice, since there the short loops (plaquettes) are abundant. Nonetheless, the average case solution gives a reasonably good approximation to the single instance results in 2D and 3D, as shown in the next Section.

Iv Results on 2D EA model

Message-passing algorithms work fine in the high temperature regime () of models defined on random topologies: this is the reason why these methods have been successfully applied in random constraint satisfaction problems, like random-SAT or random-Coloring KSAT (); Col (). However, when used on regular finite-dimensional lattices, they can experience difficulties even in the paramagnetic phase, because the presence of short loops spoils message-passing convergence.

Figure 4: Convergence probability of BP (Bethe approximation) and GBP on a 2D square lattice, as a function of inverse temperature. Data points are averages over 100 systems with random bimodal interactions. System sizes are with and a damping factor has been used in the iteration of the message-passing equations. The Bethe spin glass transition is expected to occur at () for a random graph with the same connectivity as the 2D square lattice. Notably, that temperature also marks the convergence threshold for BP equations in the 2D square lattice. GBP, on the contrary, reaches lower temperatures, but eventually stop converging.

It is well known that on a random graph of fixed degree (connectivity) the cavity approximation gives a paramagnetic result above (i.e. ) with all cavity fields . Below the Bethe critical temperature, this solution becomes unstable to perturbations, and we expect many solutions to appear with non trivial messages . The presence of many solutions in the messages passing equations is connected to the existence of many thermodynamical states in the Gibbs-Boltzmann measure, or, equivalently, to the presence of replica symmetry breaking. The appearance of such a spin glass phase is also responsible for the lack of convergence of message-passing equations, since the intrinsic locality of the message-passing equations fails to coordinate distant regions of the graph (which are now long-range correlated). As a consequence, the application of BP to the 2D EA model (that also has fixed degree ) still finds the paramagnetic phase at high temperatures, but below , the Bethe instability takes the message-passing iteration away from the solution and does not allow the messages to convergence to a fixed point (i.e. the algorithm wanders forever). In Figure 4 we show the convergence probability for the BP message-passing equations in the 2D EA model.

On the other hand, a straightforward GBP Parent-to-Child implementation does not fully overcome this problem. At high temperatures, the Parent-to-Child equations converge to a paramagnetic solution with all and non trivial , which turns out to be the same solution found by our Dual algorithm. When going down in temperature, the convergence properties of the algorithm worsen, and are sensitive to tricks like damping and bounding in the fields. A thorough discussion of these properties is left for a future work, but let us summarize that typically the algorithm stop converging at low temperatures, somewhere below , as shown in Figure 4.

So, in general, BP and GBP equations are not simple to use in finite dimensional systems at low enough temperatures: this warning was already reported in Refs. pelizzola (); HAK03 (); MooKap07 (). Indeed a different method for extremizing the constrained free energy named Double Loop algorithm HAK03 (); Yuille () was developed to overcome such difficulties. As mentioned earlier, Double Loop guarantees convergence of the beliefs, on any topology, with or without short loops. Given the convergence problems in GBP, researchers typically resort to Double Loop algorithms to extremize region graph approximations to the free energy, below the Bethe critical temperature.

In order to make a fair comparison with our Dual algorithm, we have used an optimized code for GBP and Double Loop algorithms: the open source LibDai library written in C++ libdai0.2.3 ().

The first interesting result of our work is that our Dual algorithm converges at all temperatures, just as Double Loop does. The reason why it converges is that there are no -messages, so the Bethe instability will not affect our message-passing iteration.

The second relevant result of our Dual algorithm, is the fact that it finds the same solution found by the Double Loop algorithm at all temperatures. In other words, the direct extremization of the region graph approximation to free energy Eq. (II) via a Double Loop algorithm finds a paramagnetic solution characterized by the beliefs and ; and the effective interactions found by the Double Loop algorithm are exactly equal to those found with our Dual algorithm, . This means that beliefs and correlations found by the two algorithms are identical: .

The third result is that the running times of our Dual algorithm are nearly four orders of magnitude smaller than those required by the Double Loop implementation in libdai, at least in a wide range of temperatures (see figure 5). More precisely, the convergence time of the Dual algorithm growth exponentially with , but still, in the relevant range of temperatures where the region graph approximation is a good approximation (not too low temperatures), the running time is always roughly a factor smaller than Double Loop.

Figure 5: Running times of the Double Loop algorithm HAK03 (); libdai0.2.3 () (libdai) and the Dual algorithm averaged over 10 realizations of a 2D EA model with Gaussian interactions. Generally the Double Loop algorithm requires a time 4 orders of magnitude larger than that used by the Dual algorithm. Three different precision goals where used for the Dual algorithm , while the precision of the Double Loop algorithm is . The inset shows the behavior of the running times for both algorithm versus the system size .

iv.1 Dual approximation vs Monte Carlo simulations

The fact that our Dual algorithm provides the same results (and much faster) than the Double Loop algorithm is a very good news. Essentially is telling us that we are not loosing anything by restricting the space of possible messages, as far as the region graph approximation is concerned. However, the ultimate comparison for the approximation has to be done with the exact marginals and correlations. In figure 6 we show a comparison between the exact correlations of neighboring spins obtained with a Parallel Tempering (PT) Monte Carlo simulation, and the Dual approximation estimate for the same two-spins correlations. The coincidence between and is essentially perfect at high temperatures, and it becomes weaker as the temperature is decreased. The reason for the discrepancies is obviously the fact that we are using an approximation in which collective behaviors of spins is accounted exactly only until the plaquette level; more distant correlations are approximated and these correlations become more important at low temperatures.

Figure 6: Comparison between the correlations obtained by the Dual algorithm and the nearly exact correlations obtained by a Parallel Tempering simulation. We used a EA model with Gaussian interactions. At each temperature the data correlation coefficient is reported.
Figure 7: Energy as a function of the inverse temperature for a 2D EA model, with both types of interactions, Gaussian and bimodal. Full lines represent the exact thermodynamical energy as obtained by a Monte Carlo simulation, points are the energies obtained under the Dual approximation, and dashed lines are the average case energies.

Given such a good correspondence between the correlations under the Dual approximation and the true correlations, we expect a very good estimate for the energy too. In Figure 7 we show with points the energy under the Dual approximation and with full lines the Monte Carlo exact energy: the data are indeed very close. The dashed lines show the average case energy for the Dual approximation, Eq. (17). In spite of the fact that the average case does not take into account the local structure of the lattice, the average case energy is quite close to the single instance one.

iv.2 Ground State Configuration in 2D

The good agreement between the correlations found by the Dual algorithm, and those found in a Monte Carlo simulation, for the 2D EA model, compels us to push this correspondence down to . More precisely, using the correlations obtained by our Dual algorithm at low temperatures, we try to compute a ground state configuration by the following procedure. The idea is to freeze iteratively the relative position of those interacting spins that are more strongly correlated, which is done by setting , and re-running the Dual algorithm until convergence every time one pair of spins is frozen. Note that freezing the relative position of spins is equivalent to freeze the dual variable . The freezing procedure is very simple, but for the fact one has to check that frozen links must be consistent with a spin configuration. More precisely, frozen variables must satisfy the requirement that on any closed loops the product is one,


For very short loops the satisfaction of this condition is automatically induced by the Dual algorithm: for example if three links on a plaquette freeze, the fourth link is immediately frozen to a value satisfy condition in Eq. (18). However, for longer loops (as the one shown in Fig. 8), the propagation of these constraints by the Dual algorithm is not perfect, since the information degrades with distance beyond the plaquette level. Then we need to enforce the constraints of Eq. (18) by a proper algorithm. At each stage of the freezing process, we define the clusters of frozen links as follows: if two frozen links share a spin, then they belong to the same cluster. In Figure 8 a cluster of frozen links is represented by bold lines. Notice that, once a spin is fixed in a cluster, all other spins are fixed as well by the frozen correlations. On the other hand, different clusters of spins can have arbitrary relative orientations.

Figure 8: Even if the link marked by the arrow is not the most polarized link according to the marginals provided by the Dual algorithm, the spins it connects are fully correlated by the fact that they belong to a cluster of frozen links (bold lines). Therefore, the marked link must be immediately fixed accordingly.

Consider now the situation depicted in Figure 8 and focus on the value of the correlation between the two spins connected by the link marked by an arrow. From the fact these two spins belong to the same cluster of frozen links (shown as bold link in Figure 8) we know they are perfectly correlated, however by running the Dual algorithm we could get a weak value for this correlation and then proceed by freezing a different link. A set of sub-optimal choices of this kind may finally produce a configuration of frozen links where the constraints in Eq. (18) are not all satisfied. In order to avoid these constraint violations we force any link whose spins are already part of the same cluster to be polarized accordingly. The freezing algorithm, therefore, works as follows.

1:  repeat
2:              run the Dual algorithm until convergence (at a low enough temperature)
3:              find the link with largest finite
4:              freeze that link by setting
5:              if link is connected to clusters and of frozen links then
6:                          merge clusters , and link in a unique cluster
7:              else
8:                          if link is connected to a single cluster of frozen links then
9:                                      add link to cluster
10:                          else
11:                                      create a new cluster with link
12:                          end if
13:              end if
14:              for all non-frozen links at the boundaries of a cluster of frozen links do
15:                          if link shares both spins with the same cluster then
16:                                      freeze link accordingly {To avoid violations of constraint}
17:                          end if
18:              end for
19:  until all links are frozen
20:  return  the spins configuration obtained by setting one spin and fixing the rest according to frozen links

The results obtained with this freezing procedure are quite good. In figure 9 we compare the resulting ground state energies with the exact solutions obtained using a web service running an exact solving algorithm juergengroundwww (). We used an ensemble of 100 EA models on the 2D square lattice with Gaussian interactions (so the ground state is not degenerate) and with bimodal interactions. Most points are along the bisecting line, meaning that the ground state found by either methods are the same. The relative error for the ground state energy is for Gaussian systems, and for bimodal systems. Looking at how many links are frustrated in one of the solutions and unfrustrated in the other, we found that the Dual+Freezing algorithm returns of the correct link correlation signs with respect to the true ground state solution for the Gaussian models. For the bimodal system, given the degeneracy of the ground state configuration, it is more probable to find the actual ground state energy but for the same reason the link overlap between the exact Ground State and the one found with Dual+Freezing is significantly lower ().

Figure 9: Correlation of the Dual+Freezing ground state energy with the exact ground state energy in systems. The top left points correspond to 100 bimodal systems , while the right bottom points correspond to 100 systems with Gaussian interactions. For bimodal interactions, the degeneracy of the ground state improves the probability of actually finding the correct ground state energy (), and conversely reduces the expected link correlation overlap with the exact ground state solution (). For Gaussian interactions, the ground state is not degenerated, and only in of the samples the Dual+Freezing method finds the actual ground state; however the average link overlap is very high (). The line is shown to guide the eye. Kindly note that two set of axes are being used.

In general we believe these results on the ground states to be very encouraging, considering that the Dual algorithm is very fast and not restricted to the 2D case (at variance to fast exact algorithms for computing ground states). They provide evidence that the marginals obtained by this Dual algorithm are reliable even at very low temperatures.

V Generalization to other dimensions

Let us now consider the region graph based approximation to the free energy for a generic -dimensional (hyper-)cubic lattice, using the same hierarchy of regions: square plaquettes, links and spins. After computing the counting numbers for a general -dimensional lattice, see Eq. (3), the free energy approximation becomes

Plaquettes are still the biggest regions considered at so have counting number 1, but now each link is contained in plaquettes, and each spin is in links and plaquettes. The message passing equations for the Dual algorithm in dimensions are then


where (resp.  and ) are the plaquettes containing the link (resp.  and ) excluding plaquette .

In the high temperature phase, this Dual approximation with all should be still a valid approach for any dimensionality . At low temperatures, however, the EA model in more than two dimensions have a spin glass phase transition and, therefore, we expect the Dual approximation to become poorer, as it can not account for a non trivial order parameter.

By running the Dual algorithm for the 3D EA model we have found a divergence of -fields around for bimodal couplings and around for Gaussian couplings. This divergence is due to the fact the -fields get too much self-reinforced under iteration. This divergence does not come as a surprise given that it happens also when studying the simpler pure ferromagnetic Ising model. However in the ferromagnetic model the temperature at which -fields diverge is always below the critical temperature and so the Dual algorithm still provides a very good description of the entire paramagnetic phase.

Unfortunately in the 3D EA the divergence of -fields takes place well above the critical temperature (which is for bimodal coupling and at for Gaussian couplings, see Ref. KKY, for a summary of critical temperatures in 3D spin glasses) and this would make the Dual algorithm of very little use. We have studied the origin of this divergence and we have found a general principle for reducing the divergence of -fields due to self-reinforcement, thus improving the convergence properties of the Dual algorithm. The idea is the following. When writing the Dual approximation as a constraint satisfaction problem with a non uniform prior, see Eq. (15), the constraints may be redundant. This is the case for the 3D cubic lattice: indeed both the number of links (i.e. variables in the dual problem) and the number of plaquettes (i.e. constraints in the dual problem) are . So, if constraints were independent, the entropy would be null at and negative for (and this is clearly absurd). The solution to the apparent paradox is that constraints are not independent: actually only of these are independent, and the remaining third is uniquely fixed by the value of the former. In this way the correct entropy is recovered at , given that a problem with unbiased binary variables subject to independent parity-check constraints has entropy . The dependence among constraints can be easily appreciated by looking at the 6 plaquette around a cube: if 5 of the 6 constraints are satisfied, then the sixth one is automatically satisfied and redundant.

The general rule for improving the convergence of the Dual algorithm is to remove redundant constraints (this principle is similar to the maxent-normal property of region based free energy approximations yedidia ()). Redundant constraints have no role in determining the fixed point values for the beliefs (since they are redundant), but during the iterations they provide larger fluctuations to messages and may be responsible for the lack of convergence. In practice, on a 3D cubic lattice, we may remove redundant constraints in many different ways: the basic rule states that one constraint (i.e. a plaquette) should be removed for each elementary cube, otherwise if a cube remains with its 6 plaquettes at least one redundant constraint will exist. We have used two different ways of removing one constraint per cube and we have found the same results in the entire paramagnetic phase. So, for simplicity, we are going to present data obtained by removing all constraints corresponding to plaquettes in the plane.

Figure 10: Running times of the Double Loop algorithm HAK03 (); libdai0.2.3 () (libdai) and of the Dual algorithm on a EA model with bimodal interactions (). The Dual algorithm is generally several orders of magnitude faster and returns the same solution as the Double Loop algorithm.
Figure 11: Comparison between the correlations obtained with the Dual algorithm and the (nearly) exact correlations obtained with a Parallel Tempering simulation in a 3D EA model of size with random bimodal interactions . At each temperature the correlation coefficient is reported. For the lowest temperature shown, , we also report the fraction w-s of pairs of spins such that .
Figure 12: The energy predicted by the Dual approximation in 3D EA model, compared to the average case energy, and the Monte Carlo simulation. We used a system with both types of random interactions, bimodal () and Gaussian distributed.

The Dual algorithm for the 3D EA model on the cubic lattice with no redundant constraints converges for any temperature above and so we can use it to study the entire paramagnetic phase. The lack of convergence deep in the spin glass phase is to be expected. Just as in the 2D case, the Dual algorithm (when converges) still finds the same solution obtained by a Double Loop algorithm, and again, it finds the solution nearly 100 times faster (see Fig. 10). Double Loop has the apparent advantage of converging at any temperature even at very low ones. However, deep in the spin glass phase, where the underlying paramagnetic approximation is clearly inaccurate, we believe that an algorithm (like the Dual one) that stops converging, is providing an important warning that something wrong is probably happening. Such a warning would be lacking by using a Double Loop algorithm.

In Figure 11 the correlations predicted by the Dual approximation, and those obtained by a Parallel Tempering Monte Carlo simulation are compared. At high temperatures the correspondence is quite good, but not as good as in 2D. However, it is important to stress that the 3D EA model is much more difficult to simulate than the 2D case: there is no exact method for computing the thermodynamics (at variance to the 2D case) and Monte Carlo methods require huge thermalization times, while the Dual algorithm runs in linear time with the system size.

In Figure 12 we show the estimates for the energy obtained from the Monte Carlo and the Dual algorithm (both on a single sample and on the average case). The very strong agreement between the Dual algorithm results on single samples and on the average case is telling us that -messages arriving at a given point on the lattice are uncorrelated to a very large extent. In other words, the effect of short loops in the lattice is not manifestly present in correlations between messages. On the contrary, the comparison between Dual algorithm results and Monte Carlo results is good only at high temperatures, and it degrades when approaching the critical temperature. This discrepancy can be understood as due to a growing correlation length in the EA model that diverges at the critical temperature: our Dual approximation does not account for correlations beyond the plaquette level and so it becomes inevitably poorer when the correlation length diverges. However, given the extremely fast converging times of the Dual algorithm, it can be viewed as a very effective algorithm for sampling the high temperature paramagnetic phase and as a reasonable approximate algorithm when approaching the critical point.

Vi Conclusions

We have introduced a novel Dual algorithm to compute marginals probabilities in the paramagnetic phase of frustrated spin models (e.g. spin glasses) on finite dimensional lattices. Inspired by the fact that in a paramagnetic phase with no external field each variable is unbiased (i.e. local magnetizations are null), the Dual algorithm is derived by adding such paramagnetic constraints in the GBP equations. While BP (i.e. Bethe approximation) and GBP algorithms have serious convergence problems at low temperatures even in the paramagnetic phase, the Dual algorithm converges very fast in a much wider range thanks to these constraints. The Dual algorithm can also be seen as BP on the dual lattice, where the interactions act as external fields on dual variables, thus improving convergence properties of the message passing algorithm.

We have tested the Dual algorithm for the Edwards Anderson spin glass model with bimodal and Gaussian couplings on 2D (square) and 3D (cubic) lattices. The results are very encouraging, showing convergence in the whole paramagnetic phase (and even slightly in the frozen phase for the 3D EA model) and comparing very well with exact correlations measured in Monte Carlo simulations. A comparison with a Double Loop algorithm (which is the state-of-the-art among general purpose inference algorithms) shows that both algorithms found the same result, but our Dual algorithm runs roughly 100 times faster. We also tried to push the Dual approximation to the limit, and we used the correlations inferred by the Dual algorithm to compute ground states configurations in the 2D EA model by a freezing procedure. Again, we showed that the ground states obtained in this way compare very well with exact computations.

The success of our proposal clearly shows that as long as variables are not long range correlated, the computation of correlations in a generic spin model can be done in a very fast way by means of message passing algorithms, based on mean-field like approximations. This kind of inference algorithms do not provide in general an exact answer (unless one uses it at very high temperatures or on locally tree-like topologies), and so they can not be seen as substitutes for a Monte Carlo (MC) sampling. However there are many situations where a fast and approximate answer is required more that a slow and exact answer. Let us just make a couple of examples of these situations. On the one side, if one need to sample from very noisy data, an approximated inference algorithm whose level of approximation is smaller than data uncertainty is as valid as a perfect MC sampler. On the other side, if one need to use the inferred correlations as input for a second algorithm (as for the freezing algorithm in Section IV.2) that will eventually modify/correct these correlations, a fast and reasonably good inference is enough.

The promising results shown in the present work naturally ask for an improvement in several directions. For example, in the paramagnetic phase of a model defined on a 3D lattice, our inference algorithm could be improved by using the cube as the elementary region, instead of the plaquette. An even more important improvement would be to extend the applicability range of the algorithm to the low temperature phase: but this requires a rather non trivial modification, since in low temperatures phase the assumption of zero local magnetizations needs to be broken.

F. Ricci-Tersenghi acknowledges financial support by the Italian Research Minister through the FIRB project RBFR086NN1 on “Inference and optimization in complex systems: from the thermodynamics of spin glasses to message passing algorithms”.


  • (1) M. Mézard and A. Montanari, Information, physics, and computation, Oxford University Press (2009).
  • (2) A. Pelizzola, J. Phys. A 38, R309 (2005).
  • (3) A.P. Young ed., Spin glasses and random fields, World Scientific (1998).
  • (4) Y. Kabashima and D. Saad, Europhys. Lett. 44, 668 (1998).
  • (5) R. Kikuchi, Phys. Rev. 81, 988 (1951).
  • (6) J. Yedidia, W. T. Freeman, and Y. Weiss, IEEE Transactions on Information Theory 51, 2282 (2005).
  • (7) T. Heskes, C. A. Albers, and H. J. Kappen, Proceedings of UAI-2003, 313 (2003).
  • (8) Y. Kabashima, J. Phys. Soc. Jpn. 74, 2133 (2005).
  • (9) T. Rizzo, A. Lage-Castellanos, R. Mulet, and F. Ricci-Tersenghi, J. Stat. Phys. 139, 375 (2010).
  • (10) M. Mézard and G. Parisi, Eur. Phys. J. B 20, 217 (2001); J. Stat. Phys. 111, 1 (2003). T. Castellani, F. Krzakala, and F. Ricci-Tersenghi, Eur. Phys. J. B 47, 99 (2005).
  • (11) T. Jorg, J. Lukic, E. Marinari, and O. Martin, Phys. Rev. Lett. 96, 237205 (2006).
  • (12) M. Mézard and R. Zecchina, Phys. Rev. E 66, 056126 (2002). A. Montanari, F. Ricci-Tersenghi, and G. Semerjian, J. Stat. Mech. P04004 (2008). F. Ricci-Tersenghi and G. Semerjian, J. Stat. Mech. P09001 (2009).
  • (13) R. Mulet, A. Pagnani, M. Weigt, and R. Zecchina, Phys. Rev. Lett. 89, 268701 (2002). L. Zdeborova and F. Krzakala, Phys. Rev. E 76, 031131 (2007).
  • (14) J. M. Mooij and H. J. Kappen, IEEE Transactions on Information Theory 53, 4422 (2007).
  • (15) Y. S.-K. Eye and A. L. Yuille, Neural Computation 14, 2002 (2001).
  • (16) J. M. Mooij, Journal of Machine Learning Research 11, 2169 (2010). http://www.libdai.org/
  • (17) http://www.informatik.uni-koeln.de/ls_juenger/research/sgs/index.html
  • (18) H. G. Katzgraber, M. Koerner, and A. P. Young, Phys. Rev. B 73, 224432 (2006).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description