format: text

text:
# Survey propagation: an algorithm for satisfiability

## I Introduction

The satisfiability problem is the archetype of combinatorial optimization problems which are well known to be intractable in the worst case. However, experimental studies show that many instances of satisfiability are surprisingly easy, even for naive heuristic algorithms. In an attempt to get a better understanding of which instances are easy or hard to solve, a lot of efforts have focused in recent years on the ’random K-sat’ problemCook\_review . Instances of this problem are generated by considering $`N`$ variables and $`M=\alpha N`$ clauses, where each clause contains exactly $`K`$ distinct variables, and is picked up with uniform probability distribution from the set of $`\left(\begin{array}{c}N\\ K\end{array}\right)2^K`$ possible clauses. For a given value of $`\alpha `$, the probability $`P_N(\alpha )`$ that a randomly generated instance is SAT is a decreasing function, with $`P_N(0)=1`$ and $`lim_\alpha \mathrm{}P_N(\alpha =0)`$, which has been shown to approach, as $`N`$ increases, a step function characteristic of a zero-one law friedgut , or a ‘phase transition’. It is convenient to identify a crossover regime between ’SAT’ and ’UNSAT’ regimes using the value $`\alpha _c(N)`$ of the number of constraints per variable where $`P_N(\alpha _c(N))=1/2`$. From numerical simulations, $`\alpha _c(N)`$ is supposed to converge, in the large $`N`$ limit, to a value around $`\alpha _c4.27`$ TCS\_issue ; KirkSel ; Crawford ; MZKST , but this convergence has not yet been established rigorously. Interestingly, the performance of algorithms is found to be much worse around this value of $`\alpha =4.27`$: randomly generated instances with $`\alpha `$ near to the phase transition are particularly difficult to solve.

Rigorous lower and upper bounds have been found for this conjectured satisfiability threshold: it has been established that $`lim_N\mathrm{}P_N(\alpha )=1`$ for $`\alpha <\alpha _{lb}`$ and $`lim_N\mathrm{}P_N(\alpha )=0`$ for $`\alpha >\alpha _{ub}`$. The present best bounds for the case of the random 3-SAT problem (with $`K=3`$) are $`\alpha _{ub}=4.506`$ (from dubois , using the first moment method) and $`\alpha _{lb}=3.42`$ (from Kirousis\_lb , using algorithmic analysis). Note also the interesting algorithm-independent upper bound found in moore\_achl ; achl\_perez  using the second moment method, which becomes better for larger values of $`K`$.

Recently, some elaborate statistical physics methods have been brought to bear on the random satisfiability problem. These non-rigorous analytical calculations have put forward some interesting conjectures about what happens in the solution space of the problem as this threshold is approachedMEPAZE ; MZ\_pre  (see also previous work in MonZec ; Biroli ). They suggest the following overall picture, which should hold for a generic sample of the random satisfiability problem, in the limit $`N\mathrm{}`$, with $`\alpha `$ fixed.

1. There exists a SAT-UNSAT phase transition at a critical value $`\alpha _c`$ which can be computed by solving some (complicated) integral equation; for $`K=3`$, one gets $`\alpha _c4.267`$
2. There exists a second threshold $`\alpha _{clust}`$ separating two phases which are both ‘SAT’ (in each of them there exists a satisfiable assignment with probability 1), but with very different geometric structures:
3. For $`\alpha <\alpha _{clust}`$, a generic problem has many solutions, which tend to form one giant “cluster”; the set of all satisfying assignments forms a connected cluster in which it is possible to find a path between two solutions that requires short steps only (each pair of consecutive assignments in the path are close together in Hamming distance). In this regime, local search algorithms and other simple heuristics can relatively easily find a solution. This region is called the ’easy-SAT’ region
4. For $`\alpha _{clust}<\alpha <\alpha _c`$, there exists a ’hard SAT’ phase where the solution space breaks up into many smaller clusters. Solutions in separate clusters are generally far apart: it is not possible to transform a SAT assignment in one cluster into another one in a different cluster by changing only a finite number of variables. Because of this clustering effect, local search algorithms tend to have a very slow convergence when applied to large $`N`$ instances.

So far, the analytic method used in the most recent statistical physics analysis, named the cavity methodBethe\_cav , is non rigorous, and turning this type of approach into a rigorous theory is an open subject of current researchtalag ; FraLeo . Note however that, in the simpler case of the random K-XOR-SAT, the validity of this statistical physics analysis can be confirmed by rigorous studies xorsat1 ; xorsat2 , and the above clustering conjecture has been fully confirmed.

Interestingly, the statistical physics analysis suggests a new efficient heuristic algorithm for finding SAT assignments in the hard SAT phase, which has been put forward by two of us in MZ\_pre . The aim of this paper is to provide a detailed self-contained description of this algorithm, which does not rely on the statistical physics background. We shall limit the description to the regime where solutions exist, the so called SAT phase; some modification of the algorithm allows to address the optimization problem of minimizing the number of violated constraints in the UNSAT phase, but it will not be discussed here.

The basic building block of the algorithm, called survey propagation (SP), is a message passing procedure which resembles in some respect the iterative algorithm known as belief propagation (BP), but with some crucial differences which will be described. BP is a generic algorithm for computing marginal probability distributions in problems defined on factor graphs, which has been very useful in the context of error correcting codes Gallager  and Bayesian networks pearl .

While in simple limits we are able to give some rigorous results together with an explicit comparison with the belief propagation procedures, in general there exists no rigorous proof of convergence of the algorithm. However, we provide clear numerical evidence of its performance over benchmarks problems which appear to be far larger than those which can be handled by present state-of-the-art algorithms.

The paper is organized as follows: Sect. II describes the satisfiability problem and its graphical representation in terms of a factor graph. Sect. III explains two message passing algorithms, namely warning propagation (WP) and belief propagation (BP). Both are exact for tree factor graphs. Even if they are typically unable to find a solution for random SAT in the “interesting” hard-SAT region, they are shown here because they are in some sense the basic building blocks of our survey propagation algorithm. Sect. IV explains the survey propagation algorithm itself, a decimation procedure based on it, the ‘survey inspired decimation’ (SID), and the numerical results. In Sect.V we give some heuristic arguments from statistical physics which may help the reader to understand where the SP algorithm comes from. Sect. VI contains a few general comments.

## II The SAT problem and its factor graph representation

We consider a satisfiability problem consisting of $`N`$ Boolean variables $`\{x_i\{0,1\}\}`$ (where $`\{0,1\}\{F,T\}`$), with $`i\{1,\mathrm{},N\}`$, with $`M`$ constraints. Each constraint is a clause, which is the logical OR of the variables or of their negations. A clause $`a`$ is characterized by the set of variables $`i_1,\mathrm{},i_K`$ which it contains, and the list of those which are negated, which can be characterized by a set of $`K`$ numbers $`J_{i_r}^a\{\pm 1\}`$ as follows. The clause is written as

$$\left(z_{i_1}\mathrm{}z_{i_r}\mathrm{}z_{i_K}\right)$$
(1)
where $`z_{i_r}=x_{i_r}`$ if $`J_{i_r}^a=1`$ and $`z_{i_r}=\overline{x}_{i_r}`$ if $`J_{i_r}^a=1`$ (note that a positive literal is represented by a negative $`J`$). The problem is to find whether there exists an assignment of the $`x_i\{0,1\}`$ which is such that all the $`M`$ clauses are true. We define the total cost $`C`$ of a configuration $`𝐱=(x_1,\mathrm{},x_N)`$ as the number of violated clauses.

In what follows we shall adopt the factor graph representation factor\_graph  of the SAT problem. This representation is convenient because it provides an easy graphical description to the message passing procedures which we shall develop. It also applies to a wide variety of different combinatorial problems, thereby providing a unified notation.

The SAT problem can be represented graphically as follows (see fig.1). Each of the $`N`$ variables is associated to a vertex in the graph, called a “variable node” (circles in the graphical representation), and each of the $`M`$ clauses is associated to another type of vertex in the graph, called a “function node” (squares in the graphical representation). A function node $`a`$ is connected to a variable node $`i`$ by an edge whenever the variable $`x_i`$ (or its negation) appears in the clause $`a`$. In the graphical representation, we use a full line between $`a`$ and $`i`$ whenever the variable appearing in the clause is $`x_i`$ (i.e. $`J_i^a=1`$), a dashed line whenever the variable appearing in the clause is $`\overline{x}_i`$ (i.e. $`J_i^a=1`$). Variable nodes compose the set $`X`$ ($`|X|=N`$) and function nodes the set $`A`$ ($`|A|=M`$).

In summary, each SAT problem can be described by a bipartite graph, $`G=\left(XA;E=X\times A\right)`$ where $`E`$ is the edge set, and by the set of “couplings” $`\{J_a^i\}`$ needed to define each function node. For the K-SAT problem where each clause contains $`K`$ variables, the degree of all the function nodes is $`K`$.

Throughout this paper, the variable nodes indices are taken in $`i,j,k,\mathrm{}`$, while the function nodes indices are taken in $`a,b,c,\mathrm{}`$. For every variable node $`i`$, we denote by $`V(i)`$ the set of function nodes $`a`$ to which it is connected by an edge, by $`n_i=|V(i)|`$ the degree of the node, by $`V_+(i)`$ the subset of $`V(i)`$ consisting of function nodes $`a`$ where the variable appears un-negated (the edge $`(a,i)`$ is a full line), and by $`V_{}(i)`$ the complementary subset of $`V(i)`$ consisting of function nodes $`a`$ where the variable appears negated (the edge $`(a,i)`$ is a dashed line). $`V(i)b`$ denotes the set V(i) without a node $`b`$. Similarly, for each function node $`a`$, we denote by $`V(a)=V_+(a)V_{}(a)`$ the set of neighboring variable nodes, decomposed according to the type of edge connecting $`a`$ and $`i`$, and by $`n_a`$ the degree. Given a function node $`a`$ and a variable node $`j`$, connected by an edge, it is also convenient to define the two sets: $`V_a^u(j)`$ and $`V_a^s(j)`$, where the indices $`s`$ and $`u`$ respectively refer to the neighbors which tend to make variable $`j`$ satisfy or unsatisfy the clause $`a`$, defined as (see fig.2):

if $`J_j^a=1:V_a^u(j)=V_+(j);V_a^s(j)=V_{}(j)a`$ (2)
if $`J_j^a=1:V_a^u(j)=V_{}(j);V_a^s(j)=V_+(j)a`$ (3)
The same kind of factor graph representation can be used for other constraint satisfaction problems, where each function node $`a`$ defines an arbitrary function over the set $`X_aX`$ of variable nodes to which is connected, and could also involve hidden variables.

## III The message passing solution of SAT on a tree

In the special case in which the factor graph of a SAT problem is a tree (we shall call it a tree-problem), the satisfiability problem can be easily solved by many methods. Here we shall describe two message passing algorithms. The first one, called warning propagation (WP), determines whether a tree-problem is SAT or not; if it is SAT, WP finds one satisfiable assignment. The second algorithm, called belief propagation (BP), computes the number of satisfiable assignments, as well as the fraction of these assignments where a given variable is set to true. These algorithms are exact for tree-problems, but they can be used as heuristic in general problems, and we first give their general definition, which does not rely on the tree-like structure of the factor graph.

### III.1 Warning propagation

The basic elementary message passed from one function node $`a`$ to a variable $`i`$ (connected by an edge) is a Boolean number $`u_{ai}\{0,1\}`$ called a ‘warning’.

The update rule is defined as follows. Given a function node $`a`$ and one of its variables nodes $`i`$, the warning $`u_{ai}`$ is determined from the warnings $`u_{bj}`$ arriving on all the variables $`jV(a)i`$ according to:

