Research Article

Analytic and Algorithmic Solution of Random Satisfiability Problems

See allHide authors and affiliations

Science  02 Aug 2002:
Vol. 297, Issue 5582, pp. 812-815
DOI: 10.1126/science.1073287


We study the satisfiability of random Boolean expressions built from many clauses with K variables per clause (K-satisfiability). Expressions with a ratio α of clauses to variables less than a threshold αc are almost always satisfiable, whereas those with a ratio above this threshold are almost always unsatisfiable. We show the existence of an intermediate phase below αc, where the proliferation of metastable states is responsible for the onset of complexity in search algorithms. We introduce a class of optimization algorithms that can deal with these metastable states; one such algorithm has been tested successfully on the largest existing benchmark of K-satisfiability.

The K-satisfiability problem (Ksat) asks whether one can satisfy simultaneously a set of Mconstraints between N Boolean variables, where each constraint is a clause built as the logical OR involving Kvariables (or their negations). Ksat is at the core of combinatorial optimization theory (1) and often serves as a benchmark for search algorithms in artificial intelligence and computer science. An efficient algorithm for solving Ksat for K ≥ 3 would immediately lead to other algorithms for efficiently solving thousands of different hard combinatorial problems. The class of combinatorial problems sharing such a crucial feature is called NP-complete (2), and it is a basic conjecture of modern computer science that no such efficient algorithm exists. Algorithms that are used to solve real-world NP-complete problems display a huge variability of running times, ranging from linear to exponential. A theory for the typical-case behavior of algorithms, on classes of random instances chosen from a given probability distribution, is therefore the natural complement to the worst-case analysis (3–5). Whereas 1sat and 2sat problems are solved efficiently by polynomial time algorithms (6), K > 2 randomly generated Boolean formulae may become extraordinarily difficult to solve: It has been observed numerically (7, 8) that computationally hard random instances are generated when the problems are critically constrained [i.e., close to the satisfiable-unsatisfiable (SAT-UNSAT) phase boundary]. The study of critical instances represents a theoretical challenge toward a greater understanding of the onset of computational complexity and the analysis of algorithms. Moreover, such hard instances are a popular test bed for the performance of search algorithms (9).

The random Ksat problem has close similarities with models of complex materials such as spin glasses, which are fundamental systems in the statistical physics of disordered systems (10). Spin glasses deal with binary variables (“spins”) interacting with random exchange couplings. Each pair of interacting spins can be seen as a constraint, and finding the state of minimal energy in a spin glass amounts to minimizing the number of violated constraints. Although the constraints in spin glasses and Ksat differ with respect to their precise form, in both cases the difficulty comes from the possible existence of “frustration” (11), which makes it difficult to find the global optimal state by a purely local optimization procedure. Links between combinatorial optimization and statistical physics have been known for years (10, 12, 13). Techniques from statistical physics are particularly useful when the size of the instance is large.

Two main categories of questions can be addressed. One type is algorithmic (e.g., finding an algorithm that decides whether an instance is SAT or UNSAT, or that tries to minimize the number of violated constraints). Another one is more theoretical and deals with random instances for which one wants to predict the typical behavior (e.g., phase transitions and structure of the solution space).

We address the two types of questions in the 3sat problem. When the numbers of variables N and of clauses M both increase at a fixed value of α = M/N, random Ksat problems become generically SAT at small α and generically UNSAT at large α. The existence of a SAT-UNSAT phase transition in the infinite N limit has been established rigorously for any K (14), but the critical value αc (that separates the two phases) has been found only in the (polynomial) K = 2 problem where αc = 1 (15–17). For the NP-complete caseK ≥ 3, much less is known. The present best numerical estimate for αc at K = 3 is 4.26 (18), and the rigorous bounds are 3.42 < αc < 4.506 (19, 20); previous statistical mechanics analysis, using the replica method, has found αc(3) ∼ 4.48 (21) and αc(3) ∼ 4.396 (22) in the framework of variational approximations. The SAT-UNSAT decision problem is also known experimentally to be algorithmically harder to solve in the neighborhood of αc, depending on the characteristics of the SAT-UNSAT phase transition. Indeed, 2sat and 3sat are different in this respect (23).

