{'fraction_non_alphanumeric': 0.028212227793160218, 'fraction_numerical': 0.006286013619696176, 'mean_word_length': 4.401576394502829, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An improved Belief Propagation algorithm finds many Bethe states in the random field Ising model on random graphs

We first present an empirical study of the Belief Propagation (BP) algorithm, when run on the random field Ising model defined on random regular graphs in the zero temperature limit. We introduce the notion of maximal solutions for the BP equations and we use them to fix a fraction of spins in their ground state configuration. At the phase transition point the fraction of unconstrained spins percolates and their number diverges with the system size. This in turn makes the associated optimization problem highly non trivial in the critical region. Using the bounds on the BP messages provided by the maximal solutions we design a new and very easy to implement BP scheme which is able to output a large number of stable fixed points. On one side this new algorithm is able to provide the minimum energy configuration with high probability in a competitive time. On the other side we found that the number of fixed points of the BP algorithm grows with the system size in the critical region. This unexpected feature poses new relevant questions on the physics of this class of models.

ferromagnetic models or planted models [3]), but for a general disordered model we are not aware of any algorithm that would output at least some of the states that characterize the Gibbs probability distribution in a given sample.
The same definition of a pure state in a disordered model is a delicate issue: indeed the standard trick of imposing the right boundary conditions is not easy to implement, not only because the number of different states may be very large, but also because for a generic graph the boundary is not well defined (think e.g. in random graphs or fully connected graphs).
Given a model of N interacting variables σ = {σ i } i=1,...,N , one would like to find the minimal decomposition of the Gibbs measure such that the measure µ α (σ) within state α has the clustering property, i.e. connected correlations decay fast enough at large distance, and thus distant variables are mostly independent (and the corresponding measure approximately factorizes).
In this work we consider models defined on sparse random graphs. These models can be solved exactly under the Bethe-Peierls approximation as long as correlations decay fast enough along the edges of the graph. Moreover they show more realistic physical properties, with respect to models defined on fully-connected network: e.g. the latter have interaction couplings scaling with an inverse power of the the system size and critical lines in the field vs. temperature plane diverging in the T → 0 limit.
For models defined on sparse random graphs, it is natural to identify the measure within a state µ α (σ) with a Bethe measure [2], i.e. a measure that can be, at least locally, factorized as where the first product runs over all edges of the graph (we assume for simplicity that variables interact only pairwise, the generalization to higher order interactions being straightforward) and the second is over all vertices, while µ α i and µ α ij are marginal probabilities over single variables and pairs of neighboring variables respectively. random graph can be expressed as the superposition of a relatively small number of Bethe states, that can be put in correspondence to BP fixed points.
The aim of the present work is to introduce a new algorithm which is able to identify many different Bethe states µ α in the Random Field Ising Model (RFIM) defined on a random graph at T = 0. Our new algorithm is based on the usual BP algorithm [2], which in the T = 0 limit is also known as max-sum algorithm. BP is an iterative algorithm whose fixed points correspond to minima of the Bethe free energy [5] and thus provides the local marginal probabilities µ α ij and µ α i within a Bethe state. In the case of the RFIM it has been recently shown that the global minimum of the Bethe free energy at T = 0 does actually correspond to the model ground state, irrespective of the graph which is defined on [6].
While running the standard BP algorithm on a given sample of the RFIM one usually reaches 1 or at most 2 fixed points, without any guarantee of having found the one of lowest free energy; our new algorithm achieves many different fixed points and the probability that the ground state is among these many fixed points turns out to be practically 1. Moreover, finding also the lowest excited states, our algorithm provides more physical information about the model than what can be extracted solely from the knowledge of the ground state [7].