$$u_{ai}=\underset{jV(a)i}{}\theta \left(J_j^a\left(\underset{bV(j)a}{}J_j^bu_{bj}\right)\right),$$
(4)
where $`\theta (x)=0`$ if $`x0`$ and $`\theta (x)=1`$ if $`x>0`$. This update rule is used sequentially, resulting in the following algorithm:

WP algorithm

INPUT: the factor graph of a Boolean formula in conjunctive normal form; a maximal number of iterations $`t_{max}`$

OUTPUT: UN-CONVERGED if WP has not converged after $`t_{max}`$ sweeps. If it has converged: the set of all warnings $`u_{ai}^{}`$.

* At time $`t=0`$: For every edge $`ai`$ of the factor graph, randomly initialize the warnings $`u_{ai}(t=0)\{0,1\}`$, e.g. with probability $`1/2`$.
* For $`t=1`$ to $`t=t_{max}`$:

+ sweep the set of edges in a random order, and update sequentially the warnings on all the edges of the graph, generating the values $`u_{ai}(t)`$, using subroutine WP-UPDATE.
    + If $`u_{ai}(t)=u_{ai}(t1)`$ on all the edges, the iteration has converged and generated $`u_{ai}^{}=u_{ai}(t)`$: go to 2.
* If $`t=t_{max}`$ return UN-CONVERGED. If $`t<t_{max}`$ return the set of fixed point warnings $`u_{ai}^{}=u_{ai}(t)`$

Subroutine WP-UPDATE$`(u_{ai})`$

INPUT: Set of all warnings arriving onto each variable node $`jV(a)i`$

OUTPUT: new value for the warning $`u_{ai}`$.

* For every $`jV(a)i`$, compute the cavity field $`h_{ja}=\left(_{bV_+(j)a}u_{bj}\right)\left(_{bV_{}(j)a}u_{bj}\right)`$ (If $`V(j)a`$ is empty, then $`h_{ja}=0`$).
* Using these cavity fields $`h_{ja}`$, compute the warning $`u_{ai}=_{jV(a)i}\theta \left(h_{ja}J_j^a\right)`$ (If $`V(a)i`$ is empty, then $`u_{ai}=1`$).

The interpretation of the messages and the message-passing procedure is the following. A warning $`u_{ai}=1`$ can be interpreted as a message sent from function node $`a`$, telling the variable $`i`$ that it should adopt the correct value in order to satisfy clause $`a`$. This is decided by $`a`$ according to the messages which it received from all the other variables $`j`$ to which it is connected: if $`\left(_{bV(j)a}J_j^bu_{bj}\right)J_j^a<0`$, this means that the tendency for site $`j`$ (in the absence of $`a`$) would be to take a value which does not satisfy clause $`a`$. If all neighbors $`jV(a)i`$ are in this situation, then $`a`$ sends a warning to $`i`$. An example of the use of WP is shown in fig. 3.

The warning propagation algorithm can be applied to any SAT problem. When it converges, this dynamics defines a fixed point, which is a set of warnings $`u_{ai}^{}`$. These can be used to compute, for each variable $`i`$, the “local field” $`H_i`$ and the “contradiction number” $`c_i`$ which are two integers defined as:

$$H_i=\underset{bV(i)}{}J_i^bu_{bi}^{}$$
(5)
$`c_i`$ $`=`$ $`1\text{if}\left({\displaystyle \underset{bV_+(i)}{}}u_{bi}^{}\right)\left({\displaystyle \underset{bV_{}(i)}{}}u_{bi}^{}\right)>0`$ (6)
$`c_i`$ $`=`$ $`0\text{otherwise.}`$ (7)
The local field $`H_i`$ is an indication of the preferred state of the variable $`i`$: $`x_i=1`$ if $`H_i>0`$, $`x_i=0`$ if $`H_i<0`$. The contradiction number indicates whether the variable i has received conflicting messages.

The interest in WP largely comes from the fact that it gives the exact solution for tree-problems. This is summarized in the following simple theorem:

THEOREM 1:

Consider an instance of the SAT problem for which the factor graph is a tree. Then the WP algorithm converges to a unique set of fixed point warnings $`u_{ai}^{}`$, independently on the initial warnings. If at least one of the corresponding contradiction numbers $`c_i`$ is equal to $`1`$, the problem is UNSAT, otherwise it is SAT.

Corollary: In the case where the problem is SAT, the local fields $`H_i`$ can be used to find an assignment of the variables satisfying all the clauses, using the following algorithm called “Warning Inspired Decimation” or WID:
WID algorithm   INPUT: the factor graph of a Boolean formula in conjunctive normal form   OUTPUT: UN-CONVERGED, or status of the formula (SAT or UNSAT); If the formula is SAT: one assignment which satisfies all clauses.          While the number of unfixed variables is $`>0`$, do:    Run WP    If WP does not converge, return UN-CONVERGED. Else compute the local fields $`H_i`$ and the contradiction numbers $`c_i`$, using eqs. (5,7).    If there is at least one contradiction number $`c_i=1`$, return UNSAT. Else:    If there is at least one local field $`H_i0`$: fix all variables with $`H_i0`$ ($`H_i>0x_i=1`$ and $`H_i<0x_i=0`$ ), and clean the graph, which means: $`\{`$ remove the clauses satisfied by this fixing, reduce the clauses that involve the fixed variable with opposite literal, update the number of unfixed variables$`\}`$. GOTO label 1. Else:    Choose one unfixed variable, fix it to an arbitrary value, clean the graph. GOTO label 1        return the set of assignments for all the variables.

PROOF of theorem 1:

The convergence of message passing procedures on tree graphs is a well known result (see e.g.factor\_graph ). We give here an elementary proof of convergence for the specific case of WP, and then show how the results on $`H_i`$ and $`c_i`$ follow.

Call $``$ the set of nodes. Define the leaves of the tree, as the nodes of degree $`1`$. For any edge $`(a,i)`$ connecting a function node $`a`$ to a variable node $`i`$, define its level $`r`$ as follows: remove the edge $`(a,i)`$ and consider the remaining subgraph containing $`a`$. This subgraph $`𝒯_{ai}`$ is a tree factor graph defining a new SAT problem. The level $`r`$ is the maximal distance between $`a`$ and all the leaves in the subgraph $`𝒯_{(a,i)}`$ (the distance between two nodes of the graph is the number of edges of the shortest path connecting them). If an edge $`(a,i)`$ has level $`r=0`$ (which means that $`a`$ is a leaf of the subgraph), $`u_{ai}(t)=1`$ for all $`t1`$. If $`(a,i)`$ has level $`r=1`$, then $`u_{ai}(t)=0`$ for all $`t1`$. From the iteration rule, a warning $`u_{ai}`$ at level $`r`$ is fully determined from the knowledge of all the warnings $`u_{bj}`$ at levels $`r2`$. Therefore the warning $`u_{ai}(t)`$ at a level $`r`$ is guaranteed to take a fixed value $`u_{ai}^{}`$ for $`t1+r/2`$.

Let us now turn to the study of local fields and contradiction numbers.

We first prove the following lemma:

If a warning $`u_{ai}^{}=1`$, the clause $`a`$ is violated in the reduced SAT problem defined by the subgraph $`𝒯_{ai}`$.

This is obviously true if the edge has level $`r=0`$ or $`r=1`$. Supposing that it holds for all levels $`r2`$, one considers an edge $`ai`$ at level $`r`$, with $`u_{ai}^{}=1`$. From (4), this means that for all variable nodes $`jV(a)i`$, the node $`j`$ receives at least one message from a neighboring factor node $`bV(j)a`$ with $`u_{bj}=1`$. The edge $`bj`$ is at level $`r2`$, therefore the reduced problem on the graph $`𝒯_{bj}`$ is UNSAT and therefore the clause $`b`$ imposes the value of the variable $`j`$ to be $`1`$ (True) if $`J_j^b=1`$, or $`0`$ (False) if $`J_j^b=1`$: we shall say that clause $`b`$ fixes the value of variable $`j`$. This is true for all $`jV(a)i`$, which means that the reduced problem on $`𝒯_{(a,i)}`$ is UNSAT, or equivalently, the clause $`a`$ fixes the value of variable $`i`$.

Having shown that $`u_{ai}^{}=1`$ implies that clause $`a`$ fixes the value of variable $`i`$, it is clear from (7) that a nonzero contradiction number $`c_i`$ implies that the formula is UNSAT.

If all the $`c_i`$ vanish, the formula is SAT. One can prove this for instance by showing that the WID algorithm generates a SAT assignment. The variables with $`H_i0`$ receive some nonzero $`u_{ai}^{}=1`$ and are fixed. One then ’cleans’ the graph, which means: remove the clauses satisfied by this fixing, reduce the clauses that involve the fixed variable with opposite literal. By definition, this process has removed from the graph all the edges on which there was a nonzero warning. So on the new graph, all the edges have $`u^{}=0`$. Following the step 2.2 of WID, one chooses randomly a variable $`i`$, one fixes it to an arbitrary value $`x_i`$, and cleans the graph. The clauses $`a`$ connected to $`i`$ which are satisfied by the choice $`x_i`$ are removed; the corresponding subgraphs are trees where all the edges have $`u^{}=0`$. A clause $`a`$ connected to $`i`$ which are not satisfied by the choice $`x_i`$ may send some $`u^{}=1`$ messages (this happens if such a clause had degree $`2`$ before fixing variable $`i`$). However, running WP on the corresponding subgraph $`𝒯_{(a,i)}`$, the set of warnings can not have a contradiction: A variable $`j`$ in this subgraph can receive at most one $`u^{}=1`$ warning, coming from the unique path which connects $`j`$ to $`a`$. Therefore $`c_j=0`$: one iteration of WID has generated a strictly smaller graph with no contradiction. By induction, it thus finds a SAT assignment.$`\mathrm{}`$

One should notice that the variables which are fixed at the first iteration of WID (those with non-zero $`H_i`$) are constrained to take the same value in all satisfiable assignments.

### III.2 Belief propagation

While the WP algorithm is well adapted to finding a SAT assignment, the more complicated belief propagation (BP) algorithm is able to compute, for satisfiable problems with a tree factor graph, the total number of SAT assignments, and the fraction of SAT assignments where a given variable $`x_i`$ is true.

We consider a satisfiable instance, and the probability space built by all SAT assignments taken with equal probability. Calling $`a`$ one of the clauses in which $`x_i`$ appears, the basic ingredients of BP are the messages:

* $`\mu _{ai}(x_i)[0,1]`$, interpreted as the probability that clause $`a`$ is satisfied, given the value of the variable $`x_i\{0,1\}`$.
* $`\mu _{ia}(x_i)[0,1]`$, interpreted as the probability that the variable takes value $`x_i`$, when clause $`a`$ is absent (this is again a typical ’cavity’ definition). Notice that $`_{x_i\{0,1\}}\mu _{ia}(x_i)=1`$, while there is no such normalization for $`\mu _{ai}(x_i)`$.

The BP equations are:

$$\mu _{ia}(x_i)=C_{ia}\underset{bV(i)a}{}\mu _{bi}(x_i),$$
(8)
$$\mu _{ai}(x_i)=\underset{\{x_j(ji)\}}{}f_a(X)\underset{jV(a)i}{}\mu _{ja}(x_j),$$
(9)
where $`C_{ia}`$ is a normalization constant ensuring that $`\mu _{ia}`$ is a probability, the sum over $`\{x_j(ji)\}`$ means a sum over all values of the variables $`x_j\{0,1\}`$, for all $`j`$ different from $`i`$, and $`f_a(X)`$ is a characteristic function taking value $`1`$ if the configuration $`X=\{x_i\}`$ satisfies clause $`a`$, taking value $`0`$ otherwise.

