# Numerical Analysis of Constraint Satisfactory Problems using Belief Propagation and similar methods

Arthur Hardiagon

► To cite this version:

Arthur Hardiagon. Numerical Analysis of Constraint Satisfactory Problems using Belief Propagation and similar methods. Data Analysis, Statistics and Probability [physics.data-an]. 2016. dumas-05384362

HAL Id: dumas-05384362

https://dumas.ccsd.cnrs.fr/dumas-05384362v1

Submitted on 26 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical Analysis of Constraint Satisfactory
Problems using Belief Propagation and similar
methods

Arthur Hardiagon

June 30, 2016

Department of Physics, University of Rome La Sapienza

Supervisor : Federico Ricci-Tersenghi

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>
# Contents

<table><tbody><tr><td><strong>1 Introduction</strong></td><td><strong>3</strong></td></tr><tr><td>1.1 The Class of Constraint Satisfaction Problems</td><td>3</td></tr><tr><td>1.2 Phase Transition in CSPs</td><td>3</td></tr><tr><td>1.3 ARPE Project</td><td>4</td></tr><tr><td><strong>2 Graphical Models and Statistical Mechanics</strong></td><td><strong>4</strong></td></tr><tr><td>2.1 Definitions</td><td>5</td></tr><tr><td>2.1.1 Factor graph</td><td>5</td></tr><tr><td>2.1.2 The example of NAE-k-SAT</td><td>6</td></tr><tr><td>2.2 Bethe approximation and Belief Propagation</td><td>6</td></tr><tr><td>2.2.1 General statements</td><td>6</td></tr><tr><td>2.2.2 Case of NAE-SAT</td><td>8</td></tr><tr><td>2.3 The Cavity Method and Replica Symmetry Breaking (1RSB)</td><td>9</td></tr><tr><td>2.3.1 The Replica Symmetry Cavity Method</td><td>9</td></tr><tr><td>2.3.2 1RSB Cavity Method</td><td>11</td></tr><tr><td><strong>3 BP-inspired Search Algorithms</strong></td><td><strong>16</strong></td></tr><tr><td>3.1 Random Decimation</td><td>17</td></tr><tr><td>3.1.1 Description</td><td>17</td></tr><tr><td>3.1.2 Probability of convergence</td><td>19</td></tr><tr><td>3.1.3 Convergence of BP and precursor</td><td>20</td></tr><tr><td>3.2 Guided Decimation Algorithm</td><td>21</td></tr><tr><td>3.3 Biased Decimation Algorithm</td><td>23</td></tr><tr><td><strong>4 Conclusion</strong></td><td><strong>25</strong></td></tr><tr><td>4.1 Perspectives</td><td>25</td></tr><tr><td>4.2 Acknowledgments</td><td>26</td></tr></tbody></table>
# 1 Introduction

## 1.1 The Class of Constraint Satisfaction Problems

Constraint satisfaction problems (CSPs) is a class of problems that offers a ideal framework for testing new types of algorithms in hard computation. In different disciplines, they provide statistical tools to treat efficiently information of systems with a large number of components. An instance of a CSP is said to be satisfiable if there is an assignment of $N$ variables which simultaneously satisfy $M$ constraints. The interesting ratios $\alpha = N/M$ are close to the satisfiability threshold $\alpha_s$ above which it does not exist any assignments of variables that satisfy all constraints. They are in the core of the recent advances in information theory and artificial intelligence when problems need to be numerically solved in a reasonable time. In this report, we will focus on the Not-All-Equal-Satisfiability problem (NAE-SAT) that shares some features with other models like KSAT and XORSAT problems which have been widely studied in the mathematical, computational and physics literature. All of these problems have been shown to be NP complete for some values of the parameters [6, 17, 9]. It means that for this class of problems, only algorithms with exponential complexity are conjectured to solve them, which make them practically useless. However, some heuristics have been developed over the last few decades by physicists to design algorithms that find some solutions in a linear time [14].

## 1.2 Phase Transition in CSPs

Many of the features of the space of solutions of NAE-SAT are shared with other CSPs. We will give a brief overview of the successive transitions that have been studied in XORSAT [13], q-coloring (q-COL) and KSAT [11]. To illustrate our purpose, we reported a figure of the last article which summarizes the global situation in Fig.1. We shall now precise the key point of this paragraph by defining the term of *cluster*. Two solutions belongs to two distinct clusters if they are separated by a distance $O(N)$. For low values of the density $\alpha < \alpha_{d,+}$, the space of solution is entirely connected, all solutions belong to the same cluster. For $\alpha_{d,+} < \alpha < \alpha_d$, some clusters appear, some solutions are then separated by infinite barriers in the thermodynamic limit. Above the dynamical transition, for $\alpha_d < \alpha < \alpha_c$, an exponentially number of clusters are present, each of them contening an exponential number of solutions. Above the condensation transition, for $\alpha_c < \alpha < \alpha_s$, only a few clusters are dominant, in the sense that they contain most solutions. Finally, above the satifiability transition, for $\alpha > \alpha_s$, no solution remain.

In statistical mechanics, clusters are related to ground states of the system. The apparition of separated clusters also leads to local minima in the energy density. At the computing level, configurations retained by search algorithm are likely to be trapped in these metastable states. It explains why some algorithms like Simulated
<figure><img src="image_4.png" /><figcaption>Figure 1: Schematic Illustration of the space of solutions of a typical CSP. Each black circle stands for a cluster of solution.</figcaption></figure>

Annealing (SA) based on local Monte Carlo processes fail to find some solution up to the dynamical threshold.

An other approach consists in using message passing algorithm as *Belief Propagation* (BP) or *Survey Propagation* (SP) [3]. Both exhibit features of the underlying statistical landscape of the space of solutions. This study is largely dedicated to BP methods, although SP, which is slightly much difficult to implement has been showed to be efficient up to the condensation transition and beyond [3]. Other recent studies have focused on the analysis of the frozen clusters in NAE-SAT which provide additional features of the different phases for $\alpha_c < \alpha < \alpha_s$ [4].

## 1.3 ARPE Project

During these ten months, our first aim was to apply existing methods from statistical physics to have a global understanding of random NAE-SAT instances close to the satisfiability threshold [15]. We used population dynamic algorithm in the framework of the so called Symmetry Breaking (RS) and First Replica Symmetry Breaking (1RSB) methods in order to determine numerical values of the thresholds. We investigated the efficiency of different algorithms to search for solutions, in particular we displayed a successful BP algorithm that cross the dynamical transition. Additionally, we presented a new algorithm which aim to focus on a specific region of the space of solutions, discarding less probable configurations at early steps. We first checked the feasibility of the method, by determining the new phase diagram. It allowed us to design a search algorithm that solves random instances of the NAE-SAT problem beyond the dynamical transition. In order to complete this numerical study, we are likely to investigate analytically the decimation process used in our search algorithm as it has been done in [19] by including the new measure.

# 2 Graphical Models and Statistical Mechanics

In this first part, we will present some general features about CSPs and their mathematical formalization. We aimed to introduce in a synthetic manner the concepts of
<figure><img src="image_5.png" /><figcaption>Figure 2: Left panel : a factor graph representation of an hypergraph. Right panel: a portion of the factor graph of an instance of the NAE-3-SAT Problem, with Poissonian degree for variables nodes.</figcaption></figure>

the field, from graphical representations to statistical methods, mainly extracted from
spin-glass physics. For clarity we will focus on a specific CSPs, the Bicoloring Prob-
lem, that presents some similar characteristics to other famous problems as K-SAT or
XORSAT. We particularly refer to [12] that presents a overview of the state of the art.

2.1 Definitions

2.1.1 Factor graph

A Contraint Satisfaction Problem (CSP) is defined on a set of *N* variables σ₁, ..., σ*N*
that take values in a finite ensemble, the basic case being combinatory variables or spins
in Physics : σᵢ ∈ {-1, +1}. We denote σ = (σ₁, ..., σ*N*) a global configuration of the
variables and σ_S = {σᵢ, i ∈ S} a partial configuration of the variables in a given subset
S. Finding a solution of the problem corresponds to display a global configuration that
satisfy simultaneously *M* constraints (or clauses), each of them being specified by a
function ψₐ(σ∂aN) , where ∂ₐ denotes the subset of variables that enter in the clause a.
The function ψ can take values in {0, 1}, the value 1 (or 0) denoting that the clause is
satisfied (or unsatisfied). The choice of the function ψ shall define the specific CSP.

