# Sampling with flows, diffusion and autoregressive neural networks:
A spin-glass perspective

Davide Ghio  IdePHICS laboratory, École Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Switzerland  Yatin Dandi  IdePHICS laboratory, École
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Switzerland  SPOC
laboratory, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Switzerland  Florent Krzakala  IdePHICS laboratory, École Polytechnique
Fédérale de Lausanne (EPFL), CH-1015 Switzerland  Lenka Zdeborová  SPOC
laboratory, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Switzerland

###### Abstract

Recent years witnessed the development of powerful generative models based on
flows, diffusion or autoregressive neural networks, achieving remarkable
success in generating data from examples with applications in a broad range of
areas. A theoretical analysis of the performance and understanding of the
limitations of these methods remain, however, challenging. In this paper, we
undertake a step in this direction by analysing the efficiency of sampling by
these methods on a class of problems with a known probability distribution and
comparing it with the sampling performance of more traditional methods such as
the Monte Carlo Markov chain and Langevin dynamics. We focus on a class of
probability distribution widely studied in the statistical physics of
disordered systems that relate to spin glasses, statistical inference and
constraint satisfaction problems.

We leverage the fact that sampling via flow-based, diffusion-based or
autoregressive networks methods can be equivalently mapped to the analysis of
a Bayes optimal denoising of a modified probability measure. Our findings
demonstrate that these methods encounter difficulties in sampling stemming
from the presence of a first-order phase transition along the algorithm’s
denoising path. Our conclusions go both ways: we identify regions of
parameters where these methods are unable to sample efficiently, while that is
possible using standard Monte Carlo or Langevin approaches. We also identify
regions where the opposite happens: standard approaches are inefficient while
the discussed generative methods work well.

##  I Introduction

The field of machine learning recently witnessed the development of powerful
generative models able to produce new data-samples based on learning on
datasets of existing samples. Among the most prominent ones achieving recent
successes are flow-based models [1, 2, 3, 4, 5], diffusion-based models [6, 7,
8, 9] and generative autoregressive neural networks [10, 11, 12, 13]. These
approaches are achieving remarkable success in diverse areas such as image
generation [14], language modelling [15], generation of molecules [16, 17] or
theoretical physics [18]. A theoretical understanding of the capabilities of
these models, of their limitations and performance, remains, however, a
challenge. A major aspect of these techniques is in the learning of the
probability measure or its representation from examples. In this paper, we
focus instead on a restricted setting where the probability distribution we
aim to sample from is known beforehand. Our main goal is to contribute to
setting the theoretical understanding of the capabilities and limitations of
these powerful generative models.

When applied to parametric probability distributions, the generative models
are indeed designed with the goal to provide samples uniform at random from
the distribution. Here, we study whether they are able to efficiently sample
from types of probability distributions that are encountered in the study of
mean-field spin glasses and related statistical inference problems [19, 20,
21, 22, 23, 24].

There are several benefits to studying this class of probability
distributions. On the one hand, due to numerous studies in the statistical
physics of disordered systems, we possess a comparatively good grasp of
parameter regions where traditional sampling methods—like Monte Carlo sampling
or Langevin dynamics—are effective and where they are not [25, 26, 27, 28, 29,
30]. On the other hand, the tools available for outlining the phase diagrams
of these problems [20, 21, 31, 32] turn out to be highly effective in
analytically describing the performance of generative techniques such as flow-
based, diffusion-based, or autoregressive networks as samplers for the
respective probability measures.

The above-mentioned tools, along with their mathematically rigorous
counterparts, have recently been applied to the analysis and design of
sampling algorithms in mean-field spin glass models in the context of
stochastic localization [33, 34, 35, 36, 37, 38]. This was later found to have
a close relationship with diffusion-based models [38, 37]. In particular, [36,
37] showed how one can turn the message-passing denoising algorithms into
samplers, a technique we shall use as well in the present study. While we
could not find similar studies for flow-based models, closely related work on
autoregressive networks exists on the decimation of message-passing algorithms
used for finding solutions of the random K-satisfiability problem in [39, 40].
Here, we build on these works and bring the following contributions:

  * •

Using the formalism of stochastic interpolants [4, 41], we analyse sampling
with flow-based methods, which leads to Bayesian denoising with an additive
white Gaussian noise (AWGN). This turns out to be equivalent to what arises in
the analysis of diffusion-based sampling and stochastic localization derived
in [42, 36, 38].

  * •

In the case of autoregressive generative networks the analysis leads instead
to a Bayesian denoising problem now correcting erased variables (the Binary
Erasure Channel [43]) as in [40].

  * •

Focusing on prototypical exactly solvable models where one can perform
asymptotic analysis, we then study the phase diagrams of the corresponding
Bayesian denoising problems as a function of the strength of the signal-to-
noise ratio.

  * •

We investigate the feasibility of sampling using these techniques. We posit
that these methods are capable of efficient sampling, provided there is no
range in noise amplitude exhibiting the metastability characteristic of first-
order phase transitions. If such metastability exists, the denoising becomes
computationally hard.

  * •

We locate these metastable regions in prototypical models, specifically the
Ising and spherical $p$-spin models, the bicoloring of random hypergraphs
problem, and a matrix estimation problem with a sparse spike.

  * •

We also provide a GitHub repository with our numerical experiments in [44],
see https://github.com/IdePHICS/DiffSamp.

Figure 1: Schematic summary of the comparison of the efficiency of sampling
with flow-based, diffusion-based and autoregressive methods versus Langevin or
Monte-Carlo approaches in spin glass models with a random first-order
transition (top), and in statistical inference models with a computationally
hard phase (bottom). The computationally hard phase in inference problems
appears for For $\Delta_{\rm alg}<\Delta<\Delta_{\rm IT}$ where efficient
algorithms achieving close-to-optimal estimation error are not known and
conjectured not to exist [45].

In terms of comparison between the flow-based, diffusion-based or
autoregressive generative models with more traditional Monte Carlo Markov
chains (MCMC) and Langevin sampling procedures, the following picture emerges
(summarized in Fig. 1) for two classes of models.

  * •

Disordered models that exhibit a phase diagram of the random-first-order-
theory (RFOT) type, also often called discontinuous one-step replica symmetry
breaking [19, 46], are typical in the mean-field theory of glass transition
[28], but they also appear in a variety of random constraint satisfaction
problems [22]. For such models, there exists a so-called dynamical temperature
$T_{d}$ such that Monte Carlo Markov chains or Langevin algorithms are
predicted to sample efficiently for $T>T_{d}$ [25, 26, 27, 28, 29] while the
sampling problem is computationally hard for $T<T_{d}$. Recent work [47]
showed empirically that the region $T<T_{d}$ is also hard when sampling with
autoregressive neural networks. Our analysis of the currently used flow-based,
diffusion-based and autoregressive networks-based sampling procedures goes
further and reveals analytically that they perform worse than the traditional
techniques and that there is another temperature $T_{\rm tri}$, that depends
on the detail of the method, such that for $T_{d}<T<T_{\rm tri}$ these
generative methods do not sample efficiently while traditional procedures do.

  * •

A second class of problems are the noisy random statistical inference problems
with a hard phase and a statistical-to-computational gap [23, 24, 48, 45]. In
such problems, a good statistical inference can only be performed below a
critical value of the noise (or inverse signal-to-noise ratio) $\Delta_{\rm
IT}$, at least with infinite computational power. It turns out, however, that
the best algorithms we know for such problems (in particular message-passing
ones) are only able to work below a second threshold $\Delta_{\rm alg}$ and
fail to learn in the ”hard” region $\Delta_{\rm alg}<\Delta<\Delta_{\rm IT}$.
Existing literature predicts an even more significant gap to exist in those
cases for MCMC or Langevin-based samplers [30, 49, 50] down to a noise
amplitude $\Delta_{\rm MCMC}<\Delta_{\rm alg}$. In this setting, however, it
appears that flow-based, diffusion-based, and autoregressive network-based
methods outperform the standard approaches and sample efficiently as soon as
$\Delta<\Delta_{\rm alg}$.

###  I.1 Further related works

The present study of sampling algorithms owes a lot to recent lines of work.
In particular, we shall follow the helpful approach of continuous-time flows
based on stochastic interpolants [4, 41] as the starting point. Additionally,
our way of setting up the statistical physics equivalent problem for a
continuous-time flow method closely follows the one pioneered recently in [42,
36] for stochastic localization and generalized in [37] for diffusion models.

The difficulty in denoising we uncover turns out to be connected to an easy-
hard transition arising in statistical inference problems [23] often called
the statistical-to-computational gap [24, 48]. For these models, mean-field
algorithms from spin glass theory, such as the Belief-Propagation (BP)
equations [21] and the Approximate-Message-Passing (AMP) [31], or the
equivalent Thouless-Anderson-Palmer (TAP) [51] approach, are believed to be
among the most efficient one [52]. Using AMP for sampling via diffusion method
(or via stochastic localization) was discussed recently in [36, 37]. The
resulting algorithm turns out to have very close connections to the so-called
reinforcement BP [53], and a direct connection exists between the analysis in
the present paper and the reinforced and decimated version of message passing
algorithms [54, 55, 56, 39].

Finally, we note that since we assume that the model is known, we are not
discussing here the limitations of training denoisers models from a limited
number of data points. This difficulty was discussed recently in [57] in the
context of statistical physics models.

##  II Sampling with flow and diffusion-based models

We start by discussing continuous-time flow-based models [3, 5, 4]. The goal
is to sample a vector ${\bf x}_{0}\in{\mathbb{R}}^{N}$ from a so-called target
distribution $P_{0}({\bf x}_{0})$. Continuous-time flow-based models [3, 5,
58, 4] achieve this by constructing an ODE (flow) that continuously transforms
Gaussian noise ${\bf z}\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$ to a sample
${\bf x}_{0}\sim P_{0}$. One way to obtain such a flow is by inverting the
time-evolving law of stochastic or deterministic processes bridging from a
sample ${\bf x}_{0}\sim P_{0}$ to a Gaussian noise ${\bf
z}\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$. Such a continous-time reverse
process underlies a wide range of flow-models [5, 4] as well as the score-
based diffusion models [9]. In practice, we would only observe samples from
the distribution $P_{0}({\bf x}_{0})$, whereas in this paper we aim to focus
on the limitations of these procedures, we will thus assume that the
distribution $P_{0}({\bf x}_{0})$ is known as this can only render the task
easier than when only samples from $P_{0}({\bf x}_{0})$ are available.

To present this method clearly, we find it convenient to use the formalism of
stochastic linear interpolants introduced in [4, 41]. Starting from a vector
${\bf x}_{0}$ sampled from $P_{0}$ and a Gaussian vector ${\bf z}$, we
consider the process —or one-sided stochastic interpolant— defined by

${\bf y}(t)=\alpha(t){\bf x}_{0}+\beta(t){\bf z}\\\ ,\quad t\in[0,1]\,,$ (1)

where the functions $\alpha$ and $\beta$ are generic, but constrained by the
following relations:

$\displaystyle\alpha(0)=\beta(1)=1;\quad\alpha(1)=\beta(0)=0;$ (2)
$\displaystyle\forall t\in[0,1]:\alpha(t)\geq 0,\;\dot{\alpha}(t)\leq
0,\;\beta(t)\geq 0,\;\dot{\beta}(t)\geq 0$ (3)

such that indeed we have ${\bf y}(0)={\bf x}_{0}$ and ${\bf y}(1)={\bf z}$.
Suppose that $P_{0}({\bf x_{0}})$ admits a density $\rho_{0}$ w.r.t the
Lebesgue measure. It can be shown [41], see also Appendix E, that the
probability density $\rho({\bf y}(t))$ associated to the measure $P_{t}$ of
the random variable ${\bf y}(t)$ satisfies the following transport equation:

$\partial_{t}\rho({\bf y},t)+\nabla\cdot(b({\bf y},t)\rho({\bf y},t))=0\,,$
(4)

where we defined the velocity field

$b({\bf y},t)=\mathbb{E}[\partial_{t}{\bf y}(t)|{\bf y}(t)={\bf
y}]=\mathbb{E}[\dot{\alpha}(t){\bf x_{0}}+\dot{\beta}(t)z|{\bf y}(t)={\bf
y}]\,.$ (5)

Indeed, $b({\bf y},t)$ is simply the expected velocity of ${\bf y}(t)$
conditioned on being at ${\bf y}$ at time $t$. In Appendix E, we also provide
a formal definition of Eq. (5) based on the analysis in [4], along with a
discussion of the case when the measure $P_{0}({\bf x_{0}})$ is discrete.

Equation (1) defines a forward process interpolating from ${\bf x_{0}}$ to
${\bf z}$. Equation (4) further reveals that, in law, ${\bf y}(t)$ can be
obtained by applying the vector field (Eq. (5)) starting from ${\bf x_{0}}$ at
time $t=0$. The algorithm proposed by [4] relies on applying the velocity
field in the reverse direction of time. Concretely, it relies on approximating
the unique solution to the following ordinary differential equation starting
from a random Gaussian initial condition from $t=1$, back to $t=0$:

$\frac{d{\bf Y}(t)}{dt}=b({\bf Y}(t),t)\,.$ (6)

Eq. (4) then implies that the random variable defined by Eq. (6) solved
backwards in time from the final value ${\bf Y}_{t=1}={\bf
z}\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$ is also distributed according to
$\rho({\bf y}(t))$ (which is nothing but the continuity equation for the flow
defined by Eq. (6)) and, in particular, is distributed as the desired target
$P_{0}({\bf x}_{0})$ at time $t=0$. If we can numerically solve this ODE, then
we have a sampling algorithm for $P_{0}(x)$.

In order to do that, we discretize Eq. (6) using the forward Euler method (in
reverse time) to write

${\bf Y}_{t-\delta}={\bf Y}_{t}-\delta\,b({\bf Y}_{t},t)\,.$ (7)

Noticing that we can rewrite the vector field using Eq. (1) as

$\displaystyle b({\bf Y}_{t},t)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\mathbb{E}\left[\dot{\alpha}(t){\bf
x}_{0}+\frac{\dot{\beta}(t)}{\beta(t)}\beta(t){\bf z}\bigg{|}{\bf y}(t)={\bf
Y}_{t}\right]$ (8) $\displaystyle=$
$\displaystyle\\!\\!\\!\mathbb{E}\left[\dot{\alpha}(t){\bf
x}_{0}+\frac{\dot{\beta}(t)}{\beta(t)}\left({\bf y}(t)-\alpha(t){\bf
x}_{0}\right)\bigg{|}{\bf y}(t)={\bf Y}_{t}\right]$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{\dot{\beta}(t)}{\beta(t)}{\bf
Y}_{t}+\left(\dot{\alpha}(t)-\frac{\dot{\beta}(t)\alpha(t)}{\beta(t)}\right)\mathbb{E}\left[{\bf
x}_{0}|{\bf y}(t)={\bf Y}_{t}\right]\,,$

we put back Eq. (8) in Eq. (7) to reach

${\bf
Y}_{t-\delta}\\!=\\!\left(\\!1-\delta\frac{\dot{\beta}(t)}{\beta(t)}\\!\right)\\!{\bf
Y}_{t}-\delta\left(\\!\dot{\alpha}(t)-\frac{\dot{\beta}(t)\alpha(t)}{\beta(t)}\\!\right)\mathbb{E}\left[{\bf
x}_{0}|{\bf y}(t)={\bf Y}_{t}\right]$ (9)

which, given the initial condition ${\bf Y}_{t=1}={\bf z}$, will evolve back
to a sample from the target ${\bf Y}_{t=0}\sim P_{0}$, provided that we are
able to estimate $\mathbb{E}\left[{\bf x}_{0}|{\bf y}(t)={\bf Y}_{t}\right]$
well at each time $t$. In Algorithm 1 we report a schematic implementation of
the resulting flow-based sampling technique. In Appendix E, we further discuss
the effect of discretization and the associated sampling guarantees under
access to a perfect denoiser.

###  II.1 Diffusion-based models and SDEs

The flow-based algorithm [4] (Alg. 1) relies on replicating the law of the
interpolant ${\bf y}(t)$ defined by Eq. (1) through the deterministic ordinary
differential equation (ODE) defined by Eq. (6). Interestingly, however, the
same law for ${\bf y}(t)$ can be obtained via a stochastic differential
equation (SDE). Similarly, when ${\bf y}(t)$ is generated through a forward
diffusion process applied to the data, the law can be generated through an
associated time-reversed SDE (and an ODE) [59, 9]. It turns out that both
flow-based (ODE) methods and diffusion-based (SDE) ones associated to the
above processes (as well as stochastic localization [33, 34, 35]) require
estimating the posterior-mean $\mathbb{E}\left[{\bf x}_{0}|{\bf y}(t)={\bf
Y}_{t}\right]$, i.e. averages w.r.t a “tilted” probability measure.
Utilization of such an evolving tilted measure for sampling was done first in
the stochastic-localization-based sampling algorithm in [36, 37] and related
to generic diffusion models in [38]. The reliance of both the ODE and SDE-
based approaches on the same tilted measure is a consequence of the posterior
mean (or equivalently the score) yielding both the velocity field for the ODE
defined by Eq. (6) and the drift term for the equivalent SDE as noted in [4,
41, 9]. Since the law of ${\bf y}(t)$ remains the same with the ODE or SDE-
based approaches, the deterministic evolution of the “tilted” measure $P({\bf
x}|{\bf y}(t)={\bf Y}_{t})$ matches in law with that of the associated
stochastic localization process [38]. This is convenient for our analysis
since this means that all the phase transitions discussed in our work apply to
flow-based methods, diffusion-based ones and stochastic localization
approaches.

