# Counting and Hardness-of-Finding Fixed Points
in Cellular Automata on Random Graphs

Cédric Koller  Freya Behrens  and Lenka Zdeborová

###### Abstract

We study the fixed points of outer-totalistic cellular automata on sparse
random regular graphs. These can be seen as constraint satisfaction problems,
where each variable must adhere to the same local constraint, which depends
solely on its state and the total number of its neighbors in each possible
state. Examples of this setting include classical problems such as independent
sets or assortative/dissasortative partitions. We analyse the existence and
number of fixed points in the large system limit using the cavity method,
under both the replica symmetric (RS) and one-step replica symmetry breaking
(1RSB) assumption. This method allows us to characterize the structure of the
space of solutions, in particular, if the solutions are clustered and whether
the clusters contain frozen variables. This last property is conjectured to be
linked to the typical algorithmic hardness of the problem. We bring
experimental evidence for this claim by studying the performance of the
belief-propagation reinforcement algorithm, a message-passing-based solver for
these constraint satisfaction problems.

## Introduction

Systems with seemingly simple local interactions have been known to possibly
present complex and unexpected behavior at a macroscopic scale. Cellular
automata (CAs) are a prime example of this phenomenon. Conway’s game of life
[1], probably the most famous cellular automaton, has for example been proven
to be Turing complete [2]. Self-reproducing cellular automata have been
studied by von Neumann as early as the 1940s [3]. Even elementary cellular
automata have been shown to exhibit complex behaviors [4], and some of them
have been proven to be Turing complete [5]. CAs are usually defined on regular
grids, which impose strong correlations between the cells. Following [6], we
relax this structure by considering cellular automata evolving synchronously
on random regular graphs. This will enable analysis of their behaviour in the
large size limit.

In this paper, we study the fixed points, or stationary configurations, of
these cellular automata on random regular graphs. The stationary
configurations can be seen as solutions to constraint satisfaction problems
(CSPs) where each variable must adhere to the same local constraint.
Constraint satisfaction problems, due to their analogies with disordered
systems, have been extensively studied using methods from statistical physics.
The existence and entropy of solutions can be obtained by studying the uniform
measure over the solutions. The partition function of this measure, or its
proxy the free entropy, then counts the number of solutions. The computation
of the partition function is hard in general, but approximations of the free
entropy can be computed using the replica symmetric cavity method [7, 8, 9,
10], which is particularly suited to sparse graphs. Refining the replica
symmetric approach, one can introduce assumptions of replica symmetry breaking
within the cavity formalism [11, 12, 13]. These advancements lead to the
development of new efficient algorithms [14, 15]. The statistical physics
approach also allows the study of phase transitions in the space of solutions,
of which the clustering [16, 17, 18] and the freezing [19, 20] transitions are
of particular interest due to their conjectured link to the typical
computational algorithmic complexity of the problem [21, 22].

In this paper, we systematically investigate the class of outer-totalistic
CSPs on random $3-$regular graphs using the cavity method, with the aim to
discover novel interesting problems and structural similarities between them.
In Section 1, we define the outer-totalistic CSPs and provide a brief overview
of some outer-totalistic CSPs studied in the literature. In Section 2, we
derive the RS cavity fixed point equation, which we use to estimate the
entropy of solutions. In Section 3, we detail the classification based on the
entropy of solutions. Section 4 introduces the general one-step replica
symmetry breaking formalism, including an overview of the different phases,
and the 1RSB cavity equations. The presence of replica symmetry breaking and
its relation to the structure of the space of solutions is discussed in
Section 5. Finally, the connection between the phase and the experimental
computational hardness of the outer-totalistic CSPs is explored in Section 6.

The code to generate the data and plot the figures can be found at
https://github.com/SPOC-group/fixed_points_graph_CA.

##  1 Background and Notation

#### Outer-totalistic constraint satisfaction problems.

Consider an undirected $d$-regular graph $G$ with $N$ nodes
$i=\\{1,\ldots,n\\}$. We denote the set of edges $E(G)$ and the edge between
node $i$ and $j$ as $ij$. Each node $i$ is in a state $s_{i}\in S$ from a
discrete alphabet $S$. In this work $S=\\{0,1\\}$, and we denote node $i$ as
empty if $s_{i}=0$ and occupied if $s_{i}=1$. We refer to
$\mathbf{s}=(s_{1},...,s_{N})$ as a configuration. A configuration is said to
be homogeneous if all nodes are in the same state, i.e.
$s_{i}=s_{j}\,\,\forall i,j=1,\ldots,N$.

A configuration is said to be a solution of a CSP if all its nodes respect all
the given constraints. In our case, every node of the graph must respect the
same local constraint. We call this constraint a rule. Let us denote as
$f(s_{i},\\{s_{j}\\}_{j\in\partial i})$ the indicator function that is $1$ if
node $i$ with value $s_{i}$ respects the rule when its direct neighbors
$\partial i$ take values $\\{s_{j}\\}_{j\in\partial i}$, and $0$ if it does
not. We study outer-totalistic rules, rules that depend only on the value of a
node and the total number of occupied nearest neighbors. In the case where
$S=\\{0,1\\}$, we have that $f(s_{i},\\{s_{j}\\}_{j\in\partial
i})=f(s_{i},\sum_{j\in\partial i}s_{j})$, but we will keep the first notation
for generality.

To denote an outer-totalistic rule, we extend the notation introduced by [23].
This notation consists of an ordered list of $d+1$ symbols, indicating what
value a node can take to satisfy the constraint for $0,1,\ldots,d$ occupied
neighbors. The symbols are

  * 0

: the node must be empty to satisfy the constraint,

  * 1

: the node must be occupied to satisfy the constraint,

  * +

: the constraint is satisfied regardless of the value of the node,

  * -

: the constraint cannot be satisfied.

As an example, consider rule -+10: the constraint is not satisfied if a node
has $0$ occupied neighbors, the node can take any value for $1$ occupied
neighbor, it must be occupied for $2$ occupied neighbors, and it must be empty
if it has $3$ occupied neighbors. Note that inverting the empty and occupied
states gives rise to equivalent rules. To avoid studying twice the same
problem, we remove this equivalence with the procedure explained in Appendix
A. For example, we have only $136$ non-equivalent rules out of $4^{4}=256$
rules for $d=3$. This notation was originally used to describe the dynamics of
outer-totalistic cellular automata, but we adopt it here to describe static
CSPs. While we use the cavity method to study the static case, note that it
can also be used to study dynamical systems, for instance in [24] where the
same rule notation as here is used.

Note that this notation introduces a hierarchy between the rules: starting
from rule ++++ for which every configuration is a solution, one can add
constraints, first turning + into 0 or 1 and then into -. This gives a simple
criterion when comparing the number of solutions of a given rule: all the
rules with more constraints should have at most the same number of solutions.

#### Number of stationary configurations.

Our aim is to count the number of solutions to a given outer-totalistic rule.
The uniform probability measure over solutions of an outer-totalistic rule on
a given $d$-regular graph $G$ is

$P_{G}(\mathbf{s})=\frac{\prod_{i=1}^{N}f(s_{i},\\{s_{j}\\}_{j\in\partial
i})}{Z_{G}},$ (1)

where $f(s_{i},\\{s_{j}\\}_{j\in\partial i})$ indicates whether the rule is
satisfied by taking the value $1$, and taking the value $0$ otherwise. $Z_{G}$
is the normalization constant (i.e. the partition function) and is given by

$Z_{G}=\sum_{\mathbf{s}\in\mathcal{C}}\prod_{i=1}^{N}f(s_{i},\\{s_{j}\\}_{j\in\partial
i})$ (2)

where $\sum_{\mathbf{s}\in\mathcal{C}}$ indicates a sum over all possible
configurations $\mathcal{C}=\\{0,1\\}^{N}$. Then, $Z$ is the number of
stationary configurations. We define the free entropy density $\phi_{G}$ on a
given $d$-regular graph $G$ as

$\phi_{G}=\frac{\log{Z_{G}}}{N}.$ (3)

These quantities can be analyzed on $d$-regular graphs in the large size
limit. The cavity method (see Section 2) allows us to compute the quenched
free entropy density $\phi$ defined as

$\phi=\lim_{N\rightarrow\infty}\frac{1}{N}\mathbb{E}_{G}[\log Z_{G}]$ (4)

where $\mathbb{E}_{G}$ denotes the average over the ensemble of random
$d$-regular graphs of $N$ nodes. In statistical physics, it is standard that
the quenched free entropy density, as well as other averaged quantities of
interest, are self-averaging, i.e. concentrate around their mean. In
particular, this means that the free entropy density (3) for a typical large
$d$-regular graph $G$ is equal to the quenched free entropy density (4) with
high probability. It is instructive to also define the annealed free entropy
density $\phi_{\text{ann}}$

$\phi_{\text{ann}}=\lim_{N\rightarrow\infty}\frac{1}{N}\log\left(\mathbb{E}_{G}[Z_{G}]\right)$
(5)

that provides an upper bound to the quenched free entropy density.

#### Relation to other constraint satisfaction problems.

Multiple well-studied CSPs are encompassed within the introduced rule
notation. We present a non-exhaustive list of them, along with relevant
literature that studied them using statistical physics approaches, in
particular the cavity method.

The independent set or its complement, the vertex cover, was studied using the
cavity method in [25, 26]. The corresponding rules are of the form
+0$\ldots$0. To study the maximum independent set, a Lagrange multiplier $\mu$
can be introduced and taken to infinity.

The outer-totalistic CSPs are related to the lattice glass model proposed in
[27] and further studied by Rivoire et al. [28]. In this model, occupied nodes
can have at most $k$ occupied neighbors. This corresponds to rule
+$\ldots$+0$\ldots$0. They obtained an analytic expression for the densest
packing ($\mu\rightarrow\infty$) within the 1RSB approximation.

Another problem that can be recast in the language of outer-totalistic CSPs is
the maximal independent set problem (rule 10$\ldots$0), studied with the
cavity method in [29], or the dominating set problem (rule 1+$\ldots$+),
studied with the cavity method in [30].

The assortative/dissasortative partition problems were studied in [31]. The
corresponding rules are of the form 0$\ldots$0$\pm$$\ldots$$\pm$1$\ldots$1 and
1$\ldots$1$\pm$$\ldots$$\pm$0$\ldots$0. These problems were studied for
balanced partitions, where the number of occupied and empty nodes are the
same. This balance can be imposed by choosing an appropriate value of the
parameter $\mu$ or by forcing the $0\leftrightarrow 1$ symmetry of the
messages.

Occupation problems on hyper-graphs were studied in [20], and interestingly
the averaged entropy of their solutions correspond to outer-totalistic rules
composed only of + and -. More specifically, the free entropy obtained in the
replica symmetric case is the same even if the underlying graphical model is
not, and there is no clear mapping from the solutions of one to the solutions
of the other on a single instance. The BP equation in present outer-totalistic
CSPs seem to reduce to the BP equation in occupation problems from [20], but
we were only able to show this by assuming a structure in our messages that we
do not know how to justify. The correspondence was verified numerically for
$d=3,4$. See Appendix F for a more detailed discussion.

##  2 Replica Symmetric Solution using Belief Propagation

The free entropy density for a given graph (3) and the quenched free entropy
density (4) are in general challenging to compute. Both can be estimated on
large random graphs using the replica symmetric cavity method [9, 10]. The
approximation on a given graph $G$ is known as the Bethe-Peierls
approximation, and yields the so-called Bethe free entropy. Averaging the
Bethe free entropy over the ensemble of $d$-regular graphs gives the replica
symmetric approximation of the quenched free entropy density (4). If this
approximation is exact, we say that the problem is replica symmetric (RS). The
cavity method yields a fixed-point equation, which can be iterated as an
algorithm. This algorithm is known as belief propagation (BP), message passing
or sum-product. We use the general method of casting the problem into a
graphical model to obtain the BP equations.

#### Graphical model.

Figure 1: Left: original graph and naive factor graph. The edges forming a
short loop are drawn in red. Right: original graph and factor graph. Notice
that the nodes become the factor nodes and the edges the variable nodes.

A graphical model is the representation of a factorized probability
distribution, such as eq. (1) on a factor graph. A factor graph is a graph
with two types of nodes: variable nodes, one for each variable, and factor
nodes, one for each interaction (or constraint in our case). Variable nodes
are usually drawn as circles and factor nodes as squares. An edge is present
between a variable node and a factor node if its associated constraint depends
on the value of the variable. The factor graph is thus bipartite: edges exist
only between variable nodes and factor nodes.

The naive way of constructing the graphical model of the probability
distribution (1) is to have one variable for each state $s_{i}$, and a factor
node that checks if the constraint on $s_{i}$ is satisfied. We say that the
graph $G$ on which the problem is set is the original graph (or simply the
graph when clear from the context). An example of an original graph and
associated factor graph is illustrated in Figure 1 (left). We notice that,
even if the original graph is a tree, the factor graph contains short loops
(in red). This can be problematic, because the belief propagation algorithm
presented below is heuristically known to work better on factor graphs with
large loops [32]. Therefore, we use an alternative graphical model.

Consider the variables to be on the edges and the factor nodes to be on the
nodes of the original graph. This means that the variables are tuples
$s_{i},s_{j}$. We denote the variables nodes as $ij$ and the factor nodes as
$i$ (since they are equivalent to the nodes of the original graph). The
variable $s_{i},s_{j}$ take value in $S\times S$. For $S=\\{0,1\\}$, this
gives 4 different possible values. By construction, the factor nodes have an
edge towards each variable needed to check the constraint. An example of the
original graph and obtained factor graph is presented in Figure 1 (right).
Notice that the factor graph has exactly the same structure as the original
graph. In particular, the loops are the same as in the original graph. Similar
constructions are done e.g. in [29, 31, 24].

