Statistical physics of inference: Thresholds and algorithms
Lenka Zdeborová1,∗ , and Florent Krzakala2,∗

arXiv:1511.02476v5 [cond-mat.stat-mech] 22 Jan 2018

2

1
Institut de Physique Théorique,
CNRS, CEA, Université Paris-Saclay,
F-91191, Gif-sur-Yvette, France
Laboratoire de Physique Statistique, CNRS, PSL Universités
Ecole Normale Supérieure. Sorbonne Universités
Université Pierre & Marie Curie, 75005, Paris, France.
∗
lenka.zdeborova@cea.fr and florent.krzakala@ens.fr

Many questions of fundamental interest in today’s science can be formulated as inference problems:
Some partial, or noisy, observations are performed over a set of variables and the goal is to recover,
or infer, the values of the variables based on the indirect information contained in the measurements.
For such problems, the central scientific questions are: Under what conditions is the information
contained in the measurements sufficient for a satisfactory inference to be possible? What are the
most efficient algorithms for this task? A growing body of work has shown that often we can
understand and locate these fundamental barriers by thinking of them as phase transitions in the
sense of statistical physics. Moreover, it turned out that we can use the gained physical insight
to develop new promising algorithms. The connection between inference and statistical physics
is currently witnessing an impressive renaissance and we review here the current state-of-the-art,
with a pedagogical focus on the Ising model which, formulated as an inference problem, we call
the planted spin glass. In terms of applications we review two classes of problems: (i) inference of
clusters on graphs and networks, with community detection as a special case and (ii) estimating a
signal from its noisy linear measurements, with compressed sensing as a case of sparse estimation.
Our goal is to provide a pedagogical review for researchers in physics and other fields interested in
this fascinating topic.

2
CONTENTS

A. Organization of the manuscript
B. What is statistical inference?
1. Terminology and simple examples
2. Teacher-student scenario
3. Bayes optimal versus mismatched
4. Estimators
5. The high-dimensional limit
C. Statistical physics and inference
1. Basic dictionary
2. The teacher-student scenario defines the planted ensemble
3. Inference via Monte Carlo and variational mean-field
4. Useful statistical physics concepts

5
6
6
8
10
10
12
12
13
13
14
15

II. Planted spin glass as a paradigm of statistical inference
A. Definition of the planted spin glass
B. Nishimori’s mapping via gauge transformation
C. Bayes-optimality on the Nishimori line
D. Parameter mismatch and learning
E. Relations between planted, quenched and annealed disorders
1. The quenched and annealed average
2. Properties of the planted ensemble
3. Quiet planting
F. Replica symmetry breaking
G. No replica symmetry breaking in Bayes-optimal inference.
H. Applications of planting for studies of glasses
1. Equilibration for free in glasses
2. Following states and analyzing slow annealing
3. Melting and glass transition

18
18
20
21
22
23
23
24
25
25
26
27
27
28
29

III. Phase diagram from the cavity method
A. Random graphs
B. Belief propagation and the cavity method
C. Properties of the randomly-quenched ensemble
D. Phase diagram of the planted spin glass
E. Systems with first order phase transitions
F. Comments on the phase transitions in various problems

30
30
31
33
35
36
38

IV. From physics insight to new algorithmic ideas
A. Non-backtracking spectral method
B. On computational hardness of the hard phase
C. Spatial coupling
D. The spatially coupled Curie-Weiss model
E. A brief history of time indices in the TAP equations

40
40
41
42
43
44

V. Clustering of networks and community detection
A. Biased overview of community detection
B. Asymptotically exact analysis of the SBM
C. Examples of phase transitions in the SBM
D. Spectral redemption for clustering sparse networks
E. More on the non-backtracking operator
F. The deceptiveness of the variational Bayes method
G. Real networks are not generated by the SBM

47
47
48
50
52
54
55
55

VI. Linear estimation and compressed sensing
A. Some examples of linear estimation
B. Compressed sensing

57
57
58

3
C. Approximate message passing
1. Relaxed Belief propagation
2. From r-BP to G-AMP
3. The potential
4. Note on the convergence and parameter learning
5. Examples of priors and outputs
D. Replica free energy, state evolution and phase transitions
1. State Evolution
2. Free energy
3. Phase transitions
E. Spatial coupling
F. Non-random matrices and non-separable priors
1. Structured priors
2. Orthogonal matrices
3. Reconstruction in discrete tomography
G. An application in physics: Optics in complex media
H. Beyond linear estimation: matrix factorization
1. Learning random sparse dictionaries
2. Low rank decomposition
VII. Perspectives

59
60
62
63
64
64
65
65
67
67
68
69
69
70
70
71
71
72
72
74

Acknowledgements

76

Fundings

76

References

76

4
I.

INTRODUCTION

Our goal in this review is to describe recent developments in a rapidly evolving field at the interface between
statistical inference and statistical physics. Statistical physics was developed to derive macroscopic properties of
material from microscopic physical laws. Inference aims to discover structure in data. On the first look, these goals
seem quite different. Yet, it was the very same man, Pierre Simon, Marquis de Laplace (1749–1827), who did one of
the first sophisticated derivations of the gas laws within the caloric theory, and who also created the field of statistical
inference [208]. This suggests there may be a link after all. In fact, the methods and tools of statistical physics are
designed to describe large assemblies of small elements - such as atoms or molecules. These atoms can quite readily
be replaced by other elementary constituents - bits, nodes, agents, neurons, data points. This simple fact, and the
richness and the power of the methods invented by physicists over the last century, is at the roots of the connection,
and the source of many fascinating work in the last decades.
There is a number of books discussing this interdisciplinary connection. Without claim of exhaustivity, see for
instance [86, 111, 181, 217, 238]. Our focus here is to present some more recent contributions that were obtained
using methodology that originated in the field of spin glasses [214]. In particular the replica and cavity methods,
with the related message passing algorithms, are at the roots of our toolbox. We focus on the analogy between phase
transitions in thermodynamics, and the threshold phenomena that arise in the ability to infer something useful from
data. The most interesting outcome of these considerations are deep consequences for the possibility of development
of new algorithms, and their relations to computational feasibility. This focus is summarized in the title by the words:
“threshold and algorithms”.
As we shall see, the natural setup in which the connection between inference and the statistical physics of disordered
systems is particularly transparent is the so-called teacher-student scenario. In a nutshell, the teacher-student scenario
can be described as follows: The teacher uses some ground truth information and a probabilistic model to generate
data that he then hands to the student
who is supposed to recover the ground truth as well as he can only from the knowledge of the data and the model.
This formulation of inference problems is particularly popular in neural networks [86, 281, 310], where the teacher
uses a rule to generate output, and the student is then trying to learn the rule to reproduce the same output. The
advantage of this setting is that it is easily amenable to the Bayesian approach and strategies that recover the ground
truth in an optimal way can be analyzed.
Let us make several comments: First we do not make this Bayesian choice because we think the Bayesian approach
is more correct than other lines of thought common in statistics. We do it because it is in this setting that the
connection between statistical inference and statistical physics becomes the most transparent. Our main point is to
explain and illustrate the insight that can be obtained from this connection in terms of phase transitions. In a sense
the strategy followed in this review is closer to the approach in information theory where we know the properties
of the noise that corrupted the signal than to traditional approach in statistics where we aim to assume as little as
possible.
Second, we focus on the teacher-student scenario and consequently we do not analyze problems where one is given
data or observations and has absolutely no idea of the underlying system that generated them. Nor do we discuss
model selection where one needs to judge from data which model is more suitable. There is a range of textbooks on
statistics where such problems are studied [13, 122, 198]. In theoretical physics it is common to study abstract models
that display interesting behavior and use gained knowledge to generalize about what is happening in the real world.
Third, while the subject is certainly algorithmic in nature, we shall concern ourselves primarily, as the title shows,
with phase transitions and threshold phenomena, where the connection with statistical physics is clear. This means
that while we certainly are interested in analyzing and discussing algorithms, and ask if it is possible for them to
reach the phase transition boundaries, we do not review, for any given problem, the vast amount of literature used to
actually solve it in practice. Our goal, instead, is (in a typical statistical physics attitude) to reach some “universal”
understanding of the problem.
The typical reader we have in mind has some knowledge of statistical physics, knows e.g. about the Ising model,
phase transition and free energy. We have in mind a reader who is curious to know how methods and concepts from
statistical physics such as phase transitions translate into the context of machine learning and inference. Perhaps the
most important message we want to pass in this review is that many calculations that physics uses to understand the
world can be translated into algorithms that machine learning uses to understand data. We also have in mind readers
from statistics and computer science who might want to learn more about the statistical physics approach. While we
do explain some basic physics concepts, such readers could further read a short introduction to statistical physics, as
for instance the appendix B (“some physics”) in the book of David McKay [198], or the Chap II in [217].

5
A.

Organization of the manuscript

The present review is organized as follows.
In section I B we introduce the basic concepts of statistical inference and the terminology. Most of this section is
standard knowledge in statistics and the main purpose of this section is to set the notations and present the concepts
to readers from physics who are not particularly familiar with probability and statistics. The only part of sec. I B that
is perhaps not standard, but certainly crucial for the present discussion, is the idealized teacher-student setting of
Sec. I B 2. Traditional statistics would typically not assume knowledge of the model that generated the data and large
part of statistics is precisely concerned with determining a good model. The main focus of this review is, however, on
the connections to statistical physics. In the teacher-student setting, the analysis of exact phase diagrams and phase
transitions is possible and this is our primary goal.
In Sec. I C we present some basic concepts of statistical physics that non-physics readers might not be familiar
with. The emphasis is given on formulation of the teacher-student setting as a planted statistical ensemble – this
planted ensemble is not standard in physics but we do believe it is a really fundamental concept, that has become
popular in the last decades, and that introducing it greatly elucidates the deep connection between properties of
high-dimensional inference problems, optimization problem, and disordered systems.
Having introduced the basic concepts we then choose a bottom-up approach and present in sec. II a prototypical
inference problem that we call the planted spin glass. We discuss how the physics of this inference problem can
be deduced using the existing knowledge from the physics of spin glasses. The main highlight of this section is that
Bayes-optimal inference can be translated into thermodynamics on the Nishimori line in the underlying planted model.
The physical properties of the Nishimori line are well understood for a wide range of models, this connection hence
provides a unifying picture for a large class of inference problems. Consequently, concepts discussed in section II
bring a lot of insight into many inference problems other than the planted spin glass. We conclude section II with a
discussion on how the concept of planting turns out to be innovative in study of glasses in sec. II H. All the results of
this section are presented in a way that should be understandable without going into any calculations. The purpose
of this is twofold: Readers who are not familiar with the replica and the cavity calculations can still understand
the reasoning and the main results. Readers who are familiar with the calculations can appreciate the reasoning,
connections and consequences without having to separate them from the technicalities.
Section III, however, is a more analytical chapter, that should be read with a pencil so that one can redo the
computations. It summarizes the basic analytical technique based on belief propagation and the cavity method that
we can use to analyse the properties and phase diagrams of high-dimensional inference problems. We do not go into
many details and the related methods are presented in greater depths in other existing literature and textbooks,
e.g. [217]. Readers familiar with these techniques can skip this section, except perhaps for its last part III F where we
present the zoology of the main phase transitions that are encountered in inference problems.
Armed with the statistical physics’ understanding of inference from sections II, III, we discuss related recent
algorithmic ideas in Sec. IV. The various algorithmic ideas are discussed again on the example of the planted spin
glass model or its simple variants. This section can be read separately from the previous ones by readers interested
mainly in the algorithmic consequences stemming from physics.
The planted spin glass problem is a simple example. To illustrate how to generalize from it to other, more practically
relevant, inference problems we present in section V the problem of clustering of networks with its recent development,
and in section VI the generalized linear estimation problem with the compressed sensing as its special case. The two
examples of applications in sections V, VI are chosen with respect to author’s own work. Each of these last two
sections can be seen as a concise introduction into these practically relevant problems and is therefore mainly for
readers that are interested in various applications of physics concepts or who want to specifically start learning about
clustering or compressed sensing.
Finally, a word of caution about the subjects that we have chosen to cover (and more importantly not to cover).
Historically, one of the most prominent example of statistical inference studied in physics are the error correcting codes.
Their studies in physics were pioneered by Sourlas [290] and the results are reviewed nicely e.g. in [217, 225, 238].
We do not discuss error correction in this review, except for a couple of interesting and/or historical connections.
Another very prominent problem that can be formulated as statistical inference and that has been extensively
studied by tools of statistical physics are neural networks. The book [238] covers very nicely the early results. A very
detailed presentation, with an emphasis on teacher-student protocols, can be found in the book [86]. The reviews
[281, 310] are also very good references. While we do not focus on this application that is rather well documented,
we will see in section VI that the problem of perceptron is a special case of the generalized linear estimation setting.
The two applications, network clustering and compressed sensing, that we chose to present in this review have
beed studied only more recently (in the last decade) and there are interesting results stemming from the connection
between inference and physics that are not (yet) well covered in other reviews.

6
B.

What is statistical inference?

Inference is the act or process of deriving logical conclusions from premises known or assumed to be true. Statistical
inference is the process of deducing properties of an underlying distribution by analysis of data. A way to think about
inference problems, that is particularly relevant in the context of today’s data deluge (often referred to as big data), are
procedures to extract useful information from large amounts of data. As the boundaries between fields are breaking,
nowadays inference problems appear literally everywhere, in machine learning, artificial intelligence, signal processing,
quantitative biology, medicine, neuroscience and many others.
A more formal way to set an inference problem is the following: One considers a set of variables x = {xi }i=1,...,N
on which one is able to perform some partial observations or measurements y = {yµ }µ=1,...,M . The goal is to deduce,
or infer, the values of the variables x based on the perhaps indirect and noisy information contained in the observed
data y. To give a couple of examples, the variables may correspond for instance to image pixels, an acoustic signal
or to an assignment of data points into clusters, while the corresponding observations would correspond for instance
to a blurred photograph, a recording of a song performed on a platform while the train was arriving, or positions of
points in space.
For inference problems the two main scientific questions of interest are:
• The question of sufficient information – Under what conditions is the information contained in the observations
sufficient for satisfactory recovery of the variables?
• The question of computational efficiency – Can the inference be done in an algorithmically efficient way? What
are the optimal algorithms for this task?
These two questions are of very different nature. The first one is a question that belongs to the realm of statistics
and information theory. The second one, instead, is in the realm of computer science and computational complexity.
Before the era of big data, statistics was often not concerned with computational tractability, because the number of
parameters to be estimated was usually small. Nowadays the tradeoff between statistical and computational efficiency
has changed, making it most reasonable to study the two above questions together.
We will see that these two questions have a very clear interpretation in terms of phase transitions in physics. More
specifically, the first one is connected to a genuine phase transition while the second one is connected to the type
of the transition (first or second order), and metastability associated with first order phase transitions is related to
algorithmic hardness.
1.

Terminology and simple examples

In this review, we will follow the approach of Bayesian inference. Developed by Bayes and Laplace, it concentrates
on conditional probabilities. The probability of having an event A conditioned on the event B is denoted P (A|B) and
defined as
P (A|B) =

P (A, B)
,
P (A)

(1)

P
where P (A, B) is the joint probability of both events A and B, and P (A) = B P (A, B) is the probability of event A
(the sum is taken over all the possible realizations of event B). When events A, B are taken from a continuous domain
the P (A) denotes probability distribution function, whereas P (A)dA denotes the probability of observing an outcome
in the interval (A, A + dA). For discrete variables the probability of one event, and the probability distribution are
interchangeable. Note that in this review we discuss both models with discrete variables (sec. II or V) and with
continuous variables (sec. VI), the concepts that we build hold for both these cases.
The central point of the approach is the so-called Bayes formula that follows directly from the definition of the
conditional probability (1)
P (x|y) =

P (y|x)
P (x) .
P (y)

(2)

There is a precise terminology for all the terms entering this equation, so we shall start by discussing the Bayesian
inference vernacular:
• P (x|y) is called the posterior probability distribution for x given all the data y.
• P (y|x) is called the likelihood.

7
• P (y), also denoted as a normalisation Z(y), is called the evidence.
• P (x) is called the prior probability distribution for x. It reflects our belief of what x should be before we
received any data.
Example of decaying particles. Our first example is borrowed from Ref. [198]. Unstable particles are emitted
from a source and decay at a distance y from the source. Unfortunately, we cannot build a perfect infinite length
detector, and decay events can only be observed if they occur in a window extending from distances ymin = 1 cm to
ymax = 20 cm from the source.
We run the experiment for some time, and we observe M decays at locations y = {y1 , ..., yM }. We also know that
we should expect an exponential probability distribution with characteristic length x for the event. Writing things in
the inference language, y are the observations and x is the unknown signal (of dimension N = 1). The question we
would like to answer is: Given the observations, how to estimate x?
In this particular example physics tells us exactly what is the probability distribution of a particle decaying at
distance y given the characteristic length x. In terms of Bayesian probability this means that we know exactly the
model that created the observations. This is of course not true in general, but when it is true we can set up the
so-called Bayes-optimal inference: The likelihood for one observation is
(
1
e−y/x if 1 < y < 20 ,
(3)
P (y|x) = xM(x)
0
otherwise ,
R 20
with a normalisation M(x) = 1 dz x1 e−z/x = e−1/x − e−20/x . The observations are independent, and therefore in
this case their joint probability distribution is simply a product of the single observation distributions. From eq. (3)
we can now extract P (x|y) using the Bayes’ formula and get
P (x|y) =

P (x)

e−
M
Z(y) [xM(x)]

PM

i=1

yi /x

.

(4)

Looking at this equation, we see that we have P (x) on the right hand side. If we do not have any additional information
on x, we assume than the prior distribution P (x) is a constant that just enters in the x-independent normalization,
denoting then Z̃(y) = Z(y)/P (x), we obtain
P (x|y) =

1
M

Z̃(y) [xM(x)]

e−

PM

i=1

yi /x

.

(5)

This is the final answer from Bayesian statistics. It contains all the information that we have on x in this approach.
For a dataset consisting of several points, for instance the six points y = 1.5, 2, 3, 4, 5, 12 cm, one can compute the
most probable value of x to be around x̂ = 3.7 cm. Note that this is considerably different from the mere average of
the points which is y = 4.6 cm. The reader can convince herself (for instance performing simulation and comparing
to the true correct value) that the presented approach gives a way better result that the naive average.
Example of linear regression. Most physicists are used to fit a theoretical model to experimental data. Consider
for instance that we measure the height of an object falling due to Newton’s gravity, we assume no initial speed and
no friction. We expect that the height h(t) of the object at time t will follow h(t, g, h0 ) = h0 − gt2 /2, where g is the
acceleration of gravity. We can do an experiment and measure heights of the object yi for many values of ti , with
i = 1, . . . , M . We expect, if we plot yi against t2i /2, to get a straight line that will give us the value of the acceleration
g.
Due to small but real measurement errors, the measured points will not lie exactly on a line. Used to such a
situation, a physicist would in this case perform a linear regression to minimize the so-called sum of squares of errors.
The “best” value of g, or to use the statistics language the estimator ĝ of g, would simply be the one that minimizes
PM
the cost function i=1 [h(ti , g, h0 ) − yi ]2 .
In order to get more familiarity with the concepts of inference we discuss how this familiar procedure can be cast
into the Bayesian language. In general we may not know much about the nature of the measurement noise. In such
cases it is often reasonable to expect that the noise is a Gaussian random variable of zero mean and variance ∆. Under
this assumption the likelihood P (yi |g, ∆) of a single data point is a Gaussian with mean h(ti , g, h0 ) and variance ∆].
We assume that the noise for different measurements is independent which leads to the total likelihood
P (y|g, h0 , ∆) =

M
Y

1
(2π∆)

M
2

i=1

e−

[h(ti ,g,h0 )−yi ]2
2∆

(6)

8
As in the previous example we assume a uniform prior distribution on the acceleration constant P (g), initial height
P (h0 ) and the variance P (∆), the posterior probability distribution is then given by
P (∆, g, h0 |y) =

M
Y

1
Z(y)(2π∆)

M
2

e−

[h(ti ,g,h0 )−yi ]2
2∆

(7)

i=1

PM
Looking for the most probable value of g indeed amounts to minimizing the sum of squares
i=1 [h(ti , g, h0 ) −
yi ]2 , which is the standard solution. Note, however, that the assumption of Gaussian noise independent for every
measurement might not be justified in general. Also, depending on the precise setting, there might be better prior
distributions for the parameters g, h0 , ∆. For more advanced discussion of classical linear regression see for instance
[309]. In section VI of this review we will discuss a high-dimensional version of linear regression, when the number of
parameters to be estimated is comparable to the number of observations M .
Let us make some important comments on the above considerations, which have been at the roots of the Bayesian
approach to inference from the very start. First, we notice that probabilities are used here to quantify degrees of
belief. To avoid possible confusions, it must thus be emphasized that the true x in the real world is not a random
variable, and the fact that Bayesians use a probability distribution P (x|y) does not mean that they think of the
world as stochastically changing its nature between the states described by the different hypotheses. The notation of
probabilities is used here to represent the beliefs about the mutually exclusive hypotheses (here, values of x), of which
only one is actually true. The fact that probabilities can denote degrees of belief is at the heart of Bayesian inference.
The reader interested in these considerations is referred to the very influential work of Bruno de Finetti [70].
Another subtle point is how to choose the prior distribution P (x). In fact, discussions around this point have been
rather controversial in the beginning of the history of statistical inference, and even now there are many schools of
thought. It is not the goal of this review to present these discussions and debates. In the spirit of a reductionism
that is usually done in physics, we introduce in the next section the so-called teacher-student scenario in which we
know what the model and the prior distribution really are, so that the use of the Bayes formula is straightforward
and fully justified. Majority of this review is then set within this teacher-student framework. This can be seen as a
simplified model of the generic situation. As we will see the picture of what is happening is interestingly rich already
in this framework. A first step towards the generic situation of unknown models that generated the data and unknown
prior probabilities is the case when the teacher hides some parameters of the model or/and of the prior distribution
from the student. Analyzing this case provides insight into the generic case, but still within the realm of a solvable
simplified teacher-student framework.
In sec. I C 2 we then argue that the teacher-student scenario for a generic inference problem can be seen as a
so-called planted ensemble of the corresponding statistical physics model. This is a known but not widely used way
of presenting inference in statistical physics, but the authors feel that it is an extremely useful and insightful way. It
allows to extend the intuition gathered over several decades in studies of disordered systems such as glasses and spin
glasses into many currently studied inference problems that otherwise appear highly non-trivial. One of the goals
of the present review is to bring this insight forward and establish the way of thinking about inference problems as
planted statistical ensembles.
2.

Teacher-student scenario

The teacher-student scenario is defined as follows:
• Teacher: In a first step the teacher generates a realization of variables/parameters x∗ from some probability
distribution Ptp (x), where ’tp’ stands for teacher’s prior. In a second step he uses these ground truth values
x∗ to generate the data y from some statistical model characterized by a likelihood of y given x∗ , denoted as
Ptm (y|x∗ ), where ’tm’ stands for teacher’s model. Finally, the teacher hands the data y to the student together
with some information about the distributions Ptp (x) and Ptm (y|x).
• Student: The goal of the student is to infer as precisely as possible (and tractably, when computational complexity is included into considerations) the original values x∗ of the variables from the information provided by
the teacher, i.e. from the data y and the available information about the distributions Ptp (x) and Ptm (y|x).
The variables with a ∗ (e.g. x∗ ) denote the ground truth that was used by the teacher, but that is not available to
the student. In the present manuscript we will consider cases where both the x and y are high-dimensional vectors,
and where the teacher’s prior Ptp (x) and model Ptm (y|x) are parameterized by some scalar (or low-dimensional vector)
values θ∗ and ∆∗ , respectively. This is a case relevant to many real-world problems where we assume existence of many

9
latent variables (collected into the vector x), and we assume that the model and prior distribution is parametrized by
a small number of unknown parameters.
We will in particular be interested in inference problems where the teacher’s prior distribution is separable, i.e. can
be written as a product over components
Ptp (x) =

N
Y

Ptp (xi ) .

(8)

i=1

Further, in this review, we are interested in cases where also the observations {yµ }µ=1,...,M are independent and each
of them depends on a subset of variables denoted ∂µ. The likelihood of the teacher’s model is then written as
Ptm (y|x) =

M
Y
µ=1

Ptm (yµ |{xi }i∈∂µ ) .

(9)

Given the information available to the student, she then follows the strategy of Bayesian statistics. The prior
information she has from the teacher is denoted as P (x), where x are the variables to be inferred. The statistical
model is represented by its likelihood P (y|x). Again we assume that both these can be written in a product form as
above. She then considers Bayes’ theorem (2) and writes the posterior probability which contains all the information
available to her
P (x|y) =

N
M
Y
1 Y
P (xi )
P (yµ |{xi }i∈∂µ ) .
Z(y) i=1
µ=1

(10)

A very useful quantity to define that derives directly from the posterior distribution is the marginal probability of
one variable
Z
Y
µi (xi ) = P (x|y)
dxj .
(11)
j6=i

Examples:
• Let us set the problem of simple linear regression that we used as an example in section I B 1 in the teacherstudent setting. Consider the teacher generated a number g ∗ from some probability distribution Ptp (g), this
g ∗ is the ground-truth value we aim to find back as precisely as possible from the noisy measurements. Then
the teacher generated heights yi independently from a Gaussian probability distribution of mean −g ∗ t2i /2 and
variance ∆∗ , where ti are the known measurement times, i = 1, . . . , N . He then handed the heights y and the
distribution Ptp (g), Ptm (yi |g) to the student. Set in this way the Bayesian solution presented in sec. I B 1 has
no ambiguity in terms of assuming the noise was Gaussian and independent, nor in terms of the choice of the
prior distribution for g.
• Historically the teacher-student setting was introduced for the problem of perceptron, called model B in [103].
In perceptron the aim is to learn weights Ji , i = 1, . . . , N in such a way that the scalar product with each of P
patterns ξiµ , µ = 1, . . . , P is constraint to be either smaller (yµ = −1) of larger (yµ = 1) than some threshold κ.
Mathematically stated one requires for all µ = 1, . . . , M
N
X
yµ = sign(
Ji ξiµ − κ) .

(12)

i=1

This perceptron setting models classification of patterns into two classes. One way to set a solvable version
of the model is to consider both ξiµ and yµ to be independent identically distributed (iid) random variables.
The goal is then to find a set of weights Ji such that the constraints (12) are satisfied. Another way to set
the problem that corresponds to the teacher-student scenario is to consider some teacher weights Ji∗ , generated
from some probability distribution, e.g. ±1 uniformly at random. The teacher perceptron then generates the
PN
output as yµ = sign( i=1 Ji∗ ξiµ − κ), where ξiµ are some known iid patterns. The goal is then to recover the
teacher weights Ji∗ from the knowledge of ξiµ and yµ . A very comprehensive presentation and solution of the
perceptron that uses this teacher-student terminology appears in [238].

10
3.

Bayes optimal versus mismatched

Throughout this manuscript we will be distinguishing two main cases depending on what information about the
distributions Ptp (x) and Ptm (y|x) the teacher gave to the student
• Bayes optimal case: In this case the teacher hands to the student the full and correct form of both the prior
distribution Ptp (x) and the likelihood Ptm (y|x).
• Mismatching prior and/or model: In this case the teacher hands no or only partial information about the prior
distribution Ptp (x) and/or the statistical model Ptm (y|x). In the present manuscript we will assume that the
teacher handed the correct functional form of the prior distribution and of the statistical model, but not the
values of the parameters θ and ∆, respectively.
Let us now consider the following exercise. Take x∗ to be the ground truth values of the variables generated by the
teacher and take x, x1 , x2 to be three independent samples from the posterior probability distribution P (x|y). We
then consider some function f (a, b) of two configurations of the variables a, b. Consider the following two expectations
Z
E [f (x1 , x2 )] = f (x1 , x2 )P (y)P (x1 |y)P (x2 |y)dx1 dx2 dy ,
(13)
Z
Z
E [f (x∗ , x)] = f (x∗ , x)P (x∗ , x)dx dx∗ = f (x∗ , x)P (x∗ , x, y)dx dx∗ dy
Z
= f (x∗ , x)P (x|y, x∗ )Ptm (y|x∗ )Ptp (x∗ )dx dx∗ dy .
(14)
where we used the Bayes formula. We further observe that P (x|y, x∗ ) = P (x|y) because x is independent of x∗ when
conditioned on y. Remarkably, in the Bayes optimal case, i.e. when P (x) = Ptp (x) and P (y|x) = Ptm (y|x), we then
obtain
E [f (x1 , x2 )] = E [f (x∗ , x)] ,

Bayes optimal :

(15)

meaning that under expectations there is no statistical difference between the ground truth assignment of variables
x∗ and an assignment sampled uniformly at random from the posterior probability distribution. This is a simple yet
immensely important property that will lead to numerous simplifications in the Bayes optimal case and it will be used
in several places of this manuscript, mostly under the name Nishimori condition.
In the case of mismatching prior or model the equality (15) typically does not hold.

4.

Estimators

What is the optimal estimator x̂ for the variables x? The answer naturally depends on what the quantity we aim
to optimize is. The most commonly considered estimators are:
Maximum a posteriori (MAP) estimator. MAP simply maximizes the posterior distribution and is given by
x̂MAP = argmax P (x|y) .

(16)

x

The main disadvantage of the MAP estimator is the lack of confidence interval and related overfitting issues.
Minimum mean squared error (MMSE). Ideally we would like to minimize the squared error between x̂ and x∗
SE(x̂, x∗ ) =

N
1 X
(x̂i − x∗i )2 .
N i=1

(17)

In general, however, we do not know the ground truth values of the variables x∗ . The best we can do within Bayesian
statistics is to assume that x∗ is distributed according to the posterior probability distribution. We then want to find
an estimator x̂ that minimizes the squared error averaged over the posterior
1
MSE(x̂) =
N

Z
P (x|y)

N
X
(x̂i − xi )2 dx .
i=1

(18)

11
By a simple derivative with respect to x̂i we see that the minimum of this quadratic function is achieved when
x̂MMSE
=
i

Z

Z
P (x|y)xi dx =

µi (xi ) xi dxi ,

(19)

where the right hand side is actually the mean of the marginal on the i-th variable defined by eq. (11).
The value of the MMSE is then computed as the MSE (18) evaluated at x̂i (19).
Maximum mean overlap (MMO). In cases where the support of the variables is discrete and the variables represent
indices of types rather than continuous values, the mean-squared error is not very meaningful. In that case we would
rather count how many of the N positions obtained the correct type. We define the overlap O as the fraction of
variables for which the estimator x̂ agrees with the ground truth
O(x̂, x∗ ) =

1 X
δx∗i ,x̂i .
N i

(20)

As in the case of MMSE, when we do not have the knowledge of the ground truth, the best Bayesian estimate of the
overlap is its average over the posterior distribution
MO(x̂) =

1
N

Z
P (x|y)

X

δxi ,x̂i dx .

(21)

i

The mean overlap MO maximized over the estimator x̂ leads to the maximum mean overlap estimator
x̂MMO
= argmax µi (xi ) ,
i

(22)

xi

where µ(xi ) is the marginal probability of variable i having type xi , defined by eq. (11). We should note that the
MMO is not an estimator very commonly considered in statistics, but given that statistical physics often deals with
discrete variables, and anticipating the close relation, the MMO will turn out instrumental.
Let us compare the MAP estimator with the ones based on marginals. In most cases the MAP estimator is easier
to approach algorithmically (from the computational complexity point of view optimization is generically easier than
enumeration). Moreover, in many traditional settings, where the number of samples is much higher than the number of
variables, there is usually very little difference between the two since the marginal probabilities are basically peaked
around the value corresponding to the MAP estimator. For instance, let us consider the teacher-student scenario
with unknown parameters θ of the teacher’s prior, or ∆ of the teacher’s model. In situations where the number of
samples M is large, but the parameters θ and ∆ are scalar (or low-dimensional), the posterior probability of these
parameters is so closely concentrated around its expectation that there is no practical difference between the MAP
and the marginals-based estimators for θ and ∆. When it comes to estimating the variables x then in the setting of
high-dimensional statistics (see Sec. I B 5), where the number of samples M is comparable to the number of variables
N the MAP estimator for x will generically provide crude overfitting and we will prefer estimators based on marginals.
Another reason to favor marginalization based estimators is because they come directly with a notion of statistical
significance of the estimate. Yet another reason is that from the perspective of spin glass theory there might even be
computational advantages (ground states might be glassy, whereas Bayes-optimal computations are not, see sec II G).
Let us now investigate the MMSE (and the MMO) estimator in the Bayes optimal case in the view of equality (15).
In section I B 3 we concluded that in the Bayes optimal case (i.e. when the teacher handed to the student the precise
form of the prior distribution and of the model likelihood) and in expectation, a random sample x from the posterior
distribution can be replaced by the ground truth assignment x∗ without changing the values of the quantity under
consideration. It hence follows that
Bayes optimal :

MMSE = E[SE(x̂, x∗ )] .

(23)

In words, the MMSE computed solely using the posterior distribution is on average equal to the squared error between
the MMSE estimator and the ground truth value of the variables x∗ . This is a very nice property. An analogous
identity holds between the MMO and the average of the overlap between the MMO estimator and the ground truth
assignment
Bayes optimal :

MMO = E[O(x̂, x∗ )] .

(24)

12
5.

The high-dimensional limit

Let us denote by N the number of variables x = {xi }i=1,...,N (referred to as dimensionality in statistics), and by M
the number of observations (or samples) y = {yi }i=1,...,M . Traditionally, statistical theory would consider the case of
a finite number N of variables/parameters to estimate and a diverging number of samples M . Think for instance of
the two examples (decaying particles and linear regression) from sec. I B 1.
In today’s data deluge more than ever, it is crucial to extract as much information as possible from available data.
From the information theoretic perspective, we would expect that useful information is extractable even when the
number of variables (dimensionality) N is comparably large (or even somewhat larger) than the number of observations
(samples) M . In that case, separating the useful information about x from noise is much more challenging. This is
referred to as the high-dimensional statistical theory. It is, in particular, in this high-dimensional framework, where
the amount of data to analyze is huge, that we need to pay attention not only to solvability but also to computational
efficiency, keeping in mind both questions from section I B.

FIG. 1. Schematic representation of a typical high-dimensional inference problem.

Typically, we will consider the limit where α = M/N is fixed, whereas N → ∞. In this limit the following scenario,
illustrated schematically in Fig. 1, very often applies: For low values of α < αc , successful inference of the variables
is not possible for any algorithm, the corresponding observed information is simply insufficient. For high values of
α > αs , computationally efficient inference algorithms do exist. In an intermediate regime αc < α < αs , successful
inference is in principal possible but algorithmically considerably harder than above αs . In many settings, the values
αc and αs can be defined in a precise mathematical way as values at which some measure of performance (or its
first derivative) presents a discontinuity in the limit of large size. In such a case, we talk about thresholds or phase
transitions (see section I C 4 for a more precise meaning of the second). In many problems the intermediate hard
phase is missing and αc = αs . On the other hand, settings and examples where the hard phase exists (αc < αs ) are
theoretically very interesting and challenging many researchers to prove some formal results about this hardness. For
a general inference problem, the values of αc and αs are not known, and neither are efficient algorithms that would
work down to αs . The state-of-the-art results usually only give a lower bound on αc , and an algorithmic upper bound
on αs .
The objective of a large part of the statistical physicists interested in these directions and included in this manuscript
could be summarized as follows (see also the schema Fig. 1): Determine the threshold αc between the inferenceimpossible and inference-possible region, and the threshold αs between the inference-hard and inference-easy region.
And develop efficient algorithms that would be able to infer the variables successfully with as little information
contained in the observations as corresponds to αs .

C.

Statistical physics and inference

The connection between statistical physics and statistical inference has a long history. Arguably, the most influential
early works are due to E. T. Jaynes, starting with his Statistical mechanics and information theory in 1957 [141] where
he develops the maximum entropy principle, thus building upon the pioneering work of C. Shannon [283], who named
the amount of information as entropy inspired by mathematically analogous formulas in thermodynamics. Jaynes
actually suggests that statistical mechanics does not have to be regarded as a physical theory dependent on its
assumptions, but simply as a special case of a more general theory of Bayesian inference. An inspiring reference is a
collection of lecture notes Statistical Physics to Statistical Inference and Back from a meeting held in Cargèse in 1992
[111]. The goal of that meeting was precisely to revise the foundational connections between the two fields. This close
connection leads to a continued cross-fertilization between the two fields. Important methodological developments in

13
one field find numerous applications in the other, and a new applicational area in one field profits greatly from the
methodology more common to the other.
The present manuscript is mainly concerned with recent applications of methods from spin glasses and glasses
[214] to inference problems and associated algorithms. This direction was perhaps most influenced by the work on
simulated annealing by Kirkpatrick, Gelatt, and Vecchi [164]. Connections between spin glasses and error correcting
codes [290], or neural networks such as the perceptron [281, 310] belong among the early very influential works. This
line of research is relatively well explored and illustrated on a number of problems as seen from classical books such
as the one by H. Nishimori Statistical Physics of Spin Glasses and Information Processing An Introduction [238], or
more recently by M. Mézard and A. Montanari, Information, Physics, and Computation [217].
Data processing is, however, a very fast developing field with many emerging yet far reaching applications. On top
of that, there was an important development in the theory of spin glasses on sparse random graphs in the early 2000s
[211, 212]. This manuscript aims at exposing some of this recent development.

1.

Basic dictionary