Algorithm 1 Flow-based sampling algorithm

Input: Denoiser, parameters: $\delta,\alpha(t),\beta(t)$

${\bf Y}_{t=1}={\bf Z}\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$, $N_{\rm
steps}\equiv\lfloor 1/\delta\rfloor$

for $l=0,1,\dots,N_{\rm steps}-1$ do $t=1-l\delta$ ; $\alpha=\alpha(t)$ ;
$\beta=\beta(t)$ ; $\gamma=\alpha/\beta$ Compute $\hat{\bf
x}_{0}=\mathbb{E}\left[{\bf x}_{0}|{\bf y}(t)={\bf Y}_{t}\right]$ using the
denoiser

(assuming a channel ${\bf Y}_{t}=\alpha{\bf X}_{0}+\beta{\bf Z}$) Update the
field $\bm{Y}_{t-\delta}$ via Eq. (9)

end for

Return $\hat{\bm{x}}_{0}$

###  II.2 The posterior average and Bayes-optimal denoising

The key difficulty is thus to estimate the average $\mathbb{E}\left[{\bf
x}_{0}|{\bf y}(t)={\bf Y}_{t}\right]$. This is nothing but the Bayes-optimal
denoising where ${\bf Y}_{t}$ is a noisy observation of ${\bf x}_{0}$
corrupted by additive white Gaussian noise in the form of Eq. (1).

The Bayes-optimal denoising estimator $\hat{\bf x}$ (optimal in the sense of
minimizing the mean-squared-error between the estimator and ${\bf x}_{0}$) is
to use the posterior mean $\hat{\bf x}=\mathbb{E}\left[{\bf x}_{0}|{\bf
y}(t)={\bf Y}_{t}\right]$ [60]. Determining this average can be
computationally hard, as it involves an integral in high-dimension. In
practical machine learning setups, one learns the denoiser using a neural
network trained on data coming from the distribution. Thus, the main narrative
of the diffusion and flow-based methods is that if one can access to a good
denoiser, then one can turn it into a sampler.

In this paper, we aim to focus on the limitations of this procedure and we
shall hence assume that the task of obtaining a denoiser has been completed to
near perfection by focusing on target distributions $P_{0}({\bf x}_{0})$
stemming from problems studied in statistical physics for which the Bayes-
optimal denoising problem including its algorithmic feasibility has been
extensively investigated. Before turning to these models, we write the
denoising problem in terms that are more familiar in the statistical physics
but also in the information-theoretic context.

Let us consider again the process defined by Eq. (1). Using the Bayes theorem,
we can see that the posterior distribution of the sample conditioned on the
observation ${\bf y}(t)={\bf Y}_{t}$ is given by

$\displaystyle P({\bf x}|{\bf y}(t)={\bf Y}_{t})=\frac{P_{0}({\bf x})P({\bf
y}(t)={\bf Y}_{t}|{\bf x})}{P({\bf y}(t)={\bf Y}_{t})}$ $\displaystyle=$
$\displaystyle\frac{1}{Z({\bf
Y}_{t})}\exp\left(\frac{\alpha(t)}{\beta(t)^{2}}\langle{\bf Y}_{t},{\bf
x}\rangle-\frac{\alpha(t)^{2}}{2\beta(t)^{2}}||{\bf x}||^{2}\right)P_{0}({\bf
x})\,.$

where in the last line we put all the terms not depending on ${\bf x}$ inside
the normalization aka partition function $Z$.

We now recall that the law of ${\bf Y}_{t}$ can be obtained through an
observation of a sample from the target distribution ${\bf x}_{0}\sim P_{0}$
through the AWGN channel rescaled by the factors $\alpha(t)$ and $\beta(t)$ as
in Eq. (1).

Here we denote ${\bf x_{0}}$ as the “truth” signal and make the difference
between the dummy variable ${\bf x}$ in the integral. We can hence rewrite the
measure in Eq. (II.2) further as

$\displaystyle P_{\gamma}({\bf x}|{\bf x}_{0},{\bf z})\\!\\!\\!\\!$
$\displaystyle\propto$ $\displaystyle\\!\\!\\!\\!e^{\gamma(t)^{2}\langle{\bf
x},{\bf x}_{0}\rangle+\gamma(t)\langle z,{\bf
x}\rangle-\frac{\gamma(t)^{2}}{2}\|{\bf x}\||^{2}}P_{0}({\bf x})\,,$ (11)

where $z\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$ and where we have defined the
effective signal-to-noise ratio

$\gamma(t)=\frac{\alpha(t)}{\beta(t)}.$ (12)

We see that, by construction, $\gamma$ will be such that $\gamma(1)=0$ and
$\gamma(0)=+\infty$. We will refer to $P_{\gamma}$ as the tilted measure, and
it will be a central object in the remaining discussion. We have added the
index $\gamma$ in the notation of the distribution $P_{\gamma}$ to distinguish
it from other considered probability distributions from now on.

##  III Autoregressive networks and ancestral sampling

Let us now discuss the second type of sampling algorithm considered in this
paper – the autoregressive networks. A classical way of sampling a vector
${\bf x}\in{\mathbb{R}}^{N}$ from a target distribution $P_{0}({\bf x})$, in
computer science and statistics is sequential sampling, or ancestral sampling:
One starts by computing the marginal probability of the first coordinate
$P_{0}(x_{1})$, and samples $x_{1}$ accordingly. Subsequently, the algorithm
generates all the coordinates for a sample, by sampling each coordinate
conditioned on its “parent” coordinates [61]. In the simplest case, the
distribution of each coordinate is assumed to depend on all previous
coordinates. Thus, after sampling $x_{1}$, for each subsequent node, one looks
at the distribution $P_{0}({\bf x}|x_{1})$, one considers the marginal
distribution of $x_{2}$ etc. Of course marginalization of a high-dimensional
probability distribution is in general hard, and the strategy used in an
autoregressive networks [62] is to directly learn from data a probability
distribution written in the (autoregressive) form $P_{0}({\bf
x})=P_{0}(x_{1})P_{0}(x_{2}|x_{1})\dots P_{0}(x_{N}|x_{1},\dots,x_{N-1})$ with
each term being represented via a neural network. While this decomposition
works for any ordering of the components, in practical applications the order
may be relevant (this is an important point in ancestral sampling). In the
present paper, we will consider the order to be random.

We showed in the previous section that sampling through diffusion in our
formalism boils down to the performance of a sequence of denoising problems,
where in information theory terms the signal is observed through an additive
white Gaussian noise channel. Analogously to that, we can interpret the
autoregressive networks, or its ideal version sequential sampling, as the
estimation of the marginal when a fraction $\theta$ of entries of one
configuration sampled uniformly at random from the target $P_{0}$ is revealed
exactly. To come back to the sampling scheme, reversing this process is
equivalent to fixing one variable at a time from the marginal conditioned to
the variables previously fixed until all variables are fixed. In statistical
physics, this procedure has been studied and analysed under the name
decimation [40]. In Algorithm 2 we report a schematic implementation of the
autoregressive-based sampling technique.

In information-theoretic terms, this corresponds to a denoising problem under
the so-called binary erasure channel (BEC) in which a transmitter sends a bit
and the receiver either receives the bit correctly or with some probability
$1-\theta$ receives a message that the bit was not received but “erased”
instead.

Repeating the way of thinking we used for sampling with diffusion, we now want
to analyse the denoising for such models and write the modified measure

$P_{\theta}({\bf x}|{\bf x}_{0},S_{\theta})\propto P_{0}({\bf x})\prod_{i\in
S_{\theta}}\delta([{x}_{i}-[{\bf x}_{0}]_{i})$ (13)

where $S_{\theta}$ is the set of revealed variables of size $\theta N$. We can
also think of these variables as pinned to their ground truth value ${\bf
x}_{0}$ and consequently call the measure $P_{\theta}({\bf x}|{\bf
x}_{0},S_{\theta})$ the pinning measure. From a statistical physics point of
view, the new measure can again be thought of as the original one, but with an
”infinite” magnetic field pointing to the direction of the components of ${\bf
x}_{0}$ for an expected fraction $\theta$ of components. We are again
interested in the marginals of $P_{\theta}$ that we will denote $\hat{\rm
x}(\theta)$. The evolution of the pinning measure can also be interpreted as a
coordinate-by-coordinate stochastic localization process [35].

Algorithm 2 Autoregressive-network-based sampling algorithm

Input: BEC denoiser

for $l=1,\dots,N$ do

Compute the posterior marginal $P_{l}(x_{l}|x_{1},\dots,x_{l-1})$ with BEC
denoiser Eq. (15)

Assign $x_{l}\sim P_{l}(x_{l}|x_{1},\dots,x_{l-1})$

end for

Return $\bm{x}$

##  IV Properties of Bayesian-optimal denoising

In both the cases – diffusion and autoregressive – we see that one has to
perform an optimal denoising based on an observation ${\bf y}$. In the case of
diffusion, we have access to an observable where a sample from the target
${\bf x}_{0}\sim P_{0}$ is polluted by an AWGN channel:

${\bf y}=\alpha{\bf x}_{0}+\beta{\bf z}\,,$ (14)

with $z\sim\mathcal{N}({\bf 0},\mathbb{I}_{N})$, while for autoregressive, it
is polluted by the BEC channel where for every $i=1,\dots,N$ independently

${y}_{i}=\begin{cases}[{\bf x}_{0}]_{i}&\text{with probability}\,\,\theta\\\
*&\text{otherwise}\end{cases}$ (15)

Our goal is now to study the properties of such channels and the corresponding
Bayes-optimal denoisers. A crucial point is that these Bayes-optimal
estimation problems lead to the Nishimori identities and the single state (or
replica symmetric in the spin glass jargon) properties of the measures
$P_{\gamma}$, Eq. (11), and $P_{\theta}$, Eq. (13), see e.g. the review [23].
While these are classical properties, we remind their derivations in the
appendix, and state informally their form and main consequences here.

Concretely, we shall be interested in the evolution in time (or equivalently
in $\gamma\in[0,\infty[$ (AWGN) and $\theta\in[0,1]$ (BEC)) of the following
order parameters, or overlaps:

$\displaystyle\mu(\gamma)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\mathbb{E}[\hat{\bf x}(\gamma)\cdot{\bf x}_{0}]\,,$
(16) $\displaystyle\chi(\gamma)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\mathbb{E}[\lVert\hat{\bf x}(\gamma)\rVert^{2}]\,.$
(17)

The same definitions hold if we replace $\gamma$ by $\theta$, i.e. for the
autoregressive process instead of the diffusive. The expectations are over the
disorder ${\bf x}_{0}$ and ${\bf z}$. The single state property relates to the
self-averaging of these order parameters: Almost anywhere in $\gamma$ (AWGN)
or $\theta$ (BEC), the equilibrium measure $P_{\gamma}$ and $P_{\theta}$ is a
single ”phase” probability measure, where overlap and order parameter are
self-averaging, i.e. concentrate to a single scalar quantity as $N\to\infty$.
Note that the “almost anywhere” is not a void statement: for instance, if one
is sampling a complicated glassy problem such as, say, the Sherrington-
Kirkpatrick model [20], strictly at $\gamma=0$, the Boltzmann measure is
complex and overlaps are not self-averaging. For almost all values of $\gamma$
(that is, except at some critical threshold values) the measure then
corresponds to a much simpler single-phase problem, and in particular, there
is no phenomenon akin to a static glassy transition or replica-symmetry-
breaking, thus the replica symmetric assumption is sufficient to describe the
thermodynamic behaviour of these measures. This will be instrumental in the
exact asymptotic analysis.

While such results have a long history [63, 64], their proof in the present
context follows from the analysis of optimal Bayesian denoising, and in
particular the I-MMSE theorem and the fluctuation-dissipation one [65, 66, 32,
67], for the AWGN channel, and by the study of the so-called pinning lemma
[68, 69] in the BEC channel, that also has roots in the study of the glass
transition [29, 70].

A second property, often called the Nishimori symmetry in physics [64, 23],
follows from Bayes theorem and states that for optimal Bayesian denoising one
has $\mu(\gamma)=\chi(\gamma)$ in the diffusion setting and
$\mu(\theta)=\chi(\theta)$ in the autoregressive one.

Note that while we cannot experimentally compute $\mu(\gamma)$ in simulations,
we can instead compute numerically $\chi(\gamma)$ which is just the norm of
the estimator. The study of phase transitions in such problems is thus reduced
to the behaviour of a single scalar quantity.

##  V Prototypical exactly analysable models

We will now analyse the properties of the tilted measure $P_{\gamma}$ and the
pinning measure $P_{\theta}$ for several concrete cases of the target measure
$P_{0}$. This will be possible exactly in the thermodynamic limit
$N\to\infty$. We shall focus on several classical problems from spin glass
theory, statistical inference, and constraint optimization problems, but
analogous analysis of the tilted and pinning measures can be done for many
other problems for which the phase diagram was obtained via the replica or the
cavity method [20, 21].

First, we shall start by studying models that present a so-called random
first-order transition [71, 28], or in replica parlance, a discontinuous one-
step replica symmetry breaking phenomenology [20]. There are two crucial
temperatures in RFOT systems. The so-called Kauzmann temperature $T_{K}$ below
which the system behaves as an ideal glass and the so-called dynamical
temperature $T_{d}$ that is defined as the temperature at which the point-to-
set correlation length diverges and consequently Monte Carlo or Langevin
dynamics equilibration time diverges as well [28, 29]. It is a widely accepted
conjecture that no efficient method can sample such models in their low
temperature $T<T_{d}$ glassy phase (and in fact, one can prove such hardness
for classes of algorithms [36, 72]). We shall focus on the “paramagnetic”
phases of these models $T>T_{d}$ for which the Monte Carlo or Langevin
equilibration time is known to be finite and hence efficient sampling with
MCMC or Langevin is possible. Specifically, we shall consider the following
models:

  * •

The Ising and spherical $p$-spin glass [19], whose Hamiltonian reads (for
$p=3$):

$\displaystyle{\cal H}({\bf
x})=-\frac{\sqrt{3}}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}\,,$ (18)

with $J_{ijk}\sim\mathcal{N}(0,1)$. The Boltzmann distribution is then
$P_{0}({\bf x})\propto\exp(-\beta{\cal H}({\bf x}))$ with $\beta=1/T$. In the
Ising case, we take $x_{i}=\pm 1$, for $i=1,\dots,N$, in the spherical case
${\bf x}\in{\cal S}^{N-1}$ (even though, since we are discussing the high-
temperature phase, we shall use the equivalent Gaussian model). This is one of
the most studied models in spin glass theory.

  * •

NAE-SAT, or bicoloring: Another class of popular models arises in the context
of constraint satisfaction problems, e.g. random satisfiability problem or
random graph coloring [22]. Here we shall focus on a prototypical case from
this class, the problem of coloring random $k$-hypergraphs with two colours.
This model was studied using statistical physics techniques in [73, 74, 75].
Numerous rigorous results for this model were also established e.g. in [76,
77]. The probability distribution is the following:

$P_{0}({\bf x})\propto\prod_{a=1}^{M}\omega({\bf x}_{\partial a})\,,\quad{\bf
x}\in\\{-1,+1\\}^{N}$ (19)
$\alpha=\frac{M}{N}\,,\quad\omega(x_{1},\dots,x_{k})=\begin{cases}0&\text{if}\;\sum_{i=1}^{k}x_{i}=\pm
k\\\ 1&\text{otherwise}\end{cases}$ (20)

where by ${\bf x}_{\partial a}$ we refer to the group of $k$ variables
entering into the clause $a$. Again, this model presents a RFOT phenomenology,
with a dynamical/clustering transition $\alpha_{d}$. The literature supports
the property that MCMC is able to sample efficiently for $\alpha<\alpha_{d}$
and is not for $\alpha>\alpha_{d}$ [22, 29, 78].

In these models, by studying the tilted/pinning measures, we shall see the
limitations of flow-based, diffusion and autoregressive sampling methods that
will fail at temperature $T_{d}<T<T_{\rm tri}$ and constraint density
$\alpha_{\rm tri}<\alpha<\alpha_{d}$. These methods thus turn out to be not as
performing as MCMC/Langevin approaches above the dynamical transition. We
expect this to be the case for any model with RFOT phenomenology. Every cloud,
however, has its silver lining, and we shall also note that there is a class
of models where flow, diffusion or autoregressive approaches outperform
MCMC/Langevin. This will be the case in statistical inference problems
presenting a hard phase, i.e. presenting a sharply defined region of the noise
$\Delta_{\rm IT}>\Delta>\Delta_{\rm alg}$ that is computationally hard for
message-passing algorithms and conjectured hard for any other efficient
algorithm [45]. A recent line of work [30, 49, 79, 50, 80], argues that when
it comes to sampling algorithms such as Langevin or other algorithms walking
in the space of configurations, such as gradient descent, then the hardness
actually extends to even beyond the threshold $\Delta_{\rm alg}$ up to some
$\Delta_{\rm MCMC}$ that depends on the specific algorithm. Yet, diffusion
models were shown to be able to sample down to $\Delta_{\rm alg}$ [37]. We
will also illustrate this here by studying the tilted and pinning measures of
the following prototypical model:

  * •

Sparse rank-one matrix factorization: The sparse spiked Wigner model and its
phase diagram were discussed e.g. in [81, 82, 83, 84]. This model is a
variation of the “planted” Sherrington-Kirkpatrick model in spin glass
physics. In such models, one is given a “spiked” version of a random symmetric
matrix with a rank-one perturbation with a “planted” vector ${\bf x^{*}}$ and
aims at finding back ${\bf x^{*}}$ from the matrix. This is done by sampling
from the posterior probability, and therefore one considers the probability
distribution:

$\displaystyle P_{0}({\bf x)}$
$\displaystyle\propto\prod_{i}P_{X}(x_{i})\prod_{i<j}\exp\left(\frac{1}{\Delta\sqrt{N}}J_{ij}x_{i}x_{j}-\frac{1}{2\Delta
N}x_{i}^{2}x_{j}^{2}\right)\,,$ (21) $\displaystyle
J_{ij}=\frac{x^{*}_{i}x^{*}_{j}}{\sqrt{N}}+z_{ij}\,,\quad
z_{ij}=z_{ji}\sim\mathcal{N}(0,\Delta)\,;$ $\displaystyle
P_{X}(x)=(1-\rho)\delta_{x,0}+\frac{\rho}{2}\left(\delta_{x,+1}+\delta_{x,-1}\right)\,,\quad
x^{*}_{i}\sim P_{X}\;\forall\,i.$

Flow-based and diffusion-based

Autoregressive-based

(a) Spherical $p$-spin (b) Ising $p$-spin (c) Rank-one estimation (d)
Bicoloring

Figure 2: Phase diagrams for the tilted measure $P_{\gamma}$ (top), and the
pinning measure $P_{\theta}$ (bottom): (a) The spherical $p$-spin model, b)
the Ising $p$-spin model, c) the sparse rank-one matrix estimation d) the
bicoloring problem on hypergraphs (NAE-SAT). The x-axis is the temperature
$T=1/\beta$ in (a) and (b), the inverse-SNR $\Delta/\rho^{2}$ in (c) and the
clauses-to-variables ratio $\alpha$ in (d), while the y-axis shows the SNR
ratio $\gamma^{2}=\alpha^{2}/\beta^{2}$ (top) and the decimated ratio $\theta$
(bottom). In the green phase there is a single maxima to the free entropy. The
red and orange regions display a phase coexistence with two maxima. In the red
region efficient denoising is predicted to be algorithmically hard. In (a) the
spherical $p$-spin at $\gamma=\theta=0$ the dynamical threshold
$T_{d}=\sqrt{3/8}$, the Kauzmann transition $T_{\rm K}\approx 0.58$, while the
tri-critical point is at $T_{\rm tri}=2/3$ for diffusion and $T_{\rm
tri}=\sqrt{1/2}$ for autoregressive. In (b) the Ising $p$-spin the values are:
$T_{d}\approx 0.682$, $T_{\rm K}\approx 0.652$, $T_{\rm tri}\approx 0.741$ for
diffusion, and $T_{\rm tri}\approx 0.759$ for autoregressive. In (c) for the
sparse rank-one matrix estimation at $\rho=0.08$ we have
$\Delta_{d}/\rho^{2}\approx 1.041$, $\Delta_{\rm K}/\rho^{2}\approx 1.029$,
and $\Delta_{\rm alg}/\rho^{2}\approx 0.981$. The tri-critical points are at
$\Delta_{\rm tri}/\rho^{2}\approx 1.08$ for diffusion, and $\Delta_{\rm
tri}/\rho^{2}\approx 1.069$ for autoregressive. In (d) the bicoloring the
values are $\alpha_{d}\approx 9.465$, $\alpha_{\rm K}\approx 10.3$. The tri-
critical points are $\alpha_{\rm tri}\approx 8.4$ for both diffusion and
autoregressive. The curves for bicoloring were obtained by a polynomial fit,
while in all the other cases we represent directly the data points.