#### Belief propagation on a given graph.

From the graphical model, we can deduce a fixed point equation named the
replica symmetric cavity or belief propagation (BP) equation. The equation is
derived in detail in Appendix B for the case of trees, where the Bethe-Peierls
approximation is exact and yields the true free entropy density (7). A broader
condition for BP to be exact is the independence of the incoming messages for
each variable node. This is evidently the case if the factor graph is a tree.
However, one also heuristically expects locally tree-like factor graphs, such
as typical $d$-regular graphs for large $N$ and constant $d$, to give good
results if the correlations decay fast enough in the different loops. This
motivates the construction of a graphical model with no short loops. The
equation derived for tree graphs can be used heuristically on general graphs
and is given by

$\psi^{i\rightarrow j}_{s_{i},s_{j}}=\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus
j}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus j})\prod_{k\in\partial
i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow i}.$ (6)

where we write $\\{s_{k}\\}_{k\in\partial i}$ as
$s_{j},\\{s_{k}\\}_{k\in\partial i\setminus j}$ to emphasise that $s_{j}$ is
fixed by $\psi^{i\rightarrow j}_{s_{i},s_{j}}$. We denote a configuration of
the neighbors of node $i$ with node $j$ removed as $\\{s_{k}\\}_{k\in\partial
i\setminus j}$, and the sum is over all these possible configurations.
$\psi^{i\rightarrow j}_{s_{i},s_{j}}$ is called the message from node $i$ to
$j$, and can be used to compute the free entropy density, and thus also other
quantities of interest. $Z^{i\rightarrow j}$ is the normalization such that
$\sum_{s_{i}}\sum_{s_{j}}\psi^{i\rightarrow j}_{s_{i},s_{j}}=1$. The Bethe
free entropy density $\phi_{\text{Bethe}}$ for a given graph $G$, also derived
in Appendix B, is given by

$\phi_{\text{Bethe}}=\frac{1}{N}\sum_{i=1}^{N}\log
Z^{i}-\frac{1}{N}\sum_{ij\in E(G)}\log Z^{ij}$ (7)

with

$Z^{i}=\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i},$ (8)
$Z^{ij}=\sum_{s_{i}}\sum_{s_{j}}\psi_{s_{(ij)}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i}.$ (9)

#### Belief propagation on large random $d$-regular graphs.

In this work, we study the outer-totalistic CSPs on large random $d$-regular
graphs. Under the assumption that all the nodes are locally equivalent, which
is motivated by the tree-like structure of random regular graphs, the fixed
point equation in the $N\rightarrow\infty$ limit reduces to 4 self-consistent
equations:

$\psi_{s_{i},s_{j}}=\frac{1}{Z_{\psi}}\sum_{\\{s_{k}\\}_{k=1}^{d-1}}f(s_{i},s_{j},\\{s_{k}\\}_{k=1}^{d-1})\prod_{k=1}^{d-1}\psi_{s_{k},s_{i}}.$
(10)

Recall that $f(s_{i},s_{j},\\{s_{k}\\}_{k=1}^{d-1})$ is $1$ if the constraint
is satisfied and $0$ otherwise. $Z_{\psi}$ is the normalization and is given
by

$Z_{\psi}=\sum_{s_{i}}\sum_{s_{j}}\sum_{\\{s_{k}\\}_{k=1}^{d-1}}f(s_{i},s_{j},\\{s_{k}\\}_{k=1}^{d-1})\prod_{k=1}^{d-1}\psi_{s_{k},s_{i}}.$
(11)

The iteration scheme to obtain a fixed point is presented in Appendix B.1. As
before, the messages allow us to compute what we call the replica symmetric
free entropy density $\phi_{\rm RS}$:

$\phi_{\rm
RS}=\log\left(Z_{\text{n}}\right)-\frac{d}{2}\log\left(Z_{\text{e}}\right),$
(12)

where

$Z_{\text{n}}=\sum_{s_{i},\\{s_{k}\\}_{k=1}^{d}}f(s_{i},\\{s_{k}\\}_{k=1}^{d})\prod_{k=1}^{d}\psi_{s_{k},s_{i}},$
(13) $Z_{\text{e}}=\sum_{s_{i},s_{j}}\psi_{s_{i},s_{j}}\psi_{s_{j},s_{i}}.$
(14)

Note that in this setting the graph is no longer specified. Moreover, in the
specific case of regular graphs, it can be shown that the replica symmetric
free entropy is equal to the annealed free entropy (5) [33] so that we have
$\phi_{\text{ann}}=\phi_{\rm RS}$. In the following sections we will
investigate for which rules we also have equality $\phi_{\rm RS}=\phi$ with
the quenched free energy defined in (4).

#### Stability and validity of the RS solution.

A necessary (but not sufficient) condition for the validity of the belief
propagation and the resulting free entropy is the so-called stability [34,
35]. Stability is equivalent to the convergence of the BP equations on a
single large graph. This can be verified numerically using a population
dynamics algorithm. The algorithm is presented in Appendix D. We investigate
the stability of each rule by initializing the population dynamics with the BP
messages with added noise. We use Gaussian noise with standard deviation
$10^{-6}$, $\psi_{s_{i},s_{j}}^{\rm
noisy}=\max(\psi_{s_{i},s_{j}}+\eta,0),\,\eta\sim\mathcal{N}(0,10^{-12})$,
which we then normalize. If the population dynamics goes back to the BP fixed
point, the fixed point is stable [35].

If the fixed point is not stable, then the replica symmetric assumption is
incorrect and the problem presents replica symmetry breaking. In that case, we
analyse the problem within the one-step replica symmetry breaking (1RSB)
assumption, which is presented in Section 4. Note that the reciprocal is not
necessarily true: the replica symmetric assumption can be incorrect even if
the RS fixed point is stable.

We find that the RS solution is stable for all rules with existing solutions
($\phi\geq 0$), except for two rules: +-00 and +-0+. As will be seen in
Section 5, these rules indeed present replica symmetry breaking. More
generally, the validity of the RS solution can also be checked by doing the
1RSB calculation: if only the trivial 1RSB solutions exist (i.e. the only 1RSB
solution is the replica symmetric one), then the RS solution is correct. This
was verified for all rules marked ‘RS’ in Tables 2 and 3.

##  3 Classification of Outer-totalistic Rules for $d=3$

We classify all possible non-equivalent outer totalistic rules depending on
the number of solutions. To compute the RS fixed point, we start the BP
iteration (see Appendix B.1) with $10^{6}$ different initializations
$\psi_{s_{i},s_{j}}=0,0.01,0.02,...,1\,\forall\,s_{i},s_{j}\in S$. It is
possible that multiple fixed points exist, and that they are found numerically
depending on the initialization. In that case, the fixed point with the
highest free entropy is kept. In the following section, we present a
classification based on the number and type of solutions.

The first clear distinction is whether a rule admits homogeneous solutions or
not. The formula of the rule clearly indicates if the rule admits a
homogeneous solution. In this case, there are three possible scenarios:

  * •

Only homogeneous solution: The problem only admits homogeneous solutions. This
means that $\phi=0$ and that the fixed point messages are only compatible with
homogeneous solutions.

  * •

Subexponentially many solutions: Other solutions than the homogeneous ones
exist, but their number is subexponential in $N$, i.e. $\phi=0$.

  * •

Exponentially many solutions: An exponential number of solutions exist
($\phi>0$) in addition to the homogeneous solutions.

Table 1 and 2 list all the rules that have homogeneous solutions, Table 1
listing the rules that have only homogeneous solutions. Table 3 presents all
the rules that have no homogeneous solutions.

In the latter case, the following possibilities arise:

  * •

Locally contradictory: It is not possible to find a solution on a $d$-regular
tree of depth 2 for any boundary condition.

  * •

No solution: The rule is not locally contradictory but the free entropy
density is negative.

  * •

Subexponentially many solutions: The free entropy density is $0$, there exists
a subexponential number of solutions.

  * •

Exponentially many solutions: The free entropy density is strictly positive,
there is an exponentially number solutions.

For the rules that are locally contradictory, either the message is not
normalisable (all the entries of the messages are $0$) or there is no
convergence. In this case, we define the free entropy density to be $-\infty$
and the phase as unsatisfiable (UNSAT) as there are no solutions. The rules
with no solutions but that are not locally contradictory admit a fixed point
with negative free entropy density. Thus, the obtained fixed point allows us
to distinguish between the rules that are locally contradictory and the rules
that have no solutions but are not locally contradictory.

The classification in Tables 1, 2 and 3 are done using the replica symmetric
free entropy, which is not necessarily the same as the true free entropy. To
study the validity of the replica symmetric free entropy and correct it if
necessary, we compute the free entropy within the one-step replica symmetry
breaking framework which is presented below in Sections 4 and 5. As we will
see, for $d=3$, most rules are replica symmetric. The inclusion of one-step
replica symmetry breaking does not change the classification.

##  4 One-step Replica Symmetry Breaking

In this section, we give a brief overview of one-step replica symmetry
breaking viewed from the lens of the cavity method. The method was introduced
in [12, 13], and for more details we refer the reader to [10, 32].

#### Clustering.

When the replica symmetric assumption holds, almost all stationary solutions
are said to be in a cluster or pure state in configuration space. When the
replica symmetric assumption is not correct, we say that there is a replica
symmetry breaking. If the configuration space is correctly described by
decomposing the replica symmetric cluster into smaller clusters, we say that
there is a one-step replica symmetry breaking (1RSB). For a given graph $G$
there are two messages $\psi^{i\rightarrow j}$ and $\psi^{j\rightarrow i}$ for
each edge $ij\in E(G)$. We identify each cluster by the fixed point
$\tilde{\psi}=\\{\psi^{i\rightarrow j},\psi^{j\rightarrow i}\\}_{ij\in E(G)}$.
The probability measure $\nu$ over the cluster associated to a fixed point
$\tilde{\psi}$ is

$\nu(\tilde{\psi})=\frac{1}{Z_{\nu}}e^{Nm\phi_{\rm int}(\tilde{\psi})}$ (15)

where $m$ is the Parisi parameter, $Z_{\nu}$ is the normalization and
$\phi_{\rm int}(\tilde{\psi})$ is the Bethe free entropy density associated to
the fixed point $\tilde{\psi}$. We will call $\phi_{\rm int}$ the internal
free entropy density. The free entropy density associated with the measure
$\nu$ is called the replicated free entropy $\Psi(m)=\log(Z_{\nu})/N$. The
complexity or configurational entropy density $\Sigma$ is the entropy density
at the level of the clusters, i.e. $\mathcal{N}_{c}(\phi)=e^{N\Sigma(\phi_{\rm
int})}$ where $\mathcal{N}_{c}(\phi_{\rm int})$ is the number of cluster with
internal entropy $\phi_{\rm int}$. Then we have that

$e^{N\Psi(m)}=\sum_{\phi_{\rm int}}e^{N(m\phi_{\rm int}+\Sigma(\phi_{\rm
int}))}$ (16)

where the sum is over all the possible different internal free entropy
densities. Applying the Laplace method we obtain

$\Psi(m)=m\phi_{\rm int}+\Sigma(\phi_{\rm int})$ (17)

where $\phi_{\rm int}$ extremizes $\Psi(m)$. Thus,
$\partial\Sigma/\partial\phi_{\rm int}=-m$ from the properties of the Legendre
transform. To obtain the internal free entropy $\hat{\phi}_{\rm int}$ of the
clusters containing the typical solutions (i.e. the solutions that dominate
the measure 1), we maximize the total free entropy $\phi_{\rm
int}+\Sigma(\phi_{\rm int})$ under the constraint $\Sigma(\phi_{\rm int})\geq
0$. Indeed, assuming that solutions exist, there must be at least one cluster.
If we first ignore this constraint, the maximal total entropy is obtained for
$m=1$. If the complexity is positive, we call such a clustered phase the
dynamical one-step replica symmetry breaking phase (d1RSB). In that case,
there exists an exponential number of clusters with internal free entropy
$\hat{\phi}_{\rm int}$. However, if the corresponding complexity is negative
for $m=1$, then the total free entropy is maximized for the largest internal
free entropy density $\hat{\phi^{\prime}}_{\rm int}$ such that
$\Sigma(\hat{\phi^{\prime}}_{\rm int})\geq 0$. We write
$\hat{\phi^{\prime}}_{\rm int}=\phi_{s}$. In that case we are in the static
one-step replica symmetry breaking phase (s1RSB), there are subexponentially
many clusters with internal free entropy $\phi_{s}$. Since the RS free entropy
is equal to the annealed free entropy in the case of regular graphs [33], the
s1RSB phase indicates that the annealed free entropy is different from the
quenched free entropy. If no $\hat{\phi^{\prime}}_{\rm int}$ exist such that
$\Sigma(\hat{\phi^{\prime}}_{\rm int})\geq 0$, then the problem is said to be
unsatisfiable (note that the problem can of course also be unsatisfiable
without replica symmetry breaking). Figure 2 (left) shows a sketch of the
space of solutions for the different phases.

#### Freezing.

A variable $s_{i}$ is said to be frozen if it takes the same value for every
configuration within a cluster. If an extensive number of variables are frozen
within a cluster, we say that the cluster is frozen. Using the usual
nomenclature (summarized for example in [36]), the rigidity transition happens
when frozen variables appear in typical solutions, and such a phase is called
rigid. The transition where all clusters become frozen is called the freezing
transition, and such a phase is called frozen. If almost all the variables
within a cluster are frozen, the cluster is said to be locked. In that case
the internal free entropy of the cluster is $0$.

