Sampling with flows, diffusion and autoregressive neural networks: A
  spin-glass perspective
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.
Davide Ghio, Yatin Dandi, Florent Krzakala, Lenka Zdeborová
# Sampling with flows, diffusion and autoregressive neural networks:

###### 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 [37; 38]. 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 [36; 38; 42].
* 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](https://github.com/IdePHICS/DiffSamp).

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 [46; 19], 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_{\text{tri}}\), that depends on the detail of the method, such that for \(T_{d}<T<T_{\text{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; 45; 48]. 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_{\text{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_{\text{alg}}\) and fail to learn in the "hard"

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_{\text{alg}}<\Delta<\Delta_{\text{IT}}\) where efficient algorithms achieving close-to-optimal estimation error are not known and conjectured not to exist [45].
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}\).

### 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 [36; 42] 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 [48; 24]. 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; 39; 56].

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; 4; 5]. The goal is to sample a vector \(\mathbf{x}_{0}\in\mathbb{R}^{N}\) from a so-called target distribution \(P_{0}(\mathbf{x}_{0})\). Continuous-time flow-based models [3; 4; 5; 58] achieve this by constructing an ODE (flow) that continuously transforms Gaussian noise \(\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbb{I}_{N})\) to a sample \(\mathbf{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 \(\mathbf{x}_{0}\sim P_{0}\) to a Gaussian noise \(\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbb{I}_{N})\). Such a continous-time reverse process underlies a wide range of flow-models [4; 5] as well as the score-based diffusion models [9]. In practice, we would only observe samples from the distribution \(P_{0}(\mathbf{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}(\mathbf{x}_{0})\) is known as this can only render the task easier than when only samples from \(P_{0}(\mathbf{x}_{0})\) are available.

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

\[\mathbf{y}(t)=\alpha(t)\mathbf{x}_{0}+\beta(t)\mathbf{z},\quad t\in\left[0,1 \right], \tag{1}\]

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

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

such that indeed we have \(\mathbf{y}(0)=\mathbf{x}_{0}\) and \(\mathbf{y}(1)=\mathbf{z}\). Suppose that \(P_{0}(\mathbf{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(\mathbf{y}(t))\) associated to the measure \(P_{t}\) of the random variable \(\mathbf{y}(t)\) satisfies the following transport equation:

\[\partial_{t}\rho(\mathbf{y},t)+\nabla\cdot\left(b(\mathbf{y},t)\rho(\mathbf{y },t)\right)=0\,, \tag{4}\]

where we defined the _velocity field_

\[b(\mathbf{y},t)=\mathbb{E}[\partial_{t}\mathbf{y}(t)|\mathbf{y}(t)=\mathbf{y }]=\mathbb{E}[\dot{\alpha}(t)\mathbf{x}_{\mathbf{0}}+\hat{\beta}(t)z|\mathbf{ y}(t)=\mathbf{y}]\,. \tag{5}\]

Indeed, \(b(\mathbf{y},t)\) is simply the expected velocity of \(\mathbf{y}(t)\) conditioned on being at \(\mathbf{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}(\mathbf{x}_{\mathbf{0}})\) is discrete.
Equation (1) defines a forward process interpolating from \(\mathbf{x_{0}}\) to \(\mathbf{z}\). Equation (4) further reveals that, in law, \(\mathbf{y}(t)\) can be obtained by applying the vector field (Eq. (5)) starting from \(\mathbf{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\mathbf{Y}(t)}{dt}=b(\mathbf{Y}(t),t)\,. \tag{6}\]

Eq. (4) then implies that the random variable defined by Eq. (6) solved backwards in time from the final value \(\mathbf{Y}_{t=1}=\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbb{I}_{N})\) is also distributed according to \(\rho(\mathbf{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}(\mathbf{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

\[\mathbf{Y}_{t-\delta}=\mathbf{Y}_{t}-\delta\,b(\mathbf{Y}_{t},t)\,. \tag{7}\]

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

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

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

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

which, given the initial condition \(\mathbf{Y}_{t=1}=\mathbf{z}\), will evolve back to a sample from the target \(\mathbf{Y}_{t=0}\sim P_{0}\), provided that we are able to estimate \(\mathbb{E}\left[\mathbf{x}_{0}|\mathbf{y}(t)=\mathbf{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.

### Diffusion-based models and SDEs

The flow-based algorithm [4] (Alg. 1) relies on replicating the law of the interpolant \(\mathbf{y}(t)\) defined by Eq. (1) through the deterministic ordinary differential equation (ODE) defined by Eq. (6). Interestingly, however, the _same_ law for \(\mathbf{y}(t)\) can be obtained via a stochastic differential equation (SDE). Similarly, when \(\mathbf{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) [9; 59]. 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[\mathbf{x}_{0}|\mathbf{y}(t)=\mathbf{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; 9; 41]. Since the law of \(\mathbf{y}(t)\) remains the same with the ODE or SDE-based approaches, the deterministic evolution of the "tilted" measure \(P(\mathbf{x}|\mathbf{y}(t)=\mathbf{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.
### The posterior average and Bayes-optimal denoising

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

The Bayes-optimal denoising estimator \(\hat{\mathbf{x}}\) (optimal in the sense of minimizing the mean-squared-error between the estimator and \(\mathbf{x}_{0}\)) is to use the posterior mean \(\hat{\mathbf{x}}=\mathbb{E}\left[\mathbf{x}_{0}|\mathbf{y}(t)=\mathbf{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}(\mathbf{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 \(\mathbf{y}(t)=\mathbf{Y}_{t}\) is given by

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

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

We now recall that the law of \(\mathbf{Y}_{t}\) can be obtained through an observation of a sample from the target distribution \(\mathbf{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 \(\mathbf{x_{0}}\) as the "truth" signal and make the difference between the dummy variable \(\mathbf{x}\) in the integral. We can hence rewrite the measure in Eq. (10) further as

\[P_{\gamma}(\mathbf{x}|\mathbf{x}_{0},\mathbf{z})\!\propto\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
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}(\mathbf{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}(\mathbf{x})=P_{0}(x_{1})P_{0}(x_{2}|x_{1})\ldots P_{0}(x_{N}|x_{1}, \ldots,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}(\mathbf{x}|\mathbf{x}_{0},S_{\theta})\propto P_{0}(\mathbf{x})\prod _{i\in S_{\theta}}\delta([x_{i}-[\mathbf{x}_{0}]_{i}) \tag{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 \(\mathbf{x}_{0}\) and consequently call the measure \(P_{\theta}(\mathbf{x}|\mathbf{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 \(\mathbf{x}_{0}\) for an expected fraction \(\theta\) of components. We are again interested in the marginals of \(P_{\theta}\) that we will denote \(\hat{\mathbf{x}}(\theta)\). The evolution of the pinning measure can also be interpreted as a coordinate-by-coordinate stochastic localization process [35].

```
Input: BEC denoiser for\(l=1,\ldots,N\)do  Compute the posterior marginal \(P_{l}(x_{l}|x_{1},\ldots,x_{l-1})\) with BEC denoiser Eq. (15)  Assign \(x_{l}\sim P_{l}(x_{l}|x_{1},\ldots,x_{l-1})\) endfor Return\(x\)
```

**Algorithm 2** Autoregressive-network-based sampling algorithm

## 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 \(\mathbf{y}\). In the case of diffusion, we have access to an observable where a sample from the target \(\mathbf{x}_{0}\sim P_{0}\) is polluted by an AWGN channel:

\[\mathbf{y}=\alpha\mathbf{x}_{0}+\beta\mathbf{z}\,, \tag{14}\]

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

\[y_{i}=\begin{cases}[\mathbf{x}_{0}]_{i}&\text{with probability }\theta\\ *&\text{otherwise}\end{cases} \tag{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_:

\[\mu(\gamma) \equiv \frac{1}{N}\mathbb{E}[\dot{\mathbf{x}}(\gamma)\cdot\mathbf{x}_{ 0}]\,, \tag{16}\] \[\chi(\gamma) \equiv \frac{1}{N}\mathbb{E}[\|\dot{\mathbf{x}}(\gamma)\|^{2}]\,. \tag{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 \(\mathbf{x}_{0}\) and \(\mathbf{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 [23; 64], 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 [72; 36]). 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\)): \[\mathcal{H}(\mathbf{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}(\mathbf{x})\propto\exp(-\beta\mathcal{H}(\mathbf{x}))\) with \(\beta=1/T\). In the Ising case, we take \(x_{i}=\pm 1\), for \(i=1,\dots,N\), in the spherical case \(\mathbf{x}\in\mathcal{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}(\mathbf{x})\propto\prod_{a=1}^{M}\omega(\mathbf{x}_{\partial a})\,,\quad \mathbf{x}\in\{-1,+1\}^{N}\] (19) \[\alpha=\frac{M}{N}\,,\quad\omega(x_{1},\ldots,x_{k})=\begin{cases}0&\text{if } \sum_{i=1}^{k}x_{i}=\pm k\\ 1&\text{otherwise}\end{cases}\] where by \(\mathbf{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; 78; 29].

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_{\text{tri}}\) and constraint density \(\alpha_{\text{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_{\Gamma\Gamma}>\Delta>\Delta_{\text{alg}}\) that is computationally hard for message-passing algorithms and conjectured hard for any other efficient algorithm [45]. A recent line of work [80; 30; 49; 50; 79], 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_{\text{alg}}\) up to some \(\Delta_{\text{MCMC}}\) that depends on the specific algorithm. Yet, diffusion models were shown to be able to sample down to \(\Delta_{\text{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 \(\mathbf{x}^{*}\) and aims at finding back \(\mathbf{x}^{*}\) from the matrix. This is done by sampling from the posterior probability, and therefore one considers the probability distribution: \[P_{0}(\mathbf{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)\,,\] (20) \[J_{ij}=\frac{x_{i}^{*}x_{j}^{*}}{\sqrt{N}}+z_{ij}\,,\quad z_{ij} =z_{ji}\sim\mathcal{N}(0,\Delta)\,;\] \[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.\]

## 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 \(\mathbf{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. [36, 85]) is to notice that \(\forall T>T_{K}\) the original \(p\)-spin model is contiguous to its "planted" version, where a vector \(\mathbf{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_{\mathrm{K}}\) are given by

\[\chi^{*} =\operatorname*{argmax}_{\chi}\Phi_{\mathrm{RS}}(\chi) \tag{21}\] \[\Phi_{\mathrm{RS}}(\chi) =\frac{\widetilde{\chi}}{2}+\frac{1}{2}\log\left(\frac{2\pi}{ \widetilde{\chi}+1}\right)-\frac{1}{2T^{2}}\chi^{3}\,,\] \[\widetilde{\chi} =\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

Figure 2: Phase diagrams for the tilted measure \(P_{\mathrm{r}}\) (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 (NAESTAT). 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_{\mathrm{K}}\approx 0.58\), while the tri-critical point is at \(T_{\mathrm{tri}}=2/3\) for diffusion and \(T_{\mathrm{tri}}=\sqrt{1/2}\) for autoregressive. In (b) the Ising \(p\)-spin the values are: \(T_{d}\approx 0.682\), \(T_{\mathrm{K}}\approx 0.652\), \(T_{\mathrm{tri}}\approx 0.741\) for diffusion, and \(T_{\mathrm{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_{\mathrm{K}}/\rho^{2}\approx 1.029\), and \(\Delta_{\mathrm{alg}}/\rho^{2}\approx 0.981\). The tri-critical points are at \(\Delta_{\mathrm{tri}}/\rho^{2}\approx 1.08\) for diffusion, and \(\Delta_{\mathrm{tri}}/\rho^{2}\approx 1.069\) for autoregressive. In (d) the bicoloring the values are \(\alpha_{d}\approx 9.465\), \(\alpha_{\mathrm{K}}\approx 10.3\). The tri-critical points are \(\alpha_{\mathrm{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.
\(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) [45, 48, 52, 23]. This is nothing else than the well-known metastability problem in thermodynamics, as one is "trapped" in the wrong maxima of (21).

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) [87, 31]. 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_{ij\hat{x}}\widehat {x}_{j}^{t}\widehat{x}_{k}^{t}-\frac{3}{N}\beta^{2}\widehat{x}_{i}^{t-1} \sigma^{t}\widehat{\mathbf{x}}^{t}\cdot\widehat{\mathbf{x}}^{t-1}\\ \widehat{x}_{i}^{t+1}=\frac{B_{i}^{t}+\frac{\alpha(t)}{\beta(t)T}\left[\mathbf{ Y}_{i}\right]}{3\beta^{2}\left\|\mathbf{x}^{k}\right\|_{2}^{2}/(2N)+\gamma^{2}+1}\,, \sigma^{t+1}=\frac{1}{3\beta^{2}\left\|\mathbf{x}^{k}\right\|_{2}^{2}/(2N)+ \gamma^{2}+1}\,.\end{cases}\]

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}\,, \tag{22}\]

which is nothing but the fixed point equation of Eq. (21). 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 [90, 36], they show that the AMP iterates satisfy Lipschitz-continuity w.r.t the observation \(\mathbf{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.

### 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^{*}=\operatorname{argmax}_{\chi}\Phi_{\rm RS}(\chi)\) where the RS free entropy reads

\[\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 }}\,,\] \[\widetilde{\chi} =\frac{3\chi^{2}}{2T^{2}}\,. \tag{23}\]

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.

### 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.

### 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. [94; 70], 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.

###### Acknowledgements.

 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.
## 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}(\mathbf{x}|\mathbf{x}_{0},\mathbf{z})\propto e^{\gamma(t)^{2}( \mathbf{x},\mathbf{x}_{0})+\gamma(t)(z,\mathbf{x})-\frac{\gamma(t)^{2}}{2} \left\|\mathbf{x}\right\|^{2}}P_{0}(\mathbf{x})\,. \tag{10}\]

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 \(\mathbf{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. Sample \(\mathbf{x_{0}}\) uniformly. Then sample the interaction matrix \(J\) as \(J=\frac{\beta}{n}\mathbf{x_{0}}\mathbf{x_{0}}^{\top}+W\) where \(W\sim\text{GOE}\).
2. Sample \(J\sim\text{GOE}\). Then sample \(\mathbf{x_{0}}\sim P_{\text{SK}}(J)\),

where \(P_{\text{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 [86; 87]: 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 (10) 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; 32; 67; 103], 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 [83; 87] 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. (10) 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}(\mathbf{x}|\mathbf{x}_{0},S_{\theta})\propto P_{0}(\mathbf{x})\prod _{i\in S_{\theta}}\delta([x_{i}-[\mathbf{x}_{0}]_{i}) \tag{11}\]

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:

\[\mu(\gamma) \equiv \frac{1}{N}\mathbb{E}[\hat{\mathbf{x}}(\gamma)\cdot\mathbf{x}_{0 }]\,, \tag{12}\] \[\chi(\gamma) \equiv \frac{1}{N}\mathbb{E}[\|\hat{\mathbf{x}}(\gamma)\|^{2}]\,, \tag{13}\]

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.

### Sparse rank-one matrix factorization

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

Given a hidden vector \(\mathbf{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\,.\]

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:

\[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 \(\mathbf{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, 108, 81, 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}(\mathbf{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)\,. \tag{14}\]

With this distribution, \(\mathbf{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 \(\mathbf{x}_{0}\) is equivalent to the one where the vector \(\mathbf{x}_{0}\) is planted. Therefore, in what follows, we shall assume that \(\mathbf{x}_{0}\) corresponds to the planted configuration. We thus obtain the following tilted measure for diffusion and flow-based models:

\[P_{\gamma}(\mathbf{x})\ =\ \frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right)e^{ \gamma(t)^{2}(\mathbf{x},\mathbf{x}_{0})+\gamma(t)\langle\mathbf{z},\mathbf{ x}\rangle-\frac{\gamma(t)^{2}}{2}\|\mathbf{x}\|^{2}}\left(\prod_{i<j}e^{\frac{1}{ \Delta\sqrt{N}}\frac{z_{ij}x_{i}x_{j}}{2}+\frac{1}{2\Delta^{3}}x_{i}x_{j}[ \mathbf{x}_{0}]_{i}[\mathbf{x}_{0}]_{j}-\frac{1}{2\Delta^{3}}\pi^{2}_{i}x_{j}^ {2}}\right) \tag{10}\]

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}\)):

\[P_{\theta}(\mathbf{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)\left(\prod_{i<j}e^{\frac{1}{\Delta\sqrt{N}}\frac{z_{ij}x_{i}x_{j}+ \frac{1}{2\Delta^{3}}x_{i}x_{j}[\mathbf{x}_{0}]_{i}[\mathbf{x}_{0}]_{j}-\frac{ 1}{2\Delta^{3}}\pi^{2}_{i}x_{j}^{2}}}\right) \tag{11}\]

#### b.2.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, 84, 103, 104, 108, 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] :

\[\frac{1}{N}\mathbb{E}_{\mathbf{x}_{0},\mathbf{z},\tilde{\mathbf{ z}}}\log Z_{\gamma}\ \xrightarrow[N\to\infty]{}\ \operatorname{argmax}\Phi_{RS}(m) \tag{12}\] \[\Phi_{RS}(m)\quad=\quad\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} \tag{13}\]

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** (10), we have

\[Z_{x}(A,B;x_{0}) =\int\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2 }x^{2}/2)\exp(Bx-Ax^{2}/2) \tag{14}\] \[=\rho e^{-(A+\gamma^{2})/2}\cosh(B+\gamma w+\gamma^{2}x_{0})+(1-\rho)\]

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

\[\Phi_{\mathrm{RS}}(\chi)= \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] \tag{15}\] \[+(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** (15) leads to

\[Z_{x}(A,B;x_{0}) =\begin{cases}\int\mathrm{d}xP_{X}(x)\delta_{x,x_{0}}\exp(Bx-Ax^{2 }/2)&\text{with probability }\theta\\ \int\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)&\text{with probability }1-\theta \end{cases} \tag{16}\] \[=\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}\]

which in turn gives

\[\Phi_{\mathrm{RS}}(\chi)= \theta\left(\rho\log\rho+(1-\rho)\log(1-\rho)+\frac{\rho}{2} \widetilde{\chi}\right) \tag{17}\] \[+(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]\] \[+(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}\,.\]
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. (14) 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}}(\mathbf{Y}_{i}))}{\rho+\frac{(1-\rho)\exp(\alpha t^{t}+\gamma ^{2})/2}{\cosh(B_{i}^{t}+\frac{\alpha(t)}{\beta(t)^{2}}(\mathbf{Y}_{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}}(\mathbf{Y}_{i}))}{\left(\rho\cosh(B_{i }^{t}+\frac{\alpha(t)}{\beta(t)^{2}}(\mathbf{Y}_{i}))+(1-\rho)\exp((A+\gamma^ {2})/2)\right)^{2}}\\ A^{t}=\frac{\|\widehat{\mathbf{x}}^{t}\|_{2}^{2}}{\Delta N};\quad B_{i}^{t}= \frac{1}{\Delta\sqrt{N}}\mathbf{J}_{1}\cdot\widehat{\mathbf{x}}^{t}-\frac{1} {N\Delta}\widehat{x}_{i}^{t-1}\sum_{k}\sigma_{k}^{t}\end{cases} \tag{15}\]

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

_Pinning measure:_ When considering the pinning measure (14), 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}=[\mathbf{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. (15), at \(\gamma=0\).

The resulting algorithm is thus

\[\begin{cases}\quad\widehat{x}_{i}^{t+1}=\begin{cases}[\mathbf{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{\|\widehat{\mathbf{x}}^{t}\|_{2}^{2}}{\Delta N};\quad B_{i}^{t}= \frac{1}{\Delta\sqrt{N}}\mathbf{J}_{1}\cdot\widehat{\mathbf{x}}^{t}-\frac{1} {N\Delta}\widehat{x}_{i}^{t-1}\sum_{k}\sigma_{k}^{t}\end{cases} \tag{16}\]

#### b.2.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

\[m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\text{in}}\left(\frac{m^{t}}{\Delta}, \frac{m^{t}}{\Delta}x_{0}+\sqrt{\frac{m^{t}}{\Delta}}w\right)x_{0}\right] \tag{17}\]

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

_Tilted measure:_ For the **tilted measure** (13) we get

\[\begin{split} f_{\text{in}}(A,B;x_{0})&=\frac{\int \mathrm{d}xPx_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2}x^{2}/2)\exp( Bx-Ax^{2}/2)}{\int\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((A+\gamma^{2})/2)}{\cosh(B+\gamma w+\gamma^{2}x_{0})}}\end{split} \tag{18}\]

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) \tag{19}\]
Pinning measure:In the same way, the **pinning measure** (14) is associated to

\[f_{\text{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 dxx_{P}(x)\exp\left(Bx-Ax^{2}/2\right)}{\int dxP_{X}(x)\exp\left(Bx -Ax^{2}/2\right)}&\text{with probability }1-\theta\end{cases} \tag{15}\] \[=\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}\]

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)\,. \tag{16}\]

#### b.2.4 Phase diagrams

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_{\text{inf}}-\chi_{\text{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.

Autoregressive networks vs FlowsAs we can see from the plots, for this model the tri-critical point for flow-based sampling is at \(\Delta_{\text{tri}}/\rho^{2}\approx 1.08\), while for autoregressive-based sampling is at \(\Delta_{\text{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 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 (12), 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_{\text{tri}}/\rho^{2}\approx 1.08\) for flow-based and \(\Delta_{\text{tri}}/\rho^{2}\approx 1.069\) for autoregressive based sampling.
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.

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 [83, 112]) 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.

### 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:

\[\mathcal{H}(\mathbf{x})=-\frac{\sqrt{3}}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}\,, \tag{101}\]

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, 116, 117, 46].

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}(\mathbf{x})\propto\prod_{i}P_{X}(x_{i})\prod_{i<j<k}e^{\frac{\sqrt{3} }{N}J_{ijk}x_{i}x_{j}x_{k}}\,, \tag{102}\]

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 (100) we will exploit again the fact that considering an equilibrium configuration \(\mathbf{x}_{0}\) is equivalent to take as \(\mathbf{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}(\mathbf{x})=\frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right) e^{\gamma(t)^{2}(\mathbf{x},\mathbf{x}_{0})+\gamma(t)(z,\mathbf{x})-\frac{\gamma(t)^{2} }{2}\|\mathbf{x}\|^{2}}\prod_{i<j<k}e^{\frac{\sqrt{3}}{2}J_{ijk}x_{i}x_{j}x_{k}} \tag{103}\]

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}\)
As before, we can notice that with respect to the original measure in (111), 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 (118) the original problem becomes (denoting as \(S_{\theta}\) the pinned list):

\[P_{\theta}(\mathbf{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}\theta}{\lambda}J_{ijk}x_{i}x_{j}x_{k}} \tag{119}\]

#### 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] :

\[\frac{1}{N}\mathbb{E}_{\mathbf{x}_{0},\mathbf{z},J}\log Z_{\gamma} \xrightarrow[N\to\infty]{}\text{ argmax }\Phi_{RS}(m) \tag{120}\] \[\Phi_{RS}(m) = \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} \tag{121}\]

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** (117) we have

\[Z_{x}(A,B;x_{0}) =\int\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2 }x^{2}/2)\exp(Bx-Ax^{2}/2) \tag{122}\] \[=e^{-(A+\gamma^{2})/2}\cosh(B+\gamma w+\gamma^{2}x_{0})\]

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

\[\Phi_{\mathrm{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}\,. \tag{123}\]

_Pinning measure:_ Meanwhile, considering the **pinning measure** (118) leads to

\[Z_{x}(A,B;x_{0}) =\begin{cases}\int\mathrm{d}xP_{X}(x)\delta_{x,x_{0}}\exp(Bx-Ax^ {2}/2)&\text{with probability }\theta\\ \int\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)&\text{with probability }1-\theta \end{cases} \tag{124}\] \[=\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}\]

which in turn gives

\[\Phi_{\mathrm{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}}\,. \tag{125}\]

#### 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)}[ \mathbf{Y}_{t}]_{i}),\quad\sigma_{i}^{t+1}=\cosh(B_{i}^{t}+\frac{\alpha(t)}{ \beta(t)^{2}}[\mathbf{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}\widehat{\mathbf{x}}^{t}\cdot \widehat{\mathbf{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\end{cases} \tag{121}\]

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

_Pinning measure:_ The AMP equations for the pinning measure (109) 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. (121), provided that we fix \(\gamma=0\).

The resulting equations are:

\[\begin{cases}&\widehat{x}_{i}^{t+1}=\begin{cases}[\mathbf{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}\widehat{\mathbf{x}}^{t} \cdot\widehat{\mathbf{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\end{cases} \tag{122}\]

#### b.3.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 (107). We can use the formalism in [86] to get

\[m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\text{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] \tag{123}\]

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** (108) we have

\[\begin{split} f_{\text{in}}(A,B;x_{0})&=\frac{\int \mathrm{d}xxP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2}x^{2}/2)\exp( Bx-Ax^{2}/2)}{\int\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} \tag{124}\]

which leads to the State Evolution equations

\[\begin{split}\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) \end{split} \tag{125}\]

_Pinning measure:_ In the same way, for the **pinning measure** (109) we get

\[\begin{split} f_{\text{in}}(A,B;x_{0})&=\begin{cases} \frac{P_{X}(x_{0})x_{0}\exp(Bx_{0}-Ax_{0}^{2}/2)}{P_{X}(x_{0})\exp(Bx_{0}-Ax_ {0}^{2}/2)}&\text{with probability }\theta\\ &\frac{\int dx_{P}X(x)\exp(Bx-Ax^{2}/2)}{\int dxP_{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} \tag{126}\]

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) \tag{127}\]
#### b.3.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.

### 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 \(\mathbf{x}\in\mathcal{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:

\[\mathcal{H}(\mathbf{x})=-\frac{\sqrt{3}}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}\,, \tag{100}\]

where \(J_{ijk}\sim\mathcal{N}(0,1)\) and now \(P_{X}(x_{i})=\mathcal{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}(\mathbf{x})\propto\prod_{i}P_{X}(x_{i})\prod_{i<j<k}e^{\frac{\sqrt{3}}{ N}J_{ijk}x_{i}x_{j}x_{k}}\,, \tag{101}\]

where \(P_{X}(x)=\mathcal{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}(\mathbf{x})=\frac{1}{Z_{\gamma}}\left(\prod_{i}P_{X}(x_{i})\right)e ^{\gamma(t)^{2}\langle\mathbf{x},\mathbf{x}_{0}\rangle+\gamma(t)\langle z, \mathbf{x}\rangle-\frac{\gamma(t)^{2}}{2}\|\mathbf{x}\|^{2}}\prod_{i<j<k}e^{ \frac{\sqrt{3}g}{N}J_{ijk}x_{i}x_{j}x_{k}} \tag{102}\]

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

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 (100), 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.
For the **pinned measure** defined in (111) the original problem becomes (denoting as \(S_{\theta}\) the pinned list):

\[P_{\theta}(\mathbf{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{S_{\theta}}}{N}J_{ijk}x_{i}x_{j}x_{k}} \tag{112}\]

#### b.2.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] :

\[\frac{1}{N}\mathbb{E}_{\mathbf{x}_{0},\mathbf{z},J}\log Z_{\gamma} \xrightarrow[N\to\infty]{}\text{ argmax }\Phi_{RS}(m) \tag{113}\] \[\Phi_{RS}(m) = \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} \tag{114}\]

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** (110) we have

\[Z_{x}(A,B;x_{0}) =\int\mathrm{d}xP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2 }x^{2}/2)\exp(Bx-Ax^{2}/2) \tag{115}\] \[=\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)\]

which leads to

\[\Phi_{\text{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}\,. \tag{116}\]

_Pinning measure:_ Meanwhile, considering the **pinning measure** (111) leads to

\[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} \tag{117}\]

which in turn gives

\[\Phi_{\text{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}}\,. \tag{118}\]

#### b.2.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} \widehat{\mathbf{x}}^{t}\cdot\widehat{\mathbf{x}}^{t-1}\\ \widehat{x}_{i}^{t+1}=\frac{B_{i}^{t}+\frac{\alpha(t)}{\beta\sigma(2)}\left( \mathbf{Y}_{t}\right)_{i}}{3\beta^{2}\|\mathbf{x}^{2}\|^{2}_{2}/(2N)+\gamma^{2 }+1}\,,\sigma^{t+1}=\frac{1}{3\beta^{2}\|\mathbf{x}^{4}\|^{2}_{2}/(2N)+\gamma^ {2}+1}\end{cases} \tag{119}\]

where \(\alpha(t)\) and \(\beta(t)\) are the functions defining the interpolant process, fixed at the start, and \(\mathbf{Y}_{t}\) is the value of the noisy observation at time \(t\).
Pinning measure:When considering the pinning measure (144), 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. (130), provided we fix \(\gamma=0\).

The resulting algorithm is the following:

\[\left\{\begin{array}{ll}\widehat{x}_{i}^{t+1}=\begin{cases}[ \mathbf{x}_{0}]_{i}&\text{if }i\in S_{\theta}\\ \frac{B_{i}^{t}}{3\beta^{2}\|\mathbf{x}^{t}\|^{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}\|\mathbf{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}\widehat{\mathbf{x}}^{t} \cdot\widehat{\mathbf{x}}^{t-1}\sum_{k}\sigma_{k}^{t}/N\end{cases}\right. \tag{131}\]

#### 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. (130). Here, again we follow the results presented in [86], such that as the Ising version of the model, we have

\[m^{t+1}=\mathbb{E}_{x_{0},w}\left[f_{\text{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] \tag{132}\]

where \(x_{0}\sim P_{X}\), \(w\sim\mathcal{N}(0,1)\) and \(f_{\text{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** (144) we have

\[\begin{split} f_{\text{in}}(A,B;x_{0})&=\frac{\int \mathrm{d}xxP_{X}(x)\exp(\gamma^{2}xx_{0}+\gamma wx-\gamma^{2}x^{2}/2)\exp( Bx-Ax^{2}/2)}{\int\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} \tag{133}\]

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}\,. \tag{134}\]

Pinning measure:In the same way, for the **pinning measure** (144) we find

\[\begin{split} f_{\text{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{\mathrm{d}xP_{X}(x)\exp(Bx-Ax^{2}/2)}{\int\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} \tag{135}\]

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}} \tag{136}\]

#### 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.
### 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 \(\mathbf{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 \(\mathbf{x}_{\partial a}\). In terms of probability distribution, this translates to

\[P_{0}(\mathbf{x})=\frac{1}{Z(G)}\prod_{a=1}^{M}\omega(\mathbf{x}_{\partial a} )\,,\quad\omega(x_{1},\ldots,x_{k})=\begin{cases}0&\text{if }\sum_{i=1}^{k}x_{i}=\pm k\\ 1&\text{otherwise}\end{cases} \tag{100}\]

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\to a}=f(\{u_{b\to i}\}_{b\in\partial_{i}\setminus a})\,,\quad u_{a\to i}= g(\{h_{j\to a}\}_{j\in\partial_{a}\setminus i}) \tag{101}\]

where

\[f(u_{1},\ldots,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})} \tag{102}\]

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 (101), 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_{\text{tri}}=2/3\) for flow-based and \(T_{\text{tri}}=\sqrt{1/2}\) for autoregressive based sampling.
is the function defining the messages going from factor nodes to variable nodes and

\[g(h_{1},\ldots,h_{k-1})=\frac{\sum_{x_{1},\ldots,x_{k}}\omega(x_{1},\ldots,x_{k})x _{k}\prod_{i=1}^{k-1}(1+h_{i}x_{i})}{\sum_{x_{1},\ldots,x_{k}}\omega(x_{1}, \ldots,x_{k})\prod_{i=1}^{k-1}(1+h_{i}x_{i})} \tag{100}\]

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

\[\frac{1}{N}\ln Z(G)=\frac{1}{N}\sum_{i=1}^{N}\ln\mathcal{Z}_{0}^{\mathrm{v}}( \{u_{a\to i}\}_{a\in\partial i})+\frac{1}{N}\sum_{a=1}^{M}\ln\mathcal{Z}_{0}^ {\mathrm{c}}(\{h_{i\to a}\}_{i\in\partial a})-\frac{1}{N}\sum_{(i,a)}\ln \mathcal{Z}_{0}^{\mathrm{e}}(h_{i\to a},u_{a\to i})\, \tag{101}\]

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

\[\mathcal{Z}_{0}^{\mathrm{v}}(u_{1},\ldots,u_{d}) =\sum_{x}\prod_{i=1}^{d}\left(\frac{1+xu_{i}}{2}\right)\, \tag{102}\] \[\mathcal{Z}_{0}^{\mathrm{c}}(h_{1},\ldots,h_{k}) =\sum_{x_{1},\ldots,x_{k}}\omega(x_{1},\ldots,x_{k})\prod_{i=1}^{k }\left(\frac{1+x_{i}h_{i}}{2}\right)\,\] (103) \[\mathcal{Z}_{0}^{\mathrm{e}}(h,u) =\sum_{x}\left(\frac{1+xh}{2}\right)\left(\frac{1+xu}{2}\right). \tag{104}\]

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:

\[\mathcal{P}^{RS}(h) =\sum_{d=0}^{\infty}p_{d}\int\left(\prod_{i=1}^{d}\mathrm{d}u_{i }\widehat{\mathcal{P}}^{RS}(u_{i})\right)\,\delta(h-f(u_{1},\ldots,u_{d}))\, \tag{105}\] \[\widehat{\mathcal{P}}^{RS}(u) =\int\left(\prod_{i=1}^{k-1}\mathrm{d}h_{i}\mathcal{P}^{RS}(h_{i })\right)\,\delta(u-g(h_{1},\ldots,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.3.3 1RSB and Tree reconstruction equations:

The 1RSB cavity method aims to compute the potential

\[\Phi_{1}(m)=\lim_{N\to\infty}\frac{1}{N}\log\left(\sum_{\gamma}Z_{\gamma}^{m}\right) \tag{106}\]

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

\[Q_{+}^{(t+1)}(h) =\sum_{d=0}^{\infty}p_{d}\int\left(\prod_{i=1}^{d}\mathrm{d}u_{i} \widehat{Q}_{+}^{(t)}(u_{i})\right)\,\delta(h-f(u_{1},\ldots,u_{d}))\, \tag{111}\] \[\widehat{Q}_{+}^{(t)}(u) =\sum_{x_{1},\ldots,x_{k-1}}\tilde{p}(x_{1},...,x_{k-1}|+)\int \left(\prod_{i=1}^{k-1}\mathrm{d}h_{i}Q_{+}^{(t)}(h_{i})\right)\delta(u-g(x_{1} h_{1},\ldots,x_{k-1}h_{k-1}))\,\]

where

\[\tilde{p}(x_{1},\ldots,x_{k-1}|+)=\frac{\omega(x_{1},\ldots,x_{k-1},+)}{\sum \limits_{x_{1}^{\prime},\ldots,x_{k-1}^{\prime}}\omega(x_{1}^{\prime},\ldots, x_{k-1}^{\prime},+)}=\frac{\sum\limits_{p=0}^{k-1}\omega_{p}\operatorname{ \mathbb{I}}\left[\sum\limits_{i=1}^{k-1}x_{i}=k-1-2p\right]}{\sum\limits_{p=0 }^{k-1}\binom{k-1}{p}\omega_{p}}. \tag{112}\]

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

\[\begin{split}\Phi_{\mathrm{RS}}=&\sum_{d=0}^{\infty }p_{d}\int\left(\prod_{i=1}^{d}\mathrm{d}u_{i}\widehat{Q}_{+}(u_{i})\right) \ln\mathcal{Z}_{0}^{\mathrm{v}}(u_{1},\ldots,u_{d})+\alpha\int\left(\prod_{i=1 }^{k}\mathrm{d}h_{i}Q_{+}(h_{i})\right)\ln\mathcal{Z}_{0}^{\mathrm{c}}(h_{1}, \ldots,h_{k})\\ &-\alpha k\int\mathrm{d}h\,\mathrm{d}uQ_{+}(h)\widehat{Q}_{+}(u )\ln\mathcal{Z}_{0}^{\mathrm{e}}(h,u)\,\end{split} \tag{113}\]

where \(\mathcal{Z}_{0}^{\mathrm{c}}\), \(\mathcal{Z}_{0}^{\mathrm{e}}\) and \(\mathcal{Z}_{0}^{\mathrm{v}}\) are the local partitions reported in Eq. (110), (112) and (111) respectively.

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

\[\chi=\int\mathrm{d}hQ_{+}(h)h^{2}\,. \tag{114}\]

#### b.2.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^{\mathrm{tilt}}(u_{1},\ldots,u_{d})=\frac{e^{2\gamma(\gamma+z)}\prod\limits _{i=1}^{d}(1+u_{i})-\prod\limits_{i=1}^{d}(1-u_{i})}{e^{2\gamma(\gamma+z)} \prod\limits_{i=1}^{d}(1+u_{i})+\prod\limits_{i=1}^{d}(1-u_{i})}\,,\quad z \sim\mathcal{N}(0,1)\,, \tag{115}\]

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

\[\mathcal{Z}_{0}^{\mathrm{v}}(u_{1},\ldots,u_{d})=e^{\gamma(\gamma+z)}\prod \limits_{i=1}^{d}\left(\frac{1+u_{i}}{2}\right)+e^{-\gamma(\gamma+z)}\prod \limits_{i=1}^{d}\left(\frac{1-u_{i}}{2}\right)\,,\quad z\sim\mathcal{N}(0,1)\,. \tag{116}\]

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^{\mathrm{pinn}}(u_{1},\ldots,u_{d})=\begin{cases}1&\text{with probability $\theta$}\\ \frac{\prod\limits_{i=1}^{d}(1+u_{i})-\prod\limits_{i=1}^{d}(1-u_{i})}{\prod \limits_{i=1}^{d}(1+u_{i})+\prod\limits_{i=1}^{d}(1-u_{i})}&\text{with probability $1-\theta$}\,\end{cases} \tag{117}\]
and again the only contribution to the partition function that is modified is

\[\mathcal{Z}_{0}^{\mathrm{v}}(u_{1},\ldots,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}\,. \tag{101}\]

Phase diagrams:In Fig. 10 we present the phase diagrams for the \(k\)-NAESAT problem, considering the case \(k=5\). Compared to the plots reported in the main text, here we display directly the difference \(\chi_{\mathrm{inf}}-\chi_{\mathrm{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_{\mathrm{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. 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_{\mathrm{tri}}\), where we expect the sampling scheme to work, to a value of \(T\) in the interval \([T_{d},T_{\mathrm{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

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 (100), 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_{\mathrm{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.
spherical \(p\)-spin. We do this by analysing the following observable:

\[\text{OV}(\gamma)\equiv\frac{1}{N}\mathbb{E}\left[\mathbf{Y_{t}}\langle\mathbf{x }\rangle\right] \tag{104}\]

which is the overlap between the observation vector \(\mathbf{Y_{t}}\) and the average magnetization \(\langle\mathbf{x}\rangle\).

Indeed, by definition of the observation process we can write

\[\text{OV}(\gamma)=\frac{\alpha(t)}{N}\mathbb{E}\left[\mathbf{x_{0}}\langle \mathbf{x}\rangle\right]+\frac{\beta(t)}{N}\mathbb{E}\left[\mathbf{z}\langle \mathbf{x}\rangle\right] \tag{105}\]

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

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

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

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_{\text{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_{\text{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. 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.

## 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\left|Y\right)\right)\). Denoting \(\left\langle\cdot\right\rangle\) the "Boltzmann" expectation, that is the average with respect to the \(P\left(X=\cdot\left|Y\right)\right)\), and \(\mathbb{E}\left[\cdot\right]\) the "Disorder" expectation, that is with respect to \(\left(X^{*},Y\right)\). 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] \tag{12}\]

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

Figure 12: **Spherical p-spin: Checking the Nishimori conditions for flow-based sampling.** We compare the results for the overlap OV, defined in (12) 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.
function \(g\) as

\[\left\langle g(X^{(1)},\ldots,X^{(k)})\right\rangle_{k}=:\int\prod_{i=1}^{k} \mathrm{d}x_{i}P(x_{i}|Y)g(X^{(1)},\ldots,X^{(k)}) \tag{104}\]

we have, starting from the right-hand side

\[\mathbb{E}_{Y,X^{*}}\left[\left\langle f\left(Y,X^{(1)},\ldots,X^ {(k-1)},X^{*}\right)\right\rangle_{k-1}\right]\] \[=\int dx^{*}dyP(x^{*}|Y)P(Y)\left\langle f\left(Y,X^{(1)},\ldots,X ^{(k-1)},X^{*}\right)\right\rangle_{k-1}\] \[=\mathbb{E}_{Y}\int\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}\] \[=\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}(\mathbf{x})=\frac{1}{Z_{n}}\exp\left(\gamma^{2}\langle\mathbf{x_{0 }},\mathbf{x}\rangle+\gamma\langle\mathbf{z},\mathbf{x}\rangle-\frac{\gamma^ {2}}{2}\|\mathbf{x}\|^{2}\right)P_{0}(\mathbf{x}) \tag{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)] \tag{106}\]

_then_

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

_and_

\[\partial_{\gamma}^{2}f_{N}=N\gamma^{2}\mathbb{E}\left[\left(\langle(\frac{ \mathbf{x}(\gamma)\cdot\mathbf{x}_{0}}{N})^{2}\rangle\right)-\left(\frac{ \langle\mathbf{x}(\gamma)\rangle\mathbf{x}_{0}}{N}\right)^{2}\right]=N\gamma^ {2}\mathbb{E}[\mathrm{var}\left(\frac{\mathbf{x}\cdot\mathbf{x}_{0}}{N}\right)] \tag{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:

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

where we use Stein lemma on line 110 and Nishimori on line 110. 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}[\mathrm{var}\left(\frac{\mathbf{x}\cdot\mathbf{x}_{0}}{N}\right)] \rightarrow_{N\rightarrow\infty}0\]

Proof.: From the derivative, we have that

\[\int_{\gamma_{1}}^{\gamma_{2}}\mathrm{var}_{\gamma}\left(\frac{\mathbf{x} \cdot\mathbf{x}_{0}}{N}\right)=\frac{1}{\gamma^{2}N}(\mu(\gamma_{2})-\mu( \gamma_{1}))\leq\frac{K}{N} \tag{101}\]

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(\mathbf{x}|\mathbf{y}(t)=\mathbf{Y}_{t})=\frac{1}{Z(\mathbf{Y}_{t})}\exp \left(\frac{\alpha(t)}{\beta(t)^{2}}\langle\mathbf{Y}_{t},\mathbf{x}\rangle- \frac{\alpha(t)^{2}}{2\beta(t)^{2}}||\mathbf{x}||^{2}\right)P_{0}(\mathbf{x}). \tag{102}\]

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

### 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 interpreted as an instance of the Fluctuation-dissipation theorem:

**Lemma 2**.: _Suppose that the measure \(P_{0}\) has compact support. Then the Jacobian of the vector \(\frac{\partial b(\mathbf{y},t)}{\partial\mathbf{y}}\) is related to the covariance of the tilted measure \(P_{t,\mathbf{Y}_{t}}=P(\mathbf{x}|\mathbf{y}(t)=\mathbf{Y}_{t})\) as follows:_

\[\frac{\partial b(\mathbf{y},t)}{\partial\mathbf{y}}|_{\mathbf{y}=\mathbf{Y}_{t }}=\frac{\alpha(t)}{\beta(t)^{2}}(\dot{\alpha}(t)-\frac{\dot{\beta}(t)\alpha( t)}{\beta(t)})\,\mathrm{Covar}[P_{t,\mathbf{Y}_{t}}]+\frac{\dot{\beta}(t)}{ \beta(t)} \tag{103}\]

Proof.: We have, from Eq. (5):

\[b(\mathbf{y},t) =\mathbb{E}[\dot{\alpha}(t)\mathbf{x}_{0}+\dot{\beta}(t)z|\mathbf{ y}(t)=\mathbf{y}] \tag{104}\] \[=\mathbb{E}[\dot{\alpha}(t)\mathbf{x}_{0}+\frac{\dot{\beta}(t)}{ \beta(t)}(\mathbf{y}(t)-\alpha(t)\mathbf{x}_{0})|\mathbf{y}(t)=\mathbf{y}]\] \[=(\dot{\alpha}(t)-\frac{\dot{\beta}(t)\alpha(t)}{\beta(t)})\mathbb{ E}[\mathbf{x}_{0}|\mathbf{y}(t)=\mathbf{y}]+\frac{\dot{\beta}(t)}{\beta(t)} \mathbf{y}.\]
where we used the relation in Eq. (1). Next, we have through the expression for the posterior measure \(P_{t,\mathbf{Y}_{t}}\):

\[\mathbb{E}[\mathbf{x_{0}}|\mathbf{y}(t)=\mathbf{y}]=\mathbb{E}_{P_{0}}\left[ \frac{1}{Z(\mathbf{y})}\mathbf{x}\exp\left(\frac{\alpha(t)}{\beta(t)^{2}} \langle\mathbf{y},\mathbf{x}\rangle-\frac{\alpha(t)^{2}}{2\beta(t)^{2}}|| \mathbf{x}||^{2}\right)\right]. \tag{100}\]

Using the boundedness of the support of \(P_{0}\) and dominated convergence theorem, we may differentiate the L.H.S inside the expectation. We obtain that for all \(t\in(0,1]\), \(\mathbb{E}[\mathbf{x_{0}}|\mathbf{y}(t)=\mathbf{Y}_{t}]\) is differentiable w.r.t \(\mathbf{Y}_{t}\) with the Jacobian given by \(\frac{\alpha(t)}{\beta(t)^{2}}\operatorname{Covar}[P_{t,\mathbf{Y}_{t}}]\). Substituting into Eq. (100) completes the proof.

### Derivation of the ODE (Eq. (6))

We start by presenting an informal derivation of the continuity equation, borrowed from [124]. Recall that, by definition, \(\mathbf{y}(t)=\alpha(t)\mathbf{x}_{0}+\beta(t)\mathbf{z})\). The density \(\rho(\mathbf{y},t)\) can be expressed as:

\[\rho(\mathbf{y},t)=\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}}\delta(\mathbf{y }-\mathbf{y}(\mathbf{t}))\rho_{\mathbf{0}}(\mathbf{x_{0}})\rho_{\gamma}( \mathbf{z})\mathbf{dx_{0}dz}, \tag{101}\]

where \(\rho_{\gamma}\) denotes the density of the standard Gaussian measure \(\gamma=\mathcal{N}(\mathbf{0},\mathbb{I}_{N})\). Differentiating both sides of Equation (101) yields:

\[\begin{split}\partial_{t}\rho(\mathbf{y},t)&=-\int _{\mathbb{R}^{N}\times\mathbb{R}^{N}}\nabla\delta(\mathbf{y}-\mathbf{y}( \mathbf{t}))\cdot\partial_{t}\mathbf{y}(\mathbf{t})\rho_{\mathbf{0}}( \mathbf{x_{0}})\rho_{\gamma}(\mathbf{z})\mathbf{dx_{0}dz}\\ &=-\nabla\cdot\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}}\delta( \mathbf{y}-\mathbf{y}(\mathbf{t}))\partial_{t}\mathbf{y}(\mathbf{t})\rho_{ \mathbf{0}}(\mathbf{x_{0}})\rho_{\gamma}(\mathbf{z})\mathbf{dx_{0}dz}\\ &=-\nabla\cdot\left(\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}} \left(\delta(\mathbf{y}-\mathbf{y}(t))\partial_{t}\mathbf{y}(t)\frac{\rho_{ \mathbf{0}}(\mathbf{x_{0}})\rho_{\gamma}(\mathbf{z})}{\rho(\mathbf{y},t)}d \mathbf{x_{0}dz}\right)\rho(\mathbf{y},t)\right)\end{split} \tag{102}\]

Eq. (4) is then obtained by noticing that equals \(b(\mathbf{y},t)\) defined by (5).

We refer to [124, 41] for the complete derivation based on the above approach.

We next prove that the pushforward measure obtained after applying the flow defined by the ODE (6) to initial Gaussian noise \(\mathbf{z}\) results in a sample from the target measure \(P_{0}\). For the sake of completeness, we include an alternative derivation for Eq. (4) in our proof, based on the weak form of the continuity equation. Our result allows \(P_{0}\) to be a discrete measure, by restricting the time to \(t\in(0,1)\).

**Lemma 3**.: _Let \(\epsilon\) be an arbitrarily-fixed real in \((0,1)\). Let \(\mathbf{Y_{t}}(\mathbf{z})\) denote the flow associated with the ODE in Eq. (6) starting from \(t=1\) to \(t=\epsilon\) with \(\mathbf{z}\sim\gamma\), where \(\gamma\) denotes the standard Gaussian measure \(\gamma=\mathcal{N}(\mathbf{0},\mathbb{I}_{N})\). Suppose that \(P_{0}\) has bounded support. Then, the pushforward measure \(\mathbf{Y_{t}}_{\#}\gamma\) at any time \(t\in(0,1]\) equals the measure corresponding to the law \(P_{t}\) of the interpolant \(\mathbf{y}(t)\) defined by equation (1)._

Proof.: Let \(\psi\in C_{c}^{\infty}(\mathbb{R}^{d})\) be an arbitrary test function. Using the change of variables formula, the expectation of \(\psi\) w.r.t the measure \(\mu_{t}\) at time \(t\) can be expressed as:

\[\int_{\mathbb{R}^{N}}\psi(\mathbf{y})dP_{t}(\mathbf{y})=\int_{\mathbb{R}^{N} \times\mathbb{R}^{N}}\psi(\mathbf{y}(t))dP_{0}(\mathbf{x}_{0})d\gamma_{z}( \mathbf{z}), \tag{103}\]

where \(\mathbf{y}(t)\) is defined as a measurable function of \(\mathbf{x}_{0},\mathbf{z}\) through Eq. (1).

Therefore, \(P_{t}\) evolves in a distributional sense as follows:

\[\int_{\mathbb{R}^{N}}\psi(\mathbf{y})\partial_{t}dP_{t}(\mathbf{y})=\frac{d \int_{\mathbb{R}^{N}}\psi(\mathbf{y})dP_{t}(\mathbf{y})}{dt}=\int_{\mathbb{R}^ {N}\times\mathbb{R}^{N}}\nabla\psi(\mathbf{y}(t))\cdot\partial_{t}\mathbf{y}( t)dP_{0}(\mathbf{x}_{0})d\gamma_{z}(\mathbf{z}). \tag{104}\]

Recall that:

\[b(\mathbf{y},t)=\mathbb{E}[\partial_{t}\mathbf{y}(t)|\mathbf{y}(t)=\mathbf{y}]. \tag{105}\]
Using the definition of the conditional expectation, and the change of variables formula, we have:

\[\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}}\nabla\psi(\mathbf{y}(t)) \cdot\partial_{t}\mathbf{y}(t)dP_{0}(\mathbf{x}_{0})d\gamma_{z}(\mathbf{z}) =\int_{\mathbb{R}^{N}}\nabla\psi(\mathbf{y}).\mathbb{E}[\partial_{ t}\mathbf{y}(t)|\mathbf{y}(t)=\mathbf{y}]dP_{t}(\mathbf{y})\] \[=\int_{\mathbb{R}^{N}}\nabla\psi(\mathbf{y})\cdot b(\mathbf{y},t )dP_{t}(\mathbf{y})\] \[=-\int_{\mathbb{R}^{N}}\psi(\mathbf{y})\nabla.(b(\mathbf{y},t)P_ {t}(\mathbf{y}))d\mathbf{y},\]

where we used the compactness of the support of \(\psi(x)\) and the distributional definition of the divergence operator \(\nabla\). Substituting in Eq. (100), we obtain:

\[\int_{\mathbb{R}^{N}}\psi(\mathbf{y})\partial_{t}dP_{t}(\mathbf{y})=-\int_{ \mathbb{R}^{N}}\psi(\mathbf{y})\nabla.(b(\mathbf{y},t)P_{t}(\mathbf{y}))d \mathbf{y}. \tag{101}\]

Furthermore, for all \(t>\epsilon\) for fixed \(\epsilon>0\), we have \(\beta(t)>0\). Then, Lemma 2 implies that \(b(\mathbf{y},t)\) is Lipschitz w.r.t \(\mathbf{y}\). Thus, the probability flow \(\mathbf{Y_{t}}\) for \(t\in[\epsilon,1]\) exists and is unique. Since both the push-forward measure \(\mathbf{Y_{t}}\#\gamma\) and the law of \(\mathbf{y}(t)\) i.e. \(P_{t}\) satisfy the continuity equation with velocity field \(b(\mathbf{y},t)\) (Lemma 4.1.1. in [125]). By the uniqueness of the solution of the continuity equation (see e.g. [126]), \(b(t,x)_{\#}\mu(z)=\mu_{t}\). 

### Sampling Guarantees

In this section, we quantify the effect of discretization errors in the velocity field, leading to sampling guarantees for Algorithm 1. Again, we fix a parameter \(\epsilon\) lying in \((0,1)\) We rely on the following assumptions:

**Assumption 1**.: _The measure \(P_{0}\) has compact support lying in a sphere with radius bounded as \(\mathcal{O}(\sqrt{N})\)_

\[\mathrm{Supp}(P_{0})\subseteq\mathrm{S}^{N-1}(\sqrt{N}R), \tag{102}\]

_for some \(R\) independent of \(N\)._

**Assumption 2**.: _The spectral norm of the vector field's Jacobian \(\frac{\partial b(\mathbf{Y}_{t},t)}{\partial\mathbf{Y}_{t}}\) w.r.t \(\mathbf{Y}\) is uniformly bounded, i.e. \(\|\frac{\partial b(\mathbf{Y},t)}{\partial\mathbf{Y}}\|_{2}\leq L,\forall t \in(\epsilon,1),\forall\mathbf{Y}\in\mathbb{R}^{d}\) for some \(L\geq 0\)._

**Assumption 3**.: \(\|\frac{\partial b(\mathbf{Y},t)}{\partial t}\|\) _is uniformly bounded by \(M\sqrt{N}\) in \(\mathbf{Y},t\) for \(t\in(\epsilon,1)\), for some \(M\geq 0\)._

We have the following error bounds for the discretization error associated to the flow:

**Lemma 4**.: _Consider the ODE defined by Eq. (6) i.e \(\frac{d\mathbf{Y}}{dt}=b(\mathbf{Y},t)\), with \(\mathbf{Y}\in\mathbb{R}^{N}\). Let \(\mathbf{Y}_{t}(\mathbf{z})\) denote the flow associated to the above ODE at time \(t\in(0,1]\) starting from some fixed \(\mathbf{z}\in\mathbb{R}^{N}\) at time \(t=1\). Let \(\mathbf{Y}_{\delta,t}(\mathbf{z})\) denote the iterates of the forward Euler method applied to the above ODE (in reverse time) starting from the same initialization \(\mathbf{z}\) with step-size \(\delta\) i.e, for \(i\in\mathbb{N}\):_

\[\mathbf{Y}_{\delta,\delta(i-1)}(\mathbf{z})=\mathbf{Y}_{\delta,\delta i}( \mathbf{z})-\delta b(\mathbf{Y}_{\delta,\delta i}(\mathbf{z}),\delta i), \tag{103}\]

_with \(\mathbf{z}\) fixed. Under Assumptions 1,2,3, there exists a constant \(A(\epsilon)\) such that, for all \(k\in\mathbb{N},\mathbf{z}\in\mathbb{R}^{N}\) with \(k\delta\geq\epsilon\):_

\[\|\mathbf{Y}_{\delta,k\delta}(\mathbf{z})-\mathbf{Y}_{k\delta}(\mathbf{z})\|_ {2}\leq(\frac{M+AL}{L})\sqrt{N}e^{Lk\delta}\delta, \tag{104}\]

for small enough \(\delta\).

Proof.: We first note that Assumption 1 and Lemma 2 imply that \(b(\mathbf{Y},t)\) is uniformly bounded for any fixed \(\mathbf{z}\).Then, the above bound follows from standard analysis of the forward Euler method. See for example Chapter 7 in [127]. 

**Remark**: The Lipschitzness of the posterior mean (optimal denoiser) w.r.t \(\mathbf{Y}_{t}\) is expected to hold for the setups considered in our work in light of similar results proven in [36, 42]. Furthermore, as Lemma 2 shows, it is equivalent to the boundedness of the covariance of the tilted measure. Similarly, Lemma 1 provides control over \(\frac{\partial b(\mathbf{Y},t)}{\partial t}\) through the variance of the overlap.

We now prove that the proposed algorithm with sufficiently small step for the discretized ODE produces a distribution close to the target distribution in Wasserstein distance. The obtained bounds on the error between the ODE and the algorithm iterates can be related to the Wasserstein distance between the corresponding pushforward measures.
**Theorem 3**.: _Let \(P_{\mathrm{alg,N_{steps}}}\) denote the measure of the output produced by Algorithm 1 after \(N_{\mathrm{steps}}\). Suppose that the vector field, for any \(\eta>0\), there exist \(N_{\mathrm{steps}}\), independent of \(N\), such that the normalized Wasserstein distance \(\frac{1}{\sqrt{N}}W_{2}(P_{\mathrm{alg,N_{steps}}},P_{0})<\eta\)._

Proof.: Let \(0<\epsilon<h\) be arbitrary. By definition, \(\mathbf{Y}_{h,\epsilon}\) has law \(P_{\mathrm{alg,N_{steps}}}\) while, from Lemma 3, \(\mathbf{Y}_{\epsilon}\) equals in law \(\mathbf{y}(\epsilon)\). Setting the same initialization \(\mathbf{z}\) induces a coupling between the two measures. Recall that the 2-Wasserstein distance [125, 128] between two measures \(\mu,\nu\) on \(\mathbb{R}^{N}\) is defined as \(W_{2}^{2}(\mu,\nu)=\inf_{\Gamma(\mu,\nu)}\mathbb{E}\|\mathbf{X}-\mathbf{Y}\|^ {2}\), where \(\Gamma(\mu,\nu)\) denotes the set of couplings between \(\mu,\nu\) with \(\mathbf{X},\mathbf{Y}\) being distributed as \(\mu,\nu\) respectively. Therefore, we obtain:

\[W_{2}^{2}(P_{\mathrm{alg,N_{steps}}},P_{\epsilon})\leq\mathbb{E}\|\mathbf{Y} _{\delta,\epsilon}(\mathbf{z})-\mathbf{Y}_{\epsilon}(\mathbf{z})\|^{2}. \tag{120}\]

Lemma 4 then implies that by choosing a small enough step-size \(\delta\), or equivalently, with a large enough \(N_{\mathrm{steps}}\), \(W_{2}(P_{\mathrm{alg,N_{steps}}},P_{\epsilon})\) can be made arbitrarily small for any fixed epsilon.

We further have that \(W_{2}(P_{0},P_{\epsilon})\to 0\) as \(\epsilon\to 0\). By triangle inequality for \(W_{2}\)[125, 128], we obtain:

\[W_{2}(P_{0},P_{\epsilon})\leq W_{2}(P_{\mathrm{alg,N_{steps}}},P_{\epsilon})+ W_{2}(P_{0},P_{\epsilon}). \tag{121}\]

Now, first pick \(\epsilon>0\) such that \(\frac{1}{\sqrt{N}}W_{2}(P_{0},P_{\epsilon})\leq\eta/2\). Subsequently, we pick \(N_{\mathrm{steps}}\) such that \(\frac{1}{\sqrt{N}}W_{2}(P_{\mathrm{alg,N_{steps}}},P_{\epsilon})\leq\eta/2\) to complete the proof. 

The above result establishes that Algorithm 1 samples from a distribution approximating the target distribution up to any desired accuracy, with a finite number of forward Euler steps (\(N_{\mathrm{steps}}\)), independent of the dimension \(N\).

**Remark:** In the presence of a first-order phase transition during the interpolation path \(\|\frac{\partial\mathbf{b}(\mathbf{Y},t)}{\partial t}\|\) may grow as \(\omega(\sqrt{N})\). This is apparent from Lemma 1 which relates \(\frac{\partial\mathbf{b}(\mathbf{Y},t)}{\partial t}\) to the variance of the overlap. This corresponds to \(M\) in Lemma 4 growing with \(N\). Lemma 4 reveals that as long as this growth is polynomial in \(N\), the discretization error can still be controlled by choosing \(h\) to be \(1/\mathrm{poly}(N)\). This would lead only to an additional polynomial time factor in the time complexity of the algorithm. Therefore, with an oracle access to the denoiser, Algorithm 1 might be able to efficiently sample even in the "inefficient" phase. We believe this to be an interesting question for future research.