##  VI Phase diagrams of the tilted and pinning measures

We first discuss sampling in the spherical $p$-spin model, Eq. (18). We focus
in particular on sampling in its paramagnetic phase, i.e. $T>T_{d}$, where
MCMC and Langevin algorithms are predicted to work efficiently. In order to
analyse flow-based, diffusion-based and autoregressive sampling, we need to
consider the corresponding denoising problems. For flow-based and diffusion-
based sampling, this leads to the tilted measure Eq. (11), with $P_{0}$ now
given by Eq. (18). Taken together, this defines a variant of the $p$-spin
model with a particular random field.

While the tilted measure may look complicated at first sight, because of the
random field in the direction of a particular “equilibrium” direction ${\bf
x}_{0}$, it turns out that measures of this type were already studied in the
literature. We refer the reader to Appendix A for details and only briefly
sketch the reasoning here. The trick (also used in e.g. [85, 36]) is to notice
that $\forall T>T_{K}$ the original $p$-spin model is contiguous to its
“planted” version, where a vector ${\bf x}_{0}$ has been hidden beforehand as
an equilibrium configuration (the spike tensor model [86]). The model is thus
equivalent for all practical purposes to its planted version, with the tilted
field now acting in the direction of the planted configuration.

Computing the phase diagram of the model associated with the tilted measure
thus requires computing the free entropy of the planted model with an
additional side Gaussian information. This is the same computation as needed
in various mathematical works based on Guerra’s interpolation technique where
the planted model is observed together with an additional Gaussian channel,
e.g. in [32]. This, and the single-state property of Bayes-optimal denoising
discussed in the former section (that ensures replica symmetry), allows us to
obtain the phase diagram of the tilted measure.

For the present spherical $p$-spin model, we show in Appendix B that the
equilibrium properties for $T>T_{\rm K}$ are given by

$\displaystyle\chi^{*}$ $\displaystyle=\underset{\chi}{\rm argmax}\,\Phi_{\rm
RS}(\chi)$ $\displaystyle\Phi_{\rm RS}(\chi)$
$\displaystyle=\frac{\widetilde{\chi}}{2}+\frac{1}{2}\log\left(\frac{2\pi}{\widetilde{\chi}+1}\right)-\frac{1}{2T^{2}}\chi^{3}\,,$
(22) $\displaystyle\,\widetilde{\chi}$
$\displaystyle=\frac{3}{2T^{2}}\chi^{2}+\gamma^{2}\,.$

Solving the above maximization problem is easily done. One observes that
depending on the range of parameters $T$ and $\gamma$, up to two local maxima
of $\Phi_{\rm RS}$ can be found for this model. In Fig. 2 (a, top) we depict
in green the regions of $T,\gamma$ where $\Phi_{\rm RS}$ has a unique
maximizer. The orange region is where two maximizers co-exist with the global
one having a smaller value of the order parameter $\chi$, and in the red
region, the global maximizer has a larger value of $\chi$. Such a
phenomenology is familiar in first-order phase transitions, where the red and
orange phases in Fig. 2 correspond to the phase coexistence region.

Now, a crucial point for the follow-up discussion of the flow- and diffusion-
based sampling is that it is widely believed to be algorithmically hard to
obtain the Bayes-optimal denoiser in the red region for all efficient
algorithms (even if $P_{0}$ is known) [23, 48, 52, 48, 45]. This is nothing
else than the well-known metastability problem in thermodynamics, as one is
”trapped” in the wrong maxima of (22).

The evidence for denoising being algorithmically hard in the red region goes
beyond mere physical analogies and is an intense subject of studies in
computer science, with the study of a variety of techniques such as low degree
polynomial or message passing algorithms (see e.g. [48]). Note that while
denoising is hard in the red region, direct sampling with MCMC is hard already
in the orange one.

To explain this further, let us discuss a concrete implementation of a
denoiser. Computing the marginals for the tilted measure is a classical topic
in spin glass and estimation theory, and the best-known algorithm to do it
efficiently is the so-called mean-field Thouless-Anderson-Palmer equations
[51], or —to use their modern counterpart— the iterative approximate message
passing (AMP) [31, 87]. The AMP algorithm is an iterative update procedure on
the estimates of the posterior means $\hat{x}_{i}$ and the co-variances
$\sigma_{i}$, and for the spherical 3-spin model it reads

$\begin{cases}B_{i}^{t}=\frac{\sqrt{3}\beta}{N}\sum_{j<k}J_{ijk}\widehat{x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\widehat{x}_{i}^{t-1}\sigma^{t}{\bf\widehat{x}}^{t}\cdot{\bf\widehat{x}}^{t-1}\\\
\widehat{x}_{i}^{t+1}=\frac{B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i}}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+\gamma^{2}+1}\,,\sigma^{t+1}=\frac{1}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+\gamma^{2}+1}\,.\\\ \end{cases}$ (23)

The virtue of AMP is that its performance can be tracked rigorously over
iteration time, and in fact, one can show that the overlap $\chi^{t}$ defined
by the AMP estimates obeys the following state evolution:

$\chi^{t+1}=\frac{\widetilde{\chi^{t}}}{1+\widetilde{\chi^{t}}}\,,\quad\widetilde{\chi}^{t}\equiv\frac{3}{2T^{2}}(\chi^{t})^{2}+\gamma^{2}\,,$
(24)

which is nothing but the fixed point equation of Eq. (22). We see now how the
presence of multiple fixed points can trap the algorithm in the wrong
maximizer.

Turning an AMP-denoiser into a sampler was precisely the idea introduced in
[42, 37] under the framework of stochastic localization [33, 34, 35]. While we
derived the equation for a flow-based approach, the conclusion of [42, 37]
remains. In particular, leveraging the rigorous analysis of the asymptotic
error obtained by AMP [88, 89], they prove that AMP can approximate optimal
denoising throughout the interpolation path in the regime where the global
maximum for $\Phi_{\rm RS}$ is the first one reached by the state evolution
[37]. Furthermore, using the local convexity of the TAP free energy [36, 90],
they show that the AMP iterates satisfy Lipschitz-continuity w.r.t the
observation ${\bf Y}_{t}$. This is crucial for the SDE-based sampling in [36,
37] as well as the continuous flow-based sampler in our case, since both
require control of the discretization errors. We discuss this further in
Appendix E. However, AMP and, conjecturally [45], any other polynomial-time
denoiser (and in particular any neural-network denoiser learned from data)
fail to return the correct marginal whenever the global maxima of $\Phi_{\rm
RS}$ is not the first one encountered by the state evolution, i.e. in the red
region in Fig. 2.

###  VI.1 Failure of generative models while MCMC succeeds