#### 1RSB cavity solution.

To estimate $Z_{\nu}$, we write similar equations as BP for the probability
distribution $\nu$. Again assuming that the graph is $d$-regular and that the
messages are locally equivalent we obtain

$P(\psi)=\frac{1}{Z_{P}}\int\prod_{k=1,...,d-1}dP(\psi^{(k)})\left(\frac{Z_{\text{n}}}{Z_{\text{e}}}\right)^{m}\delta\left(\psi-\mathcal{F}(\\{\psi^{(k)}\\}_{k=1,...,d-1})\right).$
(18)

where $P$ is a continuous probability distribution over the messages. The
integral is over all possible $d-1$ messages that are normalized. $\psi$ and
$\psi^{(k)}$ indicate the 4-component message introduced in eq. (10).
$\mathcal{F}\left(\\{\psi^{(k)}\\}_{k=1,...,d-1}\right)$ indicates the message
$\psi$ obtained from eq. (10) with the substitution
$\psi_{s_{k},s_{i}}=\psi^{(k)}_{s_{k},s_{i}}$. $Z_{P}$ is the normalization
and $Z_{\text{n}}$ and $Z_{\text{e}}$ are defined in eq. (13) and (14). For
notational purposes we do not explicitly write the dependence of
$Z_{\text{n}}$ and $Z_{\text{e}}$ on the message. Supposing that the 1RSB
solution holds, the replicated free entropy is given by

$\Psi_{\rm
1RSB}(m)=\log\left(\mathcal{Z}_{n}\right)-\frac{d}{2}\log\left(\mathcal{Z}_{e}\right)$
(19)

where

$\mathcal{Z}_{n}=\int\prod_{k=1,...,d-1}dP(\psi^{(k)})(Z_{\text{n}})^{m},$
(20) $\mathcal{Z}_{e}=\int dP(\psi^{(1)})dP(\psi^{(2)})(Z_{\text{e}})^{m}.$
(21)

The internal entropy density can be obtained from

$\phi_{\rm int}=\frac{\partial\Psi(m)}{\partial m}$ (22)

which allows to compute the complexity using (17). A detailed derivation of
eq. (18), (19) and (22) is given in Appendix C. Since the fixed point equation
is now over a continuous distribution, we must usually approximate it
numerically. For this, we use the population dynamics algorithm described in
Appendix D. The fraction of frozen variables can be determined by counting the
number of messages in the population that imply that a node has to take a
fixed value.

##  5 Structure of the Space of Solutions

Figure 2: Left: Sketch of the space of configurations for the different
phases. Blue indicate typical solutions. For the 1RSB phases the complexity vs
internal entropy curves are sketch on the top left, the blue doted line
indicates the tangent line where the slope of the complexity curve is $-1$,
which is where $m=1$. Right: Complexity as function of the internal free
entropy to estimate $\phi_{s}$ for the s1RSB rule +-00.

#### Population dynamics and phases.

As mentioned in Section 4, we estimate the solution to the 1RSB equation (18)
using a population dynamics algorithm (see Appendix D). This also allows one
to determine the phase. The population dynamics is conducted with a population
of $M=10^{7}$ messages initialized on hard-fields (i.e. each message of the
population dynamics is $1$ for one of the variables $(s_{i},s_{j})$, and $0$
on the others, which fixes the state of the nodes) proportionally to the
obtained message from BP, e.g. if $\psi_{0,0}=0.3$ then we initialize $0.3M$
messages with $\psi_{0,0}=1$. The initialization on hard-fields is justified
by the problem of reconstruction on trees. Indeed, the existence of non-
trivial solutions is linked to the reconstruction on trees property [37]:
consider a tree with a valid configuration, and remove all information about
the state of the nodes except on the leaves. If the leaves contain some
information about the value of the root, then the reconstruction on trees is
possible. The initial population plays the role of the leaves and is thus
initialized on hard-fields.

Another initialization of interest corresponds to the small noise
reconstruction on trees. It is similar to the initialization on hard-fields,
but with a small noise added on the messages. It can happen that both
initializations do not give the same results, in particular that the hard-
fields initialization yields a non-trivial solution whereas the small noise
initialization converges to the RS fixed point. In that case, the solutions
are ‘fake clustered’ in the same sense as was proposed in [20] where a similar
situation arises. In these cases, we posit that the RS solution in fact
correctly describes the problem. An illustrative case of the phenomenon can be
seen for rule 0010. This rule admits as solutions all configurations that have
non-intersecting loops. Thus, we can go from one solution to another by
removing or adding one loop, which can be of subextensive length. This
phenomenon seems to happen particularly for graphs with a small average
degree. The ‘phase’ of the rules were the hard-fields initialization has a
non-trivial solution but the small noise initialization converges to the RS
fixed point are noted RS*, and the complexity and internal free entropy are
those obtained from the hard-field initialization.

The other phases are described in Sections 3, 4 and are

  * •

UNSAT: No RS fixed point or contradictory RS fixed point

  * •

RS: Replica symmetric

  * •

d1RSB: Dynamical one-step replica symmetry breaking

  * •

s1RSB: Static one-step replica symmetry breaking

They are shown in both Tables 2 and 3. The d1RSB phase is identified when the
population dynamics solution at Parisi parameter $m=1$ is not trivial (i.e.
not the same as the RS solution) and both the internal free entropy density
and complexity are positive (see Figure 2 left for a sketch). The dynamical
1RSB phase holds the same entropy as the RS solution [35]. Frozen variables
were found in all identified d1RSB rules. In brackets the standard deviation
of the sampled observables at different iterations (see Appendix D) is
indicated. This quantity is not indicated if it is smaller than $10^{-5}$.

Rules with only homogeneous solutions.
---
0000, 000+, 000-, 0001, 00-0, 00-+, 00--, 00-1, 0-00, 0-0+, 0-0-, 0-01, 0-+0,
0-+-, 0--0,
0--+, 0---, 0--1, 0-10, 0-1-, 010-, 01-0, 01-+, 01--, 011-, +-0-, +--0, +--+,
+---, +-10,
+-1-, +1-0, +1--, +11-.
Table 1: All the non-equivalent rules with $d=3$ neighbors that only have
homogeneous solutions.

Subexponentially many solutions
---
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$ | $\phi_{s}$
0101 | RS | 0.0000 |  |  |
+-+- | d1RSB | 0.0000 | 0.0000 | 0.0000 |
Exponentially many solutions
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$ | $\phi_{s}$
00+0 | RS* | 0.2046 | 0.0000 | 0.2046(1) |
00++ | RS | 0.4133 |  |  |
00+- | RS* | 0.1925 | 0.0000 | 0.1925(2) |
00+1 | RS | 0.3793 |  |  |
0010 | RS* | 0.1204 | 0.0000 | 0.1203(1) |
001+ | RS | 0.2833 |  |  |
001- | RS* | 0.1160 | 0.0000 | 0.1160(1) |
0011 | RS | 0.2350 |  |  |
0+00 | RS | 0.3671 |  |  |
0+0+ | RS | 0.3967 |  |  |
0+0- | RS | 0.3484 |  |  |
0+01 | RS | 0.3740 |  |  |
0++0 | RS | 0.5548 |  |  |
0+++ | RS | 0.6350 |  |  |
0++- | RS | 0.5134 |  |  |
0++1 | RS | 0.5829 |  |  |
0+-0 | RS | 0.2443 |  |  |
0+-+ | RS | 0.2601 |  |  |
0+-- | RS | 0.2359 |  |  |
0+-1 | RS | 0.2497 |  |  |
0+10 | RS | 0.3603 |  |  |
0+1+ | RS | 0.4265 |  |  |
0+1- | RS | 0.3348 |  |  |
0-++ | RS | 0.3456 |  |  |
0-1+ | RS | 0.2833 |  |  |
0100 | d1RSB | 0.0416 | 0.0107(2) | 0.0309(2) |
010+ | d1RSB | 0.1303 | 0.0071(2) | 0.1232(2) |
01+0 | RS | 0.3715 |  |  |
(a)

01++ | RS | 0.5219 |  |  |
---|---|---|---|---|---
01+- | RS | 0.2107 |  |  |
0110 | RS | 0.1893 |  |  |
011+ | RS | 0.4131 |  |  |
+000 | RS | 0.4354 |  |  |
+00+ | RS | 0.4354 |  |  |
+00- | RS | 0.4131 |  |  |
+0+0 | RS | 0.4354 |  |  |
+0++ | RS | 0.4616 |  |  |
+0+- | RS | 0.4131 |  |  |
+0-0 | RS | 0.2938 |  |  |
+0-+ | RS | 0.2938 |  |  |
+0-- | RS | 0.2833 |  |  |
+010 | RS | 0.2938 |  |  |
+01+ | RS | 0.2964 |  |  |
+01- | RS | 0.2833 |  |  |
++00 | RS | 0.5513 |  |  |
++0+ | RS | 0.5605 |  |  |
++0- | RS | 0.5143 |  |  |
+++0 | RS | 0.6489 |  |  |
++++ | RS | 0.6931 |  |  |
+++- | RS | 0.5975 |  |  |
++-0 | RS | 0.3610 |  |  |
++-+ | RS | 0.3662 |  |  |
++-- | RS | 0.3415 |  |  |
++10 | RS | 0.4233 |  |  |
++1- | RS | 0.3900 |  |  |
+-00 | s1RSB | 0.0448 | 0.0491(2) | -0.0044(1) | 0.0445
+-0+ | s1RSB | 0.0448 | 0.0491(1) | -0.0044(1) | 0.0445
+-+0 | RS* | 0.1304 | 0.0259(3) | 0.1045(2) |
+100 | RS | 0.2681 |  |  |
+10+ | RS | 0.2964 |  |  |
+10- | d1RSB | 0.1038 | 0.0000 | 0.1038(2) |
+1+0 | RS | 0.4584 |  |  |
+1+- | RS | 0.2782 |  |  |
+110 | RS | 0.2617 |  |  |
(b)

Table 2: Classification of all non-equivalent outer-totalistic CSPs with
homogeneous solutions for 3-regular graphs. The rules with only homogeneous
solutions are indicated in Table 1. RS* indicates that the hard-fields
initialization yields a non-trivial solution while the small-noise
initialization does not. The internal free entropy $\phi_{\rm int}$ and
complexity $\Sigma$ and are those obtained from the hard-field initialization
at Parisi parameter $m=1$. $\phi_{s}$ indicates the free entropy for the rules
presenting a static 1RSB phase.

#### Static 1RSB.

The RS estimation of the free entropy density is wrong in the case of s1RSB.
The rules +-00 and +-0+ are in the s1RSB phase. In that case, the static one-
step replica symmetry free entropy density $\phi_{s}$ is found by computing
the complexity as a function of the internal free entropy. This is done by
varying the Parisi parameter between $0$ and $1$ (in our case we used 30
values of $m$), and fitting the obtained points with a function of the form
$\Sigma(\phi_{\rm int})=a+b2^{\phi_{\rm int}}+c3^{\phi_{\rm int}}$ as done in
[19]. The intersection of this curve with $\Sigma=0$ gives the value of
$\phi_{s}$. Figure 2 (right) shows the complexity versus internal free entropy
curve for rule +-00. We also recall that the RS solution of rule +-00 and +-0+
is unstable (see the discussion on stability in Section 2), indicating the
presence of replica symmetry breaking. We recall that the annealed free
entropy is equal to the RS free entropy in the case of regular graphs [33].
Thus, rules +-00 and +-0+ are the only rules in the $d=3$ case where according
to our study the quenched free entropy is not equal to the annealed one.

The classification introduced in Section 3 does not change by considering the
1RSB solution. Thus, at least for $d=3$, the RS solution is sufficient to
determine if solutions exist.

Locally contradictory
---
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$
-000 | UNSAT | $-\infty$ |  |
-00- | UNSAT | $-\infty$ |  |
-0-0 | UNSAT | $-\infty$ |  |
-0-- | UNSAT | $-\infty$ |  |
\--00 | UNSAT | $-\infty$ |  |
\--0- | UNSAT | $-\infty$ |  |
\---0 | UNSAT | $-\infty$ |  |
\---- | UNSAT | $-\infty$ |  |
No solutions but not locally contradictory
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$
-010 | UNSAT | -0.1116 | 0.0000 | -0.1116(2)
-01- | UNSAT | -0.2027 | 0.0000 | -0.2027(2)
-+-0 | UNSAT | -0.0774 | 0.0000 | -0.0774(3)
-+-- | UNSAT | -0.1744 | 0.0000 | -0.1745(3)
\--+0 | UNSAT | -0.0303 | 0.0000 | -0.0308(2)
\--10 | UNSAT | -0.1733 | 0.0000 | -0.1733(2)
-100 | UNSAT | -0.0012 | 0.0000 | -0.0012(2)
-10- | UNSAT | -0.2027 | 0.0000 | -0.2027(3)
-1-0 | UNSAT | -0.2502 | 0.0000 | -0.2502(2)
10-0 | UNSAT | -0.0760 | 0.0000 | -0.0760(2)
1-00 | UNSAT | -0.0189 | 0.0000 | -0.0189(1)
1--0 | UNSAT | -0.3466 | 0.0000 | -0.3466(2)
(a)