Let us first establish a dictionary between statistical inference and statistical physics. The posterior probability
(10) can be written as
)
( M
N
X
X
1
1
P (x|y) =
log P (xi ) =
log P (yµ |{xi }i∈∂µ ) + β
exp β
exp −βH(x, y) .
Z(y)
Z(y)
µ=1
i=1

(25)

We have used the exponential form to note that it takes the form of the Boltzmann probability distribution on variables
(discrete or continuous) with Z being the partition function, log P (xi ) the local magnetic field, and log P (yµ |{xi }i∈∂µ )
being the interaction between variables, H(x, y) being the Hamiltonian. We also introduced the auxiliary parameter
β, that was previously β = 1, and plays the role of the inverse temperature. The whole exponent divided by the
inverse temperature is then minus the Hamiltonian H of the physical system.
Note that the distribution (25) encompasses a wide range of models, the variables x can be discrete (e.g. spins)
or continuous, the interactions can depend on few or many variables (e.g. in the planted spin glass of section II the
interactions are pair-wise, and in the linear estimation of sec. VI each of the interactions involves all of the variables).
From a physics point of view the marginal probabilities (11) are nothing else than local magnetizations obtained by
summing the Boltzmann distribution over all but one variable. The MMSE (or the MMO) estimator is then computed
using these local magnetizations (at β = 1). On the other hand the MAP estimator can be viewed as the ground
state of the physical system, i.e. the minimum of the Hamiltonian. Another way to think about the MAP estimator
is as the zero temperature limit, β → ∞, of MMSE estimator.
2.

The teacher-student scenario defines the planted ensemble

Looking at the Boltzmann distribution (25) it could represent quite a generic statistical physics model where the
observations y play the role of a quenched disorder. Keeping in mind that we want to study Bayesian inference in
the teacher-student scenario of sec. I B 2 we have to emphasize that the disorder y is generated in a very particular
way (by the teacher) related to the form of the Boltzmann distribution itself. This becomes particularly striking in
the Bayes-optimal setting in which the teacher handed to the student all the information about how he generated the
data y except the realization x∗ .
In the rest of this manuscript we call models created in the teacher-student scenario of sec. I B 2 the planted models.
We call a disorder y that results from the teacher-student scenario the planted disorder. We contrast it with the
most commonly considered case of randomly-quenched disorder where the component yµ are independent identically
distributed variables. In the planted disorder the components yµ are not independent, they are correlated because
they all depend on the ground-truth x∗ that the teacher used to generate y. We call x∗ the planted configuration.
This is the same difference between solving a random system of linear equations y = Fx (which will have no solution
if F is a M × N matrix and M > N ) and solving one where the left hand side was generated as y = Fx∗ , instead of
just being random. In the second case, we have planted a solution x∗ in the problem.
A generic strategy to define the planted statistical ensemble is to take an inference problem in the teacher-student
setting. The Bayes-optimal posterior probability distribution of that inference problem then corresponds to the
Boltzmann distribution of the planted model at temperature T = 1/β = 1. What makes the planted model particular
is the form of the quenched disorder that is not independently random but instead correlated in a very particular way.

14
Perturbing a random ensemble by hiding a solution in it is often a natural thing to do. For instance, protein folding
is often modeled, in statistical physics, by a random energy model (to take into account its glassiness) in which one
special configuration (the so-called native configuration) is artificially created with a very low energy [46]. Thinking
of Bayesian inference in terms of the planted ensemble of statistical physics models is implicitly done in a number of
works. For instance the pioneering articles of Sourlas on finite-temperature decoding in error correcting codes [290–
292] and the relation to the Nishimori line, summarized excellently in [238], were very influential. In mathematics,
planting was introduced by Jerrum [142] as a model to study hardness of finding a large clique planted in a graph, and
also by Jerrum and Sorkin [143] as a model to study graph bisection. Planting was used in many mathematical works
either as a useful model or as a proof technique. It is the use of planting in proofs for constraint satisfaction problems
[4] that brought the authors’ attention to this concept. The authors find that making the connection between planting
and inference very explicit brings a lot of insight into a generic inference problem and for this reason we build the
present review around this connection.
To clarify further what is exactly the planted ensemble we give several examples and counter-examples:
• Planted coloring. Graph coloring is a famous mathematical problem where nodes of a graph are assigned colors
with a constraint that two connected nodes should not have the same color. Coloring of random graphs is a
constraint satisfaction problem, well studied in statistical physics [234, 320]. In planted coloring, the teacher
first generates a random assignment of colors to nodes x∗ , and then he generates edges of the the graph at
random but in a way that two nodes of the same color are not connected. The edges of the graph play the role
of the planted disorder y. The graph (i.e. its set of edges) is then handed to the student, but not the planted
configuration x∗ . The role of the student is to find a configuration as correlated to the planted one as possible.
Physics-study of the planted coloring appeared in [182], and we will see that it is a special case of the stochastic
block model to which we devote section V.
• Planted perceptron. In sec. I B 2 we presented the perceptron as the problem for which the teacher-student
terminology was established. Formulated as a planted problem, the set of weights J∗ the teacher-perceptron
uses is the planted configuration. The patterns ξiµ are in the physics works considered as iid, they play the role
of the ordinary randomly-quenched disorder. The teacher generates the outputs yµ using the patterns and the
planted weights J∗ (12). The output y then plays the role of the planted disorder. The planted weights J∗ are
then hidden from the student and her goal is to recover them as precisely as possible. The associated posterior
probability distribution is the Boltzmann measure of the planted perceptron. On the other hand, when asking
how many random patterns a perceptron can classify (the so-called capacity [103]), one is studying a non-planted
problem.
• Hopfield model [135] is a well known model studied in statistical physics. The goal in the Hopfield model is to
construct an Ising spin model for which a certain pre-defined set of randomly generated patterns ξiµ ∈ {±1},
i = 1, . . . , N and µ = 1, . . . , M , serve as configurations attracting local dynamics. It is, in fact, a model of
associative memory where, when starting close to one of the patterns, the dynamics would converge back to it.
PM
The Ising model studied by Hopfield has Hebb-like interactions Jij = µ=1 ξiµ ξjµ .
We could study the inference problem for recovering the patterns ξiµ from the knowledge on the Jij . For a single
pattern, M = 1, the posterior distribution of this inference problem would indeed correspond to the Hopfield
model, i.e. Ising spin model with interactions Jij . However, for more than one pattern M > 1 the posterior
of the inference problem will be a function of N M variables, and will be very different from the Boltzmann
distribution of the Hopfield model (which is just a function of only N spins). This shows that inferring or
remembering are very different tasks. Clearly, the associative memory studied in the Hopfield model, where one
start the dynamics close to one of the patterns, is different from inference of patterns from the Jij . The Hopfield
model is not a planted model in the sense to which we adhere in this review.
Since this is important we reiterate that the planted model is a statistical physics model where the Boltzmann
distribution is at the same time the posterior probability distribution of the corresponding inference problem. There
are many models considered in the literature where the quenched disorder is created using some auxiliary variables,
but the Boltzmann distribution is not the posterior probability for inference of these auxiliary variables. Then we do
not talk about a planted model, an example is given above using the Hopfield model with more than one pattern.
3.

Inference via Monte Carlo and variational mean-field

When it comes to solving an inference problem, i.e. evaluating the corresponding MMO of MMSE estimators, see
sec. I B 4, there are many algorithmic methods studied and used in computer science, statistics and machine learning.

15
For nice and comprehensive textbooks see e.g. [198, 309]. Two of the classical inference algorithms are particularly
interesting in the present context because they originate in statistical physics. These are the Gibbs sampling of
Monte Carlo Markov chain (MCMC) based methods [104, 105, 235] and the variational mean-field inference methods
[146, 244, 251, 308]. Both these methods and their applications are well covered in the existing literature and in this
review we will not focus on them. We only describe in the present section their basic concepts and ideas and relation
to the message-passing-based algorithms on which we focus in the manuscript.
The development of Monte Carlo goes back to the well know paper of Metropolis et al. [209] and was widely used
since in physics and other scientific fields. The Metropolis-Hastings algorithm [125, 209], most commonly used in
physics, is based on the comparison of energy cost in one state of a system and a small perturbation of the state
(typically the value of one variable is changed). When the cost decreases the change is accepted, when the cost
increases the change is accepted with a probability exponentially proportional to the cost difference, p = e−β∆E ,
where β is a inverse-temperature-like parameter. Another perturbation is introduced and the procedure is repeated.
Gibbs sampling [104, 105], known as heat-bath in physics, is a variant on the Metropolis-Hastings algorithm where
given a fixed state of all but one variable, the probability distribution of the one variable is evaluated conditioning on
the fixed state of the other variables. The value of the variables is then sampled from this probability distribution.
The process is then repeated. Both Metropolis-Hastings and Gibbs sampling are designed to sample the Boltzmann
probability measure. But they can also be famously turned into a generic purpose optimization algorithms called
simulated annealing [164].
The advantage of Monte Carlo based methods it that if iterated for sufficiently long time they provide the exact
answer for a very generic and broad class of systems (only ergodicity and balance condition are required). The
disadvantage is that the time needed to obtain the exact result can be as large as exponential in the size of the system
and it is in general very hard to decide what is the running time necessary for satisfactory precision. Nevertheless,
Monte Carlo based methods are very useful in many cases, in particular when the domain of variables is discrete, and
throughout this manuscript we will be comparing to them.
Tracing the precise origins of the use of variational mean-field for inference from data is a little harder. Ideas related
to variational mean-field were used in statistical physics since Gibbs (together with Maxwell and Boltzmann) created
the field, and are certainly present in the works of Van der Waals [303] and Pierre Weiss [313] arround the 1900’s.
The modern view on variational mean-field for inference and learning is very nicely summarized in [146, 198].
The basic idea of variational inference is to approximate the posterior probability distribution P (x|y) (25) by
another probability distribution Q(x|y) such that (a) for Q(x|y) the Bayesian estimators (MMO or MMSE) can be
evaluated tractably, and (b) the Kullback-Leibler divergence
DKL (Q||P ) =

X
x

Q(x|y) log

Q(x|y)
P (x|y)

(26)

between P and Q is the smallest possible. The variational inference methods are very fast and in many settings and
practical problems they provide precise approximations. However, in some settings they provide crude approximations
and sometimes misleading results (see an example in sec. V F).
The present manuscript focuses on message passing inference algorithms, such as belief propagation (BP). In physics
the origins of these methods trace back to the Bethe approximation [34], the work of Peierls [250] and the ThoulessAnderson-Palmer equations [300]. Interestingly, these equations can be seen as “improved” variational mean-field
equations, and this is at the roots of many perturbation approaches (such as Plefka/Yedidia-Georges [107, 253] or
Kikuchi [163]).
In information theory belief propagation was invented by Gallager [100] and in Bayesian inference by Pearl [249]. In
an excellent article Yedidia, Freeman and Weiss [318] relate belief propagation and the variational mean field methods.
However, strictly speaking belief propagation is not a variational method, because the related Bethe free energy does
not have to be larger than the true free energy. Compared to variational mean field methods belief propagation is
as fast, and generally more precise. An important property on which we focus in this manuscript, is that on a class
of models, that can be called mean-field spin glasses, belief propagation method gives asymptotically exact results
(whereas variational mean field does not) and it is amenable to exact analysis that enables detailed understanding of
many systems of interest.

4.

Useful statistical physics concepts

Traditional physical systems are composed of a number of spins or molecules that is comparable to the Avogadro
number NA ≈ 1023 . It is hence almost always relevant to investigate the limit of large system size that is in statistical
physics referred to as the thermodynamic limit.

16
Self-averaging. In the analogy between inference and statistical physics of disordered systems, the observations y
play the role of the so-called quenched disorder in the Hamiltonian in eq. (25). Generically, quantities of interest, e.g.
the magnetization of variable number 10, depend on the realization of this quenched disorder.
Statistical physics, however, mostly focuses on quantities that have the so-called self-averaging property, i.e. in
the thermodynamic limit their value does not depend on the realization of the disorder y (but only on its statistical
properties). More formally, a quantity A(y) is self-averaging if for every  > 0
lim Prob {|A(y) − E [A(y)] | > } = 0 .

N →∞

(27)

One of the most important quantities that we expect to be self-averaging is the free energy density
f (β, y, N ) = −

1
log Z(y)
Nβ

(28)

that is assumed to have a thermodynamic limit
f (β) = lim f (β, y, N )
N →∞

(29)

that does not depend on the realization of y.
Among the quantities discussed previously in this manuscript, the ones for which the self-averaging property is
particularly interesting are the MMSE and the MMO. Under self-averaging, the equalities (23-24) hold not only on
average, but also when the averages on the r.h.s. are removed. This means that the typical squared distance between
the ground truth and the MMSE (MMO) estimator is equal, in the thermodynamic limit, to the one evaluated without
the knowledge of the ground truth. Clearly, focusing on the analysis of self-averaging quantities greatly simplifies the
description of the behavior in the thermodynamic limit and enables a large part of the results obtained in statistical
physics.
It should be noted, however, that from a rigorous point of view, proving that the above quantities indeed are
self-averaging poses a non-trivial technical challenge. This is a recurring discrepancy between theoretical physics and
rigorous approaches, in physics we make assumptions that are correct and enable progress and understanding of the
problems, but that are not always easy to prove.
Phase transitions. Another concept borrowed from physics that will turn out very influential in the analysis of
high-dimensional inference are phase transitions. Intuitively, a phase transition is an abrupt change in behavior as a
parameter (commonly the temperature) is tuned. Since not everything that is termed phase transition in computer
science or engineering literature is a genuine phase transition from a physics point of view, it is useful to make some
reminders.
First of all, true phase transitions do not exist at finite system size. In mathematical physics a phase transition is
defined as a non-analyticity in the free energy density f (β) (29). Since the finite N free energy density is a logarithm
of a sum of exponentials, it is always an analytical function. Only when the limit N → ∞ is taken can the nonanalyticities appear. Therefore when we talk about a phase transition, we should always make sure that it indeed
corresponds to a non-analyticity and not to a smooth (even though perhaps very rapid) change.
Secondly, a traditional phase transition is always associated to either a critical slowing down or to metastability.
According to the Ehrenfest’s classification of phase transitions, there are two main types
• 1st order phase transitions are associated to a discontinuity in the first derivative of the free energy. The first
derivative of the free energy is usually related to some measurable quantity that consequently has a discontinuity.
Since a number of quantities, including e.g. the MMSE or the magnetization, are related to this first derivative,
they also display a discontinuity at a first order phase transition. The mechanism behind a 1st order phase
transition is a competition between two stable phases of the system. For the sake of vocabulary let us call the two
phases “paramagnet” and “ferromagnet”. Around the phase transition there is a region of phase coexistence
in which the properties of the system are ruled by properties of both the phases. In systems that cannot
be embedded into a (finite dimensional) Euclidean space we can define sharp boundaries of this parameterregion, called the spinodal points. The phase with lower free energy (i.e. higher log-likelihood) is said to be
thermodynamically stable, the other one is called metastable. Metastability has profound consequence on the
behavior of dynamics and algorithms. We will come back to this in detail.
• 2nd order phase transitions are associated to a discontinuity in the second derivative of the free energy. This type
of phase transition is associated to a so-called critical slowing down, meaning that dynamical processes become
increasingly slow next to the phase transition. Diverging timescales are closely linked to diverging length-scales
and diverging strength of fluctuations. These are described by the theory of criticality and critical exponents
that ruled statistical physics in the second half of last century.

17
Genuine phase transitions are quite peculiar creatures, and this is largely why they fascinate many physicists.
To make a distinction in vocabulary in physics we use cross-over for a more generic rapid change in behavior as a
parameter is tuned. We will keep the term phase transition or threshold for what is described above. Going back to
Figure 1, a physicist’s mind becomes very intrigued when for a given inference problem we find that, as the amount
of information in the measurements is increased, there exist genuine phase transitions. For a physicist, a phase
transition is something intrinsic, so we immediately expect a profound meaning and non-trivial consequences for the
system under investigation.

18
II.

PLANTED SPIN GLASS AS A PARADIGM OF STATISTICAL INFERENCE

As biology has its drosophila, or computational complexity its K-SAT problems, in statistical physics the canonical
example is the Ising model
H(S) = −

X
(ij)∈E

Jij Si Sj −

N
X

hi Si ,

(30)

i=1

where S = {Si }i=1,...,N and the spins are Si ∈ {−1, 1}, E is the set of interacting pairs, Jij is the interaction strength
between spins and hi is the magnetic field. Statistical physics is the study of the Boltzmann probability measure on
spin configurations
µ(S) =

1
e−βH(S) .
Z(β, J, h)

(31)

In this section we study an inference problem, that we call the planted spin glass, for which, in the teacher-student
setting of sec. I B 2, the posterior probability distribution is the Boltzmann distribution of the Ising model. The
purpose of this exercise is simply that we will be able to use a large part of the knowledge existing in the literature on
spin glasses to describe in detail the properties of this inference problem. Moreover, the physical picture that arises
in this section generalizes to other inference problems in the teacher-student setting. Once we understood in detail
the planted spin glass problem, we can basically deduce without calculation what will qualitatively happen in a wide
range of other inference problems. Whereas the problem presented in this section may appear as a very special and
relatively simple case, all the concepts that we illustrate on this example are very generic and apply to many other
highly non-trivial inference problems. Sections V and VI then give two recently studied examples.
The Ising model with the Hamiltonian (30) is presented in numerous textbooks. The most commonly considered
geometries of interactions are finite dimensional lattices, fully connected graphs, random graphs, or for instance social
networks. Interactions Jij are chosen such that Jij ≥ 0 in the ferromagnetic case (most commonly Jij = J for all
(ij) ∈ E) while in the antiferromagnetic case Jij ≤ 0. In the random field Ising model the fields hi are random. Perhaps
the most intriguing case of an Ising model is the spin glass, in which the interactions Jij are taken independently at
random from some distribution including both negative and positive numbers. In most of the situations considered
in physics literature, the interaction couplings J are chosen independently of any specific configuration S.
A.

Definition of the planted spin glass

To motivate the planted spin glass we consider the following problem: We have a (large) number N of people and
a pack of cards with the number +1 or −1 on each of them. There is about the same number of each. We deal one
card to every person. In the next step we choose pairs of people and ask them whether they have the same card or
not. Their instructions are to roll a dice and according to the result answer truthfully with probability ρ and falsely
with probability (1 − ρ). We collect the answers into a matrix Jij , Jij = 1 if the pair (ij) answered that they had the
same card, Jij = −1 if the answer was that they did not have the same card, and Jij = 0 if we did not ask this pair.
We want to model this system, and one of the questions we have is whether and how we can estimate which person
has which card.
Thinking back about the teacher-student scenario of sec. I B 2 the teacher distributes the cards, i.e. generates a
latent variable Si∗ = ±1, 1 ≤ i ≤ N for each of the nodes/people independently at random. We then select the graph
of observations. Let us assume that the pairs to whom we ask questions are chosen at random, and that we ask M of
them. The set of queried pairs E can then be represented by an Erdős-Rényi (ER) random graph [87, 88] on N nodes
having average degree c = 2M/N , as defined in section III A. Graphs created in this manner are at the core of many
results in this review. The teacher collects the answers into a matrix Jij = ±1, (ij) ∈ E taken from the following
distribution:
P (Jij |Si∗ , Sj∗ ) = ρδ(Jij − Si∗ Sj∗ ) + (1 − ρ)δ(Jij + Si∗ Sj∗ ) ,

(32)

where ρ is a bias of the dice (i.e. the probability of getting the correct answer), and δ(·) is the Dirac delta function.
The latent values Si∗ are then hidden from the student and the goal is to infer from the values of J as precisely as
possible which node had which card. Depending on the situation, the teacher did or did not reveal the value ρ he
used. For the moment we will assume ρ is known.
Let us look at the simplest cases: In the noiseless case, ρ = 1, each connected component of the graph has only
two possible configurations (related by global flip symmetry). When the graph is densely connected, and hence has

19
a giant component, one can therefore easily recover the two groups. When ρ = 1/2, it is as if we chose the couplings
independently at random just as in the most commonly considered spin glass model. The values J do not contain
any information about the latent variables Si∗ , and recovery is clearly impossible. What is happening in between, for
1/2 < ρ < 1? (The behavior for 0 < ρ < 1/2 is then deduced by symmetry.)
The above model is simple to state, but its solution is very interesting, and rather non-trivial. This problem
appeared in several different contexts in the literature, and is called censored block model in [2, 3, 274]. From a
rigorous point of view this problem is still a subject of current research. In the physics literature a closely related
model is known as the Mattis model [207]. Here, we will call the above problem the planted spin-glass model to
differentiate it from the standard spin-glass model, where we chose all the interactions J randomly and uncorrelated.
It should be noted, that “spin glass” here is used as a name for the model, it does not necessarily imply that this
model must exhibit glassy behavior. The phase diagram on the planted spin glass is discussed in detail in Sec. III D.
We stress that the main goal in inference is to retrieve the planted configuration, whereas the usual goal in physics of
spin glasses is to understand the glassy behavior of the model. These two goals are of course related in some aspects,
but also complementary in others.
More formally defined, to create an instance, that is a triple (S∗ , G, J), of the planted spin glass we proceed as
follows
• We consider that the hidden assignment Si∗ , called the planted configuration of spins, is chosen uniformly at
random from all the 2N different possibilities.
• We then choose a graph of interactions G = (V, E), but not yet the associated values of Jij , independently of
the planted configuration S∗ .
• For each edge (ij) ∈ E, the value of Jij is chosen to be either +1 or −1 with a probability taken from (32). The
signs of the couplings then carry information about the planted configuration.
So far the relation with the Ising model may not be very clear. To see it very explicitly we write eq. (32) differently.
Without losing generality, we can define a parameter β ∗ by
∗

eβ
ρ=
and
2 cosh β ∗

∗

e−β
1−ρ=
with
2 cosh β ∗

1
β = log
2
∗



ρ
1−ρ


,

(33)

so that eq. (32) reads
∗

P (Jij |Si∗ , Sj∗ ) =

∗

∗

eβ Jij Si Sj
.
2 cosh β ∗

(34)

Since we want to know the values of the spins S given the knowledge of the couplings J, we write the Bayes formula
Q
β ∗ Jij Si Sj
P (J|S)P (S)
(ij)∈E e
P (S|J) =
= N
P (J)
2 (2 cosh β ∗ )M P (J)

(35)

where we used (34). The posterior probability distribution has to be normalized, therefore
P (J) =

Z(β ∗ , J)
2N (2 cosh β ∗ )M

,

(36)

where Z(β ∗ , J) is the partition function of the Ising model (31) at inverse temperature β ∗ (in zero magnetic field). In
other words, the posterior distribution is simply the standard Boltzmann measure of the Ising model at a particular
inverse temperature β ∗ . In general we shall consider this Ising model at a generic inverse temperature β 6= β ∗ . This
is because a priori we do not know the proportion of errors. In the teacher-student setting, the teacher may not have
handed to the student the value β ∗ .
The special case β = β ∗ corresponds to the Bayes optimal inference as discussed in section I B 3. One of the most
physically important properties of the Bayes-optimal case is that the planted configuration can be exchanged (under
averages, or in the thermodynamic limit) for a random configuration taken uniformly from the posterior probability
measure, which is in physics simply called an equilibrium configuration. All equilibrium configurations with high
probability share their macroscopic properties, e.g. their energy, magnetization, various overlaps etc. This is directly
implied by the Nishimori condition (15) derived in Sec. I B 3 for the Bayes-optimal inference. This property will turn
out to be crucial in the present context.

20
As discussed in Sec. I B 4, since we are dealing with a problem with discrete variables (spins), it is desirable to
minimize the number of mismatches in the estimation of the planted assignment. We hence aim to evaluate the
maximum mean overlap (MMO) estimator (22) which in the case of an Ising model is simply
Ŝi = sign (mi ) ,

(37)

where mi is the equilibrium value of the magnetization of node i in the Ising model. We reduced the inference of the
planted configuration to the problem of computing the magnetizations in a disordered Ising model as a function of
temperature.
An off-the-shelf method for this that every physicist has in her handbag is a Monte-Carlo simulation, for instance
with the Metropolis algorithm. In Sec. III we shall see that for the mean-field setting (i.e. the underlying graph
being random or fully connected) there exists an exact solution of this problem via the cavity method. However,
before turning to the mean-field setting we want to discuss concepts that are more generic and we will hence keep
Monte-Carlo in our mind as a simple and universal method to evaluate phase diagrams of Hamiltonians with discrete
degrees of freedom.
B.

Nishimori’s mapping via gauge transformation

For the planted spin glass there is an alternative way that provides understanding of its phase diagram and quite
interesting additional insights [238–240]. One takes advantage of the so-called gauge transformation
Si Si∗ → S̃i ,

Jij Si∗ Sj∗ → J˜ij .

(38)

If this transformation is applied, the Ising Hamiltonian is conserved since all Si = ±1. In the new variables the
planted configuration becomes ferromagnetic: S̃i∗ = (Si∗ )2 = 1. What happened to the couplings? The energy of each
edge in the planted configuration has not been changed by the gauge transform. Since the frustrated interactions
were chosen independently at random in the planting, after the gauge trasformation we end up with a standard spin
glass with a fraction
∗

ρ=

eβ
2 cosh β ∗

(39)

of positive Jij = +1 couplings and the rest −1, where the −1’s are chosen independently at random.
After the gauge transformation, the planted assignment has been transformed to the ferromagnetic one where all
spins are positive. The overlap between the planted and another randomly chosen equilibrium configuration becomes
simply the magnetization. The system is now a standard spin glass with iid interactions, albeit with a ferromagnetic
bias ρ. The question of the identification of the hidden assignment is now simply mapped to the question of finding the
ferromagnetic “all spins up” configuration. This is a step that from the physics point of view simplifies the discussion.
It should be stressed, however, that the theory of Bayes-optimal inference, that is the main subject of the present
manuscript, applies to a class of planted models for which a similar gauge trasformation does not always exist.
The phase diagram of the spin glass models, with a ferromagnetic bias ρ, versus a temperature 1/β, has been
determined and one has simply to look it up in the literature [56, 57, 189]. The non-zero value of magnetization in
the ferromagnetically biased model then translates into the ability to be able to infer the planted configuration better
than by a random guessing. The value of the local magnetization of spin i is then related to the MMO estimator via
(37). In a remarkable contribution, Nishimori realized that the physics on the line ρ = eβ /(2 cosh β) is particularly
simple and quite different from a generic β. This line is now called the Nishimori line and corresponds to nothing
else but the condition for Bayes optimality β = β ∗ in the planted model.
Consequently, a number of interesting properties, that do not hold in general, hold on the Nishimori line. These
are called Nishimori conditions. In general, Nishimori conditions are properties that hold in the Bayes optimal case,
such as (23) and (24). Another useful Nishimori identity concerns the average energy E = hHi/N (30). For more
Nishimori conditions see [185, 238]. Using again the property that the planted configuration can replace an equilibrium
configuration we obtain immediately for the planted spin glass on Erdős-Rényi random graph (as defined in section
III A) of average degree c
c
E = − tanh β ∗ .
2

(40)

In Fig. 2 we give the phase diagram of the spin glass at temperature T = 1/β with ferromagnetic bias ρ. As an
example we use the spin glass on a random graph in which case an exact solution exists (see sec. III). However, on

21
a generic geometry the phase diagram could be obtained using Monte-Carlo simulation and would be qualitatively
similar. Due to the Nishimori mapping this is also the phase diagram of the planted spin glass, where the ferromagnetic
phase is understood as correlated to the planted configuration. The phenomenology of a rather generic Bayesian
inference problem can be directly read off from this phase diagram.
2

P

P

1.5

T=1/β

F
Tc
1

SG
0.5
F
M

Tc
TF
Nishimori line
0
0.5

0.6

0.7

0.8

ρc

0.9

1

β*

ρ=e /2cosh β*

FIG. 2. Phase diagram of the planted spin glass as a function of temperature T = 1/β and density ρ of edge interactions that
agree with product of spins in the planted configuration. This phase diagram is for random regular graphs of degree c = 3.
The Nishimori line, ρ = eβ /(2 cosh β) (red online), corresponds to the Bayes optimal inference. There is the spin glass (SG)
phase where the system is in a glassy state, the ferromagnetic (F) phase, the paramagnetic (P) phase and the mixed (M) one
where the system is both SG and F. The boundary of the mixed phase (black dashed) is only a guide for the eye between the
multi-critical point and the zero-temperature point from [58]. The same for the boundary of the SG phase (black dotted) that
is obtained from the inexact 1RSB approximation and is hardly distinguishable from the vertical line (black full).

C.

Bayes-optimality on the Nishimori line

First, let us assume that we know (or that we learn from data, see next section) the true value of β ∗ (or equivalently
ρ) for which the couplings J have been generated. Then, our consideration is restricted to the Nishimori line. If the
temperature is large enough, i.e. β ∗ is low enough, we are in the paramagnetic phase, and we find that all marginals
(magnetizations) are zero. It would thus be difficult to find their signs. In fact, this simply means that, in the
paramagnetic phase, inference of the hidden configuration is impossible. There is just not enough information to say
anything about the hidden/planted assignment.
As the temperature is lowered, i.e. β ∗ is increased, the situation changes. If we used β ∗ corresponding to an inverse
temperature above a certain critical value βc , then a Monte-Carlo simulation would rapidly discover the ferromagnetic
phase where local magnetizations are (mostly) positive. Applying rule (37) would give at this point a decent overlap
with the true assignment (which is, after the gauge transformation, all +1).
In terms of inference of the planted configuration these two phases have the following interpretation
• Undetectable: For β ∗ < βc , ρ < ρc , the planted spin glass system is in the paramagnetic state. Meaning that the
set of observed variables J does not contain any information about the planted configuration and its inference
is hence information-theoretically impossible.
• Detectable: For β ∗ > βc , ρ > ρc , the planted spin glass is in a ferromagnetic phase where equilibrium configurations are all correlated with the planted one and inference is possible, and algorithmically tractable e.g. via
Monte-Carlo sampling.
For the planted spin glass of a random graph, the critical temperature Tc = 1/βc where one goes from impossible
inference to easy inference can be taken from the literature [56, 57, 299] (we will show how to compute it for this case

22
in Sec. III) and in the case of random Erdős-Rényi graphs (as defined in section III A) with average degree c it reads

−1
1
TcER (c) = atanh √
,
c

(41)

−1

1
.
TcREG (c) = atanh √
c−1

(42)

while for regular random graphs it reads

On an Erdős-Rényi random graph (as defined in section III A), this indicates that it is only possible to infer something
about the planted assignment if the number of edges in the graph cN/2 is such that
c > cdetect =

1
.
(2ρ − 1)2

(43)

If this is not satisfied, the system is in the paramagnetic phase and nothing can be said about the planted configuration.
If, on the other hand, (43) is satisfied, then the magnetization aligns partly with the planted configuration. This is
a very generic situation in Bayes-optimal inference when we encounter a second order phase transition. We go from
a phase where inference is impossible, to a phase where it is rather straightforward. Note that this is not only a
statistical physics prediction, but one that can be proven rigorously. The impossibility result has been shown in
[126, 191] and the possibility one (if (43) is satisfied) in [274].
D.

Parameter mismatch and learning

Let us now discuss what happens if one does not know the value β ∗ . One approach that seems reasonable is to use
β → ∞, so that we look for the most probable assignment (this would correspond to the MAP estimator) which is the
ground state of the corresponding Ising model. We observe from the phase diagram in Fig. 2 that for ρ < ρc , at T = 0
the system is in the spin glass phase. In terms of Monte-Carlo, the presence of such a phase manifests itself via a
drastic growth of the equilibration time. Computing the corresponding magnetizations becomes very hard. Moreover,
even in the undetectable phase, ρ < ρc , the zero-temperature local magnetizations will not be zero. This is a sign of
overfitting that the MAP estimator suffers from.
In the detectable phase, ρ > ρc , the ferromagnetic phase extends only up to a temperature that for an Erdős-Rényi
random graph corresponds to a certain


1
βF = atanh
.
(44)
c tanh β ∗
Hence if we chose too low value of β as a starting point we could miss the detectability phase. At low temperatures
we might again encounter a spin glass phase that may or may not be correlated to the planted configuration. In either
case, such a phase causes equilibration problems and should better be avoided. Fortunately, there is no such glassy
phase on the Nishimori line as we shall explain in sec. II G.
The above considerations highlight the necessity of learning parameters (in this case β ∗ ). Within the theory
described above there is one particularly natural strategy for learning such parameters, closely related to the well
known algorithm of expectation maximization [74]. Let us illustrate this parameter learning on the example of the
planted Ising model where we assume that the parameter β ∗ is not known.
In Bayesian inference, we always aim to write a posterior probability distribution of what we do not know conditioned
to what we do know. In the present case we need P (β|J). Using the Bayes formula we obtain
P (β|J) =

X

P (S, β|J) =

S

Z(β, J)
P (β)
.
2N P (J) (2 cosh β)M

(45)

We notice that as long as the hyper-prior P (β) is independent of N , in the thermodynamic limit this probability
distribution P (β|J) converges to a Dirac delta function on the maximum of the second fraction on the right hand
side. This maximum is achieved at
*
+
1 X
cN
E=−
Jij Si Sj = −
tanh β .
(46)
N
2
(ij)∈E

23
The corresponding iterative algorithm for learning β is as follows. We start with some initial choice of β and compute
the average energy to update the value of β according to (46) until a fixed point is reached. Depending on the initial
value of β, and on the algorithm used to estimate the average energy, in some cases we may obtain several different
fixed points of (46). In those cases we choose the fixed point that is the global maximum of (45).
Notice now a very useful property, the condition for stationarity of (45) is the same as the Nishimori condition
(40) and indeed β = β ∗ is a fixed point of (46). This strategy is very reasonable, in the Bayes-optimal case all the
parameters are the ground truth ones and a number of Nishimori conditions hold only in this case. Hence in order to
learn the right values of parameters we iteratively impose the Nishimori conditions to hold.
Finally, let us note that in the above planted Ising model the underlying phase transition is of second order, therefore
there is no particular computational issue. Inference is either impossible, or possible and tractable. The hard phase
outlined in section (I B 5) is absent. The picture is richer in problems where in the random case there is a discontinuous
phase transition. Examples of such a case are described in detail in Sec. III E.
E.

Relations between planted, quenched and annealed disorders

From the early works, it sometimes seems necessary to have a gauge transformation (and consequently a special
form of the Hamiltonian) to be able to speak about the Nishimori conditions or other properties of the Nishimori line
[106, 239, 240]. Importantly this is not the case. The properties of the Nishimori line are more general and apply to
every case where it makes sense to talk about Bayes optimal inference. In fact, a more general notion of the Nishimori
line can be identified with the notion of Bayes-optimality. A nice article on this connection was given by Iba [137].
In this section, we shall discuss properties of the planted disorder, and its relation with the more commonly considered - the randomly-quenched and the annealed disorders. This clarifies the generic connection between statistical
physics of spin glasses and Bayesian inference. We stress again that none of this requires a gauge symmetry of the
underlying Hamiltonian.
Let us consider back the Hamiltonian of the Ising system (30). For a given system, e.g. the magnet on the fridge
in your house, the set of couplings Jij is a given fixed finite set of numbers. In the planted spin glass above, this set
was the values of interactions J we observed.
In statistical physics, where in the vast majority of cases we can take advantage of self-averaging, see sec. I C 4,
we often do not think of a given set of numbers Jij but instead about a probability distribution over such sets. For
instance randomly-quenched disorder is generated following
P (Jij ) =

1
1
δ(Jij − 1) + δ(Jij + 1)
2
2

(47)

for all (ij) ∈ E, where E are the edges of the graph of interactions. In this way we do not have a given Hamiltonian,
but rather its probability of occurrence, or equivalently an ensemble where instances appear proportionally to this
probability. Properties of the system then have to be averaged over the disorder (47). This view is adopted in
statistical physics because the self-averaging is assumed to hold in most physical systems of interest and the theory
is then easier to develop. This notion of averaging over disorder is so well engraved into a physicist’s mind that going
away from this framework and starting to use the related (replica and cavity) calculations for a single realization of
the disorder has lead to a major paradigm shift and to the development of revolutionary algorithms such as the survey
propagation [215]. In this section we do, however, want to stay in the mindset of averaging over the disorder and
explain how to fit in the notion of the planted ensemble, i.e. ensemble of instances where we average over the disorder
that was planted.
1.

The quenched and annealed average

How to take the average over disorder of Jij meaningfully is a problem that was solved a long time ago by Edwards
in his work on spin glasses and vulcanization (see the wonderful book Stealing the gold [109]). If one takes a large
enough system, he suggested, then the system becomes self-averaging: all extensive thermodynamic quantities have
the same values (in densities) for almost all realizations of the Hamiltonian. Therefore, one needs to compute the
average free energy
 
F
1
fquenched =
= lim −
[log Z] ,
(48)
N →∞
N
βN
where here [·] denotes the average over the disorder. This is called the quenched average, and the corresponding
computations are the quenched computations. In fact, the self-averaging hypothesis for the free energy has now been

24
proven rigorously in many cases, in particular for all lattices in finite dimension [312] and for some mean-field models
[115]. Edwards did not stop here and also suggested (and gave credit to Mark Kac for the original idea) a way to
compute the average of the logarithm of Z, known today as the replica trick, using the identity
Zn − 1
.
n→0
n

log Z = lim

(49)

The idea here is that while averaging the logarithm of Z is usually a difficult task, taking the average of Z n might
be feasible for any integer value of n. Performing a (risky) analytic continuation to n = 0, one might compute the
averaged free energy over the disorder as
1
[Z n ] − 1
lim
.
N →∞ N β n→0
n

fquenched = − lim

(50)

The computation of the quenched free energy is in general still very difficult and remains a problem at the core of
studies in statistical physics of disordered systems. It is much easier to consider the so-called annealed average, one
simply averages the partition sum and only then takes the logarithm
fannealed = −

1
log [Z].
Nβ

(51)

This is a very different computation, and of course there is no reason for it to be equal to the quenched computation.
For the example of the Ising model with binary uniformly distributed couplings (47) we obtain (independently of
the dimension or geometry)
Ising
− βfannealed
= log 2 +

c
log(cosh β) ,
2

(52)

where c is the average coordination number (degree) of the interaction graphs.
It is important to notice that the annealed average is wrong if one wants to do physics. The point is that the free
energy is an extensive quantity, so√ that the free energy per variable should be a quantity of O(1) with fluctuations
going in most situations as ∝ (1/ N ) (the exponent can be different, but the idea is that fluctuations are going to
zero as N → ∞). The partition sum Z, however, is exponentially large in N , and so its fluctuations can be huge.
The average can be easily dominated by rare, but large, fluctuations.
Consider for instance the situation where Z is exp (−βN ) with probability 1/N and exp (−2βN ) with probability
1 − 1/N . With high probability, if one picks up a large system, its free energy should be fquenched = 2, however in the
annealed computation one finds


1
1
[Z] =
exp (−βN ) + 1 −
exp (−2βN )
(53)
N
N
and to the leading order, the annealed free energy turns out to be fannealed = 1.
One should not throw away the annealed computation right away, as we shall see, it might be a good approximation
in some cases. Moreover, it turns out to be very convenient to prove theorems. Indeed, since the logarithm is a
concave function, the average of the logarithm is always smaller or equal to the logarithm of the average, so that
fannealed ≤ fquenched .

(54)

This is in fact a crucial property in the demonstrations of many results in the physics of disordered systems, and in
computer science this is the inequality behind the “first moment method” [227].
Furthermore, there is a reason why, in physics, one should sometimes consider the annealed average instead of the
quenched one. When the disorder is changing quickly in time, on timescales similar to those of configuration changes,
then we indeed need to average both over configurations and disorder and the annealed average is the correct physical
one. This is actually the origin of the name annealed and quenched averages.
2.

Properties of the planted ensemble

In sec. II A we described the planted spin glass, where we first considered the planted configuration S∗ and a graph
G = (V, E), and then we considered the probability over interactions P (J|S∗ ) eq. (34). The ensemble of instances
(S∗ , G, J) generated by the planting, i.e. by the three steps of Sec. II A, is called the planted ensemble.
From a physics point of view the two most important properties of the planted ensemble, the two golden rules to
remember, are:

25
• The planted configuration is an equilibrium configuration of the Hamiltonian derived from the posterior distribution. By this we mean that the planted configuration has exactly the same macroscopic properties (e.g. its
energy or magnetization) as any other configuration that is uniformly sampled from the posterior distribution.
Indeed, when we view planting as an inference problem, we can write the corresponding posterior probability distribution as in eq. (35). In section I B 3 we derived that in the Bayes optimal inference the planted
configuration behaves exactly in the same way as a configuration sampled from the posterior.
Given that generating an equilibrium configuration from a statistical physics model is in general a difficult task
(this is, after all, why the Monte-Carlo method has been invented!) this is certainly a very impressive property.
It actually means that, by creating a planted instance, we are generating at the same time an equilibrium
configuration for the problem, something that would normally require a long monte-carlo simulation.
• The realization of the disorder of the planted problem is not chosen uniformly, as in the randomly-quenched
ensemble (47), but instead each planted instance appears with a probability proportional to its partition sum.
This can be seen from eq. (36) where the probability distribution over J in the planted ensemble is related to
the partition function of the planted system Z(J). A more explicit formula can be obtained using the annealed
partition function (52)
Pplanted (J) =

Z(J)
,
ΛZannealed

(55)

where Λ is the number of all possible realizations of the disorder, Λ = 2M for the Ising model. Note the 1/Λ can
also be interpreted as the probability to generate a given set of couplings in the randomly-quenched ensemble
where the couplings are chosen uniformly at random.
The average over the planted ensemble might seem related to the annealed average because for instance the planted
energy (40) is always equal to the annealed average energy, which is easily derived from the annealed free energy
(52). However, the free energy of the planted ensemble is in general different from the annealed free energy (52), as
illustrated e.g. in Fig. 3.

3.

Quiet planting

Let us now imagine we have found a system where the (quenched) free energy is exactly equal to the annealed
free energy. Then according to eq. (55) we see that for such a system Pplanted (J) = 1/Λ = Pquenched (J), meaning
that generating instances from the planted ensemble is the same thing as generating from the randomly-quenched
ensemble. Such planting is denoted as quiet. In fact, we do not even need the free energies to be exactly equal, it is
sufficient that the free energy densities are the same. This is because atypical instances are usually exponentially rare
and hence any difference in the free energies that does not show up in the leading order will not generate atypical
instances.
It so happens that in mean-field systems the (self-averaging) free energies of paramagnets are indeed equal to the
annealed free energies. The equilibrium in these systems can hence be studied using the planted ensemble and this is
greatly advantageous as we will see in section II H 1.
The idea of quiet planting comes from a rigorous work [4] where the above notion was proven rigorously. In
mathematics this property is formalized and the two ensembles termed contiguous [4, 231]. The paper that pioneered
the usage of quiet planting in physics, and that also established the name “quiet”, is Hiding quiet solutions in random
constraint satisfaction problems [182].

F.

Replica symmetry breaking

The replica and cavity method, that are among the main assets that statistical physics offers for analysis of
inference problems, were originally developed to understand behavior of glasses and spin glasses via the so-called
replica symmetry breaking (RSB). We will explain in the next section, that for Bayes-optimal inference RSB is not
needed. But to appreciate fully the results exposed in this manuscript it is useful to have a very basic idea of what
RSB is and what it implies.
RSB is usually explained together with the related (relatively involved) calculations, see e.g. [217, 238]. In this
section we try to convey the basic ideas behind RSB without going into the computations. We rely on the descriptions
of what would happen to a Monte Carlo simulation in a system that exhibits various versions of RSB.

26
0.25

0.2

splanted
0.15

s planted−cluster

0.1

s quenched

s annealed

0.05

Σ
0

12.6

12.8

13

Cd

13.2

Cc

13.4

13.6

13.8

14

Ccol
C

FIG. 3. Entropy density (which is equal to the free energy at zero temperature) for the 5-coloring problem on an Erdős-Rényi
random graph (as defined in section III A) of average degree c. The planted entropy (in red) is compared to the quenched
one (in green) and the annealed one (in black). For c < cc the three curves agree. The quenched entropy goes to −∞ at the
colorability threshold ccol . The region between cd and cc is relevant only for cases with a first order phase transition as discussed
in Sec. III E. In mean field theory of spin glasses this region is called the d1RSB phase [41]. We also show the entropy of the
planted cluster below cc (black dashed), which is less than the total entropy, the difference is called the “complexity” Σ (in
blue) and is the logarithm of the number of pure states. Bottom: Sketch of the shape of the space of all valid colorings in the
planted case. At an average degree cd the space of solutions shatters into exponentially many clusters, the planted cluster being
one of them. Beyond cc the planted cluster contains more solutions than all the others together. At cs the last non-planted
cluster disappears (this is the coloring threshold in a non planted, purely random model). Figure taken from [182].

Let us start with the concept of replica symmetry. For the purpose of the present manuscript we can say that a
probability measure P (x|y) is replica symmetric if a MCMC that is initialized at equilibrium, i.e. in a configuration
sampled uniformly at random from P (x|y), is able to sample P (x|y) (close to) uniformly in a number of iterations
that is linear in the size of the system. A probability measure P (x|y) exhibits dynamical one-step of replica symmetry
breaking (d1RSB) if MCMC initialized at equilibrium explores in linear time only configurations that belong into an
exponentially small (in N ) fraction of the equilibrium configurations. In both RS and d1RSB the distance between
two configurations sampled from the Boltzmann probability measure is the same number with probability going to
one (as N → ∞).
We speak about static one-step replica symmetry breaking (1RSB) if MCMC initialized at equilibrium is able to
explore configurations belonging to a finite fraction of all equilibrium configurations. Moreover the distance between
two randomly sampled configurations can take (with probability going to one) one of two possible values.
A probability measure P (x|y) corresponds to full-step replica symmetry breaking (FRSB) if the distance between two randomly chosen configurations is distributed according to some probability distribution with continuous
non-trivial support. In terms of behavior of MCMC, in a FRSB regime the time needed to sample configurations
uniformly starting from equilibrium has a more complex behavior: relaxations are critical and a power-law behavior
is expected [289].

G.

No replica symmetry breaking in Bayes-optimal inference.

From the methodological point of view, the Nishimori line has one absolutely crucial property: the absence of glass
phase at equilibrium. In the words of the mean-field theory of spin glasses, there is no static replica symmetry breaking
on the Nishimori line. Note that there might be, and indeed is in all the cases where the underlying transition is of
the 1st order, the dynamical one-step RSB phase (d1RSB) [41]. The marginals in the d1RSB phase are, however, still

27
exactly described by the belief propagation algorithm (as described in section III), and this is what matters to us in
the Bayes-optimal inference.
A replica symmetry breaking phase (or equivalently the static spin glass phase) can be defined by the non-selfaveraging of the overlap q between two configurations randomly sampled from the Boltzmann distribution. Depending
on the realization of the sampling, the Hamming distance between the two configurations is with nonzero probability
different from its average.
In section I B 3 we derived that, on the Nishimori line, the overlap m between a typical configuration and the planted
one is equal to the overlap between two typical configurations q. With similar arguments one can derive equality of the
distribution (over realizations of the sampling) P (q) of the overlaps between two configurations randomly sampled from
the posterior, and the distribution of the magnetizations P (m), which is the overlap between the planted configuration
and a random sample from the posterior [238]. In physics, the magnetization is always argued to be self-averaging
(the same is true for the free energy density, and the magnetization is its derivative with respect to the magnetic
field), and hence on the Nishimori line the overlap is also self-averaging. From this, one argues that P (q) is trivial,
which in mean-field spin glasses indicates the absence of a spin glass phase. However, from a rigorous point of view,
self-averaging of the magnetization remains an open problem.
Another way to see the problem is to remark that, given that the spin glass susceptibility is equal to the ferromagnetic
one for planted systems (see for instance [184, 185]), then it is hard to reconcile the presence of an equilibrium
ferromagnet phase (where the susceptibility is finite) with an equilibrium spin glass phase (where the susceptibility
diverges).
A complementary (rigorous) argument was presented by Montanari [220] who showed that in sparse systems (i.e.
each observations depends on a bounded number of variables) in the Bayes-optimal setting, two point correlations
decay. In the static spin glass phase these correlations could not decay.
This means that we can exclude the presence of a static equilibrium transition to a glass phase on the Nishimori
line. Note, however, that the dynamics can still be complicated, and in fact as we mentioned, a dynamic spin glass
(d1RSB) phase can appear. As soon as we deviate from the Bayes optimal setting and mismatch the prior, the model
or their parameters, the stage for glassiness to arrive is open. And indeed as we see in Fig. 2 there are regions with
RSB when β 6= β ∗ , typically when we take the temperature too small. When studying inference with a mismatched
prior or model we thus always have to keep this in mind before making claims about exactness.

H.

Applications of planting for studies of glasses

Now that we have explored the concept of planting and quiet planting, and its connection to inference, we will show
that it also has a number of consequences for the study of structural glasses in physics. While these are not directly
relevant to inference, it illustrates in many ways the fascinating properties of the planted ensemble and the deep
connection between inference, optimization problems and the statistical physics of glasses, and we thus summarize
these properties in this section.

1.

Equilibration for free in glasses

Let us discuss the physical consequences of the fact that in quietly planted systems the planted configuration
behaves exactly the same way as any other equilibrium configuration. We can hence literarily equilibrate for free.
This is particularly interesting in the study of glassy models. As a matter of fact, the very concept of glassiness can
be defined as super-fast growth of the equilibration time as the glass transition temperature is approached [41].
The core of the interest in glasses is concerned by finite (usually three) dimensional models. In those cases the
annealed and quenched free energies are known not to coincide. A class of theories of glasses is, however, based on
mean-field (i.e. systems living on random geometries of fully connected graphs) disordered systems that present a
discontinuous phase transition. The analog of cd in these systems is called the dynamical glass transition, and the
analog of cc is called the Kauzmann glass temperature [32]. Properties of the equilibration time are well studied and
it turns out that the equilibration time (i.e. number of Monte Carlo steps per variable in the thermodynamic limit)
diverges as the dynamical temperature is approached [222] and is exponentially large in the system size between the
dynamical and Kauzmann temperature. This in practice means that numerical simulations of these models are very
cumbersome.
In most of the mean-field models used for the studies of glasses, the quiet planting works. Numerical simulations
studying the dynamics from equilibrium can hence be speeded up considerably by initializing in the planted configuration. A particular version from of this trick was first used in [221] where it enabled for the first time precise

28
simulation of the dynamics starting from equilibrium. It was later popularized by [182] and has now been used in a
considerable number of works. Without being exhaustive let us mention Refs. [12, 60, 61, 93, 202, 203].
Even if one does not care about the equivalence to the randomly-quenched ensemble one can use planting to obtain
equilibrium configurations even in general finite dimensional systems. Indeed the property of planted configuration
being the equilibrium one for a corresponding Boltzmann measure is general. This has been used in numerical studies
of so-called pinning in glasses, see e.g. [132].

2.

Following states and analyzing slow annealing

In a very influential work, Franz and Parisi [96] considered the so called “potential” method, that aimed to analyse
the properties of a single pure state in glasses. Perhaps more interestingly, this allows to describe the adiabatic
evolution of these glassy Gibbs states as an external parameter, such as the temperature, is slowly tuned. Again,
there is a deep connection with the planted ensemble which we shall now discuss. Ref. [182] studied the state that is
created in the planted coloring around the planted configuration. Adding to this the concept of quiet planting, this
framework allows us to study the shape of states in mean-field glassy systems and is equivalent to the Franz-Parisi
construction.
When a macroscopic system is in a given phase, and if one tunes a parameter, say the temperature, very slowly
then all observables, such as the energy or the magnetization in a magnet, will be given by the equilibrium equation
of state. In a glassy system, where there could be many such states, it is very interesting to see how they evolve
upon cooling, as they tend to fall off equilibrium very fast. This is for instance what is happening in the so-called
3-XOR-SAT problem, also called 3-spins spin-glass [114, 216, 265]. This is nothing but an Ising spin glass where the
two body interaction of Hamiltonian (30) is replaced by a three spin interaction (i.e. each term is a triplet Jijk Si Sj Sk ).
A plot that illustrates this phenomena is given in Fig. 4. The blue line is the equilibrium (in technical terms, the
one-step replica symmetry breaking solution) energy of the 3-regular hyper-graph (this means that every variables
has exactly 3 triplets of interactions, chosen randomly) 3-XOR-SAT problem as a function of temperature T . For
this system quiet planting works at all temperatures. The red lines are energies that we obtain when we plant
quietly (at temperature where the red lines crosses the blue one) and then lower or increase the temperature and
monitor the corresponding energy given by the dynamics. This procedure has been called state following in [183, 323],
because physically this is a study of the corresponding state at temperature at which the state is no longer among
the equilibrium ones.
Note that such state following was often described in physics of glasses on the level of a thought-experiment.
Methods for quantitative analysis existed only for the particular case of spherical p-spin model, via an exact solution
of the dynamics [41, 67, 68]. For other systems, in particular diluted ones, state following was considered to be an open
problem, see e.g. the concept suggested but not analyzed in [172]. The method for state following is mathematically
very much related to the concept of Franz-Parisi potential at two different temperatures [22, 96].
From an algorithmic point of view, the method of state following opens the exciting possibility of analyzing limiting
energy of a very slow annealing. It is well known that in a finite size system annealing always finds the ground
state, in a time that is in general exponential in the system size. A crucially interesting question is: What is the
lowest achievable energy when the annealing rate scales linearly with the size of the system? In mean-field models of
glasses equilibration is possible in a time that is linear in the system size down to the dynamical temperature Td (or
equivalently up to cd ), after that the dynamics gets blocked in one of the exponentially many states. If we take an
equilibrium state at Td −  and follow it down to zero temperature the resulting energy is a very good candidate for
the behavior and for the limiting energy of very slow simulated annealing. Very slow annealing for the Ising p-spin
model was compared to the theoretical prediction from state following in [186] .
Paradox in state following. There is, however, a problem that appears when one tries to carry out the calculation
for following states that are at equilibrium close to Td . When following the state down in temperature, the replica
symmetry gets broken at some point (the boundary between full and dashed lines in Fig. 4). Moreover, at yet lower
temperature the solution correlated to the planted configuration disappears (in a mathematically analogous way as
when a spinodal line is encountered). Even using the 1RSB computation to describe the internal structure of the
state we were not able to follow it down to zero temperature. A possible explanation is that the correct description of
this region is the full-step replica symmetry breaking (FRSB). To test this hypothesis Ref. [294] analyzed the mixed
spherical p-spin model where the FRSB is tractable.
The findings of [294], however, raised even more questions. It is found by the Nishimori mapping to the model
with a ferromagnetic bias (as described in section II B) that there is no FRSB solution to the problem, yet there
is no magnetized 1RSB solution at very low temperature either. Yet in MCMC simulation one always confirms the
expectation that as temperature is lowered the bottom of the state is correlated to the reference configuration. In
other words, using the Nishimori mapping, the non-magnetized spin glass phase that is assumed to exist at low

29

FIG. 4. State following in regular 3-XOR-SAT. Figure taken from [323].

temperatures for ρ > ρc in Fig. 2 is not physical in the context of state following. These observations are hence an
inconsistency in the theory. Loose ends such as this one are always interesting and basically always eventually lead
to new theoretical development. This state following paradox should hence not be forgotten and should be revisited
as often as new possibilities open.
One such possibility might be related to the calculation involving a chain of reference (planted) configurations at
lower and lower temperatures as suggested in [97].
3.

Melting and glass transition

Using the planted models and extended form of the Nishimori conditions the analogy between glassy dynamics and
melting was studied in [184, 185]. Perhaps the theoretically most intriguing and very fundamental question in the
theory of glasses is: Is there a true glass transition in a finite-dimensional system (i.e. system where the graph of
interactions can be embedded into finite-dimensional Euclidean space)? Such an ideal glass transition has to come
with diverging time and length scales when approaching the transition. A definitive answer to this question is difficult
to obtain as both simulations and experiments are faced with the extremely slow dynamics. According to the random
first order theory of the glass transition [168, 169], there exist systems with time and length scales with a genuine
divergence at the ideal glass transition temperature. Moreover this theory is able to predict in which systems such
an ideal glass transition exists and in which it does not. However, the random first order theory is not free from
criticisms, see for instance [36] and references therein.
The question of existence of the ideal glass transition remains open, but in [185] it was derived that if there is a
finite-dimensional model with a first order phase transition on the Nishimori line (or equivalently in the Bayes optimal
inference) then there is an ideal glass transition as well. This result was reached by using the Nishimori mapping
from a planted model to a ferromagnet, sec. II B, and comparing their dynamical properties. In particular, if there is
a first order phase transition on the Nishimori line then it comes with two quite uncommon properties - absent latent
heat, and the melting dynamics being equivalent to the stationary equilibrium dynamics. Arguably, first order phase
transitions are easier to describe or exclude than the tricky glass transition. This is hence an interesting path towards
the fundamental question of existence of the ideal glass transition.

30
III.

PHASE DIAGRAM FROM THE CAVITY METHOD

Now we move to special cases of inference problems where the phase diagram and phase transitions can be computed
analytically exactly, and for which related algorithms more efficient than Monte-Carlo can be derived. We shall achieve
this goal using extended mean-field methods. Slightly abusively we will call mean-field systems those where these
methods lead to exact results.
Mean-field methods are often the first step in the understanding of systems in physics, and disordered systems
are no exception. In spin glasses, in particular, the mean-field model proposed in 1975 by D. Sherrington and S.
Kirkpatrick [285] has been particularly important. We refer to the classical literature [214] for the replica solution,
involving the so-called replica symmetry breaking, devised by Parisi. It is not that we are not interested in the replica
method, quite the contrary, the replica method is at the roots of many results in this review. We shall, however, focus
here in the alternative cavity method [211–213] since this leads more naturally to specific algorithms and algorithmic
ideas. The development of more efficient algorithms is the goal in inference problems.
There are two main types of lattices for which the mean-field theory is exact. The first type is the fully connected
lattice underlying the canonical Sherrington–Kirkpatrick model. Similar solutions have been derived for the p-spin
model [114] and for the Potts glass model [112], and these played a major role in the development of the mean-field
theory of the structural glass transition [165–167]. The linear estimation problem that we treat in Sec. VI is an
example of a problem on a fully connected lattice.
The second type of lattice for which the mean-field theory is exact is given by large random graphs with constant
average degree, a case commonly refereed to as “Bethe lattice” in the physics literature [304]. In 2001 Mézard and
Parisi [211, 212], using the cavity method, adapted the replica symmetry breaking scheme to solve models on such
sparse lattices. There is a number of reasons why Bethe lattices are interesting, the two major ones are: (a) because
of the finite degree, they provide a better approximation of the finite dimensional case and because the notions of
distance and neighboring can be naturally defined, and (b) the related random graphs are fundamental in many
interdisciplinary cases, for instance in computer science problems, which is our main motivation here.
The cavity method is a generalization of the Bethe and Onsager ideas [33, 242] for disordered systems. It was initially developed by Mézard, Parisi and Virasoro [213] as a tool to recover the solution of the Sherrington-Kirkpatrick
model without the use of replicas. It was subsequently developed by Mézard and Parisi to deal with the statistical
physics of disordered systems on random graphs [211]. There is a deep global connection between the Bethe approximation and the so-called belief propagation (BP) approach in computer science (in error correction [99] and Bayesian
networks [249]) as was realized in connection to error correcting codes by Kabashima and Saad [154], and was put
in a more general setting by [318]. For a more recent textbook that covers in detail the cavity method and belief
propagation see [217], for a very comprehensive explanation of belief propagation see [318].

A.

Random graphs

We will remind basic properties of sparse Erdős-Rényi (ER) [87, 88] random graphs, that we used briefly already
in sec. II A. At the core of the cavity method is the fact that such random graphs locally look like trees, i.e. there are
no short cycles going trough a typical node.
An ER random graph is taken uniformly at random from the ensemble, denoted G(N, M ), of graphs that have
N vertices and M edges. To create such a graph, one has simply to add random M edges to an empty graph.
Alternatively, one can also define the so called G(N, p) ensemble where an edge exists independently for each pair of
nodes with a given probability 0 < c/N < 1. The two ensembles are asymptotically equivalent in the large N limit,
when M = c(N − 1)/2. The constant c is called the average degree. We denote by ci the degree of a node i, i.e. the
number of nodes to which i is connected. The degrees are distributed according to Poisson distribution, with average
c.
Alternatively, one can also construct the so-called regular random graphs from the ensemble R(N, c) with N vertices
but where the degree of each vertex is fixed to be exactly c. This means that the number of edges is also fixed to
M = cN/2.
We will mainly be interested in the sparse graph limit, where N → ∞, c ∼ O(1), M ∼ N . The key point is that,
in this limit, such random graphs can be considered locally as trees [37, 139]. The intuitive argument for this result
is the following one: starting from a random site, and moving following the edges, in ` steps c` sites will be reached.
In order to have a loop, we thus need c` ∼ N to be able to come back on the initial site, and this gives ` ∼ log(N ).

31
B.

Belief propagation and the cavity method

Now that we have shown that random graphs are locally tree-like, we are ready to discuss the cavity method. We
will consider a system with Hamiltonian
X
X
H=−
Jij Si Sj −
hi Si
(56)
i

(ij)∈E

and derive the thermodynamics of this model on a tree. We have to deal with N variables Si with i = 1, . . . , N , and
M interactions, that we denote ψij (Si , Sj ), on each edge (ij) of the graph G(V, E), with ψij (Si , Sj ) = exp(βJij Si Sj ).
Further denoting ψi (Si ) = exp(βhi Si ) the “prior” information the Hamiltonian gives about Si , the total partition
sum reads


X
Y
Y
 ψi (Si )
Z=
ψij (Si , Sj )
(57)
i

{Si=1,...,N }

(ij)∈E

To compute Z on a tree, the trick is to consider instead the variable Zi→j (Si ), for each two adjacent sites i and j,
defined as the partial partition function for the sub-tree rooted at i, when excluding the branch directed towards j,
with a fixed value Si of the spin variable on the site i. We also need to introduce Zi (Si ), the partition function of the
entire complete tree when, again, the variable i is fixed to a value Si . On a tree, these intermediate variables can be
exactly computed according to the following recursion
!
Y
X
Zi→j (Si ) = ψi (Si )
Zk→i (Sk )ψik (Si , Sk ) ,
(58)
Sk

k∈∂i\j



Y X

Zi (Si ) = ψi (Si )
Zj→i (Sj )ψij (Si , Sj ) ,
j∈∂i

(59)

Sj

where ∂i \ j denotes the set of all the neighbors of i, except spin j. In order to write these equations, the only
assumption that has been made was that, for all k 6= k 0 ∈ ∂i \ j, Zk→i (Sk ) and Zk0 →i (Sk0 ) are independent. On a
tree, this is obviously true: since there are no loops, the sites k and k 0 are connected only through i and we have
“cut” this interaction when considering the partial quantities. This recursion is very similar, in spirit, to the standard
transfer matrix method for a one-dimensional chain.
In practice, however, it turns out that working with partition functions (that is, numbers that can be exponentially
large in the system size) is somehow impractical, and we can thus normalize eq. (58) and rewrite these recursions in
terms of probabilities. Denoting ηi→j (Si ) as the marginal probability distribution of the variable Si when the edge
(ij) has been removed, we have
Zi (Si )
ηi (Si ) = P
.
0
S 0 Zi (S )

Zi→j (Si )
,
ηi→j (Si ) = P
0
S 0 Zi→j (S )

(60)

So that the recursions (58-59) now read
!
X
ψi (Si ) Y
ηi→j (Si ) =
ηk→i (Sk )ψik (Si , Sk ) ,
zi→j
Sk
k∈∂i\j


Y
X
ψi (Si )

ηi (Si ) =
ηj→i (Sj )ψij (Si , Sj ) ,
zi
j∈∂i

(61)

(62)

Sj

where the zi→j and zi are normalization constants defined by:
!
zi→j =

X

ψi (Si )

Si

zi =

X
Si

ψi (Si )

Y

X

k∈∂i\j

Sk

Y
j∈∂i

ηk→i (Sk )ψik (Si , Sk )



X

ηj→i (Sj )ψij (Si , Sj ) .
Sj

,

(63)

(64)

32
The iterative equations (61,62), and their normalization (63,64), are called the belief propagation equations. Indeed,
since ηi→j (Si ) is the distribution of the variable Si when the edge to variable j is absent, it is convenient to interpret
it as the “belief” of the probability of Si in absence of j. It is also called a “cavity” probability since it is derived
by removing one node from the graph. The belief propagation equations are used to define the belief propagation
algorithm
• Initialize the cavity messages (or “beliefs”) ηi→j (Si ) randomly or following a prior information ψi (Si ) if we have
one.
• Update the messages in a random order following the belief propagation recursion eq. (61,62) until their
convergence to their fixed point.
• After convergence, use the beliefs to compute the complete marginal probability distribution ηi (Si ) for each
variable. This is the belief propagation estimate on the marginal probability distribution for variable i.
Using the resulting
marginal distributions, one can compute, for instance, the equilibrium local magnetization via
P
mi = hSi i = Si ηi (Si )Si , or basically any other local quantity of interest.
At this point, since we have switched from partial partition sums to partial marginals, the astute reader could
complain that it seems that we have lost out prime objective: the computation of the partition function. Fortunately,
one can compute it from the knowledge of the marginal distributions. To do so, it is first useful to define the following
quantity for every edge (ij):
zij =

X

ηj→i (Sj )ηi→j (Si )ψij (Si , Sj ) =

Si ,Sj

zi
zj
=
,
zj→i
zi→j

(65)

where the last two equalities are obtained by plugging (61) into the first equality and realizing that it almost gives
eq. (64). Using again eqs. (61-64), we obtain


X
Y X

zi =
ψi (Si )
ηj→i (Sj )ψij (Si , Sj )
Si

Sj

j∈∂i


=

X

ψi (Si )

Si

Y

X


j∈∂i

Sj


P
Zj→i (Sj )
Si Zi (Si )

P
P
ψij (Si , Sj ) = Q
,
0
S 0 Zj→i (S )
j∈∂i
Sj Zj→i (Sj )

(66)

and along the same steps
P

Zj→i (Sj )
P
.
k∈∂j\i
Sk Zk→j (Sk )

zj→i = Q

Sj

(67)

P
For any spin Si , the total partition function can be obtained using Z = Si Zi (Si ). We can thus start from an
arbitrary spin i




X
Y X
Y
Y X

zj→i
Z=
Zi (Si ) = zi
Zj→i (Sj ) = zi
Zk→j (Sk ) ,
(68)
Si

j∈∂i

Sj

k∈∂j\i Sk

j∈∂i

and we continue to iterate this relation until we reach the leaves of the tree. Using Eq. (65), we obtain




Q
Y
Y
Y zj Y zk
zi
zj→i

··· = Q i
.
Z = zi
zk→j · · ·  = zi
zij
zjk
(ij) zij
j∈∂i

k∈∂j\i

j∈∂i

(69)

k∈∂j\i

We thus obtain the expression of the free energy in a convenient form, that can be computed directly from the
knowledge of the cavity messages, often called the Bethe free energy
X
X
f N = −T log Z =
fi −
fij ,
i
(ij)
(70)
fi = −T log zi , fij = −T log zij ,

33
where fi is a “site term” coming from the normalization of the marginal distribution of site i, and is related to the
change in Z when the site i (and the corresponding edges) is added to the system. fij is an “edge” term that can be
interpreted as the change in Z when the edge (ij) is added. This provides a convenient interpretation of the Bethe
free energy eq. (70): it is the sum of the free energy fi for all sites but, since we have counted each edge twice we
correct this by subtracting fij (see for instance [212]).
We have now entirely solved the problem on a tree. There is, however, nothing that prevents us from applying the
same strategy on any graph. Indeed the algorithm we have described is well defined on any graph, but we are not
assured that it gives exact results nor that it will converge. Using these equations on graphs with loops is sometimes
referred to as “loopy belief propagation” in Bayesian inference literature [318].
When do we expect belief propagation to be correct? As we have have discussed, random graphs are locally tree-like:
they are trees up to any finite distance. Further assuming that we are in a pure thermodynamic state, we expect that
we have short range correlations, so that the large O(log N ) loops should not matter in a large enough system, and
these equations should provide a correct description of the model.
Clearly, this must be a good approach for describing a system in a paramagnetic phase, or even a system with a
ferromagnetic transition (where we should expect to have two different fixed points of the iterations). It could be,
however, that there exists a huge number of fixed points for these equations: how to deal with this situation? Should
they all correspond to a given pure state? Fortunately, we do not have such worries, as the situation we just described
is the one arising when there is a glass transition. In this case, one needs to use the cavity method in conjunction with
the so-called “replica symmetry breaking” approach as was done by Mézard and Parisi [211, 212, 215]. In this review,
however, we discuss optimal Bayesian inference and we are guaranteed from the Nishimori conditions discussed in
sec. II G that no such replica symmetry breaking occurs. We can assume safely that in the Bayes-optimal setting
the simple belief propagation approach is always providing the correct description (unless there are other sources of
correlations in the system, e.g. the graph is not random but a lattice or the quenched disorder is correlated in a
non-local manner).
C.

Properties of the randomly-quenched ensemble

The first system on which we will illustrate how to use the cavity method and the related belief propagation to
compute the phase diagram is the Viana-Bray (V-B) model for the so-called diluted spin glass [304]. The reason we
discuss first the solution for this well known model is that, using the theory developed in sec. II, we will be able to
read off the solution for the planted spin glass model from the phase diagram of the randomly-quenched model, see
Fig. 2. Viana and Bray studied the Ising spin glass Hamiltonian (30) with the lattice being an Erdős-Rényi random
graph with average degree c. Interactions Jij are chosen from ±1 independently and uniformly at random and, for
simplicity, we consider the absence of a magnetic field. We call the ensemble of instances generated this way the
randomly-quenched ensemble.
The description of the V-B model within this approach was originally done by [44]. Instead of working with the
iteration (61), it is convenient to rewrite it in terms of an “effective” local cavity fields ui→j defined as
ηi→j (Si ) =

eui→j Si
.
2 cosh eui→j

(71)

With a bit of algebra, the belief propagation equation (61) can be written as



X
1
 ,
ui→j
= f (Jij , {uk→i
atanh tanh(βJij ) tanh β
uk→i
t
t−1
t−1 }k∈∂i\j ) =
β

(72)

k∈∂i\j

where we have added an index t that stands for a “time” iteration. As stated, iterating this procedure allows to actually
solve the problem on a tree, which is at the root of the cavity method. Here, however, we apply this procedure on the
graph, iteratively, hence the time indices in eq. (72). The local magnetization is then computed from the BP fixed
point using the hyperbolic tangent of the local field acting on the spin, and thus reads
X
mi = tanh (β
uj→i ) .
(73)
j∈∂i

Applying the same notation to the Bethe free energy per site (70) leads to
"
#
X
X
X


i
SβJij
j→i
i→j
βN f =
(ci − 1) log 2 cosh(βh ) −
log
2e
cosh(βh
+ Sβh
) ,
i

(ij)∈E

S=±1

(74)

34
P
P
where hi = k∈∂i uk→i and hi→j = k∈∂i\j uk→i , and ci is the degree of node i. Physical properties of the system
are computed from the fixed points of the BP equations (72). At this point, there are two different strategies to solve
the problem. The first one is to investigate numerically what is happening on a given realization of the problem, i.e.
to generate a large random graph and to actually run the BP iteration on it. This is actually very instructive (as we
shall see) but no matter how large our graph will be, it will not be at the thermodynamic limit.
The second strategy is to work on an infinite graph directly via a so called population dynamics [211]. We attempt
to describe the model by the knowledge of the probability distribution P(u) of the messages u on an edge with value
J over the (infinite) graph. Consider for instance the case of a random graph with fixed degree ci = c. Clearly, within
a pure Gibbs states, P(u) satisfies the BP equations (72) so that
P(u) =

1
Z

Z

P (J)dJδ[u − f (J, {ud }d=1,...,c−1 )]

c−1
Y
d=1

P(ud )dud

(75)

where f (.) is given by eq. (72). While this appears to be a difficult equation to solve explicitly, it can be conveniently
approximated to arbitrary precision with the so-called population dynamics method [211]. In this case, one models
the distribution P(u) by a population of M values
M

Ppop (u) =

1 X
δ(u − ul ) .
M

(76)

l=1

Plugging this expression in eq. (75) leads to the population dynamics algorithm:
• Initialize all values ul randomly.
• For a given element l of the population, draw c − 1 values of u from the population, and a value J of the coupling
from P (J), compute a new value unew = f (J, {u}), and replace the element ul by unew .
• Repeat the last step until Ppop (u) has converged to a good approximation of P(u).
This method is very simple and, given the population dynamics distribution, allows to compute the value of the free
energy per site via (70).
Either way, depending on the inverse temperature β in the Viana-Bray model, one observes one of the two following
possibilities: At low β there is a stable paramagnetic fixed point ui→j = 0 for all (ij) ∈ E. The corresponding
paramagnetic free energy reads
− βf (β, J) = log 2 +

c
[log (2 cosh β) − log 2] .
2

(77)

We do notice that the same expression is obtained for the annealed free energy for this model.
At high β (low temperature) the iterations of (72) do not converge on a single graph. The last statement, of course,
can not be checked simply when performing the population dynamics algorithm. A similar effect can, however, be
observed by performing two-populations dynamics with the same random numbers starting from infinitesimally close
initial populations. At high β the two populations will diverge to entirely different values [246]. This is a clear sign of
long range correlations and of the appearance of a spin glass phase. We are, of course, interested to know the precise
location of the phase transition.
Thouless [299] analyzed the threshold value of βc by linearizing the update (72) around the uniform fixed point,
thereby obtaining
X
ui→j
= tanh(βJij )
uk→i
(78)
t
t−1 .
k∈∂i\j

Thouless then argues that averaged over the disorder Jij the above equation gives zero, and when its square is averaged
we get
hu2t i = c tanh2 β hu2t−1 i ,

(79)

where we remind that c is the average degree of the ER graph, that is also equal to the average branching factor of
the tree that approximates locally the random graph. We thus expect that the paramagnetic phase becomes unstable
when c tanh2 β ≥ 1.

35
In fact, by analyzing such perturbations, one is investigating the behavior of the susceptibility. Indeed the spin
glass susceptibility can be written as
χSG ≡ β

∞
X
X
hS0 Si i2c ≈ β
cr hS0 Sr i2c .

(80)

r=0

i

In eq. (80), · · · represents a spatial average over the whole graph, h· · ·i represents a thermal average, Sr are all spins
at distance r from S0 , and Xc is the connected version of the correlation function X: hS0 Si ic = hS0 Si i − hS0 ihSi i.
The correlation can be computed using derivatives with respect to an external magnetic field on the spin S0

βhS0 Si i2c

=

∂hSi i
∂hh0 i

2
.

(81)

Since the physical field is a function of the cavity fields, we can monitor the propagation of the response using the
chain rule
∂hSr i
∂hSr i ∂hhr i ∂hhr−1 i
∂hh1 i
=
···
∂hhc0 i
∂hhr i ∂hhr−1 i ∂hhr−2 i
∂hh0 i

(82)

where we indicate with hi the cavity field going from spin i to spin i + 1 for simplicity and hr , hr−1 , ..., h0 is the
sequence of cavity fields on the shortest path from S0 to Sr (we are sure that there exists just one shortest path
because we are locally on a tree). To check whether or not the susceptibility defined in eq. (80) diverges, one should
thus simply analyze the square of the derivative of the belief propagation equations, as we have done in eq. (79). The
conclusion is that the critical temperature is given by the condition that a perturbation would grow, so that


1
ER
√
βc = atanh
.
(83)
c
Given the branching ratio is c − 1 for a random regular graph, we get in this case:


1
.
βcreg = atanh √
c−1

(84)

It follows that for β < βc the state of the Viana-Bray model is correctly described by the paramagnetic fixed point.
For β > βc the Bethe approximation fails and for the correct physical description the so-called replica symmetry
breaking framework needs to be applied [211]. Since, in this review, we are concerned by the Bayes-optimal inference
that effectively lives on the Nishimori line and does not need the concept of replica symmetric breaking, we refer the
reader interested in this to the textbook [217].
D.

Phase diagram of the planted spin glass

Now we compute the phase diagram of the planted spin glass, where the interactions J are chosen conditionally on
the planted configuration S∗ , following eq. (34). In this case the Bethe approximation and its linearization are still
described by (72) and (78). In the planted model the stability analysis of the paramagnetic fixed point is slightly
different, since the distribution from which the interactions Jij are taken is not random but depends on the planted
configuration. This is where our mapping to the Nishimori ensemble is making the problem rather easy to analyze
anyway. Indeed, according to our analysis in the last section, we simply need to ask whenever, on the Nishimori
∗
eβ
line with a bias ρ = 2 cosh
β ∗ , the system acquires a spontaneous magnetization. Or more simply, we just need to
compute the phase diagram when a ferromagnetic bias is present. Taking the average of (78) and counting properly
the probability of Jij given Si∗ and Sj∗ we immediately see that, while the spin glass transition is unchanged, the
ferromagnetic one is now given by
hut i = c(tanh β ∗ )(tanh β) hut−1 i .

(85)

According to this computation, the system thus starts to polarize towards the ferromagnetic configuration when


1
β > βF = atanh
.
(86)
c tanh β ∗

36
For the Bayes-optimal case (or equivalently, on the Nishimori line) when β ∗ = β, this immediately gives the transition
from the no-detectability phase to the recovery phase that we gave in eq. (43). We see now how easy it was to compute
it using the cavity method and the Nishimori mapping.
A couple of subtle but interesting points can be made about the phase diagram, always keeping in mind the fact
that the inference problem can be seen as looking for the ferromagnetic state on the Nishimori line.
• As β ∗ → ∞ (ρ = 1) we recover the ordinary ferromagnetic Ising model with its well known critical temperature.
√
• It is interesting to observe that there a multi-critical point at β = atanh(1/ c), where all the three conditions
(83), (86) and Bayes-optimality, β = β ∗ , meet. This point, of course, is on the Nishimori line, and is precisely
the transition point eq. (43).
• For the Bayes-optimal case, β ∗ = β, for all temperature higher than the transition point β ∗ < βc , the paramagnetic fixed point is the only fixed point of the recursion (72). Since the paramagnetic and annealed free
energies are equal, this is a case where the planting is quiet, as defined in section II E 3. This means that in
the whole region β ∗ < βc the planted ensemble behaves exactly the same as the randomly-quenched ensemble.
This is indeed reflected in the phase diagram in Fig. 2, where everything on the left of the multi-critical point
is independent of the value of ρ.
• Out of the Bayes-optimal case, i.e. when β ∗ 6= β, there exists at very low temperature (large β) a so-called
mixed phase where equilibrium configurations are correlated to the planted one (the ferromagnetic after our
mapping) but where the equilibrium phase is actually a spin glass one.
The phase diagram of the planted spin glass is summarized in Fig. 2 where the boundary between the paramagnet
(P) and the spin glass (SG) is plotted in green, between the paramagnet (P) and ferromagnet (F) in blue, and the
Nishimori line corresponding to the Bayes-optimal inference in red. We also sketch the position of the mixed phase
(M).
E.

Systems with first order phase transitions

In the above planted Ising spin glass model the underlying phase transition is of second order. In that case there
is no particular computational issue. Either inference is impossible, or it is possible and tractable. The hard phase
outlined in section (I B 5) is absent.
The picture is richer in problems when a discontinuous phase transition is found in the randomly-quenched ensemble.
Famously, this is the case in the low density parity check error correcting codes as described e.g. in [217, 238], and
has been analyzed also in details in the pioneering work of Nishimori and Wong [241]. The examples of clustering of
networks and compressed sensing (to which we dedicate sections V and VI) also fall in this category.
In order to observe such a phenomenon, we have to change the Hamiltonian eq. (30) and move from two-body
interactions to three-body interactions. This leads to the so-called 3-spin spin glass model, which we already mentioned
in sec. II H 2 under the name 3-XOR-SAT [114, 216, 265]. The Hamiltonian now reads
X
H(S) = −
Jijk Si Sj Sk ,
(87)
(ijk)∈E

Again, we can create this model on a random hyper-graph. For instance, we can consider a L regular hyper-graph
where each variable belongs to exactly L triplets, chosen randomly at random. Again we can study this model under
a planted setting, or equivalently, on the Nishimori line. What is changing with respect to the preceding section?
First of all, we need to rewrite the belief propagation equation. It is an exercise (see for instance [216, 265] or the
appendix A in [323]) to see that they now read


X
Y
mi→a = tanh β
Jb
mj→b 
(88)
b∈∂i\a

j∈∂b\i

where a, b denote the triplets of interactions. Repeating the analysis of the previous section, one can obtain the phase
diagram which we show in Fig. 5.
At first sight, this phase diagram looks really similar to the one from Fig. 2. What has changed now is the fact
that the transition is a first order one. Let us first concentrate on the Nishimori line. As soon as T < Tspin there
are two fixed points of BP, a paramagnetic and a ferromagnetic one, and the one dominating the measure is the

37

2

Tspin.

P

T=1/β

1.5
TF

Tdyn.
1
TSG

0.5

SG
F

0
0.5

0.6

0.7
0.8
ρ=eβ*/2cosh β*

0.9

1

FIG. 5.
Phase diagram of the 3-XOR-SAT (or 3-spin glass) problem as a function of temperature T = 1/β and density
of ferromagnetic bonds ρ. This phase diagram is for random regular hyper-graphs of degree L = 5. The Nishimori line,
ρ = eβ /(2 cosh β), corresponds to the Bayes-optimal inference of a planted configuration. As in Fig. 2 there is the spin glass
(SG) phase where the system is in a glassy state, the ferromagnetic (F) phase, and the paramagnetic (P) phase (for simplicity,
the mixed phase is not shown). A major difference with respect to Fig. 2, however, is the fact that the ferromagnetic transition
is now a first order one, with a spinodal line. The spin glass transition is also a first-order one, with the so-called dynamic
transition happening before the static transition.

one with lower free energy. Only when T < TF does the free energy of the ferromagnetic fixed point become lower
than the paramagnetic one. This means that, for T < TF , a perfect sampling would give a configuration with a
positive magnetization. Using the inference-Nishimori connection, this means that inference is possible for T < TF ,
and impossible for T > TF : the equilibrium phase transition marks the information theoretic transition.
Now, let us try to actually perform inference in practice, for instance by running belief propagation or a MonteCarlo algorithm. Unfortunately, we have no idea about what is the planted configuration so at best we can initialize
randomly. The problem is that this will put us in the paramagnetic state, which exists and is locally stable at all
positive temperatures. Both BP or MCMC will simply explore this paramagnetic state, and, in any sub-exponential
time, will miss the existence of the ferromagnetic one. This is the classical picture of a first-order phase transition in
a mean-field model. Because of the first order transition, we are trapped into a high free energy minima that is stable
at all temperatures, and therefore, we cannot perform inference. While we should be able to find a configuration
correlated with the planted one (in the sense that it is information-theoretically possible), we cannot do it within the
class of existing polynomial algorithms. Related ideas, in the context of the glass transition, were originally present
in [95]. This illustrates the difference between a “possible” inference and an “easy” one in the sense of section I B 5.
XOR-SAT is actually a little bit special in the sense that the easy detectable phase (see below) is missing. There
is a spinodal transition (Tspin ) where the ferromagnetic phase becomes unstable, but not a second one where the
paramagnet become unstable. In most models, however, if the temperature is lowered a second spinodal is reached,
the paramagnetic phase becomes eventually unstable and inference is possible.
Let us thus discuss another example where such a transition arises. For this, we turn to the planted coloring
problem [182]. In planted coloring one first chooses a random assignment of q colors out of the q N possible ones.
Then, one puts M = cN/2 edges at random among all the edges that do not have the same color on their two
adjacent nodes. The Bayes optimal inference in this setting evaluates marginals over the set of all valid colorings of
the created graph. Computation of these marginals builds on previous works [234, 320], where the problem of coloring
random graphs was studied, and goes as follows: One considers belief propagation initialized in two different ways.
In the random initialization the messages are chosen at random but close to the paramagnetic ones. In the planted
initialization the initial values of the messages are set to point in the direction of the planted configurations. A fixed
point is reached from both these initializations and its free energy is evaluated. One obtains the following picture,
illustrated in Figs. 3 and 6 for planted 5-coloring
• Undetectable paramagnetic phase, for average degree c < cd . The planting is quiet and belief propagation
converges to the uniform fixed point both from random and planted initialization. We will call cd the dynamical
or the clustering phase transition [179, 182].
• Undetectable clustered phase, or dynamic one-step replica symmetry breaking (d1RSB) phase, for average degree
cd < c < cc . The planting is still quiet, but the set of solutions is split into exponentially many clusters that
are uncorrelated and all look the same as the one containing the planted solutions. We will refer to cc as the
detectability phase transition. In the randomly-quenched ensemble cc is the condensation transition [179, 182].

1

0.5

0.8

0.4
overlap, planted init
overlap, random init
fpara-fplanted

0.6

0.3

0.4

0.2
cd

cs

0.2

fpara-fplanted

normalized overlap

38

0.1

0

0
cc
12

13

14

15
c

16

17

18

FIG. 6. The phase transitions for planted 5-coloring. The overlap (normalized in such a way that for a random configuration
it is zero) achieved by BP initialized in the random and planted configuration. The difference between the corresponding free
energies is also plotted. Figure taken from [71].

• Hard detectable phase, for average degree cc < c < cs . This is arguably the most peculiar phase of the four. The
Bayes optimal marginals are correlated to the planted configuration, however, there is a metastable paramagnetic
phase to which belief propagation converges from a random initialization. A ghost of the quiet planting is visible
even for c > cc in the sense that the space of configurations looks exactly like the one of the randomly-quenched
ensemble except for the planted cluster.
• Easy detectable phase, for average degree c > cs . Belief propagation finds a configuration correlated to the
planted coloring. The phase transition cs can be located doing a local stability analysis of the uniform BP fixed
point, along the same lines as done for the planted spin glass in section III D. We will refer to cs as the spinodal
or the hard-easy phase transition. In the randomly-quenched ensemble, cs is the point starting from which belief
propagation stops converging, related to the local instability towards the spin glass phase [182, 320].
A careful reader will recognize that this corresponds to the situation we have described in the introduction when
we draw figure 1. The situation in this figure is that of a problem with a first order transition, where the possible
inference is defined by the equilibrium transition, and the easy/hard transition is defined by the presence of a spinodal
marking the instability of the paramagnet towards the planted state. This is arising in learning problems as diverses
as the planted perceptron, community detection with 4 or more communities, compressed sensing, error correction,
etc. If, however, the transition is a second order one, then these two transition coincides as in the case of the planted
spin glass (or, as we shall see, community detection with two equal sizes communities). We find it remarkable that
the century old classification of phase transition in physics, introduced by Paul Ehrenfest in 1933 [138], turns out to
be the right tool to discuss a deep phenomena of statistical inference and computational complexity.

F.

Comments on the phase transitions in various problems

One more question we want to clarify is whether the concept of quiet planting works for every inference problem.
Here the answer is both yes and no. First of all, for quiet planting to work, we need to be in a setting where belief
propagation has the uniform (sometimes confusingly called factorized) fixed point, where the messages in the fixed
point do not depend on the indices of the edge. If such a fixed point exists then the annealed free energy density equals
the quenched one [228], and consequently the argument about planting being quiet works. There are many problems
where BP on the randomly-quenched ensemble does not have a uniform fixed point, the most studied examples are
perhaps the problem of K-satisfiability (K-SAT), or the Ising model in a random external magnetic field hi .
Even in problems such as K-SAT, where canonical BP does not have a uniform fixed point [218], we can introduce
a reweighted Boltzmann measure, as done in [178], in such a way that the corresponding BP fixed point is uniform.
Planting according to this reweighted measure then leads to instances that are contiguous to the random ones. Note

39
that this reweighting and corresponding planted instances were used in the proof of [5, 6] and studied also in [23, 144].
The planted solution is not an equilibrium solution with respect to the canonical Boltzmann measure but with respect
to the reweighted measure.
Another comment is related to the following natural question: How to distinguish between problems where the
transition is continuous or discontinuous? And when it is discontinuous, how wide is the hard phase? Is there some
criteria simpler than performing the whole analysis? To distinguish between the first and second order transition, we
do not know of any other way. We only have a couple of rules of thumb - in Potts model with sufficient number of
states continuous transitions become discontinuous, and in models where more than two variables interact, e.g. the
p-spin model or K-SAT, the continuous transition becomes discontinuous for sufficiently large K and p [179]. But
there are binary pairwise problems such as the independent set problem where at sufficiently large average degree the
continuous transition also becomes discontinuous [17].
Interestingly, one can identify a class of constraint satisfaction problems (CSP) where the hard phase is very wide.
In [178] these problems were identified as those with cs → ∞. Previously, only XOR-SAT was known to have this
property [95], but in [89] a whole hierarchy of such problems was identified and is related to k-transitivity of the
constraints. A constraint is k-transitive if and only if the number of satisfying assignments compatible with any given
assignment of any k variables is the same. For all CSPs that are at least 2-transitive the hard phase extends over the
whole region of linearly many constraints, equivalently cs → ∞. Calling r the smallest number such that a constraint is
not r-transitive, the number of constraints predicted to render the problem tractable scales as N r/2 [89]. Interestingly,
this classification makes the K-XOR-SAT the hardest of all CSPs (despite the fact that it can be solved by Gaussian
elimination, which is not robust to any kind of noise). The notion of 2-transitivity, and consequently infinitely large
hard phase, extends also to continuous variables, and as it was recently pointed out, renders tensor factorization
and completion very difficult [266]. Physically, the corresponding systems are such when the paramagnetic phase is
locally stable down to zero temperature [95]. Examples of such systems indeed include the Ising p-spin (related to
K-XOR-SAT) or the spherical p-spin model related to tensor factorization.
Having examples where the hard phase is huge is very useful when thinking about formal reasons of why this phase
is hard. In [89] the authors proved hardness for a class of the so-called statistical algorithms using precisely this kind
of constraint satisfaction problems. Arguably, the class of statistical algorithms is too restrictive, but we believe that
the instances of CSP with such a large hard phase should be instrumental for other attempts in this direction. We
shall summarize in section IV B arguments for why we think there is something fundamentally hard about this phase.
In [324] we discussed the example of planted locked constraint satisfaction problems (CSP). These are particularly
interesting for a couple of reasons. First of all, the valid solutions are separated by extensive Hamming distance from
each other and therefore BP either converges to a paramagnet or points exactly towards the planted configuration.
Note that low-density parity check (LDPC) error correcting codes have precisely this property. Locked CSP can be
thought of as non-linear generalization of LDPC. A second reason why locked CSPs are interesting is that in some
cases the hard phase is very wide and moreover the planted assignment is with high probability the only one (perhaps
up to a global flip symmetry). The planted locked CSP are hence attractive candidates for a one-way function that
could find applications in cryptography. Having hard CSP instances that have exactly one valid assignment was used
in studies of quantum annealing where having non-degenerate ground state is greatly advantageous, see e.g. [127].

40
IV.

FROM PHYSICS INSIGHT TO NEW ALGORITHMIC IDEAS

In the introduction, sec. I B, we presented two main questions about inference that we aim to answer. One more
fundamental, about phase diagrams, phase transitions and related insights, the second more practical, about the
development of algorithms. This section is about algorithmic contributions obtained using the insight resulting from
the statistical physics analysis. For clarity, all is again illustrated on the planted spin glass model.
With the work of Mézard, Parisi and Zecchina on survey propagation [215], the statistical physics community
learned that the replica and cavity calculations are as powerful when used as algorithms as they are when used as
tools for analyzing the phase diagram. The connection between the Bethe approximation and belief propagation was
noticed earlier in connection to error correcting codes by Kabashima and Saad [154], and in a more general setting
by [318]. But looking at the volume of related literature, it is only after survey propagation that the whole statistical
physics community started to look for explicit algorithmic applications of the Bethe approximation.
Here we review some recent physics-related contributions to algorithmic development for Bayes-optimal inference.
A.

Non-backtracking spectral method

Notice, within the historical perspective, that the Bethe solution for the spin glass on a random graph was written
by Bowman and Levin [44], i.e. about 16 years before realizing that the same equation is called belief propagation
and that it can be used as an iterative algorithm in a number of applications.
The paper by Thouless [299], who analyzed the spin glass critical temperature, implicitly includes another very
interesting algorithm, that can be summarized by equation (78) where BP is linearized. The fact that this equation
has far-reaching implications when viewed as a basis for a spectral algorithm waited for its discovery even longer.
Notably until 2013 when the methodologically same algorithm was suggested for clustering of networks [180]. This
work will be discussed in sec. V D. We now turn back to our pedagogical example, the ±J planted spin glass, and
related spectral algorithm as described and analyzed in [274]. Related ideas also appeared in [325].
First note that a large class of inference algorithms is based on a spectral method of some kind. Probably the most
widely used is the so-called principal component analysis (PCA). Its operation can be thought of as revealing the
internal structure of the data in a way that best explains the variance in the data. In PCA the information is contained
in singular vectors (or eigenvectors for symmetric matrices) corresponding to the largest singular values. Resulting
spectral algorithms are used everyday to solve numerous real world problems. The same strategy of computing
leading singular (or eigen-) vectors of some matrix related to the data is applied in many settings including the
planted partitioning problem closely related to the planted spin glass model [64, 65], or in understanding financial
data [42, 190].
A natural question that arises is whether a spectral algorithm can find a configuration correlated to the planted
assignment in the planted spin glass in the whole region β ∗ > βc where belief propagation can asymptotically do so.
Facing this question, many experts would probably suggest to take the matrix Jij (with Jij = 0 when (ij) ∈
/ E) or
its associated Laplacian (perhaps normalized in some way) and compute the leading eigenvector. This is indeed a
good strategy, but on sparse Erdős-Rényi random graphs it does not work down to the threshold βc . The problem is
that for an Ising model living on a large sparse Erdős-Rényi graph the tail of the spectra of the commonly considered
matrices is dominated by eigenvectors localized on some small subgraphs (e.g. around high-degree nodes or on some
kind of hanging sub-trees).
The strategy in designing spectral method that does not have the problem with localized eigenvectors (at least
on a much larger class of graphs) is to try to mimic what BP is doing when departing from the paramagnetic fixed
point. Eq. (78) is precisely the linearization of BP around that fixed point and suggests a definition of a so-called
non-backtracking matrix B defined on directed edges as
Bi→j,k→l = Jij δil (1 − δjk ) .

(89)

Denoting M the number of edges in the non-directed graph, this is a 2M × 2M matrix. The term (1 − δjk ) means
that a walk that follows non-zero elements of the matrix is not going back on its steps, this motivates the name of
the matrix. Building on [40, 180], it was proven in [274] that the planted configuration is indeed correlated with the
signs of elements of the eigenvector of B corresponding to the largest (in module) eigenvalue as soon as β ∗ > βc .
Fig. 7 depicts the spectrum on the matrix B (89) for a planted spin glass model on an ER graph of average degree
c and having N = 2000 nodes. The Jij = ±1 were generated following the procedure described in sec. II A with
β ∗ = 0.55. We observe (the proof is found in [274]) that on sparse random graphs
√
• For β ∗ < βc = atanh(1/ c), corresponding to the paramagnetic phase of the planted
√ spin glass model: The
spectrum of B is the same as it would be for β ∗ = 0, bounded by a circle of radius c.

41
√
• For β ∗ > βc = atanh(1/ c), corresponding to the ferromagnetic phase of the planted spin glass model: The
spectrum of B has the same property as in the previous case, except for two eigenvalues on the real axes. One
of them is out of the circle and has value c tanh β ∗ . Its corresponding eigenvector is positively correlated with
the planted assignment.
The performance of the non-backtracking spectral algorithms will be further discussed and compared to other
algorithms for the example of clustering of networks in sec. V D. In sec. V G we will also illustrate how this method
performs on some non-random graphs. Just as all the other spectral algorithms, this method is not robust to adversarial
changes in a small part of the graph (such as planting a small dense sub-graph). This can hence be a disadvantage
with respect to methods that optimize some global cost function that is not affected by small local changes. We,
however, illustrate in Fig. 13 that spectrum of the non-backtracking matrix of many real graphs does not deviate
excessively from its behavior on random graphs.
4

4

3

c=3

3

2

2

1

1

0
-4

-3

-2

-1

c=8

0
0

1

2

3

4 -4

-3

-2

-1

0

-1

-1

-2

-2

-3

-3

-4

-4

1

2

3

4

FIG. 7. Spectrum of the non-backtracking operator
√ for the ±J planted spin glass model. Right: In the undetectable phase
where all eigenvalues are inside a circle of radius c. Left: In the detectable phase where one eigenvalue is outside the circle
on the real axe, its value is c tanh β ∗ . Figure taken from [274].

B.

On computational hardness of the hard phase

In section III E we described the four phases that we observe in the Bayes optimal setting of inference problems.
The most theoretically intriguing is the hard detectable phase that appears between the detectability transition cc ,
and the spinodal transition cs . In this phase both Monte Carlo and belief propagation, that aim to sample the Bayesoptimal posterior distribution, fail and get blocked in the metastable paramagnetic phase. The spinodal transition
cs is also a barrier for all known spectral algorithms. Some of the spectral methods including the one based on the
non-backtracking matrix saturate this threshold.
There are many other problems where a phase of the same physical origin was discovered. Among them the
best known are the LDPC error correcting codes [217] and the planted clique problem [9, 75, 142, 148]. Some
recent works prove hardness in other problems assuming the hardness of planted clique (in the corresponding region)
[31, 47, 118, 197]. There is a lot of current interest in this kind of results.
When a physicist sees a phenomenon as the one above she immediately thinks about universality and tends to
conjecture that the spinodal transition in inference problems will be a barrier for all polynomial algorithms. There is,
however, the case of planted XOR-SAT that presents the same phenomenology as described in section III E. Moreover,
the hard phase in planted XOR-SAT extends to infinity (the paramagnetic state is always locally stable). The work
of [89] actually suggests that the right scaling for the number of clauses at which planted K-XOR-SAT becomes
easy is when the number of clauses is at least O(N K/2 ). Yet, thanks to its linearity planted K-XOR-SAT can be
solved in polynomial time by Gaussian elimination. Hence the above conjecture about the hard phase is false. But
the physicist’s intuition is rarely entirely wrong, and in those rare cases when it is indeed wrong it is for absolutely
fundamental reasons. In our opinion, the question:
In which exact sense is the hard phase hard?

42
is the most fundamental question related to the interface between statistical physics, statistical inference and average
computational complexity.
A line of research that should be explored, and seems quite accesible, is proving or disproving that there is no way
to use belief propagation or Monte-Carlo sampling with a mismatched prior or model in order to improve the position
of the spinodal transition.
Let us also note here that the spinodals of first order phase transitions in Bayes-optimal inference problems create
a much cleaner algorithmic barrier than phase transitions described in optimization or in constraint satisfaction
problems such as K-satisfiability. In K-SAT the idea of hard instances being close to the satisfiability threshold goes
back to the works by [62, 219, 280]. However, even after more than 20 years of study it is still not clear what is the
limiting constraint density up to which polynomial algorithms are able to find satisfiable solutions in random K-SAT
(even at very large K). For a related discussion see [321].
C.

Spatial coupling

Spatial coupling is, from a physics point of view, a very natural strategy to beat the hard-easy spinodal transition
in systems where the structure of the graphical model can be designed. It is closely related to nucleation in physics.
An exponentially long living metastability is only possible in mean-field systems. Let us describe here why it is not
possible in systems embedded in finite Euclidean dimension d. Metastability means that there are two competing
phases. Let us call ∆f the free energy difference between them. Let us consider a large system in the metastable
state with a droplet of radius L (small compared to the size of the whole system) that by random fluctuation appears
to be in the stable state. Creating such a droplet will cost some energy due to its interface where the two phases do
not match. This is quantified by the surface tension Γ. The total energy balance can be written as
E(L) = ΓLd−1 − ∆f Ld .

(90)

Physical dynamics can be roughly thought of as gradient descent in this energy. Notice that E(L) is increasing for
L < L∗ = Γ(d − 1)/(d∆f ) and decreasing for L > L∗ . This means that if fluctuations created a droplet larger than
the critical size L∗ it will continue growing until it invades the whole system. The above argument works only in
geometries where the surface of a large droplet is much smaller than its volume. On random sparse graphs or fully
connected geometries it does not work. This is why nucleation does not exist in mean-field systems. Notice also that
nucleation is simpler for lower dimensions d, the lowest meaningful case being d = 1.
We can conclude that in order to avoid long living metastability we need to work with finite-dimensional systems.
But the situation is not so simple either. In systems where the mean-field geometries are used, they are used for
a good reason. In LDPC codes they assure large distance from spurious codewords, in compressed sensing random
projections are known to conserve well the available information, etc. To avoid metastability we need to mix the meanfield geometry on a local scale with the finite dimensional geometry on the global scale. This is exactly what spatial
coupling is doing. Originally developed for LDPC error correcting codes [188, 192] and closely related to convolutional
LDPC that have longer history [90], spatial coupling has been applied to numerous other settings. Again we will
discuss some of them in more details in section VI E. To introduce the concept here, we restrict ourselves to the
spatially coupled Curie-Weiss model as introduced in [121] and studied further in [49].
The Curie-Weiss (C-W) model in a magnetic field is arguably among the simplest models that present a first order
phase transition. The C-W model is a system of N Ising spins Si ∈ {−1, +1}, that are interacting according to the
fully connected Hamiltonian
HC−W (S) = −

N
X
J X
Si Sj − h
Si ,
N
i=1

(91)

hi,ji

where the notation h·, ·i is used to denote all unique pairs, J > 0 is the ferromagnetic interaction strength and h ∈ R
is the external magnetic field.
Let us briefly remind what is obvious to everybody trained in statistical physics. For positive magnetic field h > 0
the equilibrium magnetization m(J, h) > 0 is positive, and vice versa. There exists a critical value of the interaction
strength Jc = 1 such that: for J < Jc the limh→0+ m(J, h) = limh→0− m(J, h) = 0, and for J > Jc we have
limh→0+ m(J, h) > 0 > limh→0− m(J, h). The latter is a first order phase transition, in the “low temperature” regime
J > Jc the system keeps a non-zero magnetization even at zero magnetic field h. There exists a spinodal value of the
magnetic field
!
r
p
J −1
,
(92)
hs (J) = J(J − 1) − atanh
J

43
Chain of
C-W
systems
(w = 2)

1

=80
=160
=240
=320
=400
=480
=560
=640

0.5

z
m

+1
−1

+1

−1

+1

Single C-W
system
(n = 8)

-0.5

-1

−1

+1

0

+1

0

50

100

150

200
z

250

300

350

400

FIG. 8. Left: A schematic graphical representation of the spatially coupled C-W model. A chain of C-W models interacting
within a certain range w (w = 2 in the figure). In the zoomed part the fully-connected structure of the single C-W model is
shown. A connection between C-W models along the chain indicates connections between all spins contained in both models.
Right: The propagating wave obtained by Monte-Carlo heat-bath simulation for the following parameters: J = 1.4, L = 400,
n = 100, h = 0.05, w = 5, wseed = 10 and hseed = 0.3. Simulations were performed using random sequential update and
therefore time must be rescaled by τ = t/N . The figure is taken from [49].

with the following properties: If the magnetizations are initialized to negative values and the magnetic field is of
strength 0 < h < hs (J), then both local physical dynamics and local inference algorithms, such as the Monte-Carlo
sampling, will stay at negative magnetization m− (J, h) < 0 for times exponentially large in the size of the system.
The spinodal value of the magnetic field hs (J) acts as an algorithmic barrier to equilibration, when initialized at
Si = −1 for all i = 1, . . . , N , and hence to successful inference. For h > hs (J) it is, on the other hand, easy to reach
the equilibrium magnetization m+ (J, h).

D.

The spatially coupled Curie-Weiss model

We introduce here the construction from [121] that might not appear very transparent at first. But its idea is really
simple. We create a one-dimensional chain of C-W models, put a nucleation seed somewhere and let the sub-systems
interact in order to ensure that the nucleation proceeds to the whole chain.
As illustrated in Fig. 8, we consider a one-dimensional chain of (2L + 1) C-W systems, where each of the C-W
system has n spins (referred to as a “block”) and is labelled by the index z ∈ {−L, . . . , L}. The result is that a
configuration S of the full system is now given by the values of N = n(2L + 1) spins, each labelled by a compound
index:
S = {Siz ∈ {+1, −1} : i ∈ {1, . . . , n} , z ∈ {−L, . . . , L}} .

(93)

In a similar way, the uniform external magnetic field h for a single system is replaced by an external field profile hz .
As far as the coupling is concerned, every spin not only connects to all spins in the same location z but also all spins
within w blocks from z. The corresponding Hamiltonian is then
L
n
X
X
1 X
0
0
Hn,L (S) = −
Jzz Siz Sjz −
hz
Siz .
n
0
i=1
hiz,jz i

(94)

z=−L

The couplings between spins are Jzz0 = Jg (|z − z 0 |/w) /w, where the function
 g satisfies the following condition
P+∞
1
−1
g (|x|) = 0, ∀ |x| > 1 and we choose its normalization to be w z=−∞ g w |z| = 1.
In order to ensure that the system nucleates, we must increase the field hz at some point on the chain. We choose
the magnetic field profile in such a way that hz = hseed > hsp in some small region given by z ∈ {0, . . . , wseed }.
Everywhere else, hz = h  hsp , such that the average field strength is still small. To illustrate the nucleation in Fig. 8
we initialize to Siz = −1 everywhere and let the system evolve under Monte-Carlo dynamics. We indeed observe a
nucleation “wave” that invades the whole system in a time proportional to the number of blocks. The system can be
designed in a such a way that the spinodal regime entirely disappears [121].

44
E.

A brief history of time indices in the TAP equations

In physics of disordered system, the historically first form of iterative procedure that was used to compute marginals
(local magnetizations) in a densely connected spin glass model are the Thouless-Anderson-Palmer (TAP) equations
[300]. In the early statistical physics works on error correcting codes TAP equations were used for decoding. The
primary goal of the paper by Kabashima and Saad [154], who pointed out the connection between belief propagation
and the Bethe approximation, was to compare the performance of TAP and BP. In dense systems the fixed points of
both BP and TAP agree. However, it was reported e.g. in [150, 151] that TAP equations did not converge in some
cases even when BP does and there is no replica symmetry breaking. With the renewed interest in message passing
of dense systems that followed the work of Donoho, Maleki and Montanari on compressed sensing [81], it is definitely
worth making some up-to-date remarks on the TAP equations and their convergence.
The fixed point of the TAP equations


N
N
X
X
2
mi = tanh βh + β
Jij mj − β 2 mi
Jij
(1 − m2j )
(95)
j=1

j=1

describes the replica symmetric marginals of the following fully connected spin glass Hamiltonian
X
X
HSK (S) = −
Jij Si Sj − h
Si ,

(96)

i

hi,ji

˜
where Si are the Ising spins,
√ and Jij are independent random variables of zero mean and variance J/N . Hence a
typical value of Jij = O(1/ N ).
Fixed points of the TAP equations are stationary points of the corresponding TAP free energy [300]. Derived as
such, TAP equations do not include any time-indices and were mostly iterated in the simplest form of having time
t − 1 on all magnetizations on the right hand side, and t on the left hand side. At the same time, as we said above,
a number of authors reported that the TAP equations have problems to converge even when BP does converge and
replica symmetry is not broken [150, 151]. No satisfactory explanation for this existed in the literature until very
recently [39]. The reason behind this non-convergence was that the time indices were wrong and when done correctly,
convergence is restored, and on sufficiently dense systems TAP equations behave as nicely as BP.
For the Ising spin glass this was first remarked by Bolthausen [39], who proved the convergence of TAP with the
right time indices in a non-zero magnetic field (and out of the glassy spin glassy phase). His paper inspired the proof
technique of [26] but otherwise passed quite unnoticed. In the context of linear estimation, as discussed in Sec. VI,
the analog of the TAP equations is the approximate message passing algorithm (AMP) [82] (but see also [155]). AMP
is often thought of as the large degree limit of belief propagation and in that way of thinking the correct time indices
come up naturally. Here we re-derive TAP for the Ising spin glass as a large degree limit of BP to clarify this issue
of time indices.
P
Let us go back to the belief propagation for sparse graphs eq. (72). Define messages mti→j ≡ tanh (h + β k6=j uk→i
t−1 ).
√
Recall that Jij , and hence also u, is of order O(1/ N ) and expand to the leading order. BP then becomes


X
t−1 
mti→j = tanh βh + β
Jki mk→i
.
(97)
k6=j

In this dense-graph message passing there is not yet any surprise concerning time indices. The TAP equations are
closed on full magnetizations (not cavity ones) defined as
!
X
t−1
t
mi = tanh βh + β
Jki mk→i .
(98)
k

In order to write the TAP equations we need to write messages in terms of the marginals up to all orders that
contribute to the final equations. We get by Taylor expansion with respect to the third term
!
mti→j = tanh βh + β

X
k

t−1
Jki mt−1
k→i − βJij mj→i

≈ mti − βJij mjt−1 [1 − (mti )2 ] ,

(99)

45
which leads to the iterative TAP equations




X
X
2
mti = tanh βh + β
Jij mt−1
− β 2 mt−2
Jij
[1 − (mjt−1 )2 ] .
j
i


j

(100)

j

Observe the time index (t − 2) before the last sum, which is not intuitive on the first sight. These equations are the
ones written and analyzed by Bolthausen [39].
In his proof, Bolthausen showed that the behavior of these iterative equations can be followed using what can
be called the state evolution of the TAP algorithm, at least outside of the spin glass phase (the algorithm does not
converge inside the spin glass phase). It works as follows: we consider the interactions Jij to have zero mean and
˜ . Going back to eq. (98), the argument of the tanh is a sum of independent variables (by the assumptions
variance J/N
of belief propagation) and thus follows a Gaussian distribution. Denoting m and q the mean and the variance of this
distribution of messages, two self-consistent equations can be derived
Z


p
2
1
e−z /2 tanh βz J˜ q t−1 + βh dz
(101)
mt = √
2π
Z


p
2
1
e−z /2 tanh2 βz J˜ q t−1 + βh dz
qt = √
(102)
2π
The fixed point of these equations provides the well known replica symmetric solution for the magnetization and the
spin overlap in the Sherrington-Kirkpatrick model, as was first derived in [285]. Since this is the correct solution
outside the spin glass phase (above the so-called dAT line [69], it is reassuring to see that we have a convergent
algorithm that will find the marginals fast in the densely connected case. We will return to this observation in section
VI C when we will discuss inference in densely connected models.
Note also that eqs. (101-102) correspond to a parallel update where all magnetizations are updated at the same
time. Very often, it has been found that such update might be problematic for convergence and that one should
update spins one at a time for better convergence properties. This, will also be discussed in VI C when discussing the
AMP algorithm [201]. In that case one chooses randomly at each time step which spin should be updated, then one
can still use eq. (100) with a careful book-keeping of the time-indices. The state evolution eqs. (102) become in that
case differential equations in time. Instead of q t+1 = f (q t ) as in eqs. (102) we have now dq/dt = f (q t ) − q t in the
continuous time limit. We show an exemple of the behavior of these two approaches in Fig. 9. While in this case the
parallel update was faster, we will see that sometimes the other one could be beneficial because it is smoother and
avoids e.g. the overshoot we see in the first iteration of the parallel update.
TAP with parallel update

TAP with random sequential update

0.8

0.7

0.7

0.6

0.6

0.5

0.5

q(t)

q(t)

0.8

0.4

0.3

0.4

0.3

Algorithm

0.2

Algorithm

0.2

State evolution

State evolution

0.1

0

0.1

0

5

10

t

15

0

0

5

10

15

t

FIG. 9. State evolution versus actual simulation for the TAP equations in the Sherrington-Kirkpatrick with parallel update
(left) and random sequential one (right), for a spin glass system with N = 5000, 1/β = T = 3, J˜ = 1 and h = 4. The variance
q of the distribution of messages is shown. In both cases, the state evolution equations (with the iterative form on the left and
the differential one on the right) perfectly describe the evolution of the algorithm.

These considerations stems from the derivation of TAP equations starting from the BP equations for random
interaction. To finish this section, we point out that there have been a large body of work in trying to write TAP

46
equations for interaction that are not random. We refer the interested reader to the book collection [244] for detailed
discussion of the TAP equations in this context. Along this line of thoughts, the adaptive TAP approach [245] allows
to learn the structure of the matrix J, while the structure of the Onsager term can be exactly computed depending
on the structure of the correlations (see e.g. [245, 247, 282]).
Since these approaches do not lead to a natural time ordering, it is worth mentioning that, very recently, these have
been revisited [243] in the light of the “time indices” subtlety, leading a much more satisfying theory of the dynamics
of the TAP equations on general matrices.

47
V.

CLUSTERING OF NETWORKS AND COMMUNITY DETECTION

Clustering is one of the most commonly studied problems in machine learning. The goal is to assign data points to
q clusters in a way that the cluster is somehow representative of the properties of the points inside it. About 15 years
ago, a large number of scientists got interested in the research on so-called complex networks, i.e. sets of nodes where
some pairs are connected by edges, arising in a wide range of settings. A natural question followed: How to cluster
networks? That is, how to assign nodes of a graph in groups such that the largest possible part of the edges can be
explain by a set of affinities between the clusters (instead of between the nodes themselves)? Typically, one considers
a number of clusters much smaller than the number of nodes. Part of the motivation for studying this question came
from studies of social networks which resulted in the name for this problem: community detection.
A.

Biased overview of community detection

The literature on community detection is vast and there are a number of reviews that make a great effort in putting
the existing works in chronological and merits perspective, e.g. [94]. Here we give a perspective biased towards the
approach of statistical physics of disordered systems.
In community detection we are given a graph G(V, E) of N nodes and M edges. We search to assign the nodes
into q groups, calling the assignment si ∈ {1, . . . , q}, i = 1, . . . , N . In the most commonly considered assortative
setting, the goal is that there are many edges between nodes of the same group and few edges among the groups.
Traditionally, in statistics or machine learning, this problem would be solved by a spectral relaxation. One takes
a matrix associated to the graph, its adjacency, its Laplacian or even something else, and compute its largest (or
smallest) eigenvalues and associated eigenvectors. The elements of these eigenvectors would then reveal in a visible
way the community structure [77, 91], for more details see the nice review [307].
An influential set of works on community detection is due to Mark Newman and collaborators who introduced the
so-called modularity function [63, 108, 237]


1 X
di dj
Q=
Aij −
δsi ,sj ,
(103)
2M
2M
hi,ji

where the sum is over all distinct pairs of nodes, di is the degree of node i. Community detection is performed by
finding an assignment of nodes that maximizes this modularity. Optimizing the modularity is a NP-hard problem. In
practice, however, many efficient heuristic algorithms were suggested for this task, see e.g. [94].
Modularity-based community detection has nice connections to statistical mechanics and spin glasses. Reichardt
and Bornholdt [263] pointed out that modularity maximization can be seen as a case of finding a ground state of a
particular Potts spin glass. In order to assess the statistical significance of community structure found by maximizing
the modularity, ref. [263] estimated the maximum value of modularity on a random graph, as the null model. Indeed,
on sparse random graphs the maximum value of modularity can be very large even when no community structure
whatsoever is present. To give an example, results on random graph partitioning [322] teach us that a random 3regular graph can be divided into two groups of roughly equal size in such a way that only about 11% of edges are
between the groups. This visually looks like a very good community structure with high modularity. Yet there were
no communities in the process that created the random graph. The need to compare the value of modularity to a null
model is not very satisfying. Ideally, one would like to have a method able to tell that the underlying graph does not
contain any communities if it did not.
PLet us assume that we believe there were q communities respectively taking a fraction na of nodes, a = 1, . . . , q,
a na = 1, and the probability that a given random node from community a is connected to a given random node
from community b is cab /N , where N is the notal number of nodes and cab is a system size independent q × q matrix.
There will be roughly cab na nb N edges between groups a and b, and caa n2a N/2 edges inside group a. The above very
natural assumption is called the sparse stochastic block model (SBM) [134, 314], and is commonly considered in this
setting. As Bayesian probability theory teaches us, the best possible way of inferring the hidden communities is to
enumerate all assignments of nodes into communities that have roughly the above sizes and respective number of
edges between groups. Subsequently one averages over all such assignments (factoring away the global permutational
symmetry between groups). This way one ends up with probabilities for a given node to belong to a given group. For
a truly random graph these probabilities will not be different from the fractions na .
From a computationally complexity point of view, exactly enumerating assignments as above is even harder than
finding optimal values of a cost function such as the modularity. However, from a statistical physics point of view
evaluating the above averages is equivalent to the computation of local magnetizations in the corresponding Potts
glass model. Notice, for instance, that planted coloring [182] is simply a special case of the stochastic block model

48
(caa = 0, and cab = ccd for all a 6= b, c 6= d). We shall summarize in the rest of this section the asymptotically exact
analysis of the stochastic block model [71, 72].

B.

Asymptotically exact analysis of the SBM

We explained all the methodology needed for the asymptotic analysis of the SBM on the example of the planted
spin glass in section II. In this section, we consider again the teacher-student scenario where the teacher generates
the observed graph from the stochastic block model. We denote θ = {na , cab } the set of parameters of the SBM. We
remind that na is the fraction of nodes taken by group a, a = 1, . . . , q, and cab /N is the probability that a given node
in group a connects to a given node in group b. Again we denote by θ∗ the actual values that the teacher used to
generate the graph G. The graph is represented by its adjacency matrix Aij , with N nodes and M edges. Further
si ∈ {1, . . . , q} are the Potts spins denoting the assignments to the respective groups.
The posterior likelihood of the SBM reads
P (s|G, θ) =

N
Y
Y
A 
cs s 1−Aij
1
nsi
csi sj ij 1 − i j
Z(G, θ)
N

(104)

hi,ji

i=1

where the second product is over all distinct pairs, and we included a constant N −M into the partition function. We
need to evaluate the maximum-mean-overlap (MMO) estimator (22) in order to find a configuration that maximizes
the number of nodes that are assigned into the same group as the planted configuration. In a physics language, we
would rather think of the corresponding Potts model Hamiltonian
H(G, θ, s) =

N
X
i=1

log nsi +

Xh
hi,ji


cs s i
Aij log csi sj + (1 − Aij ) log 1 − i j
N

(105)

and the MMO estimator turns out to be computed from the local magnetizations of this Potts model.
When the parameters θ∗ are not known, we need to learn them according to the lines of section II D. We write
the posterior distribution on the parameters P (θ|G) and notice that it is concentrated around the maximum of the
partition function Z(G, θ). We then iteratively maximize this partition function by imposing the Nishimori conditions
N
1 X
E(δsi ,a ) ,
N i=1
X
1
=
E(δa,si δb,sj )
N na nb

na =
cab

(ij)∈E

caa =

X
2
E(δa,si δa,sj ) .
2
N na

(106)
when a 6= b ,

(107)
(108)

(ij)∈E

This learning algorithm was called expectation maximization learning in [71]. The expectations were computed using
belief propagation or Gibbs sampling as will be discussed subsequently. But even if the expectations were computed
exactly, this parameter learning algorithm depends strongly on the initial choice for the parameters θ. It is easily
blocked in a local maximum of the partition function. The size of the space of initial conditions that need to be
explored to find the global maximum is independent of N , but in practice this is a problem, and completely general
models cab are not so easy to learn. A considerable improvement can be obtained when the learning is initialized
based on the result of spectral clustering. This procedure was suggested in [326] and works nicely, in particular in
conjunction with the best known spectral clustering techniques [180] as will be discussed in sec. V D
Two algorithms that are (in the Bayes-optimal setting) asymptotically exact for computing the marginals and the
expectations in (106-108) were proposed in [71, 72]. One is the Gibbs sampling (or Markov Chain Monte Carlo),
another one is belief propagation (BP). We present here the belief propagation algorithm, as we have introduced it
in sec. III. It can also be used very directly to understand the behavior of the system and the phase transitions in
the thermodynamic limit. The BP algorithm for the SBM was first written by [124] but in that work the parameters
were fixed ad-hoc, not in the Bayes optimal way, nor by learning. Moreover, the algorithm was tested on a network
of only 128 nodes, which is too small to observe any phase transition.
Repeating the steps of sec. III, we now derive BP. For the posterior distribution (104) BP is written in terms
of the so-called messages ψsi→j
that are defined as probabilities that node i is assigned to group si conditioned on
i

49
the absence of edge (ij). Assuming asymptotic conditional independence of the probabilities ψsk→i
for k 6= j, i, the
k
iterative equations read
"
#

Y X
csi sk 1−Aik k→i
1
Aik
i→j
csi sk 1 −
ψsk
,
(109)
ψsi = i→j nsi
Z
N
s
k6=i,j

k

Pq
where Z i→j is a normalization constant ensuring that s=1 ψsi→j = 1. We further notice that because of the factor
1/N , messages on non-edges (ij) ∈
/ E depend only weakly on the target node and hence at the leading order the BP
equations simplify to
"
#
Y X
1
−hsi
i→j
k→i
ψsi = i→j nsi e
csk si ψsk
,
(110)
Z
s
k∈∂i\j

k

where we denote by ∂i the set of neighbors of i, and define an auxiliary external field as
hsi =

N
1 XX
csk si ψskk ,
N
s
k=1

(111)

k

where ψsi i is the BP estimate of the marginal probability of node i


Y X
1

,
ψsi i = i nsi e−hsi
csj si ψsj→i
j
Z
s
j∈∂i

(112)

j

with Z i being again a normalization. In order to find a fixed point of Eq. (110) we update the messages ψ i→j ,
recompute ψ j , update the field hsi by adding the new contribution, subtracting the old one, and repeat.
Once a fixed point of the BP equations is reached, its associated free energy (log-likelihood) fBP = − log Z(G, θ)/N
reads
1 X
c
1 X
log Z i +
log Z ij − ,
(113)
fBP (G, θ) = −
N i
N
2
(i,j)∈E

where
Z ij =

X

cab (ψai→j ψbj→i + ψbi→j ψaj→i ) +

a

a<b

Zi =

X
si

X

nsi e−hsi

YX

csj si ψsk→i
.
j

caa ψai→j ψaj→i

for

(i, j) ∈ E

(114)
(115)

j∈∂i sj

In [71, 72] the author argued that the above equations can be used to analyze the asymptotic performance of
Bayes optimal inference in the SBM, θ = θ∗ . The analysis goes as follow: Let us iterate BP (110) from two different
initializations:
• P
Random initialization, for each (ij) ∈ E we chose ψsi→j = na +si→j where i→j
is a small zero-sum perturbation,
s
i→j
= 0.
s s
• Planted initialization, each message is initialized as ψsi→j = δs,s∗i , where s∗i is the actual assignment of node to
groups (that is known only for the purpose of analysis).
The free energy (113) is computed for the two corresponding fixed points. The fixed point with smaller free energy
(larger likelihood) corresponds asymptotically to the performance of the Bayes-optimal inference. Here “asymptotically” is meant in the sense that the difference between the Bayes-optimal and the BP marginals goes to zero as
N → ∞. This is an important and nontrivial property, because in general the analysis of Bayes optimal inference in
such a high-dimensional setting is very difficult (sharp-P hard). Note that the asymptotic exactness of BP is so far
only a conjecture. Rigorous results showing the asymptotic exactness of BP are so far only partial, either at relatively
high average degrees [233], on in related problems for very low average degrees [224], or e.g. in ferromagnets [73].
Another property of belief propagation that makes it very convenient for the analysis of phase transitions and
asymptotic properties is the fact that BP, in a sense, ignores some part of the finite size effects. Let us compare BP

50
to Markov chain Monte Carlo methods. In the large size limit the performance and behaviour of BP and MCMC are
perfectly comparable, as demonstrated for the stochastic block model in [71]. But for finite sizes there are several
remarkable differences. Consider a set of parameters where the two initializations of BP do not give the same fixed
point. Even if we iterate BP for infinite time we will have two different fixed points. On the other hand, in Monte Carlo
Gibbs sampling for a system of size N , we will eventually find the lower free energy state (just as for every ergodic
Markov chain satisfying detailed balance). The time needed is in general exponentially large in N , but decreases
as the free energy barrier gets smaller and therefore locating the spinodal threshold precisely from finite system size
simulations is computationally more demanding compared to BP. Another manifestation of this property is that in the
paramagnetic phase on a graph with N nodes and M edges the free energy given by BP is always the same number.
Whereas the exact (finite N ) free energy typically has some fluctuations.
In the SBM, belief propagation provides an asymptotically exact analysis of Bayes-optimal inference. In the case
where the model or its parameters θ are mismatching the ones that generated the graph, this is not true in general.
The most striking example here is the MAP estimator, in other words the ground state of the corresponding Potts
model. Modularity maximization is a special case of this setting. The ground state could also be estimated by BP if we
introduced a temperature-like parameter 1/β in the posterior and tuned this temperature to zero. A related algorithm
was discussed recently in the context of modularity-based community detection in [327]. However, as we decrease this
temperature we will typically encounter a point βSG > 1 above which the replica symmetric approximation fails and
the corresponding belief propagation marginals are not any longer asymptotically exact. When this phase transition
at βSG is continuous this will manifest itself as non-convergence of the BP algorithm. Note that depending on the set
of parameters, this spin glass phase can be correlated to the orignal assignment (mixed phase) or not (pure spin glass
phase). This behavior is qualitatively the same as the one seen in the phase digram presented in Fig. 2. In the above
sense it is hence easier to compute the marginals at the Bayes-optimal value of the parameters (that may need to be
learned as described above) than to compute the ground state (or maximize the modularity) that typically is glassy.

C.

Examples of phase transitions in the SBM

The asymptotic analysis of the stochastic block model as described above leads to the discovery of phase transitions.
The most striking picture arises in the case when the parameters of the model are such that the average degree is the
same in every group. In a sense this is the hardest case for inference, because if groups differ in their average degree,
then a simple histogram of degrees gives some information about the correct assignment into groups.
Mathematically, the same average degree condition is expressed as
X
b

cab nb = c

∀b = 1, . . . , q ,

(116)

where c is the average degree. It is straightforward to see that when (116) holds, then
ψsi→j = ns

(117)

is always a (so-called uniform) fixed point (in general not the only one) of belief propagation (110).
Let us now return to the concept of quiet planting explained for the planted spin glass in sec. II E 3 and note
that all we needed for quiet planting to work was precisely the existence of such a uniform fixed point. It is hence
straightforward that the phase transitions described in sections II C and III E also exist in Bayes-optimal inference
in the stochastic block model. Let us give the example of assortative community structure with q = 4 equally sized
groups, na = 1/q for all a = 1, . . . , q, average degree c = 16, and affinity parameters such that caa = cin and
cab = cout for all a 6= b. We define a noise-like parameter  = cout /cin . When  = 0 the model generates four entirely
disconnected groups, and when  = 1 the generation process corresponds to an Erdős-Rényi graph with no information
on the community structure. Without having any experience with phase transition and quiet planting, one might
anticipate that whenever  < 1 the Bayes optimal inference will give a better overlap than a random assignment of
nodes into groups. This expectation is wrong and the phase transition that happens instead is illustrated in Fig. 10.
For other values of parameters we can observe the triplet of discontinuous phase transition discussed in sec. III E.
Actually, the planted coloring example, that was presented in Fig. 6, is a case of the stochastic block model with
q = 5 equally sized groups and extremely dis-assortative structure cin = 0.
For the general case of the sparse SBM, the order of the phase transition, and the locations of the dynamical cd
and detectability cc phase transitions, have to be analyzed case by case mostly numerically. However, the spinodal
transition has a neat analytical expression in the case of the same average degree in every group (116). This expression
is obtained as in sec. III C by analyzing the local stability of the uniform fixed point. Belief propagation (110) linearized

51
700

N=100k, BP
N=70k, MC

0.8
0.6

q=4, c=16

0.4
0.2

detectable

N=10k
N=100k

600
500

convergence time

normalized overlap

1

undetectable

400

q=4, c=16

300
200
100

0
0

0.2

0.4
0.6
ε= cout/cin

0.8

1

εc

0
0

0.2

0.4
0.6
ε= cout/cin

0.8

1

FIG. 10. Left: Normalized overlap between the MMO estimator and the original assignment to groups. We compare the
BP result to Monte Carlo Gibbs sampling. Both agree well, but we see the absence of finite size corrections in BP in the
undetectable phase. Right: The BP convergence time, illustrating the critical slowing down around the phase transition.
Figures are both taken from [71].

around this fixed point reads
i→j
=
a

q
X
X

Tab k→i
,
b

(118)


−1 .

(119)

b=1 k∈∂i\j

where the q × q matrix Tab is defined as
T ab = na

c

ab

c

The phase transition where the paramagnetic fixed point ceases to be locally stable happens when
cλ2T = 1 ,

(120)

where λT is the largest eigenvalue of the matrix T .
In case of equally sized groups in the assortative-diagonal case where caa = cin and cab = cout for a 6= b, the largest
eigenvalue of the T ab matrix (119) is λT = (cin − cout )/(cq) from which (120) implies for the spinodal phase transition
[72]
√
(121)
|cin − cout | = q c .
When the difference between cin and cout is larger, BP detects communities and gives asymptotically the Bayesoptimal overlap. When the difference is smaller BP fails. Depending on the parameters and the order of the phase
transition the detectability phase transition itself can either be given by eq. (120) or happen at strictly smaller
|cin − cout |, giving rise to the hard phase discussed in sec. IV B.
For the planted graph coloring (i.e. the extremely dis-assortative SBM, cin = 0), the detectability phase transition
is of second order for q ≤ 3 and of first order for q ≥ 4. The spinodal transition happens at cs = (q − 1)2 . In the
large q limit the dynamical and detectability phase transitions can be deduced from [320] while keeping in mind the
concept of quiet planting
cd = q(log q + log log q) + O(1) ,
cc = 2q log q − log q − 2 log 2 + o(1) .

(122)
(123)

For the assortative case cin > cout the phase transition is of second order for q ≤ 4 and of first order (although a
very weak one for small q) for q ≥ 5. The limit of dense graphs (average probability for an edge to exist is p, pin/out
are probabilities of connection within/between groups), and large number of groups q → ∞ was analyzed in [193] and
is related to critical temperatures in the dense Potts glass [50, 113]. The regime in which both the dynamical and
detectability phase transitions appear is
N (pin − pout )2 = γ(q)p(1 − p) ,

(124)

52
with the dynamical transitions at γd (q) = 2q log q[1 + o(1)], the detectability transition at γc (q) = 4q log q[1 + o(1)],
and the spinodal threshold at γs (q) = q 2 . In this limit we see that the hard phase between γs (q) and γc (q) is very
wide (compared to the overall scale).
The prediction and location of these phase transitions in the stochastic block model from [71, 72] was followed by a
remarkable mathematical development where the location of the threshold was made rigorous for q = 2 [205, 231, 232].
The detectable side was proven for generic q in [40].
To briefly comment on the case where the average degree of every group is not the same: BP is still conjectured to
be asymptotically exact, but the phase transitions get smeared. As recently studied in [328], the second order phase
transitions disappear, and the first order phase transitions shrink into a multi-critical point beyond which they also
disappear. The situation is similar to the semi-supervised clustering, where the assignment of a fraction of nodes is
known [329]. From a physics point of view it is relevant to note that this problem is an inference interpretation of the
particle pinning that was recently used quite extensively for studies of structural glasses [53].
D.

Spectral redemption for clustering sparse networks

Any reader that did no skip section IV A is now clearly anticipating the use of the non-backtracking matrix for
spectral clustering of networks. Indeed, in [180] the following spectral clustering method was proposed
• Construct the non-backtracking matrix
Bi→j,k→l = δil (1 − δjk ) .

(125)

and compute its largest (in module) eigenvalues until they start to have a non-zero imaginary part. Denote k
the number of those large real eigenvalues.
• Consider the k P
corresponding 2M -dimensional eigenvectors ui→j and create a set of N -dimensional vectors by
resuming ui = k∈∂i uk→i . Then view the k vectors ui as N points in k dimensional space and cluster those
points into k groups using some off-the-shell clustering algorithm such as k-means.
Note in particular that k here is the estimator for the number of groups. Up to this point we were assuming that
the number of groups was known. The performance of this non-backtracking based spectral algorithms compared
to BP, and five other spectral algorithms on clustering of sparse graphs, is illustrated in Fig. 11. All these spectral
algorithms look for the largest (or smallest) eigenvalue of the associated matrix and then assign nodes into groups
based on the sign of the corresponding eigenvector. The five matrices that are most commonly used for clustering are
the adjacency matrix A, the Laplacian L = D − A (where D is a diagonal matrix with the node degrees di on the
diagonal), normalized Laplacian L̃ = D−1/2 LD−1/2 , random walk matrix Wij = Aij /di , and the modularity matrix
Qij = Aij − di dj /(2M ). The spectra of all these matrices on sparse graphs is affected by the existence of large
eigenvalues associated to localized eigenvactors. This fact deteriorated their performance with respect to both BP
and the non-backtracking based algorithm.
Ref. [180] analyzed the performance of the above algorithm on networks generated from the stochastic block model
with equal average degree in every group. A large part of results that were not rigorous in [180] were very remarkably
made rigorous in [40].
On a random graph taken from the configurational model, i.e. at random with a given degree distribution P (d), in
the large
√ size limit, the spectrum of the non-backtracking matrix is (except one eigenvalue) included inside a circle of
radius c̃. Here c̃ is the average excess degree of the graph computed as
P 2
d P (d)
c̃ = Pd
− 1.
(126)
d dP (d)
The one eigenvalue that is not inside this circle is real and has value c̃. This random graph spectrum should be
contrasted with the spectrum of the adjacency matrix (and other related ones) which for graphs with constant
average and unbounded maximum degree is not even bounded in the thermodynamic limit.
Now, consider a graph generated by the SBM with
√ equal average degree for every group. Denote λT an eigenvalue
of the matrix T , eq. (119), then as long as λT > 1/ c̃ the non-backtracking matrix will have purely real eigenvalues
at positions cλT . This way we obtain up to q − 1 additional eigenvalues on the real axes out of the circle of radius
√
c̃. Note that in the assortative-diagonal case the matrix T has q − 1 degenerate eigenvalues λT = (cin − cout )/(qc),
and therefore the non-backtracking matrix has a (q − 1)-times degenerate eigenvalue at (cin − cout )/q which merges
into the bulk below the spinodal phase transition. The spectrum of the non-backtracking operator is illustrated in
Fig. 12.

53
1
1

Overlap

0.8

Non−backtracking
Normalized Laplacian
Modularity
Random Walk
Adjacency
Laplacian
BP

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0
2.5

3

3.5

4

4.5

cin−cout

5

5.5

0
2

6

4

6

8

10

12

Average Degree

FIG. 11. The accuracy of spectral algorithms based on different linear operators, and of belief propagation, for two groups of
equal size. On the
√ left, we vary cin − cout while fixing the average degree c = 3; the detectability transition given by occurs
at cin − cout = 2 3 ≈ 3.46. On the right, we set cout /cin = 0.3 and vary c; the detectability transition is at c ≈ 3.45. Each
point is averaged over 20 instances with N = 105 . The spectral algorithm based on the non-backtracking matrix B achieves
an accuracy close to that of BP, and both remain large all the way down to the transition. Standard spectral algorithms based
on the adjacency matrix, modularity matrix, Laplacian, normalized Laplacian, and the random walk matrix all fail well above
the transition, giving a regime where they do no better than chance. Figure is taken from [180].

1.5

1.0

0.5

-1

µc
1

2

3

-0.5
-1.0
-1.5

FIG. 12. The spectrum of the non-backtracking matrix B for a network generated by the block model with N = 4000, cin = 5,
and cout = 1. The leading eigenvalue is at c√= 3, √
the second eigenvalue is close to µc = (cin − cout )/2 = 2, and the bulk of the
spectrum is confined to the disk of radius c = 3. The spectral algorithm that labels vertices according to the sign of B’s
second eigenvector (summed over the incoming edges at each vertex) labels the majority of vertices correctly. Figure is taken
from [180].

Consequently, when the average degree of every group is the same, condition (116), the phase transition of the
spectral algorithm coincides with the spinodal phase transition of the belief propagation. In other words as long the
average degree is larger than given by eq. 120 the eigenvectors corresponding to the eigenvalues on the real axes out
of the circle contain information about he ground-truth communities.
We need to stress that the agreement between the spectral and belief propagation thresholds arises only for the
case of equal average degree in every group. When the condition (116) is not satisfied there is information about the
groups on the degree and belief propagation explores this information in an optimal way. What concerns the phase
transitions, they tend to be smeared (second order phase transitions disappear, and first order ones become weaker)
[328]. The spectral methods, on the other hand, always have a threshold where the informative eigenvalues merge
back into the bulk. At this spectral phase transition, when (116) does not hold, BP does not witness any particular
change in performance.

54
Let us compare properties of belief propagation and the non-backtracking (NB) based spectral clustering. The
positive point for BP is that it gives a slightly better overlap. The positive points for NB-spectral clustering, that are
shared with other spectral methods, e.g. all those tested in Fig. 11, are that it does not need to get any information
about the model parameters nor the number of groups. It works in time linear in the number of groups, as opposed to
quadratic for BP. It is not expected to be affected by the presence of small loops which usually cause the failure of the
Bethe approximation. An inconvenience that is shared by all the considered spectral methods is that the performance
can deteriorate considerably upon small adversarial changes to the graph (such as planting a small dense subgraph).
Approaches based on some global optimization, such as BP, are expected to be robust towards small local changes.
Another inconvenient property of the non-backtracking matrix is that it can be rather large. However, this is no
obstacle at all because the spectrum (different from ±1) of B can be computed from the spectrum of the following
2N × 2N matrix B 0 [180]


0 D−1
0
B =
,
(127)
−1
A
where A is the adjacency matrix and D is the diagonal matrix with Dii = di being the degree of node i.
The matrix B 0 is smaller but still non-symmetric. The power method for computing the largest eigenvalues works
much better for symmetric matrices. In [271], it was suggested to use the so-called Bethe Hessian matrix defined as
BH(r) = (r2 − 1)1 − rA + D ,

(128)

√
with r = c̃ (or more general the squared root of the leading eigenvalue of (127) which is usually not problematic
to compute) as a basis for spectral clustering. In particular, to take eigenvectors of all of its negative eigenvalues.
Ref. [271] showed that this has all the advantages of the non-backtracking matrix while being N × N and symmetric.
The name Bethe Hessian is motivated by the relation between the Hessian of the Bethe free energy and this matrix
[226, 264], which is particularly apparent when generalized to the weighted case.

E.

More on the non-backtracking operator

The non-backtracking operator as defined above is relatively well known in mathematics, often called the Hashimoto
directed edges adjacency matrix [120]. It is used in a number of works as a tool in proofs, see e.g. [8, 24, 98, 159,
252, 287]. However, using it as a basis for a spectral algorithm was first suggested in [180]. The idea of using the
linearization of BP was inspired by the earlier work of [66]. Since then, a considerable number of papers on promising
algorithmic applications of the non-backtracking matrix on sparse networks started appearing, see [204, 229, 325].
Here, we briefly discuss selected contributions in this direction.
Looking at Fig. 12, we observe a very interesting structure in the bulk spectral density of the non-backtracking
operator. Motivated by this, ref. [272] applied the cavity method designed to compute spectral densities of nonHermitian random matrices [268] to study the non-backtracking spectral√density. A phase transition in the spectral
density was observed to happen on the boundary of the circle of radius c̃. Such a phase transition is absent in the
spectral density of traditional matrices such as the adjacency, Laplacian, random walks etc. This is no proof since we
are only concerned with spectral density not with sub-dominantly numerous eigenvalues, but from a physics point of
view a phase transition is a strong indication that the properties of the non-backtracking matrix on a random graph
are fundamentally different from those of traditional matrices. Ref. [272] also observed that the eigenvalues that
are close to the center of the circle are not visible in the asymptotic calculation of the spectral density. Numerical
investigations revealed that eigenvectors corresponding to these eigenvalues are localized. More detailed understanding
of this localization and its physical consequences is an interesting direction of future work.
In [10], the non-backtracking spectral clustering method was generalized to clustering of sparse hypergraphs. These
are graphs where instead of connected pairs of nodes, one connects k-uples of nodes. Interestingly, using the nonbacktracking based spectral method detection can be done also e.g. in planted constraint satisfaction problems without
even knowing the form of the corresponding constraints. Moreover, for the same reasons as in clustering, the method
works down to previously computed spinodal thresholds [324]. It was also observed that the overlap (performance)
of the spectral method is always a continuous function of the noise parameter, even if the Bayes optimal belief
propagation has a strongly discontinuous transition. This further supports our previous suggestion that the spectral
methods are extremely useful to point towards the right direction, but message-passing based methods should then
be used to improve the precision. Note that the hypergraph clustering can also be formulated as tensor factorization
that became recently popular in connection with interdependent networks or networks with different kinds of edges
that are referred to as multiplexes in the literature on complex networks.

55
Another exciting application of the non-backtracking matrix, or the Bethe Hessian, is the matrix completion problem
[273]. In matrix completion one aims to estimate a low-rank matrix from a small random sample of its entries. This
problem became popular in connection to collaborative filtering and related recommendation systems. Of course, the
fewer entries are observed, the harder the problem is. Ref. [273] shows that with the use of the non-backtracking
matrix the number of entries from which the low-rank structure is detectable is smaller than with other known
algorithms. The way the matrix completion problem is approached is by realizing that a Hopfield model with a Hebb
learning rule can be viewed as a rank r structure that is able to recover the r patterns as distinct thermodynamic
phases [136]. One views the observed entries of the low-rank matrix to complete as a sparse Hopfield model and
estimate the patterns from the leading eigenvectors of the corresponding Bethe Hessian. In a second stage, we use
this as a initialization for more standard matrix completion algorithms in the spirit of [162], and achieve overall better
performance.
Another problem where the leading eigenvalue of the non-backtracking matrix is relevant is percolation on networks
[161]. In percolation, we start with a graph G and remove each of its edges with probability 1 − p. Then we look
at the largest component of the remaining graph and ask whether it covers an extensive fraction of all nodes. The
value of pc below which the giant component does not exist anymore with high probability is called the percolation
threshold. Percolation thresholds were computed for many geometries including random graphs. A remarkable result
concerns any sequence of dense graphs for which it was proven that the threshold is given by the inverse of the leading
eigenvalue of the adjacency matrix [38]. This result does not extend to sparse graphs, where a counter-example is
given by random regular graphs in which the leading eigenvalue is the degree whereas the percolation threshold is the
inverse of the degree minus one.
Based on a message-passing algorithm for percolation, [161] argues that on a large class of sparse locally tree-like
graphs the percolation threshold is given by the inverse of the leading eigenvalue of the non-backtracking matrix.
Specifying for what exact class of graphs this result holds is still an open problem. Ref. [161] also proved that the
non-backtracking matrix provides a lower bound on the percolation threshold for a generic graph (with the exception
of trees). An analogous result was obtained independently by [119]. Due to these results, the non-backtracking matrix
was used to study percolation on real networks [229, 255, 267].

F.

The deceptiveness of the variational Bayes method

Both in physics and in computer science, a very common and popular method for approximate Bayesian inference
is the so-called naive mean field method, often called variational Bayesian inference [27, 244]. For clustering of
sparse networks generated by the stochastic block model this method was studied e.g. in [59, 133]. Away from the
detectability threshold the variational Bayes method gives very sensible results. Deep in the undetectable phase (very
small |cin − cout |) the predicted magnetizations are close to zero, reflecting the undetectability. Deep in the detectable
phase the variational Bayesian method is able to estimate the correct group assignment, and is statistically consistent
as proven in [59]. The performance of variational Bayes in the neighborhood of the detectability transition was studied
in [326]. The conclusions of the comparison to belief propagation are quite interesting. It is not surprising that BP
is more precise and converges faster on sparse graphs. What is surprising is that the naive mean-field inference fails
in a very deceptive way: the marginals it predicts are very biased towards a configuration that is not correlated
with the planted (ground truth) configuration. Ref. [326] quantifies this by introducing the so-called illusive overlap
and comparing it to the actual one. This means that it can be rather dangerous to use naive mean-field for testing
statistical significance of the results it provides and this could have profound consequences beyond the problem of
clustering networks. Especially when trying to use the variational Bayes method in the vicinity of the corresponding
detectable threshold, i.e. in regimes where the number of samples is barely sufficient for estimation. Monte Carlo
methods do not have this problem, but are in general slower.

G.

Real networks are not generated by the SBM

Our study of clustering of networks is asymptotically exact for the case when the networks were generated by the
stochastic block model. In view of actual data sets and realistic applications this is never the case.
One important discrepancy between the stochastic block model and real networks is that the latter often have
communities with a broad degree distribution, whereas the degree distribution of every group generated by the SBM
is Poissonian. A variant of the stochastic block model that takes the degree distribution into account is the so-called
degree corrected stochastic block model [160]. An asymptotically exact belief propagation can also easily be written
for this model, see [317]. In Ref. [317] it was also shown that the corresponding Bethe free energy is a great tool for
model selection between the canonical and degree corrected stochastic block models, and that classical model selection

56
30

3

20

Polblogs q=2
Overlap: 0.8640

10

-40

4

2

3

Dolphins q=2
Overlap: 0.7419

1

-20

20

40

60

80

-4

-2

1
2

-10

Karate q=2
Overlap: 1

2

4

6

8

-4

-2

-1

-1

2

4

6

8

10

12

-2
-20

-2

-30

-3

-3
-4

4

4

3

4

Football q=12
Overlap: 0.9163

2

Adjnoun q=2
Overlap: 0.6250

3

2

1
-4

-2

-1
-2
-3
-4

Polbooks q=3
Overlap: 0.7571

2
1

2

4

6

8

10

12

-5

5
-2
-4

10

-4

-2

-1

2

4

6

8

10

12

-2
-3
-4

FIG. 13. Spectrum of the non-backtracking matrix in the complex plane for some real networks commonly used as benchmarks
for community detection, taken from [7, 108, 171, 196, 236, 319]. The radius of the circle is the square root of the largest
eigenvalue, which is a heuristic estimate of the bulk of the spectrum. The overlap is computed using the signs of the second
eigenvector for the networks with two communities, and using k-means for those with three and more communities. The
non-backtracking operator detects communities in all these networks, with an overlap comparable to the performance of other
spectral methods. As in the case of synthetic networks generated by the stochastic block model, the number of real eigenvalues
outside the bulk appears to be a good indicator of the number q of communities. Figure taken from [180].

criteria do not apply to this problem because of the limited information per degree of freedom. Interestingly, BP for
the degree-corrected SBM also leads to the same non-backtracking matrix upon linearization.
There are of course many more discrepancies between real networks and those generated by the (degree corrected)
SBM. However, when observing the spectrum of the non-backtracking operator on some of the commonly used
benchmarks for community detection, we conclude that the properties derived for the SBM basically stay conserved
also in these examples, see Fig. 13. To conclude this section, let us quote the British mathematician and Bayesian
statistician George E. P. Box (1919-2013): “Essentially, all models are wrong, but some are useful”.

57
VI.

LINEAR ESTIMATION AND COMPRESSED SENSING

The inference problem addressed in this section is the following. We observe a M dimensional vector yµ , µ =
1, . . . , M that was created component-wise depending on a linear projection of an unknown N dimensional vector xi ,
i = 1, . . . , N by a known matrix Fµi
!
N
X
yµ = fout
Fµi xi ,
(129)
i=1

where fout is an output function that includes noise. Examples are the additive white Gaussian noise (AWGN) there
the output is given by fout (x) = x + w with w being random Gaussian variables, or the linear threshold output where
fout (x) = sign(x − κ), with κ being a threshold value. The goal in linear estimation is to infer x from the knowledge
PN
of y, F and fout . An alternative representation of the output function fout is to denote zµ = i=1 Fµi xi and talk
about the likelihood of observing a given vector y given z
P (y|z) =

M
Y
µ=1

A.

Pout (yµ |zµ ) .

(130)

Some examples of linear estimation

To put the immense number of important applications of the above setting into perspective, let us give an illustrative
example from medical diagnostics. Consider we have data about M cancer-patients that took a given treatment. Then
yµ is a variable rating the effectiveness of the treatment for a patient µ. Our goal is to predict the outcome of the
treatment for a new patient that did not undergo treatment yet. To do that, we collect medical and personal data
about each of the M patients, their medical history, family situation, but also for instance Facebook feeds, where
they worked, what sport they did in their life, what education they have, how they rate their relationship with their
wife, and all kinds of other information that is on first sight seems totally unrelated to the eventual outcome of
the treatment. This information is collected in the vector F~µ . We aim to find a set of coefficients xi such that the
PN
observations yµ can be well explained by zµ = i=1 Fµi xi for some Pout (yµ |zµ ).
In the way we set the problem, it is reasonable to think that most patient-attributes will have no effect on the
outcome of the treatment, i.e. for many i’s the xi will be close to zero. In a method, called LASSO (which stands
for least absolute shrinkage and selection operator) [123, 301], that is used everyday for similar types of problems (in
banking or insurance industry for instance) the coefficients xi are such that


!2
M
N
N
X
X
X
LASSO :
min 
yµ −
Fµi xi
+λ
|xi | ,
(131)
x

µ=1

i=1

i=1

where λ is a regularization parameter that controls the trade-off between a good fit and the sparsity of the vector
x. The advantage of LASSO is that it is a convex problem that can be solved efficiently and with strong guarantees
using linear programming. LASSO and related linear regression problems are well studied in statistics and machine
learning, instead of reviewing the literature we refer the reader to the excellent textbook [122].
LASSO can be framed in the Bayesian probabilistic setting by considering
Pout (yµ |zµ ) = √

PN
2
i=1 Fµi xi )
2∆

(yµ −
1
e−
2π∆

(132)

and the prior of xi being the Laplace distribution
`1
PX
(xi ) = e−λ|xi | /2 .

(133)

Another sparsity inducing prior is a combination of strict zeros and fraction ρ of non-zeros
PX (xi ) = ρφ(xi ) + (1 − ρ)δ(xi ) ,

(134)

where φ(xi ) is the distribution of non-zero elements. Linear programming does not work for priors such as (134), but
the statistical-physics framework is more general in this aspect.

58
In many applications, the output variables yµ are binary (interpreted in the above example as the treatment
worked or did not), then the most commonly considered output distribution is the logistic function. The conditional
dependence then is
logistic
Pout
(yµ = 1|zµ ) =

1
.
1 + e−zµ

(135)

The example of this linear estimation setting that is the best known in statistical physics is the perceptron [269],
as introduced briefly in sec. I B 2. In the perceptron setting, components of the unknown vector x are interpreted as
synaptic weights (that are in general not sparse). The rows of the matrix F are patterns presented to the perceptron.
The measurements yµ are binary and correspond to classifications of the patterns (think of some patterns being
pictures of cats yµ = +1
dogs yµ = −1). The goal is to learn the synaptic weights
P and the others being pictures of P
xi in such a way that i Fµi xi > 0 when yµ = +1, and i Fµi xi < 0 when yµ = −1. The random version of the
problem has been considerably studied by physicists, see for instance [102, 103, 170, 286, 310].
Another example well known and studied in the statistical physics literature is the code-division multiple access
(CDMA) problem in which N users simultaneously aim to transmit their binary signal xi ∈ {±1} to one base
station. In order
√ to be able to deal with interference one designs a so-called spreading code associated to each user
Fµi ∈ {±1}/ N , µ = 1, . . . , M , that spreads the signal of each user among M samples. The base station then receives
a noisy linear combination of all the signals yµ , eq. (129). The goal is to reconstruct the signal of every user xi .
Interestingly, the two examples above are the first inference problems on which, to our knowledge, the strategy
discussed in the present chapter has been introduced in a statistical physics. The seminal work of Marc Mézard
derived the TAP equations for the perceptron problem as early as in 1989 [210]. It has then been generalized to many
other cases, mainly by Kabashima and co-workers (see [152, 155] again for the perceptron and [150] for randomly
spreading CDMA).
In this chapter, we shall follow these pioneering works and adopt the Bayesian probabilistic framework for the
above generalized linear estimation problem. In general this is algorithmically much more involved than the convex
relaxations of the problem, and consequently less widely studied. The cases where interesting theoretical analysis can
be done are limited to special classes of matrices Fµi (mostly with random iid elements, but see the discussion in
section VI F). Nevertheless the resulting algorithms can (as always) be used as heuristic tools for more general data
Fµi . One problem where the above linear estimation problem truly involves a matrix Fµi with random iid elements
is compressed sensing [78] to which we dedicate the next sections.
B.

Compressed sensing

Compressed sensing, a special case of the above linear estimation setting, triggered a major evolution in signal
acquisition. Most signals of interest are compressible and this fact is widely used to save storage place. But so far
compression is made only once the full signal is measured. In many applications (e.g. medical imaging using MRI or
tomography) it is desirable to reduce the number of measurements as much as possible (to reduce costs, time, radiation
exposition dose etc.). It is hence desirable to measure the signal directly in the compressed format. Compressed sensing
is a concept implementing exactly this idea by designing new measurement protocols. In compressed sensing one aims
to infer a sparse N -dimensional vector x from the smallest possible number M of linear projections y = F x. The
M × N measurement matrix F is to be designed.
Most commonly used techniques in compressed sensing are based on a development that took place in 2006 thanks
to the works of Candes, Tao, and Donoho [54, 55, 78]. They proposed to find the vector x satisfying the constraints
y = F x which has the smallest `1 norm (i.e., the sum of the absolute values of its components). This `1 minimization
problem can be solved efficiently using linear programming techniques. They also suggested using a random matrix
for F . This is a crucial point, as it makes the M measurements random and incoherent. Incoherence expresses the
idea that objects having a sparse representation must be spread out in the domain in which they are acquired, just as
a Dirac or a spike in the time domain is spread out in the frequency domain after a Fourier transform. Rather strong
and generic theoretical guarantees can be derived for the `1 -based sparse reconstruction algorithms [123]. These ideas
have led to a fast, efficient and robust algorithm, and this `1 -reconstruction is now widely used, and has been at the
origin of the burst of interest in compressed sensing over the last decade.
It is possible to compute exactly the performance of `1 -reconstruction in the limit N → ∞, and the analytic
study shows the appearance of a sharp phase transition [83]. For any signal with density ρ = K/N , where K is the
number of non-zero elements, the `1 -reconstruction gives the exact signal with probability one if the measurement rate,
α = M/N , is above the so-called Donoho-Tanner line, α > α`1 (ρ) > ρ, see Fig. 15 (black dashed line). Interestingly,
the very same transition can be found within the replica method [101, 157] but can also be computed, rigorously, with
the approximate message passing technics that we discuss in this chapter [81].

59
If one is willing to forget the (exponentially large) computational cost, sparse signal reconstruction from random
linear projections is possible whenever α > ρ (this can be done for instance by exhaustively checking all possible
positions of the non-zero elements and trying to invert the reduced system of linear equations). A challenging
question is whether there is a compressed sensing design and a computationally tractable algorithm that can achieve
reconstruction down to such low measurement rates. In a series of works [173, 175], complemented by an impressive
rigorous proof of [79], this question was answered positively under the assumption that the empirical distribution
of the non-zero elements of the signal is known. This contribution to the theory of compressed sensing is explained
below and uses many of the statistical physics elements from section II.

C.

Approximate message passing

Approximate message passing (AMP) is the algorithmic procedure behind the aforementioned work of [81] on linear
estimation and compressed sensing. Origins of AMP trace back to statistical physics works on the perceptron problem
[210] and on the code division multiple access (CDMA) [150].
What nowadays is known as the generalized approximate message passing (G-AMP) algorithm, i.e. the generalization to arbitrary element-wise output functions, was derived for the first time, to our knowledge, in connection
to perceptron neural network in [155, 210]. The algorithm presented in [155] had the correct time-indices and good
performance was illustrated on the CDMA problem. However, for an unknown reasons, that paper did not have much
impact and many of its results waited to be rediscovered in the connection to compressed sensing much later. The
opinion that this type of algorithm has invincible convergence problems prevailed in the literature until the authors of
[45] demonstrated impressive performance of a related message passing algorithm for learning in binary perceptron.
AMP is closely related to relaxed-belief propagation (r-BP) [256]. The algorithm of [45] can actually be seen as a
special case of r-BP.
Donoho, Maleki and Montanari introduced AMP for compressed sensing [81], they also established the name. Very
remarkably, the rigorous proof of its performance for random matrices F was done in [25, 26]. A comprehensive
generalization to arbitrary priors and element-wise output function is due to Rangan [257], who established the name
generalized approximate message passing for the version with a general output channel. Many variants of this approach
have then been considered (see for instance, without being exhaustive, [92, 158, 175, 288]).
As is clear from the earlier papers [150], AMP is actually very closely related to the TAP equations [300] presented
in section IV E. The main difference between TAP and AMP is that TAP was originally derived for a model with
Ising spins and pairwise interactions, whereas AMP is for continuous variables and N -body interactions.

FIG. 14. Factor graph of the linear estimation problem corresponding to the posterior probability (136). Circles represent the
unknown variables, whereas squares represent the interactions between variables.

60
1.

Relaxed Belief propagation

A detailed statistical physics derivation of the AMP algorithm was presented in [175] for the Gaussian output
channel, and for general output as a part of the yet more general problem of matrix factorization in [153]. Here we
shall describe the main elements necessary for the reader to understand where the resulting algorithm comes from.
We are given the posterior distribution
P (x|y, F ) =

M
N
Y
1 Y
Pout (yµ |zµ )
PX (xi ) ,
Z µ=1
i=1

where

zµ =

N
X

Fµi xi ,

(136)

i=1

where the matrix Fµi has independent random entries (not necessarily identically distributed) with mean and variance
O(1/N ). This posterior probability distribution corresponds to a graph of interactions yµ between variables (spins)
xi called the graphical model as depicted in Fig. 14.
A starting point in the derivation of AMP is to write the belief propagation algorithm corresponding to this graphical
model, as done originally in [21]. The matrix Fµi plays the role of randomly-quenched disorder, the measurements yµ
are planted disorder. As long as the elements of Fµi are independent and their mean and variance of of order O(1/N )
the corresponding system is a mean field spin glass. In the Bayes-optimal case (i.e. when the prior is matching the
true empirical distribution of the signal) the fixed point of belief propagation with lowest free energy then provides
the asymptotically exact marginals of the above posterior probability distribution.
For model such as (136), where the interaction are not pairwise, the BP equation we presented in sec. III B need
to be slightly modified, for examples see [217], but the logic is very similar. BP still implements a message-passing
scheme between nodes in the graphical model of Fig. 14, ultimately allowing one to compute approximations of the
posterior marginals. Messages mi→µ are sent from the variable-nodes (circles) to the factor-nodes (squared) and
subsequent messages mµ→i are sent from factor-nodes back to variable-nodes that corresponds to algorithm’s current
“beliefs” about the probabilistic distribution of the variables xi . It reads
Y
1
PX (xi )
mγ→i (xi )
(137)
mi→µ (xi ) =
zi→µ
γ6=µ
Z Y
X
1
mµ→i (xi ) =
[dxj mj→µ (xj )] Pout (yµ |
Fµl xl )
(138)
zµ→i
j6=i

l

While being easily written, this BP is not computationally tractable, because every interaction involves N variables,
and the resulting belief propagation equations involve (N − 1)-uple integrals.
Two facts enable the derivation of a tractable BP algorithm: the central limit theorem, on the one hand, and
a projection of the messages to only two of their moments (also used in algorithms such as Gaussian BP [284] or
non-parametric BP [293]). This results in the so-called relaxed-belief-propagation (r-BP) [256]: a form of equations
that is tractable and involves a pair of means and variances for every pair variable-interaction.
Let us go through the steps that allow
P to go from BP to r-BP and first concentrate our attention to eq. (138).
Consider the variable zµ = Fµi xi + j6=i Fµj xj . In the summation in eq. (138), all xj s are independent, and
P
2
thus we expect, according to the central limit theorem, j6=i Fµj
xj to be a Gaussian random variable with mean
P
P
2
ωµ→i = j6=i Fµj aj→µ and variance Vµ→i = j6=i Fµj vj→µ , where we have denoted
Z
Z
ai→µ ≡ dxi xi mi→µ (xi )
vi→µ ≡ dxi x2i mi→µ (xi ) − a2i→µ .
(139)
We thus can replace the multi-dimensional integral in eq. (138) by a scalar Gaussian one over the random variable z:
Z
(z−ωµ→i −Fµi xi )2
−
2Vµ→i
.
(140)
mµ→i (xi ) ∝ dzµ Pout (yµ |zµ )e
We have replaced a complicated multi-dimensional integral involving full probability distributions by a single scalar
one that involves only first two moments.
2 2
One can further simplify the equations: The next step is to rewrite (z − ωµ→i − Fµi xi )2 = (z − ωµ→i )2 + Fµi
xi −
√
2(z − ωµ→i )Fµi xi Fµi xi and to use the fact that Fµi is O(1/ N ) to expand the exponential
Z
(z−ω
)2
− 2V µ→i
µ→i
mµ→i (xi ) ∝ dzµ Pout (yµ |zµ )e
(141)


1
2 2
2 2
× 1 + Fµi
xi − 2(z − ωµ→i )Fµi xi + (z − ωµ→i )2 Fµi
xi + o(1/N ) .
2

61
At this point, it is convenient to introduce the output function gout , defined via the output probability Pout as
R
(z−ω)2
dzPout (y|z) (z − ω) e− 2V
.
(142)
gout (ω, y, V ) ≡
R
(z−ω)2
V dzPout (y|z)e− 2V
The following useful identity holds for the average of (z − ω)2 /V 2 in the above measure:
R
(z−ω)2
dzPout (y|z) (z − ω)2 e− 2V
1
2
+ ∂ω gout (ω, y, V ) + gout
(ω, y, V ) .
=
R
(z−ω)2
V
−
V 2 dzPout (y|z)e 2V

(143)

Using definition (142), and re-exponentiating the xi -dependent terms while keeping all the leading order terms, we
obtain (after normalization) that the iterative form of equation eq. (138) reads:
s
t
(Bµ→i
)2
x2
x
t
i At
√i
Atµ→i − 2N
µ→i +Bµ→i N − 2At
µ→i
mµ→i (t, xi ) =
e
(144)
2πN
with
t
t
t
Bµ→i
= gout (ωµ→i
, yµ , Vµ→i
) Fµi

t
t
2
and Atµ→i = −∂ω gout (ωµ→i
, yµ , Vµ→i
) Fµi
.

Given that mµ→i (t, xi ) can be written with a quadratic form in the exponential, we just write it as a Gaussian
distribution. We can now finally close the equations by writing (137) as a product of these Gaussians with the prior
mi→µ (xi ) ∝ PX (xi )e

−

(xi −Ri→µ )2
2Σi→µ

(145)

so that we can finally give the iterative form of the r-BP algorithm Algorithm 1 (below) where we denote by ∂ω (resp.
∂R ) the partial derivative with respect to variable ω (resp. R) and we define the “input” functions as
R
(x−R)2
dx x PX (x) e− 2Σ
fa (Σ, R) = R
,
fv (Σ, R) = Σ∂R fa (Σ, R) .
(146)
(x−R)2
dx PX (x) e− 2Σ
Algorithm 1 relaxed Belief-Propagation (r-BP)
Input: y
Initialize: ai→µ (t = 0),vi→µ (t = 0), t = 1
repeat
r-BP Update of {ωµ→i , Vµ→i }
Vµ→i (t) ←

X

2
Fµj
vj→µ (t − 1)

(147)

Fµj aj→µ (t − 1)

(148)

j6=i

ωµ→i (t) ←

X
j6=i

r-BP Update of {Aµ→i , Bµ→i }
Bµ→i (t) ← gout (ωµ→i (t), yµ , Vµ→i (t)) Fµi ,
2
Aµ→i (t) ← −∂ω gout (ωµ→i (t), yµ , Vµ→i (t))Fµi

(149)
(150)

r-BP Update of {Rµ→i , Σµ→i }
1
Aν→i (t)
X
Ri→µ (t) ← Σi→µ (t)
Bν→i (t)
Σi→µ (t) ← P

(151)

ν6=µ

(152)

ν6=µ

AMP Update of the estimated partial marginals ai→µ (t),vi→µ (t)
ai→µ (t) ← fa (Σi→µ (t), Ri→µ (t)) ,

(153)

vi→µ (t) ← fv (Σi→µ (t), Ri→µ (t)) .

(154)

t←t+1
until Convergence on ai→µ (t),vi→µ (t)
output: Estimated marginals mean and variances:

P
Bν→i
1
, Pν Aν→i
A
ν→i
ν
ν


P
Bν→i
fv P A1 ν→i , Pν Aν→i
ν
ν

ai ← fa
vi ←



P

,

(155)

.

(156)

62
2.

From r-BP to G-AMP

While one can implement and run the r-BP, it can be simplified further without changing the leading order behavior
of the marginals by realizing that the dependence of the messages on the target node is weak. This is exactly what
we have done to go from the standard belief propagation to the TAP equations in sec. IV E.
After the corresponding Taylor expansion the corrections add up into the so-called Onsager reaction terms [300].
The final G-AMP iterative equations are written in terms of means ai and their variances ci for each of the variables
xi . The whole derivation is done in a way that the leading order of the marginal ai are conserved. Given the BP was
asymptotically exact for the computations of the marginals so is the G-AMP in the computations of the means and
variances. Let us define
X
ωµt+1 =
Fµi ati→µ ,
(157)
i

Vµt+1

X

=

2 t
Fµi
vi→µ ,

(158)

i

1
t+1 ,
µ Aµ→i
P
t+1
µ Bµ→i
.
= P t
µ Aµ→i

Σt+1
=P
i

(159)

Rit+1

(160)

First we notice that the correction to V are small so that
X
X
2 t
2 t
Vµt+1 =
Fµi
vi→µ ≈
Fµi
vi .
i

(161)

i

From now on we will use the notation ≈ for equalities that are true only in the leading order.
From the definition of A and B we get
"
#−1 "
#−1
X
X
t+1
t+1
2
2
t+1
t+1
Σi = −
Fµi ∂ω gout (ωµ→i , yµ , Vµ→i )
≈ −
Fµi ∂ω gout (ωµ , yµ , Vµ )
µ

µ

#−1

"
Ri = −

X

2
Fµi
∂ω gout (ωµt+1 , yµ , Vµt+1 )

µ

×

"
X

#
t+1
t+1
Fµi gout (ωµ→i
, yµ , Vµ→i
)

µ

Then we write that
t+1
t+1
gout (ωµ→i
, yµ , Vµ→i
) ≈ gout (ωµt+1 , yµ , Vµt+1 ) − Fµi ati→µ ∂ω gout (ωµt+1 , yµ , Vµt+1 )

≈

gout (ωµt+1 , yµ , Vµt+1 )

−

Fµi ati ∂ω gout (ωµt+1 , yµ , Vµt+1 )

(162)
(163)

So that finally
#−1

"
(Σi )

t+1

= −

X

2
Fµi
∂ω gout (ωµt+1 , yµ , Vµt+1 )

(164)

µ

"
Rit+1

= (Σi )

t+1

×

#
X
µ

= ati + (Σi )t+1

Fµi gout (ωµ , yµ , Vµ ) −

X

2
Fµi
ai ∂ω gout (ωµ , yµ , Vµ )

Fµi gout (ωµt+1 , yµ , Vµt+1 )

(165)
(166)

µ

Now let us consider ai→µ :
t
t
ati→µ = fa (Ri→µ
, Σi→µ ) ≈ fa (Ri→µ
, Σi )

≈

≈

t
fa (Rit , Σi ) − Bµ→i
fv (Rit , Σi )
ati − gout (ωµt , yµ , Vµt )Fµi vit

(167)
(168)
(169)

63
So that
ωµt+1 =

X
i

Fµi ati −

X

2 t
gout (ωµt , yµ , Vµt )Fµi
vi =

X

i

i

Fµi ai − Vµt gout (ωµt , yµ , Vµt )

(170)

Algorithm 2 Generalized Approximate Message Passing (G-AMP)
Input: y
0
Initialize: a0 ,v0 , gout,µ
,t=1
repeat
AMP Update of ωµ , Vµ
Vµt ←

X

2 t−1
Fµi
vi

(171)

t−1
− Vµt gout,µ
Fµi at−1
i

(172)

i

ωµt ←

X
i

AMP Update of Σi , Ri , gout,µ

t
gout,µ
← gout (ωµt , yµ , Vµt )
"
#−1
X 2
t
t
t
Σi ← −
Fµi ∂ω gout (ωµ , yµ , Vµ )

(173)
(174)

µ

Rit ← at−1
+ Σti
i

X

t
Fµi gout,µ

(175)

µ

AMP Update of the estimated marginals ai , vi

ati ← fa (Σti , Rit )

(176)

vit ← fv (Σti , Rit )

(177)

t←t+1
until Convergence on a,v
output: a,v.

An important aspect of these equations to note is the index t − 1 in the Onsager reaction term in eq. (172) that is
crucial for convergence and appears for the same reason as in the TAP equations in sec. IV E. Note that the whole
algorithm is comfortably written in terms of matrix multiplications only, this is very useful for implementations where
fast linear algebra can be used.

3.

The potential

Using the usual Bethe approach, one can also compute the corresponding potential, and show that these message
passing equations actually aim to minimize the following Bethe free energy [174, 262, 306]:
X
X
F Bethe ({Σi }, {Ri }) =
DKL (Qi ||PX ) +
DKL (Mµ ||Pout )
i

µ


1 X
+
log 2πVµ∗ + 1 + Vµ∗ ∂ω gout (ωµ∗ , yµ , Vµ∗ ) ,
2 µ

(178)

R
with Vµ∗ and ωµ∗ satisfying their respective fixed-point conditions, and DKL (Q||P ) = [log Q(x) − log P (x)]dQ(x)
being the Kullback-Leibler divergences between two probability distributions. The two Kullback-Leibler divergences
above are taken with respect to probability distributions
Mµ (z; ωµ∗ , yµ , Vµ∗ ) =
Qi (x; Ti , Σi ) =

1
1
−
Pout (yµ |z) p
e
ZM
2πVµ∗
(x−Ti )2
1
−
PX (x)e 2Σi ,
ZQ

∗ )2
(z−ωµ
∗
2Vµ

(179)
(180)

where ZM and ZQ are the corresponding normalizations. Clearly these two probability distributions are related to
the denominators in (146) and (142).

64
4.

Note on the convergence and parameter learning

When G-AMP converges, its performance in terms of mean-squared error is usually considerably better than the
performance of other existing algorithms. The AMP equations are proven to converge in the idealized setting of large
systems, random iid matrices F with zero mean entries, and when the prior PX either corresponds to the empirical
distribution of the actual signal x∗ or to the convex relaxation of the problem [25, 26]. However, for mismatched
signal distributions, for non-random matrices, positive mean random matrices or small system-sizes they might have
unpleasant convergence issues.
Justified by the promise of a strikingly better performance, it is an important subject of current research how to
make the G-AMP algorithm the most robust as possible, while keeping its performance and speed. Several recent
works followed this line of research. In [52] it was identified why AMP has convergence problems for the simple case
of iid matrices with elements of non-zero mean. For means of the matrix elements larger than some critical value,
there is an anti-ferromagnetic-like (i.e. corresponding to a negative eigenvalue of the Hessian matrix) instability of
the Nishimori line. Finite size effects always cause a slight mismatch between the true and the empirical distribution
of the signal, and hence a slight deviation from the Nishimori line. The instability amplifies the deviation and causes
the iteration of the AMP equation to fail and go far away from the Nishimori line. Even though we know that in the
Bayes-optimal case the correct fixed point must satisfy the Nishimori conditions.
The simplest way to fix this problem for matrices of non-zero mean is the so-called mean removal as used in many
of the early implementations and discussed in [306]. However, the instability of the Nishimori line is a problem that
can cause failure of the algorithm in a wider range of settings than the non-zero mean matrices, e.g. for some of the
structured matrices. It is therefore useful to study other strategies to stabilize the Nishimori line in order to improve
the convergence. We observed that randomization of the update is doing just that. Combining the sequential update
while keeping track of the correct time-indices, Ref. [201] designed the so-called Swept-AMP (SwAMP). Ref. [175]
studied the expectation-maximization learning of parameters for the compressed sensing problem when the number
of parameters is small, typically finite while the system size grows. Since the parameter-learning is precisely imposing
the validity of the Nishimori conditions, this is another way that greatly stabilizes the Nishimori line and improves
the convergence.
An additional approach to stabilize the Nishimori line is to slow down the iterations via damping, i.e. taking a
linear combination of the new and old values when updating the messages. Such damping can, however, slow down
the iteration-process considerably. In order to avoid a drastic slow-down we take advantage of the following related
development.
The authors of [174] studied the Bethe free energy associated to G-AMP and wrote a variational form of the
Bethe free energy (178) for which the G-AMP fixed points are local minima (not only generic stationary points such
as saddles). This free energy can be optimized directly leading to a provably converging (but considerably slower)
method to find the fixed points. This variational Bethe free energy can be used to slow down the iterations in a way
to keep decreasing the free energy - this we call adaptative damping in [306]. Other groups suggested several closely
related convergence improvements [260, 262].
Combining the above implementation features - adaptive damping, sequential update, mean removal and parameter
learning - leads to the current state-of-the art implementation of the G-AMP algorithm for linear estimation problems
[1].

5.

Examples of priors and outputs

The G-AMP algorithm is written for a generic prior on the signal PX (as long as it is factorized over the elements)
and a generic element-wise output channel Pout . The algorithm depends on their specific form only trough the function
fa and gout defined by (146) and (142). It is useful to give a couple of explicit examples.
The sparse prior that is most commonly considered in probabilistic compressed sensing is the Gauss-Bernoulli prior,
that is when in (134) we have φ(x) = N (x, σ) Gaussian with mean x and variance σ. For this prior the input function
fa reads
(T −x)2

ρ e− 2(Σ+σ)
faGauss−Bernoulli (Σ, T ) =

(1 −

√
Σ
3

(xΣ + T σ)

(Σ+σ) 2
√
(T −x)2
T2
Σ
−
ρ)e 2Σ + ρ √Σ+σ
e− 2(Σ+σ)

,

(181)

The most commonly considered output channel is simply additive white Gaussian noise (AWGN) (132). The output

65
function then reads
AW GN
gout
(ω, y, V ) =

y−ω
.
∆+V

(182)

As we anticipated above, the example of linear estimation that was most broadly studied in statistical physics is
the case of the perceptron problem discussed in detail e.g. in [310]. In the perceptron problem each of the M N dimensional patterns Fµ is multiplied by a vector of synaptic weights xi in order to produce an output yµ according
to
yµ = 1

N
X

if

Fµi xi > κ ,

(183)

i=1

yµ = −1

otherwise ,

(184)

where κ is a threshold value independent of the pattern. The perceptron is designed to classify patterns, i.e. one
starts with a training set of patterns and their corresponding outputs yµ and aims to learn the weights xi in such a
way that the above relation between patterns and outputs is satisfied. To relate this to the linear estimation problem
above, let us consider the perceptron problem in the teacher-student scenario where the teacher perceptron generated
the output yµ using some ground-truth set of synaptic weights x∗i . The student perceptron knows only the patterns
and the outputs and aims to learn the weights. How many patterns are needed for the student to be able to learn the
synaptic weights reliably? What are efficient learning algorithms?
In the simplest case where the threshold is zero, κ = 0 one can redefine the patterns Fµi ← Fµi yµ in which case the
corresponding redefined output is yµ = 1. The output function in that case reads
ω2

perceptron
gout
(ω, V

e− 2V
1

,
)= √
2πV H − √ω
V

(185)

where
Z

∞

H(x) =
x

t2
dt
√ e− 2 .
2π

(186)

Note also that the perceptron problem formulated this way is closely related to what is called 1-bit compressed sensing
in the signal processing literature [43, 316].
In physics a case of a perceptron that was studied in detail is that of binary synaptic weights xi ∈ {±1} [102, 103,
117]. To take that into account in the G-AMP we consider the binary prior PX (x) = [δ(x − 1) + δ(x + 1)]/2 which
leads to the input function
 
T
fabinary (Σ, T ) = tanh
.
(187)
Σ
D.

Replica free energy, state evolution and phase transitions
1.

State Evolution

The performance of the G-AMP algorithm above can be analyzed in the large-size limit, for matrices with independent entries. The corresponding analysis was baptized state evolution in [82]. Rigorous results on the state evolution
are provided in [25, 26, 140] for both AMP and G-AMP.
A statistical physics derivation of the state evolution in the spirit of the cavity method is given in [175] and is
formally very related to the state evolution of the TAP equation in sec. IV E. The main point in the derivation of the
state evolution is not to start from the G-AMP algorithm but from its message passing form, eqs. (148-154). One
then recalls that in belief propagation we assumed conditional independence of incoming messages which translates
into independence of terms present in sums of an extensive number of terms in the message passing. One then utilizes
the central limit theorem, enabling a reduction to the mean and variance of these terms matters. One simply needs
to keep track of the evolution of this mean and variance.
In statistical physics terms, the state evolution exactly corresponds to the derivation of the replica symmetric
equations via the cavity method. Since the cavity and replica methods always yield the same result, the state
evolution equations are exactly equivalent to the equations stemming from the replica method for the corresponding

66
problem. For compressed sensing this is illustrated in [175]. For the perceptron problem all the replica symmetric
equations and related results, e.g. in [310], can be obtained as a special case of the G-AMP state evolution using the
perceptron output.
Let us sketch briefly the derivation of state evolution for G-AMP. We follow the derivation of [153] for the specific
case of linear estimation. In order to write the equations for the MMSE for a general output functions, we need to
rewrite the output function using y = h(z, w) where w is a random variable distributed according to PW (w). Then
Z
Pout (y|z) = dwPW (w)δ[y − h(z, w)] .
(188)
P 2
P
Now we consider the random variable Vµ = i Fµi
vi . On average its value is i vi /N whereas its is variance o(1).
It thus concentrates around its mean as N → ∞ and therefore
1 X
vi .
(189)
V t+1 =
N i
P
P
Let us now turn our attention to ωµ→i = j6=i Fµi ai→j and zµ→i = j6=i Fµi si . These quantities are, according the
the assumptions of the cavity method, sums of uncorrelated variables, so they must converge to (correlated) Gaussian
variables with co-variances reading
E[ω 2 ] = E[a2 ] = q ,
E[zω] = E[sa] = m ,
E[z 2 ] = E[s2 ] = q0 ,

(190)
(191)
(192)

where we defined the order parameters m and q as the average correlation between the estimator a and the groundtruth signal s, and as a self-correlation for the estimator a, respectively.
Let us now see how Ri behaves (note the difference between the letters ω and w):
X
X
Ri
=
Bµ→i =
Fµi gout (ωµ→i , yµ , Vµ→i )
Σi
µ
µ
X
X
=
Fµi gout (ωµ→i , h(
Fµj sj + Fµi si , wµ ), V )
µ

=

(194)

j6=i

X
µ

(193)

Fµi gout (ωµ→i , h(

X

Fµj sj , wµ ), V ) + si αm̂ ,

(195)

j6=i

where we define
m̂ = Eω,z,w [∂z gout (ω, h(z, w), V )] .

(196)

p
Ri
= N (0, 1) αq̂ + si αm̂
Σi

(197)

We further write

with N (0, 1) being a random Gaussian variables of zero mean and unit variance, and where
 2

q̂ = Eω,z,w gout
(ω, h(z, w), V )

(198)

while the computation of Σ simply gives, thanks to concentration,
1
= αq̂ .
Σ
Now we can close all equations by computing q and m from their definitions (190-191) as


q = Es ER,Σ fa2 (Σ, R)
m = Es [ER,Σ sfa (Σ, R)]

(199)

(200)
(201)

In the Bayes-optimal case, the Nishimori consition imposes q = m and thus provides further simplifications. In fact,
one can show that q̂ is also equal to m̂ [153]. In this case, the state evolution leads to the asymptotic Bayes-optimal
minimum mean squared error (MMSE) for the corresponding inference problem. Just as it did for the planted spin

67
glass or the clustering of networks above. Since in linear estimation we are working with a fully connected graphical
model the asymptotic analysis can be written in terms of only 2 scalar parameters [175].
Finally, writing explicitly the integrals, the state evolution reads


ξ2
Z
Z
ξ
e− 2 2
1
t+1
,x + √
,
(202)
m
= dx PX (x) dξ √ fa
αm̂t
2π
αm̂t
(z−p)2
p2
Z
Z
−
e− 2mt e 2(ρ0 −mt )
t
m̂ = − dw PW (w) dp dz p
∂p gout (p, h(z, w), ρ0 − mt ) ,
(203)
2π mt (ρ0 − mt )
where we denote ρ0 = N EF (F 2 )EPX (x2 ). For the sparse prior (134) where the non-zero elements have the second
moment equal to x2 and the matrix F having iid elements of zero mean and variance 1/N we get ρ0 = ρx2 . For the
usual AWGN output (132) the second equation simplifies considerably into
m̂t =

1
∆ + ρx2 − mt

.

(204)

The mean-squared error of the AMP estimator is then MSE = ρx2 − mt . To evaluate the performance of AMP
initialized in a way unrelated to the unknown signal, the state evolution is iterated until convergence starting from
mt=0 = ρ2 x2 , where x is the mean of the distribution φ(x) in (134).
2.

Free energy

The same procedure we use to derive the state evolution can be used to write a scalar expression for the Bethe
free energy (178). Alternatively, the same expression can be obtained using the replica method, as e.g. in [153]. The
average Bethe free energy for the generalized linear estimation problem can be written as f Bethe = minm,m̂ Φ(m, m̂)
where
R
R
√
√ 
√
√
Φ(m, m̂) =α dy Dξ DsPout (y|s ρ0 − m + ξ m) log DxPout (y|x ρ0 − m + ξ m)
hR
√
√ i
R
αm̂s2
αm̂x2
αmm̂
−
+ DξdsPX (s)e− 2 +ξs αm̂ log dxPX (x)e− 2 +ξx m̂
(205)
2
The free energy can be also written as a function of the MSE E. In the case of the AWGN channel with variance
∆, for instance, it reads f Bethe = minE Φ(E) where [175]:
α(ρx2 − E)
α α
Φ(E) = − − log (∆ + E) −
2
2
2(∆ + E)
Z

Z
Z
√
α
α
x
√
∆+E x(s− 2 )+ξx ∆+E
+ dsPX (s) Dξ log
dx e
PX (x) .

(206)

Actually, this formula is, again, very close to the one discussed in the context of CDMA by Tanaka [298] and Guo and
Verdú [116]. An equivalent form has been derived for compressed sensing by Verdú [315]. Much more can be done
with the replica method in this context (see, e.g. [261, 295]).
3.

Phase transitions

What we want to discuss here is again the appearance of the very generic phenomena we have discussed already in
the first chapter, and illustrated on Fig. 1: the presence of two different transitions when trying to solve the Bayesian
optimal problem.
In Fig. 15 we plot down to what measurement rate α = M/N AMP recovers a signal of density ρ for several distributions of non-zero signal elements φ(x). We observe a first order phase transition that depends on the distribution
of non-zeros.
Interestingly, the corresponding line is always better than the `1 transition, as shown in [80], where the distribution
maximizing the spinodal line for Bayes-optimal inference was computed. This is to be expected, since we might hope
that the spinodal for the Bayesian optimal problem is the “best” spinodal, though there is no generic proof of this.
One might wonder how a measurement noise is affecting this diagram? Ref. [175] further studied how the optimal
inference performance is influenced by measurement noise, mismatch between the prior and the actual signal distribution, how the performance changes under expectation maximization learning of the parameters, such as ρ or the mean

68

1

α=M/N

0.8
0.6
0.4
ℓ1
N(0,1)
N(1,1)
N(-1,0.1)+N(1,0.1)

0.2
0
0

0.2

0.4

0.6

0.8

1

ρ=K/N
FIG. 15. Phase diagram for the AMP reconstruction in the optimal Bayesian case when the signal model is matching the
empirical distribution of signal elements. The elements of the M × N measurement matrix F are iid variables with zero mean
and variance 1/N . The spinodal transition αs (ρ) is computed with the state evolution and plotted for the following signal
distributions: φ(x) = N (0, 1) (green), φ(x) = N (1, 1) (blue) φ(x) = [N (−1, 0.1) + N (1, 0.1)]/2 (magenta). The data are
compared to the Donoho-Tanner phase transition α`1 (ρ) (dashed) for `1 reconstruction that does not depend on the signal
distribution, and to the theoretical limit for exact reconstruction α = ρ (red). Figure taken from [175].

x and the variance σ of the signal distribution. Ref. [16] explored the role of approximately sparse signals when the
δ(xi ) in (134) is replaced by a narrow Gaussian. The influence of an uncertainty in the elements of the matrix F was
studied in [176]. As expected from statistical physics, where in the presence of disorder first order phase transitions
become weaker, second order or disappear entirely, the same happens to the first order phase transition in compressed
sensing under presence of various kinds of noise. Qualitative results are illustrated in [16, 175, 176].
Existence of phase transition is not restricted to compressed sensing. In fact, the very same phenomenology has
been discussed and identified as early as in the 90’s in the context of learning a rule via a perceptron with binary
weights [103, 210, 281, 310]. In the present setting, this corresponds to a binary prior for the signal, and a perceptron
channel given by eq. (185). In this case, it becomes possible to identify perfectly the signal when the ratio between the
number of examples M and the number of variables N is larger than α ≈ 1.245, a result first derived by Derrida and
Gardner [103]. It is only beyond α ≈ 1.49, however, that the G-AMP algorithm actually succeed in finding the hidden
rule (i.e. the planted weights). All these classic results can thus be recovered by the present, unifying, approach.
Note that, for this problem, the G-AMP (or TAP) equations has been written and proposed as early as in 1989 in
[210]. The more recent application of either relaxed belief propagation [45] or G-AMP to actual data [330] has also
been considered.
Another classical example of this approach with a similar phenomenology is given by the classical work of Tanaka
and Kabashima in the context of CDMA [150, 298].

E.

Spatial coupling

As discussed in section IV B, if nothing is changed about the graphical model corresponding to compressed sensing,
the 1st order phase transition presented in Fig. 15 causes a barrier that is conjectured to be unbeatable by any
polynomial algorithm. In an idealized setting of compressed sensing, however, the design of the measurement matrix
Fµi is entirely up to us. Thanks to this, compressed sensing belongs to the class of problems where spatial coupling
can be successfully applied. This was the main point of Ref. [173], where spatially coupled measurement matrices
were designed that allowed AMP to reconstruct the signal in compressed sensing anytime the measurement rate α
is larger than the density of non-zero elements in the signal ρ, i.e. all the way above the red line in Fig. 15. This
threshold saturation was soon after proven rigorously in [79]. We should note that spatial coupling for compressed
sensing was considered even before in [187], but in that work the observed improvement was limited and it did not

69
seem possible to achieve the information theoretical limit α = ρ.
A spatially coupled matrix suggested in [173] is illustrated in Fig. 16. The physical mechanism of why the spinodal
transition disappears with such spatially-coupled measurement matrices is exactly the same as for the spatially coupled
Currie-Weiss model discussed in section IV C. The first block of the matrix has a larger effective measurement rate and
plays the role of a nucleation seed. The interactions in the blocks next to the diagonal then ensure that the seed grows
into the whole system. Ref. [175] further discusses the choice of parameters of the spatially coupled measurement
matrices, other related designs and the effect of noise. With the goal of optimizing the spatially coupled matrices the
properties of the nucleation and necessary properties of the seed in compressed sensing were studied in [51].

0y1
B
B
B
B
@

0

C B
C B
C=B
C B
A @

F
1
J1

J2
1
J1

0

J2
1
J1

J2
1
J1

1 0s1

0
J2
1
J1

J2
1
J1

J2
1
J1

J2
1

: unit coupling
: coupling J1
: coupling J2
: no coupling (null elements)

C B
C B
C⇥ B
C B
A B
B
B
B
B
B
B
B
@

C
C
C
C
C
C
C
C
C
C
C
C
A

FIG. 16. Example of a spatially coupled measurement matrix F for compressed sensing. The elements of the matrix are taken
independently at random, their variance has the block diagonal structure as illustrated in the figure. The figure it taken from
[175].

There are several problems in the literature closely formally related to compressed sensing to which spatial coupling
was applied. One of them is, again, the CDMA defined in sec. VI A. When the spreading code is not random iid but
correspond to a spatially coupled random matrix then the performance of the CDMA system can be further improved
[277, 297].
Another problem formally very closely related to compressed sensing is the real valued error correction where the
codeword is corrupted by gross errors on a fraction of entries and a small noise on all the entries. In [18] it was shown
that the error correction and its robustness towards noise can be enhanced considerably using AMP and spatial
coupling.
A similar application of these ideas has allowed for the development of the spatially coupled sparse superposition
codes [14, 15, 19, 147, 270]. These are efficient coding schemes using exactly the compressed sensing setting with a
twist: the variables are L-dimensional and are forced to point into one corner of the L-hypercube. These codes, with
the associated AMP decoder, have been shown to be fast and reliable, and actually reach the Shannon limit for the
AWGN channel in the large L limit. It is actually expected that they are Shannon achieving for any memoryless
channel.

F.

Non-random matrices and non-separable priors
1.

Structured priors

AMP as described above uses the empirical distribution of signal elements. A natural question is whether its
performance can be improved further if more information is known about the signal. There are works that have
investigated this question and combined AMP with more complex, structured priors. Utilizing such structured priors
is the key to leveraging many of the advancements recently seen in statistical signal representation. For instance
techniques such as hybrid AMP [259, 278] have shown promising results. In those works, an additional graphical
model is introduced in order to take into account the correlation between variables. The improvement, for instance
for image denoising, can be spectacular.

70
Another possible approach to structured priors is not to model the correlations directly, but instead to utilize a
bipartite construction via hidden variables, as in the restricted Boltzmann machine (RBM) [128, 129]. If a binary
RBM can be trained to model the support pattern of a given signal class, then the statistical description of the
RBM admits its use within AMP, as was shown recently in [302]. This is particularly interesting since RBMs are
the building blocks of “deep belief networks” [29] and have recently sparked a surge of interest, partly because of
the efficient algorithms developed to train them (e.g. contrastive divergence (CD) [129]). The possibility, shown in
[302], to incorporate deep learned priors into generalized linear problems such as compressed sensing appears very
promising.
2.

Orthogonal matrices

Using non-random matrices for the linear projection is crucial to explore the range of applicability. The `1 minimization based approaches provably work well for a much larger class of matrices than random ones [123]. It
is therefore important to explore performance of the approach stemming from statistical physics for more general
matrices. There has been a collection of works in this direction, for instance in the analysis of the so-called perceptron
learning with correlated pattern (which form the sensing matrix) [152, 286].
As we have mentioned in Sec. IV E, there are many works in the context of the Ising model that have extended the
TAP approach for more generic random matrices. Similar contributions have recently emerged in the context of AMP
as well: the recent work on the so-called S-AMP algorithm that extends the approximate message-passing algorithm
to general matrix ensembles when there is a certain well-defined large system size limit, based on the so-called Stransform of the spectrum of the measurement matrix [48]. Older works, in fact, already used a similar approach for
the CDMA problem [296]. Analysis of performances in this direction is also doable [156, 195]. This shows that there
are still interesting open problems in this direction and we feel that there is still much more to be done in the context
of AMP beyond random matrices.
3.

Reconstruction in discrete tomography

Signal reconstruction in x-ray computed tomography (CT) is another problem that belongs to the class of linear estimation. The matrix Fµi corresponding to computed tomography has a very particular structure. Classical
reconstruction algorithms such as filtered back-projection require the corresponding linear system to be sufficiently
determined. Therefore, the number of measurements needed is generally close to the number of pixels to be reconstructed. Nevertheless, many applications would benefit from being able to reconstruct from a smaller number of
angles. Ref. [110] designed and tested a message passing algorithm that is adapted to the measures performed in
CT and is able to provide a reliable reconstruction with a number of measurements that is smaller than what was
required by previous techniques. This algorithm is another example of how inspiration from physics can provide a
promising algorithmic approach.
Each measurement in x-ray CT is well approximated by a sum of absorption coefficient of the material under
investigation along a straight line of the corresponding x-ray. Let us discretize the object under study into small
elements, and assign a discrete value xi ∈ {1, . . . , q} of the absorption coefficient to every small element i. One
component of the measurement in CT is then
X
yµ =
xi ,
(207)
i∈µ

where the sum is over all elements that lie on the line µ. The structural information that in general allows us to
decrease the number of measurement is that two pixels that lie next to each other mostly have the same value. This
expectation is represented in the following posterior distribution
 


M
P
Y
X
1
δ yµ −
P (x|y) =
xi  eJµ (ij)∈µ δxi ,xj  ,
(208)
Z µ=1
i∈µ
where the notation (ij) ∈ µ means the elements i and j are next to each other on the line corresponding to measurement
µ, Jµ is the interaction constant.
Ref. [110] writes the belief propagation algorithm for estimating marginals of this posterior distribution. However,
one update of this belief propagation involves intractable sums over all configurations of variables on one line. We
realize that these sums correspond to the partition functions of one-dimensional Potts model with fixed magnetization.

71
A one-dimensional line is a tree and for trees partition functions can be computed exactly with belief propagation.
We use this to design an algorithm that uses BP within BP. The resulting implementation is comparably fast to the
competing convex relaxation based algorithms and provides better reconstruction and robustness to noise. We believe
that overall the algorithm suggested in [110] provides a promising perspective for reconstruction in x-ray computed
tomography.

G.

An application in physics: Optics in complex media

We would like to close the loop with physics and mention one explicit physics experiment in the field of compressed
optics where AMP and related tools have been especially useful.
Wave propagation in complex media is a fundamental problem in physics, be it in acoustics, optics, or electromagnetism. In optics, it is particularly relevant for imaging applications. Indeed, when light passes through a
multiple-scattering medium, such as a biological tissue or a layer of paint, ballistic light is rapidly attenuated, preventing conventional imaging techniques, and random scattering events generate a so-called speckle pattern that is
usually considered useless for imaging.
Recently, however, wavefront shaping using spatial light modulators has emerged as a unique tool to manipulate
multiply scattered coherent light, for focusing or imaging in scattering media [230]. This makes it possible to measure
the so-called transmission matrix of the medium [254], which fully describes light propagation through the linear
medium, from the modulator de vice to the detector. Interestingly, these transmission matrices, due to the fact that
they account for the multiple scattering and interference, are essentially iid random matrices [254], which is exactly
what AMP techniques have been designed for.
An important issue in most such recent experiments lies in accessing the amplitude and phase of the output field,
that in optics usually requires a holographic measurement, i.e., the use of a second reference beam. The phase and
amplitude of the measured field can then be extracted by simple linear combinations of interference patterns with a
phase shifted or off-axis reference. This however poses the unavoidable experimental problem of the interferometric
stability of the reference arm. In order to be able to avoid entirely the use of the reference beam, we need to solve a
problem that looks like compressed sensing: finding x (sparse or binary) such that y = |F x| where F is complex and
iid. Combining the G-AMP and the SwAMP approach for phase retrieval [201, 279], it has been possible to precisely
do this, to calibrate efficiently, and to perform a reference-less measurement of the transmission matrix of a highly
scattering material. Within this framework, the work of [84] showed that the transfer matrix can still be retrieved
for an opaque material (a thin layer of white paint), and can then be used for imaging or focusing, using the opaque
complex material as a lens. This was experimentally validated.
As a result the full complex-valued transmission matrix of a strongly scattering material can be estimated, up to a
global phase, with a simple experimental setup involving only real-valued inputs and outputs in combination with the
techniques presented in this review. We believe this illustrates the power of the techniques presented in this review,
and shows that statistical physics does not only provide interesting perspective to inference problems, but also that
inference techniques can be useful for physics experiments in return.

H.

Beyond linear estimation: matrix factorization

The technics that have been used in the context of (generalized) linear estimation are not limited to this setting.
In fact, AMP can also be applied to the more challenging setting of matrix factorization.
Decomposing a matrix into a product of two matrices of given properties is another problem with an huge range
of applications. Mathematically stated, we are given a noisy measurement of a matrix Y of dimension M × N and
we aim to decompose it into a product Y = DX, where D is a M × R matrix, and X is a R × N matrix. Depending
on the application the problem is subject to various constraints such as sparsity of D or X, low rank R, or specific
models for the noise. Such constraints turn the matrix factorization into an algorithmically extremely challenging
problem. Let us list several examples of potential applications
• Dictionary learning: Many modern signal processing devices take a great advantage of the fact that almost every
class of signals of interest has a sparse representation, i.e. there exists a basis such that the observed signal can
be written as a sparse linear combination of atoms (i.e. columns of the corresponding matrix) from this basis.
Widely known, studied and used examples include the wavelet basis for images, of Fourier basis for acoustic
signals. For other examples of interest, especially those which utilize uncommon or novel data types, a suitable
basis for sparse representation may not yet be known. The goal of dictionary learning [85, 199, 200] is to learn
such a basis purely from samples of the data. Naturally the more samples are available the easier is this task.

72
The challenge is hence to design an algorithm that would be able to learn the dictionary efficiently from the
smallest possible number of samples.
In the mathematical statement of the problem the N columns of matrix Y represent the N samples of the data
for which we aim to find a sparse representation. M is then the dimension of each of the data points. The
dictionary learning problem is a decomposition of the matrix Y into a product Y = DX + W , where D is the
dictionary of R so-called atoms, and X is the sparse representation of the data, having only a fraction ρ of
non-zero elements in each column. W is a possible noise that may be taking into account the approximative
nature of the sparse representation.
• Feature learning in neural network: Arguably one of the most impressive engineering success of the last decade
is the design of the so called “deep-learned” neural networks [30, 130, 131]. Using these, computers are now
able to recognize people on video, to tell a dog from a cat on a picture, to process speech efficiently and even
to answer complicated questions. At the roots of this efficiency is the ability to learn features. The simplest
building block of deep neural networks can be seen as a matrix factorization problem with an output Pout (y|z)
corresponding to the activation functions.
• Blind source separation: A typical example of blind source separation [28, 149] are M sensors recording noise
at a party during time N . We then aim to separate the speech of the R participants of the party based on these
recordings. Recording at one of the sensors is an unknown linear combination of the speech of all the present
people. A particularly challenging, but also very interesting application-wise, is when the number of sensors is
smaller than the number of people, M < R. The individual speeches may then be in principle recovered only if
they can be represented as sparse in some known basis.
• Sparse PCA: In PCA the aim is to approximate the matrix Y by a low-rank matrix, i.e. columns of Y are
written as a linear combination of only few elements. This reduces the information needed to describe the data
contained in Y . This is usually done via the singular value decomposition (SVD). In some applications it is
desirable that the linear combination to be sparse [145, 331], thus further reducing the necessary information.
The sparsity constraint cannot be straightforwardly included in the SVD and hence alternative algorithmic
approaches are needed.
1.

Learning random sparse dictionaries

The matrix factorization problem can again be studied within the teacher-student scenario: The teacher generates
a the matrices F and X with random iid entries, from some probability distribution (possibly sparse) and provides
to the student element-wise measurement of the product F X via an output function Pout (Y |F X). Let us consider
the dimensions M, N, R to be all of the same order. Both the AMP [153, 177, 248, 305] and the replica approach
[153, 177, 275, 276] can be applied to such problems. In this the problem of identifiability of F and X from Y , the
corresponding sample complexity and the associated Bayes-optimal mean squared error can be derived. We refer to
these works for the phase diagrams for special cases of the problem such as the sparse dictionary learning, blind source
separation, sparse PCA, robust PCA or matrix completion.
These results are very striking because in a number of relevant cases (including the dictionary learning and source
separation) it is observed that the ideal performance predicted by this approach is far better than the one achieved
by current algorithms. This suggests that there is a lot to gain in algorithmic efficiency in such problems. Making
the algorithm of [153, 248] more robust along the lines described above for AMP for linear estimation and applying
it to more realistic data is a very promising subject of future work.
2.

Low rank decomposition

In some part of the applications mentioned above (e.g. the PCA) the dimension R, corresponding to the rank
of matrix Y , is considered as very small. In those cases it is reasonable to consider the limit where M and N
are comparable and both N, M → ∞ whereas R = O(1). The graphical model corresponding to this case has Rdimensional variables and pairwise interactions. Problems that can be represented via this graphical model include
variations of PCA, submatrix localization problem, biclustering, clustering.
Again AMP type of algorithm can be written for a general element-wise output function Pout (Y |F X), and the
associated state evolution leads to an analysis of the phase diagram, phase transitions and optimal mean squared
errors. Interestingly in this approach the prior distribution can be any R-variate distribution, thus allowing for
instance for priors representing membership to one or more clusters.

73
A problem corresponding to the rank one matrix factorization, clustering in the Gaussian mixture model, was studied
long time ago using the replica method [20, 35, 311]. It was presented as a prototypical example of unsupervised
learning, see e.g. chapter 8 in the textbook [86]. Note, however, that many more settings closer to current unsupervised
learning architectures should be analyzed. The revival of interest in low-rank matrix factorization was linked with the
development of the corresponding approximate message passing algorithm by Rangan and Fletcher for rank R = 1
[258], and Matsushita and Tanaka for general rank [206] followed by rigorous work of [76]. Recently, a version for
generic rank and generic output (as in G-AMP) has been presented and analysed in [193, 194]. The resulting algorithm
when tested numerically has good convergence properties and these works hence open the stage to further studies and
applications along the line discussed in this review.
As a final comment, it is to be noted that the phase transitions observed in the matrix factorization setting are
again of the same kind that we have described all along this review, schematically in Fig. 1. In many problems, the
presence of distinct algorithmic and information theoretic thresholds is observed. In this setting, it is particularly
interesting to characterize when spectral methods [11], the most widely used approach for low rank factorization, are
optimal, and when they are not. As shown in, e.g. [75, 193, 194, 223], there are a number of situations where they are
sub-optimal. This is a valuable information one can get from the physics approach, that should motivate the creation
of new algorithms approaching the best possible performance.

74
VII.

PERSPECTIVES

This review discusses recent progress in the statistical physics approach to understanding of different problems of
statistical inference. Just as theoretical physics often focuses on understanding of idealized models that represent
in a simplified way the salient features of realistic settings, we focus on the teacher-student setting under which the
Bayes-optimal inference can be analyzed without ambiguity, and the phase transition can be identified and discussed.
The main concept that we presented in detail is the mathematical relation between the Bayes-optimal inference and
properties of disordered systems on the co-called Nishimori line. We discussed in detail various phase transitions that
arrise in the Bayes-optimal inference and related them to statistical and computational thresholds in the corresponding
inference problems. The resulting picture, that we explain in detail on the simple example of the planted spin glass, is
relevant for a large class of inference problems. In later sections we illustrate this picture on the problem of clustering
of sparse networks and on the problem of generalized linear estimation. In today’s scientific landscape the number of
very interesting data-driven problems is considerable and, in the authors’ opinion, in a considerable number of them
the strategy presented in this review will bring exciting new results.
Let us recall here the main theme of the present review. For all the presented problems, the central scientific questions that one can attempt to answer are: (1) Under what conditions is the information contained in the measurements
sufficient for a satisfactory inference to be possible? (2) What are the most efficient algorithms for this task? In the
probabilistic setting discussed here, what is striking is that all of these problems have a very similar phenomenology,
that also appeared earlier when physicists studied error correcting codes and perceptron learning, and that often these
two questions can be precisely answered. When the amount of information is too low a successful inference of the
signal is not possible for any algorithm: the corresponding information is simply insufficient. On the other hand, for
large enough amount of data, inference is possible, and the two regime are separated by a sharp phase transition.
Finally, and perhaps more importantly, there is often (as in first order transition in physics), an intermediate regime
where a successful inference is in principal possible but algorithmically hard. We have shown many example where
these transitions can be computed. A first order phase transition is always linked to appearance of computational
hardness, whereas a second order phase transition is not.
There is a deep value in this knowledge: If these transitions were not known, one could not know how good are the
state-or-the-art algorithms. To give an example, it is the very knowledge of the precise Shannon bound that is driving
research in information theory: if we know that we could be better in principle, when it is worth working on better
algorithms. As we have discussed in the very last chapter, for instance, present day matrix factorization algorithms
are often very far from optimality.
It is thus natural that, next to the theoretical understanding, an important outcome of the research reviewed in
this manuscript are new algorithmic ideas that are applicable beyond the idealized setting for which the presented
theoretical results hold. This translation of theoretical physics calculations into algorithmic strategies is very fruitful,
and promising direction of research. What is, perhaps, more challenging is to push such algorithmic ideas into actual
applications. From a point of view of a physicist, some of the problems presented here may seem very applied, but
in reality the research presented here is positioned on the theoretical side of the fields of statistics, machine learning
or signal processing. There is still a large step towards actual applications. Surely this is a very exciting direction to
explore and many of the researchers active in this field do explore it.
Let us now discuss some of the limitations and challenges related to the results and concepts presented in this
review.
On the algorithmic side a large part of the algorithms that stem from the research presented in this review involve
some kind of message passing derived under various assumptions of correlation decay. These assumptions rarely
hold of real datasets and hence sometimes this causes a failure of the corresponding algorithm. In our opinion the
theoretical potential of these message passing algorithms is very large, but much work needs to be done in order to
render these algorithms more robust and suitable for a generic problem.
Computer science and related fields mostly aim at understanding the worst possible setting or to prove something
under the weakest possible conditions. Analysis in a setting such as the teacher-student scenario usually does not
seem generic enough. In machine learning and more applied parts of computer science the focus is on a couple of
benchmark-datasets (such as e.g. the MNIST database). Simple models that aim to represent reasonably well the
generic situation are a domain of physical sciences. This is how progress in physics is made, and in our opinion this
is transferable to other fields as well.
Physicists can use with large effectiveness tools and methods from statistical physics of disordered systems such as
glasses and spin glasses. The advantage of these methods is that, unlike any other existing generic method, they are
able to describe the Bayes-optimal probabilistic inference. The disadvantage is that they are not fully rigorous and
are therefore sometimes viewed with dismissal or skepticism by mathematically oriented researchers. On the other
hand, rigorous proofs of results obtained with the cavity and replica method are a great mathematical challenge and
many partial results were obtained that have led to new mathematical tools making the whole probability theory more

75
powerful. One of the goals of the presented research is to transfer the physics reasoning and insight to mathematicians
in order to help them to craft plausible proof strategies.
abbreviation definition section
meaning
MMO
IB4
maximum mean overlap
MMSE
IB4
minimum mean squared error
MAP
IB4
maximum a posterior estimator
MCMC
IC3
Monte Carlo Markov chain
BP
IC3
belief propagation
RS
II F
replica symmetry
RSB
II F
replica symmetry breaking
1RSB
II F
one-step replica symmetry breaking
d1RSB
II F
dynamical one-step replica symmetry breaking
FRSB
II F
full-step replica symmetry breaking
XOR-SAT
II H 2
XOR-satisisfiability problem
ER
III A
Erdős-Rényi
V-B
III C
Viana-Bray model
LDPC
III F
low-density parity-check codes
CSP
III F
constraint satisfaction problem
K-SAT
III F
K-satisfiability problem
PCA
IV A
principal component analysis
C-W
IV C
Curie-Weiss model
TAP
IV E
Thouless-Anderson-Palmer
SBM
VA
stochastic block model
AWGN
VI
additive white Gaussian noise
CDMA
VI A
code-division multiple access
LASSO
VI A
least absolute shrinkage and selection operator
AMP
VI C
approximate message passing
G-AMP
VI C
generalized approximate message passing
r-BP
VI C 1
relaxed belief propagation
RBM
VI F
restricted Boltzmann machine
CT
VI F 3
computed tomography
TABLE I. Glossary of abbreviations with the number of the section where the term is defined.

76
ACKNOWLEDGEMENTS

We would like to express our gratitude to all our colleagues with whom many of the results presented here have been
obtained. In particular, we wish to thank Maria-Chiara Angelini, Jean Barbier, Francesco Caltagirone, T. Castellani,
Michael Chertkov, Andrea Crisanti, Leticia F. Cugliandolo, Laurent Daudet, Aurelien Decelle, Angelique Drémeau,
Laura Foini, Silvio Franz, Sylvain Gigan, Emmanuelle Gouillart, Jacob E. Jensen, Yoshyiuki Kabashima, Brian Karrer,
Lukas Kroc, Srinivas Kudekar, Jorge Kurchan, Marc Lelarge, Antoine Liutkus, Thibault Lesieur, Luca Leuzzi, André
Manoel, David Martina, Marc Mézard, Andrea Montanari, Cristopher Moore, Richard G. Morris, Elchanan Mossel,
Joe Neeman, Mark Newman, Hidetoshi Nishimori, Boshra Rajaei, Sundeep Rangan, Joerg Reichardt, Federico RicciTersenghi, Alaa Saade, David Saad, Ayaka Sakata, Fran cois Sausset, Christian Schmidt, Christophe Schülke, Guilhem
Semerjian, Cosma R. Shalizi, David Sherrington, Alan Sly, Phil Schniter, Eric W. Tramel, Rudiger Urbanke, Massimo
Vergassola, Jeremy Vila, Xiaoran Yan, Sun Yifan, Francesco Zamponi, Riccardo Zecchina and Pan Zhang.
A part of this manuscript has been written during a visit to UC Berkeley as part of the semester long program
on ”Counting Complexity and Phase Transitions”, at the Simons Institute for the Theory of Computing, which we
thanks for the kind hospitality.
Finally, we also wish to thank Alena and Julie for their continuous support and for letting us work, if only from
time to time, on this review.

FUNDINGS

Part of the research leading to the presented results has received funding from the European Research Council
under the European Union’s 7th Framework Programme (FP/2007-2013/ERC Grant Agreement 307087-SPARCS).

[1] Gamp in matlab. http://gampmatlab.wikia.com/wiki/Generalized Approximate Message Passing.
[2] Emmanuel Abbe, Afonso S Bandeira, Annina Bracher, and Amit Singer. Decoding binary node labels from censored
edge measurements: Phase transition and efficient recovery. Network Science and Engineering, IEEE Transactions on,
1(1):10–22, 2014.
[3] Emmanuel Abbe and Andrea Montanari. Conditional random fields, planted constraint satisfaction and entropy concentration. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 332–346.
Springer, 2013.
[4] Dimitris Achlioptas and Amin Coja-Oghlan. Algorithmic barriers from phase transitions. In Foundations of Computer
Science, 2008. FOCS’08. IEEE 49th Annual IEEE Symposium on, pages 793–802. IEEE, 2008.
[5] Dimitris Achlioptas and Cristopher Moore. The asymptotic order of the random k-SAT threshold. In Foundations of
Computer Science, 2002. Proceedings. The 43rd Annual IEEE Symposium on, pages 779–788. IEEE, 2002.
[6] Dimitris Achlioptas and Yuval Peres. The threshold for random K-SAT is 2k log 2 − o(k). Journal of the American
Mathematical Society, 17(4):947–973, 2004.
[7] Lada A Adamic and Natalie Glance. The political blogosphere and the 2004 us election: divided they blog. In Proceedings
of the 3rd international workshop on Link discovery, pages 36–43. ACM, 2005.
[8] Noga Alon, Itai Benjamini, Eyal Lubetzky, and Sasha Sodin. Non-backtracking random walks mix faster. Communications
in Contemporary Mathematics, 9(04):585–603, 2007.
[9] Noga Alon, Michael Krivelevich, and Benny Sudakov. Finding a large hidden clique in a random graph. Random Structures
and Algorithms, 13(3-4):457–466, 1998.
[10] Maria Chiara Angelini, Francesco Caltagirone, Florent Krzakala, and Lenka Zdeborová. Spectral detection on sparse
hypergraphs. arXiv preprint arXiv:1507.04113, 2015.
[11] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample
covariance matrices. Annals of Probability, pages 1643–1697, 2005.
[12] Victor Bapst, Guilhem Semerjian, and Francesco Zamponi. Effect of quantum fluctuations on the coloring of random
graphs. Physical Review A, 87(4):042322, 2013.
[13] David Barber. Bayesian reasoning and machine learning. Cambridge University Press, 2012.
[14] Jean Barbier and Florent Krzakala. Replica analysis and approximate message passing decoder for superposition codes.
In Information Theory (ISIT), 2014 IEEE International Symposium on, pages 1494–1498. IEEE, 2014.
[15] Jean Barbier and Florent Krzakala. Approximate message-passing decoder and capacity-achieving sparse superposition
codes. arXiv preprint arXiv:1503.08040, 2015.
[16] Jean Barbier, Florent Krzakala, Marc Mézard, and Lenka Zdeborová. Compressed sensing of approximately-sparse signals:
Phase transitions and optimal reconstruction. In Communication, Control, and Computing (Allerton), 2012 50th Annual
Allerton Conference on, pages 800–807. IEEE, 2012.

77
[17] Jean Barbier, Florent Krzakala, Lenka Zdeborová, and Pan Zhang. The hard-core model on random graphs revisited. In
Journal of Physics: Conference Series, volume 473, page 012021. IOP Publishing, 2013.
[18] Jean Barbier, Florent Krzakala, Lenka Zdeborová, and Pan Zhang. Robust error correction for real-valued signals via
message-passing decoding and spatial coupling. In Information Theory Workshop (ITW), 2013 IEEE, pages 1–5. IEEE,
2013.
[19] Jean Barbier, Christophe Schülke, and Florent Krzakala. Approximate message-passing with spatially coupled structured
operators, with applications to compressed sensing and sparse superposition codes. Journal of Statistical Mechanics:
Theory and Experiment, 2015(5):P05013, 2015.
[20] N Barkai and Haim Sompolinsky. Statistical mechanics of the maximum-likelihood density estimation. Physical Review
E, 50(3):1766, 1994.
[21] Dror Baron, Shriram Sarvotham, and Richard G Baraniuk. Bayesian compressive sensing via belief propagation. Signal
Processing, IEEE Transactions on, 58(1):269–280, 2010.
[22] A. Barrat, S. Franz, and G. Parisi. Temperature evolution and bifurcation of metastables states in mean-field spin glasses,
with connections with structural glasses. J. Phys. A, 30:5593–5612, 1997.
[23] Wolfgang Barthel, Alexander K Hartmann, Michele Leone, Federico Ricci-Tersenghi, Martin Weigt, and Riccardo
Zecchina. Hiding solutions in random satisfiability problems: A statistical mechanics approach. Physical Review Letters,
88(18):188701, 2002.
[24] Hyman Bass. The Ihara-Selberg zeta function of a tree lattice. International Journal of Mathematics, 3(06):717–797,
1992.
[25] Mohsen Bayati, Marc Lelarge, Andrea Montanari, et al. Universality in polytope phase transitions and message passing
algorithms. The Annals of Applied Probability, 25(2):753–822, 2015.
[26] Mohsen Bayati and Andrea Montanari. The dynamics of message passing on dense graphs, with applications to compressed
sensing. Information Theory, IEEE Transactions on, 57(2):764–785, 2011.
[27] Matthew James Beal. Variational algorithms for approximate Bayesian inference. University of London, 2003.
[28] Adel Belouchrani, Karim Abed-Meraim, Jean-François Cardoso, and Eric Moulines. A blind source separation technique
using second-order statistics. Signal Processing, IEEE Transactions on, 45(2):434–444, 1997.
[29] Yoshua Bengio. Learning deep architectures for AI. Foundations and Trends R in Machine Learning, 2(1):1–127, 2009.
[30] Yoshua Bengio, Pascal Lamblin, Dan Popovici, Hugo Larochelle, et al. Greedy layer-wise training of deep networks.
Advances in neural information processing systems, 19:153, 2007.
[31] Quentin Berthet and Philippe Rigollet. Computational lower bounds for sparse pca. arXiv preprint arXiv:1304.0828,
2013.
[32] Ludovic Berthier and Giulio Biroli. Theoretical perspective on the glass transition and amorphous materials. Reviews of
Modern Physics, 83(2):587, 2011.
[33] H. A. Bethe. Statistical physics of superlattices. Proc. Roy. Soc. London A, 150:552–575, 1935.
[34] Hans A Bethe. Statistical theory of superlattices. Proceedings of the Royal Society of London. Series A, Mathematical
and Physical Sciences, 150(871):552–575, 1935.
[35] Michael Biehl and Andreas Mietzner. Statistical mechanics of unsupervised structure recognition. Journal of Physics A:
Mathematical and General, 27(6):1885, 1994.
[36] Giulio Biroli and Jean-Philippe Bouchaud. The random first-order transition theory of glasses: a critical assessment.
Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications, pages 31–113, 2012.
[37] B. Bollobás. Random graphs. Cambridge University Press, second edition, 2001.
[38] Béla Bollobás, Christian Borgs, Jennifer Chayes, Oliver Riordan, et al. Percolation on dense graph sequences. The Annals
of Probability, 38(1):150–183, 2010.
[39] Erwin Bolthausen. An iterative construction of solutions of the TAP equations for the Sherrington–Kirkpatrick model.
Communications in Mathematical Physics, 325(1):333–366, 2014.
[40] Charles Bordenave, Marc Lelarge, and Laurent Massoulié. Non-backtracking spectrum of random graphs: community
detection and non-regular ramanujan graphs. In Foundations of Computer Science (FOCS), 2015 IEEE 56th Annual
Symposium on, pages 1347–1357. IEEE, 2015.
[41] J.-P. Bouchaud, L. Cugliandolo, J. Kurchan, and M Mézard. Out of equilibrium dynamics in spin glasses and other glassy
systems. In A. P. Young, editor, Spin Glasses and Random Fields, pages 161–223. World Scientific, Singapore, 1998.
[42] Jean-Philippe Bouchaud and Marc Potters. Theory of financial risk and derivative pricing: from statistical physics to
risk management. Cambridge university press, 2003.
[43] Petros T Boufounos and Richard G Baraniuk. 1-bit compressive sensing. In Information Sciences and Systems, 2008.
CISS 2008. 42nd Annual Conference on, pages 16–21. IEEE, 2008.
[44] DR Bowman and K Levin. Spin-glass theory in the bethe approximation: Insights and problems. Physical Review B,
25(5):3438, 1982.
[45] A. Braunstein and R. Zecchina. Learning by Message Passing in Networks of Discrete Synapses. Physical Review Letters,
96(3):030201, 2006.
[46] Joseph D Bryngelson and Peter G Wolynes. Spin glasses and the statistical mechanics of protein folding. Proceedings of
the National Academy of Sciences, 84(21):7524–7528, 1987.
[47] T Tony Cai, Tengyuan Liang, and Alexander Rakhlin. Computational and statistical boundaries for submatrix localization
in a large noisy matrix. arXiv preprint arXiv:1502.01988, 2015.
[48] Burak Cakmak, Ole Winther, and Bernard H Fleury. S-amp: Approximate message passing for general matrix ensembles.
In Information Theory Workshop (ITW), 2014 IEEE, pages 192–196. IEEE, 2014.

78
[49] Francesco Caltagirone, Silvio Franz, Richard G Morris, and Lenka Zdeborová. Dynamics and termination cost of spatially
coupled mean-field models. Physical Review E, 89(1):012102, 2014.
[50] Francesco Caltagirone, Giorgio Parisi, and Tommaso Rizzo. Dynamical critical exponents for the mean-field potts glass.
Physical Review E, 85(5):051504, 2012.
[51] Francesco Caltagirone and Lenka Zdeborová. Properties of spatial coupling in compressed sensing. arXiv preprint
arXiv:1401.6380, 2014.
[52] Francesco Caltagirone, Lenka Zdeborová, and Florent Krzakala. On convergence of approximate message passing. In
Information Theory (ISIT), 2014 IEEE International Symposium on, pages 1812–1816. IEEE, 2014.
[53] Chiara Cammarota and Giulio Biroli. Ideal glass transitions by random pinning. Proceedings of the National Academy of
Sciences, 109(23):8850–8855, 2012.
[54] Emmanuel J Candè and Michael B Wakin. An introduction to compressive sampling. Signal Processing Magazine, IEEE,
25(2):21–30, 2008.
[55] Emmanuel J Candes and Terence Tao. Near-optimal signal recovery from random projections: Universal encoding
strategies? Information Theory, IEEE Transactions on, 52(12):5406–5425, 2006.
[56] JM Carlson, JT Chayes, L Chayes, JP Sethna, and DJ Thouless. Bethe lattice spin glass: the effects of a ferromagnetic
bias and external fields. i. bifurcation analysis. Journal of statistical physics, 61(5-6):987–1067, 1990.
[57] JM Carlson, JT Chayes, James P Sethna, and DJ Thouless. Bethe lattice spin glass: the effects of a ferromagnetic
bias and external fields. ii. magnetized spin-glass phase and the de almeida-thouless line. Journal of statistical physics,
61(5-6):1069–1084, 1990.
[58] Tommaso Castellani, Florent Krzakala, and Federico Ricci-Tersenghi. Spin glass models with ferromagnetically biased
couplings on the Bethe lattice: analytic solutions and numerical simulations. The European Physical Journal B-Condensed
Matter and Complex Systems, 47(1):99–108, 2005.
[59] Alain Celisse, Jean-Jacques Daudin, Laurent Pierre, et al. Consistency of maximum-likelihood and variational estimators
in the stochastic block model. Electronic Journal of Statistics, 6:1847–1899, 2012.
[60] Patrick Charbonneau, Yuliang Jin, Giorgio Parisi, Corrado Rainone, Beatriz Seoane, and Francesco Zamponi. Numerical
detection of the gardner transition in a mean-field glass former. Physical Review E, 92(1):012316, 2015.
[61] Patrick Charbonneau, Yuliang Jin, Giorgio Parisi, and Francesco Zamponi. Hopping and the stokes–einstein relation
breakdown in simple glass formers. Proceedings of the National Academy of Sciences, 111(42):15025–15030, 2014.
[62] Peter Cheeseman, Bob Kanefsky, and William M. Taylor. Where the Really Hard Problems Are. In Proc. 12th IJCAI,
pages 331–337, San Mateo, CA, USA, 1991. Morgan Kaufmann.
[63] Aaron Clauset, Mark EJ Newman, and Cristopher Moore. Finding community structure in very large networks. Physical
Review E, 70(6):066111, 2004.
[64] Amin Coja-Oghlan. Graph partitioning via adaptive spectral techniques. Combinatorics, Probability and Computing,
19(02):227–284, 2010.
[65] Amin Coja-Oghlan and André Lanka. Finding planted partitions in random graphs with general degree distributions.
SIAM Journal on Discrete Mathematics, 23(4):1682–1714, 2009.
[66] Amin Coja-Oghlan, Elchanan Mossel, and Dan Vilenchik. A spectral approach to analysing belief propagation for 3colouring. Combinatorics, Probability and Computing, 18(06):881–912, 2009.
[67] A Crisanti, H Horner, and H-J Sommers. The sphericalp-spin interaction spin-glass model. Zeitschrift für Physik B
Condensed Matter, 92(2):257–271, 1993.
[68] L. F. Cugliandolo and J. Kurchan. Analytical solution of the off-equilibrium dynamics of a long-range spin glass model.
Phys. Rev. Lett., 71:173, 1993.
[69] J. R. L. de Almeida and D. J. Thouless. Stability of the Sherrington-Kirkpatrick solution of a spin-glass model. J. Phys.
A, 11:983–990, 1978.
[70] Bruno De Finetti. Theory of Probability. A critical introductory treatment. John Wiley & Sons, 1979.
[71] Aurelien Decelle, Florent Krzakala, Cristopher Moore, and Lenka Zdeborová. Asymptotic analysis of the stochastic block
model for modular networks and its algorithmic applications. Physical Review E, 84(6):066106, 2011.
[72] Aurelien Decelle, Florent Krzakala, Cristopher Moore, and Lenka Zdeborová. Inference and phase transitions in the
detection of modules in sparse networks. Physical Review Letters, 107(6):065701, 2011.
[73] Amir Dembo, Andrea Montanari, et al. Ising models on locally tree-like graphs. The Annals of Applied Probability,
20(2):565–592, 2010.
[74] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm.
Journal of the royal statistical society. Series B (methodological), pages 1–38, 1977.
[75] Yash Deshpande and Andrea Montanari. Finding hidden cliques of size\ sqrt {N/e} in nearly linear time. Foundations
of Computational Mathematics, pages 1–60, 2015.
[76] Yateendra Deshpande and Andrea Montanari. Information-theoretically optimal sparse PCA. In Information Theory
(ISIT), 2014 IEEE International Symposium on, pages 2197–2201. IEEE, 2014.
[77] William E Donath and Alan J Hoffman. Lower bounds for the partitioning of graphs. IBM Journal of Research and
Development, 17(5):420–425, 1973.
[78] David L Donoho. Compressed sensing. Information Theory, IEEE Transactions on, 52(4):1289–1306, 2006.
[79] David L Donoho, Adel Javanmard, and Andrea Montanari. Information-theoretically optimal compressed sensing via
spatial coupling and approximate message passing. Information Theory, IEEE Transactions on, 59(11):7434–7464, 2013.
[80] David L Donoho, Iain Johnstone, and Andrea Montanari. Accurate prediction of phase transitions in compressed sensing
via a connection to minimax denoising. Information Theory, IEEE Transactions on, 59(6):3396–3433, 2013.

79
[81] David L Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. Proceedings
of the National Academy of Sciences, 106(45):18914–18919, 2009.
[82] David L. Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. Proc. Natl.
Acad. Sci., 106(45):18914–18919, 2009.
[83] David L Donoho and Jared Tanner. Sparse nonnegative solution of underdetermined linear equations by linear programming. Proceedings of the National Academy of Sciences of the United States of America, 102(27):9446–9451, 2005.
[84] Angélique Drémeau, Antoine Liutkus, David Martina, Ori Katz, Christophe Schülke, Florent Krzakala, Sylvain Gigan,
and Laurent Daudet. Reference-less measurement of the transmission matrix of a highly scattering material using a dmd
and phase retrieval techniques. Optics express, 23(9):11898–11911, 2015.
[85] Michael Elad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries.
Image Processing, IEEE Transactions on, 15(12):3736–3745, 2006.
[86] Andreas Engel and Christian P. L. Van den Broeck. Statistical Mechanics of Learning. Cambridge University Press, New
York, NY, USA, 2001.
[87] P. Erdős and A. Rényi. On random graphs. Publ. Math. Debrecen, 6:290–297, 1959.
[88] P. Erdős and A. Rényi. On the evolution of random graphs. Publ. Math. Inst. Hungar. Acad. Sci., 5:17–61, 1960.
[89] Vitaly Feldman, Will Perkins, and Santosh Vempala. On the complexity of random satisfiability problems with planted
solutions. In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, pages 77–86. ACM,
2015.
[90] A. J. Felstrom and K.S. Zigangirov. Time-varying periodic convolutional codes with low-density parity-check matrix.
IEEE Trans. Inf. Theory, 45(6):2181 –2191, 1999.
[91] Miroslav Fiedler. Algebraic connectivity of graphs. Czechoslovak mathematical journal, 23(2):298–305, 1973.
[92] Alyson K Fletcher, Sundeep Rangan, Lav R Varshney, and Aniruddha Bhargava. Neural reconstruction with approximate
message passing (neuramp). In Advances in neural information processing systems, pages 2555–2563, 2011.
[93] Laura Foini, Florent Krzakala, and Francesco Zamponi. On the relation between kinetically constrained models of
glass dynamics and the random first-order transition theory. Journal of Statistical Mechanics: Theory and Experiment,
2012(06):P06013, 2012.
[94] Santo Fortunato. Community detection in graphs. Physics Reports, 486(3):75–174, 2010.
[95] S. Franz, M. Mézard, F. Ricci-Tersenghi, M. Weigt, and R. Zecchina. A ferromagnet with a glass transition. Europhys.
Lett., 55:465, 2001.
[96] Silvio Franz and Giorgio Parisi. Phase diagram of coupled glassy systems: A mean-field study. Phys. Rev. Lett.,
79(13):2486–2489, Sep 1997.
[97] Silvio Franz and Giorgio Parisi. Quasi-equilibrium in glassy dynamics: an algebraic view. Journal of Statistical Mechanics:
Theory and Experiment, 2013(02):P02003, 2013.
[98] Joel Friedman. A proof of Alon’s second eigenvalue conjecture and related problems. American Mathematical Soc., 2008.
[99] R. G. Gallager. Information theory and reliable communication. John Wiley and Sons, New York, 1968.
[100] Robert G. Gallager. Low-density parity check codes. IEEE Trans. Inform. Theory, 8:21–28, 1962.
[101] Surya Ganguli and Haim Sompolinsky. Statistical mechanics of compressed sensing. Physical Review Letters,
104(18):188701, 2010.
[102] E Gardner and B Derrida. Optimal storage properties of neural network models. Journal of Physics A: Mathematical
and General, 21(1):271, 1988.
[103] E Gardner and B Derrida. Three unfinished works on the optimal storage capacity of networks. Journal of Physics A:
Mathematical and General, 22(12):1983, 1989.
[104] Alan E Gelfand and Adrian FM Smith. Sampling-based approaches to calculating marginal densities. Journal of the
American statistical association, 85(410):398–409, 1990.
[105] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.
IEEE Transactions on pattern analysis and machine intelligence, (6):721–741, 1984.
[106] A Georges, D Hansel, P Le Doussal, and J-P Bouchaud. Exact properties of spin glasses. ii. nishimori’s line: new results
and physical implications. Journal de Physique, 46(11):1827–1836, 1985.
[107] Antoine Georges and Jonathan S Yedidia. How to expand around mean-field theory using high-temperature expansions.
Journal of Physics A: Mathematical and General, 24(9):2173, 1991.
[108] Michelle Girvan and Mark EJ Newman. Community structure in social and biological networks. Proceedings of the
National Academy of Sciences, 99(12):7821–7826, 2002.
[109] Paul M Goldbart, Nigel Goldenfeld, and David Sherrington. Stealing the Gold-A Celebration of the Pioneering Physics
of Sam Edwards. Oxford University Press, 2005.
[110] Emmanuelle Gouillart, Florent Krzakala, Marc Mézard, and Lenka Zdeborová. Belief-propagation reconstruction for
discrete tomography. Inverse Problems, 29(3):035003, 2013.
[111] Peter Grassberger and Jean-Pierre Nadal. From statistical physics to statistical inference and back. Kluwer Academic,
1994.
[112] D. J. Gross, I. Kanter, and H. Sompolinsky. Mean-field theory of the potts glass. Phys. Rev. Lett., 55(3):304–307, Jul
1985.
[113] DJ Gross, I Kanter, and Haim Sompolinsky. Mean-field theory of the potts glass. Physical Review Letters, 55(3):304,
1985.
[114] D.J. Gross and M. Mézard. The simplest spin glass. Nucl. Phys. B, 240:431, 1984.

80
[115] Francesco Guerra and Fabio Lucio Toninelli. The thermodynamic limit in mean field spin glass models. Communications
in Mathematical Physics, 230(1):71–79, 2002.
[116] Dongning Guo and Sergio Verdú. Randomly spread cdma: Asymptotics via statistical physics. Information Theory, IEEE
Transactions on, 51(6):1983–2010, 2005.
[117] Géza Györgyi. First-order transition to perfect generalization in a neural network with binary synapses. Physical Review
A, 41(12):7097, 1990.
[118] Bruce Hajek, Yihong Wu, and Jiaming Xu. Computational lower bounds for community detection on random graphs.
arXiv preprint arXiv:1406.6625, 2014.
[119] Kathleen E Hamilton and Leonid P Pryadko. Tight lower bound for percolation threshold on an infinite graph. Physical
Review Letters, 113(20):208701, 2014.
[120] Ki-ichiro Hashimoto. Zeta functions of finite graphs and representations of p-adic groups. Automorphic forms and
geometry of arithmetic varieties., pages 211–280, 1989.
[121] S Hamed Hassani, Nicolas Macris, and Ruediger Urbanke. Chains of mean-field models. J. Stat. Mech.: Theor. and Exp.,
page P02011, 2012.
[122] Trevor Hastie, Robert Tibshirani, Jerome Friedman, and James Franklin. The elements of statistical learning: data
mining, inference and prediction. The Mathematical Intelligencer, 27(2):83–85, 2005.
[123] Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical learning with sparsity: the lasso and generalizations.
CRC Press, 2015.
[124] Matthew B Hastings. Community detection as an inference problem. Physical Review E, 74(3):035102, 2006.
[125] W Keith Hastings. Monte carlo sampling methods using markov chains and their applications. Biometrika, 57(1):97–109,
1970.
[126] Simon Heimlicher, Marc Lelarge, and Laurent Massoulié. Community detection in the labelled stochastic block model.
arXiv preprint arXiv:1209.2910, 2012.
[127] Itay Hen and AP Young. Exponential complexity of the quantum adiabatic algorithm for certain satisfiability problems.
Physical Review E, 84(6):061152, 2011.
[128] Geoffrey E. Hinton. Training products of experts by minimizing contrastive divergence. Neural Computation, 14(8):1771–
1800, August 2002.
[129] Geoffrey E. Hinton. A practical guide to training restricted boltzmann machines. Technical Report UTML TR 2010-003,
University of Toronto, Toronto, Canada, 2010.
[130] Geoffrey E Hinton, Simon Osindero, and Yee-Whye Teh. A fast learning algorithm for deep belief nets. Neural computation,
18(7):1527–1554, 2006.
[131] Geoffrey E Hinton and Ruslan R Salakhutdinov. Reducing the dimensionality of data with neural networks. Science,
313(5786):504–507, 2006.
[132] Glen M Hocky, Ludovic Berthier, and David R Reichman. Equilibrium ultrastable glasses produced by random pinning.
The Journal of chemical physics, 141(22):224503, 2014.
[133] Jake M Hofman and Chris H Wiggins. Bayesian approach to network modularity. Physical Review Letters, 100(25):258701,
2008.
[134] Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt. Stochastic blockmodels: First steps. Social networks,
5(2):109–137, 1983.
[135] J. J. Hopfield. Neural networks and physical systems with emergent collective computational properties. Proc. Nat. Acad.
Sci. USA, 79:2554–2558, 1982.
[136] John J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of
the national academy of sciences, 79(8):2554–2558, 1982.
[137] Yukito Iba. The nishimori line and bayesian statistics. Journal of Physics A: Mathematical and General, 32(21):3875,
1999.
[138] Gregg Jaeger. The ehrenfest classification of phase transitions: Introduction and evolution. Archive for History of Exact
Sciences, 53(1):51–81, 1998.
[139] S. Janson, T. Luczak, and A. Rucinski. Random graphs. Wiley, New-York, 2000.
[140] Adel Javanmard and Andrea Montanari. State evolution for general approximate message passing algorithms, with
applications to spatial coupling. Information and Inference, page iat004, 2013.
[141] Edwin T Jaynes. Information theory and statistical mechanics. Physical Review, 106(4):620, 1957.
[142] Mark Jerrum. Large cliques elude the metropolis process. Random Structures & Algorithms, 3(4):347–359, 1992.
[143] Mark Jerrum and Gregory B Sorkin. Simulated annealing for graph bisection. In Foundations of Computer Science, 1993.
Proceedings., 34th Annual Symposium on, pages 94–103. IEEE, 1993.
[144] Haixia Jia, Cristopher Moore, and Doug Strain. Generating hard satisfiable formulas by hiding solutions deceptively. In
Proceedings of the National Conference on Artificial Intelligence, volume 20, page 384. Menlo Park, CA; Cambridge, MA;
London; AAAI Press; MIT Press; 1999, 2005.
[145] Iain M Johnstone and Arthur Yu Lu. Sparse principal components analysis. 2004.
[146] Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul. An introduction to variational methods
for graphical models. Machine learning, 37(2):183–233, 1999.
[147] Antony Joseph and Andrew R Barron. Least squares superposition codes of moderate dictionary size are reliable at rates
up to capacity. Information Theory, IEEE Transactions on, 58(5):2541–2557, 2012.
[148] Ari Juels and Marcus Peinado. Hiding cliques for cryptographic security. Designs, Codes and Cryptography, 20(3):269–280,
2000.

81
[149] Tzyy-Ping Jung, Scott Makeig, Colin Humphries, Te-Won Lee, Martin J Mckeown, Vicente Iragui, and Terrence J
Sejnowski. Removing electroencephalographic artifacts by blind source separation. Psychophysiology, 37(02):163–178,
2000.
[150] Yoshiyuki Kabashima. A CDMA multiuser detection algorithm on the basis of belief propagation. Journal of Physics A:
Mathematical and General, 36(43):11111, 2003.
[151] Yoshiyuki Kabashima. Propagating beliefs in spin-glass models. Journal of the Physical Society of Japan, 72(7):1645–1649,
2003.
[152] Yoshiyuki Kabashima. Inference from correlated patterns: a unified theory for perceptron learning and linear vector
channels. In Journal of Physics: Conference Series, volume 95, page 012001. IOP Publishing, 2008.
[153] Yoshiyuki Kabashima, Florent Krzakala, Marc Mézard, Ayaka Sakata, and Lenka Zdeborová. Phase transitions and sample
complexity in bayes-optimal matrix factorization. to appear in IEEE Trans. Inf. Theory, preprint arxiv:1402.1298, 2016.
[154] Yoshiyuki Kabashima and David Saad. Belief propagation vs. tap for decoding corrupted messages. EPL (Europhysics
Letters), 44(5):668, 1998.
[155] Yoshiyuki Kabashima and Shinsuke Uda. A BP-based algorithm for performing Bayesian inference in large perceptrontype networks. In Algorithmic Learning Theory, pages 479–493. Springer, 2004.
[156] Yoshiyuki Kabashima and Mikko Vehkapera. Signal recovery using expectation consistent approximation for linear
observations. In Information Theory (ISIT), 2014 IEEE International Symposium on, pages 226–230. IEEE, 2014.
[157] Yoshiyuki Kabashima, Tadashi Wadayama, and Toshiyuki Tanaka. A typical reconstruction limit for compressed sensing
based on lp-norm minimization. Journal of Statistical Mechanics: Theory and Experiment, 2009(09):L09003, 2009.
[158] Ulugbek Kamilov, Sundeep Rangan, Michael Unser, and Alyson K Fletcher. Approximate message passing with consistent
parameter estimation and applications to sparse learning. In Advances in Neural Information Processing Systems, pages
2438–2446, 2012.
[159] David R Karger, Sewoong Oh, and Devavrat Shah. Iterative learning for reliable crowdsourcing systems. In Advances in
neural information processing systems, pages 1953–1961, 2011.
[160] Brian Karrer and Mark EJ Newman. Stochastic blockmodels and community structure in networks. Physical Review E,
83(1):016107, 2011.
[161] Brian Karrer, MEJ Newman, and Lenka Zdeborová. Percolation on sparse networks. Physical Review Letters,
113(20):208702, 2014.
[162] Raghunandan H Keshavan, Andrea Montanari, and Sewoong Oh. Matrix completion from a few entries. Information
Theory, IEEE Transactions on, 56(6):2980–2998, 2010.
[163] R. Kikuchi. A theory of cooperative phenomena. Phys. Rev., 81:988–1003, 1951.
[164] S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi. Optimization by simulated annealing. Science, 220:671–680, 1983.
[165] T.R. Kirkpatrick and D. Thirumalai. Dynamics of the structural glass transition and the p-spin-interaction spin-glass
model. Phys. Rev. Lett., 58:2091, 1987.
[166] T.R. Kirkpatrick and D. Thirumalai. p-spin-interaction spin-glass models: Connections with the structural glass problem.
Phys. Rev. B, 36:5388, 1987.
[167] T.R. Kirkpatrick, D. Thirumalai, and P.G. Wolynes. Scaling concepts for the dynamics of viscous liquids near an ideal
glassy state. Phys. Rev. A, 40:1045, 1989.
[168] TR Kirkpatrick and PG Wolynes. Connections between some kinetic and equilibrium theories of the glass transition.
Physical Review A, 35(7):3072, 1987.
[169] TR Kirkpatrick and PG Wolynes. Stable and metastable states in mean-field potts and structural glasses. Physical Review
B, 36(16):8552, 1987.
[170] W. Krauth and M. Mézard. Storage capacity of memory networks with binary couplings. J. Phys., 50:3057–3066, 1989.
[171] Valdis Krebs. Books about us politics. unpublished, http://www. orgnet. com, 2004.
[172] F. Krzakala and J. Kurchan. A landscape analysis of constraint satisfaction problems. Phys. Rev. E, 76:021122, 2007.
[173] F Krzakala, M Mézard, F Sausset, YF Sun, and L Zdeborová. Statistical-physics-based reconstruction in compressed
sensing. Physical Review X, 2(2):021005, 2012.
[174] Florent Krzakala, Andre Manoel, Eric W Tramel, and Lenka Zdeborová. Variational free energies for compressed sensing.
In Information Theory (ISIT), 2014 IEEE International Symposium on, pages 1499–1503. IEEE, 2014.
[175] Florent Krzakala, Marc Mézard, Francois Sausset, Yifan Sun, and Lenka Zdeborová. Probabilistic reconstruction in
compressed sensing: algorithms, phase diagrams, and threshold achieving matrices. Journal of Statistical Mechanics:
Theory and Experiment, 2012(08):P08009, 2012.
[176] Florent Krzakala, Marc Mézard, and Lenka Zdeborová. Compressed sensing under matrix uncertainty: Optimum thresholds and robust approximate message passing. In Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on, pages 5519–5523. IEEE, 2013.
[177] Florent Krzakala, Marc Mézard, and Lenka Zdeborová. Phase diagram and approximate message passing for blind
calibration and dictionary learning. In Information Theory Proceedings (ISIT), 2013 IEEE International Symposium on,
pages 659–663. IEEE, 2013.
[178] Florent Krzakala, Marc Mézard, and Lenka Zdeborová. Reweighted belief propagation and quiet planting for random
k-sat. Journal on Satisfiability, Boolean Modeling and Computation, 8:149, 2014.
[179] Florent Krzakala, Andrea Montanari, Federico Ricci-Tersenghi, Guilhem Semerjian, and Lenka Zdeborová. Gibbs states
and the set of solutions of random constraint satisfaction problems. Proc. Natl. Acad. Sci. U.S.A, 104:10318, 2007.
[180] Florent Krzakala, Cristopher Moore, Elchanan Mossel, Joe Neeman, Allan Sly, Lenka Zdeborová, and Pan Zhang. Spectral
redemption in clustering sparse networks. Proceedings of the National Academy of Sciences, 110(52):20935–20940, 2013.

82
[181] Florent Krzakala, Federico Ricci-Tersenghi, Lenka Zdeborová, Riccardo Zecchina, Eric W. Tramel, and Leticia F. Cugliandolo. Statistical Physics, Optimization, Inference and Message-Passing Algorithms. Oxford University Press, 2015.
[182] Florent Krzakala and Lenka Zdeborová. Hiding quiet solutions in random constraint satisfaction problems. Physical
Review Letters, 102(23):238701, 2009.
[183] Florent Krzakala and Lenka Zdeborová. Following gibbs states adiabatically - the energy landscape of mean-field glassy
systems. EPL (Europhysics Letters), 90(6):66002, 2010.
[184] Florent Krzakala and Lenka Zdeborová. On melting dynamics and the glass transition. i. glassy aspects of melting
dynamics. The Journal of chemical physics, 134(3):034512, 2011.
[185] Florent Krzakala and Lenka Zdeborová. On melting dynamics and the glass transition. ii. glassy dynamics as a melting
process. The Journal of chemical physics, 134(3):034513, 2011.
[186] Florent Krzakala and Lenka Zdeborová. Performance of simulated annealing in p-spin glasses. In Journal of Physics:
Conference Series, volume 473, page 012022. IOP Publishing, 2013.
[187] S. Kudekar and H.D. Pfister. The effect of spatial coupling on compressive sensing. In Communication, Control, and
Computing (Allerton), pages 347–353, 2010.
[188] Shrinivas Kudekar, Thomas J Richardson, and Rüdiger L Urbanke. Threshold saturation via spatial coupling: Why
convolutional ldpc ensembles perform so well over the bec. IEEE Transactions on Information Theory, 57(2):803–834,
2011.
[189] C. Kwon and D. J. Thouless. Ising spin glass at zero temperature on the bethe lattice. Phys. Rev. B, 37(13):7649–7654,
1988.
[190] Laurent Laloux, Pierre Cizeau, Jean-Philippe Bouchaud, and Marc Potters. Noise dressing of financial correlation matrices. Physical Review Letters, 83(7):1467, 1999.
[191] M. Lelarge, L. Massoulié, and Jiaming Xu. Reconstruction in the labeled stochastic block model. In Information Theory
Workshop (ITW), 2013 IEEE, pages 1–5, Sept 2013.
[192] Michael Lentmaier, Arvind Sridharan, Daniel J Costello, and K Sh Zigangirov. Iterative decoding threshold analysis for
ldpc convolutional codes. IEEE Trans. Inf. Theory, 56(10):5274–5289, 2010.
[193] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Mmse of probabilistic low-rank matrix estimation: Universality
with respect to the output channel. In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing
(Allerton), pages 680–687. IEEE, 2015.
[194] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Phase transitions in sparse pca. In Information Theory (ISIT),
2015 IEEE International Symposium on, pages 1635–1639. IEEE, 2015.
[195] Ting Liu, Chao-Kai Wen, Shi Jin, and Xiaohu You. Generalized turbo signal recovery for nonlinear measurements and
orthogonal sensing matrices. arXiv preprint arXiv:1512.04833, 2015.
[196] David Lusseau, Karsten Schneider, Oliver J Boisseau, Patti Haase, Elisabeth Slooten, and Steve M Dawson. The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations. Behavioral Ecology
and Sociobiology, 54(4):396–405, 2003.
[197] Zongming Ma, Yihong Wu, et al. Computational barriers in minimax submatrix detection. The Annals of Statistics,
43(3):1089–1116, 2015.
[198] D. J. C. MacKay. Information theory, inference, and learning algorithms. Cambridge University Press, Cambridge, 2003.
[199] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online dictionary learning for sparse coding. In Proceedings of the 26th annual international conference on machine learning, pages 689–696. ACM, 2009.
[200] Stéphane G Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dictionaries. Signal Processing, IEEE
Transactions on, 41(12):3397–3415, 1993.
[201] Andre Manoel, Florent Krzakala, Eric Tramel, and Lenka Zdeborovà. Swept approximate message passing for sparse
estimation. In Proceedings of the 32nd International Conference on Machine Learning (ICML-15), pages 1123–1132,
2015.
[202] Romain Mari and Jorge Kurchan. Dynamical transition of glasses: From exact to approximate. The Journal of chemical
physics, 135(12):124504, 2011.
[203] Manuel Sebastian Mariani, Giorgio Parisi, and Corrado Rainone. Calorimetric glass transition in a mean-field theory
approach. Proceedings of the National Academy of Sciences, 112(8):2361–2366, 2015.
[204] Travis Martin, Xiao Zhang, and MEJ Newman. Localization and centrality in networks. Physical Review E, 90(5):052808,
2014.
[205] Laurent Massoulié. Community detection thresholds and the weak Ramanujan property. In Proceedings of the 46th
Annual ACM Symposium on Theory of Computing, pages 694–703. ACM, 2014.
[206] Ryosuke Matsushita and Toshiyuki Tanaka. Low-rank matrix reconstruction and clustering via approximate message
passing. In Advances in Neural Information Processing Systems, pages 917–925, 2013.
[207] DC Mattis. Solvable spin systems with random interactions. Physics Letters A, 56(5):421–422, 1976.
[208] Sharon Bertsch McGrayne. The Theory That Would Not Die: How Bayes Rule Cracked the Enigma Code, Hunted Down
Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy. Yale University Press, 2011.
[209] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of State Calculations by Fast
Computing Machines. The Journal of Chemical Physics, 21:1087–1092, June 1953.
[210] M Mézard. The space of interactions in neural networks: Gardner’s computation with the cavity method. Journal of
Physics A: Mathematical and General, 22(12):2181, 1989.
[211] M. Mézard and G. Parisi. The bethe lattice spin glass revisited. Eur. Phys. J. B, 20:217, 2001.
[212] M. Mézard and G. Parisi. The cavity method at zero temperature. J. Stat. Phys., 111:1–34, 2003.

83
[213] M. Mézard, G. Parisi, and M. A. Virasoro. SK model: The replica solution without replicas. Europhys. Lett., 1:77–82,
1986.
[214] M. Mézard, G. Parisi, and M. A. Virasoro. Spin-Glass Theory and Beyond, volume 9 of Lecture Notes in Physics. World
Scientific, Singapore, 1987.
[215] M. Mézard, G. Parisi, and R. Zecchina. Analytic and algorithmic solution of random satisfiability problems. Science,
297:812–815, 2002.
[216] M. Mézard, F. Ricci-Tersenghi, and R. Zecchina. Alternative solutions to diluted p-spin models and XORSAT problems.
J. Stat. Phys., 111:505, 2003.
[217] Marc Mézard and Andrea Montanari. Information, physics, and computation. Oxford University Press, 2009.
[218] R. Monasson and R. Zecchina. Entropy of the K-satisfiability problem. Phys. Rev. Lett., 76:3881–3885, 1996.
[219] Rémi Monasson, Riccardo Zecchina, Scott Kirkpatrick, Bart Selman, and Lidror Troyansky. Determining computational
complexity from characteristic “phase transitions”. Nature, 400(6740):133–137, 1999.
[220] A. Montanari. Estimating random variables from random sparse observations. European Transactions on Telecommunications, 19(4):385–403, 2008.
[221] A. Montanari and G. Semerjian. On the dynamics of the glass transition on bethe lattices. J. Stat. Phys., 124:103–189,
2006.
[222] A. Montanari and G. Semerjian. Rigorous inequalities between length and time scales in glassy systems. J. Stat. Phys.,
125:23, 2006.
[223] Andrea Montanari, Daniel Reichman, and Ofer Zeitouni. On the limitation of spectral methods: From the gaussian hidden
clique problem to rank-one perturbations of gaussian tensors. In Advances in Neural Information Processing Systems,
pages 217–225, 2015.
[224] Andrea Montanari and Devavrat Shah. Counting good truth assignments of random k-SAT formulae. In Proceedings of
the eighteenth annual ACM-SIAM symposium on Discrete algorithms, pages 1255–1264. Society for Industrial and Applied
Mathematics, 2007.
[225] Andrea Montanari and Rudiger Urbanke. Modern coding theory: The statistical mechanics and computer science point
of view. Complex Systems, Les Houches lecture notes, page 69, 2007.
[226] Joris M Mooij, Hilbert J Kappen, et al. Validity estimates for loopy belief propagation on binary real-world networks. In
NIPS, 2004.
[227] Cristopher Moore and Stephan Mertens. The nature of computation. Oxford University Press, 2011.
[228] T. Mora. Géométrie et inférence dans l’optimisation et en théorie de l’information. PhD thesis, Université Paris-Sud,
2007. http://tel.archives-ouvertes.fr/tel-00175221/en/.
[229] Flaviano Morone and Hernán A Makse. Influence maximization in complex networks through optimal percolation. Nature,
524:6568, 2015.
[230] Allard P Mosk, Ad Lagendijk, Geoffroy Lerosey, and Mathias Fink. Controlling waves in space and time for imaging and
focusing in complex media. Nature photonics, 6(5):283–292, 2012.
[231] Elchanan Mossel, Joe Neeman, and Allan Sly. Stochastic block models and reconstruction. arXiv preprint arXiv:1202.1499,
2012.
[232] Elchanan Mossel, Joe Neeman, and Allan Sly. A proof of the block model threshold conjecture. arXiv preprint
arXiv:1311.4115, 2013.
[233] Elchanan Mossel, Joe Neeman, and Allan Sly. Belief propagation, robust reconstruction and optimal recovery of block
models. In Proceedings of The 27th Conference on Learning Theory, pages 356–370, 2014.
[234] R. Mulet, A. Pagnani, M. Weigt, and R. Zecchina. Coloring random graphs. Phys. Rev. Lett., 89:268701, 2002.
[235] Radford M Neal. Probabilistic inference using markov chain monte carlo methods. Technical report, Department of
Computer Science, University of Toronto Toronto, Ontario, Canada, 1993.
[236] Mark EJ Newman. Finding community structure in networks using the eigenvectors of matrices. Physical Review E,
74(3):036104, 2006.
[237] Mark EJ Newman. Modularity and community structure in networks. Proceedings of the National Academy of Sciences,
103(23):8577–8582, 2006.
[238] H. Nishimori. Statistical Physics of Spin Glasses and Information Processing: An Introduction. Oxford University Press,
Oxford, UK, 2001.
[239] Hidetoshi Nishimori. Exact results and critical properties of the ising model with competing interactions. Journal of
Physics C: Solid State Physics, 13(21):4071, 1980.
[240] Hidetoshi Nishimori. Internal energy, specific heat and correlation function of the bond-random ising model. Progress of
Theoretical Physics, 66(4):1169–1181, 1981.
[241] Hidetoshi Nishimori and KY Michael Wong. Statistical mechanics of image restoration and error-correcting codes. Physical
Review E, 60(1):132, 1999.
[242] Lars Onsager. Crystal statistics. i. a two-dimensional model with an order-disorder transition. Physical Review, 65(34):117, 1944.
[243] Manfred Opper, Burak Cakmak, and Ole Winther. A theory of solving tap equations for ising models with general
invariant random matrices. Journal of Physics A: Mathematical and Theoretical, 49(11):114002, 2016.
[244] Manfred Opper and David Saad. Advanced mean field methods: Theory and practice. MIT press, 2001.
[245] Manfred Opper and Ole Winther. Tractable approximations for probabilistic models: The adaptive thouless-andersonpalmer mean field approach. Physical Review Letters, 86(17):3695, 2001.

84
[246] A. Pagnani, G. Parisi, and M. Ratiéville. Near optimal configurations in mean field disordered systems. Phys. Rev. E,
68:046706, 2003.
[247] Giorgio Parisi and Marc Potters. Mean-field equations for spin models with orthogonal interaction matrices. Journal of
Physics A: Mathematical and General, 28(18):5267, 1995.
[248] Jason T Parker, Philip Schniter, and Volkan Cevher. Bilinear generalized approximate message passingpart i: Derivation.
Signal Processing, IEEE Transactions on, 62(22):5839–5853, 2014.
[249] J. Pearl. Reverend bayes on inference engines: A distributed hierarchical approach. In Proceedings American Association
of Artificial Intelligence National Conference on AI, pages 133–136, Pittsburgh, PA, USA, 1982.
[250] R Peierls. Statistical theory of superlattices with unequal concentrations of the components. In Proceedings of the Royal
Society of London A: Mathematical, Physical and Engineering Sciences, volume 154, pages 207–222. The Royal Society,
1936.
[251] Carsten Peterson and James R. Anderson. A mean field theory learning algorithm for neural networks. Complex Systems,
1:995 – 1019, 1987.
[252] Henry D Pfister and Pascal O Vontobel. On the relevance of graph covers and zeta functions for the analysis of SPA
decoding of cycle codes. In Information Theory Proceedings (ISIT), 2013 IEEE International Symposium on, pages
3000–3004. IEEE, 2013.
[253] T Plefka. Convergence condition of the tap equation for the infinite-ranged ising spin glass model. Journal of Physics A:
Mathematical and general, 15(6):1971, 1982.
[254] SM Popoff, G Lerosey, R Carminati, M Fink, AC Boccara, and S Gigan. Measuring the transmission matrix in optics:
an approach to the study and control of light propagation in disordered media. Physical Review Letters, 104(10):100601,
2010.
[255] Filippo Radicchi. Predicting percolation thresholds in networks. Physical Review E, 91(1):010801, 2015.
[256] Sundeep Rangan. Estimation with random linear mixing, belief propagation and compressed sensing. In Information
Sciences and Systems (CISS), 2010 44th Annual Conference on, pages 1–6. IEEE, 2010.
[257] Sundeep Rangan. Generalized approximate message passing for estimation with random linear mixing. In Information
Theory Proceedings (ISIT), 2011 IEEE International Symposium on, pages 2168–2172. IEEE, 2011.
[258] Sundeep Rangan and Alyson K Fletcher. Iterative estimation of constrained rank-one matrices in noise. In Information
Theory Proceedings (ISIT), 2012 IEEE International Symposium on, pages 1246–1250. IEEE, 2012.
[259] Sundeep Rangan, Alyson K Fletcher, Vivek K Goyal, and Philip Schniter. Hybrid approximate message passing with
applications to structured sparsity. arXiv preprint arXiv:1111.2581, 2011.
[260] Sundeep Rangan, Alyson K Fletcher, Philip Schniter, and Ulugbek S Kamilov. Inference for generalized linear models
via alternating directions and bethe free energy minimization. In Information Theory (ISIT), 2015 IEEE International
Symposium on, pages 1640–1644. IEEE, 2015.
[261] Sundeep Rangan, Vivek Goyal, and Alyson K Fletcher. Asymptotic analysis of map estimation via the replica method
and compressed sensing. In Advances in Neural Information Processing Systems 22, pages 1545–1553. 2009.
[262] Sundeep Rangan, Philip Schniter, and Alyson Fletcher. On the convergence of approximate message passing with arbitrary
matrices. In Information Theory (ISIT), 2014 IEEE International Symposium on, pages 236–240. IEEE, 2014.
[263] Jörg Reichardt and Stefan Bornholdt. Statistical mechanics of community detection. Physical Review E, 74(1):016110,
2006.
[264] F. Ricci-Tersenghi. The bethe approximation for solving the inverse ising problem: a comparison with other inference
methods. J. Stat. Mech.: Th. and Exp., page P08015, 2012.
[265] F. Ricci-Tersenghi, M. Weigt, and R. Zecchina. Simplest random k-satifiablitity problem. Phys. Rev. E, 63:026702, 2001.
[266] Emile Richard and Andrea Montanari. A statistical model for tensor PCA. In Advances in Neural Information Processing
Systems, pages 2897–2905, 2014.
[267] Tim Rogers. Assessing node risk and vulnerability in epidemics on networks. EPL (Europhysics Letters), 109(2):28005,
2015.
[268] Tim Rogers and Isaac Pérez Castillo. Cavity approach to the spectral density of non-Hermitian sparse matrices. Physical
Review E, 79(1):012101, 2009.
[269] Frank Rosenblatt. The perceptron: a probabilistic model for information storage and organization in the brain. Psychological Review, 65(6):386, 1958.
[270] Cynthia Rush, Adam Greig, and Ramji Venkataramanan. Capacity-achieving sparse regression codes via approximate
message passing decoding. In Information Theory (ISIT), 2015 IEEE International Symposium on, pages 2016–2020.
IEEE, 2015.
[271] Alaa Saade, Florent Krzakala, and Lenka Zdeborová. Spectral clustering of graphs with the bethe hessian. In Advances
in Neural Information Processing Systems, pages 406–414, 2014.
[272] Alaa Saade, Florent Krzakala, and Lenka Zdeborová. Spectral density of the non-backtracking operator on random
graphs. EPL (Europhysics Letters), 107(5):50005, 2014.
[273] Alaa Saade, Florent Krzakala, and Lenka Zdeborová. Matrix completion from fewer entries: Spectral detectability and
rank estimation. In Advances in Neural Information Processing Systems, pages 1261–1269, 2015.
[274] Alaa Saade, Marc Lelarge, Florent Krzakala, and Lenka Zdeborová. Spectral detection in the censored block model. In
2015 IEEE International Symposium on Information Theory (ISIT), pages 1184–1188. IEEE, 2015.
[275] Ayaka Sakata and Yoshiyuki Kabashima. Sample complexity of bayesian optimal dictionary learning. In Information
Theory Proceedings (ISIT), 2013 IEEE International Symposium on, pages 669–673. IEEE, 2013.

85
[276] Ayaka Sakata and Yoshiyuki Kabashima. Statistical mechanics of dictionary learning. EPL (Europhysics Letters),
103(2):28008, 2013.
[277] Christian Schlegel and Dmitri Truhachev. Multiple access demodulation in the lifted signal graph with spatial coupling.
In IEEE International Symposium on Information Theory Proceedings (ISIT), pages 2989–2993. IEEE, 2011.
[278] Philip Schniter. Turbo reconstruction of structured sparse signals. In Information Sciences and Systems (CISS), 2010
44th Annual Conference on, pages 1–6. IEEE, 2010.
[279] Philip Schniter and Sundeep Rangan. Compressive phase retrieval via generalized approximate message passing. Signal
Processing, IEEE Transactions on, 63(4):1043–1055, 2015.
[280] Bart Selman, David G. Mitchell, and Hector J. Levesque. Generating hard satisfiability problems. Artif. Intell., 81(12):17–29, 1996.
[281] HS Seung, Haim Sompolinsky, and N Tishby. Statistical mechanics of learning from examples. Physical Review A,
45(8):6056, 1992.
[282] Maoz Shamir and Haim Sompolinsky. Thouless-anderson-palmer equations for neural networks. Physical Review E,
61(2):1839, 2000.
[283] C. E. Shannon. A mathematical theory of communication. Bell System Tech. Journal, 27:379–423, 623–655, 1948.
[284] Ori Shental, Paul H Siegel, Jack K Wolf, Danny Bickson, and Danny Dolev. Gaussian belief propagation solver for systems
of linear equations. In Information Theory, 2008. ISIT 2008. IEEE International Symposium on, pages 1863–1867. IEEE,
2008.
[285] D. Sherrington and S. Kirkpatrick. Solvable model of a spin-glass. Phys. Rev. Lett., 35:1792–1796, 1975.
[286] Takashi Shinzato and Yoshiyuki Kabashima. Perceptron capacity revisited: classification ability for correlated patterns.
Journal of Physics A: Mathematical and Theoretical, 41(32):324013, 2008.
[287] Sasha Sodin. Random matrices, nonbacktracking walks, and orthogonal polynomials. Journal of Mathematical Physics,
48(12):123503, 2007.
[288] Subhojit Som and Philip Schniter. Compressive imaging using approximate message passing and a markov-tree prior.
Signal Processing, IEEE Transactions on, 60(7):3439–3448, 2012.
[289] H. Sompolinsky and Annette Zippelius. Relaxational dynamics of the edwards-anderson model and the mean-field theory
of spin-glasses. Phys. Rev. B, 25:6860–6875, Jun 1982.
[290] N. Sourlas. Spin-glass models as error-correcting codes. Nature, 339:693–694, 1989.
[291] N. Sourlas. Spin-glasses, error-correcting codes and finite-temperature decoding. Europhys. Lett., 25:159–164, 1994.
[292] N. Sourlas. Statistical mechanics and error-correcting codes. In P. Grassberger and J.-P. Nadal, editors, From Statistical
Physics to Statistical Inference and Back, pages 195–204. Kluwer, 1994.
[293] Erik B Sudderth, Alexander T Ihler, Michael Isard, William T Freeman, and Alan S Willsky. Nonparametric belief
propagation. Communications of the ACM, 53(10):95–103, 2010.
[294] YiFan Sun, Andrea Crisanti, Florent Krzakala, Luca Leuzzi, and Lenka Zdeborová. Following states in temperature in
the spherical s+ p-spin glass model. Journal of Statistical Mechanics: Theory and Experiment, 2012(07):P07002, 2012.
[295] Koujin Takeda and Yoshiyuki Kabashima. Statistical mechanical analysis of compressed sensing utilizing correlated
compression matrix. In Information Theory Proceedings (ISIT), 2010 IEEE International Symposium on, pages 1538–
1542. IEEE, 2010.
[296] Koujin Takeda, Shinsuke Uda, and Yoshiyuki Kabashima. Analysis of cdma systems that are characterized by eigenvalue
spectrum. EPL (Europhysics Letters), 76(6):1193, 2006.
[297] Keigo Takeuchi, Toshiyuki Tanaka, and Tsutomu Kawabata. Improvement of BP-based CDMA multiuser detection by
spatial coupling. In IEEE International Symposium on Information Theory Proceedings (ISIT), pages 1489–1493. IEEE,
2011.
[298] Toshiyuki Tanaka. A statistical-mechanics approach to large-system analysis of cdma multiuser detectors. Information
Theory, IEEE Transactions on, 48(11):2888–2910, 2002.
[299] D. J. Thouless. Spin-glass on a Bethe lattice. Phys. Rev. Lett., 56:1082–1085, 1986.
[300] D. J. Thouless, P. W. Anderson, and R. G. Palmer. Solution of ‘solvable model of a spin-glass’. Phil. Mag., 35:593–601,
1977.
[301] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B
(Methodological), pages 267–288, 1996.
[302] Eric W Tramel, Angélique Drémeau, and Florent Krzakala. Approximate message passing with restricted boltzmann
machine priors. Journal of Statistical Mechanics: Theory and Experiment, 2016(7):073401, 2016.
[303] Johannes D Van der Waals. The equation of state for gases and liquids. Nobel lectures in Physics, 1:254–265, 1910.
[304] L. Viana and A. J. Bray. Phase diagrams for dilute spin-glasses. J. Phys. C, 18:3037–3051, 1985.
[305] Jeremy Vila, Philip Schniter, and Joseph Meola. Hyperspectral image unmixing via bilinear generalized approximate
message passing. In SPIE Defense, Security, and Sensing, pages 87430Y–87430Y. International Society for Optics and
Photonics, 2013.
[306] Jeremy Vila, Philip Schniter, Sundeep Rangan, Florent Krzakala, and Lenka Zdeborová. Adaptive damping and mean
removal for the generalized approximate message passing algorithm. 2015 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP), pages 2021 – 2025, 2015.
[307] Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4):395–416, 2007.
[308] Martin J Wainwright and Michael I Jordan. Graphical models, exponential families, and variational inference. Foundations
and Trends R in Machine Learning, 1(1-2):1–305, 2008.
[309] Larry Wasserman. All of statistics: a concise course in statistical inference. Springer Science & Business Media, 2013.