II. THE MODEL
The random filed Ising model is well known to the statistical physics community as one of the simplest disordered systems [8]. It is defined by the following Hamiltonian where E and V are respectively the edge set and the vertex set of the interaction graph (e.g. a complete graph, a random graph, or a finite dimensional regular lattice). Variables s i = ±1 are N = |V | Ising spins and we choose to work with random fields {h i } extracted from a Gaussian distribution of zero mean and unitary variance, that ensures to have a second order phase transition in the whole field vs. temperature plane, including on the T = 0 axis. Varying J the model undergoes a second order phase transition at J c between a paramagnetic phase for J < J c and a ferromagnetic phase for J > J c . Other kinds of long range order, e.g. a spin glass phase, dominating the thermodynamics have been excluded [9,10] for the model in Eq. (2), that is the one with pairwise interactions between Ising spins (although replica symmetry breaking effects have been numerically detected in the p-spin version [11] and in the pairwise model with continuous spins [12]). The absence of replica symmetry breaking (RSB) for the states dominating the thermodynamics does not imply the absence of many metastable states [9], which may affect dramatically the dynamical behavior in the out of equilibrium regime [13].
Renormalization group studies indicate that thermal fluctuations are subdominant and the model can be studied directly at zero temperature [14]. This has a great numerical advantage since a min-cut algorithm exists, that provides the ground state (GS) configuration in polynomial time [15]. It is worth noticing, that this min-cut algorithm provides a single GS, even in case of strong GS degeneracy: in other words it is not a good sampler of the Gibbs measure, neither it can give any information about the gap between the GS and the first excited state.
In this work we study the RFIM on random 4-regular graphs (RRG). This ensemble consists of graphs where each vertex has exactly 4 neighbors randomly chosen (we avoid self-connections and double edges between a pair of nodes). These graphs are locally tree-like, in the sense that the local neighborhood of a randomly chosen vertex converges with high probability to a tree in the N → ∞ limit [2] and the typical size of loops is log(N ).
In the last two decades a lot of attention has been devoted to the RFIM defined on random regular graphs. Mainly the interest was focused on the non-equilibrium physics of the RFIM, as it emerged as a very effective model for Barkhausen noise and hysteresis in magnets (see [16] and references therein). Much efforts over the years have been made to characterize analytically the outof-equilibrium magnetization [16], correlation functions [17], hysteresis loops [18], and expansions toward fully connected models [19], just to mention few. The out-of-equilibrium Glauber dynamics too was solved in [20]. The vast number of metastable non-equilibrium states was numerically analyzed in [21], where different techniques were used to provide a clear description of the model complexity.
For what concerns the equilibrium physics of the model, beside the classic paper [22] where the bimodal random fields version was studied, the only modern analytic study of the RFIM on a random regular graph is the one present in [23], where the model is studied at zero and finite temperature using the cavity method approach [24,25]. Moreover on such graphs the loops corrections can be handled analytically [7,26,27].

