Solving Constraint Satisfaction Problems through Belief Propagation-guided decimation

I Introduction

An instance of a constraint satisfaction problem consists of $n$ variables $\underline{x}=(x_{1},\dots,x_{n})$ and $m$ constraints among them. Solving such an instance amounts to finding an assignment of the variables that satisfies all the constraints, or proving that no such assignment exists. A remarkable example in this class is provided by $k$-satisfiability, where variables are binary, $x_{i}\in\{{\tt 0},{\tt 1}\}$, and each constraint requires $k$ of the variables to be different from a specific $k$-uple. Explicitly, the $a$-th constraint (clause), $a\in[m]\equiv\{1,\dots,m\}$ is specified by $k$ variable indexes $i_{1}(a),\dots,i_{k}(a)\in[n]$, and $k$ bits $z_{1}(a),\dots,z_{k}(a)\in\{{\tt 0},{\tt 1}\}$. Clause $a$ is satisfied by assignment $\underline{x}$ if and only if $(x_{i_{1}(a)},\dots,x_{i_{k}(a)})\neq(z_{1}(a),\dots,z_{k}(a))$.

A constraint satisfaction problem admits a natural factor graph representation, cf. Fig. 1. Given an instance, each variable can be associated to a variable node, and each constraint to a factor node. Edges connect factor node $a\in F\equiv[m]$ to those variable nodes $i\in V\equiv[n]$ such that the $a$-th constraint depends in a non-trivial way on variable $x_{i}$. For instance, in the case of $k$-satisfiability, clause $a$ is connected to variables $i_{1}(a),\dots,i_{k}(a)$. If the resulting graph is sparse, fast message passing algorithms can be defined on it.

Although constraint satisfaction problems are generally NP-hard, a large effort has been devoted to the development of efficient heuristics. Recently, considerable progress has been achieved in building efficient ‘incomplete solvers’ . These are algorithms that look for a solution but, if they do not find one, cannot prove that the problem is unsolvable. A particularly interesting class is provided by message passing-guided decimation procedures. These consist in iterating the following steps:

Use the result to choose a variable index $i\in V$, and a value $x_{i}^{*}$ for the corresponding variable.

Replace the constraint satisfaction problem with the one obtained by fixing $x_{i}$ to $x_{i}^{*}$.

The iteration may stop for two reasons. In the first case a contradiction is produced: the same variable $x_{i}$ appears in two constraints whose other arguments have already been fixed, and that are satisfied by distinct values of $x_{i}$. If this does not happen, the iteration stops only when all the variables are fixed and a solution is found. Notice that earlier algorithms, such as unit clause propagation (UCP) did not used message passing in step 2, and were not nearly as effective.

Random constraint satisfaction problems are a useful testing ground for new heuristics. For instance, random $k$-satisfiability is the distribution over $k$-SAT formulae defined by picking a formula uniformly at random among all the ones including $m$ clauses over $n$ variables. Decimation procedures of the type sketched above proved particularly successful in this context. In particular survey propagation-guided decimation outperformed the best previous heuristics based on stochastic local search . More recently belief propagation-guided decimation was shown empirically to have good performances as well .

Unfortunately, so far there exists no analysis of message-passing guided decimation. Our understanding almost entirely relies on simulations, even for random instances. Consequently the comparison among different heuristics, as well as the underpinnings of their effectiveness are somewhat unclear. In this paper we define a simple class of randomized message passing-guided decimation algorithms, and present a technique for analyzing them on random instances. The technique is based on the identification of a process on infinite trees that describes the evolution of the decimation algorithm. The tree process is then analyzed through an appropriate generalization of density evolution . Our approach is close in spirit to the one of . While it applies to a large class of random constraint satisfaction problems (including, e.g. coloring of random graphs), for the sake concreteness, we will focus on random $k$-SAT.

We expect the tree process to describe exactly the algorithm behavior in the limit of large instances, $n\to\infty$. While we could not prove this point, numerical simulations convincingly support this conjecture. Further, non-rigorous predictions based on tree calculations have been repeatedly successful in the analysis of random $k$-satisfiability. This approach goes under the name of ‘cavity method’ in statistical mechanics .