It is convenient to parameterize $`\mu _{ia}(x_i)`$ by introducing the number $`\gamma _{ia}[0,1]`$ which is the probability that the variable $`x_i`$ is in the state which violates clause $`a`$, in a problem where clause $`a`$ would be absent (writing for instance $`\mu _{ia}(x_i)=\gamma _{ia}\delta (x_i,0)+(1\gamma _{ia})\delta (x_i,1)`$ in the case where $`J_i^a=1`$).

Let us denote by

$$\delta _{ai}\underset{jV(a)i}{}\gamma _{ja}$$
(10)
the probability that all variables in clause $`a`$, except variable $`i`$, are in the state which violates the clause.

The BP algorithm amounts to an iterative update of the messages $`\delta _{ai}`$ according to the rule:

BP algorithm: INPUT: the factor graph of a Boolean formula in conjunctive normal form; a maximal number of iterations $`t_{max}`$; a requested precision $`ϵ`$.

OUTPUT: UN-CONVERGED if BP has not converged after $`t_{max}`$ sweeps. If it has converged: the set of all messages $`\delta _{ai}^{}`$.

* At time $`t=0`$: For every edge $`ai`$ of the factor graph, randomly initialize the messages $`\delta _{ai}(t=0)[0,1]`$
* For $`t=1`$ to $`t=t_{max}`$:

+ sweep the set of edges in a random order, and update sequentially the warnings on all the edges of the graph, generating the values $`\delta _{ai}(t)`$, using subroutine BP-UPDATE.
    + If $`|\delta _{ai}(t)\delta _{ai}(t1)|<ϵ`$ on all the edges, the iteration has converged and generated $`\delta _{ai}^{}=\delta _{ai}(t)`$: go to 2.
* If $`t=t_{max}`$ return UN-CONVERGED. If $`t<t_{max}`$ return the set of fixed point warnings $`\delta ai^{}=\delta _{ai}(t)`$

Subroutine BP-UPDATE$`(\delta _{ai})`$

INPUT: Set of all messages arriving onto each variable node $`jV(a)i`$

OUTPUT: new value for the message $`\delta _{ai}`$.

* For every $`jV(a)i`$, compute the cavity field

$$\gamma _{ja}=\frac{P_{ja}^u}{P_{ja}^u+P_{ja}^s}$$
(11)
where

$`P_{ja}^u`$ $`=`$ $`{\displaystyle \underset{bV_a^s(j)}{}}(1\delta _{bj}),`$
$`P_{ja}^s`$ $`=`$ $`{\displaystyle \underset{bV_a^u(j)}{}}(1\delta _{bj}).`$ (12)
If an ensemble is empty, for instance $`V_a^s(j)=\mathrm{}`$, the corresponding $`P_{ja}^u`$ takes value $`1`$ by definition.
* Using these numbers $`\gamma _{ja}`$, compute the new message: $`\delta _{ai}_{jV(a)i}\gamma _{ja}`$. If a factor node $`a`$ is a leaf (unit clause) with a single neighbor $`i`$, the corresponding $`\delta _{ai}`$ takes value $`1`$ by definition.

As WP, the BP algorithm is exact on trees (see for instance factor\_graph ). In fact it gives a more accurate results than WP since it allows to compute the exact probabilities $`\mu _i`$ (while WP identifies the variables which are fully constrained, and gives a zero local field on the other variables). A working example of BP is shown in fig.4. In this example and more in general for trees, BP also provides the exact number $`𝒩`$ of SAT assignments, as given by the following theorem:

THEOREM 2: Consider an instance of the SAT problem for which the factor graph is a tree, and there exist some SAT assignments. Then:

a) The BP algorithm converges to a unique set of fixed point messages $`\delta _{ai}^{}`$.

b)The probability $`\mu _i`$ that the variable $`x_i=1`$ is given by:

$$\mu _i=\frac{_{aV_{}(i)}(1\delta _{ai}^{})}{_{aV_{}(i)}(1\delta _{ai}^{})+_{aV_+(i)}(1\delta _{ai}^{})}.$$
(13)
c) The number of SAT assignments is $`𝒩=\mathrm{exp}(S)`$, where the entropy $`S`$ is given by:

$`S`$ $`=`$ $`{\displaystyle \underset{aA}{}}\mathrm{ln}\left[{\displaystyle \underset{iV(a)}{}}\left({\displaystyle \underset{bV_a^s(i)}{}}(1\delta _{bi}^{})+{\displaystyle \underset{bV_a^u(i)}{}}(1\delta _{bi}^{})\right){\displaystyle \underset{iV(a)}{}}\left({\displaystyle \underset{bV_a^u(i)}{}}(1\delta _{bi}^{})\right)\right]`$ (14)
$`+{\displaystyle \underset{iX}{}}(1n_i)\mathrm{ln}\left[{\displaystyle \underset{bV_+(i)}{}}(1\delta _{bi}^{})+{\displaystyle \underset{bV_{}(i)}{}}(1\delta _{bi}^{})\right]`$
PROOF:

The proof of convergence is simple, using the same strategy as the proof in sect.III: messages at level 0 and 1 are fixed automatically, and a message at level $`r`$ is fixed by the values of messages at lower levels.

The probability $`\mu _i`$ is computed from the same procedure as the one giving the BP equations (8,9), with the difference that one takes into account all neighbors of the site $`i`$.

The slightly more involved result is the one concerning the entropy. We use the probability measure $`P(X)`$ on the space of all assignments which has uniform probability for all SAT assignments and zero probability for all the assignments which violate at least one clause:

$$P(X)=\frac{1}{𝒩}\underset{aA}{}f_a(X).$$
(15)
From $`P`$ one can define the following marginals:

* the ’site marginal’ $`p_i(x_i)`$ is the probability that variable $`i`$ takes value $`x_i\{0,1\}`$
* the ’clause marginal’ $`p_a(X_a)`$ is the probability that the set of variables $`x_i,iV(a)`$, takes a given value, denoted by $`X_a`$ (among the $`2^{n_a}`$ possible values).

For a tree factor graph, one easily shows by induction on the size of the graph that the full probability can be expressed in terms of the site and clause marginals as:

$$P(X)=\underset{aA}{}p_a(X_a)\underset{iX}{}p_i(x_i)^{1n_i}.$$
(16)
The entropy $`S=\mathrm{ln}(𝒩)`$ is then obtained as

$$S=\underset{X}{}P(X)\mathrm{ln}P(X)=\underset{a}{}\underset{X_a}{}p_a(X_a)\mathrm{ln}[p_a(X_a)]\underset{i}{}(1n_i)\underset{x_i}{}p_i(x_i)\mathrm{ln}[p_i(x_i)]$$
(17)
Let us now derive the expression of this quantity in terms of the messages used in BP. One has:

$$p_i(x_i)=c_i\underset{bV(i)}{}\mu _{bi}(x_i)$$
(18)
and

$$p_a(X_a)=c_af_a(X_a)\underset{iV(a)}{}\mu _{ia}(x_i),$$
(19)
where $`c_i`$ and $`c_a`$ are two normalization constants. From (18) one gets after some reshuffling:

$$\underset{i}{}(n_i1)\underset{x_i}{}p_i(x_i)\mathrm{ln}p_i(x_i)=\underset{i}{}(n_i1)\mathrm{ln}c_i+\underset{a}{}\underset{iV(a)}{}\underset{X_a}{}p_a(X_a)\mathrm{ln}\left[\underset{bV(i)a}{}\mu _{bi}(x_i)\right];$$
(20)
Using the BP equation (8), this gives:

$$\underset{i}{}(n_i1)\underset{x_i}{}p_i(x_i)\mathrm{ln}p_i(x_i)=\underset{i}{}(n_i1)\mathrm{ln}c_i+\underset{a}{}\underset{X_a}{}p_a(X_a)\mathrm{ln}\left[\underset{iV(a)}{}\mu _{ia}(x_i)f_a(X_a)\right]\underset{a}{}\underset{iV(a)}{}\mathrm{ln}C_{ia},$$
(21)
where the term $`f_a(X_a)`$ inside the logarithm has been added, taking into account the fact that, as $`f_a(X_a)\{0,1\}`$, one always has $`p_a(X_a)\mathrm{ln}f_a(X_a)=0`$. Therefore:

$$S=\underset{a}{}\mathrm{ln}c_a+\underset{i}{}(n_i1)\mathrm{ln}c_i\underset{a}{}\underset{iV(a)}{}\mathrm{ln}C_{ia}$$
(22)
In the notations of (12), one has

$$c_a=\frac{1}{1_{iV(a)}\gamma _{ia}}=\frac{1}{1_{iV(a)}P_{ia}^u/(P_{ia}^u+P_{ia}^s)},$$
(23)
$$C_{ia}=\frac{1}{P_{ia}^u+P_{ia}^s},$$
(24)
and

$$c_i=\frac{1}{_{bV_+(i)}(1\delta _{bi})+_{bV_{}(i)}(1\delta _{bi})}.$$
(25)
Substitution into (22) gives the expression (14) for the entropy. $`\mathrm{}`$

## IV Survey Propagation

### IV.1 The algorithm

The WP and BP algorithms have been shown to work for satisfiability problems where the factor graph is a tree. In more general cases where the factor graph has loops, they can be tried as heuristics, but there is no guarantee of convergence. In this section we present a new message passing algorithm, survey propagation (SP), which is also a heuristic, without any guarantee of convergence. It reduces to WP for tree-problems, but it turns out to be more efficient than WP or BP in experimental studies of random satisfiability problems. SP has been discovered using concepts developed in statistical physics under the name of ’cavity method’. Here we shall first present SP, then give some experimental results, and in the end expose the qualitative physical reasoning behind it.

A message of SP, called a survey, passed from one function node $`a`$ to a variable $`i`$ (connected by an edge) is a real number $`\eta _{ai}[0,1]`$. The SP algorithm uses exactly the same main procedure as BP (see III.2), but instead of calling the BP-UPDATE, it uses a different update rule, SP-UPDATE, defined as:

Subroutine SP-UPDATE$`(\eta _{ai})`$

INPUT: Set of all messages arriving onto each variable node $`jV(a)i`$

OUTPUT: new value for the message $`\eta _{ai}`$.

* For every $`jV(a)i`$, compute the three numbers:

$`\mathrm{\Pi }_{ja}^u`$ $`=`$ $`\left[1{\displaystyle \underset{bV_a^u(j)}{}}\left(1\eta _{bj}\right)\right]{\displaystyle \underset{bV_a^s(j)}{}}\left(1\eta _{bj}\right)`$
$`\mathrm{\Pi }_{ja}^s`$ $`=`$ $`\left[1{\displaystyle \underset{bV_a^s(j)}{}}\left(1\eta _{bj}\right)\right]{\displaystyle \underset{bV_a^u(j)}{}}\left(1\eta _{bj}\right)`$
$`\mathrm{\Pi }_{ja}^0`$ $`=`$ $`{\displaystyle \underset{bV(j)a}{}}\left(1\eta _{bj}\right)`$ (26)
if a set like $`V_a^s(j)`$ is empty, the corresponding product takes value $`1`$ by definition.
* Using these numbers , compute and return the new survey:

$$\eta _{ai}=\underset{jV(a)i}{}\left[\frac{\mathrm{\Pi }_{ja}^u}{\mathrm{\Pi }_{ja}^u+\mathrm{\Pi }_{ja}^s+\mathrm{\Pi }_{ja}^0}\right].$$
(27)
If $`V(a)i`$ is empty, then $`\eta _{ai}=1`$.