Setting out the statistical physics problem.

The Ksat problem deals with N Boolean variablesxi , i ∈ {1, ..., N}. Each clause a∈ {1, ..., M} involvesK variables {xi 1 (a), ...,x iK (a)}. Each such variable can be negated or not, and the clause is built as the OR function of the K resulting variables. In physical terms, the variable xi can be represented by a “spin” si = ±1 through the one-to-one mapping si = −1 if xi is false and si = +1 if xi is true. For each variablex ir (a)appearing in clause a, one introduces a “coupling”J a r = −1 if the variable appears negated in the clause; otherwise the coupling isJ a r = 1. The sets of indicesi 1(a), ...,iK (a) and of “couplings”Ja = {J a 1,...,J a K} define an instance of the problem under study. Given a spin configuration, the “energy”ɛ Ja(si 1 (a),...,s ik (a)) of clause a is equal to 0 if the clause is satisfied, or equal to 1 if it is violated (24). The total energyE equals the number of violated clauses.

In statistical physics, one assigns to each of the 2N spin configurations a Boltzmann probability exp(−βE)/Z, where β is an auxiliary parameter playing the role of the inverse of temperature, andZ is a normalization term; here we are interested in the β → ∞ “zero-temperature” limit, where Boltzmann's law selects optimal states.

The spin glass approach.

We first study the large N limit of the random 3sat problem, where the indices in each clause are chosen randomly, as well as the sign of each coupling, with uniform distributions. Our approach to these problems uses a general strategy initiated years ago in spin glass theory (10). The first concept we need to introduce is that of a state. Roughly speaking, states correspond to connected regions of configurations, such that one must cross energy barriers that diverge when N → ∞ to go from one state to another. The archetype of such a situation is the ferromagnetic transition where the spins collectively polarize, either toward an “up” state or toward a “down” state. In frustrated systems such as satisfiability problems, many states can exist: The number of states with energyE behaves as exp[N Σ(e)], whereeE/N; the function Σ(e), called the complexity, is a crucial concept in studies of structural glasses. The ground-state energy densitye can be found by the condition Σ(e) = 0. Here we choose a restricted zero-temperature definition that applies to random Ksat: A state is simply a cluster of configurations of equal energy related by single spin-flip moves, such that the energy cannot be decreased by any sequence of single spin flips (25). Generalizing the approach of (26), one can develop a whole “zero-temperature thermodynamics” of the states by introducing a “free energy” function Φ(y) defined fromEmbedded Image(1)The reweighting y is a Lagrange parameter (similar to an inverse temperature) that allows the energy of the states to be fixed. Larger reweighting selects states of lower energies until one reaches y = y*, corresponding to the lowest energy states (y* = ∞ in the SAT region).

The cavity method: Message-passing procedures.

To compute Φ(y), we use the zero-temperature cavity method (27, 28), in which the basic ingredients are the cavity fields and the cavity biases, which are defined in each state. The cavity field hi ameasures the tendency of spin i to be up, when one of the clauses, a, to which i belongs, has been disconnected (Fig. 1). It is equal to the sum of cavity biasesub i, sent to sitei from all the other clauses b to which it belongs. In computer science terminology, cavity fields are messages sent from a variable node to a function node, whereas cavity biases are messages sent from a function node to a variable node (29). The cavity biases are determined by a local optimization procedure. Consider one clause a, involving K variabless 1, ...,sK , and a penalty functionɛ J(s 1, ..., sK ). The optimization on the variables s 2, ..., sK , expressed as Embedded Image Embedded Image Embedded Image(2)defines the mean energy shiftaJ and the cavity biasua →1 =uJ (h 2, ..., hK ) propagated from this clause to the variable s 1 (30).

Figure 1