Subexponentially many solutions
---
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$
1010 | d1RSB | 0.0000 | 0.0000 | 0.0000(2)
Exponentially many solutions
rule | phase | $\phi_{\rm RS}$ | $\phi_{\rm int},m=1$ | $\Sigma,m=1$
-0+0 | RS* | 0.1016 | 0.0000 | 0.1016(2)
-0+- | RS* | 0.0523 | 0.0000 | 0.0523(2)
-+00 | RS | 0.2191 |  |
-+0- | RS | 0.1541 |  |
-++0 | RS | 0.4880 |  |
-++- | RS | 0.4055 |  |
-+10 | RS* | 0.2149 | 0.0817(3) | 0.1331(3)
-1+0 | RS | 0.3551 |  |
-110 | RS | 0.1893 |  |
1000 | RS | 0.2617 |  |
10+0 | RS* | 0.2845 | 0.1460(4) | 0.1384(4)
1+00 | RS | 0.4414 |  |
1++0 | RS | 0.5829 |  |
1+-0 | d1RSB | 0.1038 | 0.0000 | 0.1038(2)
1100 | RS | 0.2350 |  |
(b)

Table 3: Classification of all non-equivalent outer-totalistic CSPs without
homogeneous solutions for 3-regular graphs. RS* indicates that the hard-fields
initialization yields a non-trivial solution while the small-noise
initialization does not. The internal free entropy $\phi_{\rm int}$ and
complexity $\Sigma$ and are those obtained from the hard-field initialization
at Parisi parameter $m=1$. In the case of rules with no homogeneous solutions
for $d=3$, no rule presents a static 1RSB phase.

#### Stability of the 1RSB solution

Checking the validity of the 1RSB solution is more delicate than in the RS
case. We only numerically check the so-called type II stability [38, 35],
which is done by adding noise to the obtained messages of the population
dynamics. We add a noise $|\eta|,\eta\sim\mathcal{N}(0,10^{-8})$ and normalize
the messages. If the population dynamics starting from the noisy messages
converges back to the previous non-noisy messages, we say that the solution is
stable. If the noisy messages go to the RS fixed point, we are in the case
where the small noise reconstruction is not possible, and we are in the ‘fake
cluster’ case discussed above. However, if the noisy messages converge to
another non-trivial fixed point, there is an instability of type II which
indicates that further steps of replica symmetry breaking are required. This
last case did not arise, giving an indication that the 1RSB result is correct.

##  6 Freezing and Computational Hardness

It has been proven that a large class of algorithm fail to find solutions to
problems presenting the overlap gap property (OGP) [39, 40]. Informally, a
problem presents OGP if the clusters of solutions are well separated, in the
sense that there are no pairs of solutions at a certain intermediate distance.
Works in statistical physics put forward a conjecture that solutions within
frozen clusters (a weaker condition than OGP) are algorithmically difficult to
find (e.g. require exponential time) [20, 41, 21, 22]. It is also believed
that solutions to problems that are replica symmetric can be found in
polynomial time. We provide experimental evidence for these conjectures by
solving outer-totalistic CSPs using the belief propagation reinforcement
algorithm. The algorithm is presented below.

It is important to keep in mind that we only characterized the
thermodynamically relevant clusters, which contain typical solutions. However,
it is possible that non-dominating clusters have very different properties in
terms of the hardness of finding solutions that belong to them. It has been
shown that in the presence of atypical clusters with large local entropy,
solutions can still be found in polynomial time [42, 43]. One such case arises
for a rule of degree $d=4$ and is discussed below.

#### Belief propagation reinforcement

The belief propagation reinforcement algorithm, introduced in [15], is a
solver based on BP. Given a graph and a rule, the solver tries to produce a
solution. The main idea is to reinforce, i.e. to bias, the variables in the
direction of the marginals obtained by BP. The BP-reinforcement equation reads

$\psi^{i\rightarrow j}_{s_{i},s_{j}}=\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus
j}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus j})\prod_{k\in\partial
i\setminus j}b^{(k)}_{s_{k}}\psi_{s_{k},s_{i}}^{k\rightarrow i}.$ (23)

where $b^{(k)}_{s_{k}}$ is the bias. To update the bias, we use the same
heuristics as in [20]:

$\begin{split}b_{1}^{(i)}=\pi,\,b_{0}^{(i)}=1-\pi,\,\text{if}\,\chi_{0}^{(i)}>\chi_{1}^{(i)}\\\
b_{1}^{(i)}=1-\pi,\,b_{0}^{(i)}=\pi,\,\text{if}\,\chi_{0}^{(i)}\leq\chi_{1}^{(i)}\end{split}$
(24)

where $0\leq\pi\leq 1/2$ is a hyperparameter and $\chi_{s_{i}}^{(i)}$ is the
BP estimate of the marginal probability that node $i$ is in state $s_{i}$.
$\chi_{s_{i}}^{(i)}$ is given by

$\chi_{s_{i}}^{(i)}=\frac{1}{Z^{\text{site}}_{\chi}}\prod_{k\in\partial
i}\left(\chi_{0,s_{i}}^{k\rightarrow i}+\chi_{1,s_{i}}^{k\rightarrow
i}\right)$ (25)

where $Z^{\text{site}}_{\chi}$ is the normalization.
$\chi_{s_{i},s_{j}}^{i\rightarrow j}$ are the BP estimate of the marginal
probability that the variable node $ij$ in the factor graph takes values
$(s_{i},s_{j})$. This is given by

$\chi_{s_{i},s_{j}}^{i\rightarrow j}=\frac{\psi_{s_{i},s_{j}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i}}{Z_{\chi}}$ (26)

where $Z_{\chi}$ is the normalization. Each bias $b_{s_{i}}^{(i)}$ is updated
following (24), but the update is done only with probability

$p(t)=1-(1+t)^{-\gamma}$ (27)

where $t$ is the iteration step and $\gamma$ a hyperparameter. We fix
$\gamma=0.1$ in our case. To obtain a possible solution $\mathbf{s}$ from the
algorithm, we take

$s_{i}=\operatorname{argmax}_{s_{l}}b_{s_{l}}^{(i)}\,\,\forall i=1,...,N.$
(28)

The pseudocode is given in Algorithm 1. In practice, $\pi$ is found by doing a
grid search starting with $\pi=0.5$ and decreasing it until a solution is
found or $\pi=0$ is reached. In that last case, the algorithm failed to find a
solution. We use $25$ values of $\pi$ equidistant between $0$ and $0.5$, and a
maximum number of iterations $T=10^{4}$. Each time the algorithm is restarted,
the bias and the messages are initialized at random.

Algorithm 1 BP-REINFORCEMENT $(G,T,\gamma,\pi,\epsilon)$

1: Initialize the bias $b_{s_{i}}^{i}$ and the messages
$\psi_{s_{i},s_{j}}^{i\rightarrow j}$ randomly on the factor graph of $G$ and
normalize them.

2: $t\leftarrow 0$

3: Compute the current configuration
$s_{i}=\operatorname{argmax}_{s_{l}}b_{s_{l}}^{i}$

4: while $\mathbf{s}$ is not a solution and $t\leq T$ do

5: Update all the messages $\psi_{s_{i},s_{j}}^{i\rightarrow j}$ according to
(23) with dampening $\epsilon$.

6: Update every bias $b_{s_{i}}^{(i)}$ with probability $p(t)$ according to
(24)

7: Update $s_{i}=\operatorname{argmax}_{s_{l}}b_{s_{l}}^{(i)}$

8: $t\leftarrow t+1$

9: end while

#### d=3 rules

We focus on the problems that have $\phi\geq 0$ with no homogeneous solutions,
since the rules that have a homogeneous solution are trivial to solve. In the
case $d=3$, we tried to find solutions on $10$ different graphs of size
$N=10^{5}$ for each rule. Solutions were found for all $10$ graphs for all
rules except for the 1RSB rules with frozen variables 1+-0, where a solution
on only one graph was found (which we attribute to a finite size effect). For
the rule 1010 no solutions were found. This is shown in Table 4 (left). These
results are in agreement with the conjecture that frozen solutions are hard to
find for robust algorithms, and are also coherent with the hypothesis that the
‘fake clustered’ rules (marked RS*) are algorithmically easy. We also found no
solution for the rules with negative free entropy density but that are not
locally contradictory.

Note the special case of rule 1010, which can be mapped to a system of linear
equations modulo $2$ and can thus be solved in polynomial time. Using Gaussian
elimination, we found a solution for each of the $10$ graphs of size $8\cdot
10^{4}$ that we tested. Rules +---, -+--, +-+-, \--10, -10-, 0101 and 0-10 can
also be mapped to a system of linear equations over finite fields, but are of
lesser interest presenting either a trivial homogeneous solution or no
solutions at all. The 1010 rule is a parity check and is known to be solvable
in polynomial time. A similar situation arises in XORSAT problems, where the
space of solutions also exhibits clusters and frozen variables [44, 45] while
the problem is solvable in polynomial time [46, 47]. However, note that the
algebraic structure allowing Gaussian elimination disappears if one adds noise
to the problem [46]. If we focus on algorithms robust to noise the problem
0101 will also be algorithmically hard.

|  Satisfiable $d=3$ rules with no homogeneous solutions |
---|---|---
Rule | Phase | Rate of success
-0+0 | RS* | 1
-0+- | RS* | 1
-+00 | RS | 1
-+0- | RS | 1
-++0 | RS | 1
-++- | RS | 1
-+10 | RS* | 1
-1+0 | RS | 1
-110 | RS | 1
1000 | RS | 1
10+0 | RS* | 1
1+00 | RS | 1
1++0 | RS | 1
1+-0 | d1RSB | 0.1
1100 | RS | 1
1010 | d1RSB | 0
(a)

|  Rigid d1RSB $d=4$ rules with no homogeneous solution
---|---
Rule | Rate of success
-0+-0 | 0.5
-0+-- | 0
-+-+0 | 0
-+-+- | 0
-+10- | 0
\--+00 | 0
-1000 | 0
-10+0 | 0
-1+-0 | 0
-1-+0 | 0
10100 | 0
1+-00 | 1
(b)

Table 4: Left: Rate of success of the belief propagation reinforcement
algorithm on random 3-regular graphs with $N=10^{5}$ nodes for 10 different
graphs. The studied rules have $\phi\geq 0$ and no homogeneous solutions.
Right: Rate of success using the belief propagation reinforcement algorithm on
random 4-regular graphs with $N=10^{5}$ nodes for 10 different graphs for the
stable and frozen d1RSB rules with no homogeneous solutions.

#### d=4 rules and local entropy

We investigate the case $d=4$, for which we identify $12$ stable d1RSB rules
with no homogeneous solutions and containing frozen variables out of the
possible 528 non-equivalent rules. Using again BP reinforcement on $10$ graphs
of size $N=10^{5}$, we find no solution except for rule -0+-0 for which we
found solutions on 5 out of the 10 graphs, and rule 1+-00 for which a solution
was found on every tested graph, see Table 4 (right).

As discussed before, we expect the found solutions to belong to non-dominating
clusters with large local entropy. To investigate this, [48] or more recently
[49] study the local entropy density $s_{LE}(\mathbf{s},d)$, which is the
logarithm of the number of solutions at a given (normalized) Hamming distance
$d$ from a given solution $\mathbf{s^{\prime}}$ (note the difference with the
entropy of the Franz-Parisi potential [50], which is equivalent to the local
entropy in the case where the given solution $\mathbf{s}^{\prime}$ is
typical). They argue that solutions in wide and flat connected regions of
solutions can be algorithmically easy to find. These regions are characterized
by the absence of a gap in their local entropy. More precisely, we say that
there is a gap if there exist two distances $\alpha>0$ and $\beta>0$ with
$\alpha<\beta$ such that $\phi_{LE}(\mathbf{s},\alpha)\geq
0,\,\phi_{LE}(\mathbf{s},\beta)\geq 0$ and for some $\gamma\in[\alpha,\beta]$
we have $\phi_{LE}(\mathbf{s},\gamma)<0$. We compute the local entropy within
the RS approximation using belief propagation on a given graph by adding a
Lagrange multiplier $\kappa$ tuning the distance to a given solution. The
studied partition function then is

$Z_{LE}(\mathbf{s^{\prime}},d)=\sum_{\mathbf{s}\in\mathcal{C}}\prod_{i=1}^{N}e^{\kappa\mathds{1}(s_{i}\neq
s^{\prime}_{i})}f(s_{i},\\{s_{j}\\}_{j\in\partial
i})f(s^{\prime}_{i},\\{s^{\prime}_{j}\\}_{j\in\partial i})$ (29)

where $\mathds{1}$ is the indicator function that is $1$ if the argument is
true and is $0$ otherwise. The local free entropy density is
$\phi_{LE}(\mathbf{s^{\prime}},d)=\log(Z_{LE}(\mathbf{s^{\prime}},d))/N$. The
local entropy is obtained by a Legendre transform and is given by
$s_{LE}=\phi_{LE}-\kappa d$. The BP equation, the distance and the local
entropy density are obtained following a similar procedure as in Appendix B.
To obtain the local entropy density at a chosen distance, the Lagrange
multiplier $\kappa$ is tuned dynamically during the BP iterations. The
derivation of the BP equation and the iteration algorithm is presented in more
details in Appendix E.

Figure 6 (left) presents the local entropy density as a function of the
distance to a solution found by using BP reinforcement on a graph of size
$N=10^{5}$ for the rule 1+-00. We notice that there is no gap, indicating a
flat and wide minima. This also indicates that the solution is atypical, since
typical solutions are clustered and thus present a gap. However, recall that
the RS free entropy is an upper bound of the true free entropy, which is not
necessarily reached, so a gap might emerge by considering further steps of
replica symmetry breaking.