Qualitatively, the statistical physics interpretation of the survey $`\eta _{ai}`$ is a probability that a warning is sent from $`a`$ to $`i`$ (see section V for details). Therefore, whenever the SP algorithm converges to a fixed-point set of messages $`\eta _{ai}^{}`$, one can use it in a decimation procedure in order to find a satisfiable assignment, if such an assignment exists. This procedure, called the survey inspired decimation (SID), is a generalization of the WID algorithm III.1, defined by:

SID algorithm

INPUT: The factor graph of a Boolean formula in conjunctive normal form. A maximal number of iterations $`t_{max}`$ and a precision $`ϵ`$ used in SP

OUTPUT: One assignment which satisfies all clauses, or ’SP UNCONVERGED’, or ’probably UNSAT’

* Random initial condition for the surveys
* Run SP. If SP does not converge, return ’SP UNCONVERGED’ and stop (or restart, i.e.go to 0.). If SP converges, use the fixed-point surveys $`\eta _{ai}^{}`$ in order to:
* Decimate:

+ If non-trivial surveys ($`\{\eta 0\}`$) are found, then:

    - Evaluate, for each variable node $`i`$, the three ’biases’ $`\{W_i^{(+)},W_i^{()},W_i^{(0)}\}`$defined by:

        $`W_i^{(+)}`$ $`=`$ $`{\displaystyle \frac{\widehat{\mathrm{\Pi }}_i^+}{\widehat{\mathrm{\Pi }}_i^++\widehat{\mathrm{\Pi }}_i^{}+\widehat{\mathrm{\Pi }}_i^0}}`$ (28)
        $`W_i^{()}`$ $`=`$ $`{\displaystyle \frac{\widehat{\mathrm{\Pi }}_i^{}}{\widehat{\mathrm{\Pi }}_i^++\widehat{\mathrm{\Pi }}_i^{}+\widehat{\mathrm{\Pi }}_i^0}}`$ (29)
        $`W_i^{(0)}`$ $`=`$ $`1W_i^{(+)}W_i^{()}`$ (30)
        where $`\widehat{\mathrm{\Pi }}_i^+,\widehat{\mathrm{\Pi }}_i^{},\widehat{\mathrm{\Pi }}_i^0`$ are defined by

        $`\widehat{\mathrm{\Pi }}_i^+`$ $`=`$ $`\left[1{\displaystyle \underset{aV_+(i)}{}}\left(1\eta _{ai}^{}\right)\right]{\displaystyle \underset{aV_{}(i)}{}}\left(1\eta _{ai}^{}\right)`$
        $`\widehat{\mathrm{\Pi }}_i^{}`$ $`=`$ $`\left[1{\displaystyle \underset{aV_{}(i)}{}}\left(1\eta _{ai}^{}\right)\right]{\displaystyle \underset{aV_+(i)}{}}\left(1\eta _{ai}^{}\right)`$
        $`\widehat{\mathrm{\Pi }}_i^0`$ $`=`$ $`{\displaystyle \underset{aV(i)}{}}\left(1\eta _{ai}^{}\right)`$ (31)
        - fix the variable with the largest $`|W_i^{(+)}W_i^{()}|`$ to the value $`x_i=1`$ if $`W_i^{(+)}>W_i^{()}`$, to the value $`x_i=0`$ if $`W_i^{(+)}<W_i^{()}`$. Clean the graph, which means: $`\{`$ remove the clauses satisfied by this fixing, reduce the clauses that involve the fixed variable with opposite literal, update the number of unfixed variables$`\}`$.
    + If all surveys are trivial ($`\{\eta =0\}`$) , then output the simplified sub-formula and run on it a local search process (e.g. walksat).
* If the problem is solved completely by unit clause propagation, then output “SAT” and stop. If no contradiction is found then continue the decimation process on the smaller problem (go to 1.) else (if a contradiction is reached) stop.

There exist several variants of this algorithm. In the code which is available at web , for performance reasons we update simultaneously all $`\eta `$ belonging to the same clause. The clauses to be updated are chosen in a random permutation order at each iteration step. The algorithm can also be randomized by fixing, instead of the most biased variables, one variable randomly chosen in the set of the x percent variables with the largest bias. This strategy allows to use some restart in the case where the algorithm has not found a solution. A fastest decimation can also be obtained by fixing in the step 2.1(b), instead of one variable, a fraction $`f`$ of the $`N_t`$ variables which have not yet been fixed (going back to $`1`$ variable when $`fN_t<1`$).

### IV.2 Experimental study of the SP algorithm

In order to get some concrete information on the behaviour of SP for large but finite $`N`$, we have experimented SP and SID on single instances of the random 3-SAT problem with many variables, up to $`N10^7`$. In this section we summarize these (single machine) experiments and their results.

Instances of the 3-SAT problem were generated with the pseudo random number generator "Algorithm B" on p.32 of Knuth Knuth . However we found that results are stable with respect to changes in the random number generators. Formulas are generated by choosing k-tuples of variable indices at random (with no repetitions) and by negating variables with probability $`0.5`$.

We first discuss the behaviour of the SP algorithm itself. We have used a precision parameter $`ϵ=10^3`$ (smaller values don’t seem to increase performance significantly). Depending on the range of $`\alpha `$, we have found the following behaviours, for large enough $`N`$:

* For $`\alpha <\alpha _d3.9`$, SP converges towards the set of trivial messages $`\eta _{ai}=0`$, for all $`ai`$ edges. All variables are under-constrained.
* For $`3.9<\alpha <4.3`$, SP converges to a unique fixed-point set of non-trivial messages, independently from the initial conditions, where a large fraction of the messages $`\eta _{ai}`$ are in $`]0,1[`$.

Notice that, for ‘small’ values of $`N`$, around $`N=1000`$, one often finds some instances in which SP does not converge. But the probability of convergence, at a given $`\alpha <4.3`$, increases with $`N`$. This is exemplified by the following quantitative measure of the performance of the SID algorithm (which uses SP). We have solved several instances of the random 3-SAT problem, for various values of $`\alpha `$ and $`N`$, using the SID algorithm in which we fix at each step the fraction $`fN_t`$ of variables with largest $`\left|W_i^{(+)}W_i^{()}\right|`$. Table 5 gives in each case the fraction of samples which are solved by SID, in a single run of decimation (without any restart). The algorithm fails when, either SP does not converge, or the simplified sub-formula found by SID is not solved by walksat. The performance of SID improves when $`N`$ increases and when $`f`$ decreases. Notice that for $`N=10^5`$ we solve all the $`50`$ randomly generated instances at $`\alpha =4.24`$. For larger values of $`\alpha `$ the algorithm often fails. Notice that in such cases it does not give any information on whether the instance is UNSAT. Some failures may be due to UNSAT instances, others are just real failures of the SID for SAT instances. Few experiments on even larger instances ($`N=10^6,10^7`$) have also been succesfully run (the limiting factor being the available computer memory needed to store formulas).

As shown by the data on the average total number of SP iterations along the successful solution process in table 5, the convergence time of the SP algorithm basically does not grow with $`N`$ (a growth like $`\mathrm{log}N`$, which could be expected from the geometrical properties of the factor graph, is not excluded). Therefore the process of computing all the SP messages $`\eta _{ai}^{}`$ takes $`\mathrm{\Theta }(N)`$, or maybe $`\mathrm{\Theta }(N\mathrm{ln}N)`$, operations. If SID fixes at each step only one variable, it will thus converge in $`\mathrm{\Theta }(N^2\mathrm{log}N)`$ operations (the time taken by walksat to solve the simplified sub-formula seems to grow more slowly). When we fix a fraction of variables at a time, we get a further reduction of the cost to $`O(N(\mathrm{ln}N)^2)`$ (the second $`\mathrm{ln}`$ comes from sorting the biases).

A very basic yet complete version of the code which is intended to serve only for the study on random 3-SAT instances is available at the web site web . Generalization of the algorithm to other problems require some changes which are not implemented in the distributed code.

## V Heuristic arguments

Survey propagation has been invented using powerful concepts and methods developed in the statistical physics of disordered systems, notably the cavity method for diluted problems Bethe\_cav . In this section we want to give some short background on these methods, in order to help the reader understand where SP comes from, and maybe develop similar algorithms in other contexts. Unfortunately so far there is no rigorous derivation of the cavity method, so this whole section only contains heuristic arguments.

### V.1 The physical picture underlying the SP construction: clustering of configurations

Let us start with a discussion of the validity of BP. As we saw in sect.III, BP aims at computing the marginal probability distribution of a variable $`x_i`$, within the probability space built by all SAT assignments, each being given equal probability. The message $`\mu _{ai}(x_i)`$ used in BP can be computed exactly if one knows the joint probability distribution $`P^{(a)}(X)`$ of the variables in $`X=\{x_j,jV(a)i\}`$, in the graph where clause $`a`$ is absent. Using the same notations as in (9), one has:

$$\mu _{ai}(x_i)=\underset{\{x_j(ji)\}}{}f_a(X)P^{(a)}(X)$$
(32)
Comparing this eq. (32) to the eqs. (8,9) of BP, one sees that BP uses an approximation, namely the fact that the joint probability $`P^{(a)}(X)`$ factorizes: $`P^{(a)}(X)_{jV(a)i}\mu _{ja}(x_j)`$. This amounts to assuming that the variables $`x_j`$, for $`jV(a)i`$, are uncorrelated in the absence of clause $`a`$. This assumption is obviously correct when the factor graph is a tree, which confirms the validity of BP in that case. In a general problem, for such an assumption to hold, we need two conditions to be fulfilled:

* the variables $`x_j,jV(a)i`$ should be far from each other, in the factor graph where clause $`a`$ is absent.
* there should be a single phase in the problem, or the probability measure should be reduced to one single pure phase.

The first condition is easily understood, and one can see that it is generically fulfilled when one considers a random satisfiability problem. In random K-sat with $`M=\alpha N`$ clauses, the factor graph is a random bipartite graph, where function nodes have degree $`K`$ and variable nodes have fluctuating degrees, with a distribution which becomes, in the large $`N`$ limit, a Poisson law of mean $`K\alpha `$. Locally such a graph is tree like. Taking a clause $`a`$ at random, the minimal distance between two of the variables $`x_j,jV(a)i`$ is generically of order $`\mathrm{log}N`$.

The second condition is more subtle. In general a pure phase is defined in statistical physics as an extremal Gibbs measure georgii ; however the standard construction of Gibbs measures deals with infinite systems. Here we need to work with $`N`$ variables where $`N1`$ but is finite, and the corresponding construction has not been worked out yet. For the random satisfiability problem, a heuristic description of a pure phase is a cluster of SAT assignments, defined as follows. Consider all SAT assignments. Define the distance between two assignments $`\{x_i\}`$ and $`\{y_i\}`$ as $`_i(x_iy_i)^2`$. If the distance between two SAT assignment is smaller than a number $`q`$, they are said to belong to the same $`q`$-cluster. This allows to partition the set of SAT assignments into $`q`$-clusters. One is interested in the ’clusters’ obtained in large $`N`$ – large $`q`$ limit of the $`q`$-clusters, where the large $`N`$ limit is taken first. Heuristic statistical physics arguments indicate that, in the random K-satisfiability problem, for $`\alpha <\alpha _{clust}(K)`$, there should exist one single such cluster of SAT assignments: this means that one can move inside the space of SAT assignment, from any assignment to another one, by a succession of moves involving each a number of flips of variables which is $`<<N`$. In such a case the BP factorization approximation is expected to be correct. On the other hand, for $`\alpha >\alpha _{clust}(K)`$ the space of SAT assignment separates into many distant clusters, and the BP factorization does not hold globally, but it would hold if one could restrict the probability space to one given cluster $`\alpha `$. Within such a restricted space, BP would converge to a set of messages $`\mu _{ai}^\alpha (x_i)`$ which depends on the cluster $`\alpha `$.