Recall now how the flow- or diffusion-based methods use the denoiser (the
marginal of the tilted measure) to produce samples from $P_{0}$. They start
with a Gaussian noise at $\gamma=0$ and increase $\gamma$ gradually while
transforming the Gaussian noise towards the direction of the marginal of the
tilted measure. To do this for a specific temperature $T$, one needs to be
able to denoise optimally on a vertical line at all values of
$\gamma\in[0,\infty[$. However, for temperatures below the tri-critical points
$T_{\rm tri}=2/3$ (this value is for diffusion, in the spherical $p$-spin
model) we see that we encounter the first-order phase transition, and the
metastable region (red in Fig. 2) where optimal denoising is computationally
hard and hence the uniform sampling using this strategy is as well. Another,
less critical, problem is the fact that the measure will change drastically at
the phase transition (as the value $\chi^{*}$ jumps discontinuously) which
means one has to be very careful with the discretization of the diffusion
process.

Based on this reasoning, we interpret the phase diagram of the tilted measure
as a representation of the presence of a fundamental barrier for sampling by
flow- and diffusion-based methods: The sampling scheme corresponds to starting
from $\gamma=0$ and going upwards. If this path intersects the hard phase
(red) at any time, it means that the Bayesian denoising cannot be performed
optimally in an efficient way, thus preventing the sampling scheme from
working efficiently in consequence. In particular, for the spherical $p$-spin
model, we computed the curves analytically, see Appendix B, and found that
this hurdle is present at temperatures up to the tri-critical point $T_{\rm
tri}=2/3$, strictly larger than the dynamical temperature $T_{d}=\sqrt{3/8}$
(threshold between the orange and red region at $\gamma=0$) down to which MCMC
and Langevin algorithms are predicted to sample efficiently in the literature.

Indeed, one can write exact equations describing the Langevin dynamics [25]
(and prove them [27]). The analysis shows that Langevin is efficient at
sampling for all $T>T_{d}$ [91]. On the other hand, for $T<T_{d}$, we are in
the so-called ”dynamical” spin glass phase and Langevin fails to sample in
linear time, this is the ageing regime [25, 26]. In fact, all polynomial
algorithms are conjectured to fail to sample efficiently, as can be proven for
any ”stable” algorithm [72].

The situation for sampling with autoregressive network, analysed via the phase
diagram of the pinning measure, is very analogous, see the lower row of Fig.
2. The pinning measure again defines a variation of the original $p$-spin
model, and as shown in the appendix we can compute properties at equilibrium
by solving $\chi^{*}={\rm argmax}_{\chi}\,\Phi_{\rm RS}(\chi)$ where the RS
free entropy reads

$\displaystyle\Phi_{\rm RS}(\chi)$
$\displaystyle=\frac{\widetilde{\chi}}{2}+\frac{1-\theta}{2}\log\left(\frac{2\pi}{\widetilde{\chi}+1}\right)-\theta\log(\sqrt{2\pi}e)-\frac{\chi^{3}}{2T^{2}}\,,$
$\displaystyle\widetilde{\chi}$ $\displaystyle=\frac{3\chi^{2}}{2T^{2}}\,.$
(25)

The same reasoning can be applied to this phase diagram (and in fact the so-
called decimated versions of message-passing algorithms were proposed as early
as in [54, 55]). Concerning the efficiency of sampling, the same phenomenology
appears again when interpreting this phase diagram. In fact, for the spherical
$p$-spin model, it turns out to be even worse because the temperature where
the difficulty arises for the autoregressive method is $T_{\rm
tri}=\sqrt{1/2}$, which is larger than the one for flow-based and diffusion
models $T_{\rm tri}=2/3$. So in this case, autoregressive-based sampling
algorithms perform worse than flow- or diffusion-based.

###  VI.2 Other models

Fig. 2 then evaluates the phase diagrams of the tilted and pinning measures
for three other models – the Ising $p$-spin (b), the rank-one matrix
estimation (c), and the bicoloring problem on sparse hypergraphs (d). For the
bicoloring problem, defined on random sparse hypergraphs, we can use the
belief propagation equations as Bayesian denoiser [73, 74, 75]. The resulting
equations are quite long, and we defer their presentation in the appendix,
along with their derivation using the cavity method.

We observe the same phenomenology with tri-critical points causing hurdles for
flow-, diffusion- and autoregressive methods reaching out to the phase where
traditional approaches based on MCMC or Langevin work efficiently. In fact, we
expect this picture to always appear for any model with the RFOT phenomenology
where the dynamical temperature is distinct from the ideal glass or Kauzmann
temperature. Such a phenomenology was described in many problems far beyond
those we picked to study, and we can hence anticipate follow-up studies
identifying analogous phase diagrams in many other problems of interest.

Finally, we notice that depending on the model, the position of the tri-
critical point for flow- and diffusion-based methods is better (e.g. for the
spherical $3$-spin, Ising 3-spin) or worse (e.g. for the sparse rank-one
matrix estimation) than for the autoregressive methods. In any case, the
position of the tri-critical point does depend on the noise channel to which
the generative model maps. This leaves open the question of which channel one
should use, for each model, to minimize the range of values for which
generative model-based sampling is suboptimal compared to MCMC techniques that
fail at the dynamical threshold $T_{d}$. It is not inconceivable that it is
possible to reduce $T_{\rm tri}$ very close (or maybe even up to) $T_{d}$ by
optimizing over different distributions in linear interpolant [41], or using
non-linear maps [37]. This is left for future studies.

###  VI.3 Outperforming MCMC in inference models with a hard phase

We now discuss a situation more advantageous for modern techniques. Column (c)
of Fig. 2 depicts the phase diagram of the sparse rank-one matrix estimation
problem that presents an additional interesting feature: here there is a
planted “hidden” signal that we seek to recover. At large values of the noise
the signal is hidden and the RFOT-type phenomenology reappears so that the
high noise behaviour is identical to the high-temperature models of the other
models.

At low noise values, however, there is another phase transition at
$\gamma=\theta=0$ denoted as $\Delta_{\rm alg}$ below which the AMP algorithm
solves the estimation problems optimally and above which it does not up to the
value $\Delta_{\rm IT}$ (the hard region, in red, is delimited by these two
values). Going vertically up in the phase diagram in $\gamma$ or $\theta$ for
$\Delta<\Delta_{\rm alg}$ does not cause any encounter of the hard (red)
phase, and thus sampling based on flow or diffusion or autoregressive networks
works. This has been proven rigorously recently in [92] (with some technical
regularity assumptions on the denoiser [90]).

Yet existing literature collects evidence that in inference problems that
present such a hard phase, local dynamics algorithms such as MCMC and Langevin
are not able to sample efficiently until some yet lower values of noise
$\Delta_{\rm MCMC}$. In particular, this suggestion was put forward indirectly
in [30] by arguing that the metastable phase in the hard region is glassy, and
this glassy nature extends well beyond the region that is hard for message
passing algorithms such as AMP. This was then shown explicitly in follow-up
works starting with an analysis of the dynamics in a mixed spiked matrix-
tensor model in [49], in the phase retrieval problem [79], the planted
coloring [50] and on a rigorous basis in the planted clique problem [80].
Works such as [93] suggest that over-parametrization observed in modern neural
networks mitigates those hurdles, and this may be one of the reasons why over-
parametrization is beneficial.

In light of these works, it is interesting to note that sampling based on
flow, diffusion or autoregressive networks also avoids hurdles stemming from
the glassiness of the hard phase and rather effortlessly so by working in the
space of marginals rather than configurations directly. The phase diagram
presented in Fig. 2 indicates that both diffusion and autoregressive networks
sample the $P_{0}$ efficiently for any $\Delta<\Delta_{\rm alg}$. This poses
an intriguing question for future work of whether over-parametrization of
neural networks that are learning the denoisers from data would still be so
beneficial in these methods.

Finally, we comment on the relevance of these findings beyond the specific
model discussed here, and in particular for the study of physical objects in
finite dimension with short-range interaction. Do we expect a problem embedded
in finite dimension to suffer the same fate as the one discussed here? While
similar phase diagrams as Fig. 2 have been observed in finite dimension in
e.g. [70, 94], the phenomenology of first-order transition is different. From
nucleation arguments, the exponentially hard denoising phase is not expected
to exist in finite dimension [95]. Indeed, it can be proven rigorously that
for any graph that can be embedded in a finite-dimensional lattice, an
efficient algorithm exists [96]. In this case, the analysis of whether a good
denoiser can be learned with a neural network will require a finer study,
depending on the discretization of the ODE close to the transition, and on the
number of points in the data set (perhaps in the vein of [57]). These are
interesting potential new directions of research.

## Conclusion

Our investigation into the efficiency of sampling with modern generative
models, in comparison with traditional methods, reveals distinct strengths and
weaknesses for each. By examining a specific class of probability
distributions from statistical physics, we identified parameters where either
method excels or falls short. Significantly, our approach highlighted
challenges stemming from a first-order discontinuous transition for generative
models-based sampling techniques even in regions of parameters where
traditional samplers work efficiently. While generative models have shown
promise across various applications, it is crucial to understand their
potential pitfalls and advantages in specific contexts, and our paper makes a
key step in this direction.

## Acknowledgments

We acknowledge funding from the Swiss National Science Foundation grant SNFS
OperaGOST and SMArtNet. We also thank Ahmed El Alaoui, Hugo Cui, and Eric
Vanden-Eijnden for enlightening discussions on these problems.

## Appendices

##  Appendix A Computing free entropies: the planting trick

The main technical difficulty is the study of the tilted or pinned measure and
its correlation with an equilibrium configuration. Here, we explain how this
difficulty is avoided. Given a probability measure $P_{0}$, the tilted measure
reads

$P_{\gamma}({\bf x}|{\bf x}_{0},{\bf z})\propto e^{\gamma(t)^{2}\langle{\bf
x},{\bf x}_{0}\rangle+\gamma(t)\langle z,{\bf
x}\rangle-\frac{\gamma(t)^{2}}{2}\|{\bf x}\|^{2}}P_{0}({\bf x})\,.$ (26)

We shall restrict ourselves to “planted” models with a hidden assignment in
the following. While this is the case for actual inference statistical models
(such as the sparse Wigner model we consider), it turns out that the $p$-spin
[82] and the NAESAT model [97] are contiguous to their planted version as long
as $T>T_{K}$ (for $p$-spin) or $\alpha<\alpha_{K}$ for NAESAT. For the
$p$-spin models, such contiguity has also been rigorously established in
certain setups [98, 86, 99].

For all practical purposes, we shall thus study planted models. In this case,
the difficulty is greatly simplified, as the joint distribution of the planted
vector and the disorder equals the joint distribution of the disorder and an
equilibrium configuration. This equivalence, along with the contiguity with
the planted models for non-inference problems, allows us to replace the
equilibrium configuration ${\bf x_{0}}$ in the tilted measure with the planted
vector.

For the case of the Sherrington-Kirkpatrick (SK) model, such an equivalence
was rigorously proven and utilized in [36] using the contiguity between the
planted and unplanted models. Concretely, Proposition 4.2 in [36] shows the
equivalence between the following two methods for generating the tilted
measure:

  1. 1.

Sample ${\bf x_{0}}$ uniformly. Then sample the interaction matrix $J$ as
$J=\frac{\beta}{n}{\bf x_{0}}{\bf x_{0}}^{\top}+W$ where $W\sim\text{GOE}$.

  2. 2.

Sample $J\sim\text{GOE}$. Then sample ${\bf x_{0}}\sim P_{\rm SK}(J)$,

where $P_{\rm SK}(J)$ denotes the Boltzmann Sherrington-Kirkpatrick measure
with interaction matrix $J$. Similarly, based on the contiguity between
planted and unplanted models, we assume that the above equivalence holds for
$T>T_{K}$ (for $p$-spin) and for $\alpha<\alpha_{K}$ (for NAESAT). We note
that for inference problems, i.e. when $J$ is planted for both points (1), (2)
above, the equivalence holds directly based on the definition of the posterior
measure.

Additionally, the tilted measure can be then seen, in terms of free entropy,
as the measure associated now to an inference problem. Consider for
concreteness the planted $p$-spin, also called the spike-tensor model [87,
86]: One extracts a random vector $X^{0}\in\mathbb{R}^{N}$ from a Gaussian or
a Rademacher distribution, and then one aims at recovering $X^{0}$ from the
measure of a) a Gaussian measurement $Y=\alpha X^{0}+\beta Z$ and b) a noisy
tensor measurement
$J_{ijk}=X^{0}_{i}X^{0}_{j}X^{0}_{k}+\frac{1}{T}\xi_{ijk}\forall i<j<k$. The
Bayesian posterior of this model is nothing but the tilted measure (26)
applied to the tensor spike model.

Interestingly, such considerations are not new. Nishimori [63] and later Iba
[64] already used a similar trick, and [29, 78, 85] used it to discuss the
dynamics starting from equilibrium conditions, while [100] used it in the
context of error correction, all for the $p$-spin model. Recently, the same
technique was used to prove the clustering property in the $p$-spin model
[72].

Adding a Gaussian measurement to an inference problem is also a classical
trick used when proving free energies in the mathematical physics literature,
especially in the context of Guerra interpolation [101, 102] for Bayes optimal
models, see e.g. [66, 103, 32, 67], and thus such free energies have been
solved rigorously as well.

How does this change for the decimated problem? In this case, the additional
Gaussian channel is replaced by an erasure channel. This is precisely one of
the alternatives used as well in the mathematical physics literature! Indeed,
the pinning lemma [68] (see also [69]) is often used instead of the Gaussian
channel.

These considerations are really helpful, as we can now solve these problems as
a simple variant of problems already solved in the literature, often
rigorously (in the case of the $p$-spin we refer to [83] and [49], and for the
spike Wigner model to [81, 84, 104].

For the sparse case, a rigorous control is harder, and thus we shall simply
stay at the level of rigour of the cavity method [20] and use the results of
[74]. Note however that the method is trustworthy [97].

Finally, since these are variations of known inference problems, we can
leverage on the existing work on approximate message passing [31] for the
$p$-spin model [87, 83] and the spike model [81], see in particular [86] for a
detailed presentation.

##  Appendix B Asymptotic solutions and Phase Diagrams

In this section, we present how each of the phase diagrams we presented in the
main text can be produced. Specifically, we first remind our definition of
tilted and pinning measures and of the order parameters we are going to use
for the analysis. After defining the probability distribution associated to
each problem, we give the expression of the RS free entropy, from which one
can compute the values of the parameters of the problem at which the spinodal
points are located, i.e. the points at which the potential develops a second
maxima by continuous deformation, and also the point corresponding to the IT
transition, i.e. where the two maxima exchange the role of global and local
maxima.

We then report the expression of the denoisers (AMP/BP) used in the sampling
schemes illustrated in the main text, with their associated self-consistent
asymptotic equations (state evolution for AMP and Cavity equations for BP).
Looking at the fixed points of these equations, starting from an uninformed
and informed initialization, we can plot the difference between the values of
the order parameter reached at the fixed point in these two cases, which
allows detecting the phases in which multiple fixed points are present.

In Eq. (26) we have recalled the expression for the tilted measure
characterizing the flow-based sampling scheme. As we already mentioned,
$\gamma$ is nothing but a rescaled sampling time, such that studying the
properties of the tilted measure varying $\gamma$ allows us to characterize
the properties of the Bayesian denoising problem at all times during sampling.

In the same way, let us remind the pinning measure, which we use to analyse
the autoregressive-based sampling procedure:

$P_{\theta}({\bf x}|{\bf x}_{0},S_{\theta})\propto P_{0}({\bf x})\prod_{i\in
S_{\theta}}\delta([{x}_{i}-[{\bf x}_{0}]_{i})$ (27)

where $\theta$ is the fraction of pinned variables.

Finally, let us remind that we shall study the evolution in time (or
equivalently in $\gamma\in[0,\infty[$ (AWGN) and $\theta\in[0,1]$ (BEC)) of
the following order parameters:

$\displaystyle\mu(\gamma)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\mathbb{E}[\hat{\bf x}(\gamma)\cdot{\bf x}_{0}]\,,$
(28) $\displaystyle\chi(\gamma)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\mathbb{E}[\lVert\hat{\bf x}(\gamma)\rVert^{2}]\,,$
(29)

and analogously we can define $\mu(\theta)$ and $\chi(\theta)$. Concretely, we
will consider only cases in which the Nishimori identities hold, such that
these two quantities always coincide, and thus we will restrict our analysis
to $\chi$.

Let us now go through each one of the models mentioned in the main.

###  B.1 Sparse rank-one matrix factorization

We consider the Bayes-Optimal rank-one matrix estimation (or rank-one matrix
factorization) problem:

Given a hidden vector ${\bf x^{*}}$, sampled from the so-called Rademacher-
Bernoulli prior distribution

$P_{X}(x)=(1-\rho)\delta_{x,0}+\frac{\rho}{2}\left(\delta_{x,+1}+\delta_{x,-1}\right)\,,\quad
x^{*}_{i}\sim P_{X}\;\forall\,i\,.$ (30)

one has access to noisy observations, that is a matrix $J_{ij}$ is composed by
a rank-one spike plus i.i.d. Gaussian noise:

$\displaystyle
J_{ij}=\frac{x^{*}_{i}x^{*}_{j}}{\sqrt{N}}+\widetilde{z}_{ij}\,,\quad\widetilde{z}_{ij}=\widetilde{z}_{ji}\sim\mathcal{N}(0,\Delta)\,;$

and the goal is to infer ${\bf x^{*}}$ in the best way possible. There are
many important problems in statistics and machine learning that can be
expressed in this way [105], and this model has been the subject of many works
both from the statistics [106, 107, 81, 108, 104] and the statistical physics
communities [109, 83]. The presentation of this problem closely follows the
study presented in [83].

From the Bayesian point of view, the problem amounts to sampling from the
posterior. One way to introduce the model is through the following probability
distribution:

$P_{0}({\bf
x)}\propto\prod_{i}P_{X}(x_{i})\prod_{i<j}\exp\left(\frac{1}{\Delta\sqrt{N}}J_{ij}x_{i}x_{j}-\frac{1}{2\Delta
N}x_{i}^{2}x_{j}^{2}\right)\,.$ (31)

With this distribution, ${\bf x}$ will be a random vector with, on average, a
fraction $\rho$ of components that are Ising spin variables (i.e. each $x_{i}$
takes values $\pm 1$) and the rest of the entries that are put to zero, in
such a way that the parameter $\rho$ controls the sparsity of the vector we
want to retrieve.

As discussed in Appendix A, the tilted measure with an equilibrium vector
${\bf x}_{0}$ is equivalent to the one where the vector ${\bf x}_{0}$ is
planted. Therefore, in what follows, we shall assume that ${\bf x}_{0}$
corresponds to the planted configuration. We thus obtain the following tilted
measure for diffusion and flow-based models:

$\displaystyle P_{\gamma}({\bf x)}$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right)e^{\gamma(t)^{2}\langle{\bf
x},{\bf x}_{0}\rangle+\gamma(t)\langle{\bf z},{\bf
x}\rangle-\frac{\gamma(t)^{2}}{2}\|{\bf
x}\||^{2}}\left(\prod_{i<j}e^{\frac{1}{\Delta\sqrt{N}}\tilde{z}_{ij}x_{i}x_{j}+\frac{1}{\Delta{\mathbb{N}}}x_{i}x_{j}[{\bf
x}_{0}]_{i}[{\bf x}_{0}]_{j}-\frac{1}{2\Delta N}x_{i}^{2}x_{j}^{2}}\right)$
(32)

With respect to the original model, the tilted measure simply includes in
addition, a field in the planted direction, a random field and a
renormalization of the constant in front of the quadratic part. As mentioned
in Appendix A, this corresponds to an equivalent inference problem with an
additional Gaussian measurement.

The pinned measure is slightly different. In this case, it modifies the
original problems as (denoting the pinned list as $S_{\theta}$):

$\displaystyle P_{\theta}({\bf x)}$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{\theta}}\left(\prod_{i\notin
S_{\theta}}P_{X}(x_{i})\right)\left(\prod_{i\in
S_{\theta}}\delta(x_{i}-x_{i}^{*})\right)\left(\prod_{i<j}e^{\frac{1}{\Delta\sqrt{N}}\tilde{z}_{ij}x_{i}x_{j}+\frac{1}{\Delta{\mathbb{N}}}x_{i}x_{j}[{\bf
x}_{0}]_{i}[{\bf x}_{0}]_{j}-\frac{1}{2\Delta N}x_{i}^{2}x_{j}^{2}}\right)$
(33)

####  B.1.1 Replica free entropy

The replica formula for such problems can be found in many places and we refer
to the mathematical literature for the detailed rigorous statements: [81, 108,
84, 104, 103, 110]. In particular, [103] gives a generic proof using the
adaptive interpolation methods in the presence of a Gaussian channel is added,
which turns out to give the same measure as the tilted one, while [104] uses
instead a pinned measure. In both case, we can thus adapt the results in the
literature: The asymptotic free entropy is given by the maximum of the so-
called replica symmetric potential[20] :

$\displaystyle\frac{1}{N}{\mathbb{E}}_{{\bf x}_{0},{\bf z},\tilde{\bf z}}\log
Z_{\gamma}$ $\displaystyle\xrightarrow[N\to\infty]{}$ $\displaystyle{\rm
argmax}\;\Phi_{RS}(m)$ (34) $\displaystyle\Phi_{RS}(m)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{w,x_{0}}\left[\log
Z_{x}\left(\frac{m}{\Delta},\frac{m}{\Delta}x_{0}+\sqrt{\frac{m}{\Delta}}w\right)\right]-\frac{m^{2}}{4\Delta}$
(35)

where $m$ is the order parameter of the problem, $x_{0}\sim P_{X}$,
$w\sim\mathcal{N}(0,1)$ while $Z_{x}$ depends on the specific measure
considered. The same is valid for the pinned measure, as long as one
substitutes $Z_{\gamma}$ with $Z_{\theta}$.

##### Tilted measure:

In the case of the tilted measure (26), we have

$\begin{split}Z_{x}(A,B;x_{0})&=\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)\\\ &=\rho
e^{-(A+\gamma^{2})/2}\cosh(B+\gamma w+\gamma^{2}x_{0})+(1-\rho)\end{split}$
(36)

that we can put into Eq. (35) to get

$\displaystyle\Phi_{\rm RS}(\chi)=$
$\displaystyle\rho{\mathbb{E}}_{w}\left[\log\left((1-\rho)+\rho
e^{-\widetilde{\chi}/2}\cosh\left(\widetilde{\chi}+\sqrt{\widetilde{\chi}}w\right)\right)\right]$
(37) $\displaystyle+(1-\rho){\mathbb{E}}_{w}\left[\log\left((1-\rho)+\rho
e^{-\widetilde{\chi}/2}\cosh\left(\sqrt{\widetilde{\chi}}w\right)\right)\right]-\frac{\chi^{2}}{4\Delta}\,,\quad\widetilde{\chi}=\frac{\chi}{\Delta}+\gamma^{2}\,.$

##### Pinning measure:

Meanwhile, considering the pinning measure (27) leads to

$\begin{split}Z_{x}(A,B;x_{0})&=\begin{cases}\int\mathop{}\\!\mathrm{d}xP_{X}(x)\delta_{x,x_{0}}\exp(Bx-
Ax^{2}/2)&\text{with probability }\theta\\\
\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)&\text{with probability
}1-\theta\end{cases}\\\
&=\begin{cases}P_{X}(x_{0})\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)&\text{with
probability }\theta\\\ \rho\exp(-A/2)\cosh(B)+(1-\rho)&\text{with probability
}1-\theta\end{cases}\end{split}$ (38)