On the other hand, Figure 6 (right) shows the local entropy density for rule
1+-0. Since finding solutions for this rule on large graphs is hard, we use
$N=1000$ where solutions could be consistently found. We notice the presence
of a gap, as expected for a solution within a cluster. However, note that the
gap is small and the obtained curve can slightly change depending on the
generated graph, probably due to its small number of nodes. We also checked
the unique solution found for $N=10^{5}$, and it does not present a gap in its
local entropy. We suppose that this solution was found from an atypical
generation of the graph due to finite $N$. The solutions of rule -0+-0 also do
not present a gap in their local entropy. However, solutions are found more
frequently in this case which means that the explanation of atypical graphs is
unlikely. It is possible that, in this case, sub-dominant dense clusters (in
which solutions have a large local entropy) exist, but that the solutions are
not always found by our choice of hyperparameters in the BP reinforcement
algorithm.

Figure 6: Left: Local entropy $S_{LE}$ as function of the normalized distance
$d$ for an empirical solution for rule 1+-00, $N=10^{5}$. The grey line in the
inset is the entropy with no constraints. While typical solutions are frozen,
the found solution is in a subdominant cluster with large local entropy (for
small $d$ the local entropy is close to the entropy with no constraints) and
contains no frozen variables. Right: Local entropy as a function of the
normalized distance for an empirical solution containing frozen variables for
rule 1+-0, $N=1000$. Notice the presence of a gap in the zoomed-in region
close to $0$.

## Conclusion

In this work, we provided the cavity equation under both the replica symmetric
and 1-step replica symmetry-breaking assumptions for all outer-totalistic
constraint satisfaction problems on $d$-regular graphs. By numerically solving
these equations, we determined the existence and number of solutions. The 1RSB
calculation also allowed to characterize the structure of the space of typical
solutions. Additionally, the computational hardness of the problems was
explored experimentally using the belief propagation reinforcement algorithm.
This analysis uncovered the following interesting phenomena:

  * •

For the $d=3$ rules, we find that the annealed free entropy is equal to the
quenched free entropy with the exception of rules +-00 and +-0+ which present
a static 1-step replica symmetry breaking. Thus, for these two rules, the
quenched free entropy is strictly lower than the annealed one.

  * •

For the $d=3$ rules, we uncovered one rule, 1+-0, where solutions are hard to
find algorithmically. This rule presents clustering and frozen variables.

  * •

In the $d=4$ case, multiple rules were found to present clustering and frozen
variables and no trivial homogeneous solution. While no solutions were found
for most of them, solutions to rule 1+-00 were easy to find. An analysis of
the local entropy near the found solutions indicates that these solutions
belong to a subdominant dense cluster, which does not contain frozen
variables.

Several extensions of this work could be interesting to consider. While we
limited our study to a binary alphabet, expanding our approach to more states
is straightforward using our formalism, which would then also include problems
such as graph-coloring. Additionally, investigating rules that go beyond the
nearest neighbors could also yield interesting insights. Concerning the
typical algorithmic complexity, the overlap gap property could be studied
using the cavity method. Indeed, the overlaps can be estimated using this
method, which could potentially provide additional insights into the
algorithmic hardness of certain problems. Finally, the outer-totalistic CSPs
were introduced in this work through the lens of stationary solutions of
cellular automata on graphs. Going beyond stationarity to study the dynamics
of graph cellular automata is feasible with the cavity method as shown in
recent works [24, 6]. Extending this line of work could prove an exciting
avenue into the dynamic properties of complex systems.

## Acknowledgements

We thank Vittorio Erba and Barbora Hudcová for insightful and helpful
discussions.

## References

  * [1] Martin Gardner.  Mathematical games - the fantastic combinations of john conway’s new solitaire game ‘life’.  Scientific American 223: 120-123, 1970.
  * [2] Paul Rendell.  Turing Machine Universality of the Game of Life.  Springer Cham, 2015.
  * [3] John Von Neumann and Arthur W. (Arthur Walter) Burks.  Theory of self-reproducing automata.  Urbana, University of Illinois Press, 1966.
  * [4] Stephen Wolfram.  Universality and complexity in cellular automata.  Physica D: Nonlinear Phenomena, 10(1):1–35, 1984.
  * [5] Matthew Cook.  Universality in elementary cellular automata.  Complex Systems, 15, 2004.
  * [6] Freya Behrens, Barbora Hudcová, and Lenka Zdeborová.  Dynamical phase transitions in graph cellular automata.  Phys. Rev. E, 109:044312, 2024.
  * [7] H. A. Bethe.  Statistical physics of superlattices.  Proc. Roy. Soc. London A, 150 :552–575, 1935.
  * [8] Judea Pearl.  Reverend Bayes on Inference Engines: A Distributed Hierarchical Approach.  Association for Computing Machinery, 1982.
  * [9] Jonathan S. Yedidia, William T. Freeman, and Yair Weiss.  Understanding belief propagation and its generalizations.  In Exploring artificial intelligence in the new millennium, pages 239–269. Morgan Kaufmann Publishers Inc., 2003.
  * [10] Marc Mézard and Andrea Montanari.  Information, Physics and Computation.  Oxford University Press; Illustrated edition, 2009.
  * [11] M. Mézard, G. Parisi, and M. A. Virasoro.  SK model: The replica solution without replicas.  Europhysics Letters, 1(2):77, 1986.
  * [12] M. Mézard and G. Parisi.  The bethe lattice spin glass revisited.  The European Physical Journal B - Condensed Matter and Complex Systems, 20(2):217–233, 2001.
  * [13] Marc Mézard and Giorgio Parisi.  The cavity method at zero temperature.  Journal of Statistical Physics, 111(1):1–34, 2003.
  * [14] Marc Mézard and Riccardo Zecchina.  Random k-satisfiability problem: From an analytic solution to an efficient algorithm.  Physical Review E, 66(5):056126, 2002.  Publisher: American Physical Society.
  * [15] Joël Chavas, Cyril Furtlehner, Marc Mézard, and Riccardo Zecchina.  Survey-propagation decimation through distributed local computations.  J Stat Mech, 2005, 2005.
  * [16] Giulio Biroli, Remi Monasson, and Martin Weigt.  A variational description of the ground state structure in random satisfiability problems.  The European Physical Journal B, 14(3):551–568, 2000.
  * [17] Marc Mézard, Giorgio Parisi, and Riccardo Zecchina.  Analytic and algorithmic solution of random satisfiability problems.  Science (New York, N.Y.), 297:812–5, 2002.
  * [18] Marc Mézard, Matteo Palassini, and Olivier Rivoire.  Landscape of solutions in constraint satisfaction problems.  Physical Review Letters, 95(20), 2005.
  * [19] Lenka Zdeborová and Florent Krzakala.  Phase transitions in the coloring of random graphs.  Physical Review E, 76(3):031131, 2007.
  * [20] Lenka Zdeborová and Marc Mézard.  Constraint satisfaction problems with isolated solutions are hard.  Journal of Statistical Mechanics: Theory and Experiment, 2008(12):P12004, 2008.
  * [21] Lenka Zdeborová and Florent Krzakala.  Statistical physics of inference: Thresholds and algorithms.  Advances in Physics, 65(5):453–552, 2016.
  * [22] David Gamarnik, Cristopher Moore, and Lenka Zdeborová.  Disordered systems insights on computational hardness.  Journal of Statistical Mechanics: Theory and Experiment, 2022(11), 2023.
  * [23] Carsten Marr and Marc-Thorsten Hütt.  Outer-totalistic cellular automata on graphs.  Physics Letters A, 373(5):546–549, 2009.
  * [24] Freya Behrens, Barbora Hudcová, and Lenka Zdeborová.  Backtracking dynamical cavity method.  Phys. Rev. X, 13:031021, 2023.
  * [25] Martin Weigt and Alexander K Hartmann.  Number of guards needed by a museum: A phase transition in vertex covering of random graphs.  Physical review letters, 84(26):6118, 2000.
  * [26] Jean Barbier, Florent Krzakala, Lenka Zdeborová, and Pan Zhang.  The hard-core model on random graphs revisited.  Journal of Physics: Conference Series, 473:012021, 2013.
  * [27] Giulio Biroli and Marc Mézard.  Lattice glass models.  Physical review letters, 88(2):025501, 2001.
  * [28] Olivier Rivoire, Giulio Biroli, Olivier Martin, and Marc Mézard.  Glass models on bethe lattices.  Eur Phys J B, 37:55–78, 2003.
  * [29] Luca Dall’Asta, Paolo Pin, and Abolfazl Ramezanpour.  Statistical mechanics of maximal independent sets.  Physical Review E, 80(6):061136, 2009.
  * [30] Jin-Hua Zhao, Yusupjan Habibulla, and Hai-Jun Zhou.  Statistical mechanics of the minimum dominating set problem.  Journal of Statistical Physics, 159:1154–1174, 2015.
  * [31] Freya Behrens, Gabriel Arpino, Yaroslav Kivva, and Lenka Zdeborová.  (dis)assortative partitions on random regular graphs.  Journal of Physics A: Mathematical and Theoretical, 55, 2022.
  * [32] Florent Krzakala and Lenka Zdeborová.  Statistical physics methods in optimization and machine learning, 2021.  Lecture notes.
  * [33] Thierry Mora.  Géométrie et inférence dans l’optimisation et en théorie de l’information.  phdthesis, Université Paris Sud - Paris XI, 2007.
  * [34] J. R. L. de Almeida and D. J. Thouless.  Stability of the sherrington-kirkpatrick solution of a spin glass model.  Journal of Physics A: Mathematical and General, 11(5):983, 1978.
  * [35] Lenka Zdeborová.  Statistical Physics of Hard Optimization Problems.  Theses, Université Paris Sud - Paris XI ; Univerzita Karlova v Praze, June 2008.
  * [36] Marylou Gabrié, Varsha Dani, Guilhem Semerjian, and Lenka Zdeborová.  Phase transitions in the q-coloring of random hypergraphs.  Journal of Physics A: Mathematical and Theoretical, 50(50), 2017.  Publisher: IOP Publishing.
  * [37] Marc Mézard and Andrea Montanari.  Reconstruction on trees and spin glass transition.  Journal of Statistical Physics, 124(6):1317–1350, 2006.
  * [38] Andrea Montanari, Giorgio Parisi, and Federico Ricci-Tersenghi.  Instability of one-step replica-symmetry-broken phase in satisfiability problems.  Journal of Physics A: Mathematical and General, 37(6):2073–2091, 2004.
  * [39] David Gamarnik and Madhu Sudan.  Limits of local algorithms over sparse random graphs.  The Annals of Probability, 45(4):2353–2376, 2017.  Publisher: Institute of Mathematical Statistics.
  * [40] David Gamarnik.  The overlap gap property: a geometric barrier to optimizing over random structures.  Proceedings of the National Academy of Sciences, 118(41):e2108492118, 2021.
  * [41] Haiping Huang and Yoshiyuki Kabashima.  Origin of the computational hardness for learning with binary synapses.  Physical Review E, 90(5):052813, 2014.  Publisher: American Physical Society.
  * [42] L. Dall’Asta, A. Ramezanpour, and R. Zecchina.  Entropy landscape and non-gibbs solutions in constraint satisfaction problems.  Physical Review E, 77(3), 2008.
  * [43] Emmanuel Abbe, Shuangping Li, and Allan Sly.  Binary perceptron: efficient algorithms can find solutions in a rare well-connected cluster.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2022, pages 860–873. Association for Computing Machinery, 2022.
  * [44] S. Cocco, O. Dubois, J. Mandler, and R. Monasson.  Rigorous decimation-based construction of ground pure states for spin-glass models on random lattices.  Physical Review Letters, 90(4):047205, 2003.  Publisher: American Physical Society.
  * [45] Thierry Mora and Marc Mézard.  Geometrical organization of solutions to random linear boolean equations.  Journal of Statistical Mechanics: Theory and Experiment, 2006(10), 2006.
  * [46] W. Barthel, A. K. Hartmann, M. Leone, F. Ricci-Tersenghi, M. Weigt, and R. Zecchina.  Hiding solutions in random satisfiability problems: A statistical mechanics approach.  Physical Review Letters, 88(18):188701, 2002.  Publisher: American Physical Society.
  * [47] Harri Haanpää, Matti Järvisalo, Petteri Kaski, and Ilkka Niemelä.  Hard satisfiable clause sets for benchmarking equivalence reasoning techniques.  Journal on Satisfiability, Boolean Modeling and Computation, 2(1):27–46, 2006.  Publisher: IOS Press.
  * [48] Carlo Baldassi, Alessandro Ingrosso, Carlo Lucibello, Luca Saglietti, and Riccardo Zecchina.  Subdominant dense clusters allow for simple learning and high computational performance in neural networks with discrete synapses.  Physical Review Letters, 115(12):128101, 2015.
  * [49] Carlo Baldassi, Enrico M. Malatesta, Gabriele Perugini, and Riccardo Zecchina.  Typical and atypical solutions in nonconvex neural networks with discrete and continuous weights.  Physical Review E, 108(2):024310, 2023.  Publisher: American Physical Society.
  * [50] Silvio Franz and Giorgio Parisi.  Recipes for metastable states in spin glasses.  Journal de Physique I, 5(11):1401–1415, 1995.
  * [51] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.

##  Appendix A Non-equivalent Outer-totalistic Rules