In the random 3sat problem, the graph of clauses is locally isomorphic to a tree. Variables (spins) are depicted by a circle, and clauses by a square. The cavity biasua →1 sent from the red clausea to the variable s 1 summarizes the effects of optimizing clause a on s 2,s 3, taking into account all the blue + green (top) part of the graph, when the yellow (bottom) part has been taken away. In the absence of the red clause, the cavity fieldh 2→a sums up all cavity biasesw 1, ...,wq arriving onto s 2 from the blue clauses, and the cavity fieldh 3→a sums up all cavity biasesv 1, ...,vq arriving ontos 3 from the green clauses.

The advantage of cavity biases and cavity fields in large (N≫ 1) random Ksat and spin glass problems is the special structure of the interaction graph: It is locally tree-like, and the connectivity fluctuates from site to site with a Poisson distribution of meanKα (see Fig. 1). On a more global scale, these random graphs have loops with a typical length growing as log(N). As the cavity fields h 2, ..., hK are defined in the absence of the clause, they correspond to faraway variables [with a distance of order log(N)]. The “clustering property,” valid inside each state, implies that their correlations go to zero at large N (on a real tree they would be fully uncorrelated). The topology of the graph implies that the cavity equations are exact on finite subgraphs.

To determine the statistical properties of the set of cavity biases and how they change from state to state, we introduce “surveys,” which are histograms of cavity biases. For each state ω, there is one cavity bias u a→1 ωpropagated from one clause a to site 1; it can be computed from Eq. 2, where the cavity fields are those corresponding to the state ω. For a given value of the reweighting, the survey propagated to spin s 1 in Fig. 1 is defined asQ a→1 (y)(u) = CΣωδ(u u a→1 ω), whereC is a normalization constant ensuring thatQ a→1 (y)(u) is a probability distribution, and the sum over ω is restricted to the states having the energy density e selected by the reweighting y.

The survey propagation rules on the graph of Fig. 1 take the precise formEmbedded Image Embedded Image Embedded Image(3)where W =w 1 + ... +wq , V =v 1 + ... +vq′ , and C′ is a normalization constant. The exponential term in Eq. 3 takes care of the energy-level crossings induced by the propagation. Once the surveys are known, the free energy Φ(y) can be computed using the formulae of (27, 28), and the complexity can be deduced from Eq. 1. The order parameter of the theory is the collection of the surveys.

Phase diagram.

In the zero-temperature 3sat problem, one sees from the definition (2) that a given cavity biasua i takes either the values {0,1} (if the variable xi appears negated in clause a) or the values {0,−1}. The corresponding surveyQ ai (y)(u) is thus characterized by a single number, the probability thatu = 0. Using this simplification, we have been able to compute (31) the statistical distribution of surveys in random graphs in the infinite volume limit. We find two critical values of α at αd ∼ 3.921 and αc ∼ 4.256. For α < αd, the solution is of a paramagnetic type [all the surveys equal δ(u)], a generic instance is satisfiable, and the solution can be found even by a simple zero-temperature Metropolis algorithm (ZTMA) (32). For αd < α < αc, the space of configurations breaks up into many states, and there exists a nontrivial complexity (33). Some of the states have zero energy; therefore, we are still in the SAT phase. It can be argued that algorithms like ZTMA will generically get trapped into the most numerous states, which have an extensive (proportional to N) energyE th.

At α = 4.2 we find analytically E th∼ 0.0036N, and we have checked that ZTMA converges to a similar value of energy. The fact that e th= E th/N is small explains the good performance of smarter algorithms on instances involving a few thousand variables. At α > αc, the system is in its UNSAT phase, and the lowest possible energy is positive. The phase diagram is summarized in Fig. 2.

Figure 2

The phase diagram of the random 3sat problem. Plotted is e 0, the number of violated clauses per variable (red), versus the control parameter α, which is the number of clauses per variable. The SAT-UNSAT transition occurs at α = αc ∼ 4.256. The green line is e th, the threshold energy per variable, where local algorithms get trapped. The blue line is the complexity Σ of satisfiable states, equal to 1/N times the logarithm of their number.

Survey propagation algorithm.