We call *instance* a set of *M* functions ψₐ, a ∈ {1, ..., M}. Then, for a given CSP, we
can generate differents instances. For this purpose, it is simpler to represent the CSP
in a factor graph as in Fig.2, which is a bipartite graph where empty squares (function
nodes) encode the hyperedges a and empty circles (variables nodes) denote the vertices
i; an edge being drawn if and only if a ∈ ∂i (or equivalent i ∈ ∂a). First, we impose a
number N of variables and a number of constraints M, playing on the main parameter
α = M/N, the density of clauses. Given both variables nodes and function nodes, we
obtain an instance of a problem by linking each function node with a given number of
variables nodes.
## 2.1.2 The example of NAE-k-SAT

In the following, we will present our work on a specific problem involving binary variables where the constraint function is defined as:

$$ \psi_a(\underline{\sigma}_{\partial a}) = \begin{cases} 0 & \text{if } \underline{\sigma}_{\partial a} = (+1, \dots, +1) \text{ or } (-1, \dots, -1) \\ 1 & \text{else} \end{cases} $$

In other words, a clause is satisfied only if at least two variables that enter in the clause have different values. We call this problem Not-All-Equal-Satisfiability Problem (NAE-SAT), which is also refered as the Bicoloring Hypergraph Problem in Computer Science, since it corresponds to coloring nodes of a graph avoiding any monocromatic hyperedge. The underlying structure of the graph shall be defined too. We distinguish the k-regular model and the Poisson-distributed model. The first case that corresponds to Erdös-Renyi random graphs consists in arranging k variables nodes for each function node and simultaneously setting at $L = \alpha k$ the number of function nodes attached to each variable nodes. This particular model has been rigorously studied by [8] [1]. For simplicity of construction, we have focused on the Poisson case, where we still have k variables entering each clause, and the distribution of the number of neighbors of variable nodes follow a Poissoninan distribution of intensity $\alpha k$, see Fig.2.

## 2.2 Bethe approximation and Belief Propagation

### 2.2.1 General statements

We can define an energy for the system as the total number of unsatisfied clauses :

$$ E = (\{\sigma_i, i \in 1, \dots, N\}) = \psi_a(\underline{\sigma}_{\partial a}) \quad (1) $$

Thus, finding a solution for this kind of problem is equivalent to a quest of 0-energy. It is partially why statistical methods coming from statistical mechanics have been largely employed to study analytically CSPs. First, we define a joint probability function on the configurational space, that gives the way we sample the system. If we want to consider all solutions as equivalent ones, we define the uniform measure:

$$ \mu(\underline{\sigma}) = \frac{1}{Z} \prod_{a=1}^{M} \psi_a(\underline{\sigma}_{\partial a}) \quad (2) $$

where Z is the partition function of the system and is also the total number of solution of the problem. As physicists, we would like to compute the free-entropy of the system defined by $\phi[\mu] = \log(Z)$. We also are interested in computing marginals $\mu_i(\sigma_i)$ of the probability distribution in eqn (2). In the general case, it is not possible to get an analytical form for $\mu_i(\sigma_i)$, but we can use the Bethe approximation to give an estimate
value. This approximation suggests to compute local quantities using the locally tree-like structure of random graphs. It seems reasonable since the typical length of loops in random graph scales as $log(N)$, which diverges at the thermodynamic limit.

Using this approximation, marginals can be expressed according to messages (or beliefs) that carry information from variable nodes to function nodes (and vice versa) [12]. Regarding the factor graph of the CSP shown in Fig.3, one can also suggest a probabilistic interpretation of these messages. In a bipartite graph, given an edge (ia) one shall consider two kind of messages $\mu_{i \to a}$ and $\hat{\mu}_{a \to i}$, which are probability distributions on the single-variable $\sigma_i \in \{-1, +1\}$. Indeed, we have, for instance, $\mu_{i \to a}(\sigma_i) > 0$ and $\sum_{\sigma_i} \mu_{i \to a}(\sigma_i) = 1$.

<figure><img src="image_6.png" /><figcaption>Figure 3: Illustration of Belief Propagation in a factor graph</figcaption></figure>

The message $\mu_{i \to a}(\sigma_i)$ is seen as the probability that assigning the value $\sigma_i$ to the variable $i$ would satisfy the clause $a$ in the absence of the edge (ia) and therefore when the branch attached to clause $a$ is not considered. In the other hand, the message $\hat{\mu}_{b \to i}(\sigma_i)$ can be seen as the probability that assigning the value $\sigma_i$ to the variable $i$ would satisfy the clause $b$ knowing that it is not satisfied by the other variables in the absence of the edge (ib). Using independence of branches in a tree-model, one can write equations of evolution of these messages, that we will refer as *Belief Propagation* equations, in the following way:

$$ \mu_{i \to a}^{(t+1)}(\sigma_i) \propto \prod_{b \in \partial j \setminus a} \hat{\mu}_{b \to i}^{(t)}(\sigma_i) \quad (3) $$

$$ \hat{\mu}_{a \to i}^{(t+1)}(\sigma_i) \propto \sum_{\sigma_{\partial a \setminus i}} \prod_{j \in \partial a \setminus i} \mu_{j \to a}^{(t)}(\sigma_j) \psi_a(\sigma_{\partial a}) \quad (4) $$
where exact values have to be calculated according to normalization. In a few paragraphs, we will discuss on the time-dependency of messages when we will focus on their numerical implementation of eqns (3) and (4). Assuming the existence of fixed points of these equations, and simplifying the notations with $(\mu_{i \to a}^{(\infty)}(\sigma_i), \hat{\mu}_{a \to i}^{(\infty)}(\sigma_i)) \equiv (\mu_{i \to a}(\sigma_i), \hat{\mu}_{a \to i}(\sigma_i))$, we can write marginal distribution of variable $i$ using the set of incoming fixed-points messages:

$$ \mu_i(\sigma_i) \propto \prod_{a \in \partial i} \mu_{a \to i}(\sigma_i) \qquad (5) $$

Once we have a method to compute marginals, we can employ it to search for a solution. Several search algorithms inspired by marginal computation have shown good results in NAE-SAT [7] or other CSPs [19]. That will be the object of further paragraphs. In addition, this approach allows us to compute statistical observables in order to study the structure of the space of solutions. For instance, in the Bethe approximation, defining $\underline{\mu} = \{\mu_{i \to a}(\cdot), \hat{\mu}_{a \to i}(\cdot); a = \{1, ..., M\}; i \in \partial a\}$, one can write the Bethe free-entropy as follows:

$$ \mathbb{F}[\underline{\mu}] = \langle H[\underline{\mu}] - E[\underline{\mu}] \rangle_{\mu} = \sum_{a=1}^{M} \mathbb{F}_{a}(\underline{\mu}) + \sum_{i=1}^{N} \mathbb{F}_{i}(\underline{\mu}) + \sum_{(ia)} \mathbb{F}_{ia}(\underline{\mu}) \qquad (6) $$

where $H[\underline{p}]$ is the entropy of the uniform distribution defined in (2) rewritten in terms of messages, $E[\underline{p}]$ is the energy of the system rewritten in terms of messages and we define:

$$ e^{\mathbb{F}_a[\underline{\mu}]} = \sum_{\underline{\sigma}_{\partial a}} \psi_a(\underline{\sigma}_{\partial a}) \prod_{j \in \partial a} \mu_{j \to a}(\sigma_j) \qquad (7) $$

$$ e^{\mathbb{F}_i[\underline{\mu}]} = \sum_{\sigma_i} \prod_{b \in \partial i} \hat{\mu}_{b \to i}(\sigma_i) \qquad (8) $$

$$ e^{\mathbb{F}_{ia}[\underline{\mu}]} = \sum_{\sigma_i} \hat{\mu}_{a \to i}(\sigma_i) \mu_{i \to a}(\sigma_i) \qquad (9) $$