III. BELIEF PROPAGATION EQUATIONS AND MAXIMAL SOLUTIONS
Belief Propagation (BP) is an iterative message-passing algorithm for solving the self-consistency equations that determines the minima of the Bethe free-energy [28]: these minima do actually corresponds to the physical states of the system and the one with lowest free-energy is the one dominating the thermodynamics. It is known that the Bethe approximation is exact on trees [29,30] as it relies on the independence of the neighborhood of a given spin when this spin is removed from the graph (hence the name cavity method in physics). When the graph has loops the width of the approximation made is generally unclear, thou applying BP to graphs with loops is a consolidate practice, specially in computer science studies: in this case is called loopy belief propagation. On random graphs, which are locally tree-like, the failure of BP is directly related to the birth of long range correlations and thus phase transitions.
In this context the RFIM presents a lot of features that makes it a very interesting model to approach with the BP algorithm. According to results in [31] the model is replica symmetric, and thus BP should provide the exact thermodynamics in this case, a conjecture that was recently proven in [32]. However one must be aware of the fact that these conclusions only apply to the case where only one BP fixed point exists, a condition that is far from being satisfied by most of the models that have been studied with the BP algorithm. And indeed we shall see that for the RFIM has more than one fixed point. Nevertheless in [6] it was shown that the zero temperature global minimum of the Bethe free energy coincide with the ground state of the problem, for every graph topology. Strictly speaking, this means that among the (possibly many) fixed points of the BP algorithm there must be the GS solution. This peculiar feature of the RFIM stimulated a lot of research tending to systematize BP as a provably exact ground state solver for the problems with an energy function as the one in Eq. (2), see for example [33]. A nice example with this taste is the tree-reweighted message passing scheme [34,35]. A very interesting recent result is present in [36] where it is demonstrated that with a proper calibration BP can be made equivalent to the min-cut/max-flow algorithm, which is an exact solver.
All these results points out that indeed BP can be good even if loops are present. However to our knowledge an empirical study of the real performance of BP on single instances is lacking. The main advantage of a direct approach, as we will see, is to face the presence of many fixed points.
On one side this obviates most of our analytic understanding of BP. On the other side we shall see that the organization of these BP solutions seems to be relevant for the physics of the RFIM.
Given an energy function of the type (2) with pairwise interactions and binary variables, the BP update rules reads (see [2]): for every edge (i, j) ∈ E. By taking the limit T → 0 we obtain: These equations are also known as min-sum equations. Once an initial value is assigned to the i→j }, we let them evolve with the update rule (4) until a stopping criterion is met. In this work we use the following prescription: if at a time t * , max (i,j)∈E |h i→j | < for some accuracy , then the messages have converged to a fixed point. We found that such a procedure is highly stable and we used = 10 −10 . Once the equations have converged, one can easily compute the spins value with the fixed point messages via from which physical observable can be computed.
Given that the RFIM at zero temperature displays a second order phase transition at a critical value of the coupling strength J c = 0.54404(5), one expects to find one stable fixed point in the paramagnetic region J < J c , and two stable fixed points in the ferromagnetic regime J > J c , where the one with the magnetization of the same sign of the sum of the random fields will be the dominating one, i.e. the one with lower energy [7] . As we will see this picture is over-simplified, and turns out that in a broad region near the critical point there is a large number of stable fixed points, depending on the disorder realization.
From physical considerations, a natural choice for the initial condition (IC) of the BP messages is to bias them towards the sign of the sum of local random fields. Specifically we start with the two IC: In this way the (+) IC is equivalent to set all the messages equal to their maximum value with the same sign of the sum of the random fields of the specific instance, and analogously for the (−) IC.
This choice is motivated from the following observation: when J is small we expect that all the spins will prefer to align to their local random field, so that only one fixed point exists with zero average magnetization. In this situation the particular choice of the IC should play no role. On the other side, when the coupling constant is large enough the spins will tend to align all together in the same direction, irrespective of the random fields, so we expect the (+) IC to converge on the ferromagnetic minimum with lower energy and the (−) IC to converge on the sub dominating one. As we will see in a while this intuition turns out to be true only away from the critical region.
As the model is replica symmetric, for all J we expect the BP algorithm to converge for every IC. In Fig. 1 we report the number of BP steps (where a single step corresponds to the update of all the 2|E| messages) to reach the convergence criterion, t conv , as a function of the interaction strength for the (+) IC and we tested that the conclusion is the same for other choices for the IC.
As expected, the BP dynamics gets slower in the critical region and the convergence time seems to diverge at the critical point as t conv ∼ N 1/3 . In this respect the BP algorithm is competitive with the latest versions of the min-cut algorithm [37], at least on locally tree-like graphs.
However the BP algorithm is not guaranteed to converge on the minimum energy configuration (the GS) of a given graph. Although we know that the GS of the model should be a fixed point of the BP algorithm [6] it is not obvious at all how to initialize the messages in such a way to converge on the GS. Moreover we have no way of asserting that a fixed point configuration is the one with the minimum possible energy, although in the next section we shall give a criterion valid at least in the paramagnetic phase.
In order to understand the behavior of BP under the (±) IC we measured the probability that it converges on the GS configuration of the problem, obtained with the min-cut algorithm. The results are reported in Fig. 2 for graphs of sizes ranging from 2 11 to 2 16 . The results are averaged over a number of realization for the disorder (graphs and random fields) that goes from 5·10 4 for the smallest size to 5 · 10 3 for the biggest one. As expected, the probability that the (+) IC converges to the GS tends to one in the paramagnetic phase (J < J c ). In the ferromagnetic phase we found that there is a small but finite probability that the dominating minimum is not reached with the (+) IC. This happens when the two IC converge on two configurations with the magnetization of the same sign, where the one with the smallest modulus is reached with the (−) IC and has the lower energy. Beside this phenomenon we can conclude that as far as we are away from the critical region the two IC (±) guarantee the convergence on the lowest energy configuration.
The situation drastically changes in the neighborhood of J c . As the reader can see in the right panel of Fig. 2 the probability that none of the (±) IC converges on the GS is non-zero and grows with N at J c . As we will see in more detail in Section V, in the critical region many solutions of the BP equations appear and a level crossing phenomenon is at work, such that the GS can not be obtained by starting from an IC that corresponds to the GS at a nearby value of J. Beside posing some question marks on the real nature of the phase transition, this picture makes the optimization problem very hard to solve with the standard BP algorithm. Although we know that the iterative BP equations must converge for every initial condition, we do not have any suggestion on how to initialize the BP messages in order to converge to the GS with high probability very close to J c .
We tested that the remaining trivial initial condition, i.e. the one with all the messages equal to zero, does not increases significantly the probability of success as it almost always converge on one of the two (±) fixed points. The same holds for random initializations that achieve the GS with a smaller probability that starting from (±) IC.