The paper is organized as follows. Section II contains some necessary background and notation on random $k$-SAT as well as a synthetic discussion of related work. In Section III we define the decimation procedure that we are going to analyze. We further provide the basic intuition behind the definition of the tree model. The latter is analyzed in Section IV, and the predictions thus derived are compared with numerical simulations in Section V. Finally, some conclusions and suggestions for future work are presented in Section VI. Proofs of several auxiliary lemmas are omitted from this extended abstract and deferred to technical appendices.

II Random kk-SAT and message passing: Background and related work

As mentioned above, random $k$-SAT refers to the uniform distribution over $k$-SAT instances with $m$ constraints over $n$ variables. More explicitly, each constraint is drawn uniformly at random among the $2^{k}\binom{n}{k}$ possible ones. We are interested here in the limit $n,m\to\infty$ with $m/n=\alpha$ fixed.

We let $\alpha_{\rm s}(k)$ be the largest value of $\alpha$ such that random $k$-SAT instances admit with high probability a solution. It is known that $\alpha_{\rm s}(k)=2^{k}\log 2-O(k)$. A sharp conjecture on the value of $\alpha_{\rm s}(k)$ has been put forward in on the basis of statistical physics calculations, implying $\alpha_{\rm s}(k)\approx 4.267$, $9.93$, $21.12$ for (respectively) $k=3,4,5$ and $\alpha_{\rm s}(k)=2^{k}\log 2-\frac{1}{2}(1+\log 2)+O(2^{-k})$ for large $k$ .

Simple heuristics have been analyzed thoroughly and proved to find a solution with probability bounded away from $0$ if $\alpha\leq{\rm const}\,2^{k}/k$. Here the proportionality constant depends on the specific heuristic.

To the best of our knowledge, the first application of message passing algorithms to $k$-satisfiability is reported in . In this early study BP was mostly applied in a one-shot fashion (as in iterative decoding of sparse graph codes ), without decimation. By this we mean that belief propagation is run, and resulting marginal probabilities are used to guess the values of all variables at once. However the probability of success of the one-shot algorithm is exponentially small: there are $\Theta(n)$ isolated constraints, whose variables have non-trivial marginal probabilities, each of them is hence violated with finite probability in the one-shot assignment.

Statistical mechanics methods allowed to derive a very precise picture of the solution set . This inspired a new message passing algorithm dubbed survey propagation . In conjunction with decimation, this algorithm allowed to solve random instances of unprecedentedly large sizes, in difficult regimes of $\alpha$ and $k$.

In it was proved that for $\alpha\leq(2\log k)/k\,[1+o(1)]$, BP computes good approximations of the marginals of $\mu$, irrespective of its initialization. It is clear from empirical studies that the ‘worst case’ argument used in this estimate (and in other papers on belief propagation ) is far too pessimistic.

In Ref. a simple message passing algorithm, warning propagation (see below), was analyzed for a modified (‘planted’) ensemble of random formulae. The algorithm was proved to converge and find solutions for large enough density $\alpha$ (see also ). Both the ensemble and the algorithm are quite different from the ones treated in this paper.

Further, the definition and analysis of a ‘Maxwell decoder’ in , is closely related to the approach in this paper. Let us recall that the Maxwell decoder was a (mostly conceptual) algorithm for implementing maximum likelihood decoding of LDPC codes over the erasure channel. The treatment in applies almost verbatim to a simple constraint satisfaction problem known as XORSAT. The generalization in the present paper is analogous to the one from the erasure to a general binary memoryless symmetric channel.

Finally, let us mention that BP decimation can be an interesting option in engineering applications, as demonstrated empirically in the case of lossy source coding .

III A simple decimation procedure