which in turn gives

$\displaystyle\Phi_{\rm RS}(\chi)=$
$\displaystyle\theta\left(\rho\log\rho+(1-\rho)\log(1-\rho)+\frac{\rho}{2}\widetilde{\chi}\right)$
(39)
$\displaystyle+(1-\theta)\bigg{[}\rho{\mathbb{E}}_{w}\left[\log\left((1-\rho)+\rho
e^{-\widetilde{\chi}/2}\cosh\left(\widetilde{\chi}+\sqrt{\widetilde{\chi}}w\right)\right)\right]$
$\displaystyle+(1-\rho){\mathbb{E}}_{w}\left[\log\left((1-\rho)+\rho
e^{-\widetilde{\chi}/2}\cosh\left(\sqrt{\widetilde{\chi}}w\right)\right)\right]\bigg{]}-\frac{\chi^{2}}{4\Delta}\,,\quad\widetilde{\chi}=\frac{\chi}{\Delta}\,.$

####  B.1.2 Message-Passing algorithm

The derivation of the AMP algorithm for this problem has a long history, and
is connected to the Thouless-Anderson-Palmer (TAP) equations [51]. For this
problem, the introduction of TAP as an iterative algorithm is due to
Bolthausen [111] and has been adapted to the present situation in [107, 81].

##### Tilted measure:

The equivalence between the tilted and the planted measure with external field
allows us to reduce the AMP iterations for the tilted measure to the ones for
an associated inference problem. In [83], the authors provided a framework for
deriving the AMP iterates for such an inference problem involving pair-wise
interactions between spins. While the tilted measure defined by Eq. (32)
involves additional random and planted fields, the generality of the
derivation in [83] allows them to be straightforwardly incorporated into the
single-site factors $P_{X}(x_{i})$ in [83]. Through an adaptation of the
derivation of Equations 66, 67 in [83], we obtain:

$\begin{cases}\widehat{x}_{i}^{t+1}=\frac{\rho\tanh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i})}{\rho+\frac{(1-\rho)\exp((A^{t}+\gamma^{2})/2)}{\cosh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i})}},\quad\sigma_{i}^{t+1}=\rho\frac{\rho+(1-\rho)e^{(A^{t}+\gamma^{2})/2}\cosh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i})}{\left(\rho\cosh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i})+(1-\rho)\exp((A+\gamma^{2})/2)\right)^{2}}\\\
A^{t}=\frac{\|{\bf\widehat{x}^{t}}\|_{2}^{2}}{\Delta N}\,;\quad
B_{i}^{t}=\frac{1}{\Delta\sqrt{N}}{\bf
J_{i}}\cdot{\bf\widehat{x}^{t}}-\frac{1}{N\Delta}\widehat{x}_{i}^{t-1}\sum_{k}\sigma_{k}^{t}\end{cases}$
(40)

where $\alpha(t)$ and $\beta(t)$ are the functions defining the interpolant
process, fixed at the start, and ${\bf Y}_{t}$ is the value of the noisy
observation at time $t$.

##### Pinning measure:

When considering the pinning measure (27), the AMP equations are only a slight
variation of the ones presented for the flow-based case.

Specifically, in autoregressive-based sampling we choose a fraction $\theta$
of the variables, for which we fix $\widehat{x}_{i}^{t}=[{\bf
x}_{0}]_{i},\;\sigma_{i}^{t}=0$; this is due to the fact that their posterior
means are completely polarized on the solution. For the rest of the variables,
a fraction $1-\theta$, the AMP equations are exactly the ones for diffusion in
Eq. (40), at $\gamma=0$.

The resulting algorithm is thus

$\begin{cases}&\widehat{x}_{i}^{t+1}=\begin{cases}[{\bf x}_{0}]_{i}&\text{if
}i\in{S_{\theta}}\\\
\frac{\rho\tanh(B_{i}^{t})}{\rho+\frac{(1-\rho)\exp(A^{t}/2)}{\cosh(B_{i}^{t})}},&\text{otherwise}\end{cases}\,,\quad\sigma_{i}^{t+1}=\begin{cases}0&\text{if
}i\in{S_{\theta}}\\\
\rho\frac{\rho+(1-\rho)e^{A^{t}/2}\cosh(B_{i}^{t})}{\left(\rho\cosh(B_{i}^{t})+(1-\rho)\exp(A/2)\right)^{2}}&\text{otherwise}\\\
\end{cases}\\\ &A^{t}=\frac{\|{\bf\widehat{x}^{t}}\|_{2}^{2}}{\Delta
N}\,;\quad B_{i}^{t}=\frac{1}{\Delta\sqrt{N}}{\bf
J_{i}}\cdot{\bf\widehat{x}^{t}}-\frac{1}{N\Delta}\widehat{x}_{i}^{t-1}\sum_{k}\sigma_{k}^{t}\end{cases}$
(41)

####  B.1.3 State Evolution equations

The advantage of AMP is that it can be rigorously tracked by the State
Evolution equations [111, 107, 31] that turn out to be nothing but the fixed
point equations of the associated replica free entropy. Again, we can use the
generic results reported in [83] to get

$\displaystyle m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\rm
in}\left(\frac{m^{t}}{\Delta},\frac{m^{t}}{\Delta}x_{0}+\sqrt{\frac{m^{t}}{\Delta}}w\right)x_{0}\right]$
(42)

where $x_{0}\sim P_{X}$, $w\sim\mathcal{N}(0,1)$ and $f_{\rm in}$ is the input
channel and depends on the specific problem.

##### Tilted measure:

For the tilted measure (26) we get

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)}\\\ &=\frac{\rho\tanh(B+\gamma
w+\gamma^{2}x_{0})}{\rho+\frac{(1-\rho)\exp\left((A+\gamma^{2})/2\right)}{\cosh(B+\gamma
w+\gamma^{2}x_{0})}}\end{split}$ (43)

which leads to the state evolution equations:

$\chi^{t+1}=\rho^{2}{\mathbb{E}}_{w}\left[\frac{\tanh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)}{\rho+\frac{(1-\rho)\exp(\widetilde{\chi}^{t}/2)}{\cosh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)}}\right]\,,\quad\widetilde{\chi}^{t}\equiv\frac{\chi^{t}}{\Delta}+\gamma^{2}\,,\quad
w\sim\mathcal{N}(0,1)$ (44)

##### Pinning measure:

In the same way, the pinning measure (27) is associated to

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\begin{cases}\frac{P_{X}(x_{0})x_{0}\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}{P_{X}(x_{0})\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}&\text{with
probability }\theta\\\ \frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)}&\text{with
probability }1-\theta\end{cases}\\\ &=\begin{cases}x_{0}&\text{with
probability }\theta\\\
\rho\tanh(B)/(\rho+(1-\rho)\exp\left(A/2\right)/\cosh(B))&\text{with
probability }1-\theta\end{cases}\end{split}$ (45)

and consequently to the fixed point equations

$\frac{\chi^{t+1}}{\rho}=\theta+(1-\theta){\mathbb{E}}_{z}\left[\frac{\rho\tanh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)}{\rho+\frac{(1-\rho)\exp(\widetilde{\chi}^{t}/2)}{\cosh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)}}\right]\,,\quad\widetilde{\chi}^{t}\equiv\frac{\chi^{t}}{\Delta}\,,\quad
w\sim\mathcal{N}(0,1)\,.$ (46)

####  B.1.4 Phase diagrams

Figure 3: Phase diagrams for flow-based sampling (left) and autoregressive-
based sampling (right) for the sparse rank-one model, with Rademacher-
Bernoulli prior and sparsity $\rho=0.08$. On the x-axis we put the rescaled
signal-to-noise-ratio $\Delta/\rho^{2}$ and on the y-axis the ratio
$\gamma^{2}=\alpha^{2}/\beta^{2}$ (left) and the decimated ratio $\theta$
(right). We compute the order parameter $\chi/\rho$, defined in (29), both
from an uninformed and an informed initialization, and we plot the difference
between the two. The dashed white lines are the spinodal lines, while the
dashed black one is the IT threshold. Note that for the flow-based case (left
panel), we show explicitly in Fig. 4 and Fig. 5 the behaviour of the free
entropy functional for $\Delta/\rho^{2}=1.05$ and $\Delta/\rho^{2}=0.98$. Here
in both plots (left and right) we have that the dynamical transition is at
$\Delta_{d}/\rho^{2}\approx 1.041$, the IT/Kauzmann transition is at
$\Delta_{IT}/\rho^{2}\approx 1.029$, while the tri-critical points are at
$\Delta_{\rm tri}/\rho^{2}\approx 1.08$ for flow-based and $\Delta_{\rm
tri}/\rho^{2}\approx 1.069$ for autoregressive based sampling.

In Fig. 3 we present the phase diagrams for the sparse rank-one matrix
factorization problem, choosing $\rho=0.08$ as value for sparsity. We remind
that $\rho$ must be small enough to observe a first-order phenomenology [108].

The section of the parameter space displayed is the same as in the plots
presented in the main text, but here we display directly the difference
$\chi_{\rm inf}-\chi_{\rm uninf}$, so that the coloured zones of the plots are
the ones displaying multiple fixed points, as opposed to the black ones. We
furthermore draw as white dashed lines the spinodal points, and as a black
dashed line the IT threshold, both defined at the beginning of Appendix B. For
the flow-based plot, we also show explicitly in Fig. 4 and Fig. 5 the
behaviour of the free entropy functional for $\Delta/\rho^{2}=1.05$ and
$\Delta/\rho^{2}=0.98$. For Fig. 4 the first order transition is apparent,
while for Fig. 5 the transition is found to be continuous.

Spike Wigner model: $\Phi_{\rm RS}(\chi)$ for
$\Delta_{d}/\rho^{2}<\Delta/\rho^{2}=1.05<{\Delta}_{\rm tri}/\rho^{2}$

(a) $\gamma^{2}=0$ (b) $\gamma^{2}=0.03$ (c) $\gamma^{2}=0.05$ (d)
$\gamma^{2}=0.08$

Figure 4: The free entropy function $\Phi_{\rm RS}(\chi)$ in the region where
the flow-based model fails because of the jump around panel(c).

Spike Wigner model: $\Phi_{\rm RS}(\chi)$ for
$\Delta/\rho^{2}=0.98<\Delta_{\rm alg}/\rho^{2}$

(a) $\gamma^{2}=0$ (b) $\gamma^{2}=0.02$ (c) $\gamma^{2}=0.05$ (d)
$\gamma^{2}=0.1$

Figure 5: The free entropy function $\Phi_{\rm RS}(\chi)$ in the region where
the flow-based model succeeds as the position of the maxima is always unique.

##### Autoregressive networks vs Flows

As we can see from the plots, for this model the tri-critical point for flow-
based sampling is at $\Delta_{\rm tri}/\rho^{2}\approx 1.08$, while for
autoregressive-based sampling is at $\Delta_{\rm tri}/\rho^{2}\approx 1.069$,
meaning that the gap with MCMC and Langevin sampling is smaller in this latter
case. In other words, this means that there is a range of values of inverse-
SNR $\Delta$ for which autoregressive networks based sampling (along with MCMC
and Langevin) is efficient, while flow-based sampling is not.

This appears not to be specific to this particular value of $\rho$, but in all
the range of sparsity in which the model presents a first order phase
transition, autoregressive-based sampling appears to be more efficient than
flow-based sampling, in the sense explained above, as shown in Fig. 6.

Figure 6: Flows vs autoregressive in the sparse rank-one problem. We plot the
values of the tri-critical point for flow-based (orange dots) and
autoregressive-based (blue crosses) sampling on the sparse rank-one model,
when varying the sparsity parameter $\rho$. Specifically, we plot $\Delta_{\rm
tri}/\rho$ for the two cases, comparing them also to the Dynamical transition
value $\Delta_{d}/\rho$ (red line) and the IT transition value $\Delta_{\rm
IT}/\rho$ (green line). Finally, we also put a black cross at the point
$(\rho_{\rm max},\Delta_{\rm max}/\rho_{\rm max})$, taken from [83], which
corresponds to the maximum value of $\rho$ at which there is a first order
phase transition.

##### The easy phase (and a subtle point)

For the expert reader, it is worth pointing a subtle difference between what
happens at $\Delta=\rho^{2}$ and $\Delta=\Delta_{\rm alg}<\rho^{2}$ (see also
[83]). Indeed, we use here the definition of the easy phase as the absence of
a metastable maxima that is trapping the dynamics. In other words, the state
evolution of AMP initialized at the uninformed fixed point finds the global
maxima. This is indeed what is going on for $\Delta=\Delta_{\rm alg}$, as
illustrated in Figure 5.

However, it is worth mentioning that already below $\Delta=\rho^{2}$, the
fixed point at zero is unstable, so that there are two fixed point: the
”correct” one at large $\chi$, and the one found by AMP at low but non-zero
$\chi$. This phenomenon, usually called the Baik-Ben Arous-Peche (BBP)
transition, is illustrated in Figure 7. In this phase, while AMP, for instance
(but also a standard spectral method such as PCA [112, 83]) would find an
estimator correlated with the ground truth (and thus $\chi^{*}>0$ even at
$\gamma=0$), there is still a first-order phase transition and AMP would not
be optimal. Since there is a discontinuity, the flow-based method suffers from
the same problem, and fails to sample as well.

The BBP transition at $\gamma^{2}=0$: $\Phi_{\rm RS}(\chi)$ for
$\Delta/\rho^{2}=0.99$ (left) and $\Delta/\rho^{2}=1.01$ (right)

(a) $\Delta/\rho^{2}=0.99$ (zoom) (b) $\Delta/\rho^{2}=0.99$ (unzoom) (c)
$\Delta/\rho^{2}=1.01$ (zoom) (d) $\Delta/\rho^{2}=1.01$ (unzoom)

Figure 7: Illustration of the BBP[112] transition around $\Delta=\rho^{2}$
where the spurious minima stop to be at zero (see the difference between the
zoomed (a) and (c)), but there the correct non-spurious minima is still at
larger value of $\chi$ (see the unzoomed (b) and (d)). This is still a hard
phase for inference and sampling, and the problem remains so until
$\Delta<\Delta_{\rm alg}$

###  B.2 Ising $p$-spin model

We consider now the $p$-spin model (here with $p=3$), which is one of the most
important statistical physics problems in spin glass theory. First, let us
look at the Ising version, defined by the following Hamiltonian:

$\displaystyle{\cal H}({\bf
x})=-\frac{\sqrt{3}}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}\,,$ (47)

where $J_{ijk}\sim\mathcal{N}(0,1)$ and the variables $x_{i}$ are constrained
to be $\pm 1$. This version of the model was introduced in [113], and solved
with the replica method in [19]. It is hard to underestimate its importance in
the spin glass and glass theory as the prototype of the mean-field random
first-order model [114, 115, 70, 46, 116].

For temperatures higher than the spin glass, or Kauzmann, temperature, the
model can be proven [98, 86, 99] to be contiguous to its “planted version”,
the tensor factorization problem [87]. This is nothing but the tensor
generalization of the former matrix estimation problem of the preceding
section. We shall report the derivations presented in [86], which considers
the planted version of the problem, but that (thanks to contiguity) also
describe the unplanted model, and the mapping between the two can be done
using $\Delta=\frac{2}{3}T^{2}$.

As before, from a Bayesian perspective, the problem boils down to sampling the
following posterior

$P_{0}({\bf
x})\propto\prod_{i}P_{X}(x_{i})\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}\,,$
(48)

where the prior $P_{X}(x)=\delta_{x,-1}/2+\delta_{x,+1}/2$ constrains the
variables to be $\pm 1$.

As we have reminded for the previous model, when considering the tilted
measure (26) we will exploit again the fact that considering an equilibrium
configuration ${\bf x}_{0}$ is equivalent to take as ${\bf x}_{0}$ the planted
configuration, as long as we are beyond the Spin Glass temperature.

We shall thus consider the following tilted measure for diffusion and flow-
based sampling:

$P_{\gamma}({\bf
x})=\frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right)e^{\gamma(t)^{2}\langle{\bf
x},{\bf x}_{0}\rangle+\gamma(t)\langle z,{\bf
x}\rangle-\frac{\gamma(t)^{2}}{2}\|{\bf
x}\||^{2}}\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}$ (49)

As before, we can notice that with respect to the original measure in (48),
the tilting adds a field in the planted direction, a random field and a
normalization factor depending on the l2 norm.