In this situation, the cavity method uses a statistical approach. It considers all the clusters of SAT assignments, and attributes to each cluster a probability proportional to the number of configuration that it contains. Then one introduces, on each edge $`ai`$, a survey which gives the probability, when a cluster $`\alpha `$ is chosen randomly with this probability, that the message $`\mu _{ai}^\alpha (x)`$ is equal to a certain function $`P(x)`$.

This object is a probability of a probability and it is thus difficult to use in practical algorithms. For this reason, SP departs from the usual cavity method and uses a simpler object which is a survey of warnings, interpreted as follows: Consider one cluster $`\alpha `$ and an edge $`ai`$ of the factor graph. If, in every SAT assignments of the cluster $`\alpha `$, all the variables $`x_j,jV(a)i`$ don’t satisfy clause $`a`$, then a warning $`u_{ai}^\alpha =1`$ is passed along the edge from $`a`$ to $`i`$. The SP message along this edge is the survey of these warnings, when one picks up a cluster $`\alpha `$ at random: $`\eta _{ai}=_\alpha u_{ai}^\alpha /(_\alpha 1)`$. So basically the SP message gives the probability that there is a warning sent from $`a`$ to $`i`$ in a randomly chosen cluster. With respect to the full-fledged cavity method, this is a much simplified object, which focuses onto the variables which are constrained.

The experimental results on random 3-satisfiability discussed in sect. IV.2 confirm the theoretical analysis of MEPAZE ; MZ\_pre  which indicate that all $`\eta _{ai}`$ vanish for $`\alpha <\alpha _d3.91`$. This can be interpreted as the fact that, in this range of $`\alpha `$, there are no ’constrained clusters’ (meaning clusters in which some of the variables are constrained). For $`\alpha _d<\alpha <\alpha _c=4.267`$ the theory predicts the existence of non-trivial messages, meaning that there exist constrained clusters. This is the region where SP and SID are able to outperform existing algorithms. One should notice that the SID does not seem to converge up to the conjectured threshold $`\alpha _c=4.267`$: although it is very small, there remains a region in the ’Hard SAT’ phase, close to $`\alpha _c`$ where the simple version of SP/SID does not give the result. Recent work shows that this gap can be reduced by using some backtracking strategy in SIDParisi\_backtrack . Whether it will be possible to close the gap with generalized versions of SP/SID, while keeping a typically polynomial running time, is an interesting open issue.

### V.2 The “don’t-care” state

In a given cluster $`\alpha `$, a variable $`x_i`$ can be thought of being in three possible states: either it is constrained equal to $`0`$ (this means that $`x_i=0`$ in all SAT assignments of the cluster $`\alpha `$), or it is constrained equal to $`1`$, or it is not constrained. In this last situation we attribute it the value $``$. Therefore we can describe a cluster by the values of the $`N`$ generalized variables $`x_i^\alpha \{0,1,\}`$, where $``$ will be denoted as the don’t-care state. Such a description associates to each cluster a single point in $`\{0,1,\}^N`$. It discards a lot of information: $`x_i^\alpha =`$ has lost all the information on the fraction of assignments in the cluster $`\alpha `$ where $`x_i=0`$. But it gives a simplified description of the cluster and it focuses onto the constrained variables.

It is interesting to notice that the SP equations can be interpreted as BP equations in the presence of this extra don’t-care state. This can be seen as follows. Borrowing the notations of the BP equations (8, 9) we denote by $`\gamma _{ia}[0,1]`$ the probability that the variable $`x_i`$ is in the state which violates clause $`a`$, in a problem where clause $`a`$ would be absent, and by

$$\delta _{ai}=\underset{jV(a)i}{}\gamma _{ja}$$
(33)
the probability that all variables in clause $`a`$, except variable $`i`$, are in the state which violates the clause.

Let us compute $`\mu _{ia}(x_i)`$. This depends on the messages sent from the nodes $`b\{V(i)a\}`$ to variable $`i`$. The various possibilities for these messages are:

* No warning arriving from $`bV_a^s(i)`$, and no warning arriving from $`bV_a^u(i)`$. This happens with a probability

$$\mathrm{\Pi }_{ia}^0=\underset{bV(i)a}{}\left(1\delta _{bi}\right)$$
(34)
* No warning arriving from $`bV_a^s(i)`$, and at least one warning arriving from $`bV_a^u(i)`$. This happens with a probability

$$\mathrm{\Pi }_{ia}^u=\left[1\underset{bV_a^u(i)}{}\left(1\delta _{bi}\right)\right]\underset{bV_a^s(i)}{}\left(1\delta _{bi}\right)$$
(35)
* No warning arriving from $`bV_a^u(i)`$, and at least one warning arriving from $`bV_a^s(i)`$. This happens with a probability

$$\mathrm{\Pi }_{ia}^s=\left[1\underset{bV_a^s(i)}{}\left(1\delta _{bi}\right)\right]\underset{bV_a^u(i)}{}\left(1\delta _{bi}\right)$$
(36)
* At least one warning arriving from $`bV_a^u(i)`$, and at least one warning arriving from $`bV_a^s(i)`$. This happens with a probability

$$\mathrm{\Pi }_{ia}^c=\left[1\underset{bV_a^s(i)}{}\left(1\delta _{bi}\right)\right]\left[1\underset{bV_a^u(i)}{}\left(1\delta _{bi}\right)\right]$$
(37)

As we work only with SAT configurations, the contradictory messages must be excluded. Therefore, the probability $`\gamma _{ia}[0,1]`$ that the variable $`x_i`$ is in the state which violates clause $`a`$, given that there is no contradiction, is:

$$\gamma _{ia}=\mathrm{\Pi }_{ia}^u/(\mathrm{\Pi }_{ia}^u+\mathrm{\Pi }_{ia}^s+\mathrm{\Pi }_{ia}^0)$$
(38)
The above equations in the enlarged space including the null message, given in (33-38) are identical to the SP equations(27,26), with the identification $`\eta _{ai}=\delta _{ai}`$.

### V.3 Complexity

The above interpretation of SP using the don’t-care state suggest a method to estimate the complexity, defined as the normalized logarithm of the number of constrained clusters of SAT assignments. As each constrained cluster of SAT assignments is associated with a point in $`\{0,1,\}^N`$, the complexity should be given by the corresponding entropy, which can be estimated with the usual BP formula (14) in the presence of the don’t-care state (see ref. BZ\_jstat  for a rigorous derivation). The result, originally derived in MZ\_pre  from a direct statistical physics analysis without the use of the don’t-care state, is as follows.

The total complexity $`\mathrm{\Sigma }`$ can be decomposed into contributions associated with every function node and with every variable, and reads:

$$\mathrm{\Sigma }=\underset{a=1}{\overset{M}{}}\mathrm{\Sigma }_a\underset{i=1}{\overset{N}{}}(n_i1)\mathrm{\Sigma }_i$$
(39)
where

$$\mathrm{\Sigma }_a=+\mathrm{log}\left[\underset{jV(a)}{}\left(\mathrm{\Pi }_{ja}^u+\mathrm{\Pi }_{ja}^s+\mathrm{\Pi }_{ja}^0\right)\underset{jV(a)}{}\mathrm{\Pi }_{ja}^u\right]$$
(40)
$$\mathrm{\Sigma }_i=+\mathrm{log}\left[\widehat{\mathrm{\Pi }}_i^++\widehat{\mathrm{\Pi }}_i^{}+\widehat{\mathrm{\Pi }}_i^0\right]$$
(41)
In Fig. 6 we report the data for the complexity of random 3-SAT formulae of size $`N=10^6`$ and $`\alpha `$ in the clustering range. This has been obtained with the following procedure. One generates first a ’random’ 3-SAT formula with $`N=10^6`$ and $`M=\mathrm{4.27\; 10}^6`$, using a pseudo-random number generator. The SP algorithm is run on this formula, and the complexity is evaluated from (39). Then a new formula is generated from the previous one by eliminating $`10^4`$ (pseudo-)randomly chosen clauses, and the algorithm is run again, etc… It turns out that, for $`N=10^6`$, the resulting curve is very ‘reproducible’, in the sense that the fluctuations of the curve from one random instance to the next are small (typically below $`1\%`$).

### V.4 Interpretation of the SID algorithm: categories of variables

Once SP has reached convergence, we can compute the total biases $`\{W_i^{(\pm )},W_i^{(0)}\}`$. According to the previous discussion, these numbers should be interpreted as giving the fraction of constrained clusters where the variable $`x_i`$ is respectively (frozen positive)/ (frozen negative)/ unconstrained. Having computed these weights, we may distinguish three reference types of variable nodes (of course all the intermediate cases will also be present): the under-constrained ones with $`W_i^{(0)}1`$, the biased ones with either $`W_i^{(+)}1`$ or $`W_i^{()}1`$, and the balanced ones with $`W_i^{(+)}W_i^{()}`$ and $`W_i^{(0)}`$ small.

Fixing a variable of each of these types produces different effects, consistently with the interpretation of the surveys. The following behaviours can be easily checked in numerical experiments:

Fixing a biased variable does not alter the structure of the clusters and the complexity changes smoothly (few constrained clusters are eliminated). This is the strategy used by the SID procedure.

Fixing an under-constrained variable has no effect on the complexity.

As expected, fixing a balanced variable produces a decrease very close to $`\mathrm{ln}2`$ in the complexity.

### V.5 A summary of the main conjectures

The whole interpretation relies on the existence, in certain ‘hard SAT’ regions of parameters (here in a window of $`\alpha `$ below $`\alpha _c`$), of clusters of SAT assignments, which are very far apart (one cannot reach a cluster from another one unless one flips a finite fraction of the $`N`$ variables). It would be very interesting to prove this statement.

When there exist such clusters, some of them may be ‘constrained clusters’, in which some variables are constrained to take the same value for all the assignments of the cluster. The introduction of the don’t-care state is an attempt at identifying all these clusters (in which each cluster is characterized by a single point in the enlarged assignment space including the don’t-care state). SP, interpreted as BP in this enlarged assignment space, is an attempt at obtaining a statistical description of these constrained clusters.

Two conjectures arise naturally from the heuristic statistical physics approach and the numerical experiments on SP. They should hold for the random satisfiability problem, in the large $`N`$ limit, in some hard SAT window of $`\alpha `$ just below $`\alpha _c`$ <sup>1</sup><sup>1</sup>1After the completion of this work, two groups have arrived to the conclusion that the ’simple’ clustering scenario described here should hold in the region $`\alpha [4.15,\alpha _c4.2667]`$ MMZ ; MPR :

* 1) With probability one (on the set of initial messages), SP converges to a unique set of fixed point messages.
* 2) These fixed point messages contain the correct information about the constrained clusters, and in particular the number of constrained clusters can be computed as in (39).

Note that for finite $`N`$, there are obvious counterexamples to these conjectures, as shown in fig. 7.

## VI Comments and perspectives