Let us recall the definition of BP for our specific setup ( are general references). BP is a message passing algorithm: at each iteration messages are sent from variable nodes to neighboring clause nodes and vice versa. To describe the message update equations, we need some more notation. As in the case of factor nodes, we shall call ${\partial i}$ the set of factors that depends on the variable $x_{i}$. If $i\in{\partial a}$, say $i=i_{l}(a)$, we denote $z(i,a)=z_{l}(a)$ the value of $x_{i}$ which does not satisfy the $a$-th clause. For a pair of adjacent variable ($i$) and factor ($a$) nodes (i.e. $i\in{\partial a}$), let us call ${\partial_{+}i}(a)$ (resp. ${\partial_{-}i}(a)$) the set of factor nodes adjacent to $i$, distinct from $a$, that agrees (resp. disagrees) with $a$ on the satisfying value of $x_{i}$. In formulae, ${\partial_{+}i}(a)=\{b\in{\partial i}\setminus a|z(i,b)=z(i,a)\}$ and ${\partial_{-}i}(a)=\{b\in{\partial i}|z(i,b)=1-z(i,a)\}$.

It is convenient to use log-likelihood notations for messages as is done in iterative decoding , with two caveats: $(1)$ We introduce a factor $1/2$ to be consistent with physics notation; $(2)$ The message from variable node $i$ to factor node $a$ corresponds to the log-likelihood for $x_{i}$ to satisfy/not-satisfy clause $a$ (rather than to be ${\tt 0}/{\tt 1}$).

Let $\{h^{(r)}_{i\to a}\}$, $\{u^{(r)}_{a\to i}\}$ denote the messages that are passed at time $r$ along the directed edges $i\to a$ and $a\to i$, for $i\in V$, and $a\in F$. The update equations read

with $\epsilon=0$ (this parameter is introduced for the discussion in Sec. V).

For $i\in V$, let ${\partial_{+}i}$ be the subset of clauses that are satisfied by $x_{i}={\tt 0}$, and ${\partial_{-}i}$ the subset satisfied by $x_{i}={\tt 1}$. Then the BP estimate for the marginal of $x_{i}$ under the measure $\mu(\,\cdot\,)$ is $\nu^{(r)}_{i}(x_{i})$, where

III-B Unit clause and warning propagation

During the decimation procedure a subset $U$ of the variables are fixed to specific values, collectively denoted as $\underline{x}_{U}^{*}$. This has some direct implications. By this we mean that for some other variables $x_{j}$, $j\in V\setminus U$, it follows from ‘unit clause propagation’ (UCP) that they take the same value in all of the solutions compatible with the partial assignment $\underline{x}_{U}^{*}$. We will say that these variables are directly implied by the condition $\underline{x}_{U}=\underline{x}_{U}^{*}$. Let us recall that unit clause propagation corresponds to the following deduction procedure. For each of the fixed variables $x_{i}$, and each of the clauses $a$ it belongs to, the value $x_{i}^{*}$ can either satisfy clause $a$, or not. In the first case clause $a$ can be eliminated from the factor graph. In the second a smaller clause with one less variable is implied. In both cases variable $x_{i}$ is removed. It can happen that the size of a clause gets reduced to $1$, through this procedure. In this case the only variable belonging to the clause must take a definite value in order to satisfy it. We say that such a variable is directly implied by the fixed ones. Whenever a variable is directly implied, its value can be substituted in all the clauses it belongs to, thus allowing further reductions.

The process stops for one of two reasons: $(1)$ All the fixed or directly implied variables have been pruned and no unit clause is present in the reduced formula. In this case we refer to all variables that appeared at some point in a unit clause as directly implied variables. $(2)$ Two unit clauses imply different values for the same variable. We will say that a contradiction was revealed in this case: no solution $\underline{x}$ of the formula can verify the condition $\underline{x}_{U}=\underline{x}_{U}^{*}$.

A key element in our analysis is the remark that UCP admits a message passing description. The corresponding algorithm is usually referred to as warning propagation (WP) . The WP messages (to be denoted as ${\mathfrak{h}}_{i\to a}^{(r)}$, ${\mathfrak{u}}_{a\to i}^{(r)}$) take values in $\{{\tt I},{\tt 0}\}$. The meaning of ${\mathfrak{u}}_{a\to i}^{(r)}={\tt I}$ (respectively ${\mathfrak{u}}_{a\to i}^{(r)}={\tt 0}$) is: ‘variable $x_{i}$ is (resp. is not) directly implied by clause $a$ to satisfy it.’ For variable-to-factor messages, the meaning of ${\mathfrak{h}}_{i\to a}^{(r)}={\tt I}$ (respectively ${\mathfrak{h}}_{i\to a}^{(r)}={\tt 0}$) is: ‘variable $x_{i}$ is (resp. is not) directly implied, through one of the clauses $b\in{\partial i}\setminus a$, not to satisfy clause $a$.’