We now consider one given instance (31), that is, one fixed large graph. We have seen experimentally that in the glassy region α > αd, the standard (y = 0) iteration of cavity biases either ceases to converge or converges to the trivial paramagnetic solution where allua i = 0. Ifi is the rth site connected to the function nodea, we introduce a surveyQ ai (y)(u) = ηaiδ(u) + (1 − ηai)δ(u+ J a r) that is characterized by the single number ηai. The survey propagation of Eq. 3 performed with random sequential updating is a message-passing procedure that defines a dynamical process in the space of the KN variables ηai. We have implemented it on large random instances in the hard part of the SAT phase, with α ∼ 4.2 to 4.25, using a sufficiently large value of y(typically y ∼ 4 to 6). The process is found to converge to a unique nontrivial solution. We expect that this survey propagation technique can be of interest in many problems of statistical inference.

The set of all surveysQ ai (y)(u) found after convergence provides a nontrivial information on the structure of the states. From all the surveys sent onto one sitei, we reconstruct through a reweighted convolution (34) the probability distribution of local fields on this site, Pi (H). This is a distribution on integers [Pi (H) = Σrδ(H r)w i r]. The total weight wi + = Σr=1 w i r of Pi (H) on positive integers gives the fraction of zero-energy states wheresi = 1; similarly, the total weightwi = Σr=−∞ −1 w i r ofPi (H) on negative integers gives the fraction of zero-energy states where si = −1. We have checked numerically, on single instances with N= 10,000, that these fractions predicted from survey propagation agree with those obtained by averaging on a few hundreds of ground states.

A decimation algorithm.

This information can be exploited to invent new types of algorithms (31) or to improve existing ones. We have worked out one such application, the survey inspired decimation (SID), which shows promising performance, but other algorithms probably could be found using the same type of information. Given an instance, we first compute all the surveys by the survey propagation algorithm with a sufficiently large value of y (e.g., y = 6). Then we deduce the distribution of local fields, and in particular their weights wi ± on positive and negative integers. We then fix the variable i with largest |wi +wi | to the valuesi =Sign(wi +wi ). Satisfied clauses are eliminated, and unsatisfied K-clauses involving i are transformed into K − 1 clauses, leading to a new instance with a reduced number of variables (and of clauses). The surveys can be propagated again on this new instance (starting from the previous ones) until convergence, and the procedure is iterated. Whenever a paramagnetic state is found (signaled by all ηai = 1) or at some intermediate steps, a rapid search process like simulated annealing at a fixed cooling rate is run.

This SID algorithm has been tested successfully on the largest (up to N = 2000) existing benchmarks (9) of random 3sat instances in the hard regime. Satisfying assignments have been found for all benchmarks. We have applied the SID to much larger instances, increasing N up to N = 105 at a fixed α = 4.2. The algorithm is very efficient: It always finds a SAT configuration, and its apparent complexity scales like N 2, although more systematic studies with higher statistics will be necessary to establish this behavior. For the very same large instances, the only existing algorithm able to find solutions, at a considerable computational cost, is a highly optimized version of the walksat algorithm (9, 35).


We have proposed an analytical method that predicts quantitatively the phase diagram of the random 3sat problem in the limit of infinite number of clauses and opens the way to other types of algorithms. The existence of an intermediate phase with many metastable states close to the SAT-UNSAT transition explains the slowing down of algorithms in this region. We would like to stress that the solution we propose is typical of a “one-step replica symmetry-breaking” solution, as it is called in spin glasses (10). All the consistency checks of the analytic results lead us to believe that this solution is exact for the 3sat problem. From the strict mathematical point of view, the phase diagram we propose should be considered as a conjecture, as for the great majority of the theoretical results in statistical physics. Our computation implies that a way to provide a fully rigorous proof of the transition behavior in random Ksat problems could be based on the study of the decomposition of the probability measure into states endowed with the clustering property (36). On the other hand, the predictions of our theory can be compared with numerical experiments, and our first such tests have confirmed its validity. On the basis of the analytical study, our algorithm looks promising in that it can solve large instances exploring a rather small number of spin configurations. It will be interesting to explore its application to other optimization problems.

  • * To whom correspondence should be addressed. E-mail: zecchina{at}


View Abstract

Navigate This Article