When the solutions of a SAT (or more generally of a constrained satisfaction problem) tend to cluster into well separated regions of the assignment space, it often becomes difficult to find a global solution because (at least in all local methods), different parts of the problem tend to adopt locally optimal configurations corresponding to different clusters, and these cannot be merged to find a global solution. In SP we proceed in two steps. First we define some elementary messages which are warnings, characteristic of each cluster, then we use as the main message the surveys of these warnings. The warning that we used is a rather simplified object: it states which variables are constrained and which are not constrained. Because of this simplification the surveys can be handled easily, which makes the algorithm rather fast. As we saw in sect.V, one might also aim at a finer description of each cluster, where the elementary messages would give the probability that a variable is in a given state ($`\{0,1\}`$), when considering the set of all SAT configurations inside this given cluster. In this case the survey will give the probability distribution of this probability distribution, when one chooses a cluster randomly. This is more cumbersome algorithmically, but there is no difficulty of principle in developing this extension. One could also think of generalizing further the messages (to probability distributions of probability distributions of probability distributions), if some problems have a structure of clusters within other clusters (this is known to happen for instance in spin glasses MPV ), but the cost in terms of computer resources necessary to implement it might become prohibitive.

We have presented here a description of a new algorithmic strategy to handle the SAT problem. This is a rather general strategy which can also be applied in principle, to all Constraint Satisfaction Problems coloring ; MZ\_pre . At the moment the approach is very heuristic, although it is based on a rather detailed conjecture concerning the phase structure of random 3-SAT. The validity of a similar conjecture has been checked exactly in the random XOR-SAT problem. It would be interesting to study this algorithmic strategy in its own, independently from the random 3-SAT problem and the statistical physics study. One can expect progress on SP to be made in the future in various directions, among which: Rigorous results on convergence, different ways of using the information contained in the surveys, generalization of the algorithm to deal with complex graphs which are not typical random graphs, use of the generalized SP with penalty to provide UNSAT certificates. The generalization of SP to generic Constraint Satisfaction Problems is discussed in a separate publication BMWZ .

## Acknowledgments

We are indebted to G. Parisi and M. Weigt for very fruitful discussions. We also thank C. Borgs, J. Chayes, B. Hayes, S. Kirkpatrick, S. Mertens, O. Martin, P. Purdom, M. Safra, B. Selman, J. Yedidia for stimulating exchanges. This work has been supported in part by the European Community’s Human Potential Programme under contract HPRN-CT-2002-00319, ’STIPCO’.


text:
# 1 Introduction

## 1 Introduction

Recently, Dijkgraaf and Vafa  have proposed a technique for calculating the effective superpotential for the glueball superfield in an $`𝒩=1`$ gauge theory through the planar diagrams of matrix model. Several tests have been performed for this proposal in -. In particular, the theories with matter fields in the fundamental representation have been found in .

In this letter, we compute the matrix integral over tree level superpotential characterized by light meson fields (mass deformation) in electric $`𝒩=1`$ $`SO(N_c)`$ gauge theory with $`N_f`$ flavors of quarks in the vector (fundamental) representation found by Intriligator-Seiberg  and we reproduce the exact effective superpotential in the gauge theory side. Moreover, we do similar analysis in magnetic theory. The matrix descriptions of both electric and magnetic theories are coincident with each other implying Seiberg duality in the gauge theory side. In , there are some relevant works on the gauge group $`SO(N_c)`$ in the view point of matrix model. For $`U(N_c)`$ gauge theory with $`N_f`$ flavors of quarks in the fundamental and anti-fundamental representations, the matrix descriptions of both electric and magnetic theories have been found in .

## 2 Matrix model description of electric theory

Let us deform our $`𝒩=1`$ supersymmetric $`SO(N_c)`$ gauge theory with $`N_f`$ flavors of quarks $`Q_a^j(j=1,2,\mathrm{},N_f,a=1,2,\mathrm{},N_c)`$ in the vector (fundamental) representation( $`N_c4,N_fN_c5`$) , by adding the mass terms of gauge invariant meson superfields $`M^{ij}=Q^iQ^j`$

$`W_{\text{tree}}={\displaystyle \frac{1}{2}}\text{Tr}mM={\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{N_f}{}}}m_j{\displaystyle \underset{a=1}{\overset{N_c}{}}}Q_a^jQ_a^j`$ (2.1)
where it is understood that quark superfields are represented in an $`N_f\times N_c`$ matrix form in terms of the flavor and color indices. In the gauge theory side, one can make this superfield matrix diagonal by using the gauge and global rotation.

The idea is to use this tree level superpotential as the potential for the matrix model. At first, let us consider $`SO(N)`$ matrix model at large $`N`$ by replacing the gauge theory fields with $`N_f\times N`$ matrices in order to calculate the contributions to the free energy. Then the partition function can be written as

$`Z={\displaystyle \frac{1}{\text{Vol}(SO(N))}}\left({\displaystyle \frac{\mathrm{\Lambda }}{2\pi g_s}}\right)^{\frac{1}{2}N_fN}{\displaystyle \underset{j=1}{\overset{N_f}{}}[dQ_a^j]e^{\frac{1}{2g_s}_{j=1}^{N_f}m_j_{a=1}^NQ_a^jQ_a^j}}`$
where we put the factor $`\left(\frac{\mathrm{\Lambda }}{2\pi g_s}\right)^{\frac{1}{2}N_fN}`$ in order to make the integration measure dimensionless111In fact we can put any dimensionful term with the right dimension inside the parenthesis, but this particular form gives the right form of the superpotential for the gauge theory. and $`\mathrm{\Lambda }`$ is the scale of $`SO(N_c)`$ gauge theory and $`m_j`$ is the mass of $`j`$-th quark. In this gauge theory there is no difference between left- and right-handed quarks which is different from $`SU(N_c)`$ gauge theory where the right-handed quark is the left-handed anti-quark. Recall that the superpotential should have $`R`$-charge 2 and dimension 3.

The large $`N`$-limit behavior of the volume for $`SO(N)`$ group, the normalization of the matrix path integral, can be read off from

$`\mathrm{log}\text{Vol}(SO(N))`$ $`=`$ $`\mathrm{log}{\displaystyle \frac{\sqrt{2}\left(2\pi \right)^{N^2/4}}{\left(N3\right)!\left(N5\right)!\mathrm{}3!1!\left(N/21\right)!}}`$
$`=`$ $`{\displaystyle \frac{1}{4}}N^2\mathrm{log}N+{\displaystyle \frac{1}{4}}N^2\left({\displaystyle \frac{3}{2}}+\mathrm{log}\pi +\mathrm{log}2\right)+{\displaystyle \frac{N}{4}}\mathrm{log}N`$
$`+{\displaystyle \frac{1}{24}}\mathrm{log}N+{\displaystyle \frac{N}{4}}\left(1+\mathrm{log}2\mathrm{log}\pi \right)+𝒪(1).`$
The usual matrix integral for flavor part is

$`{\displaystyle \underset{j=1}{\overset{N_f}{}}[dQ_a^j]e^{\frac{1}{2g_s}_{j=1}^{N_f}m_j_{a=1}^NQ_a^jQ_a^j}}=\left({\displaystyle \frac{\pi ^{N_f}g_s^{N_f}}{\text{det}m}}\right)^{N/2}.`$
In the large $`N`$-limit we are interested in, the glueball superfield $`S`$ can be identified with $`g_sN`$ as a second step and one can write the log of the partition function as follows:

$`\mathrm{log}Z`$ $`=`$ $`{\displaystyle \frac{1}{4}}N^2\mathrm{log}\left({\displaystyle \frac{N}{2\pi e^{3/2}}}\right)+{\displaystyle \frac{N}{2}}N_f\mathrm{log}\mathrm{\Lambda }{\displaystyle \frac{N}{2}}\mathrm{log}\text{det}m`$
$`=`$ $`{\displaystyle \frac{S^2}{4g_s^2}}\mathrm{log}\left({\displaystyle \frac{S}{2\pi e^{3/2}g_s}}\right)+{\displaystyle \frac{S}{2g_s}}N_f\mathrm{log}\mathrm{\Lambda }{\displaystyle \frac{S}{2g_s}}\mathrm{log}\text{det}m`$
$``$ $`{\displaystyle \frac{1}{g_s^2}}_2{\displaystyle \frac{1}{g_s}}_1.`$
Note that all the effects of the matter is in $`_1`$.

Then the effective superpotential for the glueball field $`S`$ can be computed as the derivative of the contribution to free energy plus the contribution from flavors

$`W`$ $`=`$ $`\left(N_c2\right){\displaystyle \frac{_2}{S}}+_1`$ (2.2)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)\left(SS\mathrm{log}\left({\displaystyle \frac{S}{\mathrm{\Lambda }^3}}\right){\displaystyle \frac{S}{N_c2}}N_f\mathrm{log}\mathrm{\Lambda }+{\displaystyle \frac{S}{N_c2}}\mathrm{log}\text{det}m\right)`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)\left(SS\mathrm{log}{\displaystyle \frac{S}{\left(\mathrm{\Lambda }^{3(N_c2)N_f}\text{det}m\right)^{1/\left(N_c2\right)}}}\right).`$
Here we have used the results that the derivative of $`_2`$ gives the Veneziano-Yankielowicz superpotential.

Solving the F-flatness condition $`_SW=0`$ (minimizing $`W`$ with respect to a glueball superfield $`S`$) one gets, by the scale normalization $`\mathrm{\Lambda }^{3\left(N_c2\right)N_f}=16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}`$,

$`S=\left(16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}\text{det}m\right)^{1/\left(N_c2\right)}ϵ_{N_c2},ϵ_{N_c2}=e^{2\pi ik/\left(N_c2\right)},k=1,\mathrm{},\left(N_c2\right)`$
with the phase factor $`e^{2\pi ik/\left(N_c2\right)},k=1,\mathrm{},\left(N_c2\right)`$ reflecting the $`(N_c2)`$ supersymmetric vacua but physically equivalent vacua of the theory coming from the spontaneous breaking of a discrete symmetry. Then the exact superpotential by plugging $`S`$ into (2.2)( integrating out the $`S`$) leads to