We want to apply WP to the case in which a part of the variables have been fixed, namely $x_{i}=x^{*}_{i}$ for any $i\in U\subseteq V$. In this case the WP rules read

In the following we shall always assume that WP is initialized with ${\mathfrak{h}}_{i\to a}^{(0)}={\tt 0}$, ${\mathfrak{u}}_{i\to a}^{(0)}={\tt 0}$ for each edge $(ia)\in E$. It is then easy to prove that messages are monotone in the iteration number (according to the ordering ${\tt 0}<{\tt I}$). In particular the WP iteration converges in at most $O(n)$ iterations. We denote $\{{\mathfrak{u}}_{a\to i}^{(\infty)}\}$ the corresponding fixed point messages, and say that $i\in V\setminus U$ is WP-implied by the fixed variables if there exist $a\in{\partial i}$ such that ${\mathfrak{u}}_{a\to i}^{(\infty)}={\tt I}$. Then the equivalence between UCP and WP can be stated in the form below.

Assume a partial assignment $\underline{x}_{U}^{*}$ to be given for $U\subseteq V$. Then

The fixed point WP messages $\{{\mathfrak{u}}_{a\to i}^{(\infty)}\}$ do not depend on the order of the WP updates (as long as any variable is updated an a priori unlimited number of times).

$i\in V\setminus U$ is directly implied iff it is WP-implied.

UCP encounters a contradiction iff there exists $i\in V$, $a\in{\partial_{+}i}$, $b\in{\partial_{-}i}$ such that ${\mathfrak{u}}_{a\to i}^{(\infty)}={\mathfrak{u}}_{b\to i}^{(\infty)}={\tt I}$.

For the clarity of what follows let us emphasize the terminology of fixed variables (those in $U$) and of directly implied variables (not in $U$, but implied by $\underline{x}_{U}^{*}$ though UCP or WP). Finally, we will call frozen variables the union of fixed and directly implied ones, and denote the set of frozen variables by $W\subseteq V$.

III-C Decimation

The BP-guided decimation algorithm is defined by the pseudocode of Table I.

There are still a couple of elements we need to specify. First of all, how the BP equations (2), (3) are modified when a non-empty subset $U$ of the variables is fixed. One option would be to eliminate these variables from the factor graph, and reduce the clauses they belong to accordingly. A simpler approach consists in modifying Eq. (2) when $i\in U$. Explicitly, if the chosen value $x_{i}^{*}$ satisfies clause $a$, then we set $h_{i\to a}^{(r+1)}=+\infty$. If it does not, we set $h_{i\to a}^{(r+1)}=-\infty$.

Next, let us stress that, while WP is run until convergence, a not-yet defined ‘stopping criterion’ is used for BP. This will be precised in Section V. Here we just say that it includes a maximum iteration number $r_{\rm max}$, which is kept of smaller order than $O(n)$.

The algorithm complexity is therefore naively $O(n^{3}r_{\rm max})$. It requires $n$ cycles, each involving: $(1)$ at most $r_{\rm max}$ BP iterations of $O(n)$ complexity and $(2)$ at most $n$ WP iterations of complexity $O(n)$. It is easy to reduce the complexity to $O(n^{2}r_{\rm max})$ by updating WP in sequential (instead of parallel) order, as in UCP. Finally, natural choice (corresponding to the assumption that BP converges exponentially fast) is to take $r_{\rm max}=O(\log n)$, leading to $O(n^{2}\log n)$ complexity.

In practice WP converges after a small number of iterations, and the BP updates are the most expensive part of the algorithm. This could be reduced further by using the fact that fixing a single variable should produce only a small change in the messages. Ref. uses this argument for a similar algorithm to argue that $O(n\log n)$ time is enough.

III-D Intuitive picture