We want to remove the “empty $\leftrightarrow$ occupied” symmetry that
produces equivalent rules. Extending the procedure of [23], the operator
$\mathcal{T}:\xi\rightarrow 1-\xi,\,\xi\in\\{0,1\\}$, which switches 0s and
1s, is introduced. Abusing the notation, $\mathcal{T}$ can be applied to a
configuration $\mathbf{s}$ using the notation $\mathcal{T}\mathbf{s}$, e.g.
$\mathcal{T}(0,1,0,1,1)=(1,0,1,0,0)$. $\mathcal{T}$ can also be applied to a
rule $\alpha_{0}\alpha_{1}\ldots\alpha_{d}$ with
$\alpha_{0},\alpha_{1},\ldots,\alpha_{d}\in\\{\texttt{0,1,+,-}\\}$ using the
notation $\mathcal{T}(\alpha_{0}\alpha_{1}\ldots\alpha_{d})$, where the
elements $\\{\texttt{+,-}\\}$ are invariant under $\mathcal{T}$. For example,
$\mathcal{T}(\texttt{0+1-})=\texttt{1+0-}$. We consider here that the rule is
describing the dynamics for the next time step: + indicates that the node
stays in the same state, - that it switches state, 0 that is becomes $0$ and 1
that it becomes 1. One time-step update of a configuration $\mathbf{s}$ under
rule $\alpha_{0}\alpha_{1}\ldots\alpha_{d}$ is denoted
$\alpha_{0}\alpha_{1}\ldots\alpha_{d}\cdot\mathbf{s}$. It can be observed that
$\alpha_{0}\alpha_{1}\ldots\alpha_{d}\cdot\mathbf{s}=\mathcal{T}\left(\mathcal{T}(\alpha_{d}\ldots\alpha_{1}\alpha_{0})\cdot\mathcal{T}\mathbf{s}\right)$
where $\alpha_{d}\ldots\alpha_{1}\alpha_{0}$ is
$\alpha_{0}\alpha_{1}\ldots\alpha_{d}$ in reverse order. Thus, it can be said
that the rule $\mathcal{T}(\alpha_{d}\ldots\alpha_{1}\alpha_{0})$ is symmetric
to the rule $\alpha_{0}\alpha_{1}\ldots\alpha_{d}$. Notice that some rules are
self-symmetric, e.g. rule 0-1.

If each node has two neighbors, a total of $4^{3}=64$ outer-totalistic rules
exist. 8 of them are self-symmetric: 0+1, 1+0, 0-1, 1-0, +++, +-+, -+- and
\---. All the other rules have exactly one different symmetric rule. Thus, a
total of $8+(64-8)/2=36$ nonequivalent outer-totalistic rules are
obtained111Note that [23] only identifies 34 nonequivalent outer-totalistic
rules using the same method. The authors were contacted about this
discrepancy, and came to the conclusion that it was a typo.. This procedure
can easily be extended to identify the nonequivalent rules for $d>2$.

##  Appendix B Derivation of the BP Equation and Observables on a Tree

In this appendix, we give a detailed derivation of eq. (10) and (12) starting
from the graphical model presented in Section 2. We will treat a more general
case than the probability distribution defined in eq. (1):

$P_{G}(\mathbf{s})=\frac{\prod_{i=1}^{N}e^{\mu
s_{i}}f(s_{i},\\{s_{j}\\}_{j\in\partial i})}{Z_{G}}.$ (30)

where $\mu$ is the chemical potential, a Lagrange multiplier allowing us to
tune the density of occupied nodes $r=\sum_{i=1}^{N}s_{i}/N$. Note that
additional Lagrange multipliers can easily be added to compute other
quantities of interest, such as the local entropy as done in Section 6. To
obtain eq.(10) and (12), simply set $\mu=0$.

#### BP equation.

Suppose that the original graph is a tree. From our construction of the
graphical model, this implies that the factor graph is also a tree. We define
the auxiliary partition function $V_{s_{i},s_{j}}^{i\rightarrow j}$ as being
the partition function of the system above $ij$ with the values of
$s_{i},s_{j}$ fixed. This corresponds to all the variables in the branch of
the tree that include $i$ but with the edge $ij$ removed. Let us name this
sub-graph $G^{i\rightarrow j}$. $G^{i\rightarrow j}$ is shown in red in Figure
7. The auxiliary partition function reads

$V_{s_{i},s_{j}}^{i\rightarrow j}=\sum_{\\{s_{l}\\}_{l\in G^{i\rightarrow
j}}}\prod_{l\in G^{i\rightarrow j}}e^{\mu
s_{l}}f(s_{l},\\{s_{m}\\}_{m\in\partial l}).$ (31)

Figure 7: Factor graph with the subgraph $G^{i\rightarrow j}$ in red and the
subgraphs $\\{G^{k\rightarrow i}\\}_{k\in\partial i\setminus j}$ in blue.

Note that we indicate the neighborhood on the original graph for simplicity.
Supposing that the factor graph is a tree, fixing the value of variable
$s_{i},s_{j}$ makes the branches above $i$ independent of each other. Let us
name these branches $G^{k\rightarrow i}$ with $k\in\partial i\setminus j$.
They are shown in blue in Figure 7. Thus, we can split the sum according to
the branch they belong to into a recursive relation:

$\begin{split}V_{s_{i,j}}^{i\rightarrow j}&=\sum_{\\{s_{k}\\}_{k\in\partial
i\setminus j}}e^{\mu s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\\\ &\,\,\,\,\,\,\,\,\,\,\,\,\times\prod_{k\in\partial i\setminus
j}\sum_{\\{s_{l}\\}_{l\in G^{k\rightarrow i}}}\prod_{l\in G^{k\rightarrow
i}}e^{\mu s_{l}}f(s_{l},\\{s_{m}\\}_{m\in\partial l})\\\
&=\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}V_{s_{k,i}}^{k\rightarrow i}.\end{split}$
(32)

Using (32) starting from the leaves, one can obtain all the auxiliary
partition functions of the tree.

Define the message $\psi^{i\rightarrow j}_{s_{i,j}}$ as the normalized
auxiliary partition function

$\psi^{i\rightarrow j}_{s_{i},s_{j}}:=\frac{V_{s_{i},s_{j}}^{i\rightarrow
j}}{\sum_{s_{i},s_{j}}V_{s_{i},s_{j}}^{i\rightarrow j}},$ (33)

where $\sum_{s_{i},s_{j}}$ is the sum over all possible values of $s_{i}$ and
$s_{j}$, e.g. in our case $s_{i}=s_{j}=0$, $s_{i}=0$ and $s_{j}=1$, $s_{i}=1$
and $s_{j}=0$, $s_{i}=s_{j}=1$. Messages are more convenient to work with than
the auxiliary partition function, since they do not scale exponentially in $N$
and represent the probability that variable $ij$ takes value $s_{i},s_{j}$
when restricted to the sub-graph $G^{i\rightarrow j}$ (the red part of Figure
7). Using the recursion relation (32) and the definition (33), we can also
obtain a recursion relation for the messages:

$\begin{split}\psi^{i\rightarrow
j}_{s_{i},s_{j}}&=\frac{\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}V_{s_{k},s_{i}}^{k\rightarrow
i}}{\sum_{s_{i},s_{j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}V_{s_{k},s_{i}}^{k\rightarrow i}}\\\
&\,\,\,\,\,\,\,\,\times\overbrace{\frac{\prod_{k\in\partial i\setminus
j}\sum_{s_{k},s_{i}}V_{s_{k},s_{i}}^{k\rightarrow i}}{\prod_{k\in\partial
i\setminus j}\sum_{s_{k},s_{i}}V_{s_{k},s_{i}}^{k\rightarrow i}}}^{=1}\\\
&=\frac{\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow
i}}{\sum_{s_{i},s_{j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow
i}}\end{split}$ (34)

Noticing that

$Z^{i\rightarrow j}:=\sum_{s_{i},s_{j}}\sum_{\\{s_{k}\\}_{k\in\partial
i\setminus j}}e^{\mu s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{(ki)}}^{k\rightarrow i}$ (35)

is a normalization, we obtain the following self-consistent equation for the
message:

$\psi^{i\rightarrow j}_{s_{i},s_{j}}=\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow i}.$ (36)

This is the Belief Propagation equation given in eq. (6). As for the auxiliary
partition functions, one can obtain recursively each message by starting from
the leaves of the tree.

#### Observables.

Starting from the leaves, we can collect recursively the contribution of each
branch to the partition function of the original problem to obtain the full
partition function:

$Z_{G}=\sum_{s_{r}}\sum_{\\{s_{k}\\}_{k\in\partial r}}e^{\mu
s_{r}}f(s_{r},\\{s_{k}\\}_{k\in\partial r})\prod_{k\in\partial
r}V^{k\rightarrow r}_{s_{k},s_{r}}$ (37)

where $r$ indicates a root of the tree. We can write this expression as a
function of the messages instead of the auxiliary partition function:

$\begin{split}Z_{G}&=Z_{G}\overbrace{\frac{\prod_{i=1,i\neq
r}^{N}\sum_{s_{i}}\sum_{s_{j}}V_{s_{j},s_{i}}^{j\rightarrow
i}V_{s_{i},s_{j}}^{i\rightarrow j}}{\prod_{(ij)\in
E(G)}\sum_{s_{i}}\sum_{s_{j}}V_{s_{j},s_{i}}^{j\rightarrow
i}V_{s_{i},s_{j}}^{i\rightarrow j}}}^{=1}\\\ &=Z_{G}\frac{\prod_{i=1,i\neq
r}^{N}\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}V_{s_{k},s_{i}}^{k\rightarrow i}}{\prod_{ij\in
E(G)}\sum_{s_{i}}\sum_{s_{j}}V_{s_{j},s_{i}}^{j\rightarrow
i}V_{s_{i},s_{j}}^{i\rightarrow j}}\\\
&=\frac{\prod_{i=1}^{N}\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}V_{s_{k},s_{i}}^{k\rightarrow i}}{\prod_{ij\in
E(G)}\sum_{s_{i}}\sum_{s_{j}}V_{s_{j},s_{i}}^{j\rightarrow
i}V_{s_{i},s_{j}}^{i\rightarrow j}}\\\
&=\left(\prod_{i=1}^{N}\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}e^{\mu s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}V_{s_{k},s_{i}}^{k\rightarrow i}}{\prod_{k\in\partial
i}\sum_{s_{k},s_{i}}V_{s_{k},s_{i}}^{k\rightarrow i}}\right)\\\
&\,\,\,\,\,\,\,\,\,\,\times\left(\prod_{ij\in
E(G)}\frac{\sum_{s_{i}}\sum_{s_{j}}V_{s_{j},s_{i}}^{j\rightarrow
i}V_{s_{i},s_{j}}^{i\rightarrow
j}}{\sum_{s_{i},s_{j}}V_{s_{i},s_{j}}^{i\rightarrow
j}\sum_{s_{j},s_{i}}V_{s_{j},s_{i}}^{j\rightarrow i}}\right)^{-1}\\\
&=\frac{\prod_{i=1}^{N}\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i}}{\prod_{ij\in
E(G)}\sum_{s_{i}}\sum_{s_{j}}\psi_{s_{i},s_{j}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i}}\end{split}$ (38)

where $E(G)$ denotes the set of edges of $G$. The property of the first line
can be seen if one considers a tree with $i\rightarrow j$ going in the
direction of the root. To go from the 1st to the 2nd line, we used the
recursion equation (32). From the 2nd to 3rd line, we see from (37) that
$Z_{G}$ can be absorbed in the product over the nodes. We multiplied line 3 by
$(\prod_{k\in\partial i}\sum_{s_{k},s_{i}}V_{s_{k},s_{i}}^{k\rightarrow
i})(\sum_{s_{i},s_{j}}V_{s_{i},s_{j}}^{i\rightarrow
j}\sum_{s_{j},s_{i}}V_{s_{j},s_{i}}^{j\rightarrow i})^{-1}$ (which is equal to
$1$) to get to line 4. In the last equality we simply used the definition of
the message (33).

Defining

$Z^{i}=\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i},$ (39)
$Z^{ij}=\sum_{s_{i}}\sum_{s_{j}}\psi_{s_{i},s_{j}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i},$ (40)

we can more succinctly write

$Z_{G}=\frac{\prod_{i=1}^{N}Z^{i}}{\prod_{ij\in E(G)}Z^{(ij)}}.$ (41)

From definition (3), the expression (7) for the free entropy density as a
function of the messages is obtained:

$\phi_{G}=\frac{1}{N}\sum_{i=1}^{N}\log Z^{i}-\frac{1}{N}\sum_{ij\in E(G)}\log
Z^{ij}.$ (42)

The mean density of occupied nodes $\rho$ as a function of the message is
given by

$\rho=\frac{\partial\phi}{\partial\mu}=\frac{1}{N}\sum_{i=1}^{N}\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}s_{i}e^{\mu s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i}}{Z^{i}}.$ (43)

###  B.1 Belief Propagation Iterations on $d$-regular Graphs

We now consider the case of $d$-regular graphs. We will assume that the
messages are locally equivalent, i.e. that they are the same independently of
the nodes. This means that we reduce from $2|E(G)|$ messages to only one. The
fixed point equation (6) then reduces to (10), and we can iterate it:

$\tilde{\psi}_{s_{i},s_{j}}(t+1)=\frac{1}{\tilde{Z}}\sum_{\\{s_{k}\\}_{k=1}^{d-1}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k=1}^{d-1})\prod_{k=1}^{d-1}\psi_{s_{k},s_{i}}(t).$
(44)

where $t$ indicates the iteration and $\tilde{Z}$ is the normalization. In
practice, the message is updated with a dampening $\epsilon\in[0,1)$:

$\psi_{s_{i},s_{j}}(t+1)=\epsilon\psi_{s_{i},s_{j}}(t)+(1-\epsilon)\tilde{\psi}_{s_{i},s_{j}}(t+1)$
(45)

The message can be represented by a $2\times 2$ matrix $\psi$ indexed by
$s_{i}$ and $s_{j}$. We stop the iterations when the Frobenius norm of
$\psi(t+1)-\psi(t)$ is smaller than a given tolerance (or after a maximum
number of steps). We used a tolerance of $10^{-12}$ and a damping
$\epsilon=0.8$. The estimated free entropy (7) simplifies to

$\phi_{\rm RS}=\log\left(\sum_{s_{i},\\{s_{k}\\}_{k=1}^{d}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k=1}^{d})\prod_{k=1}^{d}\psi_{s_{k},s_{i}}\right)-\frac{d}{2}\log\left(\sum_{s_{i},s_{j}}\psi_{s_{i},s_{j}}\psi_{s_{j},s_{i}}\right)$
(46)

which we call the replica symmetric free entropy density. The $d/2$ factor
comes from the fact that a random d-regular graph of $N$ nodes has
$\frac{Nd}{2}$ edges. This is the expression given in eq. (12) where we write
$Z^{i}\rightarrow Z_{\text{n}}$ and $Z^{ij}\rightarrow Z_{\text{e}}$ to avoid
referencing the graph.

##  Appendix C One-step Replica Symmetry Breaking Cavity Equation

In this annex we recover the 1RSB eq. (18), the expression for the replicated
free entropy density (19) and an explicit formula of the internal free entropy
as function of the messages.

Suppose again that the graph is a tree. Let us write the self-consistent
equation (36) as

$\begin{split}\psi^{i\rightarrow j}_{s_{i},s_{j}}&=\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow i}\\\
&=:\mathcal{F}(\\{\mathbf{\psi}_{s_{k},s_{i}}^{k\rightarrow
i}\\}_{k\in\partial i\setminus j}).\end{split}$ (47)

If this is valid for every $s_{i},s_{j}\in S$, then we write

$\psi^{i\rightarrow j}=\mathcal{F}\left(\\{\psi^{k\rightarrow
i}\\}_{k\in\partial i\setminus j}\right)$ (48)

where a message $\psi^{a\rightarrow b}$ without the subscript indicates all
the components of the message (e.g. in the form of a matrix).

Supposing that the internal free entropy density $\phi_{\rm
int}(\tilde{\psi})=\phi_{\text{Bethe}}(\tilde{\psi})$ (i.e. that the Bethe-
Peierls approximation holds within each cluster) and using the expression (42)
for $\phi_{\text{Bethe}}$, we can write the partition function of (15)

$\begin{split}Z_{\nu}&=\sum_{\\{\mathbf{\tilde{\psi}}\\}}e^{Nm\phi_{\rm
int}(\mathbf{\tilde{\psi}})}\\\ &\stackrel{{\scriptstyle
N\rightarrow\infty}}{{=}}\int\prod_{ij\in E(G)}[d\psi^{i\rightarrow
j}]e^{m\left(\sum_{i=1}^{N}\log Z^{i}-\sum_{ij\in E(G)}\log Z^{ij}\right)}\\\
&\,\,\,\,\,\,\,\,\,\,\,\times\prod_{i=1}^{N}\prod_{j\in\partial
i}\delta\left(\psi^{i\rightarrow j}-\mathcal{F}(\\{\mathbf{\psi}^{k\rightarrow
i}\\}_{k\in\partial i\setminus j})\right)\\\ &=\int\prod_{ij\in
E(G)}[d\psi^{i\rightarrow j}]\prod_{i=1}^{N}\left(Z^{i}\right)^{m}\prod_{ij\in
E(G)}\left(Z^{ij}\right)^{-m}\prod_{i=1}^{N}\prod_{j\in\partial
i}\delta\left(\psi^{i\rightarrow j}-\mathcal{F}(\\{\mathbf{\psi}^{k\rightarrow
i}\\}_{k\in\partial i\setminus j})\right)\end{split}$ (49)

where the integral is over all possible messages with the condition that they
are normalized (i.e. that they are indeed summing to $1$) and $\delta(\cdot)$
is the Dirac delta. For notational purposes, we do not write the dependence of
$Z^{i}$, $Z^{(ij)}$ to $\psi^{i\rightarrow j}$. We see that eq. (49) is the
partition function of the probability distribution

$\begin{split}\tilde{P}\left(\\{\psi^{i\rightarrow j}\\}_{ij\in
E(G)}\right)=&\frac{1}{Z_{\mu}}\prod_{i=1}^{N}\left(Z^{i}\right)^{m}\prod_{ij\in
E(G)}\left(Z^{(ij)}\right)^{-m}\\\
&\,\,\,\,\times\prod_{i=1}^{N}\prod_{j\in\partial
i}\delta\left(\psi^{i\rightarrow j}-\mathcal{F}(\\{\mathbf{\psi}^{k\rightarrow
i}\\}_{k\in\partial i\setminus j})\right).\end{split}$ (50)

If the graph is a tree we can take the product over $ij$ such that the i’s are
never the same. Using the fact that $\delta(x)\delta(x)=\delta(x)$, equation
(50) becomes

$\tilde{P}\left(\\{\psi^{i\rightarrow j}\\}_{ij\in
E(G)}\right)=\frac{1}{Z_{\mu}}\prod_{ij\in
E(G)}\left(Z^{i}\right)^{m}\left(Z^{(ij)}\right)^{-m}\delta\left(\psi^{i\rightarrow
j}-\mathcal{F}(\\{\mathbf{\psi}^{k\rightarrow i}\\}_{k\in\partial i\setminus
j})\right).$ (51)

We can apply the cavity formalism on this probability distribution. Following
the same procedure as in B, the obtained self-consistent equation is

$\begin{split}P^{i\rightarrow j}(\psi^{i\rightarrow
j})=\frac{1}{\mathcal{Z}^{i\rightarrow j}}&\int\prod_{k\in\partial i\setminus
j}[d\psi^{k\rightarrow
i}]\left(\frac{Z^{i}}{Z^{ij}}\right)^{m}\delta\left(\psi^{i\rightarrow
j}-\mathcal{F}(\\{\mathbf{\psi}^{k\rightarrow i}\\}_{k\in\partial i\setminus
j})\right)\\\ &\times\prod_{k\in\partial i\setminus j}[P^{k\rightarrow
i}(\psi^{k\rightarrow i})]\end{split}$ (52)

where $\mathcal{Z}^{i\rightarrow j}$ is the normalization. We will from here
on write

$\prod_{k\in\partial i\setminus j}[d\psi^{k\rightarrow
i}]\times\prod_{k\in\partial i\setminus j}[P^{k\rightarrow
i}(\psi^{k\rightarrow i})]=\prod_{k\in\partial i\setminus j}dP^{k\rightarrow
i}.$ (53)

In the population dynamics algorithm (see Appendix D), we use the fact that
$Z^{i}/Z^{ij}=Z^{i\rightarrow j}$. Indeed,

$\begin{split}\frac{Z^{i}}{Z^{ij}}&=\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}e^{\mu s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow
i}}{\sum_{s_{i}}\sum_{s_{j}}\psi_{s_{i},s_{j}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i}}\\\
&=\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow
i}}{\sum_{s_{i}}\sum_{s_{j}}\left(\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus j}}e^{\mu
s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow
i}\right)\psi_{s_{j},s_{i}}^{j\rightarrow i}}\\\
&=\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial i}}e^{\mu
s_{i}}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i}}{\frac{1}{Z^{i\rightarrow
j}}\sum_{s_{i}}\sum_{s_{j}}\left(\sum_{\\{s_{k}\\}_{k\in\partial i\setminus
j}}e^{\mu s_{i}}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\psi_{s_{j},s_{i}}^{j\rightarrow i}\prod_{k\in\partial i\setminus
j}\psi_{s_{k},s_{i}}^{k\rightarrow i}\right)}\\\ &=Z^{i\rightarrow
j}\end{split}$ (54)

where we used the definitions (39), (40) of $Z^{i}$, $Z^{ij}$ in the first
line and the BP eq. (36) to go from the first to the second line.

The replicated free entropy density is obtained similarly as (42):

$\Psi(m)=\frac{1}{N}\sum_{i=1}^{N}\log\mathcal{Z}^{i}-\frac{1}{N}\sum_{ij\in
E(G)}\log\mathcal{Z}^{ij}$ (55)

with

$\mathcal{Z}^{i}:=\int\prod_{k\in\partial i}dP^{k\rightarrow
i}\left(Z^{i}\right)^{m},$ (56) $\mathcal{Z}^{ij}:=\int dP^{i\rightarrow
j}dP^{j\rightarrow i}\left(Z^{ij}\right)^{m}.$ (57)

To obtain the free entropy density of a cluster, we use

$\begin{split}\phi_{\rm int}=\frac{\partial\Psi(m)}{\partial
m}=&\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\mathcal{Z}^{i}}\int\prod_{k\in\partial
i}dP^{k\rightarrow i}\log(Z^{i})\left(Z^{i}\right)^{m}\\\
&-\frac{1}{N}\sum_{ij\in E(G)}\frac{1}{\mathcal{Z}^{ij}}\int dP^{i\rightarrow
j}dP^{j\rightarrow i}\log(Z^{ij})\left(Z^{ij}\right)^{m}.\end{split}$ (58)

The density is given by

$\begin{split}\rho&=\frac{\partial\phi}{\partial{\mu}}=\frac{1}{N}\sum_{i=1}^{N}\int\prod_{k\in\partial
i}dP^{k\rightarrow
i}\frac{\partial}{\partial\mu}\left(\frac{\log(Z^{i})(Z^{i})^{m}}{\mathcal{Z}^{i}}\right)\\\
&=\frac{1}{N}\sum_{i=1}^{N}\int\prod_{k\in\partial i}dP^{k\rightarrow
i}\left(\frac{\partial_{\mu}Z^{i}}{Z^{i}}\frac{(Z^{i})^{m}}{\mathcal{Z}^{i}}+\frac{\log(Z^{i})^{2}(Z^{i})^{m}\partial_{\mu}Z^{i}}{\mathcal{Z}^{i}}-\frac{\log(Z^{i})(Z^{i})^{m}}{(\mathcal{Z}^{i})^{2}}\partial_{\mu}\mathcal{Z}^{i}\right)\\\
&=\frac{1}{N}\sum_{i=1}^{N}\int\prod_{k\in\partial i}dP^{k\rightarrow
i}\frac{(Z^{i})^{m-1}\partial_{\mu}Z^{i}}{\mathcal{Z}^{i}}\end{split}$ (59)

where we notice that $Z^{ij}$, $\mathcal{Z}^{ij}$ do not depend on the
Lagrange multiplier $\mu$ and the two last terms of the second line cancel
each other out.

Moving to d-regular graphs, we will suppose again that the messages are
locally equivalent. Thus, there is no more dependence on the graph and we use
the RS free entropy density (12) instead of the Bethe free entropy density
(7). Eq. (52) becomes eq. (18) using the notation $dP^{a\rightarrow
b}\rightarrow dP(\psi)$. Similarly, the replicated free entropy density (55)
becomes (19) following the same procedure as for eq. (46) with
$\mathcal{Z}^{i}\rightarrow\mathcal{Z}_{n}$ and
$\mathcal{Z}^{ij}\rightarrow\mathcal{Z}_{e}$. The internal free entropy
density (58) becomes

$\phi_{\rm
int}=\frac{1}{\mathcal{Z}_{n}}\int\prod_{k=1,...,d}dP(\psi^{(k)})\log(Z_{\text{n}})\left(Z_{\text{n}}\right)^{m}-\frac{d}{2}\frac{1}{\mathcal{Z}_{e}}\int
dP(\psi^{(1)})dP(\psi^{(2)})\log(Z_{\text{e}})\left(Z_{\text{e}})\right)^{m}.$
(60)

The complexity can be obtained from (17) and is given by

$\Sigma(\phi_{\rm int})=\Psi(m)-m\phi_{\rm int}.$ (61)

The density (59) becomes

$\rho=\int\prod_{k=1,...,d}dP(\psi^{(k)})\frac{(Z_{\text{n}})^{m-1}\partial_{\mu}Z_{\text{n}}}{\mathcal{Z}_{n}}.$
(62)

##  Appendix D Population Dynamics

The pseudocode for population dynamics algorithms can be found for instance in
[35] or [10]. Below we describe it informally:

  1. 1.

Initialize a list of $M$ messages $\psi$. This list is called the population

  2. 2.

Choose uniformly at random $L=\lceil(1-\epsilon)M\rceil$ messages from the
population. We will call them the subpopulation. $\epsilon$ is the damping
parameter and plays a similar role as in the BP algorithm.

  3. 3.

Update each message of the subpopulation by doing one BP iterations. The BP
iteration is done by sampling $d-1$ neighbors uniformly at random from the
population for each message of the subpopulation.

  4. 4.

Compute the normalization $Z_{l}$ for each message of the updated
subpopulation.

  5. 5.

Sample $L$ messages from the updated subpopulation from the distribution
$P(l)=\frac{Z_{l}^{m}}{\sum_{i=1}^{L}Z_{i}^{m}}$.

  6. 6.

Replace the subpopulation from the sampled messages from step 5.

  7. 7.

Repeat from step 2.