On trees, we get exactly the free-entropy of the system by evaluating the Bethe free-entropy at the fixed point of BP: $\mathbb{F}[\underline{\mu}^*] = \phi = \log(Z)$. Still, on loopy graphs, the Bethe free-entropy will provide an approximate value of the exact free-entropy.

### 2.2.2 Case of NAE-SAT

In the following, we will specialize this formalism to NAE-SAT. Since the messages are defined for two values of the variables, we can parametrize the problem with a pair of messages $(\hat{\mu}_{ai}, \mu_{ia}) \in [0, 1]^2$ defined as $\hat{\mu}_{ai} \equiv \hat{\mu}_{a \to i}(\sigma_i = +1)$ and $\mu_{ia} \equiv \mu_{a \to i}(\sigma_i = -1)$. By normalization, we retrieve the message taken on the opposite value. BP equations will be refered as:
$$ \mu_{ia}^{(t+1)} = f(\{\hat{\mu}_{bi}^{(t)}\}) \qquad (10) $$

$$ \hat{\mu}_{ai}^{(t+1)} = \hat{f}(\{\mu_{ja}^{(t)}\}) \qquad (11) $$

where

$$ f(\{\hat{\mu}_{bi}^{(t)}\}) = \frac{\prod_{b \in \partial j \setminus a} \hat{\mu}_{bi}}{\prod_{b \in \partial j \setminus a} \hat{\mu}_{bi} + \prod_{b \in \partial j \setminus a} (1 - \hat{\mu}_{bi})} \qquad (12) $$

$$ \hat{f}(\{\mu_{ja}\}) = \frac{1 - \prod_{j \in \partial a \setminus i} \mu_{ja}}{2 - \prod_{j \in \partial a \setminus i} \mu_{ja} - \prod_{j \in \partial a \setminus i} (1 - \mu_{ja})} \qquad (13) $$

## 2.3 The Cavity Method and Replica Symmetry Breaking (1RSB)

In this paragraph, we will apply the cavity method to the NAE-SAT problem in order to have an analytical description of the space of solutions. A detailed description of the various phase transitions that occur when the density of clauses increases has already been done in K-SAT [11] and XORSAT [12]. Here, we will develop a similar reasoning found in [12][7] to account for the statistical landscape of NAE-SAT.

### 2.3.1 The Replica Symmetry Cavity Method

The Replica Symmetry solution is the more intuitive solution that gives a family of fixed points of BP eqs (10)(11). Instead of considering the update rule for messages, it directly focuses on the fixed points and its properties. We can characterize this solution regarding the free-entropy per variable that can be calculated by taking the expectation of eqn (5) :

$$ f^{RS} = f_{v}^{RS} + \alpha f_{c}^{RS} - \alpha k f_{e}^{RS} \qquad (14) $$

with

$$ f_{v}^{RS} = \mathbb{E}_{l, \{\hat{\mu}_{bi}\}} \log(z_i) = \mathbb{E}_{l, \{\hat{\mu}_{bi}\}} \log[\prod_{b \in \partial i} \hat{\mu}_{bi} + \prod_{b \in \partial i} (1 - \hat{\mu}_{bi})] \qquad (15) $$

$$ f_{c}^{RS} = \mathbb{E}_{\{\mu_{ja}\}} \log(z_a) = \mathbb{E}_{\{\mu_{ja}\}} \log[1 - \prod_{j \in \partial a} \mu_{ja} + \prod_{j \in \partial a} (1 - \mu_{ja})] \qquad (16) $$

$$ f_{e}^{RS} = \mathbb{E}_{\{\mu_{ia}, \hat{\mu}_{ai}\}} \log(z_{ia}) = \mathbb{E}_{\{\mu_{ia}, \hat{\mu}_{ai}\}} \log[\hat{\mu}_{ai}\mu_{ia} + (1 - \hat{\mu}_{ai})(1 - \mu_{ia})] \qquad (17) $$
where the $\mathbb{E}$ denotes expectation with respect to the random variables in subscript. Note that the average over the degree $l$ of variable nodes account for the choice of the random graph model we defined in sec. 2.1.2, namely a Poissonian distribution of intensity $\alpha K$.

The numerical method consists in running BP equations without considering the microscopic random graph underlying structure. At each step of updating rules of eqs (10)(11), we randomly choose a set of incoming messages from an entire population and the outcome shall replace a randomly chosen element of the other population. This method also called **Population Dynamics** can be seen as a **Monte Carlo method** for sampling a probability density. The description of the algorithm is the following :

<caption>**Algorithm 1 RS Population Dynamics (density $\alpha$, degree $k$, size $N$, iterations $T$)**</caption>

1: Initialize $\{\mu_i\}$
2: **for** $t = 1, ..., T$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, .., N$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $i(1), ..., i(k-1)$ uniformly in $\{1, ..., N\}$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\hat{\mu}_j = \hat{f}(\mu_{i(1)}, ..., \mu_{i(k-1)})$
6: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
7: Draw $l$ from $\text{Poisson}(\alpha k)$
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, .., N$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $i(1), ..., i(l)$ uniformly in $\{1, ..., N\}$
10: &nbsp;&nbsp;&nbsp;&nbsp;Set $\mu_j = f(\hat{\mu}_{i(1)}, ..., \hat{\mu}_{i(l)})$
11: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
12: **end for**
13: **return** $\{\mu_i, \hat{\mu}_i\}$

After running a certain time the algorithm, we can compute the free-entropy using the
outcoming distribution of messages according to eqs (15) - (17). We used a population
of $N = 10^5$ and we run the population algorithm for a few values of the density $\alpha$. The
results for $k = 4$ are exposed in Fig.4. For a certain value $\alpha = \alpha_{s,num}^{RS}$, the average
free-entropy vanishes. A negative free-entropy does not have a physical meaning which
can be interpreted as the non-existence of RS solutions in this region.

Using the symmetry of the system, one can easily derive the value of the theoretical
threshold from eqs (14)-(17) which gives :

$$
\alpha_{s,th}^{RS} = - \frac{\log(2)}{\log(1 - 2^{1-k})} \qquad (18)
$$

Our results for different values of $k$ have been put on Tab. 1; one can check they
perfectly fit with the analytical calculation in (18). Above this threshold, no algorithm
is able to find SAT assignments, we will see in the next part that it is confirmed by our
simulations. Still, the RS description of solutions has been shown to be wrong above a
certain threshold, since the entropic landscape follows a much complex scheme, which
have been first found in spin-glasses models. This is the subject of the next paragraph.
<figure><img src="image_7.png" /><figcaption>Figure 4: RS internal free-entropy in NAE-4-SAT. The RS-SAT threshold is located when the free-entropy evaluated at long times vanishes. A refined analysis estimates $\alpha_{s,num}^{RS} = 5.19$.</figcaption></figure>

### 2.3.2 1RSB Cavity Method

In the last section, we used the Bethe approximation to derive the value of the free-entropy of the system. We assumed that random graphs are locally tree-like. This assumption stands as long as correlations along the graph structure decay sufficiently fast in the thermodynamic limit. What happen if this assumption turns to be wrong? We need to introduce the notion of clusters that is detailed in spin-glass literature [15]. In each cluster, we assume that such correlations decay sufficiently fast. It corresponds to sample solutions that are independent in each cluster. According to this hypothesis, one shall replace the Bethe measure $\mu$ defined in (2) by a decomposition into an exponentially large number of Bethe measures restricted to each cluster:

$$ \mu(\underline{\sigma}) = \sum_{n} w_n \mu^n(\underline{\sigma}) \qquad (19) $$

where the index $n$ counts states and the weights $w_n$ are related to Bethe free entropies by $w_n = e^{\mathbb{F}_n/\Phi}$, $\Phi = \sum_n e^{\mathbb{F}_n}$. In order to describe the statistical landscape of these clusters, we need to assume that the number of measures with free entropy $\phi$ is exponentially large $N(\phi) \propto \exp(N\Sigma(\phi))$ where the complexity $\Sigma(\phi)$ characterises the density of states with free entropy $\phi$.

