## Abstract

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 *M*constraints between *N* Boolean variables, where each constraint is a clause built as the logical OR involving *K*variables (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 case*K* ≥ 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 variables*x _{i}
*,

*i*∈ {1, ...,

*N*}. Each clause

*a*∈ {1, ...,

*M*} involves

*K*variables {

*x*

_{i}_{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

*x*can be represented by a “spin”

_{i}*s*= ±1 through the one-to-one mapping

_{i}*s*= −1 if

_{i}*x*is false and

_{i}*s*= +1 if

_{i}*x*is true. For each variable

_{i}*x*

_{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 is

*J*

_{a}

^{r}= 1. The sets of indices

*i*

_{1}(

*a*), ...,

*i*(

_{K}*a*) and of “couplings”

**J**= {

_{a}*J*

_{a}

^{1},...,

*J*

_{a}

^{K}} define an instance of the problem under study. Given a spin configuration, the “energy”

*ɛ*

_{Ja}(

*s*

_{i}_{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 energy

*E*equals the number of violated clauses.

In statistical physics, one assigns to each of the 2^{N} spin configurations a Boltzmann probability exp(−β*E*)/*Z*, where β is an auxiliary parameter playing the role of the inverse of temperature, and*Z* 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 energy*E* behaves as exp[*N* Σ(*e*)], where*e* ≡ *E*/*N*; the function Σ(*e*), called the complexity, is a crucial concept in studies of structural glasses. The ground-state energy density*e* 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 from(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 *h _{i}
*

_{→a}measures 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 biases

*u*

_{b}_{→i}, sent to site

*i*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*variables

*s*

_{1}, ...,

*s*, and a penalty function

_{K}*ɛ*

_{J}(

*s*

_{1}, ...,

*s*). The optimization on the variables

_{K}*s*

_{2}, ...,

*s*, expressed as (2)defines the mean energy shift

_{K}*a*and the cavity bias

_{J}*u*

_{a}_{→1}=

*u*(

_{J}*h*

_{2}, ...,

*h*) propagated from this clause to the variable

_{K}*s*

_{1}(30).

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 mean*K*α (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}, ..., *h _{K}
* 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 as*Q*
_{a→1}
^{(y)}(*u*) = *C*Σ_{ω}δ(*u *−*u*
_{a→1}
^{ω}), where*C* is a normalization constant ensuring that*Q*
_{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 form
(3)where *W* =*w*
_{1} + ... +*w _{q}
*,

*V*=

*v*

_{1}+ ... +

*v*, and

_{q′}*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 bias*u _{a}
*

_{→i}takes either the values {0,1} (if the variable

*x*appears negated in clause

_{i}*a*) or the values {0,−1}. The corresponding survey

*Q*

_{a→i}

^{(y)}(

*u*) is thus characterized by a single number, the probability that

*u*= 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*) energy

*E*

_{th}.

At α = 4.2 we find analytically *E*
_{th}∼ 0.0036*N*, 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.

#### 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 all*u _{a}
*

_{→i}= 0. If

*i*is the

*r*th site connected to the function node

*a*, we introduce a survey

*Q*

_{a→i}

^{(y)}(

*u*) = η

_{a→i}δ(

*u*) + (1 − η

_{a→i})δ(

*u*+

*J*

_{a}

^{r}) that is characterized by the single number η

_{a→i}. 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 η

_{a→i}. 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 surveys*Q*
_{a→i}
^{(y)}(*u*) found after convergence provides a nontrivial information on the structure of the states. From all the surveys sent onto one site*i*, we reconstruct through a reweighted convolution (34) the probability distribution of local fields on this site, *P _{i}
*(

*H*). This is a distribution on integers [

*P*(

_{i}*H*) = Σ

_{r}δ(

*H*−

*r*)

*w*

_{i}

^{r}]. The total weight

*w*

_{i}^{+}= Σ

_{r=1}

^{∞}

*w*

_{i}

^{r}of

*P*(

_{i}*H*) on positive integers gives the fraction of zero-energy states where

*s*= 1; similarly, the total weight

_{i}*w*

_{i}^{−}= Σ

_{r=−∞}

^{−1}

*w*

_{i}

^{r}of

*P*(

_{i}*H*) on negative integers gives the fraction of zero-energy states where

*s*= −1. We have checked numerically, on single instances with

_{i}*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 *w _{i}
*

^{±}on positive and negative integers. We then fix the variable

*i*with largest |

*w*

_{i}^{+}−

*w*

_{i}^{−}| to the value

*s*=

_{i}*Sign*(

*w*

_{i}^{+}−

*w*

_{i}^{−}). 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 η

_{a→i}= 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* = 10^{5} 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).

#### Conclusions.

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}ictp.trieste.it