For the pinned measure defined in (27) the original problem becomes (denoting
as $S_{\theta}$ the pinned list):

$P_{\theta}({\bf x})=\frac{1}{Z_{\theta}}\left(\prod_{i\notin
S_{\theta}}P_{X}(x_{i})\right)\left(\prod_{i\in
S_{\theta}}\delta(x_{i}-x_{i}^{*})\right)\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}$
(50)

####  B.2.1 Replica free entropy

The replica formula for this model has been studied extensively in the
literature, and can also be proven rigorously, see [86]. Again, the asymptotic
free entropy is given by the maximum of the so-called replica symmetric
potential [20] :

$\displaystyle\frac{1}{N}{\mathbb{E}}_{{\bf x}_{0},{\bf z},J}\log Z_{\gamma}$
$\displaystyle\xrightarrow[N\to\infty]{}$ $\displaystyle{\rm
argmax}\;\Phi_{RS}(m)$ (51) $\displaystyle\Phi_{RS}(m)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{w,x_{0}}\left[\log
Z_{x}\left(\frac{m^{2}}{\Delta},\frac{m^{2}}{\Delta}x_{0}+\sqrt{\frac{m^{2}}{\Delta}}w\right)\right]-\frac{m^{3}}{3\Delta}$
(52)

where $m$ is the order parameter of the problem, $x_{0}\sim P_{X}$,
$w\sim\mathcal{N}(0,1)$ while $Z_{x}$ depends on the specific measure
considered. The same is valid for the pinned measure, as long as one
substitutes $Z_{\gamma}$ with $Z_{\theta}$.

In the following, since we will be interested on the unplanted model, we will
use the temperature $T=\sqrt{3\Delta/2}$ as signal-to-noise parameter.

##### Tilted measure:

In the case of the tilted measure (26) we have

$\begin{split}Z_{x}(A,B;x_{0})&=\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)\\\
&=e^{-(A+\gamma^{2})/2}\cosh(B+\gamma w+\gamma^{2}x_{0})\end{split}$ (53)

that can be put into Eq. (51) to get

$\Phi_{\rm
RS}(\chi)=-\frac{\widetilde{\chi}}{2}+{\mathbb{E}}_{w}\left[\log\cosh\left(\widetilde{\chi}+\sqrt{\widetilde{\chi}}w\right)\right]-\frac{\chi^{3}}{2T^{2}}\,,\quad\widetilde{\chi}=\frac{3\chi^{2}}{2T^{2}}+\gamma^{2}\,.$
(54)

##### Pinning measure:

Meanwhile, considering the pinning measure (27) leads to

$\begin{split}Z_{x}(A,B;x_{0})&=\begin{cases}\int\mathop{}\\!\mathrm{d}xP_{X}(x)\delta_{x,x_{0}}\exp(Bx-
Ax^{2}/2)&\text{with probability }\theta\\\
\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)&\text{with probability
}1-\theta\end{cases}\\\
&=\begin{cases}P_{X}(x_{0})\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)&\text{with
probability }\theta\\\ \exp(-A/2)\cosh(B)&\text{with probability
}1-\theta\end{cases}\end{split}$ (55)

which in turn gives

$\displaystyle\Phi_{\rm
RS}(\chi)=\frac{2\theta-1}{2}\widetilde{\chi}+(1-\theta){\mathbb{E}}_{w}\left[\log\cosh\left(\widetilde{\chi}+\sqrt{\widetilde{\chi}}w\right)\right]-\frac{\chi^{3}}{2T^{2}}\,,\quad\widetilde{\chi}=\frac{3\chi^{2}}{2T^{2}}\,.$
(56)

####  B.2.2 Message-passing algorithm

As for the SK model and the sparse rank-one matrix factorization problem, the
AMP iterates and the associated Thouless-Anderson-Palmer (TAP) equations [51]
have been widely studied for the $p$-spin models [117]. Here we consider the
formalism of [86], which presents the AMP iterates for the Spike Tensor model,
and we use contiguity to derive equations valid for the unplanted model with
the tilting (or pinning) field. We remind again that the mapping between the
two models is given by $\Delta=\frac{2T^{2}}{3}$.

##### Tilted measure:

Due to contiguity with the planted model and the equivalence of the tilted
measure and an associated measure with a planted field described in Section A,
the AMP iterations for the tilted measure can be obtained through the ones for
the associated inference problem described in [83].This yields to

$\begin{cases}\widehat{x}_{i}^{t+1}=\tanh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i}),\quad\sigma_{i}^{t+1}=\cosh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i})^{-2}\\\
B_{i}^{t}=\frac{\sqrt{3}\beta}{N}\sum_{j<k}J_{ijk}\widehat{x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\widehat{x}_{i}^{t-1}{\bf\widehat{x}}^{t}\cdot{\bf\widehat{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\\\
\end{cases}$ (57)

where $\alpha(t)$ and $\beta(t)$ are the functions defining the interpolant
process, fixed at the start, and ${\bf Y}_{t}$ is the value of the noisy
observation at time $t$.

##### Pinning measure:

The AMP equations for the pinning measure (27) are a slight variation of the
ones just presented for the tilted measure.

Specifically, in the autoregressive-based sampling scheme for a fraction
$\theta$ of the variables we fix
$\widehat{x}_{i}^{t}=x_{0},\;\sigma_{i}^{t}=0$, since their posterior means
are totally polarized on the solution. For the rest of the variables, a
fraction $1-\theta$, the AMP equations are the same as the one presented for
diffusion in Eq. (57), provided that we fix $\gamma=0$.

The resulting equations are:

$\begin{cases}&\widehat{x}_{i}^{t+1}=\begin{cases}[{\bf x}_{0}]_{i}&\text{if
}i\in{S_{\theta}}\\\
\tanh(B_{i}^{t}),&\text{otherwise}\end{cases}\,,\quad\sigma_{i}^{t+1}=\begin{cases}0&\text{if
}i\in{S_{\theta}}\\\ \cosh(B_{i}^{t})^{-2}&\text{otherwise}\\\ \end{cases}\\\
&B_{i}^{t}=\frac{\sqrt{3}\beta}{N}\sum_{j<k}J_{ijk}\widehat{x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\widehat{x}_{i}^{t-1}{\bf\widehat{x}}^{t}\cdot{\bf\widehat{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\end{cases}$
(58)

####  B.2.3 State evolution equations

The advantage of AMP is that it can be rigorously tracked by the State
Evolution equations [31], that can be proven to be nothing but the fixed point
of the replica potential (51). We can use the formalism in [86] to get

$\displaystyle m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\rm
in}\left(\frac{(m^{t})^{2}}{\Delta},\frac{(m^{t})^{2}}{\Delta}x_{0}+\sqrt{\frac{(m^{t})^{2}}{\Delta}}w\right)x_{0}\right]$
(59)

where $x_{0}\sim P_{X}$, $w\sim\mathcal{N}(0,1)$ and $f_{in}$ is the input
channel and depends on the specific problem. Again, we will state our results
using the temperature $T=\sqrt{3\Delta/2}$.

##### Tilted measure:

For the tilted measure (26) we have

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)}\\\ &=\tanh(B+\gamma
w+\gamma^{2}x_{0})\end{split}$ (60)

which leads to the State Evolution equations

$\chi^{t+1}={\mathbb{E}}_{w}\left[\tanh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)\right]\,,\quad\widetilde{\chi}^{t}\equiv\frac{3(\chi^{t})^{2}}{2T^{2}}+\gamma^{2}\,,\quad
w\sim\mathcal{N}(0,1)$ (61)

##### Pinning measure:

In the same way, for the pinning measure (27) we get

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\begin{cases}\frac{P_{X}(x_{0})x_{0}\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}{P_{X}(x_{0})\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}&\text{with
probability }\theta\\\ \frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)}&\text{with
probability }1-\theta\end{cases}\\\ &=\begin{cases}x_{0}&\text{with
probability }\theta\\\ \tanh(B)&\text{with probability
}1-\theta\end{cases}\end{split}$ (62)

and thus the fixed point equations

$\chi^{t+1}=\theta+(1-\theta){\mathbb{E}}_{w}\left[\tanh(\widetilde{\chi}^{t}+\sqrt{\widetilde{\chi}^{t}}w)\right]\,,\quad\widetilde{\chi}^{t}\equiv\frac{3(\chi^{t})^{2}}{2T^{2}}\,,\quad
w\sim\mathcal{N}(0,1)$ (63)

####  B.2.4 Phase diagrams

The phase diagrams for the Ising $p$-spin are presented in Figure 8, in the
same style as the previous section. In this case, notably, we observe that the
flow-based method is advantageous with respect to the autoregressive one, in
the sense that the former shows a smaller gap $T_{\rm tri}-T_{d}$ compared to
the latter. This is the opposite situation compared to the previous problem.

Figure 8: Phase diagrams for flow-based sampling (left) and autoregressive-
based sampling (right) for the Ising p-spin model. On the x-axis we put the
temperature $T$ and on the y-axis the ratio $\gamma^{2}=\alpha^{2}/\beta^{2}$
(left) and the decimated ratio $\theta$ (right). We compute the order
parameter $\chi$, defined in (29), both from an uninformed and an informed
initialization, and we plot the difference between the two. The dashed white
lines are the spinodal lines, while the dashed black one is the IT threshold,
both defined at the beginning of the section. In both plots we have that the
dynamical transition is at $T_{d}\approx 0.682$, the Kauzmann transition is at
$T_{K}\approx 0.652$, while the tri-critical points are at $T_{\rm tri}\approx
0.741$ for flow-based and $T_{\rm tri}\approx 0.759$ for autoregressive-based
sampling.

###  B.3 Spherical p-spin model

The spherical $p$-spin model is a variation of the Ising model where the
configurations are constrained to lie on the sphere ${\bf x}\in{\cal S}^{N-1}$
[117]. Again, we shall use the planted version, which is asymptotically
equivalent at high temperature with the unplanted problem having independent
Gaussian entries. The Hamiltonian thus still reads:

$\displaystyle{\cal H}({\bf
x})=-\frac{\sqrt{3}}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}\,,$ (64)

where $J_{ijk}\sim\mathcal{N}(0,1)$ and now $P_{X}(x_{i})={\cal N}(0,1)$.

Same as its Ising version, for temperatures higher than the Kauzmann
temperature the model can be proven [86] to be contiguous to its “planted
version”, the tensor factorization problem [87]. Here we report the
derivations presented in [86], and the mapping between the two is again
$\Delta=\frac{2}{3}T^{2}$.

From a Bayesian point of view, the posterior measure we want to sample is the
following:

$P_{0}({\bf
x})\propto\prod_{i}P_{X}(x_{i})\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}\,,$
(65)

where $P_{X}(x)={\cal N}(0,1)$ in high-dimension constrains the variables to
be on the sphere.

We first consider the following tilted measure, arising in diffusion and flow-
based sampling:

$P_{\gamma}({\bf
x})=\frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right)e^{\gamma(t)^{2}\langle{\bf
x},{\bf x}_{0}\rangle+\gamma(t)\langle z,{\bf
x}\rangle-\frac{\gamma(t)^{2}}{2}\|{\bf
x}\||^{2}}\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}$ (66)

Again, with respect to the original measure in Eq. (65), this measure presents
an additional field in the planted direction, a random field and a constant
factor depending on the l2 norm.

For the pinned measure defined in (27) the original problem becomes (denoting
as $S_{\theta}$ the pinned list):

$P_{\theta}({\bf x})=\frac{1}{Z_{\theta}}\left(\prod_{i\notin
S_{\theta}}P_{X}(x_{i})\right)\left(\prod_{i\in
S_{\theta}}\delta(x_{i}-x_{i}^{*})\right)\prod_{i<j<k}e^{\frac{\sqrt{3}\beta}{N}J_{ijk}x_{i}x_{j}x_{k}}$
(67)

####  B.3.1 Replica free entropy:

The replica free entropy for this model has been studied extensively in the
spin-glass literature, and can also be proven rigorously, see e.g. [86].

As for the previous models, the asymptotic free entropy is given by the
maximum of the so-called replica symmetric potential [20] :

$\displaystyle\frac{1}{N}{\mathbb{E}}_{{\bf x}_{0},{\bf z},J}\log Z_{\gamma}$
$\displaystyle\xrightarrow[N\to\infty]{}$ $\displaystyle{\rm
argmax}\;\Phi_{RS}(m)$ (68) $\displaystyle\Phi_{RS}(m)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{w,x_{0}}\left[\log
Z_{x}\left(\frac{m^{2}}{\Delta},\frac{m^{2}}{\Delta}x_{0}+\sqrt{\frac{m^{2}}{\Delta}}w\right)\right]-\frac{m^{3}}{3\Delta}$
(69)

where $m$ is the order parameter of the problem, $x_{0}\sim P_{X}$,
$w\sim\mathcal{N}(0,1)$ and $Z_{x}$ depends on the specific measure
considered. The same is valid for the pinned measure, as long as one
substitutes $Z_{\gamma}$ with $Z_{\theta}$.

While using the same derivation, in the following we will use
$T=\sqrt{3\Delta/2}$ as signal-to-noise parameter, since we are interested in
studying the unplanted model.

##### Tilted measure:

In the case of the tilted measure (26) we have

$\begin{split}Z_{x}(A,B;x_{0})&=\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)\\\
&=\sqrt{\frac{2\pi}{A+\gamma^{2}+1}}\exp\left(\frac{(B+\gamma^{2}x_{0}+\gamma
w)^{2}}{2(A+\gamma^{2}+1)}\right)\end{split}$ (70)

which leads to

$\Phi_{\rm
RS}(\chi)=\frac{\widetilde{\chi}}{2}+\frac{1}{2}\log\left(\frac{2\pi}{\widetilde{\chi}+1}\right)-\frac{1}{2T^{2}}\chi^{3}\,,\quad\widetilde{\chi}=\frac{3\chi^{2}}{2T^{2}}+\gamma^{2}\,.$
(71)

##### Pinning measure:

Meanwhile, considering the pinning measure (27) leads to

$\displaystyle
Z_{x}(A,B;x_{0})=\begin{cases}P_{X}(x_{0})\exp\left(Bx_{0}-\frac{A}{2}x_{0}^{2}\right)&\text{with
probability }\theta\\\
\sqrt{\frac{2\pi}{A+1}}\exp\left(\frac{B^{2}}{2(A+1)}\right)&\text{with
probability }1-\theta\end{cases}$ (72)

which in turn gives

$\Phi_{\rm
RS}(\chi)=\frac{\widetilde{\chi}}{2}+\frac{1-\theta}{2}\log\left(\frac{2\pi}{\widetilde{\chi}+1}\right)-\theta\log(\sqrt{2\pi}e)-\frac{\chi^{3}}{2T^{2}}\,,\quad\widetilde{\chi}=\frac{3\chi^{2}}{2T^{2}}\,.$
(73)

####  B.3.2 Message-passing algorithm

As already mentioned for the Ising version of the model, the Thouless-
Anderson-Palmer (TAP) equations [51] have been widely studied for the $p$-spin
models [117]. We use the results of [86], which presents the AMP algorithm for
the Spike Tensor model, and we use contiguity to derive equations valid for
the unplanted model with the tilting (or pinning) field. The mapping between
the two models is again given by $\Delta=\frac{2T^{2}}{3}$.

##### Tilted measure:

The equivalence between the tilted and the planted measure with external field
allows us to map the AMP iterations for the tilted measure to the ones for an
associated inference problem. Here we consider the formalism of [86], such
that the resulting equations are:

$\begin{cases}B_{i}^{t}=\frac{\sqrt{3}\beta}{N}\sum_{j<k}J_{ijk}\widehat{x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\sigma^{t}\widehat{x}_{i}^{t-1}{\bf\widehat{x}}^{t}\cdot{\bf\widehat{x}}^{t-1}\\\
\widehat{x}_{i}^{t+1}=\frac{B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}[{\bf
Y}_{t}]_{i}}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+\gamma^{2}+1}\,,\sigma^{t+1}=\frac{1}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+\gamma^{2}+1}\\\ \end{cases}$ (74)

where $\alpha(t)$ and $\beta(t)$ are the functions defining the interpolant
process, fixed at the start, and ${\bf Y}_{t}$ is the value of the noisy
observation at time $t$.

##### Pinning measure:

When considering the pinning measure (27), the AMP equations are only a slight
variation of the ones presented for the flow-based case.

Specifically, in autoregressive-based sampling we choose a fraction $\theta$
of the variables, for which we fix
$\widehat{x}_{i}^{t}=x_{0},\;\sigma_{i}^{t}=0$, which stems from the fact that
their posterior means are completely polarized on the solution. For the rest
of the variables, a fraction $1-\theta$, the AMP equations are exactly the
ones reported in Eq. (74), provided we fix $\gamma=0$.

The resulting algorithm is the following:

$\begin{cases}&\widehat{x}_{i}^{t+1}=\begin{cases}[{\bf x}_{0}]_{i}&\text{if
}i\in{S_{\theta}}\\\ \frac{B_{i}^{t}}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+1},&\text{otherwise}\end{cases}\,,\quad\sigma_{i}^{t+1}=\begin{cases}0&\text{if
}i\in{S_{\theta}}\\\ \frac{1}{3\beta^{2}\|\hat{\bf
x^{t}}\|_{2}^{2}/(2N)+1}&\text{otherwise}\\\ \end{cases}\\\
&B_{i}^{t}=\frac{\sqrt{3}\beta}{N}\sum_{j<k}J_{ijk}\widehat{x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\widehat{x}_{i}^{t-1}{\bf\widehat{x}}^{t}\cdot{\bf\widehat{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\end{cases}$
(75)

####  B.3.3 State evolution equations

As mentioned earlier, AMP iterations possess the salient property of being
able to be rigorously tracked by the State Evolution equations, that turn out
to be the fixed point of the replica potential in Eq. (68). Here, again we
follow the results presented in [86], such that as the Ising version of the
model, we have

$\displaystyle m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\rm
in}\left(\frac{(m^{t})^{2}}{\Delta},\frac{(m^{t})^{2}}{\Delta}x_{0}+\sqrt{\frac{(m^{t})^{2}}{\Delta}}w\right)x_{0}\right]$
(76)

where $x_{0}\sim P_{X}$, $w\sim\mathcal{N}(0,1)$ and $f_{in}$ is the input
channel and depends on the specific problem. We also remind the mapping
$T=\sqrt{3\Delta/2}$, which we will use in the following presentation.

##### Tilted measure:

For the tilted measure (26) we have

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma
wx-\gamma^{2}x^{2}/2)\exp(Bx-Ax^{2}/2)}\\\ &=\frac{B+\gamma^{2}x_{0}+\gamma
w}{A+\gamma^{2}+1}\end{split}$ (77)

which leads to the State Evolution equations

$\chi^{t+1}=\frac{\widetilde{\chi^{t}}}{1+\widetilde{\chi^{t}}}\,,\quad\widetilde{\chi}^{t}\equiv\frac{3}{2T^{2}}(\chi^{t})^{2}+\gamma^{2}\,.$
(78)

##### Pinning measure:

In the same way, for the pinning measure (27) we find

$\begin{split}f_{\rm
in}(A,B;x_{0})&=\begin{cases}\frac{P_{X}(x_{0})x_{0}\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}{P_{X}(x_{0})\exp\left(Bx_{0}-Ax_{0}^{2}/2\right)}&\text{with
probability }\theta\\\ \frac{\int\mathop{}\\!\mathrm{d}xxP_{X}(x)\exp(Bx-
Ax^{2}/2)}{\int\mathop{}\\!\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)}&\text{with
probability }1-\theta\end{cases}\\\ &=\begin{cases}x_{0}&\text{with
probability }\theta\\\ \frac{B}{A+1}&\text{with probability
}1-\theta\end{cases}\end{split}$ (79)

and thus the fixed point equations are

$\chi^{t+1}=\theta+(1-\theta)\frac{\widetilde{\chi^{t}}}{1+\widetilde{\chi^{t}}}\,,\quad\widetilde{\chi}^{t}\equiv\frac{3(\chi^{t})^{2}}{2T^{2}}$
(80)

####  B.3.4 Phase diagrams

The phase diagrams for the spherical $p$-spin are presented in Figure 9, and
we plot the same quantities as for the previous models. We observe again that
the flow-based method is advantageous with respect to the autoregressive one.

Figure 9: Phase diagrams for flow-based sampling (left) and autoregressive-
based sampling (right) for the Spherical p-spin model. On the x-axis we put
the temperature $T$ and on the y-axis the ratio
$\gamma^{2}=\alpha^{2}/\beta^{2}$ (left) and the decimated ratio $\theta$
(right). We compute the order parameter $\chi$, defined in (29), both from an
uninformed and an informed initialization, and we plot the difference between
the two. The dashed white lines are the spinodal lines, while the dashed black
one is the IT threshold, both defined at the beginning of the section. In both
plots we have that the dynamical transition is at $T_{d}=\sqrt{3/8}$, the
Kauzmann transition is at $T_{K}\approx 0.58$, while the tri-critical points
are at $T_{\rm tri}=2/3$ for flow-based and $T_{\rm tri}=\sqrt{1/2}$ for
autoregressive based sampling.

###  B.4 NAE-SAT model

####  B.4.1 Target model definition

The $k$-hypergraph bicoloring (or $k$-NAESAT) problem is a prototypical model
of constrained satisfaction problems defined on hypergraphs. An instance of
the problem is determined by an hypergraph $G=(V,E)$, where $V$ is the set of
$N$ vertices and $E$ the set of $M$ hyperedges, each containing exactly $k$
vertices.

Each vertex $i\in V$ is associated to an Ising-spin variable $x_{i}=\pm 1$,
while each hyperedge $a\in E$ is associated to a constraint involving the $k$
vertices entering in $a$ (in the following, we will use ${\bf x}_{\partial a}$
to indicate this set of variables).

For bicoloring, the $a$-th constraint is satisfied if there is at least one
$+1$ and one $-1$ among the $k$ variables of ${\bf x}_{\partial a}$. In terms
of probability distribution, this translates to

$P_{0}({\bf x})=\frac{1}{Z(G)}\prod_{a=1}^{M}\omega({\bf x}_{\partial
a})\,,\quad\omega(x_{1},\dots,x_{k})=\begin{cases}0&\text{if}\;\sum_{i=1}^{k}x_{i}=\pm
k\\\ 1&\text{otherwise}\end{cases}$ (81)

In the following, we will focus on the asymptotic limit where both $N$ and $M$
go to infinity, with constant rate $\alpha=M/N$.

This model has again the advantage to be contiguous to its planted version
[69]. We refer to [76, 77, 69] for rigorous results, and to [73, 74, 75, 118]
for results within the cavity approach.

####  B.4.2 BP equations and Bethe free entropy:

To analyse the properties of the model, one can use the cavity method [21]
from statistical physics, which allows to derive the BP equations for the
problem (see e.g. [118]), which can be written in terms of cavity messages as

$h_{i\rightarrow a}=f(\\{u_{b\rightarrow i}\\}_{b\in\partial_{i}\backslash
a})\,,\quad u_{a\rightarrow i}=g(\\{h_{j\rightarrow
a}\\}_{j\in\partial_{a}\backslash i})$ (82)

where

$f(u_{1},\dots,u_{d})=\frac{\prod_{i=1}^{d}(1+u_{i})-\prod_{i=1}^{d}(1-u_{i})}{\prod_{i=1}^{d}(1+u_{i})+\prod_{i=1}^{d}(1-u_{i})}$
(83)

is the function defining the messages going from factor nodes to variable
nodes and

$g(h_{1},\dots,h_{k-1})=\frac{\sum_{x_{1},\dots,x_{k}}\omega(x_{1},\dots,x_{k})x_{k}\prod_{i=1}^{k-1}(1+h_{i}x_{i})}{\sum_{x_{1},\dots,x_{k}}\omega(x_{1},\dots,x_{k})\prod_{i=1}^{k-1}(1+h_{i}x_{i})}$
(84)

is the one defining messages going from variable nodes to factor nodes.

Once a fixed point of the BP equations is reached, one can compute the Free
entropy from the resulting BP marginals. In its general form, it can be
written as

$\displaystyle\frac{1}{N}\ln Z(G)=\frac{1}{N}\sum_{i=1}^{N}\ln{\cal
Z}_{0}^{\rm v}(\\{u_{a\to i}\\}_{a\in{\partial
i}})+\frac{1}{N}\sum_{a=1}^{M}\ln{\cal Z}_{0}^{\rm c}(\\{h_{i\to
a}\\}_{i\in{\partial a}})-\frac{1}{N}\sum_{(i,a)}\ln{\cal Z}_{0}^{\rm
e}(h_{i\to a},u_{a\to i})\ ,$ (85)

where the last sum runs over the edges of the factor graph, and the local
partition functions are defined as:

$\displaystyle{\cal Z}_{0}^{\rm v}(u_{1},\dots,u_{d})$
$\displaystyle=\sum_{x}\prod_{i=1}^{d}\left(\frac{1+xu_{i}}{2}\right)\ ,$ (86)
$\displaystyle{\cal Z}_{0}^{\rm c}(h_{1},\dots,h_{k})$
$\displaystyle=\sum_{x_{1},\dots,x_{k}}{\omega}(x_{1},\dots,x_{k})\prod_{i=1}^{k}\left(\frac{1+x_{i}h_{i}}{2}\right)\
,$ (87) $\displaystyle{\cal Z}_{0}^{\rm e}(h,u)$
$\displaystyle=\sum_{x}\left(\frac{1+xh}{2}\right)\left(\frac{1+xu}{2}\right)\
.$ (88)

##### Replica symmetric cavity equations:

The simplest version of the cavity method implies the assumption of Replica
Symmetry (RS). This is rigorously justified in our case thanks to the fact
that we are considering a Bayes-Optimal model.

In such a case, the resulting self-consistent equations are relatively
simpler:

$\displaystyle\mathcal{P}^{RS}(h)$
$\displaystyle=\sum_{d=0}^{\infty}p_{d}\int\left(\prod_{i=1}^{d}\mathop{}\\!\mathrm{d}u_{i}\widehat{\mathcal{P}}^{RS}(u_{i})\right)\,\delta(h-f(u_{1},\dots,u_{d}))\
,$ (89) $\displaystyle\widehat{\mathcal{P}}^{RS}(u)$
$\displaystyle=\int\left(\prod_{i=1}^{k-1}\mathop{}\\!\mathrm{d}h_{i}\mathcal{P}^{RS}(h_{i})\right)\,\delta(u-g(h_{1},\dots,h_{k-1}))\
.$

where $\mathcal{P}^{RS}$ and $\widehat{\mathcal{P}}^{RS}$ are probability
distributions defined on the messages $h$ and $u$ respectively.

Being self-consistency equations defined on probability distributions makes
the computation of an analytical solution possible only in very restricted
cases, while in practise one needs to use numerical techniques to find
approximate solutions.

Still, population dynamics techniques (see for example [21]) have been shown
to provide very good approximations when one is interested in observables
defined as average quantities over $\mathcal{O}(N)$ cavity messages.

When the density of interactions $\alpha$ becomes large, the hypothesis
underlying the Replica Symmetric assumption breaks down, and we have to
consider the Replica Symmetry Breaking (RSB) phenomenon [119].

####  B.4.3 1RSB and Tree reconstruction equations:

The 1RSB cavity method aims to compute the potential

$\displaystyle\Phi_{1}(m)=\lim_{N\rightarrow\infty}\frac{1}{N}\log\left(\sum_{\gamma}Z_{\gamma}^{m}\right)$
(91)

where $m$ is the so-called Parisi parameter.

In their general formulation, the 1RSB equations are self-consistent equations
defined on distributions over probability distributions. Focusing on the case
$m=1$ allows simplifying considerably the equations, and the resulting
formulas correspond to the so-called Tree Reconstruction equations [120].
Moreover, the problem we are considering has a global spin-flip symmetry which
allows us to simplify further the equations, which can be finally written as

$\displaystyle Q_{+}^{(t+1)}(h)$
$\displaystyle=\sum_{d=0}^{\infty}p_{d}\int\left(\prod_{i=1}^{d}\mathop{}\\!\mathrm{d}u_{i}\widehat{Q}_{+}^{(t)}(u_{i})\right)\,\delta(h-f(u_{1},\dots,u_{d}))\
,$ (92) $\displaystyle\widehat{Q}_{+}^{(t)}(u)$
$\displaystyle=\sum_{x_{1},...,x_{k-1}}\tilde{p}(x_{1},...,x_{k-1}|+)\int\left(\prod_{i=1}^{k-1}\mathop{}\\!\mathrm{d}h_{i}Q_{+}^{(t)}(h_{i})\right)\delta(u-g(x_{1}h_{1},\dots,x_{k-1}h_{k-1}))\
,$

where

$\displaystyle\tilde{p}(x_{1},\dots,x_{k-1}|+)=\frac{{\omega}(x_{1},\dots,x_{k-1},+)}{\underset{x^{\prime}_{1},\dots,x^{\prime}_{k-1}}{\sum}{\omega}(x^{\prime}_{1},\dots,x^{\prime}_{k-1},+)}=\frac{\overset{k-1}{\underset{p=0}{\sum}}{\omega}_{p}\,\mathbb{I}\left[\overset{k-1}{\underset{i=1}{\sum}}x_{i}=k-1-2p\right]}{\overset{k-1}{\underset{p=0}{\sum}}\binom{k-1}{p}{\omega}_{p}}\
.$ (93)

Similar as before, we can compute the RS free entropy for the planted problem
as

$\displaystyle\Phi_{\rm RS}=$
$\displaystyle\sum_{d=0}^{\infty}p_{d}\int\left(\prod_{i=1}^{d}\mathop{}\\!\mathrm{d}u_{i}\widehat{Q}_{+}(u_{i})\right)\ln{\cal
Z}_{0}^{\rm
v}(u_{1},\dots,u_{d})+\alpha\int\left(\prod_{i=1}^{k}\mathop{}\\!\mathrm{d}h_{i}Q_{+}(h_{i})\right)\ln{\cal
Z}_{0}^{\rm c}(h_{1},\dots,h_{k})$ (94) $\displaystyle-\alpha
k\int\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}uQ_{+}(h)\widehat{Q}_{+}(u)\ln{\cal
Z}_{0}^{\rm e}(h,u)\ ,$

where ${\cal Z}_{0}^{\rm c}$, ${\cal Z}_{0}^{\rm e}$ and ${\cal Z}_{0}^{\rm
v}$ are the local partitions reported in Eq. (87), (88) and (86) respectively.

Finally, the order parameter for the problem, corresponding to the definition
in Eq. (29), can be computed from the population dynamics simulations as

$\chi=\int\mathop{}\\!\mathrm{d}hQ_{+}(h)h^{2}\,.$ (95)

####  B.4.4 Tilted and Pinned measures

In the preceding section, we presented the classical $k$-NAESAT model and how
its properties can be studied using the cavity method. Now, let us consider
the tilted and pinning measures and see how this changes the previous
equations.

##### Tilted measure:

Starting with the tilted measure, the added tilting field results directly in
the function defining the messages going from the factor nodes to the variable
nodes, which is modified to

$f^{\rm
tilt}(u_{1},\dots,u_{d})=\frac{e^{2\gamma(\gamma+z)}\underset{i=1}{\overset{d}{\prod}}(1+u_{i})-\underset{i=1}{\overset{d}{\prod}}(1-u_{i})}{e^{2\gamma(\gamma+z)}\underset{i=1}{\overset{d}{\prod}}(1+u_{i})+\underset{i=1}{\overset{d}{\prod}}(1-u_{i})}\,,\quad
z\sim\mathcal{N}(0,1)\,,$ (96)

while the function $g$, defined in (84), is not modified. In terms of the
partition function, the only local term that is modified is given by

${\cal Z}_{0}^{\rm
v}(u_{1},\dots,u_{d})=e^{\gamma(\gamma+z)}\prod_{i=1}^{d}\left(\frac{1+u_{i}}{2}\right)+e^{-\gamma(\gamma+z)}\prod_{i=1}^{d}\left(\frac{1-u_{i}}{2}\right)\,,\quad
z\sim\mathcal{N}(0,1)\,.$ (97)

while the other two terms remain the same.

##### Pinning measure:

For the pinning measure, again the function $g$ remains the same, but now we
have

$f^{\rm pinn}(u_{1},\dots,u_{d})=\begin{cases}1&\text{with probability
}\theta\\\
\frac{\underset{i=1}{\overset{d}{\prod}}(1+u_{i})-\underset{i=1}{\overset{d}{\prod}}(1-u_{i})}{\underset{i=1}{\overset{d}{\prod}}(1+u_{i})+\underset{i=1}{\overset{d}{\prod}}(1-u_{i})}&\text{with
probability }1-\theta\end{cases}\,,$ (98)

and again the only contribution to the partition function that is modified is

${\cal Z}_{0}^{\rm
v}(u_{1},\dots,u_{d})=\begin{cases}\prod_{i=1}^{d}\left(\frac{1+u_{i}}{2}\right)&\text{with
probability }\theta\\\
\prod_{i=1}^{d}\left(\frac{1+u_{i}}{2}\right)+\prod_{i=1}^{d}\left(\frac{1-u_{i}}{2}\right)&\text{with
probability }1-\theta\end{cases}\,.$ (99)

##### Phase diagrams:

In Fig. 10 we present the phase diagrams for the $k$-NAESAT problem,
considering the case $k=5$.

Figure 10: Phase diagrams for flow-based sampling (left) and autoregressive-
based sampling (right) for the $k$-NAESAT model. On the x-axis we put the
constraints to variables ratio $\alpha=M/N$ and on the y-axis the ratio
$\gamma^{2}=\alpha^{2}/\beta^{2}$ (left) and the decimated ratio $\theta$
(right). We compute the order parameter $\chi$, defined in (95), both from an
uninformed and an informed initialization, and we plot the difference between
the two. The dashed white lines are the spinodal lines, while the dashed black
one is the IT threshold, both defined at the beginning of the section. In both
plots we have that the dynamical transition is at $\alpha_{d}\approx 9.465$,
the Kauzmann (or condensation) transition is at $\alpha_{K}\approx 10.3$,
while the tri-critical points are at $\alpha_{\rm tri}\approx 8.4$ for both
flow-based and autoregressive based sampling. Thus, within our numerical
precision, the two methods seems to be equally efficient.