Analyzing the dynamics of BP-decimation seems extremely challenging. The problem is that the procedure is not ‘myopic’ , in the sense that the value chosen for variable $x_{i}$ depends on a large neighborhood of node $i$ in the factor graph. By analogy with myopic decimation algorithms one expects the existence of a critical value of the clause density $\alpha_{{\mbox{\footnotesize BPd}}}(k)$ such that the algorithm finds a solution with probability bounded away from $0$ for $\alpha<\alpha_{{\mbox{\footnotesize BPd}}}(k)$, while it is unsuccessful with high probability for $\alpha>\alpha_{{\mbox{\footnotesize BPd}}}(k)$. Notice that, if the algorithm finds a solution with positive probability, restarting it a finite number of times should A caveat: here we are blurring the distinction between probability with respect to the formula realization and the algorithm realization. yield a solution with probability arbitrarily close to $1$.

We shall argue in favor of this scenario and present an approach to analyze the algorithm evolution for $\alpha$ smaller than a spinodal point $\alpha_{{\mbox{\footnotesize spin}}}(k)$. More precisely, our analysis allows to compute the asymptotic fraction of ‘directly implied’ variables after any number of iterations. Further, the outcome of this computation provides a strong indication that $\alpha_{{\mbox{\footnotesize spin}}}(k)\leq\alpha_{{\mbox{\footnotesize BPd}}}(k)$. Both the analysis, and the conclusion that $\alpha_{{\mbox{\footnotesize spin}}}(k)\leq\alpha_{{\mbox{\footnotesize BPd}}}(k)$ are confirmed by large scale numerical simulations.

Our argument goes in two steps. In this section we show how to reduce the description of the algorithm to a sequence of ‘static’ problems. The resolution of the latter will be treated in the next section. Both parts rely on some assumptions on the asymptotic behavior of large random $k$-SAT instances, that originate in the statistical mechanics treatment of this problem . We will spell out such assumptions along the way.

As a preliminary remark, notice that the two message passing algorithms play different roles in the BP-decimation procedure of Table I. BP is used to estimate marginals of the uniform measure $\mu(\,\cdot\,)$ over solutions, cf. Eq. (1), in the first repetition of the loop. In subsequent repetitions, it is used to compute marginals of the conditional distribution, given the current assignment $\underline{x}_{U}=\underline{x}^{*}_{U}$. These marginals are in turn used to choose the values $\{x_{i}^{*}\}$ of variables to be fixed. WP is on the other hand used to check a necessary condition for the current partial assignment to be consistent. Namely it checks if it induces a contradiction on directly implied variables. In fact, it could be replaced by UCP, and, in any case, it does not influence the evolution of the partial assignment $\underline{x}^{*}_{U}$.

Let us introduce some notation: $(i(1),i(2),\dots,i(n))$ is the order in which variables are chosen at step $5$ in the algorithm, $U_{t}=\{i(1),\dots,i(t)\}$ the set of fixed variables at the beginning of the $t+1$-th repetition of the loop, and $W_{t}$ the frozen variables at that time (i.e. the union of $U_{t}$ and the variables directly implied by $\underline{x}^{*}_{U_{t}}$).

We begin the argument by considering an ‘idealized’ version of the algorithm where BP is replaced by a black box, that is able to return the exact marginal distribution of the measure conditioned on the previous choices, namely $\nu_{i}(\,\cdot\,)=\mu_{i|U}(\,\cdot\,|\underline{x}_{U}^{*})$. Let us point out two simple properties of this idealized algorithm. First, it always finds a solution if the input formula is satisfiable (this will be the case with high probability if we assume $\alpha<\alpha_{\rm s}(k)$). In fact, assume by contradiction that the algorithm fails. Then, there has been a last time $t$, such that the $k$-SAT instance has at least one solution consistent with the condition $\underline{x}_{U_{t-1}}=\underline{x}^{*}_{U_{t-1}}$, but no solution under the additional constraint $x_{i}=x_{i}^{*}$ for $i=i(t)$. This cannot happen because it would imply $\mu_{i|U_{t-1}}(x^{*}_{i}|\underline{x}^{*}_{U_{t-1}})=0$, and if this is the case, we would not have chosen $x^{*}_{i}$ in step $8$ of the algorithm.

The second consequence is that the algorithm output configuration $\underline{x}^{*}$ is a uniformly random solution. This follows from our assumption since