IV. PERCOLATION OF FROZEN VARIABLES
In this section we prove some properties of the maximal solutions, i.e. the fixed point solutions obtained with the (±) IC defined in Eq. (7). These results are preliminary to the definition of the algorithm that finds many different BP fixed points and rely on the following no-passing rule (NPR), first introduced by Middleton [38] in the context of charge density waves and later extended to the Glauber dynamics [16] and to the GS evolution [39] in the zero temperature RFIM.
Let us adopt the convention that two vectors v and v are partially ordered (to be indicated Thus, if the initial messages satisfy then the same order must hold between messages at any time t > 0. Moreover when the fixed point of BP is reached, every spin will be computed as a non-decreasing function of the fixed points messages, such that the spin configurations will satisfy s (A) s (B) .

From Eq. (4) it is immediate to derive the following bound on the BP messages
The (±) initial conditions, see Eq. (7), do correspond to set the BP messages all equal and taking the largest (positive or negative) value allowed by the above bound. This in turn implies that at any time on every edge of the graph a BP message cannot assume a value greater/lower than the value of the corresponding (+)/(−) message for every initial condition ( * ). Moreover the two initial condition (±) will converge on the fixed points whose configurations are the one with the lowest and highest magnetization, as this is due to the NPR. For this reason we shall call them the maximal solutions.
then spin s i must take the same value in all BP fixed points (i.e. in all the free-energy minima, including the GS) and we call it a frozen spin.
Thanks to this simple property we can claim to have found the GS in case the two maximal The mean fraction of free spins is shown in Fig. 3 for a random 4-regular graph.
The data in Fig. 3 strongly suggest the presence of a phase transition in n free exactly at the critical point J c . Moreover for J > J c the data seem to approach in the large size limit the

V. EXPLORATION OF THE BETHE FREE ENERGY LANDSCAPE
We present here a heuristic modification of the BP algorithm that is able to find several new BP fixed points (i.e. Bethe free energy minima) by relying on BP fixed points already found. We always start with the two maximal IC (±) defined in Eq. (7), that converge to maximal solutions.
If the maximal solutions do coincide, then they provide the GS, and BP has no other fixed points: thus the algorithm stops. On the contrary, if maximal solutions are different, the frozen spins and the corresponding messages are fixed for the rest of the run, and the algorithm proceeds looking for more BP fixed points. We saw in Section III (see right panel in Fig. 2) that in the critical region the GS does not always correspond to a maximal solution, so it worth continuing the search for the possible GS.
The key question is how to initialize the non-frozen BP messages in order to find new BP fixed points (and hopefully the GS) without wasting time in random initializations. Since the fixed point already found may have a rather large basin of attraction under the BP iteration, a reasonable initialization that is more likely to flow to a different fixed point (if any) is the following: In this way if u The mean number of BP fixed points (i.e. solutions to BP equations) N sol that our algorithm outputs is reported in Fig. 4. Far from the critical point, for |J − J c | 1, typically the algorithm does not provide additional information than running BP with the maximal initial conditions (see  Table I we report such a number of samples, that are statistically compatible with a probability of missing the true GS independent of N and close to 3 · 10 −5 . In the very rare samples where our algorithm does not find the true GS, the lowest free-energy minimum found lays above the true GS by O(10 −4 ). Moreover we believe that with a proper calibration of the algorithm (e.g. by using some appropriate damping) the probability of finding the GS can be made even higher. Here, however, our primary interest is on exploring efficiently the Bethe free-energy minima and so we are not going deeper with the possible use of our improved BP algorithm for solving the associated optimization problem, thou we believe that this could be a promising direction to follow.
At this point a comment is mandatory. It is worth stressing that we have no guarantee that our improved BP algorithm finds all the BP fixed points. Moreover, it was demonstrated that while the BP fixed points are Bethe free energy minima, the converse needs not to be true [40]. That said, the high probability with which we find the GS for each sample makes us confident that with The best fit exponent γ increases with N and goes above 1 for the largest N values; however a precise extrapolation to the N → ∞ is difficult. We have estimated the value γ ∼ 3/2 (used to rescale data in the right panel of Fig. 5 by studying N sol ∝ ξ 2−γ N versus N 2 sol / N sol ∝ ξ N .

VI. CONCLUSIONS
We have presented an improved BP algorithm that is able to find many stable solutions to BP equations (i.e. free-energy minima or Bethe states) in the RFIM defined on a random graph. While the standard implementation of BP is effective only far from the critical region, close to the critical point the choice of the initial condition plays a crucial role: indeed the number of free-energy minima grows and their basins of attraction shrink, such that a randomly chosen initial condition is unlikely to find the lowest free-energy minima. To partially overcome this problem, we have proposed a new way of recursively initialize BP with a proper interpolation of BP fixed points already found. This algorithm returns the ground state with high probability, together with many other BP fixed points, corresponding to metastable states of low free-energy.
A careful analysis of the number of BP fixed points in the critical region reveals a slow divergence of its mean value N sol ∼ log(N ) and a probability distribution decaying, in the large N limit, sol . The existence of a diverging number of states is a key feature of disordered systems (e.g. is at the heart of the replica symmetry breaking theory [1]). However, to the best of our knowledge, this is the first case where many different Bethe states are explicitly found in a model. It is still unclear what are the physical consequences of the existence of a diverging number of Bethe states in the critical region of the RFIM; although the thermodynamics of the model can be solved within a replica symmetric ansatz [9], one is tempting to interpret the occurrence a diverging number of Bethe states as a reminiscence of a replica symmetry breaking phase.
It is important to remark that the Bethe states found by our improved BP algorithm are much more stable than typical minima of the energy potential function; the latter are called one-spinflip stable configurations and are found in a much broader range [41,42], however their relevance for the out-of-equilibrium dynamics and in general for the determination of physical properties of the the model is unclear. On the contrary, it has been shown [43] that BP fixed points at zero temperature correspond to configurations that are stable under the flip of spins in any subset of vertices forming a subgraph with at most one loop. This property of BP fixed points is known as maximal neighborhood stability property and makes them much better candidates to describe also the physics at finite temperature and the behavior of thermal algorithm, that may get trapped in these Bethe states even if their evolution is stochastic.
It has been recently shown rigorously [4] that the Gibbs measure of any random graphical models can be decompose into a moderate number, e.g. O(log(N )), of Bethe states, and that the probability marginals in these Bethe states can be obtained from the corresponding BP fixed point. This result implies that, being able to find the relevant BP fixed points, one could in principle compute exactly any marginal in a random graphical model. Our improved BP algorithm allows one to perform this program for the case of the RFIM on a random graph.
In the light of the present results, supporting the existence of a large number of BP fixed points in the RFIM on a random graph, we believe it would be useful to reconsider the proofs that were derived assuming the existence of a unique BP fixed point [32].
Let us conclude with a remark on the optimization problem of finding the lowest energy configuration. For the case of the RFIM on a random graph, our improved BP algorithm finds almost certainly the ground state in a time which is competitive with algorithms, like min-cut, that provably return the exact ground state. However there are applications where having at hand many low energy configurations, as those retuned by our improved BP algorithm, allows for a much better final choice [44]. Moreover this improved algorithm, finding several different low-energy configurations, allows one to study the fundamental excitations of the system, without the use of methods, like e.g. the -coupling algorithm [45], that require to modify the Hamiltonian with somehow ad-hoc and not fully justified perturbations. Finally the algorithm presented here is quite robust with respect to small changes in the Hamiltonian (e.g. the introduction of a small fraction of negative coupling), while the min-cut algorithm can not be used as soon as any small amount of frustration is introduced in the interaction couplings.