$`W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)ϵ_{N_c2}\left(16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}\text{det}m\right)^{1/\left(N_c2\right)}`$ (2.3)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)ϵ_{N_c2}\widehat{\mathrm{\Lambda }}_{N_cN_f,0}^3`$
where $`\widehat{\mathrm{\Lambda }}_{N_cN_f,0}^3=\left(16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}\text{det}m\right)^{1/\left(N_c2\right)}`$ is the low energy scale  for the $`SO(N_cN_f)`$ Yang-Mills theory. This is nothing but ADS superpotential below for pure $`SO(N_c)`$ Yang-Mills gauge theory. In the low energy effective theory, the classical vacuum degeneracy was lifted by quantum effects which is represented by a dynamically generated superpotential for the light meson fields $`M^{ij}`$. The superpotential  generated by gaugino condensation leads to

$`W_{ADS}={\displaystyle \frac{1}{2}}\left(N_cN_f2\right)ϵ_{N_cN_f2}\left({\displaystyle \frac{16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}}{\text{det}M}}\right)^{1/\left(N_cN_f2\right)}`$
where $`ϵ_{N_cN_f2}`$ is the $`(N_cN_f2)`$-th root of unity. This theory has no vacuum but by adding the mass terms (2.1) to the superpotential $`W_{ADS}`$, the theory has $`(N_c2)`$ supersymmetric vacua. If not all of the matter fields are massive, in the gauge theory side, one can integrate out massive quarks and get the effective superpotential at low energy for the massless ones. It is the same form as the above $`W_{ADS}`$ but with the scale replaced by the low energy one.

## 3 Matrix model description of magnetic theory

In the IR theory of electric theory in previous section, the magnetic theory is described by an $`SO(\stackrel{~}{N}_c=N_fN_c+4)`$ gauge theory ( $`N_f>N_c,N_c4`$ ) with $`N_f`$ flavors of dual quarks $`q_a^j(j=1,\mathrm{},N_f,a=1,\mathrm{},\stackrel{~}{N}_c)`$ and the additional gauge singlet fields $`M^{ij}`$ which is an elementary field of dimension 1 at the UV-fixed point . The matter field variables in the magnetic theory are the original electric variables $`M^{ij}`$ and magnetic quarks $`q_i`$ with the superpotential together with mass term

$`W={\displaystyle \frac{1}{2\mu }}M^{ij}q_iq_j+{\displaystyle \frac{1}{2}}\text{Tr}mM.`$
After integrating $`X`$ first, the partition function can be written as a delta function along the line of :

$`Z`$ $`=`$ $`{\displaystyle \frac{1}{\text{Vol}(SO(N))}}\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Lambda }}}{2\pi g_s}}\right)^{\frac{1}{2}N_f\stackrel{~}{N}}{\displaystyle [dX]\underset{j=1}{\overset{N_f}{}}[dq_j]\delta \left(X^{ij}M^{ij}\right)e^{\frac{1}{2g_s}\left(\text{Tr}mX+\frac{1}{\mu }_{i,j=1}^{N_f}X^{ij}q_iq_j\right)}}`$
$`=`$ $`{\displaystyle \frac{1}{\text{Vol}(SO(N))}}\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Lambda }}}{2\pi g_s}}\right)^{\frac{1}{2}N_f\stackrel{~}{N}}{\displaystyle \underset{j=1}{\overset{N_f}{}}[dq_j]\delta \left(\mu m+q_iq_j\right)}.`$
In the gauge theory side, this behavior of delta function is equivalent to the $`M^{ij}`$ equations of motion led by $`<q_iq_j>=\mu m`$. The constrained matrix integral over $`N_f`$ flavors of length $`\stackrel{~}{N}`$ can be obtained Wishart random matrices and the result  for this is

$`{\displaystyle \underset{j=1}{\overset{N_f}{}}[dq_j]\delta \left(\mu m+q_iq_j\right)}=c\times \left(\text{det}\left(\mu m\right)\right)^{\left(\stackrel{~}{N}N_f1\right)/2}`$
where the coefficient $`c`$ behaves like as

$`e^{\frac{1}{2}N_f\stackrel{~}{N}\mathrm{log}\stackrel{~}{\frac{N}{2}}}`$
for large $`\stackrel{~}{N}`$-limit.

As we did for electric case, the effective superpotential from the log of partition function can be expressed as

$`W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\stackrel{~}{N}_c2\right)\left(SS\mathrm{log}{\displaystyle \frac{S}{\stackrel{~}{\mathrm{\Lambda }}^3}}\right){\displaystyle \frac{1}{2}}N_f\left(SS\mathrm{log}{\displaystyle \frac{S}{\stackrel{~}{\mathrm{\Lambda }}^3}}\right)+SN_f\mathrm{log}\stackrel{~}{\mathrm{\Lambda }}{\displaystyle \frac{1}{2}}S\mathrm{log}\text{det}\left(\mu m\right)`$ (3.1)
$`=`$ $`{\displaystyle \frac{1}{2}}(\stackrel{~}{N}_cN_f2)(SS\mathrm{log}S+S\mathrm{log}\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Lambda }}^{3\left(\stackrel{~}{N}_c2\right)N_f}}{\text{det}\left(\mu m\right)}}\right)^{1/\left(\stackrel{~}{N}_cN_f2\right)}).`$
Solving the F-flatness condition $`_SW=0`$ one gets, with the phase factor $`e^{2\pi ik/\left(\stackrel{~}{N}_cN_f2\right)},k=1,\mathrm{},\left(\stackrel{~}{N}_cN_f2\right)`$ reflecting the $`(\stackrel{~}{N}_cN_f2)`$ supersymmetric vacua,

$`S=\left({\displaystyle \frac{16\stackrel{~}{\mathrm{\Lambda }}_{N_fN_c+4,N_f}^{3\left(\stackrel{~}{N}_c2\right)N_f}}{\text{det}\left(\mu m\right)}}\right)^{1/\left(\stackrel{~}{N}_cN_f2\right)}ϵ_{\stackrel{~}{N}_cN_f2},\stackrel{~}{\mathrm{\Lambda }}^{3\left(\stackrel{~}{N}_c2\right)N_f}=16\stackrel{~}{\mathrm{\Lambda }}_{N_fN_c+4,N_f}^{3\left(\stackrel{~}{N}_c2\right)N_f}.`$
Then the exact superpotential by plugging this $`S`$ into (3.1) leads to

$`W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\stackrel{~}{N}_cN_f2\right)ϵ_{\stackrel{~}{N}_cN_f2}\left({\displaystyle \frac{16\stackrel{~}{\mathrm{\Lambda }}_{N_fN_c+4,N_f}^{3\left(\stackrel{~}{N}_c2\right)N_f}}{\text{det}\left(\mu m\right)}}\right)^{1/\left(\stackrel{~}{N}_cN_f2\right)}`$ (3.2)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c+2\right)ϵ_{N_c2}\left({\displaystyle \frac{16\stackrel{~}{\mathrm{\Lambda }}_{N_fN_c+4,N_f}^{3\left(N_fN_c+2\right)N_f}}{\left(1\right)^{N_f}\mu ^{N_f}\text{det}m}}\right)^{1/\left(N_c2\right)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)ϵ_{N_c2}\left(16\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}\text{det}m\right)^{1/\left(N_c2\right)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(N_c2\right)ϵ_{N_c2}\widehat{\mathrm{\Lambda }}_{N_cN_f,0}^3`$
which is exactly the same as the one in (2.3). Here we used the fact that $`\stackrel{~}{N}_c=N_fN_c+4`$, $`\text{det}\left(\mu m\right)=\left(1\right)^{N_f}\mu ^{N_f}\text{det}m`$ and $`ϵ_{\stackrel{~}{N}_cN_f2}=ϵ_{N_c+2}=ϵ_{N_c2}`$. The scale of the magnetic theory in the gauge theory side was related to that of electric theory by

$`2^8\mathrm{\Lambda }_{N_c,N_f}^{3\left(N_c2\right)N_f}\stackrel{~}{\mathrm{\Lambda }}_{N_fN_c+4,N_f}^{3\left(N_fN_c+2\right)N_f}=\left(1\right)^{N_fN_c}\mu ^{N_f}`$
where the normalization factor $`1/2^8`$ was chosen to get the consistent low energy behavior under large mass deformation and along the flat directions. Note that the factor $`(1)^{N_c/(N_c2)}=1`$ in (3.2) is cancelled exactly by the overall $`1`$ factor. If not all of the matter fields are massive, the remaining low energy magnetic theory is a magnetic $`SO(N_fN_c+4M)`$ gauge theory with $`(N_fM)`$ flavors and the low energy superpotential. This low energy magnetic theory is dual to the low energy $`SO(N_c)`$ gauge theory with $`(N_fM)`$ massless quarks. One can arrive at the scale relation which connects the $`SO(N_c)`$ electric theory with $`(N_fM)`$ flavors to the $`SO(N_fN_c+4M)`$ magnetic theory with $`(N_fM)`$ flavors.

## 4 Discussions

In this paper, we apply the matrix model for fundamental flavors without any adjoint matter to check the Seiberg duality when we deform mass terms. By explict matrix path integral in both electric and magnetic theories, we have found the complete agreement with the field theory result.

One can consider and generalize degenerate mass deformation in which some eigenvalues of the mass matrix $`m`$ are vanishing. We expect to have the following superpotential, effectively $`SO(N_c)`$ theory with $`K`$ flavors,

$`W={\displaystyle \frac{1}{2}}\left(N_cK2\right)ϵ_{N_cK2}\left({\displaystyle \frac{16\mathrm{\Lambda }_{N_c,K}^{3\left(N_c2\right)K}\text{det}m}{\text{det}M}}\right)^{1/\left(N_cK2\right)}`$
where $`M^{ij}(i,j=1,2,\mathrm{},K)`$ are the meson from the massless flavors and the index $`i`$ of mass matrix eigenvalues $`m_i`$ runs $`i=K+1,K+2,\mathrm{},N_f`$. We have to add the terms from the delta function constraint for massless flavor part. For $`K=0`$, we can reproduce the result of section 2 and 3.

Acknowledgments

This research of CA was supported by Korea Research Foundation Grant(KRF-2002-015-CS0006). This research of SN was supported by Korea Research Foundation Grant KRF-2001-041-D00049.


text:
# Combining first KamLAND results with solar neutrino data

## I Introduction

In a recent paper the first results of the KamLAND collaboration became public kamlandPRL . These data contain precious information on the neutrino oscillation hypothesis which has been advocated to account for a number of neutrino experiments involving solar and atmospheric neutrinos and which indicate that neutrinos are massive and that neutrino flavor mixing occurs Fukuda:2001nj ; Fukuda:2002pe ; Fukuda:1998fd ; Fukuda:1998rq ; Cleveland:nv ; Davis:jw ; Abdurashitov:1999zd ; Hampel:1998xg ; Altmann:2000ft ; Cattadori:rd ; sno ; Fukuda:1998mi . The KamLAND experiment is a reactor neutrino experiment with its detector located at the Kamiokande site. Most of the $`\overline{\nu }_e`$ flux incident at KamLAND comes from plants at distances of $`80350`$ km from the detector, making the average baseline of about 180 kilometers, long enough to provide a sensitive probe of the LMA solution of the solar neutrino problem MSW ; Bahcall:1999ed ; Gonzalez-Garcia:1999aj . The KamLAND collaboration has for the first time measured the disappearance of neutrinos traveling to a detector from a power reactor. They observe a strong evidence for the disappearance of neutrinos during their flight over such distances, giving the first terrestrial confirmation of the solar neutrino anomaly and also establishing the oscillation hypothesis with man-produced neutrinos. Moreover the parameters that describe this disappearance in terms of the oscillations of the electron neutrino type to another, are consistent with latest pre-KamLAND determinations Maltoni:2002ni ; Bahcall:2002hv ; Bandyopadhyay:2002xj ; Barger:2002iv ; deHolanda:2002pp ; Creminelli:2001ij ; Fogli:2002pt  of solar neutrino parameters.

In this note we analyze the implications of these fundamental results by combining the KamLAND data with data from solar neutrino experiments. We will assume CPT conservation and for simplicity we consider a two-flavor massive neutrino oscillation framework. In Sec. II we analyze the impact of the KamLAND results by including the full information on the spectral distribution of the observed events. Subsequently, in Sec. III we perform a global fit that combines the full KamLAND and Chooz reactor data sample Apollonio:1999ae  with the full solar neutrino data as included in Ref. Maltoni:2002ni . In Sec. IV we check the stability of the results with respect to changes in the statistical analysis, and we summarize in Sec. V.

## II Simulation and analysis of KamLAND data

In KamLAND the target for the $`\overline{\nu }_e`$ flux consists of a spherical transparent balloon filled with 1000 tons of non-doped liquid scintillator. The anti-neutrinos are detected via the inverse neutron $`\beta `$-decay

$$\overline{\nu }_e+pe^++n.$$
(1)
In Fig. 5 of Ref. kamlandPRL  the spectral data are given in 13 bins of prompt energy above 2.6 MeV. To simulate the KamLAND data we calculate the expected number of events in each bin for given oscillation parameters as