The 1RSB method consists in studying a new system composed by $m$ replicas of the original system, $m$ being also called *Parisi 1RSB parameter*. The new partition function
of the system is then given by

$$ \Phi[m] = \sum_n e^{m \mathbb{F}_n} = \int e^{N[m\phi + \Sigma(\phi)]} d\phi \quad (20) $$

One can retrieve the original problem at $m = 1$, nevertheless computing the partition function $\Phi[m]$ for each $m$ allows to access to the complexity $\Sigma$ by a Legendre transform. Without entering into details, we will refer to the general method invoking auxiliary graphical model described in [12]. Essentially, considering the introduction of the Boltzmann factor $w_n^m$, one can write equivalent equations as (10)-(11) that monitor the evolution of the probability density distributions of each kind of messages, where for each edges (ia) of the original factor graph, we associate two distributions functions $\{Q_{ia}[\mu_{ia}], \hat{Q}_{ai}[\hat{\mu}_{ai}]\}$. Since messages are probability distributions on $\{0, 1\}$, one shall now think in terms of probability distributions of probability distributions. However, we will see that simplifications, either theoretical and numerical ones, can be achieved for special cases $m = 0$ and $m = 1$ which gives already significant information. The most general equations that rule these distributions are then given by:

$$ Q_{ia}[\mu_{ia}] = \frac{1}{\mathcal{Z}_{ia}} \int \prod_{b=1}^{l} d\hat{Q}_{bi}[\hat{\mu}_{bi}] \delta(\mu_{ia} - f(\{\hat{\mu}_{bi}\}) z(\{\hat{\mu}_{bi}\})^m \quad (21) $$

$$ \hat{Q}_{ai}[\hat{\mu}_{ai}] = \frac{1}{\hat{\mathcal{Z}}_{ai}} \int \prod_{j=1}^{k-1} dQ_{ja}[\mu_{ja}] \delta(\hat{\mu}_{ai} - \hat{f}(\{\mu_{ja}\}) \hat{z}(\{\mu_{ja}\})^m \quad (22) $$

where $\mathcal{Z}_{ia}, \hat{\mathcal{Z}}_{ai}$ are normalization constants and $z(\{\hat{\mu}_{bi}\}), \hat{z}(\{\mu_{ja}\})$ are respectively normalization factors that appears in eqs (12)(13). It is remarkable that eqs (21)(22) are strongly similar to BP equations in their structure. Nevertheless solving them numerically in the general case of $m \in \mathbb{R}$ is quite painful, since we have to generate for each edge of a graph, a couple of populations arrays. Details about this type of computation are emphasized in [7, 12].

**Case m=0** This particular case has two main assets. First, choosing $m = 0$ is equivalent to sample solutions uniformly in each clusters. Since, at the SAT/UNSAT threshold, the number of solutions of the problem goes to 0 in the thermodynamic limit, the solutions sampled at $m = 0$ should disappear, but also the other ones that are not sampled. Then, the satisfactory threshold would be given by the value of clause density where complexity disappears: $\Sigma^{m=0}(\alpha_s) = 0$. Secondly, the dynamical threshold that we put forward in the last section should be retrieved by the appearance of a non-zero complexity: $\alpha_d = \min_{\alpha}(\sum^{m=0}(\alpha) > 0)$. All these results have been studied in the case of $k = 3$ in [5]. We will only generalized here their approach to the general $k \ge 3$ in order to determine the density $\alpha_d^+$ above which the first clusters appear.

In order to take the limit $m \to 0$ in eqs (21)(22), one should be careful and consider fully polarized distributions, that corresponds to hard messages in the literature. Using the symmetry of the system $\mu \leftarrow 1-\mu$, one can reparametrized the problem with respect
<figure><img src="image_8.png" /><figcaption>Figure 5: Average value of a population of soft-messages for $k = 4$; the vertical line indicates the location of the transition.</figcaption></figure>

to the distributions of soft-messages $\eta$ arriving to each function node. The most general form of the distribution of clause-to-variable message is given by:

$$ \mu^{\eta}(r) = \eta \nu(r) + \frac{1-\eta}{2} [\delta(r-1) + \delta(r)] \quad (23) $$

where $\mu^{\eta}(r)$ is from now a probability distribution over the space $\{0, 1\}$ parametrized by the distribution function $\mu$ and $\nu$ is the distribution of soft messages that verify the normalization condition $\int \nu(r)dr = 1$. If we define $\rho(\eta)$ the probability functions of $\eta$, we can derived from (21)(22) a new self-consistency equation as follows:

$$ \rho(\eta) = \sum_{l_1=0}^{\infty} \mathcal{P}_{\alpha k}(l_1) \dots \sum_{l_{k-1}=0}^{\infty} \mathcal{P}_{\alpha k}(l_{k-1}) \int \prod_{i_1=1}^{l_1} d\rho(\eta_{i_1}) \dots \prod_{i_{k-1}=1}^{l_{k-1}} d\rho(\eta_{i_{k-1}}) \quad (24) \\ \times \delta \left[ \eta - 1 + \frac{1}{2^{k-2}} \left(1 - \frac{\prod_{i_1=1}^{l_1} \eta_{i_1}}{A_{i_1}}\right) \dots \left(1 - \frac{\prod_{i_{k-1}=1}^{l_{k-1}} \eta_{i_{k-1}}}{A_{i_{k-1}}}\right) \right] $$

where $\mathcal{P}_{\beta}$ stands for the Poissonian distribution of intensity $\beta$ and $A_i = 2 \prod_{j=1}^{i_1} \frac{1+\eta_i}{2} - \prod_{j=1}^{i_1} \eta_i$. Taking the limit $m \to 0$ corresponds to count the contribution of the whole set of clause-to-variable messages distributions except the two that are all hard-messages with same polarization and the ones that contain contradictions, namely at least two differently polarized hard-messages. We can solve this equation using a population dynamics algorithm. We start from a population $\eta$ randomly distributed in $[0, 1]$ and we proceed as follows:

* draw $l_1$ from a Poissonian distribution, take $l_1$ elements at random in the current population, compute $g_1 = \prod_{i_1=1}^{l_1} \eta_{i_1}$ and $A_{i_1}$.
<table><caption>Table 1: Various thresholds of the bicoloring problem. The first three columns have been obtained by methods explicited in section 2. The fourth column shows results for RDBP, that are discussed in 3.1. The last column have been adapted from [4], since the k-regular hypergraphs and Poissonian hypergraphs present the same unstability of the paramagnetic solution.</caption><thead><tr><th>k</th><th>α<sub>d,+</sub></th><th>α<sub>d</sub></th><th>α<sup>RS</sup><sub>s,num</sub></th><th>α<sub>num</sub></th><th>α<sub>stab</sub></th></tr></thead><tbody><tr><td>3</td><td>1.91</td><td>-</td><td>2.41</td><td>1.91</td><td>1.5</td></tr><tr><td>4</td><td>4.06</td><td>-</td><td>5.19</td><td>4.47</td><td>4.08</td></tr><tr><td>5</td><td>7.95</td><td>9.06</td><td>10.74</td><td>8.85</td><td>11.25</td></tr><tr><td>6</td><td>15.09</td><td>17.05</td><td>21.83</td><td>-</td><td>32.03</td></tr></tbody></table>

* repeat the previous step $k-1$ times

* replace a random element of the population by:

$$ \eta - 1 + \frac{1}{2^{k-2}} \left(1 - \frac{g_1}{A_{i_1}}\right) \dots \left(1 - \frac{g_{k-1}}{A_{i_{k-1}}}\right) $$

* repeat the three last steps until convergence of the distribution

For different values of $k$, we run this program for $N = 10^5$ varying the density $\alpha$. The program has to stop when the empirical distribution is stabilized. One can observe that, for $\alpha < \alpha_{d,+}$, the distribution is a delta function on $\eta = 1$ and the average value is $\langle\eta\rangle = 1$ since for $\alpha > \alpha_{d,+}$, the distribution is non trivial and $\langle\eta\rangle < 1$. In Fig.5 we display this behaviour by plotting the mean value of the distribution in the case $k = 4$. We run this algorithm for different values of $k$ and the corresponding threshold $\alpha_{d,+}$ for each are shown in Tab. 1.

**Case m=1** In the other hand, the dynamical threshold cannot be obtained by the Parisi computation at $m = 0$ since we need an information about the relative number of solutions between clusters, which is given by the entropic term $e^{Nm\phi}$ in eq. (20). At $m = 1$, we are preferentially sampling solutions inside clusters that dominate in the total number of solution, the ones with free-entropy that verify $\Sigma'(\phi) = -1$. In other words, under the dynamical threshold, the term $e^{N(\phi+\Sigma(\phi))}$ remain exponentially smaller than the total number of solutions $Z$. Thus, the appearance of a point $\phi^*$ with $\Sigma'(\phi^*) = -1$ determine the value of the dynamical threshold $\alpha_d$. One can show that equivalently, a non-trivial solution appears for eqs (21)(22) with $m = 1$ for $\alpha = \alpha_d$.

In the case $m = 1$, we can simplify eqs (21)(22) following the method described in [12, 16]. Therefore, one shall define new distributions functions on $[0, 1]$ as:

$$ R_{\sigma}(r) = (1 + \sigma(1 - 2r))\mathbb{E}\{Q_{ia}[\mu_{ia}]\} \quad (25) $$

$$ \hat{R}_{\sigma}(\hat{r}) = (1 + \sigma(1 - 2\hat{r}))\mathbb{E}\{\hat{Q}_{ai}[\hat{\mu}_{ai}]\} \quad (26) $$
where $\sigma \in \{-, +\}$ and $\mathbb{E}$ stands for the average over distributions functions on each edge and over graph structure. These functions then follow the following evolution:

$$ \hat{R}_{\sigma}(\hat{r}) = \int \sum_{\{\sigma_j\}} \pi(\{\sigma_j|\sigma\} \prod_{j=1}^{k-1} dR_{\sigma_j}(r_j)\delta(\hat{r} - \hat{f}(\{r_j\})) \quad (27) $$

$$ R_{\sigma}(r) = \sum_{l=0}^{\infty} \mathcal{P}_{\alpha k}(l) \int \prod_{i=1}^{l} dR_{\sigma}(i) \delta(r - f(\{\hat{r}_i\})) \quad (28) $$

where $\pi(\{\sigma_j|\sigma\} = \frac{1-\mathbb{I}(\sigma_1=\dots=\sigma_{k-1}=\sigma)}{2^{k-1}-1})$. These equations has the advantage to be much more easy to compute numerically since we now deal with simple distributions that require each only one population instead of one population for each edge in eqs (21)(22). The computational procedure is synthesized in the Algo. 2.

**Algorithm 2** 1RSB Population Dynamics for Parisi parameter $m = 1$ (density $\alpha$, degree $k$, size $N$, iterations $T$).

1: Initialize $(r_i^+, r_i^-)$ by (0,1)

2: **for** $t = 1, \dots, T$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1, \dots, N$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $\sigma = +, -$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $\sigma_1, \dots, \sigma_{k-1}$ from $\pi(\{\sigma_j|\sigma\})$

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $i(1), \dots, i(k-1)$ uniformly in $\{1, \dots, N\}$

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\hat{r}_i^\sigma = \hat{f}(r_{i(1)}^{\sigma_1}, \dots, r_{i(k-1)}^{\sigma_{k-1}})$

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: **end for**

10: **for** $i = 1, \dots, N$ **do**

11: &nbsp;&nbsp;&nbsp;&nbsp;**for** $\sigma = +, -$ **do**

12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $l$ from $\text{Poisson}(\alpha k)$

13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $i(1), \dots, i(l)$ uniformly in $\{1, \dots, N\}$

14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $r_i^\sigma = f(\hat{r}_{i(1)}^\sigma, \dots, \hat{r}_{i(l)}^\sigma)$

15: **end for**

16: **end for**

17: **end for**

18: $\{r_i^+, r_i^-, \hat{r}_i^+, \hat{r}_i^-\}$

In principle, one should run this program until the population is time-independent. Nevertheless, we took $T = 2000$, not too large, so as to prevent fluctuations due to finite-size effects to deteriorate our measurements. We also use damping to help the relaxation of the system; this technique will be briefly explicated in the next part. This algorithm can be easily implemented and since its complexity is $O(N)$, we choosed a population $N = 10^6$ to minimize fluctuations.

In the 1RSB picture, one can define the inter and intra overlaps as
$$q_0 \equiv \mathbb{E}\left[(\int dQ_{ia}[\mu_{ia}](2\mu_{ia}-1))^2\right] \qquad q_1 \equiv \mathbb{E}\left[\int dQ_{ia}[\mu_{ia}](2\mu_{ia}-1)^2\right] \quad (29)$$

Then using definitions (22)(23), we can rewrite these numbers as

$$q_0 = \left[ \left\langle r - \frac{1}{2} \right\rangle_{R_+}^2 + \left\langle r - \frac{1}{2} \right\rangle_{R_-}^2 \right] \qquad q_1 = \left[ \left\langle \left(r - \frac{1}{2}\right)^2 \right\rangle_{R_+} + \left\langle \left(r - \frac{1}{2}\right)^2 \right\rangle_{R_-} \right] \quad (30)$$

where $\langle \rangle_{R_\sigma}$ refers to the expectation of the distribution $R_\sigma$. One can compute these numbers using the empirical averages over the populations obtained by the upper algorithm. We plotted in Fig.6 (right panel) the difference between overlaps obtained for $k=5$. One can observe a transition phase at $\alpha_d = 9.07$ (with a good precision considering $N = 10^6$) since $q_1 - q_0 \to 0$ for $\alpha < \alpha_d$ and this difference stabilizes to a positive value for $\alpha > \alpha_d$. We can think about this typical 1RSB transition phase as a discontinuous transition with respect to the parameter order $q_1 - q_0$.

In the other hand, in the cases $k=3,4$, the system seems to display a continuous transition. Indeed, in the Fig.6 (left panel), one can see that the non-trivial solution around $\alpha = 4.1$ (rough estimation considering finite size effects) for $k=4$ is underscored by a non zero overlaps difference but very close to 0. It suggests an other mechanism than for $k=5$. This singular transition that happen below the dynamical threshold for the cases $k=3,4$ have been already highlighted in NAE-3-SAT (at $\alpha = 1.5$) in [5]. It corresponds to the instability of the paramagnetic RS solution $\{(\mu_{ia}, \hat{\mu}_{ai}) = (1/2, 1/2)\}$ through the 1RSB solution. It is also refered as the stability threshold $l_{stab}$ in the k-regular graphs model [4], which values have been adapted with our notations in Tab.1. Hence, the present method allows to determine $\alpha_d$ only for $k \ge 5$, where the paramagnetic instability arises above the dynamical threshold.

Let us now explain what happen in terms of clusters. The inter overlap $q_1$ is related to the typical distance between two solutions in the same cluster when $q_0$ is related to the typical distance between two solutions in different clusters. Let us consider that a variable is fixed to a given value (for instance, $\sigma = -1$, which corresponds to the initial condition $\{(r, r_i^-) = (0, 1), \forall i\}$ in our algorithm). Then, iterating $T$ times eqs (27)(28) is equivalent to infer the value of variables located at a distance $T$ of the fixed variable in a tree graph. Once it is possible to infer values of variables given a specific boundary condition, it means it selects clusters that are compatible with this boundary condition. If the system is above its dynamical threshold, when such a reconstruction is possible, we observe a transition phase according to the limit of $q_1 - q_0$. This transition phase has some theoretical implications in term of the free-boundary measure defined in (2); indeed it is associated to a point-to-set correlation built on this measure [16, 11] and which corresponds to $q_1 - q_0$ in the spin glass formalism.

## 3 BP-inspired Search Algorithms

In this section, we intend to describe the algorithms we used to look for solutions of the NAE-k-SAT problem. The commun setting is a message passing procedure that allows
<figure><img src="image_9.png" /><figcaption>Figure 6: Difference of intra and inter overlaps in NAE-SAT. Left panel: $k = 4$ over a population $N = 3.10^4$; scale is logarithmic. Right panel: $k = 5$ over a population $N = 10^6$.</figcaption></figure>

to compute approximated values of the marginal probabilities instead of the exact ones. After presenting the decimation algorithm that have been investigated both analytically and numerically in [19], we will present a similar algorithm which finds solutions above the dynamical threshold. Finally, we will present an other algorithm that is specific to NAE-$k$-SAT, that is still under investigation.

# 3.1 Random Decimation

## 3.1.1 Description

In the decimation procedure, we start with a random instance and we look for a solution of this given instance, it means a set of $N$ boolean variables that verify the constraints of the random formula. The idea is to fix an set of variables to their most probable value with respect to the BP equations marginals evaluated by eq (5). Basically once we fixed a few variables, the factor graph is reduced to a new factor graph in which new estimates can be computed by running the BP equations between variables nodes that are undecimated and function nodes that constrain them.

At each decimation step, one shall select a set of variables that will be assigned to their most probable value. One can think about different way to do the latter. The first consists in choosing uniformly these variables at random. This process, that we will refer as *Random Decimation Belief Propagation* (RDBP) is known to be less efficient than others as we will see in the next section. However the analytical description is feasible in this case, as the authors of [19] did for the K-SAT model. In the other hand, it brings numerical evidence of some features of the statistical landscape that we presented in the first part. The algorithm is synthetized in the following scheme :

We wrote the algorithm in C because it is quite apropriate to run fast calculation on graphs. Indeed, each variable node and each function node is represented in C by
**Algorithm 3 Random Decimation BP (size N, degree k, density α)**

1: generate a random instance of the NAE-SAT problem with parameter k and α

2: initialize messages {$\mu_{ja}^{(0)}, \hat{\mu}_{bi}^{(0)}$} with uniform distribution on $[0, 1]$

3: **while** the number of unfixed variables is > 0 **do**

4: **for** t = 1, ..., T (BP iterations) **do**

5: **if** $\Delta\mu < \delta$ exit

6: compute messages {$\mu_{ja}^{(t)}$} according to (12) for all variable nodes i not fixed

7: compute messages {$\hat{\mu}_{bi}^{(t)}$} according to (12) for all function nodes b

8: compute the update shift $\Delta\mu = \min_{(ia)}(|\mu_{ia}^{(t)} - \mu_{ia}^{(t-1)}|, |\hat{\mu}_{ai}^{(t)} - \hat{\mu}_{ai}^{(t-1)}|))$

9: **end for**

10: choose uniformly a fraction $r_d$ of unfixed variables and assign spin values according to the marginal estimates

11: set all messages sent by these variables to hard messages {$+1, -1$} according to the spins

12: **end while**

13: number of unsat clauses

a *class* that keeps in memory every input/output messages which is concerned by this node, the memory that stores a message between two nodes is allocated once and we can rapidly access to this memory from both variable and function nodes that exchange this message.

Since we first noticed that BP was not converging for large value of $\alpha$, we decided to stop the BP iterations after $T = 2000$ steps. In the other hand, we set the criterium of convergence of BP to $\delta = 10^{-4}$, asking our programm to be fast. We did not find useful to lower the value of this criterium. Indeed we checked that after a few decimation steps, messages converge exponentially by using *damping*. It means that when we update a message $m$, we do not assign to the new value $m_{new}$ but rather to a weighted sum $m \leftarrow 0.4m_{new} + 0.6m$. Note that for convenience, we used message damping but it has been shown in [18] that marginal damping leads faster to convergence for a certain value of the damping. In addition, we can reasonably assign errors during estimation of marginals to the region where BP does not converge, by opposition to the region where $\Delta\mu < \delta$.

The main advantage of this algorithm is its complexity in $O(N)$ as long as we choose a large decimation rate $r_d$. One could think that decimate variables in a simultaneous way would induce errors in the estimate of marginals for the successive reduced factor graphs. However, after trying different decimation rates, we already observed good performance of the algorithm for $r_d = 2\%$. Eventually, one could improve its performances either by increasing the number of variables $N$, either by decreasing the decimation rate $r_d$. The best achievement would be to use one after one variable to decimate, but in this case the algorithm is $O(N^2)$ which is quite painful in time and not clever for designing a research algorithm.
<figure><img src="image_10.png" /><figcaption>Figure 7: Probability of finding a solution for $k = 4$ (left panel) and $k = 5$ (right panel).</figcaption></figure>

### 3.1.2 Probability of convergence

The previous scheme of the algorithm has to be performed for each instance we generated. In order to determine the performance of the algorithm, we run the program on 100 instances and computed the frequency of finding a solution at the end of the decimation process (which corresponds to 0 unsat clauses). Our results for $k = 4$ and $k = 5$ are shown in Fig.7. We observe that for $\alpha < \alpha_{num}$, the algorithmic finds easily a solution of the problem, since for $\alpha > \alpha_{num}$, in the large $N$ limit, the probability to solve an instance goes to 0. We can determine with a good precision the value of $\alpha_{num}$ by taking the crossing point between all curves that are plotted for different $N$. For $k = 4$, we get $\alpha_{num}^{k=4} = 4.465 \pm 0.003$, the precision being given by statistical errors in computing the frequency. For $k = 5$, we get an approximated value of the algorithmic threshold $\alpha_{num}^{k=5} \approx 8.85 \pm 0.05$. One has to precise that we use only 50 instances for the latter case, and then the statistical fluctuations are quite important.

In 2.3.2., we located the dynamical threshold $\alpha_d$ of NAE-SAT that corresponds to the appearance of a large number of clusters containing an exponentially large number of solutions. In $k = 5$, we effectively noticed that the previous search algorithm stops finding solutions at a value of clause density close to the dynamical threshold $\alpha_d = 9.09$ (but still smaller), which would suggest that Belief Propagation is not anymore efficient above this threshold. Note that the algorithmic threshold stands far above the sub-dynamical threshold $\alpha_{d,+}$ (see Tab. 1) which indicates that it is not the apparition of newt clusters that make RDBP fail. It simply suggests that BP iterations selects solutions that are typical, thus in the most numerous cluster, where no isolated metastable states exist, that could trap the process.

For $k = 3, 4$, we cannot compare this algorithmic threshold with the result of the same 1RSB population dynamics for the reasons we already discussed in 2.3.2. In order to compare numerical thresholds found with random decimation and dynamical threshold, we would need a more detailed study of the generalized partition function (20), see 2.3.2. It could be part of a future work to confirm our expectations on the cases $k = 3, 4$.
### 3.1.3 Convergence of BP and precursor

The detailed analysis of the random decimation have been done in the case of K-SAT [19]. Inspired from it, we have numerically investigated the lack of convergence during decimation. For this purpose, we computed the statistics over the number of instances for which BP converges (before $T = 2000$ steps) at a given clause density. The results for $k = 4$ are exposed in Fig.8. We did not display the results in the range $\theta > 0.5$, since the BP iterations always converge in this region. We observe in the upper panels that the frequency of non convergence takes a maximum value for $\theta^*(\alpha)$ which corresponds to the decimation fraction for which the algorithm cumulates more errors in terms of marginal estimation.

Note that the behaviour of the BP iterations at the beginning of the decimation process is highlighted in the lower panel; a fixed points at $\theta = 0$ is obtained for low values of the density (for instance $\alpha = 3.9$), and convergence is lost above $\alpha \approx 4.0$. We believe it is due to the paramagnetic instability that appears at $\alpha_{stab} = 4.08$, see Tab. 1. One could check this hypothesis and find the precise location of this transition by increasing the number of variables. In addition, for $k = 5$, the same study has revealed that the BP iterations would always converge in the $\theta \to 0$ limit, which is in agreement with $\alpha_{stab}^k > \alpha_d^k$, see Tab. 1.

Moreover, one can denote another transition that has to be precised in a later study. Indeed, for low density, the peak at $\theta^*$ vanishes at the thermodynamic limit (left upper panel) while it remains at higher density (center and right upper panel). One can already notice that this transition around $\alpha \approx 4.0$ occurs far below the algorithmic threshold $\alpha_{num} \approx 4.46$ and cannot account for the lack of solutions above this latter. However these preliminary results suggest the same scenario found in [19] applies to NAE-SAT as well. They showed that $\theta^*$ is the precursor of the algorithmic barrier that occured slightly below the dynamical threshold. Therefore, it would be interesting to derive a theoretical value for the algorithmic barrier in RDBP in order to confirm these expectations.
<figure><img src="image_11.png" /><figcaption>Figure 8: Upper panels : Probability of non convergence of BP iterations during the decimation process, plotted respectively from left to right for $\alpha = 4.0, 4.1, 4.2$. Lower panel : Non convergence of BP at low values of $\theta$, for $N = 5000$. All panels refers to $k = 4$. The upper results have been obtained using a constant decimation rate $r_d = 2\%$, while the lower results have been obtained varying $r_d$ to access the small values of $\theta$.</figcaption></figure>

In conclusion of this part, we have shown in our simulations that the appearance of clusters, i.e. regions well separated in the space of solutions, prevent the random decimation algorithm to be efficient after the dynamical threshold. It indicates that this algorithm choose sequentially marginals in a way that, for example marginals of two distinct nodes are not compatible in the 1RSB solutions of the same cluster, prohibiting the merging of a global solution.

## 3.2 Guided Decimation Algorithm

The dynamical threshold can be overtaken by different algorithms. One of them called *Survey Propagation* has been demonstrated to be very efficient for this purpose in several CSPs [3][5]. It is largely inspired by the clustering picture that we have already invoked in the introduction. Basically, each cluster can be associated to a set of hard-messages, also refered as *warnings*. Then, the SP algorithm computes the marginals from the distributions of the WP messages which allow to select independently different clusters. However, we did not probe the NAE-SAT problem with this algorithm since we had
first used an improvement of the latter decimation algorithm that allow to cross the dynamical barrier too. The program we called *Guided Decimated Belief Propagation* (GDBP) consists only to change the way variables are fixed during the decimation process. The step 11 of the Algo. 3 has to be replaced by a new step in which the a fraction $r_d$ of unfixed variables are pinned according to their most biased marginals. Basically, we had to sort the unfixed variables with respect to $|\mu_i - \frac{1}{2}|$. It turns to be equivalent to consider variables with highest biased marginals as sources of warnings.

<figure><img src="image_12.png" /><figcaption>Figure 9: Ground state energy of decimated BP in NAE-3-SAT, $N = 10^4$. Vertical bars stands for dynamical ($\alpha_d$) and SAT/UNSAT thresholds ($\alpha_s$), values being taken from [5].</figcaption></figure>

The complexity of the algorithm is now $O(N \log(N))$, the log coming from the sorting step. For our study, we have used $N = 10000$ which is enough for the preliminary study. For a future work on powerful machines involving $N \approx 10^6$ variables, we have already integrated a *partial sorting algorithm*. Indeed one only need to sort $r_d N$ on each decimation step, which would reduce the complexity to $O(N + r_d N \log(r_d N))$.

In the Fig.9, we plotted the internal energy averaged over 100 instances obtained at the end of the guided decimation for $k=3$. As soon as the minimum energy reached by the algorithm is close to 0, it tells us that most of random instances are solvables. GDBP turns out to be more efficient than RDBP since it allows to solve instances for $\alpha_d < \alpha < \alpha_s$. Still, there is no clear connection between the performance of this algorithm and the solution space phase transitions, since no analytical solutions exist. In order to understand from which kind of cluster this algorithm is able to select solutions, one could suggest the study of their entropic location as it has been done for an other successful algorithm known as *Belief Propagation Reinforcement (BPR)* [7].
## 3.3 Biased Decimation Algorithm

In this last section, we will present an other approach we developed in order to improve the efficiency of the previous search algorithms. Based on the change of the measure (2), we aimed to find solution on hard instances at $\alpha > \alpha_d$ without using a sorting step like in GDBP, but restricting the space of solutions.

Moreover, several recent studies have emphasized the structure of clusters into two classes, *frozen clusters* and *unfrozen clusters* [7]; the frozen clusters being those in which some variables are assigned to the same value in all solutions of the same cluster. A general method that can be extended to all CSPs, called *whitening process*, have been extensively studied in [4], it allows to characterize a solution as being frozen or not. In NAE-SAT and KSAT, this dynamics has been shown to be equivalent to the *Warning Propagation* (WP) where BP messages are replaced by hard-messages. We will not enter into technical descriptions but we can summarize the process as follows. Given a solution $\sigma$ of the problem, one can initialize the WP equations, and iterate them until convergence. If some of the variables have stayed in the same polarization during the process, the solution is frozen. It is already known that any search algorithm for CSPs can only display unfrozen solutions. Indeed, the presence of frozen variables, for which the estimation of the marginals has to be done without any error, make impossible the quest for such a solution. In fact, marginals are estimated according to heuristic methods like RS and 1RSB and therefore are always approximated.

In this context we believe that search algorithm could be modified such that it discards marginal estimations that are likely to be close to hard fields. For instance, in a clause involving $k$ variables, we recall that the two configurations of the constrained variables where spins are all polarized in the same direction are not allowed. Moreover, in a given clause, the $2k$ configurations displaying a unique variable different from the others should be disadvantaged. One shall define a new measure as :

$$ \mu_{\epsilon}(\sigma_i) = \frac{1}{Z(\epsilon)} \prod_{a} \psi_a(\sigma_{\partial a}) \lambda_a^{\epsilon}(\sigma_{\partial a}) \quad (31) $$

where $Z(\epsilon)$ is a normalization constant that depends on a parameter $\epsilon$ that rules the reweighting functions $\lambda_a^{\epsilon}$ defined as :

$$ \lambda_a^{\epsilon}(\sigma) = \begin{cases} e^{-\epsilon} & \text{if } \exists! i \in \partial a, \forall j \in \partial a \setminus i, \sigma_i = -\sigma_j \\ 1 & \text{else} \end{cases} \quad (32) $$

Use of similar reweightings in the uniform measure has been introduced to evaluate rigorous upper bound of KSAT [2]. The authors of [10] generalized this approach and studied the remarkable properties of the generated instances. Inspiring from this latter, we derived a modified set of BP equations by replacing eq. (13) by :

$$ \hat{f}_{\epsilon}(\{\mu_{ja}\}) = \frac{P_{+}}{P_{+} + P_{-}} \quad (33) $$
where

$$P_+ = 1 - \prod_{j \in \partial a \setminus i} \mu_{ja} + (e^{-\epsilon} - 1) \left[ \prod_{j \in \partial a \setminus i} (1 - \mu_{ja}) + \sum_{j \in \partial a \setminus i} (1 - \mu_{ja}) \prod_{k \in \partial a \setminus \{i,j\}} \mu_{ka} \right] \quad (34)$$

$$P_{-} = 1 - \prod_{j \in \partial a \setminus i} (1 - \mu_{ja}) + (e^{-\epsilon} - 1) \left[ \prod_{j \in \partial a \setminus i} \mu_{ja} + \sum_{j \in \partial a \setminus i} \mu_{ja} \prod_{k \in \partial a \setminus \{i,j\}} (1 - \mu_{ka}) \right] \quad (35)$$

Note that we retrieve the original eqs (10)(11) taking $\epsilon = 0$.

At this stage, we can already investigate the effect of the bias $\epsilon$ on the behaviour of the population dynamics close to the dynamical threshold. For this end, it is useful to repeat the work of the paragraph 2.3.2 by computing the new effective dynamical threshold $\alpha_d(\epsilon)$, replacing the function $\hat{f}$ in line 7 of Algo. 2 by its new expression in (28). The results are reported in the phase diagram of Fig.10. It turns out that non trivial solutions appear for $\alpha_d(\epsilon) > \alpha(0)$ for a certain range $\epsilon \in [0, \epsilon_{max}]$. Out of this range, we expect the search of solutions being worse than in the unbiased case $\epsilon = 0$.

<figure><img src="image_13.png" /><figcaption>Figure 10: Left panel: Phase diagram of the biased NAE-SAT. The hard region is located on the right of the delimited line $\alpha_d(\epsilon)$. Right panel: Success Probability of the biased-BP, for $N = 1000$.</figcaption></figure>

Additionally, considering the shift of the hard region to upper values of $\alpha$ for $0 < \epsilon < \epsilon_{max}$, we run a modified RDBP algorithm, replacing the BP equations by the new ones, hoping that it would break the algorithm limitation due to clustering. For a matter of time, we chose $N = 1000$ and we averaged over 50 instances. We set $\epsilon = 0.12$ where the corresponding threshold is close to its maximum in the phase diagram. The probability of finding a solution is plotted in Fig.10, where error bars take into account the number of tested instances for each program. Therefore, one can clearly distinguish the two lines, the one corresponding to RDBP in 3.3 ($\epsilon = 0$) being overtaken by the one of the biased-BP ($\epsilon = 0.12$). Actually, we first run a program for $k = 4$, with a larger
$N$ but the influence of the bias was barely appreciable, that is why we did not present here these results.

Obviously, the use of random decimation search algorithm as a way to probe the efficiency of this method has still an academic interest, since we have already shown that another algorithm based on the original BP equations are also efficient to break the dynamical threshold. Nevertheless, the study of the decimation process is analytically doable [19], hence we aimed in the future to compare our numerical results to the theoretical ones.

One should be tempted to combine our two last algorithms and run guided decimation with bias. Still, we do not have tested this possibility since we denote that the guided algorithm already selected solutions that were exempt of fully-polarized configurations. This interesting feature need to be clarified to understand how fixing most polarized variables allow to access to solutions in the hard region.

# 4 Conclusion

In summary, we have reported a numerical study of the NAE-SAT problem to put forward the relation between heuristics based on the Cavity Method and the performance of a few BP-inspired message passing algorithms. We have described different methods to access to transitions occuring when the density of clause increases. The first method, also refered as d1RSB (Parisi parameter $m = 0$) allows to recover the location of the apparition of the first clusters. Note that there is no evidence that this transition affects our search algorithms. However, we located the dynamical threshold by a second method using Parisi calculation for $m = 1$. We showed that this threshold is related to algorithmic barrier in the search procedure. Indeed, we run a first search algorithm, the RDBP, that fails to find solutions above this threshold. A modifying version of this algorithm have been presented, using a reweighting of the empirical measure. Its average problem have been investigated under 1RSB method, displaying an interesting phase diagram where hard regions are shifted to higher density values. In order to break the algorithmic barrier, we presented two other BP-inspired algorithms, one is based on the previous biased measure and the other one revises the decimation procedure.

## 4.1 Perspectives

The phase diagram obtained in the last part seems to be very encouraging for the future design of algorithms. We believe this reweighting method to be much effective for larger $k$, since $k = 5$ works better than $k = 4$. First, future simulations would be performed by increasing the number of BP steps, reducing the decimation rate, and increasing the number of tested instances. The analysis of the convergence of BP equations have revealed to be essential for the understanding of the behaviour close to the threshold and need to be performed in detail in the GDBP and the biased-BP. A related aspect concerns the treatment of hard-messages since it gives some quantitative information over the accumulation of errors during the decimation process. It would be interesting to investigate that direction to understand why both GDBP and biased-BP are efficient
above the dynamical threshold. Eventually, an analytical work which allows to compute thermodynamic quantities will complete this numerical approach and this would provide sufficient matter to publish our results in a paper.

## 4.2 Acknowledgments

I would like to thank all my colleagues of the Physics Department of La Sapienza University in Rome. Thanks to the students who were graduating too, for our long discussions and their lovely roman accent. Thanks to the whole team of the Ciclofficina Populare Don Chisciotte for their enthusiasm and their knowledge in mechanics (not always statistical ones). I would keep a very good impression of the working conditions here in Rome. Actually I am very greatful to the ARPE program for giving me the opportunity to do this internship. I am still wondering if I will continue this work here in Rome or in Paris where I will attend to a conference in the beginning of July, hoping to have stimulating discussion with other influent scientists in that domain.

I think it is always positive to move into an other country for a research intership, learn a new language, meet some great scientists. In the other hand, I also wonder myself about my perspectives of career in research, especially in theoretical physics. I first though that it would demand too much work to spend three years of PhD doing something very specific and quite difficult but I think it is worth to reach this level of personal achievement. As soon as I will find a good collaboration, a few people I liked to work with and an economical situation which is sufficient, I would be very pleasant to proceed in this way. I will always defend the public research as a institution where hierarchical constraints do not exist and where individuals can fully enjoy their activity as being part of their own construction.

In that manner, I would give a special thank to my supervisor Federico Ricci-Tersenghi for its precious advice, and its kindness. I really appreciated our rich exchanges and his enthusiasm, even when I thought my work was not satisfying. I look forward to work with him in my future collaborations and thereby maintain a great professional relation.

At last, I will never thank enough my girlfriend Marina. You always supported me in the tough moments, even through a ten thousand kilometers barrier that separate us. I dedicated this work to you, I know you will love to read it for entertainment !

## References

[1] D. Achlioptas and A. Coja-Oghlan. Algorithm barriers from phase transitions. *Proceedings of FCOS 2008*, page 793, 2008.

[2] D. Achlioptas and Y. Peres. The threshold for random k-sat is $2^k \log(2) - O(k)$. *Journal of the AMS*, 17(4):947-973, 2004.

[3] A. Braunstein, M. Mézard, and R. Zecchina. Survey propagation: an algorithm for satisfiability. *Random Structures and Algorithms*, 27(2):201-226, 2005.
[4] A. Braunsteina, L. Dall'Asta, G. Semerjian, and L. Zdeborova. The large deviations of the whitening process in random constraint satisfaction problem. *J. Stat. Mech.*, 5:053401, 2016.

[5] T Castellani, V. Napoleano, F. Ricci-Tersenghi, and R. Zecchina. Bicoloring random hypergraphs. *J. Physics A : Math. Gen.*, 36:11037-11053, 2003.

[6] Stephen A. Cook. The complexity of theorem-proving procedures. *Proc. 3rd STOC*, pages 151-158, 1971.

[7] L. Dall'Asta, A. Ramezanpour, and R. Zecchina. Entropy landscape and non-gibbs solutions in constraint satisfaction problems. *Phys. Rev. E*, 77(031118), 2008.

[8] J. Ding, A. Sly, and N. Sun. Satisfiability threshold for random regular nae-sat. *Commun. Math. Phys.*, (341):435-489, 2016.

[9] P. Karp. *Complexity of Computer Computations*. Number p85. 1972.

[10] F. Krzakala, M. Mezard, and Lenka Zdeborova. Reweighted belief propagation and quiet planting for random k-sat. *Journal on Satisfiability, Boolean Modeling and Computation*, 8:149-171, 2014.

[11] F. Krzakala, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborová. Gibbs states and the set of solutions of random constraint satisfaction problems. *Proceedings of the National Academy of Sciences*, 104(25):10318-10323, 2007.

[12] M. Mézard and A. Montanari. *Physics, Information, Computation*. Oxford Press, 2009.

[13] M. Mézard and A. Montanari. *Physics, Information, Computation*. Oxford Press, 2009, p 407.

[14] M. Mézard, G. Parisi, and R. Zecchina. Analytic and algorithmic solution of random satisfiability problems. *Science*, 297(5582):812-815, 2002.

[15] M. Mézarda and G. Parisi. The cavity method at zero temperature. *J. Stat. Mech.*, 111(1/2), 2003.

[16] A. Montanari, F. Ricci-Tersenghi, and G. Semerjian. Clusters of solutions and replica symmetry breaking in random k-satisfiability. *J. Stat. Mech.*, 2008.

[17] C. H. Papadimitriou. *Computational complexity*. 1994.

[18] M. Pretty. A message-passing algorithm with damping. *J. Stat. Mech. : Theory and Experiment*, 2005.

[19] F. Ricci-Tersenghi and G. Semerjian. On the cavity method for decimated random constraint satisfaction problems and the analysis of belief propagation guided decimation algorithms. *J. Stat. Mech. : Theory and Experiment*, 2009.