Therefore, the distribution of the state of the idealized algorithm after any number $t$ of decimation steps can be described as follows. Pick a uniformly random solution $\underline{x}^{*}$, and a uniformly random subset of the variable indexes $U_{t}\subseteq V$, with $|U_{t}|=t$. Then fix the variables $i\in U_{t}$ to take value $x_{i}=x_{i}^{*}$, and discard the rest of the reference configuration $\underline{x}^{*}$ (i.e. the bits $x_{j}^{*}$ for $j\notin U_{t}$).

We now put aside the idealized algorithm and consider the effect of fixing the $t$-th variable $i(t)$ to $x_{i(t)}^{*}$. Three cases can in principle arise: $(i)$ $x_{i(t)}$ was directly implied to be equal to $1-x_{i(t)}^{*}$ by $\underline{x}^{*}_{U_{t-1}}$ and a contradiction is generated. We assume that BP is able to detect this direct implication and avoid such a trivial contradiction; $(ii)$ $x_{i(t)}$ was directly implied to $x_{i(t)}^{*}$ by $\underline{x}^{*}_{U_{t-1}}$. The set of frozen variables remains the same, $W_{t}=W_{t-1}$, as this step is merely the actuation of a previous logical implication; $(iii)$ $i(t)$ was not directly implied by $\underline{x}^{*}_{U_{t-1}}$. This is the only interesting case that we develop now.

Let us call $Z_{t}\equiv W_{t}\setminus W_{t-1}$ the set of newly frozen variables after this fixing step. A moment of reflection shows that $Z_{t}$ contains $i(t)$ and that it forms a connected subset of $V$ in $G$. Consider now the subgraph $G_{t}\subseteq G$ induced by $Z_{t}$ (i.e. $G_{t}=(Z_{t},F_{t},E_{t})$ where $F_{t}$ is the set of factor nodes having at least one adjacent variable in $Z_{t}$, and $E_{t}$ is the set of edges between $Z_{t}$ and $F_{t}$). A crucial observation is the following:

If $G_{t}$ is a tree, no contradiction can arise during step $t$.

From this lemma, and since the factor graph of a typical random formula is locally tree-like, one is naturally lead to study the size of $Z_{t}$, i.e. of the cascade of newly implied variables induced by fixing the $t$-th variable. If this size remains bounded as $n\to\infty$, then $G_{t}$ will typically be a tree and, consequently, contradictions will arise with vanishingly small probability during one step. If on the other hand the size diverges for infinitely large samples, then $G_{t}$ will contain loops and opens the possibility for contradictions to appear.

Of course $\phi$ will be an increasing function of $\theta$. The argument above implies that, as long as its derivative remain finite for $\theta\in$, then the algorithm finds a solution. When the derivative diverges at some point $\theta_{*}$, then the number of direct implications of a single variable diverges as well. The spinodal point $\alpha_{{\mbox{\footnotesize spin}}}(k)$ is defined as be the smallest value of $\alpha$ such that this happens.

The expectation in the definition of $\phi(\theta)$ is with respect to the choices made by the real BP algorithm in the first $n\theta$ steps, including the small mistakes it necessarily makes. Our crucial hypothesis is that the location of $\alpha_{{\mbox{\footnotesize spin}}}(k)$ does not change (in the $n\to\infty$ limit) if $\phi(\theta)$ is computed along the execution of the idealized decimation algorithm. In other words we assume that the cumulative effect of BP errors over $n$ decimation steps produces only a small bias in the distribution of $\underline{x}^{*}$. For $\alpha\geq\alpha_{{\mbox{\footnotesize BPd}}}(k)$ this hypothesis is no longer consistent, as the real BP algorithm fails with high probability.

IV The tree model and its analysis

Let us consider a $k$-satisfiability formula whose factor graph is a finite tree, and the uniform measure $\mu$ over its solutions (which always exist) defined in Eq. (1). It follows from general results that the recursion equations (2,3) have a unique fixed-point, that we shall denote $\{h_{i\to a},u_{a\to i}\}$. Further the BP marginals $\nu_{i}(\,\cdot\,)$, cf. Eq. (5), are the actual marginals of $\mu$. Drawing a configuration $\underline{x}$ from the law $\mu$ is most easily done in a recursive, broadcasting fashion. Start from an arbitrary variable node $i$ and draw $x_{i}$ with distribution $\nu_{i}$. Thanks to the Markov property of $\mu$, conditional on the value of $x_{i}$, $\underline{x}_{V\setminus i}$ can be generated independently for each of the branches of the tree rooted at $i$. Namely, for each $a\in{\partial i}$, one draws $\underline{x}_{{\partial a}\setminus i}$ from