86
[310] Timothy LH Watkin, Albrecht Rau, and Michael Biehl. The statistical mechanics of learning a rule. Reviews of Modern
Physics, 65(2):499, 1993.
[311] TLH Watkin and J-P Nadal. Optimal unsupervised learning. Journal of Physics A: Mathematical and General, 27(6):1899,
1994.
[312] J. Wehr and M. Aizenman. Fluctuations of extensive functions of quenched random couplings. J. Stat. Phys., 60:287–306,
1990.
[313] P Weiss and G Foëx. Magnetization of ferromagnetic substances above the curie point. Archives des sciences physiques
et naturelles, 31:5–19, 1911.
[314] Harrison C White, Scott A Boorman, and Ronald L Breiger. Social structure from multiple networks. i. blockmodels of
roles and positions. American journal of sociology, pages 730–780, 1976.
[315] Yihong Wu and Sergio Verdú. Optimal phase transitions in compressed sensing. Information Theory, IEEE Transactions
on, 58(10):6241–6263, 2012.
[316] Yingying Xu, Yoshiyuki Kabashima, and Lenka Zdeborová. Bayesian signal reconstruction for 1-bit compressed sensing.
Journal of Statistical Mechanics: Theory and Experiment, 2014(11):P11015, 2014.
[317] Xiaoran Yan, Cosma Shalizi, Jacob E Jensen, Florent Krzakala, Cristopher Moore, Lenka Zdeborová, Pan Zhang, and
Yaojia Zhu. Model selection for degree-corrected block models. Journal of Statistical Mechanics: Theory and Experiment,
2014(5):P05007, 2014.
[318] J.S. Yedidia, W.T. Freeman, and Y. Weiss. Understanding belief propagation and its generalizations. In Exploring
Artificial Intelligence in the New Millennium, pages 239–236. Science & Technology Books, 2003.
[319] Wayne W Zachary. An information flow model for conflict and fission in small groups. Journal of anthropological research,
pages 452–473, 1977.
[320] L. Zdeborová and F. Krzakala. Phase transitions in the coloring of random graphs. Phys. Rev. E, 76:031131, 2007.
[321] Lenka Zdeborová. Statistical physics of hard optimization problems1. Acta Physica Slovaca, 59(3):169–303, 2009.
[322] Lenka Zdeborová and Stefan Boettcher. A conjecture on the maximum cut and bisection width in random regular graphs.
Journal of Statistical Mechanics: Theory and Experiment, 2010(02):P02020, 2010.
[323] Lenka Zdeborová and Florent Krzakala. Generalization of the cavity method for adiabatic evolution of gibbs states.
Physical Review B, 81(22):224205, 2010.
[324] Lenka Zdeborová and Florent Krzakala. Quiet planting in the locked constraint satisfaction problems. SIAM Journal on
Discrete Mathematics, 25(2):750–770, 2011.
[325] Pan Zhang. Nonbacktracking operator for the Ising model and its applications in systems with multiple states. Physical
Review E, 91(4):042120, 2015.
[326] Pan Zhang, Florent Krzakala, Jörg Reichardt, and Lenka Zdeborová. Comparative study for inference of hidden classes
in stochastic block models. Journal of Statistical Mechanics: Theory and Experiment, 2012(12):P12021, 2012.
[327] Pan Zhang and Cristopher Moore. Scalable detection of statistically significant communities and hierarchies, using message
passing for modularity. Proceedings of the National Academy of Sciences, 111(51):18144–18149, 2014.
[328] Pan Zhang, Cristopher Moore, and MEJ Newman. Community detection in networks with unequal groups. Physical
Review E, 93(1):012303, 2016.
[329] Pan Zhang, Cristopher Moore, and Lenka Zdeborová. Phase transitions in semisupervised clustering of sparse networks.
Physical Review E, 90(5):052802, 2014.
[330] Justin Ziniel, Philip Schniter, and Per Sederberg. Binary linear classification and feature selection via generalized approximate message passing. In Information Sciences and Systems (CISS), 2014 48th Annual Conference on, pages 1–6.
IEEE, 2014.
[331] Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. Journal of computational and
graphical statistics, 15(2):265–286, 2006.

