• (i) k-satisfiability (k-SAT) with k ≥ 3. In this case = {0, 1}. The constraints are defined by fixing a k-tuple [z_a(1), …, z_a(k)] for each a, and setting ψ_a(x_ia(1), …, x_ia(k)) = 0 if (x_ia(1), …, x_ia(k) = (z_a(1), …, z_a(k)) and = 1 otherwise.
• (ii) q-coloring (q-COL) with q ≥ 3. Given a graph G with N vertices and M edges, one is asked to assign colors x_i {1, …, q} to the vertices in such a way that no edge has the same color at both ends.

The optimization (maximize the number of satisfied constraints) and counting (count the number of satisfying assignments) versions of this problems are defined straightforwardly. It is also convenient to represent CSP instances as factor graphs (2), i.e., bipartite graphs with vertex sets [N], [M] including an edge between node i [N] and a [M] if, and only if, the ith variable is involved in the ath constraint (compare Fig. 1). This representation allows one to define naturally a distance d(i, j) between variable nodes.

Fig. 1. The factor graph of a small CSP allows to define the distance d(i, j) between variables xi and xj (filled squares are constraints and empty circles variables). Here, for instance, d(6, 1) = 2 and d(3, 5) = 1.

Ensembles of random CSPs (rCSPs) were introduced (see e.g., ref. 3) with the hope of discovering generic mathematical phenomena that could be exploited in the design of efficient algorithms. Indeed several search heuristics, such as Walk-SAT (4) and “myopic” algorithms (5) have been successfully analyzed and optimized over rCSP ensembles. The most spectacular advance in this direction has probably been the introduction of a new and powerful message passing algorithm (survey propagation, SP) (6). The original justification for SP was based on the (nonrigorous) cavity method from spin glass theory. Subsequent work proved that standard message passing algorithms (such as belief propagation, BP) can indeed be useful for some CSPs (7–9). Nevertheless, the fundamental reason for the (empirical) superiority of SP in this context remains to be understood and is a major open problem in the field. Building on a refined picture of the solution set of rCSP, this article provides a possible (and testable) explanation. We consider two ensembles that have attracted the majority of work in the field: (i) random k-SAT: each k-SAT instance with N variables and M = Nα clauses is considered with the same probability; (ii) q-COL on random graphs: the graph G is uniformly random among the ones over N vertices, with uniform degree l (the number of constraints is therefore M = Nl/2).

It is well known that rCSPs may undergo phase transitions as the number of constraints per variable α is varied.g The best known of such phase transitions is the SAT-UNSAT one: as α crosses a critical value α_s(k) (that can, in principle, depend on N), the instances pass from being satisfiable to unsatisfiable with high probabilityh (10). For k-SAT, it is known that α_s(2) = 1. A conjecture based on the cavity method was put forward in ref. 6 for all k ≥ 3 that implied in particular the values presented in Table 1 and α_s(k) = 2^k log 2 − (1 + log 2)/2 + O (2^−k) for large k (11). Subsequently, it was proved that α_s(k) = 2^k log 2 − O(k), confirming this asymptotic behavior (12). An analogous conjecture for q-COL was proposed in ref. 13 yielding, for regular random graphs (14), the values reported in Table 1 and l_s(q) = 2q log q − log q − 1 + o(1) for large q [according to our convention, random graphs are with high probability uncolorable if l ≥ l_s(q)]. It was proved in refs. 12 and 15 that l_s(q) = 2q log q − O(log q).

Table 1. Critical connectivities for the dynamic, condensation, and satisfiability transitions in k-SAT and q-COL

Even more interesting and challenging are phase transitions in the structure of the set ^N of solutions of rCSP's (“structural” phase transitions). Assuming the existence of solutions, a convenient way of describing is to introduce the uniform measure over solutions μ(x):

We shall now introduce a few possible “global” characterizations of the measure μ. Each one of these properties has its counterpart in the theory of Gibbs measures, and we shall partially adopt that terminology here (17).

To define the first of such characterizations, we let i [N] be a uniformly random variable index, denote as x the vector of variables whose distance from i is at least , and by μ(x_i|x) the marginal distribution of x_i given x. Then we say that the measure (Eq. 1) satisfies the uniqueness condition if for any given i [N],

Notice that the uniqueness threshold is far below the SAT-UNSAT threshold. Furthermore, several empirical studies (19, 20) pointed out that BP [as well as many other heuristics (4, 5)] is effective up to much larger values of the clause density. In a remarkable series of papers (6, 21), statistical physicists argued that a second structural phase transition is more relevant than the uniqueness one. Following this literature, we shall refer to this as the “dynamic phase transition” and denote the corresponding threshold as α_d(k) [or l_d(q)]. To precise this notion, we provide here two alternative formulations corresponding to two distinct intuitions. According to the first one, above α_d(k) the variables (x₁, …, x_N) become globally correlated under μ. The criterion in 2 is replaced by one in which far apart variables x are themselves sampled from μ (“extremality” condition):

According to the second intuition, above α_d(k), the measure (Eq. 1) decomposes into a large number of disconnected “clusters.” This means that there exists a partition {A_n}_{n=1 …} of ^N (depending on the instance) such that: (i) one cannot find n such that μ(A_n) → 1; (ii) denoting by _ϵA the set of configurations x ^N\A whose Hamming distance from A is at most Nϵ, we have μ(_ϵA_n)/μ(A_n)(1 − μ(A_n)) → 0 exponentially fast in N for all n and ϵ small enough. Notice that the measure μ can be decomposed as:

The above ideas are obviously related to the performance of algorithms. For instance, the correlation decay condition in Eq. 3 is likely to be sufficient for approximate correctness of BP on random formulae. Also, the existence of partitions as above implies exponential slowing down in a large class of Monte Carlo Markov chain sampling algorithms.i

Recently, some important rigorous results were obtained supporting this picture (22, 23). However, even at the heuristic level, several crucial questions remain open. The most important concern the distribution of the weights {w_n}: are they tightly concentrated (on an appropriate scale) or not? A (somewhat surprisingly) related question is: can the absence of decorrelation above α_d(k) be detected by probing a subset of variables bounded in N?

SP (6) can be thought as an inference algorithm for a modified graphical model that gives unit weight to each cluster (20, 24), thus tilting the original measure toward small clusters. The resulting performances will strongly depend on the distribution of the cluster sizes w_n. Further, under the tilted measure, α_d(k) is underestimated because small clusters have a larger impact. The correct value was never determined (but see ref. 16 for coloring). Mézard et al. (25) undertook the technically challenging task of determining the cluster size distribution without, however, clarifying several of its properties.

In this article we address these issues and unveil at least two unexpected phenomena. Our results are described in Results and Discussion and summarized in Fig. 2. Finally, we discuss the connection with the performances of SP. Some technical details of the calculation are collected in Cavity Formalism, Tree Reconstruction, and SP.

Fig. 2. Pictorial representation of the different phase transitions in the set of solutions of a rCSP. At αd,+ some clusters appear, but for αd,+ < α < αd they comprise only an exponentially small fraction of solutions. (more ...)

The formulation in terms of extremality condition (see Eq. 3) allows for an heuristic calculation of the dynamic threshold α_d(k). Previous attempts were based instead on the cavity method, which is an heuristic implementation of the definition in terms of pure state decomposition (see Eq. 4). Generalizing the results of ref. 16, it is possible to show that the two calculations provide identical results. However, the first one is technically simpler and under much better control. As mentioned above we obtain, for all k ≥ 4 a value of α_d(k) larger than the one quoted in refs. 6 and 11.

Further we determined the distribution of cluster sizes w_n, thus unveiling a third “condensation” phase transition at α_c(k) ≥ α_d(k) (strict inequality holds for k ≥ 4 in SAT and q ≥ 4 in coloring, see below). For α < α_c(k) the weights w_n concentrate on a logarithmic scale [namely, −log w_n is θ(N) with θ(N^1/2) fluctuations]. Roughly speaking, the measure is evenly split among an exponential number of clusters.

For α > α_c(k) [and < α_s(k)] the measure is carried by a subexponential number of clusters. More precisely, the ordered sequence {w_n} converges to a well known Poisson-Dirichlet process {w_n}, first recognized in the spin glass context by Ruelle (26). This is defined by w_n = x_n/Σx_n, where x_n > 0 are the points of a Poisson process with rate x^−1−m(α) and m(α) (0, 1). This picture is known in spin glass theory as one-step replica symmetry breaking (1RSB) and has been proven in ref. 27 for some special models. The Parisi 1RSB parameter m(α) is monotonically decreasing from 1 to 0 when α increases from α_c(k) to α_s(k) (see Fig. 3).

Fig. 3. The Parisi 1RSB parameter m(α) as a function of the constraint density α. In the Inset, the complexity Σ(s) as a function of the cluster entropy for α = αs(k) − 0.1 [the slope at Σ(s) = 0 is − (more ...)

Remarkably, the condensation phase transition is also linked to an appropriate notion of correlation decay. If i(1), …, i(n) [N] are uniformly random variable indices, then, for α < α_c(k) and any fixed n:

Dynamic Phase Transition and Gibbs Measure Extremality. A rigorous calculation of α_d(k) along any of the two definitions provided above (compare Eqs. 3 and 4) remains an open problem. Each of the two approaches has, however, an heuristic implementation that we shall now describe. It can be proved that the two calculations yield equal results as further discussed in the last section.

The approach based on the extremality condition in Eq. 3 relies on an easy-to-state assumption and typically provides a more precise estimate. We begin by observing that, because of the Markov structure of μ, it is sufficient for Eq. 3 to hold that the same condition is verified by the correlation between x_i and the set of variables at distance exactly from i, that we shall keep denoting as x. The idea is then to consider a large yet finite neighborhood of i. Given ≥ , the factor graph neighborhood of radius around i converges in distribution to the radius- neighborhood of the root in a well defined random tree factor graph T.

For coloring of random regular graphs, the correct limiting tree model T is coloring on the infinite l-regular tree. For random k-SAT, T is defined by the following construction. Start from the root variable node and connect it to l new function nodes (clauses), l being a Poisson random variable of mean kα. Connect each of these function nodes with k − 1 new variables and repeat. The resulting tree is infinite with nonvanishing probability if α > 1/k(k− 1). Associate a formula to this graph in the usual way, with each variable occurrence being negated independently with probability 1/2.

The basic assumption within the first approach is that the extremality condition in Eq. 3 can be checked on the correlation between the root and generation- variables in the tree model. On the tree, μ is defined to be a translation invariant Gibbs measure (17) associated to the infinite factor graphj T (which provides a specification). The correlation between the root and generation- variables can be computed through a recursive procedure (defining a sequence of distributions , see Eq. 15 below). The recursion can be efficiently implemented numerically yielding the values presented in Table 1 for k (resp. q) = 4, 5, 6. For large k (resp. q) one can formally expand the equations on P and obtain:

The second approach to the determination of α_d(k) is based on the “cavity method” (6, 25). It begins by assuming a decomposition in pure states of the form 4 with two crucial properties: (i) if we denote by W_n the size of the nth cluster (and hence w_n = W_n/Σ W_n), then the number of clusters of size W_n = e^Ns grows approximately as e^NΣ(s); (ii) for each single-cluster measure μ_n, a correlation decay condition of the form 3 holds.

The approach aims at determining the rate function Σ(s), complexity: the result is expressed in terms of the solution of a distributional fixed point equation. For the sake of simplicity we describe here the simplest possible scenariok resulting from such a calculation (compare Fig. 4). For α < α_d,−∞(k) the cavity fixed point equation does not admit any solution: no clusters are present. At α_d,−∞(k) a solution appears, eventually yielding, for α > α_d,+ a non-negative complexity Σ(s) for some values of s ₊. The maximum and minimum such values will be denoted by s_max and s_min. At a strictly larger value α_d,0(k), Σ(s) develops a stationary point (local maximum). It turns out that α_d,0(k) coincides with the threshold computed in refs. 6, 11, and 14. In particular, α_d,0(4) ≈ 8.297, α_d,0(5) ≈ 16.12, α_d,0(6) ≈ 30.50 and l_d,0(4) = 9, l_d,0(5) = 13, l_d,0(6) = 17. For large k (resp. q), α_d,0(k) admits the same expansion as in Eqs. 6 and 7 with γ_d,0 = 1 − log 2. However, up to the larger value α_d(k), the appearance of clusters is irrelevant from the point of view of μ. In fact, within the cavity method it can be shown that e^N[s+Σ(s)] remains exponentially smaller than the total number of solutions Z: most of the solutions are in a single cluster. The value α_d(k) is determined by the appearance of a point s with Σ′(s) = −1 on the complexity curve. Correspondingly, one has Z ≈ e^N[Σ(s)+s]: most of the solutions are comprised in clusters of size about e^Ns. The entropy per variable = lim_N→∞ N⁻¹ log Z remains analytic at α_d(k).

Fig. 4. The complexity function [the number of clusters with entropy density s is eNΣ(s)] for the six-colorings of l-regular graphs with l {17, 18, 19, 20}. Circles indicate the dominating states with entropy s; the dashed lines have (more ...)

Condensation Phase Transition. As α increases above α_d, Σ(s) decreases: clusters of highly correlated solutions may no longer satisfy the newly added constraints. In Fig. 5Inset, we show the α dependency of Σ(s) for 4-SAT. In the large k limit, with α = ρ2^k we get Σ(s) = log 2 − ρ − log 2 e^−kρ + O(2^−k), and s = log 2e^−kρ + O(2^−k).

Fig. 5. Correlation function (3) between the root and generation variables in a random k-SAT tree formula. Here k = 4 and (from bottom to top) α = 9.30, 9.33, 9.35, 9.40 [recall that αd(4) ≈ 9.38]. In the Inset, the complexity (more ...)

The condensation point α_c(k) is the value of α such that Σ(s) vanishes: above α_c(k), most of the measure is contained in a subexponential number of large clustersl. Our estimates for α_c(k) are presented in Table 1 [see also Fig. 4 for Σ(s) in the six-coloring] while in the large-k limit we obtain α_c(k) = 2^k log 2 − 3/2 log 2 + O(2^−k) [recall that the SAT-UNSAT transition is at α_s(k) = 2^k log 2 − (1 + log 2)/2 + O(2^−k)] and l_c(q) = 2q log q − log q − 2 log 2 + o(1) [with the COL-UNCOL transition at l_s(q) = 2q log q − log q − 1 + o(1)]. Technically, the size of dominating clusters is found by maximizing Σ(s) + s over the s interval on which Σ(s) ≥ 0. For α [α_c(k), α_s(k)], the maximum is reached at s_max, with Σ(s_max) = 0 yielding = s_max. It turns out that the solutions are comprised within a finite number of clusters, with entropy e^Nsmax+Δ, where Δ = θ(1). The shifts Δ are asymptotically distributed according to a Poisson point process of rate e^−m(α)Δ with m(α) = −Σ′(s_max). This leads to the Poisson Dirichlet distribution of weights discussed above. Finally, the entropy per variable is nonanalytic at α_c(k).

Let us conclude by stressing two points. First, we avoided the 3-SAT and three-coloring cases. These cases [as well as the three-coloring on Erdös-Rényi graphs (25)] are particular in that the dynamic transition point α_d is determined by a local instability [a Kesten-Stigum (28, 29) condition, see also ref. 21], yielding α_d(3) ≈ 3.86 and l_d(3) = 6 (the case l = 5, q = 3 being marginal). Related to this is the fact that α_c = α_d: throughout the clustered phase, the measure is dominated by a few large clusters [technically, Σ(s) < 0 for all α > α_d]. Second, we did not check the local stability of the 1RSB calculation. By analogy with ref. 30, we expect that an instability can modify the curve Σ(s) but not the values of α_d and α_c.

Algorithmic Implications. Two message passing algorithms were studied extensively on random k-SAT: BP and SP (mixed strategies were also considered in refs. 19 and 20). A BP message ν_u→v(x) between nodes u and v on the factor graph is usually interpreted as the marginal distribution of x_u (or x_v) in a modified graphical model. An SP message is instead a distribution over such marginals P_u→v(ν). The empirical superiority of SP is usually attributed to the existence of clusters (6): the distribution P_u→v(ν) is a survey of the marginal distribution of x_u over the clusters. As a consequence, according to the standard wisdom, SP should outperform BP for α > α_d(k).

This picture, however, has several problems. Let us list two of them. First, it seems that essentially local algorithms (such as message passing ones) should be sensitive only to correlations among finite subsets of the variablesm, and these remain bounded up to the condensation transition. Recall in fact that the extremality condition in Eq. 3 involves a number of variables unbounded in N, while the weaker in Eq. 5 is satisfied up to α_c(k).

Second, it would be meaningful to weight uniformly the solutions when computing the surveys P_u→v(ν). In the cavity method jargon, this corresponds to using a 1RSB Parisi parameter r = 1 instead of r = 0 as is done in ref. 6. It is a simple algebraic fact of the cavity formalism that for r = 1 the means of the SP surveys satisfy the BP equations. Since the means are the most important statistics used by SP to find a solution, BP should perform roughly as SP. Both arguments suggest that BP should perform well up to the condensation point α_c(k). We tested this conclusion on 4-SAT at α = 9.5 [α_d(4), α_c(4)], through the following numerical experiment (compare Fig. 6). (i) Run BP for t_max iterations. (ii) Compute the BP estimates ν_i(x) for the single-bit marginals and choose the one with largest bias. (iii) Fix x_i = 0 or 1 with probabilities ν_i(0), ν_i(1). (iv) Reduce the formula accordingly (i.e., eliminate the constraints satisfied by the assignment of x_i and reduce the ones violated). This cycle is repeated until a solution is found or a contradiction is encountered. If the marginals ν_i were correct, this procedure would provide a satisfying assignment sampled uniformly from μ. In fact, we found a solution with finite probability (≈0.4), despite the fact that α > α_d(4). The experiment was repeated at α = 9 with a similar fraction of successes.

Fig. 6. Performance of BP heuristics on random 4-SAT formulae. The residual entropy per spin N−1 log Z (here we estimate it within Bethe approximation) as a function of the fraction of fixed variables. tmax = 20 in these experiments.

Above the condensation transition, correlations become too strong and the BP fixed point no longer describes the measure μ. Indeed the same algorithm proved unsuccessful at α = 9.7 [α_c(4), α_s(4)]. As mentioned above, SP can be regarded as an inference algorithm in a modified graphical model that weights preferentially small clusters. More precisely, it selects clusters of size e^N with maximizing the complexity Σ(s). With respect to the new measure, the weak correlation condition in Eq. 5 still holds and allows one to perform inference by message passing.

Within the cavity formalism, the optimal choice would be to take r ≈ m(α) [0, 1). Any parameter corresponding to a non-negative complexity r [0, m(α)] should, however, give good results. SP corresponds to the choice r = 0 that has some definite computational advantages, since messages have a compact representation in this case (they are real numbers).

This section provides some technical elements of our computation. The reader not familiar with this topic should consult refs. 6, 11, 25, and 32 for a more extensive introduction. The expert reader will find a new derivation and some hints of how we overcame technical difficulties.

On a tree factor graph, the marginals of μ, Eq. 1 can be computed recursively. The edge of the factor graph from variable node i to constraint node a (respectively from a to i) carries “message” _i→a (_a→i), a probability measure on defined as the marginal of x_i in the modified graphical model obtained by deleting constraint node a (resp. all constraint nodes around i apart from a). The messages are determined by the equations:

_a→i

_a→i

The log number of solutions, log Z, can be expressed as a sum of contributions that are local functions of the messages that solve Eqs. 8 and 9:

To move from trees to loopy graphs, we first consider an intermediate step in which the factor graph is still a tree but a subset of the variables, x_B = {x_j : j B} is fixed. We are therefore replacing the measure μ (compare Eq. 1), with the conditional one μ( · |x_B). In physics terms, the variables in x_B specify a boundary condition.

Notice that the measure μ( · |x_B) still factorizes according to (a subgraph of) the original factor graph. As a consequence, the conditional marginals μ(x_i|x_B) can be computed along the same lines as above. The messages η_i→a^x_B and ν_a→i^x_B obey Eq. 10, with an appropriate boundary condition for messages from B. Also, the number of solutions that take values x_B on j B [call it Z(x_B)] can be computed by using Eq. 11.

Next, we want to consider the boundary variables themselves as random variables. More precisely, given r , we let the boundary to be x_B with probability

The log-partition function Φ(r) = log (r) admits an expression that is analogous to Eq. 11,

As mentioned, the above derivation holds for tree factor graphs. Nevertheless, the local recursion equations 10 and 13 can be used as an heuristics on loopy factor graphs as well. Further, although we justified Eq. 13 through the introduction of a random boundary condition x_B, we can take B = and still look for nondegenerate solutions of such equations.

Starting from an arbitrary initialization of the messages, the recursions are iterated until an approximate fixed point is reached. After convergence, the distributions P_i→a, Q_a→i can be used to evaluate the potential Φ(r) (compare Eq. 14). From this we compute the complexity function Σ(r) Φ(r) − rΦ′(r) that gives access to the decomposition of μ in pure states. More precisely, Σ(r) is the exponential growth rate of the number of states with internal entropy s = Φ′(r). This is how curves such as in Fig. 4 are traced.

In practice, it can be convenient to consider the distributions of messages P_i→a, Q_a→i with respect to the graph realization. This approach is sometimes referred to as density evolution in coding theory. If one considers a uniformly random directed edge i → a (or a → i) in a rCSP instance, the corresponding message will be a random variable. After t parallel updates according to Eq. 13, the message distribution converges (in the N → ∞ limit) to a well defined law _t (for variable to constraint messages) or _t (for constraint to variable). As t → ∞, these converge to a fixed point , that satisfies the distributional equivalent of Eq. 13.

To be definite, let us consider the case of graph coloring. Since the compatibility functions are pairwise in this case (i.e., k = 2 in Eq. 1), the constraint-to-variable messages can be eliminated and Eq. 13 takes the form:

In general, finding messages distributions , that satisfy the distributional version of Eq. 13 is an extremely challenging task, even numerically. We adopted the population dynamics method (32), which represents the distributions by samples (this is closely related to particle filters in statistics). For instance, one represents by a sample of Ps, each encoded as a list of ηs. Since computer memory drastically limits the samples size, and thus the precision of the results, we worked in two directions: (i) we analytically solved the distributional equations for large k (in the case of k-SAT) or q (q-Col); and (ii) we identified and exploited simplifications arising for special values of r.

Let us briefly discuss the second point. Simplifications emerge for r = 0 and r = 1. The first case corresponds to SP. Refs. 6 and 11 showed how to compute efficiently Σ(r = 0) through population dynamics. Building on this, we could show that the clusters internal entropy s(r = 0) can be computed at a small supplementary cost.

The value r = 1 corresponds instead to the tree reconstruction problem (33): In this case (x_B) (compare Eq. 12), coincides with the marginal of μ. Averaging Eq. 13 (and the analogous one for Q_a→i) one obtains the BP equations (8 and 9), e.g., ʃ dP_i→a(η) η = _i→a. These remarks can be used to show that the constrained averages:

We thank J. Kurchan and M. Mézard for stimulating discussions. This work has been partially supported by the European Commission under Contracts EVERGROW and STIPCO.