Here $z$ is a normalization factor and $\nu_{i\to a}(\cdot)$ denotes the marginal of the variable $x_{i}$ in the amputated factor graph where factor node $a$ has been removed (this is easily expressed in terms of the message $h_{i\to a}$). Once all variables $j$ at distance $1$ from $i$ have been generated, the process can be iterated to fix variables at distance $2$ from $i$, and so on. It is easy to realize that this process indeed samples a solution uniformly at random.

Following the program sketched in the previous Section, we shall study the effect of fixing a subset of the variables to the value they take in one of the solutions. We first state the following lemma.

Suppose $U$ is a subset of the variables of a tree formula, and let $\underline{x}^{*}$ be a uniformly random solution. The probability that a variable $i\notin U$ is directly implied by $\underline{x}_{U}^{*}$ reads

where the new messages $\{\widehat{u}_{a\to i},\widehat{h}_{i\to a}\}$ are solutions of

We consider now a random tree factor graph and a random set of fixed variables $U$.

We expect the tree computation to provide the correct prediction for the actual curve $\phi(\theta)$ (i.e. $\phi^{\mbox{\footnotesize\rm tree}}(\theta)=\phi(\theta)$) for a large range of the satisfiable regime, including $[0,\alpha_{{\mbox{\footnotesize spin}}}^{\mbox{\footnotesize\rm tree}}(k)]$. As a consequence, we expect $\alpha_{{\mbox{\footnotesize spin}}}(k)=\alpha_{{\mbox{\footnotesize spin}}}^{{\mbox{\footnotesize\rm tree}}}(k)$ and BP decimation to be successful up to $\alpha_{{\mbox{\footnotesize spin}}}(k)$. Similar tree computations are at the basis of a number of statistical mechanics computations in random $k$-SAT and have been repeatedly confirmed by rigorous studies.

Local weak limits of combinatorial models on random graphs were recently considered in . For a generalized conjecture in the regime $[\alpha_{\rm c}(k),\alpha_{\rm s}(k)]$ see .

A slightly stronger version of this conjecture would imply that $\phi(\theta)=\phi^{{\mbox{\footnotesize\rm tree}}}(\theta)$. As a consequence (following the discussion in previous section) the tree model would correctly describe the algorithm evolution.

V Numerical simulations

In order to test the validity of our analysis we performed numerical simulations of the pseudo-code of Table I. Let us give a few further details on its implementation. The BP messages are stored as $\{\tanh h_{i\to a},\tanh u_{a\to i}\}$. Ambiguities in the update rule (3) arises when $\tanh u_{b\to i}=\tanh u_{c\to i}=1$ with $b\in{\partial_{+}i}(a)$ and $c\in{\partial_{-}i}(a)$. Because of numerical imprecisions this situation can occur even before a contradiction has been detected by WP; such ambiguities are resolved by recomputing the incoming messages $\tanh u_{b\to i}$ using the regularized version of Eq. (4) with a small positive value of $\epsilon$ (in practice we used $\epsilon=10^{-4}$).

As for the stopping criterion used in step 5, we leave the BP iteration loop if either of the two following criteria is fulfilled: $(1)$, $\sup_{i}|\tanh h_{i}^{(r)}-\tanh h_{i}^{(r-1)}|<\delta$, i.e. BP has converged to a fixed-point within a given accuracy; $(2)$ A maximal number of iterations $r_{\rm max}$ fixed beforehand has been reached. In our implementation we took $\delta=10^{-10}$ and $r_{\rm max}=200$.