$$\begin{array}{cc}\hfill N_i^{\mathrm{th}}(\mathrm{\Delta }m^2,\theta )& =fdE_\nu \sigma (E_\nu )\times \hfill \\ \hfill \underset{j}{}\varphi _j(E_\nu )& P_j(E_\nu ,\mathrm{\Delta }m^2,\theta )_i𝑑E_eR(E_e,E_e^{}).\hfill \end{array}$$
(2)
Here $`R(E_e,E_e^{})`$ is the energy resolution function and $`E_e,E_e^{}`$ are the observed and the true positron energy, respectively, and we use an energy resolution of $`7.5\%/\sqrt{E(\mathrm{MeV})}`$ kamlandPRL . The neutrino energy is related to the positron energy by $`E_\nu =E_e^{}+\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the neutron-proton mass difference. The integration interval over $`E_e`$ is determined by the prompt energy interval in each bin. The neutrino spectrum $`\varphi (E_\nu )`$ from nuclear reactors is well known, we are using the phenomenological parameterization given in Refs. Vogel:iv ; Murayama:2000iq . We adopt the average fuel composition for the nuclear reactors given in Ref. kamlandPRL . Note that possible effects due to time variations in the fuel composition have been shown to be small Murayama:2000iq . The sum over $`j`$ in Eq. (2) runs over 16 nuclear plants, taking into account the different distances from the detector and the power output of each reactor (see Table 3 of Ref. kamlandproposal ). The relevant detection cross section $`\sigma (E_\nu )`$ is given in Ref. Vogel:1999zy . In the 2-neutrino framework the survival probability for the neutrinos coming from the reactor $`j`$ is given by

$$P_j(E_\nu ,\mathrm{\Delta }m^2,\theta )=1\mathrm{sin}^22\theta \mathrm{sin}^2\frac{\mathrm{\Delta }m^2L_j}{4E_\nu }.$$
(3)
The normalization factor $`f`$ in Eq. (2) is determined in such a way that for the case of no oscillations we obtain a total number of events of 86.8, as expected from the Monte-Carlo simulation used in Ref. kamlandPRL .

For the statistical analysis we use the $`\chi ^2`$-function

$$\chi ^2=\underset{i,j}{}(N_i^{\mathrm{th}}N_i^{\mathrm{obs}})S_{ij}^1(N_j^{\mathrm{th}}N_j^{\mathrm{obs}}).$$
(4)
The observed number of events $`N_j^{\mathrm{obs}}`$ in each bin can be read off from Fig. 5 of Ref. kamlandPRL . In the covariance matrix $`S`$ we include the experimental error in each bin $`\sigma _i`$ (obtained from the same figure), which we assume to be uncorrelated, and the systematic error $`\sigma _{\mathrm{syst}}=0.0642`$ kamlandPRL  implied by the uncertainty on the total number of events expected for no oscillations:

$$S_{ij}=\sigma _i^2\delta _{ij}+\sigma _{\mathrm{syst}}^2N_i^{\mathrm{th}}N_j^{\mathrm{th}}.$$
(5)
This $`\chi ^2`$ definition assumes Gaussian distribution of the data. For the discussion of an alternative analysis based on poisson distributed data see Sec. IV.

## III Results and Discussion

Our results are summarized in Figs. 1, 2 and 3. In Fig. 1 we show the allowed regions of the oscillation parameters obtained from our re-analysis of the KamLAND data. It is in good agreement with the analysis performed by the KamLAND group, shown in Fig. 6 of Ref. kamlandPRL . This gives us confidence on our simulation of the KamLAND data and therefore encourages us to use it in a full analysis combining also with the solar data sample. Figs. 2 and 3 show the corresponding results obtained in a combined fit of the full KamLAND data sample with the global sample of solar neutrino data, as well as the Chooz result. The solar data we are using and the details of our solar neutrino analysis are given in Ref. Maltoni:2002ni .

First of all, we have quantified the rejection of non-LMA solutions and found that it is now more robust. For example, for the LOW solution we have $`\mathrm{\Delta }\chi _{\text{LOW–LMA}}^2=26.9`$, which for 2 d.o.f. ($`\mathrm{\Delta }m^2`$ and $`\theta `$) corresponds to a relative probability of $`1.4\times 10^6`$, assuming Gaussian errors. A similar result is also found for the VAC solution. Apart from selecting out LMA as the unique solution of the solar neutrino problem we find, however, that the new reactor results have little impact on the location of the best fit point:

$$\mathrm{tan}^2\theta =0.46,\mathrm{\Delta }m^2=6.9\times 10^5\text{eV}^2.$$
(6)
In particular the solar neutrino mixing remains significantly non-maximal, a point which is rather important for model-building. Indeed bi-maximal mixing models are disfavored Chankowski:2000fp  while models where the solar mixing can be non-maximal Babu:2002dz  are preferred, as before. This is not in contradiction with the fact that KamLAND data alone prefer maximal mixing kamlandPRL , since such preference has no statistical significance. Indeed, one can see from the right panel in Fig. 3 that $`\mathrm{\Delta }\chi ^2`$ is rather flat with respect to the mixing angle for $`\mathrm{tan}^2\theta 0.4`$. This explains why the addition of the KamLAND data has no impact whatsoever in the determination of the solar neutrino oscillation mixing. The allowed $`3\sigma `$ region we find for $`\theta `$ is:

$$0.29\mathrm{tan}^2\theta 0.86,$$
(7)
practically identical to the pre-KamLAND range given in Eq. (4) of Ref. Maltoni:2002ni .

On the other hand, the new data do have a strong impact in narrowing down the allowed range of $`\mathrm{\Delta }m^2`$. From the left panel of Fig. 3 one can read off that KamLAND data alone provides the bound $`\mathrm{\Delta }m^2>8\times 10^6\text{eV}^2`$, whereas the CHOOZ experiment gives $`\mathrm{\Delta }m^2<10^3\text{eV}^2`$, both at $`3\sigma `$. Hence, global reactor neutrino data provide a robust allowed interval for $`\mathrm{\Delta }m^2`$, based only on terrestrial experiments using artificial neutrino sources. However, combining this information from reactors with the solar neutrino data leads to a significant reduction of the allowed range: As clearly visible in Fig. 2, the original LMA region is now split into two sub-regions. From Fig. 3 we obtain at $`3\sigma `$ (1 d.o.f.)

$$5.1\times 10^5\text{eV}^2\mathrm{\Delta }m^29.7\times 10^5\text{eV}^2,$$
(8)
$$1.2\times 10^4\text{eV}^2\mathrm{\Delta }m^21.9\times 10^4\text{eV}^2.$$
(9)
The local minimum in the region (9) occurs for

$$\mathrm{tan}^2\theta =0.42,\mathrm{\Delta }m^2=1.4\times 10^4\text{eV}^2$$
(10)
with a $`\mathrm{\Delta }\chi ^2=4.5`$ with respect to the best fit point given in Eq. (6). This ambiguity might be resolved when more KamLAND data have been collected (see e.g. Refs. Murayama:2000iq ; deGouvea:2001su ; Barger:2000hy ).

## IV Stability of the statistical analysis

The current KamLAND data sample consists of 54 anti-neutrino events, which are distributed over the 13 energy bins. This leads to rather small numbers of events in each bin. The 5 bins with highest energies contain no event at all. In such a case the use of a $`\chi ^2`$-function based on Poisson statistics might be appropriate. In order to check the stability of our results we have performed also an analysis by using pdg

$$\chi ^2=2\underset{i}{}\left[\alpha N_i^{\mathrm{th}}N_i^{\mathrm{obs}}+N_i^{\mathrm{obs}}\mathrm{ln}\left(\frac{N_i^{\mathrm{obs}}}{\alpha N_i^{\mathrm{th}}}\right)\right]+\left(\frac{1\alpha }{\sigma _{\mathrm{syst}}}\right)^2$$
(11)
where the term containing the logarithm is absent in bins with no events. We minimize with respect to $`\alpha `$ in order to take into account the overall uncertainty of the theoretical predictions.

The analysis of KamLAND data using Eq. (11) is shown in Fig. 1 as the hollow lines. We observe that this analysis is somewhat less constraining compared to the analysis based on the Gaussian $`\chi ^2`$ of Eq. (4). One notices that smaller values of the mixing angle are allowed, especially at high convidence level. Let us note, however, that the allowed regions from the Gaussian analysis are in better agreement with the analysis done by the KamLAND group. This is the reason why we prefer to use this method for analyzing KamLAND data. The better agreement with the original KamLAND analysis might be related to the fact that the inclusion of the information on the experimental errors provided by the KamLAND collaboration in Fig. 5 of Ref. kamlandPRL  can only be included by means of a Gaussian $`\chi ^2`$-function, as in Eq. (4). In this way it is possible to take into account the asymmetric errors and the error bars in bins where the number of events is zero.

However, we note that the determination of $`\mathrm{\Delta }m^2`$ is rather stable, only the constraint on the mixing angle is somehow affected. Since the bound on the mixing angle in the combined analysis is dominated by solar data, we expect the difference between the two methods to be small after combining KamLAND with solar data. The results of this exercise are shown in Fig. 4. Comparing this figure with Fig. 2 we find indeed, that the result is very similar. The location of the best fit point and the 90%, 95% and 99% C.L. regions around the best fit point are nearly identical. However, the local minimum does not appear at the 90% C.L., though its location is, again, very stable. Some small differences are visible for the 99.73% C.L. contour.

To summarize, although there are some notable differences between the allowed regions obtained by assuming Gaussian or Poisson $`\chi ^2`$-functions for the KamLAND data taken alone, the differences are very small when combined with solar data. This illustrates the robustness of our results against variations in the statistical analysis.

## V Conclusions

We find that among all previous oscillation solutions to the solar neutrino anomaly, the new reactor results from the KamLAND experiment single out the LMA solution, rejecting all other oscillation solutions at a significant level. Furthermore, we find that already these first 145.1 days of KamLAND data lead to a significant improvement in narrowing down the allowed range of $`\mathrm{\Delta }m^2`$ when combined with solar neutrino data. The original LMA region now is split into two relatively narrow islands around the values of $`\mathrm{\Delta }m^2=6.9\times 10^5`$ eV2 (best fit point) and $`\mathrm{\Delta }m^2=1.4\times 10^4`$ eV2 (local minimum). However, our full analysis indicates that the new data have little impact on the determination of the mixing angle. In particular the solar neutrino mixing remains significantly non-maximal ($`3\sigma `$).

Before closing, let us note that we have considered here only the simplest case of two neutrinos. Analyzing in detail the impact of the KamLAND results on three-neutrino oscillation scenarios Gonzalez-Garcia:2000sq  and the resulting constraints is beyond the scope of this short note. The improved determination of $`\mathrm{\Delta }m^2`$ can also play an interesting role in probing fine details of solar physics Burgess:2002we , matter effects Fogli:2002hb , probing electro-magnetic neutrino properties Grimus:2002vb  (see also em\_others ) or testing CPT invariance in the neutrino sector Bahcall:2002ia . Similarly, the nailing down of LMA as the solution has also implications for non-oscillation solutions to the neutrino anomaly Valle:2002tm  in terms of spin-flavor precession Barranco:2002te ; Friedland:2002pg , non-standard interactions Guzzo:2001mi  or neutrino decay decay . Clearly none can now be the leading explanation to the solar neutrino anomaly Pakvasa:2003zv , although a detailed evaluation must be performed to decide, in each case, to what extent these solutions are now rejected.

###### Acknowledgements.

This work was supported by Spanish grant BFM2002-00345, by the European Commission RTN grant HPRN-CT-2000-00148, by the ESF *Neutrino Astrophysics Network* and by the *Sonderforschungsbereich 375-95 für Astro-Teilchenphysik* der Deutschen Forschungsgemeinschaft (T.S.). M.M. is supported by the Marie Curie contract HPMF-CT-2000-01008.