Compared to the plots reported in the main text, here we display directly the
difference $\chi_{\rm inf}-\chi_{\rm uninf}$, so that the coloured zones of
the plots are the ones displaying multiple fixed points, as opposed to the
black ones. We furthermore draw as white dashed lines the spinodal points, and
as a black dashed line the IT threshold, both defined at the beginning of
Appendix B.

Compared to the previous models, defined on dense graphs, here we don’t have
some fixed point equations defined on an $O(1)$ number of scalar parameter,
but instead some self-consistent equations on probability distributions. This
clearly makes it harder to have a precise estimate of the tri-critical points,
both for the tilted and the pinned measures. For the precision we were able to
achieve, both for flow-based sampling and for autoregressive-based sampling
the tri-critical point appears to be around $\alpha_{\rm tri}\approx 8.4$, and
we are not able to state if one of the two methods has a smaller gap compared
to the other and thus if there is a range in $\alpha$ where it is able to
sample efficiently where the other can not. A more careful analysis is needed
to clarify this point.

##  Appendix C Sampling simulations for Algorithm 1

We now show how the phenomena described through the previous phase diagram
influences the performance of the sampling algorithm in practice.

Figure 11: Spherical p-spin: $\chi(\gamma)$ for flow-based sampling. We
compare the results for the order parameter $\chi$ computed from the State
Evolution equations (black lines) to finite-size implementations of the
sampling algorithm in two regimes. Left: $T=\sqrt{\frac{3}{4}}>\Delta_{\rm
tri}$, for all values of $\gamma^{2}$, the SE equations have a unique fixed
point, and thus the initialization plays no role. The resulting curve (black
continuous line) is compared to algorithmic implementations of the flow-based
sampling algorithm, for sizes $N=200,400,800$. These curves, shown in
different colours above, match the asymptotic prediction. Right:
$T_{d}<T=\frac{9\sqrt{2}}{20}<T_{\rm tri}$, there is a range of values of
$\gamma^{2}$, for which the SE equations have two distinct fixed points. More
precisely, for all values between the IT point (reported as a dotted line) and
the informed spinodal point the model presents an algorithmically hard phase.
The uninformed/informed state evolution curves (black continuous/dashed lines
respectively) are compared to algorithmic implementations of the flow-based
sampling algorithm, for sizes $N=200,400,800$. Importantly, we show that in
this regime there is an evident mismatch with the asymptotic prediction.

In Fig. 11 we report some numerical simulations where we compare the
asymptotic curves derived through state evolution with some empirical
simulations implementing the flow-based sampling for the spherical $p$-spin
model.

Specifically, we compare what happens at a high value of temperature $T>T_{\rm
tri}$, where we expect the sampling scheme to work, to a value of $T$ in the
interval $[T_{d},T_{\rm tri}]$, where, as previously explained, we predict the
algorithm to fail.

In the first case, shown in the left part of the plot, we see that finite-size
simulation follows very well the theoretical prediction, already at a
considerably low number of variables. Conversely, in the second case, reported
in the right part, the situation is different, and the curves show a gap.

Moreover, we check when the Nishimori conditions are satisfied for different
values of temperature $T$ in the spherical $p$-spin. We do this by analysing
the following observable:

$\displaystyle\text{OV}(\gamma)\equiv\frac{1}{N}\mathbb{E}\left[{\bf
Y_{t}}\langle{\bf x}\rangle\right]$ (100)

which is the overlap between the observation vector ${\bf Y_{t}}$ and the
average magnetization $\langle{\bf x}\rangle$.

Indeed, by definition of the observation process we can write

$\displaystyle\text{OV}(\gamma)=\frac{\alpha(t)}{N}\mathbb{E}\left[{\bf
x_{0}}\langle{\bf x}\rangle\right]+\frac{\beta(t)}{N}\mathbb{E}\left[{\bf
z}\langle{\bf x}\rangle\right]$ (101)

and furthermore by using Stein’s lemma [121] we can write it as

$\displaystyle\begin{aligned}
\text{OV}(\gamma)&=\frac{\alpha(t)}{N}\mathbb{E}\left[{\bf x_{0}}\langle{\bf
x}\rangle\right]+\frac{\beta(t)}{N}\mathbb{E}\left[\partial_{\bf z}\langle{\bf
x}\rangle\right]\\\ &=\frac{\alpha(t)}{N}\mathbb{E}\left[{\bf
x_{0}}\langle{\bf
x}\rangle\right]+\frac{\beta(t)}{N}\mathbb{E}\left[\frac{\alpha(t)}{\beta(t)}\langle\|{\bf
x}\|_{2}^{2}\rangle-\frac{\alpha(t)}{\beta(t)}\langle{\bf
x}\rangle^{2}\right]\\\ &=\alpha(t)\frac{1}{N}\mathbb{E}\left[{\bf
x_{0}}\langle{\bf x}\rangle+\langle\|{\bf x}\|_{2}^{2}\rangle-\langle{\bf
x}\rangle^{2}\right]\approx\alpha(t)+\alpha(t)\frac{1}{N}\mathbb{E}\left[{\bf
x_{0}}\langle{\bf x}\rangle-\langle{\bf x}\rangle^{2}\right]\,.\end{aligned}$
(102)

where in the last step we used the fact that $\lim_{N\rightarrow\infty}\|{\bf
x}\|_{2}^{2}/N=1$. Now, if the Nishimori conditions are satisfied, we have
that $\mathbb{E}\left[{\bf x_{0}}\langle{\bf
x}\rangle\right]=\mathbb{E}\left[\langle{\bf x}\rangle^{2}\right]$. Putting
this equivalence back into (102), we see that in this case OV coincides with
the function $\alpha(t)$ that defines the interpolant process.

In Fig. 12, we use this equivalence and compare with $\alpha(t)=1-t$, the
behaviour of OV for $T=\frac{9\sqrt{2}}{20}>T_{\rm tri}$ and for
$T_{d}<T=\sqrt{\frac{3}{4}}<T_{\rm tri}$, showing that in the first case we
are Bayes optimal for all the values of $\gamma$, and thus at each sampling
step. Instead, in the second case, even if at the beginning the two curves
coincide, around the IT threshold they develop a gap, and after this point
they take two different paths. This behaviour is consistent with what we
observed from the study of the order parameter $\chi$, for which also a gap
with the state evolution equations was developed around the IT threshold.

Figure 12: Spherical p-spin: Checking the Nishimori conditions for flow-based
sampling. We compare the results for the overlap OV, defined in (100) computed
from finite-size $(N=200,400,800)$ implementations of the sampling algorithm
to the behaviour with $\gamma$ of the function $\alpha(t)=1-t$ in two regimes.
Left: $T=\sqrt{\frac{3}{4}}>\Delta_{\rm tri}$. The simulations, shown in
different colours above, are very close to the asymptotic prediction, even if
the sizes are relatively limited. Right: $T_{d}<T=\frac{9\sqrt{2}}{20}<T_{\rm
tri}$. Here, after the IT threshold there is an evident mismatch between the
two curves, due to the departure of the algorithm from the Bayes-optimality
regime.

##  Appendix D Bayes optimal inference: Concentration, Replica Symmetry and
Nishimori identities

We shall briefly recall here some of the important properties of the measure
associated with optimal Bayesian denoising, that we are using in the main
text. All of these properties are well known in the literature, but we recall
them for completeness

The first one are the well-known Nishimori identities [63, 64, 23], which are
valid for any Bayesian posterior estimation problem where one observes a
variable $Y$ sampled from $P(Y|X)$, and attempt to reconstruct $X$ by
computing the posterior average. We reproduce the theorem and the proof here

######  Theorem 1 (Nishimori Identity).

Let $X^{(1)},\ldots,X^{(k)}$ be $k$ i.i.d. samples (given $Y$) from the
distribution $P\left(X=\cdot\,\middle|\,Y\right)$. Denoting
$\left\langle\cdot\right\rangle$ the “Boltzmann” expectation, that is the
average with respect to the $P\left(X=\cdot\,\middle|\,Y\right)$, and
$\mathbb{E}\left[\cdot\right]$ the “Disorder” expectation, that is with
respect to $(X^{*},Y)$. Then for all continuous bounded function $f$ we can
switch one of the copies for $X^{*}$:

$\mathbb{E}\left[\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{(k)}\right)\right\rangle_{k}\right]=\mathbb{E}\left[\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{*}\right)\right\rangle_{k-1}\right]$
(103)

###### Proof.

The proof is a consequence of Bayes theorem and of the fact that both $x^{*}$
and any of the copy $X^{(k)}$ are distributed from the posterior distribution.
Denoting more explicitly the Boltzmann average over $k$ copies for any
function $g$ as

$\displaystyle\left\langle
g(X^{(1)},\ldots,X^{(k)})\right\rangle_{k}=\vcentcolon\int\prod_{i=1}^{k}\mathop{}\\!\mathrm{d}{x}_{i}P(x_{i}|Y)g(X^{(1)},\ldots,X^{(k)})\,$
(104)

we have, starting from the right-hand side

$\displaystyle\mathbb{E}_{Y,X^{*}}\left[\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{*}\right)\right\rangle_{k-1}\right]$
$\displaystyle=\int dx^{*}dyP(x^{*}|Y)P(Y)\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{*}\right)\right\rangle_{k-1}$
$\displaystyle=\mathbb{E}_{Y}\int\mathop{}\\!\mathrm{d}{x}^{k}P(x^{k}|y)\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{k}\right)\right\rangle_{k-1}$
$\displaystyle=\mathbb{E}_{Y}\left[\left\langle
f\left(Y,X^{(1)},\ldots,X^{(k-1)},X^{(k)}\right)\right\rangle_{k}\right]$

∎

In particular, we have the relation $\mu=\chi$, as stated in the main text.

We now move to the specific case of Gaussian denoising, and more specifically
to the measure:

$P_{\gamma}({\bf x})=\frac{1}{Z_{n}}\exp\left(\gamma^{2}\langle{\bf
x_{0}},{\bf x}\rangle+\gamma\langle{\bf z},{\bf
x}\rangle-\frac{\gamma^{2}}{2}\|{\bf x}\|^{2}\right)P_{0}({\bf x})\,$ (105)

The first identity is that the derivative of the free entropy associated with
this problem is simply the expected overlap $\mu$ (This is often called in
slightly different context the I-MMSE theorem [122]) and the second derivative
the (Boltzmann) variance of the overlap (this is often called the Fluctuation-
Dissipation theorem in statistical mechanics):

######  Lemma 1 (First(I-MMSE theorem[65]) and second (FDT theorem [123])
derivative of the free entropy).

Consider the free entropy density associated with the measure (105):

$f_{N}=\frac{1}{N}{\mathbb{E}}[\log Z_{n}(\gamma)]$ (106)

then

$\partial_{\gamma}f_{n}=\gamma\mu(\gamma)\equiv\frac{\gamma}{N}\mathbb{E}[\langle{\bf
x}(\gamma)\rangle{\bf x}_{0}]$ (107)

and

$\partial^{2}_{\gamma}f_{N}=N\gamma^{2}\mathbb{E}\left[\left(\langle(\frac{{\bf
x}(\gamma)\cdot{\bf x}_{0}}{N})^{2}\rangle\right)-\left(\frac{\langle{\bf
x}(\gamma)\rangle{\bf
x}_{0}}{N}\right)^{2}\right]=N\gamma^{2}{\mathbb{E}}[{\rm var}\left(\frac{{\bf
x}\cdot{\bf x}_{0}}{N}\right)]$ (108)

###### Proof.

The proof is a direct application of Nishimori identities together with Stein
lemma (which states that ${\mathbb{E}}_{Z}zg(Z)={\mathbb{E}}g^{\prime}(z)$ for
a Gaussian random variable $Z$), which we reproduce here:

$\displaystyle\partial_{\gamma}f_{n}$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\int d{\bf x}P_{\gamma}({\bf x})\frac{2\gamma{\bf
x}_{0}\cdot{\bf x}+{\bf x}\cdot{\bf z}-\gamma{\bf x}\cdot{\bf x}}{N{\bf
z}(\gamma)}$ (109) $\displaystyle=$ $\displaystyle
2\gamma\mu-\gamma{\mathbb{E}}\langle\frac{{\bf x}\cdot{\bf
x}}{N}\rangle-{\mathbb{E}}\int d{\bf x}P_{\gamma}(X)\frac{X\cdot Z}{N{\bf
z}(\gamma)}$ (110) $\displaystyle=$ $\displaystyle
2\gamma\mu-\gamma{\mathbb{E}}\langle\frac{{\bf x}\cdot{\bf
x}}{N}\rangle+\gamma{\mathbb{E}}\langle\frac{{\bf x}\cdot{\bf
x}}{N}\rangle-\gamma\frac{1}{N}\langle{\bf x}\rangle\cdot\langle{\bf
x}\rangle$ (111) $\displaystyle=$ $\displaystyle 2\gamma\mu-\chi=\gamma\mu\ $
(112)

where we use Stein lemma on line 111 and Nishimori on line 112. The second
derivative identity is obtained along the same line by deriving twice with
respect to $\gamma$. ∎

This lemma in turn implies that the concentration of the overlap $\mu$, a
trick often used in the literature of mathematical physics when proving the
replica equation (see e.g. [67]):

######  Theorem 2 (Concentration of overlaps).

Almost everywhere in $gamma$, we have for some $K(\gamma)$, that

${\mathbb{E}}[{\rm var}\left(\frac{{\bf x}\cdot{\bf
x}_{0}}{N}\right)]\to_{N\to\infty}0$

###### Proof.

From the derivative, we have that

$\displaystyle\int_{\gamma_{1}}^{\gamma_{2}}{\rm var}_{\gamma}\left(\frac{{\bf
x}\cdot{\bf
x}_{0}}{N}\right)=\frac{1}{\gamma^{2}N}(\mu(\gamma_{2})-\mu(\gamma_{1}))\leq\frac{K}{N}$
(113)

where we have assumed that $\mu$ is bounded by some constant (which is the
case for the discrete variables discussed in this paper where $-1<\mu<1$). As
a consequence, almost everywhere in $\gamma$, the variance of the overlap must
vanish. ∎

Additionally, one can also prove the concentration of the variance with
respect to the disorder, see [67].

Finally, similar results exist in the case of “pinning” [68, 69].

##  Appendix E Analysis of Algorithm 1

In this section, we provide a theoretical analysis of the performance of
Algorithm 1 for the “efficient” regime in Figure 1. In these regimes, we
conjecture being able to approximate perfect denoising during the
interpolation path. Therefore, to simplify our analysis, we assume access to a
perfect denoiser and primarily focus on the validity of the continuous-time
limit and the effect of discretization. Additional approximation and
discretization errors due to the denoiser, such as AMP can be incorporated in
our analysis through a straightforward manner. For a rigorous analysis of the
approximation error and Lipschitzness of the AMP iterates in the case of the
SK and spiked matrix models, we refer the reader to [36, 37].

Lastly, throughout the analysis, we shall assume that Algorithm 1 is run from
time $t=1$ to $t=\epsilon$ for some $\epsilon>0$. This ensures that we can
avoid singularities at $t=0$, while choosing $\epsilon$ arbitrarily close to
$0$ allows us to approximate the target measure up to arbitrary accuracy.

Recall the definition of the “tilted” posterior measure at time $t$:

$P({\bf x}|{\bf y}(t)={\bf Y}_{t})=\frac{1}{Z({\bf
Y}_{t})}\exp\left(\frac{\alpha(t)}{\beta(t)^{2}}\langle{\bf Y}_{t},{\bf
x}\rangle-\frac{\alpha(t)^{2}}{2\beta(t)^{2}}||{\bf x}||^{2}\right)P_{0}({\bf
x}).$ (114)

Here the factor $\frac{1}{Z({\bf
Y}_{t})}\exp\left(\frac{\alpha(t)}{\beta(t)^{2}}\langle{\bf Y}_{t},{\bf
x}\rangle-\frac{\alpha(t)^{2}}{2\beta(t)^{2}}||{\bf x}||^{2}\right)$ is
interpreted as the Radon-Nikodym derivative of $P({\bf x}|{\bf y}(t)$ w.r.t
$P_{0}$. Throughout the present section, we shall denote $P({\bf x}|{\bf
y}(t)={\bf Y}_{t})$ by $P_{t,{\bf Y}_{t}}$. We shall assume that
$\alpha(t),\beta(t)$ are continuously-differentiable.

###  E.1 Fluctuation-dissipation and Lipschitzness

A crucial quantity related to the well-posedness of the continuity equation
(Eq. (4)) and the validity of Algorithm 1 is the Lipschitz constant of the
vector field. We therefore start by present a preliminary result, which can be