A first numerical check is presented in Fig. 2. The two dashed curves represent the fraction of frozen variables along the execution of the BP guide decimation algorithm, for two formulas of the $4$-sat ensemble, of moderate size ($n=4000$). The first formula had a ratio of constraints per variable $\alpha=7<\alpha_{{\mbox{\footnotesize spin}}}$. In agreement with the picture obtained from the analytical computation, the algorithm managed to find a solution of the formula (no contradiction encountered) and the measured fraction of frozen variables follows quite accurately the tree model prediction. The second formula was taken in the regime $\alpha_{\mbox{\footnotesize spin}}<\alpha=8.4<\alpha_{\rm c}$. The algorithm halted because a contradiction was found, after roughly the fraction $\theta_{*}$ (computed from the tree model) of variables has been fixed. The portion of the curve before this event exhibits again a rather good agreement between the direct simulation and the model.

Figure 3 shows the probability of success of BP decimation in a neighborhood of $\alpha_{{\mbox{\footnotesize spin}}}(4)$ for random formulae of size $n=500$, $1000$, $2000$. Each data point is obtained by running the algorithm on $1000$ to $3000$ formulae. The data strongly suggest that the success probability is bounded away from $0$ for $\alpha<\alpha_{\mbox{\footnotesize spin}}(k)$, in agreement with our argument.

Finally, in Figure 4 we consider the number of variables $t_{*}$ fixed by BP decimation before a contradiction is encountered. According to the argument in Section III-D, $t_{*}/n$ should concentrate around the location $\theta_{*}$ of the discontinuity in $\phi(\theta)$. This is in fact the point at which the number of variables directly implied by a fixed one is no longer bounded. The comparison is again encouraging. Notice that for $\alpha<\alpha_{\mbox{\footnotesize spin}}(k)$ we do not have any prediction, and the estimate of $t_{*}$ concerns only a small fraction of the runs.

To summarize, our simulations support the claim that, for $\alpha<\alpha_{\mbox{\footnotesize spin}}(k)$ the success probability is strictly positive and the algorithm evolution follows the tree model. For $\alpha>\alpha_{{\mbox{\footnotesize spin}}}(k)$ the main failure mechanism is indeed related to unbounded cascades of directly implied variables, after about $n\theta_{*}$ steps.

VI Conclusions and future directions

Let us conclude by highlighting some features of this work and proposing some directions for future research. It is worth mentioning that, as was also found in , random 3-sat has a qualitatively different behavior compared to random $k$-sat with $k\geq 4$. In particular we did not found any evidence for the existence of a vertical tangent point in the $k=3$ function $\phi(\theta)$ in the regime we expect to control through the tree computation, namely $\alpha<\alpha_{\rm c}(3)\approx 3.86$.

Our analysis suggests that BP guided decimation is successful with positive probability for $\alpha\leq\alpha_{{\mbox{\footnotesize spin}}}(k)$. Further we argued that this threshold can be computed through a tree model and evaluated via density evolution. Despite these conclusions are based on several assumptions, it is tempting to make a comparison with the best rigorous results on simple decimation algorithms. For $k=4$ the best result was obtained by Frieze and Suen who proved SCB (shortest clause with limited amount of backtracking) to succeed for $\alpha<5.54$. This is far from the conjectured threshold of BP decimation that is $\alpha_{{\mbox{\footnotesize spin}}}(4)\approx 8.05$. For large $k$, an asymptotic expansion suggests that

whereas SCB is known from to reach clause densities of $c_{k}2^{k}/k$, with $c_{k}\to 1.817$ as $k\to\infty$. A rigorous version of our analysis would lead to a constant factor improvement. On the other hand, the quest for an algorithm that provably solves random $k$-SAT in polynomial time beyond $\alpha=O(2^{k}/k)$, is open.

From a practical point of view the decimation strategy studied in this paper is not the most efficient one. A seemingly slight modification of the pseudo-code of Table I consists in replacing the uniformly random choice of the variable to be fixed, privilegiating the ones with the most strongly biased marginals. The intuition for this choice is that these marginals are the less subject to the ’small errors’ of BP. The numerical results reported in suggest that this modification improves significantly the performances of the decimation algorithm; unfortunately it also makes its analysis much more difficult.

This work was partially supported by EVERGROW, integrated project No. 1935 in the complex systems initiative of the Future and Emerging Technologies directorate of the IST Priority, EU Sixth Framework.

References