The algorithm is stopped when the obtained distribution stabilizes (this is
checked by comparing distances between histograms of the population), or after
a maximal number of iterations. In the case were the algorithm is used to
determine the stability of the RS solution, the messages of the population are
initialize on noisy copies of the RS solutions. To approximate the solution of
the 1RSB equation, the messages of the population are initialized on hard-
fields, i.e. one component of the message is $1$ and the others are $0$, as
prescribed by the reconstruction on trees presented in Section 5.

The observables are then sampled from the population. This is done by sampling
uniformly at random $K$ messages from the population and computing the
observables by approximating the integrals (55) and (58). This is repeated
many times, and the average is then taken. Since the population can still
fluctuate after many iterations, the computation of the observables can be
done during the iterations of the population dynamics, after a stabilization
of the population. However, to avoid correlations, sampling should not be done
at each iteration. The average over the observables at different iterations is
then taken, and the standard deviation is used to indicate the error.

This task was vectorized and implemented using pytorch [51] to allow for large
populations (we use $M=10^{7}$) to be run in acceptable time. We use
$\epsilon=0.8$ and we consider that the populations stabilize after $8000$
iterations or when the distribution of the messages in the population ceases
to change significantly. The sampling of the observables is then done with
$K=2\cdot 10^{8}$ every $50$ iteration for a total of $2000$ iterations.

For the rules that are locked (when $\phi_{\rm int}=0)$, the numerical result
is not exactly zero due to the fact that a variable can be frozen but the
message is not a hard-field (e.g. for $\psi_{0,0}=\psi_{0,1}=0.5$). This can
lead to numerical errors due to the estimation of the integrals (20), (21) and
the fluctuations of the population. To circumvent this problem, we check if
each message of the population implies the value of a node, and put the
internal free entropy to $0$ in Tables 2 and 3 if this is the case.

##  Appendix E Local Entropy

Following the procedure of Appendix B with the partition function given by
(29), we obtain the following BP equations on a graph of size $N$:

$\psi^{i\rightarrow j}_{s_{i},s_{j}}=\frac{1}{Z^{i\rightarrow
j}}\sum_{\\{s_{k}\\}_{k\in\partial i\setminus
j}}e^{\kappa\mathds{1}(s^{\prime}_{i}\neq
s_{i})}f(s_{i},s_{j},\\{s_{k}\\}_{k\in\partial i\setminus
j})\prod_{k\in\partial i\setminus j}\psi_{s_{k},s_{i}}^{k\rightarrow i}.$ (63)

where we did not write the dependence to the given configuration
$\mathbf{s}^{\prime}$, and suppose that this configuration is a solution so as
to avoid an additional indicator function. $Z^{i\rightarrow j}$ is the
normalization such that $\sum_{s_{i}}\sum_{s_{j}}\psi^{i\rightarrow
j}_{s_{i},s_{j}}=1$.

The free entropy density is given by

$\phi_{LE}=\frac{1}{N}\sum_{i=1}^{N}\log Z^{i}-\frac{1}{N}\sum_{(ij)\in
E(G)}\log Z^{ij}.$ (64)

where

$Z^{i}=\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}e^{\kappa\mathds{1}(s^{\prime}_{i}\neq
s_{i})}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{k},s_{i}}^{k\rightarrow i}$ (65)
$Z^{ij}=\sum_{s_{i}}\sum_{s_{j}}\psi_{s_{i},s_{j}}^{i\rightarrow
j}\psi_{s_{j},s_{i}}^{j\rightarrow i}.$ (66)

The (normalized) distance $d(\mathbf{s},\mathbf{s}^{\prime})$ between two
configurations $\mathbf{s}$ and $\mathbf{s}^{\prime}$ is defined as

$d(\mathbf{s},\mathbf{s}^{\prime})=\frac{\sum_{i=1}^{N}\mathds{1}(s_{i}\neq
s^{\prime}_{i})}{N},$ (67)

where the indicator function $\mathds{1}$ is $1$ if the argument is true, and
$0$ otherwise. In the large $N$ limit, the average distance is very peaked and
is equal to the distance. The distance to $\mathbf{s}^{\prime}$ as function of
the messages is given by

$d=\frac{\partial\phi_{LE}}{\partial\kappa}=\frac{1}{N}\sum_{i=1}^{N}\frac{\sum_{s_{i}}\sum_{\\{s_{k}\\}_{k\in\partial
i}}e^{\kappa\mathds{1}(s^{\prime}_{i}\neq
s_{i})}f(s_{i},\\{s_{k}\\}_{k\in\partial i})\prod_{k\in\partial
i}\psi_{s_{(ki)}}^{k\rightarrow i}}{Z^{i}}.$ (68)

Let $\mathcal{N}(\mathbf{s}^{\prime},d)$ be the exponential leading term of
the number of solutions at distance $d$ of $\mathbf{s}^{\prime}$. The local
entropy density $s_{LE}(\mathbf{s},d)$ is defined as
$\mathcal{N}(\mathbf{s}^{\prime},d)=e^{Ns_{LE}(\mathbf{s}^{\prime},d)}$.
Writing

$e^{N\phi_{LE}}=\sum_{d=i/N,i=0,\ldots,N}e^{\kappa
Nd(\mathbf{s}^{\prime},d)}\mathcal{N}(\mathbf{s}^{\prime},d)$ (69)

and applying the Laplace method we obtain

$s_{LE}=\phi_{LE}-\kappa d.$ (70)

This standard reasoning of statistical physics is also used to obtain eq.
(17).

The Legendre transform only gives the concave envelope of the local entropy.
To obtain non-concave behavior, we can fix the distance $d^{*}$ and adapt
$\kappa$ to obtain this distance. In practice, at every $t_{\kappa}$
iterations of the BP eq. (63) we update
$\kappa\leftarrow\kappa+\gamma(d^{*}-d)$, where $d$ is the distance computed
from (68) at the current BP messages. $\gamma$ is a factor determining how
much $\kappa$ is changed. We use $\gamma=5$ and $t_{\kappa}=1$ to obtain
Figures 6. The pseudocode is given in Algorithm 2. After convergence or after
$T$ steps, the free local entropy, the distance and the local entropy can be
computed from the messages using eqs. (64), (68) and (70).

Algorithm 2 BP-LOCAL ENTROPY
$(G,\mathbf{s}^{\prime},d^{*},T,\gamma,t_{\kappa},\epsilon,\epsilon_{BP},\epsilon_{d})$

1: Initialize the messages $\psi_{s_{i},s_{j}}^{i\rightarrow j}$ randomly on
the factor graph of $G$ and normalize them, we denote $\psi^{i\rightarrow
j}=\begin{pmatrix}\psi^{i\rightarrow j}_{0,0}&\psi^{i\rightarrow j}_{0,1}\\\
\psi^{i\rightarrow j}_{1,0}&\psi^{i\rightarrow j}_{1,1}\end{pmatrix}$ a
$2\times 2$ matrix and $\psi$ the list of all the messages.

2: $\psi^{i\rightarrow j}_{old}\leftarrow\psi^{i\rightarrow j}$

3: $\kappa\leftarrow 0$

4: $t\leftarrow 0$

5: while $\frac{1}{2|E|}\sum_{(i\rightarrow j)}||\psi^{i\rightarrow
j}-\psi^{i\rightarrow j}_{old}||_{F}>\epsilon_{BP}$ and
$|d(\psi)-d^{*}|>\epsilon_{d}$ and $t\leq T$ do

6: $\psi^{i\rightarrow j}_{old}\leftarrow\psi^{i\rightarrow j}$

7: if $t_{\kappa}$ divides $t$ then

8: $\kappa\leftarrow\kappa+\gamma(d^{*}-d(\psi))$

9: end if

10: Update all the messages $\psi_{s_{i},s_{j}}^{i\rightarrow j}$ according to
(63) with dampening $\epsilon$.

11: $t\leftarrow t+1$

12: end while

13: return $\psi$, $\kappa$.

##  Appendix F Equivalence between the Occupation Problems and the +- Rules

The occupation problems studied in [20] can be seen as the rules containing
only + and - (the rules does not depend on the state of the middle node). It
was found that, even though the graphical model is not the same, the entropy
matches.

Fixing $L=K=d$, [20] gives the following BP fixed point equation:

$\begin{split}\psi_{0}&=\frac{1}{Z_{reg}}\sum_{r=0}^{d-1}A_{r}\binom{d-1}{r}\psi_{1}^{(d-1)r}\psi_{0}^{(d-1)(d-1-r)}\\\
\psi_{1}&=\frac{1}{Z_{reg}}\sum_{r=0}^{d-1}A_{r+1}\binom{d-1}{r}\psi_{1}^{(d-1)r}\psi_{0}^{(d-1)(d-1-r)}\end{split}$
(71)

where $A_{r}$ is $1$ if the symbol at position $r+1$ of the rule is a + and
$0$ otherwise.

Our eq. (10) simplifies in this case to

$\begin{split}\psi_{0,0}&=\frac{1}{Z_{\psi}}\sum_{r=0}^{d-1}A_{r}\binom{d-1}{r}\psi_{0,0}^{d-1-r}\psi_{1,0}^{r}\\\
\psi_{1,0}&=\frac{1}{Z_{\psi}}\sum_{r=0}^{d-1}A_{r}\binom{d-1}{r}\psi_{0,1}^{d-1-r}\psi_{1,1}^{r}\\\
\psi_{0,1}&=\frac{1}{Z_{\psi}}\sum_{r=0}^{d-1}A_{r+1}\binom{d-1}{r}\psi_{0,0}^{d-1-r}\psi_{1,0}^{r}\\\
\psi_{1,1}&=\frac{1}{Z_{\psi}}\sum_{r=0}^{d-1}A_{r+1}\binom{d-1}{r}\psi_{0,1}^{d-1-r}\psi_{1,1}^{r}.\end{split}$
(72)

$\psi_{s_{i}}$ indicates the probability that a node takes the value $s_{i}$
when we have a cavity. We verified numerically that
$\psi_{s_{i}}=\psi_{0,s_{i}}+\psi_{1,s_{i}}$ on all the occupation rules for
$d=3,4$.

Summing $\psi_{0,0}+\psi_{1,0}$ and $\psi_{0,1}+\psi_{1,1}$ we obtain

$\begin{split}\psi_{0,0}+\psi_{1,0}&=\sum_{r=0}^{d-1}A_{r}\binom{d-1}{r}\frac{1}{Z_{\psi}}\overbrace{(\psi_{0,0}^{d-1-r}\psi_{1,0}^{r}+\psi_{0,1}^{d-1-r}\psi_{1,1}^{r})}^{:=B(r)}\\\
\psi_{0,1}+\psi_{1,1}&=\sum_{r=0}^{d-1}A_{r+1}\binom{d-1}{r}\frac{1}{Z_{\psi}}(\psi_{0,0}^{d-1-r}\psi_{1,0}^{r}+\psi_{0,1}^{d-1-r}\psi_{1,1}^{r}).\end{split}$
(73)

Let $\psi_{0}=\psi_{0,0}+\psi_{1,0}$ and $\psi_{1}=\psi_{0,1}+\psi_{1,1}$. The
BP equations (71) become

$\begin{split}\psi_{0}=\sum_{r=0}^{d-1}A_{r}\binom{d-1}{r}\frac{1}{Z_{reg}}\overbrace{\left(\sum_{r^{\prime}=0}^{d-1}\binom{d-1}{r^{\prime}}\psi_{0,1}^{d-1-r^{\prime}}\psi_{1,1}^{r^{\prime}}\right)^{r}\left(\sum_{r^{\prime\prime}=0}^{d-1}\binom{d-1}{r^{\prime\prime}}\psi_{0,0}^{d-1-r^{\prime\prime}}\psi_{1,0}^{r^{\prime\prime}}\right)^{d-1-r}}^{:=C(r)}\\\
\psi_{1}=\sum_{r=0}^{d-1}A_{r+1}\binom{d-1}{r}\frac{1}{Z_{reg}}\left(\sum_{r^{\prime}=0}^{d-1}\binom{d-1}{r^{\prime}}\psi_{0,1}^{d-1-r^{\prime}}\psi_{1,1}^{r^{\prime}}\right)^{r}\left(\sum_{r^{\prime\prime}=0}^{d-1}\binom{d-1}{r^{\prime\prime}}\psi_{0,0}^{d-1-r^{\prime\prime}}\psi_{1,0}^{r^{\prime\prime}}\right)^{d-1-r}\end{split}$
(74)

where we used the binomial expansion of $(\psi_{0,0}+\psi_{1,0})^{d-1}$ and
$(\psi_{0,1}+\psi_{1,1})^{d-1}$.

We would like to show that the BP equations are equivalent, which is the case
if $B(r)Z_{reg}=C(r)Z_{\psi}$ for every $r$. We verified numerically that the
found fixed point satisfies this equation for each occupation problem for
$d=3,4$. Additionally, inspecting the numerical values of the fixed point we
notice that

$\psi_{0,1}\psi_{1,0}=\psi_{0,0}\psi_{1,1}$ (75)
$\psi_{1,0}^{d-1}=\psi_{1,1}^{d-2}\psi_{0,1}.$ (76)

However, we do not know how to justify these relations. From these equations,
we obtain the substitutions

$\psi_{0,1}=\frac{\psi_{1,0}^{d-1}}{\psi_{1,1}^{d-2}}$ (77)

and

$\psi_{0,0}=\frac{\psi_{0,1}\psi_{1,0}}{\psi_{1,1}}=\frac{\psi_{1,0}^{d}}{\psi_{1,1}^{d-1}}.$
(78)

Using these substitutions we obtain analytically the desired relation
$B(r)Z_{reg}=C(r)Z_{\psi}\,\forall r$ for $d=3,4,5,6,7$.