††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.††thanks: These authors contributed
equally to this work.

# The planted XY model: thermodynamics and inference

Siyu Chen  École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland  Guanhao Huang  École Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland  Giovanni Piccioli  École Polytechnique Fédérale
de Lausanne (EPFL), Statistical Physics of Computation laboratory, CH-1015
Lausanne, Switzerland  Lenka Zdeborová  École Polytechnique Fédérale de
Lausanne (EPFL), Statistical Physics of Computation laboratory, CH-1015
Lausanne, Switzerland

###### Abstract

In this paper we study a fully connected planted spin glass named the planted
XY model. Motivation for studying this system comes both from the spin glass
field and the one of statistical inference where it models the angular
synchronization problem [1]. We derive the replica symmetric (RS) phase
diagram in the temperature, ferromagnetic bias plane using the approximate
message passing (AMP) algorithm and its state evolution (SE). While the RS
predictions are exact on the Nishimori line (i.e. when the temperature is
matched to the ferromagnetic bias), they become inaccurate when the parameters
are mismatched, giving rise to a spin glass phase where AMP is not able to
converge. To overcome the defects of the RS approximation we carry out a one-
step replica symmetry breaking (1RSB) analysis based on the approximate survey
propagation (ASP) algorithm [2]. Exploiting the state evolution of ASP, we
count the number of metastable states in the measure, derive the 1RSB free
entropy and find the behavior of the Parisi parameter throughout the spin
glass phase.

##  I Introduction

In statistical physics, XY models describe a system in which the spins are
unit norm vectors living in $\mathbb{R}^{2}$, or equivalently unit norm
complex numbers $x_{i}=e^{i\theta_{i}}$. They are also characterized by an
Hamiltonian which is invariant under global phase shifts, thus displaying a
symmetry identifiable with $U(1)$. Early studies have focused on ferromagnetic
lattice models with short range interactions, described by the Hamiltonian
$H_{XY}(\theta)=J\sum_{\langle i,j\rangle}\cos(\theta_{i}-\theta_{j})$. This
line of work culminated with the discovery of the Kosterlitz-Thouless
transition [3]. Disordered models have also received attention [4]. The most
studied way of introducing disorder is through Gaussian couplings, giving a
Hamiltonian $H_{DXY}(\theta)=\sum_{i<j}J_{ij}\cos(\theta_{i}-\theta_{j})$ [4].
The model studied in this paper is a different version of disordered XY model
motivated by the statistical inference problem of angular synchronization,
which consists of retrieving a vector of angles $\theta_{i}\in[0,2\pi]$ from
measurements of their offsets $\theta_{i}-\theta_{j}$. The full definition is
given in section II. This problem arises in many applications, for example
time synchronization in distributed networks [5][6], alignment in signal
processing [7], computer vision [8] and optics [9]. Angular synchronization
was first introduced in [1] and solved using spectral algorithms and
semidefinite relaxation techniques [10]. In [11], the authors found the
replica symmetric solution via a replica and cavity computation. From the
algorithmic point of view, a version of approximate message passing (AMP)[12]
was developed for a general class of problems including the angular
synchronization [13]. Formulating the angular synchronization problem in the
language of XY Hamiltonians, see eq. (4), one obtains a planted version of the
XY model. In planted systems, the couplings between spins depend on a special
configuration $\theta^{\star}$, and the Gibbs measure is nothing but the
conditional distribution of $\theta^{*}$ given the couplings [14]. The
inference problem then corresponds to recovery of $\theta^{*}$ from the
couplings. Another variant of the XY model studied in the literature which
admits a planted interpretation is the Gauge glass, an XY spin glass with
Hamiltonian $H_{GG}=J\sum_{i,j\in E}\cos(\theta_{i}-\theta_{j}-\phi_{ij})$,
where $E$ is the set of edge interactions. The randomness is contained in the
angles $\\{\phi_{ij}\\}$ and possibly in $E$. This model has been first
studied in the random setting with $\phi$ drawn i.i.d. from a uniform
distribution in [15] and later generalized to the planted case, where
$\\{\phi_{ij}\\}$ are drawn from a zero mean von Mises distribution [16].
Previous works on the gauge glass have mainly studied the model on the
Nishimori line [17] [18] i.e. a line in the temperature-coupling diagram,
where calculations greatly simplify. This model was further studied in its
discretized version in [19] for a mixture distribution that interpolates
between ferromagnetic and uniformly distributed couplings and with
interactions on a sparse random graph. Finally, the short range gauge glass
model has also been extended to the quantum setting in [20] and has further
physical relevance [21]. The model considered in this work, Hamiltonian (4),
is closely related to $H_{GG}$. Our choice fell on model defined by eq. (4)
because of the connection to angular synchronization. Thanks to the invariance
of the Hamiltonian under a joint transformation of the couplings and the
spins, we are able to map the partial recovery transition in the inference
case, into an order-disorder transition. Moreover, for appropriate choices of
the parameters (specifically for $\lambda<1$, defined below), the model has a
fully random behavior, thus connecting with the literature on spin glasses.

Our work draws a bridge between the angular synchronization and the studies of
disordered XY models in the statistical physics literature. We consider the
fully connected disordered model associated with angular synchronization, and
we investigate the properties of the Gibbs measure given by the posterior. In
the inference case, one is usually interested in characterizing the error in
the retrieval of the signal. Instead, we focus on the phase diagram spanned by
the inverse temperature and the ferromagnetic bias. Specifically, we
concentrate on the region outside the Nishimori line, which is not studied in
related models [16]. Furthermore, we go beyond the RS analysis of [13][11] and
perform a one step replica symmetry breaking study. Our theoretical analysis
is complemented by algorithms which provide an instrumental way to study
single instances of the model. We use this as an opportunity to characterize
the behaviour of AMP, and its 1RSB version Approximate Survey Propagation. The
respective state evolutions give us the RS and 1RSB phase diagrams. The paper
will maintain a schizophrenic attitude, aiming to connect with both the
physical XY Hamiltonian and the inference problem that motivated it.

The rest of the article is structured as follows: in section II we introduce
the model from the inference point of view and we show the equivalence between
the inference and the order/disorder formulation. In section III we study the
RS phase diagram, with special attention to the outside of the Nishimori line.
In section IV we perform a 1RSB analysis, using the ASP algorithm and its
state evolution. Section V is dedicated to discussion and conclusions.

##  II Definition of the model

We start by introducing the models in the language of statistical inference,
where it is known as $U(1)$ synchronization, an instance of angular
synchronization [1, 10, 13]. The problem consists of retrieving a complex
signal $x^{\star}=(x^{\star}_{1},x^{\star}_{2},\dots,x^{\star}_{N})$, where
each $x^{\star}_{i}$ is uniformly distributed on the unit circle,
independently of other coordinates. In other words,
$x^{\star}_{i}=e^{i\theta_{i}^{\star}}$ with
$\theta_{i}^{\star}\sim\text{Uniform}([0,2\pi])\;i.i.d.$. We will refer to
$x^{\star}$ as the ground truth or the planted signal. A set of $N^{2}$
complex measurements $\\{Y_{ij}\\}_{1\leq i,j\leq N}$ is produced according to
the rule

$Y_{ij}=\sqrt{\frac{\lambda}{N}}x^{\star}_{i}\overline{x^{\star}_{j}}+W_{ij},$
(1)

where $W$ is a Hermitian matrix (i.e. $W_{ji}=\overline{W}_{ij}$) whose
elements above the diagonal are all independent and distributed as
$W_{ij}\sim\mathcal{N}(0,1/2)+i\mathcal{N}(0,1/2)$. We also set
$Y_{ii}=0,\;\forall i\in[N]$. The parameter $\lambda$ plays the role of the
signal to noise ratio, while the scaling $1/\sqrt{N}$ ensures that the problem
of recovering $x^{\star}$ is neither trivial (very large signal-to-noise
ratios) nor impossible (very small signal-to-noise ratios) [22]. The goal of
the inference problem is to recover $x^{\star}$ from the knowledge of
$Y_{ij}$.

We can write the posterior probability of $x$ given $Y$. In doing so, we
assume that the parameter $\lambda$ is unknown, hence we study the family of
probability measures with varying parameter $\hat{\lambda}$ possibly different
from $\lambda$. We stress that $Y$ is always generated using $\lambda$. When
$\lambda=\hat{\lambda}$ computing the marginals of the posterior leads to
Bayes-optimal inference; in the statistical physics language we say that the
Nishimori condition is met. The consequences of this condition are extensively
studied in [14]. We first write the likelihood

$\displaystyle
P(Y|x)=\prod_{i<j}\frac{1}{\pi}\exp\left[-\left|Y_{ij}-\sqrt{\frac{\hat{\lambda}}{N}}x_{i}\overline{x}_{j}\right|^{2}\right].$
(2)

Then, by applying Bayes theorem, we get the posterior measure

$\displaystyle P(x|Y)$ $\displaystyle=$
$\displaystyle\frac{P(Y|x)P(x)}{P(Y)}=\frac{1}{\pi^{\binom{N}{2}}(2\pi)^{N}P(Y)}\exp\left[-\sum_{i<j}\left|Y_{ij}-\sqrt{\frac{\hat{\lambda}}{N}}x_{i}\overline{x}_{j}\right|^{2}\right]$
(3) $\displaystyle=$
$\displaystyle\frac{1}{Z(Y,\hat{\lambda})}\exp\left[2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{i<j}\operatorname{Re}(Y_{ij}\overline{x_{i}}x_{j})\right].$

In order to obtain the final expression, the prior as well as the terms in the
expression of the likelihood that do not depend on $x$, were absorbed into the
normalization. To reconnect with the statistical physics setting, we write the
posterior in terms of a Hamiltonian

$P(x|Y)=\frac{1}{Z(Y,\hat{\lambda})}e^{H(x,Y,\hat{\lambda})},\qquad
H(x,Y,\hat{\lambda})=2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{i<j}\operatorname{Re}(Y_{ij}\overline{x_{i}}x_{j}).$
(4)

To highlight the connection with the gauge glass model $H_{GG}=J\sum_{i,j\in
E}\cos(\theta_{i}-\theta_{j}-\phi_{ij})$, the Hamiltonian can also be written
as
$H=2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{i<j}|Y_{ij}|\cos\left(\theta_{i}-\theta_{j}-\phi_{ij}\right)$,
where $Y_{ij}=|Y_{ij}|e^{i\phi_{ij}}$ and $x_{i}=e^{i\theta_{i}}$. The model
defined above is characterized by two main symmetries. The first being the
$U(1)$ symmetry, from which the model takes its name. It consists of the
invariance of the Hamiltonian under a global phase shift, that is
$H(x_{1},\dots,x_{N})=H(e^{i\phi}x_{1},\dots,e^{i\phi}x_{N})$ for any
$\phi\in[0,2\pi]$. As a consequence, we are able to recover the planted signal
only up to a global phase. The second symmetry is more subtle but allows us to
transform the planted model into an ordered one. This feature is not unique to
$U(1)$ synchronization: for example, the planted SK model enjoys the same
property, and can be transformed into a ferromagnetic model where the
ferromagnetic bias is proportional to the signal to noise ratio in the
original problem [14]. The $U(1)$ synchronization Hamiltonian is invariant
following simultaneous transformation of $x$ and $Y$. Given an arbitrary
vector $z=(z_{1},\dots,z_{N})=(e^{i\phi_{i}},\dots,e^{i\phi_{N}})$, we
transform

$\displaystyle
x_{i}^{\prime}=x_{i}\overline{z_{i}}=e^{i(\theta_{i}-\phi_{i})},$ (5)
$\displaystyle Y^{\prime}_{ij}=Y_{ij}\overline{z_{i}}z_{j}.$ (6)

To obtain a ferromagnetic model in the variables $x^{\prime}$, we pick
$z_{i}=x^{\star}_{i}$. The planted configuration is then transformed into an
ordered one $x^{\star\prime}_{i}=x^{\star}_{i}\overline{x}^{\star}_{i}=1$. For
large $\lambda$, configurations $x$ sampled from the measure will align with
$x^{\star}$. Thus, $x^{\prime}$ will align with $(1,\dots,1)$ (always up to a
global phase shift).

Thanks to this symmetry we can study without loss of generality the problem
with $x^{\star}=(1,1,\dots,1)$ and our results will extend to the case where
$x^{\star}$ is sampled uniformly over the unit circle. Therefore, without loss
of generality, we can restrict our analysis to the problem with measure (4)
and random couplings

$Y_{ij}=\sqrt{\frac{\lambda}{N}}+W_{ij},$ (7)

where $W$ has the same distribution as in (1). We name this particular
instance the planted XY model. The parameter $\lambda$ plays the role of a
ferromagnetic bias, while the parameter $\hat{\lambda}$ is an inverse
temperature.

##  III The RS phase diagram

In this section, we derive the RS phase diagram of the planted XY model. This
is the phase diagram under the approximation that the Gibbs measure can be
represented as a Bethe measure. To obtain the phase diagram we use AMP, a
generalization of the Thouless-Anderson-Palmer equations [23, 12], and its
state evolution, which is equivalent to the replica symmetric solution. Our
analysis also provides a case to characterize the strengths and limitations of
AMP. AMP is a general purpose iterative message passing algorithm which
outputs an approximation to the marginals of the desired probability
distribution: in our case the posterior $P(x|Y)$. While we expect AMP to give
exact results on the Nishimori line ($\lambda=\hat{\lambda}$), in the case of
mismatched parameters ($\lambda\neq\hat{\lambda}$), it can be inaccurate and
fail to converge due to the emergence of replica symmetry breaking (RSB). The
derivation of AMP and its state evolution from Belief propagation are
presented in appendix B. Here, we present the final AMP equations

$\displaystyle h_{i}^{(t)}=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\hat{x}_{k}^{(t)}-\frac{\hat{\lambda}}{N}\hat{x}_{i}^{(t-1)}\sum_{k\neq
i}|Y_{ik}|^{2}\frac{\partial\eta}{\partial h}\left(h_{k}^{(t-1)}\right)$ (8)
$\displaystyle\hat{x}_{i}^{(t+1)}=\eta\left(h_{i}^{(t)}\right),\qquad\eta(h)=\frac{h}{|h|}\frac{I_{1}(|2h|)}{I_{0}(|2h|)},$
(9)

where $I_{k}$ is the modified Bessel function of the first kind of order $k$
and $\frac{\partial\eta}{\partial h}(h)=1-|\eta(h)|^{2}$. $\hat{x}_{i}$ is the
estimator of the mean of the marginal; that is, $\hat{x}_{i}$ estimates
$\operatorname{\mathbb{E}}_{x|Y}[x_{i}]$. One of the elements which
distinguishes AMP from other iterative algorithms is the ability (in the
$N\to\infty$ limit) to track its dynamics through the state evolution
equations. In particular, we have closed evolution equations for the two
observables

$\displaystyle m=\frac{1}{N}\sum_{i=1}^{N}\hat{x}_{i}\qquad
q=\frac{1}{N}\sum_{i=1}^{N}|\hat{x}_{i}|^{2},$ (10)

representing respectively the alignment of the marginals with the planted
configuration and how concentrated each marginal is. When doing inference, one
is interested in the mean square error (MSE) of the estimator. The MSE can be
expressed as $\text{MSE}=1+q-2m$. In the language of statistical physics, $m$
is the order parameter that represents how biased the system is towards the
ordered state. The SE equations read:

$\displaystyle
m^{(t+1)}=\operatorname{\mathbb{E}}_{z}\left[\eta\left(h^{(t)}\right)\right]$
(11) $\displaystyle
q^{(t+1)}=\operatorname{\mathbb{E}}_{z}\left[\left|\eta\left(h^{(t)}\right)\right|^{2}\right],\quad\text{with}$
(12) $\displaystyle
h^{(t)}=\sqrt{\lambda\hat{\lambda}}m^{(t)}+\sqrt{\hat{\lambda}q^{(t)}/2}z,$
(13)

where $z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$ is a complex normal random
variable. Last, we’re also able to compute the Bethe approximation to the free
entropy $f=\frac{1}{N}\log Z$. The Bethe free entropy is derived in Appendix
G, its expression is

$\Phi_{\text{RS}}(\hat{\lambda},\lambda)=-\sqrt{\lambda\hat{\lambda}}m^{2}+\frac{\hat{\lambda}}{2}q^{2}-\hat{\lambda}q+\mathbb{E}_{z}\log[I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\hat{\lambda}q/2}z|)],$
(14)

where $z$ is defined as above and $m,q$ are determined by iterating SE until
convergence. In appendix G, we also show that the stationary points of
$\Phi_{\text{RS}}$ with respect to $\lambda,\hat{\lambda}$ are the fixed
points of the SE equations. This confers another interpretation of the SE
equations, the one of an iterative method to find the stationary points of the
free entropy. We perform extensive numerical experiments with the goal of
studying our model through the lenses of SE and AMP.

###  III.1 On the Nishimori line

We start by restricting ourselves to the Nishimori line
$\lambda=\hat{\lambda}$. For a large class of models, including ours, it was
proven that the RS solution is exact in the large size limit [24]. From the
inference point of view, this corresponds to the case where we know how the
data $Y$ is generated, and we can perform Bayes Optimal inference, in the
sense that $\hat{\lambda}$ matches $\lambda$. As a consequence of this fact,
one can establish the relation $m=q$ [18] [14]. AMP’s analysis on the
Nishimori line has been partially conducted in [13] for a class of models that
includes ours. Moreover, a free entropy equivalent to ours has been obtained
in [13] and in [11] via the replica method and proven to be correct in a more
general setting in [24].

Figure 1 illustrates the behavior of SE and AMP on the Nishimori line. On the
left panel, the converged values of $m$ and $q$, from both SE and AMP are
plotted as a function of $\lambda$. First we observe the agreement between AMP
and SE, i.e. SE’s fixed points have the same $m,q$ as the points to which AMP
converges. Next we see that, as one would expect, $m=q$ at convergence. In the
center panel, SE is run from initial conditions $m^{t=0}=0,q^{t=0}=10^{-2}$,
corresponding AMP initialized randomly near $0$ (in principle one would like
to use $m^{t=0}=0,q^{t=0}=0$ but numerical errors arise when initializing with
too small $q$), and $m^{t=0}=1,q^{t=0}=1$ (i.e. initializing
$\hat{x}^{t=0}=(1,1,\dots,1)$ in AMP), known as informed initialization. The
two iterations converge to the same value of $m$, showing that the SE fixed
point is unique in $m$. The right panel provides a free entropy interpretation
of this phenomenon. The $m$ dependent free entropy clearly has only one
maximum, hence SE inevitably lands on it. In other models with multimodal free
energies [22], the local maximum to which SE (and hence AMP) converges might
not be the absolute one. Being the RS free entropy $\Phi_{\text{RS}}$ exact on
the Nishimori line, the uniqueness allows us to conclude that AMP computes the
true marginals when $N\to\infty$. Beside AMP’s behavior, Figure 1 shows that
$m$ undergoes a second order phase transition at $\lambda=1$: for $\lambda<1$
we say the model is in the paramagnetic phase, because $m=0$ and there is no
correlation with the planted signal. Instead for $\lambda>1$ the system
develops a ferromagnetic order.

Figure 1: Numerical results for AMP and SE on the Nishimori line. All plots
show the values of $m,q$ at convergence. Left panel: comparing $m$ and $q$
from AMP to verify that $m=q$; also verifying that they both agree with SE
(continuous line). The initial condition for AMP is uninformative
($m^{t=0}=0$). Center panel: the two curves represent the value of $m$ to
which SE converges with respective initial condition $m^{t=0}=0,q^{t=0}=1$ and
$m^{t=0}=1,q^{t=0}=1$; the fact that they’re equal shows that there is a
unique fixed point of SE. Right Panel: Bethe free entropy as a function of
$m$; the red cross on each curve marks the unique stationary point to which SE
converges; the uniqueness of the SE fixed point is a direct consequence of the
free entropy having exactly one stationary point.

###  III.2 Contiguity to the random XY model

One important consequence of the previous analysis is the existence of a phase
where the effect of the planted signal disappears. When $\lambda<1$, the
signal is completely washed out by the noise, and the data $Y$ is
indistinguishable from random noise; in other words, it is as if $Y=W$. This
in turn implies that for $\lambda<1$, all high-probability quantities are
independent of $\lambda$ (this is referred to as contiguity of probability
distributions in the statistics literature [25]). To put it differently,
whenever $\lambda<1$, the planted nature of our model is lost and we look at a
spin glass with purely random couplings. In the rest of the article, we will
refer to the $\lambda<1$ case as the fully random phase.

###  III.3 Replica symmetric phase diagram

Moving to the general case where $\hat{\lambda}$ is possibly different from
$\lambda$, we aim to explore the full phase diagram painted by SE. We begin by
verifying again the agreement between the fixed points of AMP and SE, also
during the dynamics. In Figure 2, AMP is initialized from a random
configuration $x^{t=0}$ with entries of unit norm (thus $q^{t=0}=1$ and
$m^{t=0}=0$). It’s evident that SE accurately tracks AMP, apart for some
finite size effects.

Figure 2: Comparison between AMP ($N=1000$) and SE. The red line represents
the behavior of SE while each of the 40 gray lines is an independent AMP run.
State evolution tracks accurately both $m$ and $q$. All the runs are
initialized with a configuration where each spin’s value is picked uniform
over the unit circle.

In Figure 3, we plot the full phase diagram in the $\hat{\lambda},\lambda$
plane. We can identify several phases:

  * •

RS unstable phase The yellow area in the right panel is the region where AMP
does not converge. The non-convergence of AMP is synonym of the replica
symmetric solution being unstable and the RSB being required to correctly
model the measure. The equivalence between the convergence of AMP and the
stability of the RS phase is further discussed in section IV.2.

  * •

Paramagnetic phase corresponding to $m=0,\;q=0$. The boundary of this region
can be found analytically (see Appendix C) and is given by the curve
$\lambda=\min(1,\hat{\lambda}^{-1})$. Throughout the paramagnetic phase, AMP
will output the non-informative estimator $\hat{x}=(0,0,\dots,0)$. This
corresponds to the estimated marginals being uniform on the circle.

  * •

Ferromagnetic phase defined the intersection of the region where $m>0,\;q>0$
and the RS stable phase. The marginals produced by AMP are partially aligned
with the planted state (or partially ordered). Moreover the RS solution
correctly describes the structure of the Gibbs measure.

  * •

Mixed phase where $m>0,\;q>0$ but the RS solution is unstable. This region
corresponds to the slice between the ferromagnetic and spin glass phase.

  * •

Spin glass phase where $m=0,\;q>0$. In the spin glass phase, AMP’s marginals
are polarized towards a random value which is uncorrelated with the planted
configuration. In this phase, AMP also encounters convergence problems, and
the RS solution is unstable. In Appendix C we show that the upper boundary of
the spin glass phase, dividing the $m=0$ from the $m>0$ region, can also be
expressed analytically in an implicit form. Notice that since we are in the RS
unstable phase, the RS prediction for this boundary is not reliable, thus RSB
is needed to evaluate the correct boundary between the spin glass and the
mixed phase.

More quantitative information about the values of $q$ and $m$ is found in
Figure 4. In this figure, the top row represents the phase diagram obtained
via AMP, while the bottom row was obtained from SE. First, we notice that
across all transition lines, both $m$ and $q$ are continuous. Looking at the
panels showing $q$, we see that $q$ is always increasing in $\hat{\lambda}$.
This can be explained by interpreting $\hat{\lambda}^{1/2}$ as an inverse
temperature, then it’s clear that spins should be more and more polarized with
decreasing temperature. From a mathematical perspective according to (8),
$\hat{\lambda}^{1/2}$ controls the norm of $h^{t}$ and hence that of
$\hat{x}$.

Figure 3: Phase diagram obtained from SE. The dashed line in both panels is
the Nishimori line. Left panel: the light blue region is paramagnetic phase,
where $m=q=0$. It is delimited by the curve
$\lambda=\min(1,\hat{\lambda}^{-1})$. The pink bottom right region is the spin
glass phase (evaluated by the RS solution): here $m=0$ while $q>0$. The upper
curve delimiting the spin glass phase (evaluated on the RS level) converges to
$\lambda=4/\pi$ for $\hat{\lambda}\to\infty$. In the top right we find the
ferromagnetic region where both $m$, $q$ are strictly positive and the RS
solution is stable. The remaining slice between the ferromagnetic and spin
glass phase is called the mixed phase. Right panel: The yellow area, which
encompasses the spin glass phase and the mixed phase, is the region where the
RS stability parameter $c$ is negative and thus AMP does not converge, or
equivalently the RS solution is unstable.

Figure 4: Heatmap representing $m$ and $q$ for AMP (upper panels, $N=500$) and
SE (lower panels) as a function of $\lambda$, $\hat{\lambda}$.

###  III.4 Convergence of AMP

It is known that on the Nishimori line the RS ansatz is exact [14, 24], hence
AMP estimates the marginals of $P(x|Y)$ exactly in the large size limit. The
same cannot be said about the rest of the phase diagram. Therefore we need to
distinguish between the true behavior of the model and that of the algorithm.
For example, AMP’s shortcomings are evident when the iterations fail to
converge. In the left panel of Figure 5, we plot the convergence time of the
algorithm across the phase diagram. For small enough $\lambda$, when
$\hat{\lambda}$ is increased, we always encounter a phase in which AMP does
not converge. AMP’s convergence has important links to RS stability. In the
Sherrington Kirkpatrick (SK) model, it was proved in [26] that the RS
stability line delimits the region where AMP converges. This property is
general and also in our case, the 1RSB analysis will confirm that convergence
of AMP and stability of RS coincide. We can thus state that AMP converges iff
the RS solution is stable. Analytically, we study AMP’s convergence by looking
at stability under a random perturbation $h^{(t)}\mapsto h^{(t)}+\delta
h^{(t)}$, with $\delta h^{(t)}_{k}=\epsilon e^{i\phi_{k}^{(t)}}$, with
$\phi_{k}^{(t)}\sim\text{Unif}([0,2\pi])$. By propagating the perturbation
through the AMP equations (more details are provided in Appendix D), we obtain
that the perturbation norm grows according to the law
$\operatorname{\mathbb{E}}_{\delta h}\left[\left\lVert\delta
h^{(t+1)}\right\rVert^{2}\right]=(1-c_{\text{AMP}}^{(t)})\operatorname{\mathbb{E}}_{\delta
h}\left[\left\lVert\delta h^{(t)}\right\rVert^{2}\right]$, with

$c_{\text{AMP}}^{(t)}=1+\frac{\hat{\lambda}}{2}\sum_{i=1}^{N}\left[\frac{1}{N}\sum_{k}|Y_{ik}|^{2}\right]\left[\eta_{r}(|h_{i}^{(t)}|)^{2}+\left(\eta_{r}(|h_{i}^{(t)}|)+|h_{i}^{(t)}|\eta_{r}^{\prime}(|h_{i}^{(t)}|)\right)^{2}\right]$
(15)

where the expectation is with respect to the randomness in the perturbation
and $\eta_{r}(r)=\frac{1}{r}\frac{I_{1}(2r)}{I_{0}(2r)}$. AMP will converge if
the norm of the perturbation decreases in time ($c_{\text{AMP}}>0$), and will
oscillate otherwise. This quantity can also be tracked using state evolution

$\displaystyle
c^{(t)}_{\text{SE}}=1+\frac{\hat{\lambda}}{2}\operatorname{\mathbb{E}}_{z}\left[\eta_{r}(|h|)^{2}+\left(\eta_{r}(|h|)+|h|\eta_{r}^{\prime}(|h|)\right)^{2}\right]$
(16) $\displaystyle
h=\sqrt{\lambda\hat{\lambda}m^{(t)}}+\sqrt{\frac{\hat{\lambda}q^{(t)}}{2}}z,\quad
z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$ (17)

where $m^{(t)},q^{(t)}$ are obtained by running SE. Since $c_{\text{AMP}}$
converges to $c_{\text{SE}}$ in probability in the $N\to\infty$ limit, we will
refer to both quantities as $c$. The right panel in Figure 3 depicts in yellow
the region where AMP is not convergent. The area where $c<0$ coincides with
the union of the mixed and the spin glass phases. We conducted further
experiments about AMP’s convergence properties. In the left panel in Figure 5
we plot the number of iterations (capped at 300) after which AMP converges. We
verify that the region of non convergence coincides with the one predicted
from $c_{\text{AMP}}$ and $c_{\text{SE}}$, displayed in the center plots.
Finally, the right panel shows the quantity
$\Delta\hat{x}^{t}=\frac{1}{N}\sqrt{\sum_{i}\big{|}\hat{x}^{(t)}_{i}-\hat{x}_{i}^{(t-1)}\big{|}^{2}}$,
representing the rate of change of the estimator. We see that in the $c<0$
region, $\Delta\hat{x}^{t}$ does not decay to zero, because the dynamics keeps
oscillating. Contrary to AMP, SE is not affected by convergence problems and
correctly tracks the observables $m$ and $q$, even when $\hat{x}^{(t)}$ does
not converge.

Figure 5: Convergence properties of AMP. To measure the rate of change of the
estimator, we introduce the quantity
$\Delta\hat{x}^{t}=\frac{1}{N}\sqrt{\sum_{i}\big{|}\hat{x}^{(t)}_{i}-\hat{x}_{i}^{(t-1)}\big{|}^{2}}$.
We consider AMP’s iterations to have converged when
$\Delta\hat{x}^{t}<10^{-5}$. Left panel: convergence time of AMP over the
$\lambda,\hat{\lambda}$ plane. The convergence time is capped at $300$, each
pixel represents the average convergence time over 25 independent runs. Center
panel: (top) heatmap of $c_{\text{AMP}}$ computed at the final time, each
pixel represents the average convergence time over 25 independent runs,
(bottom) heatmap of $c_{\text{SE}}$ computed when SE has converged. Right
panel: we pick two points in the $\lambda,\hat{\lambda}$ plane (marked by
white crosses in the center panel) and plot $\Delta\hat{x}^{t}$ for 40
independent AMP runs. When $\Delta\hat{x}^{t}$ does not decay to zero, AMP
does not converge.

In conclusion, AMP and SE correspond to a replica symmetric approach in
solving the planted XY model. Under this assumption the Gibbs distribution is
well described by a single Bethe measure (or state). In the RS unstable phase,
this approximation may not hold and an analysis based on replica symmetry
breaking is needed.

##  IV 1RSB analysis

To overcome the shortcomings of the RS approach and AMP, we must carry out a
more refined analysis which takes into account replica symmetry breaking.
Methods such as BP or AMP basically ‘fit’ the Gibbs measure onto a Bethe
measure. When this ansatz turns out to be correct, we say the model is in a
replica symmetric phase and AMP converges, giving an accurate estimate of the
marginals. Otherwise, the Gibbs measure $\mu$ can break into a multitude of
Bethe measures which we index by $\alpha$, i.e.
$\mu=\frac{1}{Z}\sum_{\alpha}Z_{\alpha}\mu_{\alpha}$. The total partition
function is then $Z=\sum_{\alpha}Z_{\alpha}$, where $Z_{\alpha}$ is the
partition function of a single Bethe state (computable by exponentiating the
$f_{\text{Bethe}}$ of the single state). Message passing algorithms will
exhibit a multitude of fixed points, each corresponding to one of the states.
Replica symmetry breaking allows us to account for this structure of the Gibbs
measure.

In the first step of the construction, called 1-step RSB (1RSB), we postulate
the existence of a function $\Sigma(f)$, called the complexity function, with
the property that the number of states with a free entropy
$f_{\alpha}=\frac{1}{N}\log Z_{\alpha}$ close to a value of $f$ is, at the
leading order, $e^{N\Sigma(f)}$. The best approximation to the free entropy of
the system, at the 1RSB level then reads

$\tilde{\Phi}_{\text{1RSB}}=\frac{1}{N}\log\sum_{\alpha}e^{Nf_{\alpha}}=\frac{1}{N}\log\int
df\,e^{N(\Sigma(f)+f)}=\sup_{f:\Sigma(f)\geq
0}\left[\Sigma(f)+f\right]=\Sigma(\tilde{f}^{*})+\tilde{f}^{*}$ (18)

Here, $\tilde{f}^{*}$ is the free entropy of each of the equilibrium states
and it’s determined by the condition

$\tilde{f}^{*}=\operatorname*{\arg\\!\max}_{f:\Sigma(f)>0}\left(\Sigma(f)+f\right)=\begin{cases}f:\Sigma^{\prime}(f)=-1&\text{
if $\Sigma(f)>0$}\\\ f_{0}&\text{ otherwise }\end{cases}$ (19)

with $f_{0}$ the largest root of the equation $\Sigma(f)=0$. We will call
$\tilde{\Phi}_{\text{1RSB}}$ the _true_ 1RSB free entropy. By introducing the
positivity constraint on $\Sigma$ in (19), we are discarding the unphysical
solutions with $\Sigma(f)<0$. The negative complexity would in fact mean that
there is an exponentially (in $N$) small probability of finding a state with
the given free entropy.

From ASP we will obtain the related quantity, which we call _replicated_ free
entropy

$\Phi_{\text{1RSB}}(s)=\sup_{f}\left[\Sigma(f)+sf\right]=\Sigma(f^{*}(s))+sf^{*}(s)$
(20)

where $f^{*}$ satisfies $\Sigma^{\prime}(f^{*})=-s$. Notice that from (20) we
have the characterization $f^{*}(s)=\frac{\partial\Phi_{\text{1RSB}}}{\partial
s}$. We also remark that we can access $\Sigma(f)$ by computing it
parametrically in $s$.

$\Sigma(f^{*}(s))=\Phi_{\text{1RSB}}(s)-sf^{*}(s)$ (21)

The next goal is to find $\tilde{\Phi}_{\text{1RSB}}$, starting from the newly
found $\Sigma$ and $f^{*}$. One difference between $\Phi_{\text{1RSB}}$ and
$\tilde{\Phi}_{\text{1RSB}}$ is the relaxation of the positivity constraint on
$\Sigma$. This implies that naively setting $s=1$ might not give
$\tilde{\Phi}_{\text{1RSB}}$. Instead we obtain
$\tilde{\Phi}_{\text{1RSB}}=\Sigma(f^{*}(s_{\star}))+f^{*}(s_{\star})$ for a
well chosen $s_{\star}$

$s_{\star}=\begin{cases}1&\text{ if $\Sigma(f^{*}(1))\geq 0$}\\\
-\Sigma^{\prime}(f_{0})&\text{ if $\Sigma(f^{*}(1))<0$ }\end{cases}$ (22)

and justify this choice in Appendix E. Basically, $s_{\star}$ is the value of
$s$ for which the replicated free entropy best approximates the Gibbs measure.

###  IV.1 Approximate survey propagation

In this section, we derive the ASP algorithm. Survey propagation (SP) is a
message passing algorithm developed originally in the context of random
constraint satisfaction problems [27]. The approach has then been extended to
several other inference and optimization problems [28][29]. ASP, through its
state evolution, also allows us to compute the 1RSB replicated free entropy
exactly in the $N\to\infty$ limit. Appendix F provides the details of the
derivation, here we only go over the key steps.

In this work we follow the derivation of ASP introduced in [2]. We introduce a
replicated system composed of $s$ independent replicas. We indicate with
$x_{i}=(x_{i}^{(1)},\dots,x_{i}^{(s)})$, the replicated variables. First we
write BP for the replicated system $\\{x^{(a)}\\}_{a=1}^{s}$

$\displaystyle{m}_{{i}\rightarrow{}{ij}}(x_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{i\to ij}}\prod_{k\neq
i,j}{m}_{{ki}\rightarrow{}{i}}(x_{i})$ (23)
$\displaystyle{m}_{{ij}\rightarrow{}{i}}(x_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{ij\to i}}\int
dx_{j}{m}_{{j}\rightarrow{}{ij}}(x_{j})\exp\left[2\sqrt{\frac{\lambda}{N}}\sum_{a=1}^{s}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right],$
(24)

Then we relax BP by parametrizing the messages with their means and
covariances

$\displaystyle\langle x_{i}^{(a)}\rangle=\hat{x}_{i\rightarrow ij}$ (25)
$\displaystyle\langle
x_{i}^{(a)}\overline{x}_{i}^{(b)}\rangle=\begin{cases*}\left|\hat{x}_{i\rightarrow
ij}\right|^{2}+\Delta_{i\rightarrow ij}&if $a\neq b$\\\ 1&if
$a=b$.\end{cases*}$ (26)

with $\langle\cdot\rangle$ being the average with respect to $m_{i\rightarrow
ij}$. $\Delta$ is a measure of how coupled the replicas are: when $\Delta=0$
the replicas are independent and we recover the BP equations for $s$
independent replicas. The crucial assumption of the 1RSB approximation is that
all pairs of replicas $(a,b):\;a\neq b$ have the same $\Delta$. $\hat{x}$
instead plays the same role as in AMP. Finally, we remove the dependence of
the messages on the target node and correct for it by introducing the Onsager
term. Once again, this is only possible thanks to the fully connected nature
of the model. After accomplishing these steps, we arrive to the ASP equations

$\displaystyle{T}_{i}^{t}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k=1}^{N}Y_{ik}{\hat{x}}_{k}^{t}-\frac{\hat{\lambda}}{N}{\hat{x}}_{i}^{t-1}\sum_{k}|Y_{ik}|^{2}\frac{d{\hat{x}}_{k}}{dT_{k}}$
(27) $\displaystyle{V}_{i}^{t}$ $\displaystyle=$
$\displaystyle\frac{\hat{\lambda}}{N}\sum_{k=1}^{N}|Y_{ik}|^{2}{\Delta}_{k}^{t}$
(28) $\displaystyle{\hat{x}}_{i}^{t+1}$ $\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s-1}\left[\frac{T_{i}+V_{i}h}{|T_{i}+V_{i}h|}I_{1}(2|T_{i}+V_{i}h|)\right]}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s}}$
(29) $\displaystyle{\Delta}_{i}^{t+1}$ $\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s-2}\left[I_{1}(2|T_{i}+V_{i}h|)\right]^{2}}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s}}-|{\hat{x}}_{i}^{t+1}|^{2}.$
(30)

with $\frac{d{\hat{x}}_{k}}{dT_{k}}$, computed numerically via finite
differences.

It can be shown that by setting $s=1$, one recovers AMP. The ASP equations are
equipped with their SE, which allows us to track the scalar quantities $m$,
$q$ and $\Delta\coloneqq\frac{1}{N}\sum_{i}\Delta_{i}^{t}$. SE reads

$\displaystyle T^{t}$ $\displaystyle=$
$\displaystyle\sqrt{\hat{\lambda}\lambda}m^{t}+\sqrt{\hat{\lambda}q^{t}/2}z$
(31) $\displaystyle V^{t}$ $\displaystyle=$
$\displaystyle\hat{\lambda}\Delta^{t}$ (32) $\displaystyle m^{t+1}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[\hat{x}(T^{t},V^{t})\right]$ (33)
$\displaystyle q^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[|\hat{x}(T^{t},V^{t})|^{2}\right]$ (34)
$\displaystyle\Delta^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[\Delta(T^{t},V^{t},q^{t})\right].$ (35)

The functions $\hat{x}(T,V)$ and $\Delta(T,V,q)$ are the same as in (27), but
without indices and with $|\hat{x}_{i}|^{2}$ replaced by $q$. Finally, we
compute the 1RSB replicated free entropy for ASP (see Appendix G for the
derivation)

$\displaystyle\Phi_{\mathrm{1RSB}}(\lambda,\hat{\lambda},s)$ $\displaystyle=$
$\displaystyle-s\sqrt{\hat{\lambda}\lambda}m^{2}+\frac{s^{2}\hat{\lambda}}{2}q^{2}-s\hat{\lambda}(q+\Delta)-\frac{s(s-1)\hat{\lambda}}{2}(\Delta+q)^{2}$
$\displaystyle+\mathbb{E}_{z}\left[\log\left(\frac{\hat{\lambda}\Delta}{\pi}\int_{\mathbb{C}}dh\,\exp(-\hat{\lambda}\Delta|h|^{2})[I_{0}(2|T+\hat{\lambda}\Delta
h|)]^{s}\right)\right]$

and the corresponding free entropy of the states selected by $s$

$\displaystyle f^{*}(s,\lambda,\hat{\lambda})$ $\displaystyle=$
$\displaystyle\frac{\partial\Phi_{\mathrm{\text{1RSB}}}(s,\lambda,\hat{\lambda})}{\partial
s}=-\sqrt{\lambda\hat{\lambda}}m^{2}+s\lambda
q^{2}-\lambda(q+\Delta)-\frac{(2s-1)\lambda}{2}(\Delta+q)^{2}$
$\displaystyle+$
$\displaystyle\mathbb{E}_{z}\left[\frac{\int_{\mathbb{C}}dh\exp(-\lambda\Delta|h|^{2})\log[I_{0}(2|T+\lambda\Delta
h|)][I_{0}(2|T+\lambda\Delta
h|)]^{s}}{\int_{\mathbb{C}}dh\exp(-\lambda\Delta|h|^{2})[I_{0}(2|T+\lambda\Delta
h|)]^{s}}\right].$

In both the free entropy and SE, $z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$.
Analogously to AMP, the stationary points of the replicated free entropy
$\Phi_{\text{1RSB}}$ with respect to $m,q,\Delta$ are fixed points of the ASP
equations; a derivation of this fact is provided in Appendix G. By
manipulating the equations, it can be shown that there are two ways to recover
the RS solution: either by setting $s=1$ or by having $\Delta=0$.

###  IV.2 Numerical results

Iterating equations (27) and (31) presents some challenges due to the multiple
integrations involved at every time step. Nonetheless we manage to compute, to
satisfactory numerical precision, all the key quantities in the problem: the
complexity $\Sigma$, the 1RSB free entropy $\tilde{\Phi}_{1RSB}$, the
equilibrium free entropy $\tilde{f}^{*}$ and $s_{\star}$. We start by
verifying that ASP and SE behave as expected. Figure 6 illustrates some ASP
numerical experiments conducted at $\lambda=0.5$,
$\hat{\lambda}=2,s=s_{\star}=0.166$, where the model is fully random (no
ferromagnetic bias), and at $\hat{\lambda}=2,\lambda=1.08,s=s_{\star}=0.221$,
located in the mixed phase. Both points are located in the RS unstable region.
First, we remark that SE tracks ASP apart from some finite size effects.
Moreover we observe that $\Delta\hat{x}$ does not always go to zero, thus ASP
doesn’t always converge in the RS unstable phase. Figure 6 also confirms the
existence of the RSB phase, characterized by $\Delta>0$, contrary to the RS
phase where $\Delta=0$. For the rest of the analysis, we will present results
obtained exclusively from SE.

Figure 6: ASP and its state evolution at two different points in the
$\lambda,\hat{\lambda}$ plane. For each run we plot $m,\,q\,,\Delta$ and
$\Delta\hat{x}^{t}$ as a function of iteration time. The red curves indicate
the evolution predicted by SE. Each of the 50 transparent lines is an
independent ASP run with $N=2000$. Left panel: we run ASP deep in the spin
glass phase $\hat{\lambda}=2$, $\lambda=0.5,s=s_{\star}=0.166$; recall that
here our model is equivalent to a fully random one with $\lambda=0$. Right
panel: ASP is run at $\hat{\lambda}=2,\lambda=1.08,s=s_{\star}=0.221$. This
point is located in the mixed phase. In both cases SE tracks ASP accurately
enough, the finite discrepancy in the case of $\lambda=1.05$ is caused by the
finite size effect near the RS instability transition located at
$\lambda=1.105$.

To capture the behavior of SE, we study the algorithm along two trajectories
in the $\hat{\lambda},\lambda$ plane. Figure 7 fixes $\hat{\lambda}=2$ and
plots several quantities as a function of $\lambda$. First we notice that for
$\lambda<1$, $s_{\star},\Delta,m$ and also $q$ (not shown) are constant. In
fact, in this region, the model is equivalent to a random one.

Recall that the AMP convergence threshold is at
$\lambda_{\text{conv}}(\hat{\lambda}=2)=1.105$, indicated by the upper limit
of the grey band. We indeed observe that $\Delta$ undergoes a second order
transition (upper center plot) exactly at $\lambda_{\text{conv}}$; moreover
for $\lambda>\lambda_{\text{conv}}$ the 1RSB free entropy becomes equal to the
RS free entropy (bottom left plot). These results confirm that AMP convergence
and RS stability are equivalent.

We then see that in the whole range of $\lambda$, $s_{\star}<1$. This
indicates that there is no dynamical-1RSB phase, where the measure would be
dominated by an exponential number of states. On the 1RSB level, the measure
is either dominated by a single Bethe state in the RS phase
($\lambda>\lambda_{\text{conv}}$) or, in the 1RSB phase, by a sub exponential
number of Bethe states. This fact implies that
$\tilde{\Phi}_{\text{1RSB}}(\lambda,\hat{\lambda})=f^{*}(s_{\star},\lambda,\hat{\lambda})$.
Put differently, in the 1RSB phase, the free entropy of the system will be
given by the point where the complexity curve $\Sigma(f)$ touches zero.
Nonetheless, below $\lambda_{\text{conv}}$, $\Sigma$ (center panel) is
positive on a interval of values of $s$, attesting the presence of an
exponential number of metastable states. Approaching $\lambda_{\text{conv}}$
from below, $\Sigma(f)$ shrinks until it becomes a point at the RS stability
transition, correspondingly $\Sigma(s)$ approaches a zero function. The
complexity curves also have an unphysical branch, only shown for
$\lambda=0.9$. The unphysical branch begins when $f^{*}$ becomes decreasing in
$s$, and continues down to $s=0$.

The behavior of $m$ (top right) is also interesting: the lower margin of the
gray stripe corresponds to the value of $\lambda$ at which AMP starts
correlating with the ordered state, i.e. when $m_{\text{AMP}}>0$. We see that
ASP achieves positive $m$ even when AMP is not able to, almost achieving the
theoretical threshold at $\lambda=1$. From an inference point of view we can
say that ASP recovers the signal at smaller SNR.

Let us shift our attention to Figure 8. We fix $\lambda=0.5$ (fully random
phase) and vary $\hat{\lambda}$. Recall that for $\hat{\lambda}\leq 1$ the
measure is RS. For increasing $\hat{\lambda}$, $\Delta$ increases, signaling
that the states’ width decreases. As in the previous figure, when approaching
the RS region the complexity curves $\Sigma(s)$ become flatter, becoming the
identically zero function at $\hat{\lambda}=1$, and $\Sigma(f)$ becomes a
point at $f^{*}$. When this happens, the free entropy becomes independent of
$s$, giving back the RS solution, which corresponds to $s=1$. From the point
of view of $\Sigma(f)$, smaller values of $\hat{\lambda}$ translate into fewer
metastable states in the measure, approaching eventually the RS picture with
only one state, the paramagnetic one.

Finally we discuss the stability of ASP in both Figures 7 and 8. In a way
analogous to AMP we study the convergence of ASP by analyzing the stability of
its iterations under a random perturbation. Given the complexity of the update
functions, we perform this analysis numerically by perturbing the vector
$T^{t}$ by a small quantity. We introduce $c_{\text{ASP}}$ such that if
$c_{\text{ASP}}>0$ the perturbation norm shrinks in time. If instead
$c_{\text{ASP}}<0$, the perturbation grows in time and ASP is unstable. The
bottom right panels in 7 and 8 show where $c_{\mathrm{ASP}}$ is positive.
Notice for example that in figure 7 ASP converges in the mixed phase, where
AMP failed to converge.

Figure 7: $\hat{\lambda}=2$ and variable $\lambda$. The grey band represents
the mixed phase from the RS phase diagram (i.e. goes from the $m=0$
instability up to the AMP convergence threshold). Top panel: from left to
right we plot, as a function of $\lambda$, the value of $s$ that best
approximates the true free entropy ($s_{\star}$), the overlap between two
replicas in the same state ($\Delta$) and the overlap with the ordered
configuration ($m$). The point at which $\Delta(s_{\star})$ becomes nonzero
marks the RS instability transition, this is seen to coincide with the AMP
convergence threshold. As shown in the top right plot, ASP attains $m>0$ for
smaller values of $\lambda$ compared to AMP. Center Plots: Plotting the
complexity curves $\Sigma(f^{*})$ and $\Sigma(s)$. The curves are
characterized by a physical and an unphysical branch (shown in transparency
only for $\lambda=0.9$) which extends up to $s=0$. At the RS instability
threshold $\Sigma(f)$ collapses to a point and $\Sigma(s)$ becomes constant.
Bottom Plots: On the left we plot the difference between the RS and 1RSB free
energies; the two become equal at the RS instability threshold. On the right
we plot the stability coefficient $c_{\text{ASP}}$ of ASP. Figure 8:
$\lambda=0.5$ and variable $\hat{\lambda}$. For $\hat{\lambda}\leq 1$ the
system is in the RS phase. Top panel: from left to right we plot, for various
values of $\hat{\lambda}$, the value of $s$ that best approximates the free
entropy ($s_{\star}$) and the overlap between two replicas in the same state
($\Delta$). The error bars reflect the error in the computation of
$s_{\star}$. Notice how $\Delta\to 0$ when approaching the RS phase. Bottom
Panel: (left and center) plotting the complexity curves for different values
of $\lambda$, as a function of both $f^{*}$ and $s$. The unphysical branch is
shown in transparent trait only for $\hat{\lambda}=3$. The inset in the left
plot magnifies the region near $f^{*}=f^{*}(s_{\star})$. Approaching the RS
phase, $\Sigma(f)$ collapses to a point, while $\Sigma(s)$ approaches a
constant. (right): plotting the ASP stability coefficient $c_{\text{ASP}}$.

##  V Discussion

In this work, we studied the planted XY model defined by (4) and (7). Our
model enjoys multiple connections, both with spin glasses and inference
problems. In the inference setting, it corresponds to the angular
synchronization problem, instead from the statistical physics point of view,
it can be seen as an XY, mean field spin glass with a ferromagnetic bias. The
two problems are related by a change of variables. We first obtain the RS
phase diagram where we recognize several regions, a paramagnetic phase where
the spins are each uniformly distributed on the circle, a ferromagnetic phase
where a partial global order arises (spins approximately point in the same
direction) but the replica symmetric solution is stable, a mixed phase where
replica symmetry is not stable and magnetization is positive, and a spin glass
phase in which each spin is partially frozen in a random direction.

To mitigate the instability of the RS solution, we resort to the ASP
algorithm. ASP is the 1RSB version of AMP, allowing to model the existence of
multiple states in the Gibbs measure, each state corresponding to a fixed
point of AMP. In the 1RSB formalism, we obtain a better estimate of the free
entropy, and we can also count the number of states with each free entropy
through the complexity function $\Sigma(f)$. All the estimates are obtained
through the state evolution of ASP.

One question remains open: Is the 1RSB approach exact or are further levels of
RSB required? There are two failure modes of 1RSB: either several 1RSB states
form a cluster together, or each 1RSB state breaks into a multitude of smaller
states [30][29][31]. Studying ASP’s convergence allows us to detect the second
kind of instability. We can then conclude that the 1RSB ansatz is incorrect in
at least large part of the phase diagram (i.e. where ASP does not converge or
equivalently where $c_{\text{ASP}}<0$) where likely the Full-RSB ansatz would
be needed to provide the exact solution. In the region where ASP converges and
$c_{\text{ASP}}>0$ one would need to evaluate the 2RSB solution to
conclusively decide about the exactness of the 1RSB, this is left for future
work.

## Acknowledgments

This work started as a part of the doctoral course Statistical Physics For
Optimization and Learning taught at EPFL in spring 2021\. The work of Siyu
Chen was supported by the Swiss National Science Foundation under contract
200021-178999. We acknowledge funding from the ERC under the European Union’s
Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe.

## Data and materials availability

The code and data used to produce the plots within this work will be released
on the repository Zenodo upon publication of this preprint.

## References

  * Singer [2011] A. Singer, Angular synchronization by eigenvectors and semidefinite programming, Applied and computational harmonic analysis 30, 20 (2011).
  * Antenucci  _et al._ [2019] F. Antenucci, F. Krzakala, P. Urbani, and L. Zdeborová, Approximate survey propagation for statistical inference, Journal of Statistical Mechanics: Theory and Experiment 2019, 023401 (2019).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Sherrington [1983] D. Sherrington, The infinite-ranged m-vector spin glass, in  _Heidelberg Colloquium on Spin Glasses_ (Springer, 1983) pp. 125–136.
  * Karp  _et al._ [2003] R. Karp, J. Elson, D. Estrin, and S. Shenker,  _Optimal and Global Time Synchronization in Sensornets_ , Tech. Rep. (2003).
  * Giridhar and Kumar [2006] A. Giridhar and P. R. Kumar, Distributed clock synchronization over wireless networks: Algorithms and analysis, in  _Proceedings of the 45th IEEE Conference on Decision and Control_ (IEEE, 2006) pp. 4915–4920.
  * Bandeira  _et al._ [2014] A. S. Bandeira, M. Charikar, A. Singer, and A. Zhu, Multireference alignment using semidefinite programming, in  _Proceedings of the 5th conference on Innovations in theoretical computer science_ (2014) pp. 459–470.
  * Agrawal  _et al._ [2006] A. Agrawal, R. Raskar, and R. Chellappa, What is the range of surface reconstructions from a gradient field?, in  _European conference on computer vision_ (Springer, 2006) pp. 578–591.
  * Rubinstein and Wolansky [2001] J. Rubinstein and G. Wolansky, Reconstruction of optical surfaces from ray data, Optical review 8, 281 (2001).
  * Bandeira  _et al._ [2017] A. S. Bandeira, N. Boumal, and A. Singer, Tightness of the maximum likelihood semidefinite relaxation for angular synchronization, Mathematical Programming 163, 145 (2017).
  * Javanmard  _et al._ [2016] A. Javanmard, A. Montanari, and F. Ricci-Tersenghi, Phase transitions in semidefinite relaxations, Proceedings of the National Academy of Sciences 113, E2218 (2016).
  * Donoho  _et al._ [2009] D. L. Donoho, A. Maleki, and A. Montanari, Message-passing algorithms for compressed sensing, Proceedings of the National Academy of Sciences 106, 18914 (2009).
  * Perry  _et al._ [2018] A. Perry, A. S. Wein, A. S. Bandeira, and A. Moitra, Message-passing algorithms for synchronization problems over compact groups, Communications on Pure and Applied Mathematics 71, 2275 (2018).
  * Zdeborová and Krzakala [2016] L. Zdeborová and F. Krzakala, Statistical physics of inference: Thresholds and algorithms, Advances in Physics 65, 453 (2016).
  * Ebner and Stroud [1985] C. Ebner and D. Stroud, Diamagnetic susceptibility of superconducting clusters: Spin-glass behavior, Phys. Rev. B 31, 165 (1985).
  * Ozeki and Nishimori [1993] Y. Ozeki and H. Nishimori, Phase diagram of gauge glasses, Journal of Physics A: Mathematical and General 26, 3399 (1993).
  * Nishimori [2001] H. Nishimori,  _Statistical physics of spin glasses and information processing: an introduction_ , 111 (Clarendon Press, 2001).
  * Iba [1999] Y. Iba, The nishimori line and bayesian statistics, Journal of Physics A: Mathematical and General 32, 3875 (1999).
  * Lupo and Ricci-Tersenghi [2017] C. Lupo and F. Ricci-Tersenghi, Approximating the xy model on a random graph with a q-state clock model, Physical Review B 95, 054433 (2017).
  * Morita  _et al._ [2006] S. Morita, Y. Ozeki, and H. Nishimori, Gauge theory for quantum spin glasses, Journal of the Physical Society of Japan 75, 014001 (2006).
  * Song and Zhang [2021] F.-F. Song and G.-M. Zhang, Tensor network approach to the two-dimensional fully frustrated xy model and a bosonic metallic phase with chirality order (2021), arXiv:2112.03550.
  * Lesieur  _et al._ [2017] T. Lesieur, F. Krzakala, and L. Zdeborová, Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications, Journal of Statistical Mechanics: Theory and Experiment 2017, 073403 (2017).
  * Thouless  _et al._ [1977] D. J. Thouless, P. W. Anderson, and R. G. Palmer, Solution of’solvable model of a spin glass’, Philosophical Magazine 35, 593 (1977).
  * Miolane [2017] L. Miolane, Fundamental limits of low-rank matrix estimation: the non-symmetric case, arXiv preprint arXiv:1702.00473 (2017).
  * Mossel  _et al._ [2015] E. Mossel, J. Neeman, and A. Sly, Reconstruction and estimation in the planted partition model, Probability Theory and Related Fields 162, 431 (2015).
  * Bolthausen [2014] E. Bolthausen, An iterative construction of solutions of the tap equations for the sherrington–kirkpatrick model, Communications in Mathematical Physics 325, 333 (2014).
  * Mézard  _et al._ [2002] M. Mézard, G. Parisi, and R. Zecchina, Analytic and algorithmic solution of random satisfiability problems, Science 297, 812 (2002).
  * Lucibello  _et al._ [2019] C. Lucibello, L. Saglietti, and Y. Lu, Generalized approximate survey propagation for high-dimensional estimation, in  _International Conference on Machine Learning_ (PMLR, 2019) pp. 4173–4182.
  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, physics, and computation_ (Oxford University Press, 2009).
  * Montanari  _et al._ [2004] A. Montanari, G. Parisi, and F. Ricci-Tersenghi, Instability of one-step replica-symmetry-broken phase in satisfiability problems, Journal of Physics A: Mathematical and General 37, 2073 (2004).
  * Rivoire  _et al._ [2004] O. Rivoire, G. Biroli, O. C. Martin, and M. Mézard, Glass models on bethe lattices, The European Physical Journal B-Condensed Matter and Complex Systems 37, 55 (2004).
  * Jammalamadaka and Sengupta [2001] S. R. Jammalamadaka and A. Sengupta,  _Topics in circular statistics_ , Vol. 5 (world scientific, 2001).
  * Bayati and Montanari [2011] M. Bayati and A. Montanari, The dynamics of message passing on dense graphs, with applications to compressed sensing, IEEE Transactions on Information Theory 57, 764 (2011).

##  Appendix A Circular distributions

In this appendix we recall some basic facts about probability measures on the
unit circle, and their connection to our setting. In the following we will
always assume that $x=e^{i\theta}$ is a complex variable of unit modulus. Let
$x\sim P(x)$, then we define the raw moments of $x$ as

$m_{n}=\oint_{\Gamma}dxe^{inx}P(x).$ (38)

In analogy with the linear case one can define the circular mean and variance
respectively as $m_{1}$ and $1-|m_{1}|$. We shall explore a family of circular
distributions that appear in the analysis of the planted XY model.

Suppose $x$ belongs to the following family of probability measures spanned by
the complex parameter $h=|h|e^{i\phi}$.

$P(x)=\frac{1}{2\pi I_{0}(|h|)}e^{\operatorname{Re}(x\overline{h})},$ (39)

where $I_{k}$ is the modified Bessel function of the first kind of order $k$.
The real part can now be written as
$\operatorname{Re}(x\overline{h})=|h|\cos(\theta-\phi)$, thus yielding a Von
Mises distribution [32] for the variable $\theta$:

$P(\theta)=\frac{1}{2\pi I_{0}(|h|)}e^{|h|\cos(\theta-\phi)}.$ (40)

The moments of $P$ are

$m_{n}=\operatorname{\mathbb{E}}_{P}[x^{n}]=\frac{I_{n}(|h|)}{I_{0}(|h|)}$
(41)

In fact

$\displaystyle 2\pi
I_{0}(|h|)m_{n}=\int_{0}^{2\pi}e^{in\theta}e^{|h|cos(\theta)}=\int_{0}^{2\pi}\left[\cos(n\theta)+i\sin(n\theta)\right]e^{|h|cos(\theta)}=\int_{0}^{2\pi}\cos(n\theta)e^{|h|cos(\theta)}=2\pi
I_{n}(|h|).$ (42)

The last equality follows from the definition of $I_{n}$.

##  Appendix B AMP and state evolution derivation

We shall now briefly recall the derivation of the AMP algorithm. While there
exist several ways to do so, we choose an approach similar to [22] based on
manipulating the Belief Propagation equations. Belief propagation for the
planted XY model reads

$\displaystyle m_{i\rightarrow ij}(x_{i})=\frac{1}{Z_{i\rightarrow
ij}}\prod_{k\neq i,j}m_{ki\rightarrow i}(x_{i})$ (43) $\displaystyle
m_{ij\rightarrow j}(x_{j})=\frac{1}{Z_{ij\rightarrow j}}\oint
dx_{i}\exp\left(2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}(Y_{ij}\overline{x_{i}}x_{j})\right)m_{i\rightarrow
ij}(x_{i}),$ (44)

where all the integrals are on the unit circle in the complex plane. This is a
set of $N^{2}$ functional equations: it would be impossible to use them in
practice on a computer. The first step to obtain AMP consists of relaxing BP:
this means finding a parametric form of the messages under which the BP
equations can be closed. We first expand the exponent in (44) equation:

$\displaystyle m_{ij\rightarrow j}(x_{j})=\frac{1}{Z_{ij\rightarrow j}}\oint
dx_{i}\left[1+2\sqrt{\frac{\hat{\lambda}}{N}}Re(Y_{ij}\overline{x_{i}}x_{j})+2\frac{\hat{\lambda}}{N}\operatorname{Re}(Y_{ij}\overline{x_{i}}x_{j})^{2}+O(N^{-3/2})\right]m_{i\rightarrow
ij}(x_{i})=$ (45) $\displaystyle=\frac{1}{Z_{ij\rightarrow j}}\oint
dx_{i}\left[1+2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}(Y_{ij}\overline{x_{i}}x_{j})+\frac{\hat{\lambda}}{N}\operatorname{Re}(Y_{ij}^{2}\overline{x_{i}}^{2}x_{j}^{2})+\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}+O(N^{-3/2})\right]m_{i\rightarrow
ij}(x_{i})=$ (46) $\displaystyle=\frac{1}{Z_{ij\rightarrow
j}}\left[1+2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}(Y_{ij}\overline{\hat{x}}_{i\rightarrow
ij}x_{j})+O(N^{-1})\right]=\frac{1}{Z_{ij\rightarrow
j}}\exp\left[2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}(Y_{ij}\overline{\hat{x}}_{i\rightarrow
ij}x_{j})\right]+O(N^{-1}).$ (47)

From the first to the second line we used that
$\operatorname{Re}(z)^{2}=\frac{1}{2}(\operatorname{Re}(z^{2})+|z|^{2})$. Then
from the second to the third line we dropped the $O(N^{-1})$ terms since these
are subdominant in the $N\to\infty$ limit and we performed the average
introducing $\hat{x}_{i\rightarrow
ij}=\operatorname{\mathbb{E}}_{m_{i\rightarrow ij}}[x_{i}]$. The normalization
constant can be computed $Z_{ij\rightarrow j}=2\pi
I_{0}\left(2\sqrt{\hat{\lambda}/N}|Y_{ij}\hat{x}_{i\rightarrow ij}|\right)$ By
substituting the last expression in (43) we get

$\displaystyle m_{i\rightarrow ij}(x_{i})=\frac{1}{Z_{i\rightarrow
ij}}\exp\left[2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}\left(x_{i}\sum_{k\neq
i,j}Y_{ki}\overline{\hat{x}}_{k\rightarrow ki}\right)\right]=\frac{1}{2\pi
I_{0}(h_{ij\rightarrow
j})}\exp\left[2\operatorname{Re}\left(x_{i}\overline{h}_{ij\rightarrow
j}\right)\right],$ (48)

Where in the last line we defined $h_{ij\rightarrow
j}=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq i,j}Y_{ik}\hat{x}_{k\rightarrow
ki}$. We’ve finally arrived to the relaxed BP-equations:

$\displaystyle h^{(t)}_{ij\rightarrow
j}=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i,j}Y_{ik}\hat{x}^{(t)}_{k\rightarrow ki}$ (49)
$\displaystyle\hat{x}_{i\rightarrow
ij}^{(t+1)}=\operatorname{\mathbb{E}}_{m^{(t)}_{i\rightarrow
ij}}[x_{i}]=\eta(h^{(t)}_{ij\rightarrow
j}),\qquad\eta(h)=\frac{h}{|h|}\frac{I_{1}(|2h|)}{I_{0}(|2h|)}$ (50)

For the computation of $\operatorname{\mathbb{E}}_{m_{i\rightarrow
ij}}[x_{i}]$ we defer to appendix A. To complete the derivation of AMP we now
remove the dependence of the target variable by expanding the relaxed-BP
equations. First define the single site fields as

$h_{i}=h_{ij\rightarrow j}+\sqrt{\hat{\lambda}/N}Y_{ij}\hat{x}_{j\rightarrow
ji}=\sqrt{\hat{\lambda}/N}\sum_{k\neq i}Y_{ik}\hat{x}_{k\rightarrow ki}$ (51)

Similarly we introduce the variables $\hat{x}$, which, at convergence,
represent AMP’s approximation to the system’s marginals.

$\hat{x}_{i}^{(t+1)}=\eta(h_{i}^{(t)})$ (52)

the goal is to replace variables indexed on edges with the new variables
indexed on vertices: to do so we have to keep track of the error

$\displaystyle\hat{x}^{(t)}_{k\rightarrow
ki}-\hat{x}_{k}^{(t)}=\eta(h^{(t-1)}_{ki\rightarrow
i})-\eta(h^{(t-1)}_{k})=-\frac{\partial\eta}{\partial
h}(h_{k}^{(t-1)})\sqrt{\frac{\hat{\lambda}}{N}}Y_{ki}\hat{x}^{(t-1)}_{i\rightarrow
ik}+O(N^{-1})=$ (53) $\displaystyle=-\frac{\partial\eta}{\partial
h}(h_{k}^{(t-1)})\sqrt{\frac{\hat{\lambda}}{N}}Y_{ki}\hat{x}_{i}^{(t-1)}+O(N^{-1})$

For a detailed explanation of the form of $\partial\eta/\partial h$ see the
paragraph B.0.1. Plugging this into $h_{i}^{(t)}$ we have

$\displaystyle h_{i}^{(t)}=2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\hat{x}_{k\rightarrow
ki}^{(t)}=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\left(\hat{x}_{k}^{(t)}-\sqrt{\frac{\hat{\lambda}}{N}}Y_{ki}\frac{\partial\eta}{\partial
h}(h_{k}^{(t-1)})\hat{x}_{i}^{(t-1)}\right)=$ (54)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\hat{x}_{k}^{(t)}-\frac{\hat{\lambda}}{N}\hat{x}_{i}^{(t-1)}\sum_{k\neq
i}|Y_{ik}|^{2}\frac{\partial\eta}{\partial h}(h_{k}^{(t-1)})$ (55)

Equation (52) together with (54) constitute the AMP algorithm.

####  B.0.1 Derivation of the Onsager term

In this section we derive the form of $\frac{\partial\eta}{\partial h}$.
Notice that $\eta$ is not an analytic function, so its derivative cannot be
expressed as a complex number, instead it takes the form of a $2\times 2$
Jacobian. Writing $\eta(h)=\eta(x+iy)=u(x,y)+iv(x,y)$, the directional
derivative of $\eta$ along $z$ is $\frac{\partial\eta}{\partial
z}=\operatorname{Re}(z)\left[\partial_{x}u(x,y)+i\partial
xv(x,y)\right]+\operatorname{Im}(z)\left[\partial_{u}u(x,y)+i\partial_{y}v(x,y)\right]$.
Decomposing $z$ along the directions respectively orthogonal and parallel to
$h$, and with the notation $r=|h|$ we obtain the following alternative
expression

$\eta(h+z)-\eta(h)=\left[z-\operatorname{Re}\left(\frac{\overline{h}}{|h|}z\right)\frac{h}{|h|}\right]\eta_{r}(r)+\operatorname{Re}\left(\frac{\overline{h}}{|h|}z\right)\frac{h}{|h|}\partial_{r}(r\eta_{r}(r))+O(|z|^{2})$
(56)

where we have defined $\eta_{r}(r)=\frac{1}{r}\frac{I_{1}(2r)}{I_{0}(2r)}$. In
the computation of the Onsager term we have

$\displaystyle h_{i}^{(t)}=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\eta(h_{k}^{(t-1)})-\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i}Y_{ik}\left[\eta(h_{k}^{(t-1)})-\eta\left(h_{k}^{(t-1)}-\sqrt{\frac{\hat{\lambda}}{N}}\overline{Y}_{ik}\hat{x}_{i}^{(t-1)}\right)\right]$
(57)

Let us fix $k$ in the second summation and look at one term: to lighten the
notation rename $h_{k}^{(t-1)}\mapsto h,\;Y_{ik}\mapsto
Y,\;\hat{x}_{i}^{(t-1)}\mapsto x$. Applying (56) we have

$\displaystyle
Y\left[\eta(h)-\eta\left(h-\sqrt{\frac{\hat{\lambda}}{N}}\overline{Y}x\right)\right]=\sqrt{\frac{\hat{\lambda}}{N}}Y\left[\overline{Y}x-\operatorname{Re}\left(\frac{\overline{h}}{|h|}\overline{Y}x\right)\frac{h}{|h|}\right]\eta_{r}(r)+\sqrt{\frac{\hat{\lambda}}{N}}Y\operatorname{Re}\left(\frac{\overline{h}}{|h|}\overline{Y}x\right)\frac{h}{|h|}\partial_{r}(r\eta_{r}(r))+O(1/N)=$
(58)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\left\\{|Y|^{2}x+Y\frac{h}{|h|}\operatorname{Re}\left(\frac{\overline{h}}{|h|}\overline{Y}x\right)\left[\partial_{r}(r\eta_{r}(r))-\eta_{r}(r)\right]\right\\}+O(1/N)=$
(59)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\left\\{\frac{1}{2}x|Y|^{2}\left(\partial_{r}(r\eta_{r}(r))+\eta_{r}(r)\right)+\frac{Y^{2}h^{2}\overline{x}}{2|h|^{2}}\left[\partial_{r}(r\eta_{r}(r))-\eta_{r}(r)\right]\right\\}+O(1/N),$
(60)

where in the last step we expanded
$\operatorname{Re}(z)=\frac{1}{2}(z+\overline{z})$. Substituting this back
into (57) and summing over $k$, one sees that the second term is negligible
because $Y^{2}$ is a zero mean random variable. Moreover the by the properties
of the Bessel functions we have the identity
$\frac{1}{2}\left(\partial_{r}(r\eta_{r}(r))+\eta_{r}(r)\right)=1-|\eta(r)|^{2}$.
To conclude (54) holds with

$\frac{\partial\eta}{\partial h}(h_{k}^{(t-1)})=1-|\eta(h_{k}^{(t-1)})|^{2}$
(61)

###  B.1 State evolution heuristic derivation

One of the elements which distinguishes AMP from other iterative algorithms is
the ability (in the $N\to\infty$ limit ) to track its dynamics through the
state evolution equations. In particular we will derive closed equations for
the two observables

$\displaystyle m=\frac{1}{N}\sum_{i=1}^{N}\hat{x}_{i}\qquad
q=\frac{1}{N}\sum_{i=1}^{N}|\hat{x}_{i}|^{2},$ (62)

representing respectively the alignment of the marginals with the planted
configuration and how concentrated each marginal is. We start by deriving an
iterative equation for $m^{(t)}$.

$\displaystyle
m^{(t+1)}=\frac{1}{N}\sum_{i=1}^{N}\eta\left(\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k=1}^{N}Y_{ik}\hat{x}_{k}^{(t)}\right)=\frac{1}{N}\sum_{i=1}^{N}\eta\left(\frac{\sqrt{\lambda\hat{\lambda}}}{N}\sum_{k=1}^{N}\hat{x}_{k}+\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k=1}^{N}W_{ik}\hat{x}_{k}\right)=$
(63)
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\eta\left(\sqrt{\lambda\hat{\lambda}}m^{(t)}+\sqrt{\frac{\hat{\lambda}}{2}q^{(t)}}z_{i}\right)=\frac{1}{2\pi}\int_{\mathbb{C}}dz\;e^{-\frac{1}{2}|z|^{2}}\eta\left(\sqrt{\lambda\hat{\lambda}}m^{(t)}+\sqrt{\hat{\lambda}q^{(t)}/2}\;z\right)$
(64)

In the second passage we defined
$z_{i}=\sqrt{\frac{2N}{q^{(t)}}}\sum_{k=1}^{N}W_{ik}\hat{x}_{k}\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$,
and in the end we replaced the sum over $i$ with an integral, since
$N\to\infty$ and $z_{i}$s are assumed to be independent. The whole derivation
revolves around the assumption of independence between $\hat{x}^{(t)}$ and
$W$. This is of course not the case, because $\hat{x}$ will depend on $W$
through previous iterations, however the Onsager term in the iterations re-
establishes asymptotic independence as explained in [12][33]. Similarly for
$q$ we have

$\displaystyle
q^{(t+1)}=\frac{1}{N}\sum_{i=1}^{N}\left|\eta\left(\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k=1}^{N}Y_{ik}\hat{x}_{k}^{(t)}\right)\right|^{2}=\frac{1}{2\pi}\int_{\mathbb{C}}dz\;e^{-\frac{1}{2}|z|^{2}}\left|\eta\left(\sqrt{\lambda\hat{\lambda}}m^{(t)}+\sqrt{\hat{\lambda}q^{(t)}/2}\;z\right)\right|^{2}.$
(65)

###  B.2 Simplification of state evolution

We can further simplify SE equations. The use of this is to reduce the number
of integrals to be done numerically. We start by simplifying $f_{m}$:

$\displaystyle
f_{m}(m,q)=\frac{1}{2\pi}\int_{\mathbb{C}}dz\;e^{-\frac{1}{2}|z|^{2}}\eta\left(\sqrt{\lambda\hat{\lambda}}m+\sqrt{\frac{\hat{\lambda}q}{2}}z\right)=\frac{1}{2\pi}\int
dx\,dy\;e^{-\frac{1}{2}(x^{2}+y^{2})}\frac{a(x)}{\sqrt{a(x)^{2}+\hat{\lambda}qy^{2}/2}}\tilde{\eta}\left(\sqrt{a(x)^{2}+\hat{\lambda}qy^{2}/2}\right).$
(66)

where we split $z=x+iy$ and we introduced
$a(x)=\sqrt{\lambda\hat{\lambda}}m+\sqrt{\hat{\lambda}q/2}x$. In principle
$f_{m}$ could be complex however the imaginary part is zero. By changing
variables according to $x^{\prime}=a(x),\;y^{\prime}=\sqrt{\hat{\lambda}q/2}y$
we get

$\displaystyle f_{m}(m,q)=\frac{1}{\pi\hat{\lambda}q}\int
dxdy\;\exp\left[-\frac{1}{\hat{\lambda}q}(y^{2}+(x-\sqrt{\lambda\hat{\lambda}}m)^{2})\right]\frac{x}{\sqrt{x^{2}+y^{2}}}\tilde{\eta}\left(\sqrt{x^{2}+y^{2}}\right)=$
(67) $\displaystyle=\frac{2}{\pi\hat{\lambda}q}e^{-\frac{\lambda
m^{2}}{q}}\int_{0}^{\infty}dr\int_{-r}^{r}dx\;\exp\left[-\frac{r^{2}}{\hat{\lambda}q}+2\frac{m}{q}\sqrt{\frac{\lambda}{\hat{\lambda}}}x\right]\frac{x}{\sqrt{r^{2}-x^{2}}}\tilde{\eta}\left(r\right)=$
(68) $\displaystyle=\frac{2}{\hat{\lambda}q}e^{-\frac{\lambda
m^{2}}{q}}\int_{0}^{\infty}dr\;r\exp\left[-\frac{r^{2}}{\hat{\lambda}q}\right]I_{1}\left(2\frac{m}{q}\sqrt{\frac{\lambda}{\hat{\lambda}}}r\right)\tilde{\eta}\left(r\right)$
(69)

where in the second line we changed variables to $x,r=\sqrt{x^{2}+y^{2}}$ and
in the last passage we used that
$\int_{-r}^{r}dx\,e^{ax}\frac{x}{\sqrt{r^{2}-x^{2}}}=\pi rI_{1}(ar)$. An
analogous procedure also yields a simplified equation for $q$:

$f_{q}(m,q)=\frac{2}{\hat{\lambda}q}e^{-\frac{\lambda
m^{2}}{q}}\int_{0}^{\infty}dr\;r\exp\left[-\frac{r^{2}}{\hat{\lambda}q}\right]I_{0}\left(2\frac{m}{q}\sqrt{\frac{\lambda}{\hat{\lambda}}}r\right)\tilde{\eta}\left(r\right)^{2}.$
(70)

##  Appendix C Fixed point analysis

We write SE equations in vectorial form as

$x^{t+1}=(m^{t+1},q^{t+1})=f(x^{t})=(f_{m}(m^{t},q^{t}),f_{q}(m^{t},q^{t})).$
(71)

We say $\tilde{x}$ is a fixed point if $f(\tilde{x})=\tilde{x}$. Let
$\tilde{x}$ be a fixed point and denote $\Delta=x-\tilde{x}$. Then at the
linear order in $f$ Delta will evolve as

$\Delta^{t+1}_{i}=\sum_{j}\frac{\partial f_{i}}{\partial
x_{j}}(\tilde{x})\Delta^{t}_{j}=(\left[J_{f}(\tilde{x})\right]^{t+1}\Delta^{0})_{i}$
(72)

In order to characterize the behavior of $\Delta$ (and hence study the
stability of $\tilde{x}$) we must then look at the Jacobian. In the following
we study this Jacobian for multiple fixed points.

###  C.1 $m=0,\;q=0$

We first evaluate the stability with respect to $m$.

$\displaystyle\frac{\partial f_{m}}{\partial
m}(m,q=0)\bigg{|}_{m=0}=\frac{\partial}{\partial
m}\frac{1}{2\pi}\int_{\mathbb{C}}dz\;e^{-\frac{1}{2}|z|^{2}}\eta\left(\sqrt{\lambda\hat{\lambda}}m\right)\bigg{|}_{m=0}=\sqrt{\hat{\lambda}\lambda}\eta^{\prime}(0)=\sqrt{\hat{\lambda}\lambda}.$
(73)

We now examine the stability with respect to $q$. For this purpose define
$\tilde{\eta}=|\eta|$.

$\displaystyle\frac{\partial f_{q}}{\partial
q}(q,m=0)\bigg{|}_{q=0}=\frac{\partial}{\partial q}\frac{1}{2\pi}\int
dz\;e^{-\frac{|z|^{2}}{2}}\tilde{\eta}\left(\sqrt{\hat{\lambda}q/2}\;z\right)^{2}\bigg{|}_{q=0}=\frac{\partial}{\partial
q}\int_{0}^{\infty}dr\;re^{-\frac{r^{2}}{2}}\tilde{\eta}\left(\sqrt{\hat{\lambda}q/2}\;r\right)^{2}\bigg{|}_{q=0}=$
(74)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{8q}}\int_{0}^{\infty}dr\;r^{2}e^{-\frac{r^{2}}{2}}2\tilde{\eta}\left(\sqrt{\hat{\lambda}q/2}\;r\right)\tilde{\eta}^{\prime}\left(\sqrt{\hat{\lambda}q/2}\;r\right)\bigg{|}_{q=0}=\sqrt{\frac{\hat{\lambda}}{2q}}\int_{0}^{\infty}dr\;r^{3}e^{-\frac{r^{2}}{2}}\tilde{\eta}^{\prime}\left(0\right)^{2}\sqrt{\hat{\lambda}q/2}\bigg{|}_{q=0}=$
(75)
$\displaystyle=\frac{\hat{\lambda}}{4}\int_{0}^{\infty}dr\;r^{3}e^{-\frac{r^{2}}{2}}=\hat{\lambda}$
(76)

In the second to last passage we use the fact that $q\to 0$ to expand
$\tilde{\eta}$ to the first order. Also remember that
$\tilde{\eta}^{\prime}(0)=1$.

So the $m=q=0$ fixed point is stable if $\hat{\lambda}<1$ and
$\hat{\lambda}\lambda<1$. This region is delimited by the curve
$\lambda=\min(1,\hat{\lambda}^{-1})$ shown in Figure 3.

###  C.2 $m=0$

In the spin glass phase we expect that $m=0$ while $q>0$. To find the boundary
with the phase where both $m$ and $q$ are positive we need to compute the
stability of $m=0$ for the value of $q$ given by the converged state
evolution.

Starting from (69) we compute $\partial f_{m}/\partial m$ for general $q$.

$\displaystyle\frac{\partial f_{m}}{\partial
m}=-2\lambda\frac{m}{q}f_{m}(m,q)+\frac{2}{q^{2}\hat{\lambda}}\sqrt{\frac{\lambda}{\hat{\lambda}}}e^{-\frac{\lambda
m^{2}}{q}}\int_{0}^{\infty}dr\;r^{2}\tilde{\eta}(r)e^{-\frac{r^{2}}{\hat{\lambda}q}}\left[I_{0}\left(2\frac{m}{q}\sqrt{\frac{\lambda}{\hat{\lambda}}}r\right)+I_{2}\left(2\frac{m}{q}\sqrt{\frac{\lambda}{\hat{\lambda}}}r\right)\right]$
(77)

By setting $m=0$ one obtains

$\displaystyle\frac{\partial f_{m}}{\partial
m}\bigg{|}_{m=0}=\frac{2\sqrt{\lambda}}{q^{2}\hat{\lambda}^{3/2}}\int_{0}^{\infty}dr\;r^{2}\tilde{\eta}(r)e^{-\frac{r^{2}}{\hat{\lambda}q}}$
(78)

##  Appendix D Convergence of AMP

In this appendix we study the regions of the $(\lambda,\hat{\lambda})$ plane
where AMP converges. First we introduce some quantities that will be required
for the analysis. For convenience we will sometimes treat complex numbers and
functions as vectors in $\mathbb{R}^{2}$.
$z=(\operatorname{Re}(z),\operatorname{Im}(z))$. Accordingly we will represent
$\eta(z)=\frac{z}{|z|}\frac{I_{1}(2|z|)}{I_{0}(2|z|)}$ as
$\eta(z)=(x\eta_{r}(r),y\eta_{r}(r))$, with $z=x+iy$, $r=\sqrt{x^{2}+y^{2}}$
and $\eta_{r}(r)=\frac{1}{r}\frac{I_{1}(2r)}{I_{0}(2r)}$. Since $\eta$ is not
an analytic function, its derivative cannot be expressed a a single complex
number, but the whole $2\times 2$ Jacobian $J\eta$ is required. We find

$J\eta(x,y)=\begin{bmatrix}\eta_{r}(r)+\frac{x^{2}}{r}\eta_{r}^{\prime}(r)&\frac{xy}{r}\eta_{r}^{\prime}(r)\\\
\frac{xy}{r}\eta_{r}^{\prime}(r)&\eta_{r}(r)+\frac{x^{2}}{r}\eta_{r}^{\prime}(r)\\\
\end{bmatrix}$ (79)

Finally we will need the following fact: let $M=\\{M_{ij}\\}_{i,j=1}^{2}$ be a
$2\times 2$ matrix and $\left\lVert\cdot\right\rVert_{2}$ be the euclidean
norm in $\mathbb{R}^{2}$. Then

$\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\left\lVert M\begin{bmatrix}\cos\theta\\\
\sin\theta\end{bmatrix}\right\rVert_{2}^{2}=\frac{1}{2}\left(M_{11}+M_{12}+M_{21}+M_{22}\right)$
(80)

We perturb the vector $h^{(t)}$ with an infinitesimal vector $\delta
h^{(t)}\in\mathbb{C}^{N}$, where coordinates are i.i.d. uniformly distributed
on a circle of radius $\epsilon$ (i.e. $\delta h^{(t)}_{k}=\epsilon
e^{i\phi_{k}^{(t)}},\;\;\phi_{k}^{(t)}\sim\text{Unif}([0,2\pi])$). Let
$\tilde{h}^{(t)}=h^{(t)}+\delta h^{(t)}$ be the perturbed vector.

If the perturbation grows in time then AMP will not converge since every fixed
point would be repulsive. Define the norm of the perturbation to be
$\left\lVert\delta h\right\rVert^{2}\coloneqq\frac{1}{N}\sum_{k=1}^{N}|\delta
h_{k}|^{2}$. In this way we have $\left\lVert\delta
h^{(t)}\right\rVert^{2}=\epsilon^{2}$. Using (8), (9) we get

$\displaystyle\delta h^{(t+1)}_{i}$
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k}Y_{ik}\left[\eta(h^{(t)}_{k}+\delta
h^{(t)}_{k})-\eta(h^{(t)}_{k})\right]+\frac{\hat{\lambda}}{N}\hat{x}_{i}^{(t-1)}\sum_{k}|Y_{ik}|^{2}\left[\frac{\partial\eta}{\partial
h}(h_{k}^{(t-1)}+\delta h^{(t-1)}_{k})-\frac{\partial\eta}{\partial
h}(h_{k}^{(t-1)})\right]=$ (81)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k}Y_{ik}J\eta(h^{(t)}_{k})\delta
h^{(t)}_{k}+\frac{\hat{\lambda}}{N}\hat{x}_{i}^{(t-1)}\sum_{k}|Y_{ik}|^{2}\left[\frac{\partial}{\partial
r}\left(\frac{\tilde{\partial\eta}}{\partial
h}(r)\right)\operatorname{Re}\left(\frac{\overline{h}_{k}^{(t-1)}}{|h_{k}^{(t-1)}|}\delta
h^{(t)}_{k}\right)\right]=$ (82)
$\displaystyle=\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k}Y_{ik}J\eta(h^{(t)}_{k})\delta
h^{(t)}_{k}+O(N^{-1/2}).$ (83)

In the second passage we introduced the notation $r=|h_{k}^{(t-1)}|$ and the
function $\frac{\tilde{\partial\eta}}{\partial
h}(|h|)\coloneqq\frac{\partial\eta}{\partial h}(h)$. Moreover we used the fact
that $|h+\delta h|-|h|=\operatorname{Re}(\frac{\overline{h}}{|h|}\delta
h)+O(|\delta h|^{2})$. In the last passage we exploited the fact that
$\operatorname{Re}(\frac{\overline{h}}{|h|}\delta h)$ has zero average with
respect to the randomness in $\delta h$, thus the Onsager term will give a
$N^{-1/2}$ contribution, which can be neglected. We now compute the norm of
the perturbation at time $t+1$

$\displaystyle\left\lVert\delta h^{(t+1)}\right\rVert^{2}$
$\displaystyle=\frac{1}{N}\sum_{i}\frac{\hat{\lambda}}{N}\sum_{k}|Y_{ik}|^{2}|J\eta(h^{(t)}_{k})\delta
h^{(t)}_{k}|^{2}+\frac{c_{1}}{N}=\frac{\hat{\lambda}}{N}\sum_{k}\left[\frac{1}{N}\sum_{i}|Y_{ik}|^{2}\right]|J\eta(h^{(t)}_{k})\delta
h^{(t)}_{k}|^{2}+\frac{c_{2}}{N}=$ (84)
$\displaystyle=\frac{\hat{\lambda}}{N}\sum_{k}|J\eta(h^{(t)}_{k})\delta
h^{(t)}_{k}|^{2}+\frac{c_{3}}{N}=\hat{\lambda}\int dh\,\int d\delta
h\,P_{\text{emp}}(h,\delta h)|J\eta(h)\delta h|^{2}+\frac{c_{3}}{N}=$ (85)
$\displaystyle=\epsilon^{2}\hat{\lambda}\int dh\,P(h)\frac{1}{2\pi}\int
d\theta|J\eta(h)e^{i\theta}|^{2}+\frac{c_{4}}{N}=$ (86)
$\displaystyle=\epsilon^{2}\frac{\hat{\lambda}}{2}\int
dh\,P(h)\left[\left(\eta_{r}(r)+\frac{x^{2}}{r}\eta_{r}^{\prime}(r)\right)^{2}+2\left(\frac{xy}{r}\eta_{r}^{\prime}(r)\right)^{2}+\left(\eta_{r}(r)+\frac{y^{2}}{r}\eta_{r}^{\prime}(r)\right)^{2}\right]+\frac{c_{4}}{N}=$
(87) $\displaystyle=\epsilon^{2}\frac{\hat{\lambda}}{2}\int
dh\,P(h)\left[\eta_{r}(r)^{2}+\left(\eta_{r}(r)+r\eta_{r}^{\prime}(r)\right)^{2}\right]+\frac{c_{4}}{N}$
(88)

Going from the first to the second line we used the fact that
$\operatorname{\mathbb{E}}_{Y}[|Y|^{2}]=1+O(1/N)$ and that the empirical
average over $i$ concentrates around the expected value. In the second line we
also introduced the empirical joint distribution of $\delta h$ and $h$:
$P_{\text{emp}}(h,\delta
h)=\frac{1}{N}\sum_{k}\delta(h-h^{(t)}_{k})\delta(\delta h-\delta
h^{(t)}_{k})$. Since $h^{(t)}$ and $\delta h^{(t)}$ are independent, in the
limit of $N\to\infty$, $P_{\text{emp}}(h,\delta h)\to P(h)P(\delta h)$. $P(h)$
is determined by SE’s prediction that
$h=\sqrt{\lambda\hat{\lambda}m}+\sqrt{\frac{\hat{\lambda}q}{2}}z$, with
$z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$, while $P(\delta h)$ is, coordinate
wise, the uniform distribution on the circle of radius $\epsilon$. To obtain
the third line we used (80) and (79), with $h=x+iy,\,r=|h|$.
$c_{1},c_{2},c_{3},c_{4}$ are all constants with respect to $N$.

To conclude, the perturbation norm obeys $\left\lVert\delta
h^{(t+1)}\right\rVert^{2}=(1-c^{(t)})\left\lVert\delta
h^{(t)}\right\rVert^{2}$, with

$\displaystyle
c^{(t)}=1+\frac{\hat{\lambda}}{2}\operatorname{\mathbb{E}}_{z}\left[\eta_{r}(|h|)^{2}+\left(\eta_{r}(|h|)+|h|\eta_{r}^{\prime}(|h|)\right)^{2}\right]$
(89) $\displaystyle
h=\sqrt{\lambda\hat{\lambda}m^{(t)}}+\sqrt{\frac{\hat{\lambda}q^{(t)}}{2}}z,\quad
z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$ (90)

with $m^{(t)},q^{(t)}$ obtained iterating SE. This result is valid with high
probability with respect to $\delta h,Y$. To study the convergence of single
instances of AMP it is useful to derive the average (with respect to $\delta
h$) growth of a perturbation, when $N$ is finite. Following an analogous
derivation we obtain $\operatorname{\mathbb{E}}_{\delta
h}\left[\left\lVert\delta
h^{(t+1)}\right\rVert^{2}\right]=(1-c_{\text{AMP}}^{(t)})\operatorname{\mathbb{E}}_{\delta
h}\left[\left\lVert\delta h^{(t)}\right\rVert^{2}\right]$, with

$c_{\text{AMP}}^{(t)}=1+\frac{\hat{\lambda}}{2}\sum_{i=1}^{N}\left[\frac{1}{N}\sum_{k}|Y_{ik}|^{2}\right]\left[\eta_{r}(|h_{i}^{(t)}|)^{2}+\left(\eta_{r}(|h_{i}^{(t)}|)+|h_{i}^{(t)}|\eta_{r}^{\prime}(|h_{i}^{(t)}|)\right)^{2}\right]$
(91)

##  Appendix E Correspondence of 1RSB free entropy and replicated free entropy

We pick

$s_{\star}=\begin{cases}1&\text{ if $\Sigma(f^{*}(1))\geq 0$}\\\
-\Sigma^{\prime}(f_{0})&\text{ if $\Sigma(f^{*}(1))<0$ }\end{cases}$ (92)

such that
$\tilde{\Phi}_{\text{1RSB}}=\Sigma(f^{*}(s_{\star}))+f^{*}(s_{\star})$ for the
chosen $s_{\star}$. Let us verify that this choice is correct: Suppose
$f^{*}(s_{\star})=\tilde{f}_{\star}$, then the argument goes through because
$\tilde{\Phi}_{\text{1RSB}}=\Sigma(\tilde{f}^{*})+\tilde{f}^{*}=\Sigma(f^{*}(s_{\star}))+f^{*}(s_{\star})$.
Hence, we only have to show $f^{*}(s_{\star})=\tilde{f}_{\star}$. If
$\Sigma(f^{*}(1))>0$, then $s_{\star}=1$ and $\tilde{f}^{*}=f^{*}(1)$.
Instead, if $\Sigma(f^{*})<0$, we will have $\tilde{f}^{*}=f_{0}$. But with
the new choice of $s_{\star}$, $f^{*}$ will satisfy
$\Sigma^{\prime}(f^{*})=\Sigma^{\prime}(f_{0})$, hence also giving
$f^{*}(s_{\star})=f_{0}$. Therefore $f^{*}(s_{\star})=\tilde{f}_{\star}$.
$\square$

##  Appendix F Derivation of Approximate Survey Propagation and its state
evolution

We derive the ASP algorithm as follows. First, we start from the Belief
Propagation (BP) equations for a replicated model consisting of $s$
independent replicas. We then put forward an ansatz for the messages,
parametrizing them by their means and covariances. Propagating the ansatz
through the BP equations, we are able to obtain a closed set of equations for
the parameters. This procedure is general and applicable to other statistical
physics models with a Boltzmann probability measure, but we will stick to the
explicit form of Hamiltonian for our problem for clarity.

The BP equations are

$\displaystyle{m}_{{i}\rightarrow{}{ij}}(\vec{x}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{i\to
ij}}\prod_{a=1}^{s}P_{X}\left(x_{i}^{(a)}\right)\prod_{k\neq
j}{m}_{{ki}\rightarrow{}{i}}(\vec{x}_{i})$ (93)
$\displaystyle{m}_{{ij}\rightarrow{}{i}}(\vec{x}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{ij\to i}}\int
d\vec{x}_{j}{m}_{{j}\rightarrow{}{ij}}(\vec{x}_{j})\exp\left[\sum_{a=1}^{s}2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right],$
(94)

where $\vec{x}_{i}=\left(x_{i}^{(1)},x_{i}^{(2)},\dots,x_{i}^{(s)}\right)$
refers to the variable $x_{i}$ in different replicas, and
$P_{X}\left(x_{i}^{(a)}\right)=\delta(|x_{i}^{(a)}|-1)/2\pi$ is the prior
distribution.

The ansatz for the messages is a Gaussian distribution parametrized by the
covariance of variables within and among replicas

${m}_{{j}\rightarrow{}{ij}}(\vec{x}_{j})=\frac{1}{Z_{j\to ij}}\int
dh\exp\left[-\frac{1}{2}\frac{|h-{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}}{\Delta_{j\to
ij}^{2}}\right]\prod_{a=1}^{s}\exp\left[-\frac{1}{2}|h-x_{j}^{(a)}|^{2}\right].$
(95)

The ansatz produces the following correlation functions:

$\displaystyle\langle x_{j}^{(a)}\rangle$ $\displaystyle=$
$\displaystyle{\hat{x}}_{{j}\rightarrow{}{ij}}$ (96) $\displaystyle\langle
x_{j}^{(a)}\overline{x}_{j}^{(b)}\rangle_{a\neq b}$ $\displaystyle=$
$\displaystyle{\Delta}_{{j}\rightarrow{}{ij}}+|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}$
(97) $\displaystyle\langle x_{j}^{(a)}\overline{x}_{j}^{(a)}\rangle$
$\displaystyle=$ $\displaystyle 1\quad\text{by the unit norm constraint}$ (98)

where $\langle\cdot\rangle$ denotes the average with respect to the
distribution described by (95). The precise form (95) is of little importance,
our analysis will only make use of the correlation identities defined above.
The next step is to expand

$\exp\left[\sum_{a=1}^{s}2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right]=1+2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{a=1}^{s}\mathrm{Re}\left(Y_{ij}\overline{x}_{i}^{(a)}x_{j}^{(a)}\right)+\frac{\lambda}{N}|Y_{ij}|^{2}\sum_{a,b}^{s}\overline{x}_{i}^{(a)}x_{i}^{(b)}x_{j}^{(a)}\overline{x}_{j}^{(b)}+O(N^{-3/2}),$
(99)

and express ${m}_{{ij}\rightarrow{}{j}}(\vec{x}_{i})$ in terms of the
estimator ${\hat{x}}_{{i}\rightarrow{}{ij}}$ and covariance
${\Delta}_{{j}\rightarrow{}{ij}}$. In the following computations the symbol
$\doteq$ will mean that the two sides are equal up to terms that vanish in the
$N\to\infty$ limit.

$\displaystyle{m}_{{ij}\rightarrow{}{i}}(\vec{x}_{i})=\frac{1}{Z^{\prime\prime}_{ij\to
i}}\int
d\vec{x}_{j}{m}_{{j}\rightarrow{}{ij}}(\vec{x}_{j})\exp\left[\sum_{a=1}^{s}2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right]\doteq$
$\displaystyle\doteq$ $\displaystyle\frac{1}{Z^{\prime}_{ij\to
i}}\left\\{1+2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}{\hat{x}}_{{j}\rightarrow{}{ij}}\sum_{a=1}^{s}\overline{x}_{i}^{(a)}\right)+\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}({\Delta}_{{j}\rightarrow{}{ij}}+|\hat{x}_{j\to
ij}|^{2})\left|\sum_{a=1}^{s}x_{i}^{(a)}\right|^{2}+\right.$ $\displaystyle+$
$\displaystyle\frac{\hat{\lambda}}{N}s\left(1-{\Delta}_{{j}\rightarrow{}{ij}}-|\hat{x}_{j\to
ij}|^{2}\right)-\left.\frac{\hat{\lambda}}{N}\operatorname{Re}\left(Y_{ij}^{2}\hat{x}_{j\to
ij}^{2}\left(\sum_{a=1}^{s}\overline{x}_{i}^{(a)}\right)^{2}\right)\right\\}\doteq$
$\displaystyle\doteq$ $\displaystyle\frac{1}{Z_{ij\to
i}}\exp\left[2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}{\hat{x}}_{{j}\rightarrow{}{ij}}\sum_{a=1}^{s}\overline{x}_{i}^{(a)}\right)-\frac{\hat{\lambda}}{N}\operatorname{Re}\left[Y_{ij}^{2}\hat{x}_{j\to
ij}^{2}\right]\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}{\Delta}_{{j}\rightarrow{}{ij}}\left|\sum_{a=1}^{s}x_{i}^{(a)}\right|^{2}\right].$
(101)

Where all constants not depending on $x_{i}$ have been absorbed into $Z_{ij\to
i}$. Moreover we dropped the last term in (F) because it is subdominant in
$N$. Moving to the first BP equation (93) we have

$\displaystyle m_{i\to ij}(\vec{x}_{i})$ $\displaystyle\doteq$
$\displaystyle\frac{1}{Z_{i\to
ij}}\left(\prod_{a=1}^{s}P_{X}(x_{i}^{(a)})\right)\exp\left[2\sqrt{\frac{\hat{\lambda}}{N}}\operatorname{Re}\left(\left(\sum_{k\neq
i,j}Y_{ik}\hat{x}_{k\to
ik}\right)\sum_{a=1}^{s}\overline{x}_{i}^{(a)}\right)+\right.$
$\displaystyle+$
$\displaystyle\left.\frac{\hat{\lambda}}{N}\left|\sum_{a=1}^{s}\overline{x}_{i}^{(a)}\right|^{2}\sum_{k\neq
i,j}|Y_{ik}|^{2}\Delta_{k\to ik}\right]$

At this point it is convenient to have only the linear term of
$\sum_{a}\overline{x}_{i}^{(a)}$ in the exponent so that we can factorise over
different replicas. To achieve this, we apply the Hubbard-Stratonovich trick

$\int_{\mathbb{C}}d{x}\,\frac{\theta}{\pi}e^{-\theta(|x|^{2}-2\mathrm{Re}[\overline{x}y])}=e^{\theta|y|^{2}}.$
(103)

Picking $y=\sum_{a=1}^{s}x_{i}^{(a)}$ and $\theta=\frac{\lambda}{N}\sum_{k\neq
i,j}|Y_{ik}|^{2}{\Delta}_{{k}\rightarrow{}{ik}}$ we get

$\displaystyle m_{i\to ij}(\vec{x}_{i})$ $\displaystyle\doteq$
$\displaystyle\frac{1}{Z_{i\to
ij}}\int_{\mathbb{C}}dh\exp\left(-|h|^{2}\frac{\hat{\lambda}}{N}\sum_{k\neq
i,j}|Y_{ik}|^{2}{\Delta}_{{k}\rightarrow{}{ik}}\right)\times$
$\displaystyle\times$
$\displaystyle\prod_{a=1}^{s}\left\\{P_{X}(x_{i}^{(a)})\exp\left[2\operatorname{Re}\left(\overline{x}_{i}^{(a)}\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
i,j}Y_{ik}\hat{x}_{k\to
ik}+\overline{x}_{i}^{(a)}h\frac{\hat{\lambda}}{N}\sum_{k\neq
i,j}|Y_{ik}|^{2}\Delta_{k\to ik}\right)\right]\right\\}$

where $h$ is a complex variable and the integral is computed over the entire
complex plane. The advantage of writing the message in the form of (F) is that
one can use it as a measure to explicitly evaluate the correlation functions,
and since now the replicas are properly factorised, evaluation then follows
readily. We define the following functions for brevity,

$\displaystyle
P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})$
$\displaystyle=$ $\displaystyle
P_{X}(x_{i})\exp\left[2\mathrm{Re}\left({T}_{{i}\rightarrow{}{ij}}\overline{x}_{i}+{V}_{{i}\rightarrow{}{ij}}h\overline{x}_{i}\right)\right]$
(105) $\displaystyle{T}_{{i}\rightarrow{}{ij}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\hat{\lambda}}{N}}\sum_{k\neq
j}Y_{ik}{\hat{x}}_{{k}\rightarrow{}{ik}}$ (106)
$\displaystyle{V}_{{i}\rightarrow{}{ij}}$ $\displaystyle=$
$\displaystyle\frac{\hat{\lambda}}{N}\sum_{k\neq
j}|Y_{ik}|^{2}{\Delta}_{{k}\rightarrow{}{ik}}$ (107)

with which the message is concisely expressed as

$\displaystyle{m}_{{i}\rightarrow{}{ij}}(\vec{x}_{i})$ $\displaystyle\doteq$
$\displaystyle\frac{1}{Z_{i\to
ij}}\int_{\mathbb{C}}d{h}\,e^{-{V}_{{i}\rightarrow{}{ij}}|h|^{2}}\prod_{a=1}^{s}P_{\mathrm{p}}(x_{i}^{(a)},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}}).$
(108)

At this stage, the equations are closed, and we have at hand the following
expressions

$\displaystyle\langle x_{i}^{(a)}\rangle$ $\displaystyle=$
$\displaystyle{\hat{x}}_{{i}\rightarrow{}{ij}}$ (109) $\displaystyle\doteq$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{{i}\rightarrow{}{ij}}|h|^{2}}\left[\int
d{x}_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})\right]^{s-1}\int
d{x}_{i}x_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{{i}\rightarrow{}{ij}}|h|^{2}}\left[\int
d{x}_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})\right]^{s}}$
$\displaystyle\langle x_{i}^{(a)}\overline{x}_{i}^{(b)}\rangle$
$\displaystyle=$
$\displaystyle{\Delta}_{{i}\rightarrow{}{ij}}+|{\hat{x}}_{{i}\rightarrow{}{ij}}|^{2}$
(110) $\displaystyle\doteq$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{{i}\rightarrow{}{ij}}|h|^{2}}\left[\int
d{x}_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})\right]^{s-2}\left|\int
d{x}_{i}x_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})\right|^{2}}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{{i}\rightarrow{}{ij}}|h|^{2}}\left[\int
d{x}_{i}P_{\mathrm{p}}(x_{i},h,{T}_{{i}\rightarrow{}{ij}},{V}_{{i}\rightarrow{}{ij}})\right]^{s}}.$

Note that ${\Delta}_{{i}\rightarrow{}{ij}}$ is strictly non-negative as
ensured by the Cauchy–Schwarz inequality. It is important to keep in mind that
${\Delta}_{{i}\rightarrow{}{ij}}$ has to be non-negative for the integral in
(108) to converge.

After the TAPyfication procedure, consisting in removing the dependence on the
target node in $m_{i\to ij}$ at the price of introducing an Onsager term, we
arrive at the iterative ASP algorithm with marginal distributions. Here we
denote the various quantities computed at iteration $t$ with a superscript
$t$,

$\displaystyle{T}_{i}^{t}$ $\displaystyle=$
$\displaystyle\sum_{k}\sqrt{\frac{\hat{\lambda}}{N}}Y_{ik}{\hat{x}}_{k}^{t}-{\hat{x}}_{i}^{t-1}\sum_{k}\frac{\hat{\lambda}}{N}|Y_{ik}|^{2}\frac{d{\hat{x}}_{k}}{dT_{k}}$
(111) $\displaystyle{V}_{i}^{t}$ $\displaystyle=$
$\displaystyle\sum_{k}\frac{\hat{\lambda}}{N}|Y_{ik}|^{2}{\Delta}_{k}^{t}$
(112) $\displaystyle{\hat{x}}_{i}^{t+1}$ $\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[\int
d{x}\,P_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})\right]^{s-1}\int
d{x}\,xP_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[\int
d{x}P_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})\right]^{s}}$ (113)
$\displaystyle{\Delta}_{i}^{t+1}+|{\hat{x}}_{i}^{t+1}|^{2}$ $\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[\int
d{x}\,P_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})\right]^{s-2}\left|\int
d{x}\,xP_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})\right|^{2}}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[\int
d{x}\,P_{\mathrm{p}}(x,h,{T}_{i}^{t},{V}_{i}^{t})\right]^{s}}.$ (114)

The integrals are not new. The prior restricts the range of integration to be
on the unit circle and then the integrals evaluate to modified Bessel
functions of the first kind $I_{k}$, defined by

$\displaystyle
I_{k}(x)=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\,\cos{(k\theta)}\exp{(x\cos{\theta})}.$
(115)

Explicitly,

$\displaystyle{\hat{x}}_{i}^{t+1}$ $\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}^{t}+V_{i}h|)\right]^{s-1}\left[\frac{T_{i}^{t}+V_{i}^{t}h}{|T_{i}^{t}+V_{i}h|}I_{1}(2|T_{i}^{t}+V_{i}^{t}h|)\right]}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}^{t}+V_{i}^{t}h|)\right]^{s}}$
(116) $\displaystyle{\Delta}_{i}^{t+1}+|{\hat{x}}_{i}^{t+1}|^{2}$
$\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}^{t}+V_{i}^{t}h|)\right]^{s-2}\left[I_{1}(2|T_{i}^{t}+V_{i}^{t}h|)\right]^{2}}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2|T_{i}^{t}+V_{i}^{t}h|)\right]^{s}}$
(117)

The form is not quite convenient for us to implement an algorithm since one of
the parameters $T_{i}$ is of complex nature, and integration with such a
parameter is more expensive computationally. We notice, however, that the
phase of $T_{i}$ can be factored out, if we write
$T_{i}=|T_{i}|e^{i\theta_{T_{i}}}$ and change the integration variable $h$
accordingly to $he^{i\theta_{T_{i}}}$, so conveniently we have

$\displaystyle{\hat{x}}_{i}^{t+1}$ $\displaystyle=$ $\displaystyle
e^{i\theta_{T_{i}^{t}}}\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]^{s-1}\left[\frac{|T_{i}^{t}|+V_{i}^{t}h}{||T_{i}^{t}|+V_{i}^{t}h|}I_{1}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]^{s}}$
(118) $\displaystyle{\Delta}_{i}^{t+1}+|{\hat{x}}_{i}^{t+1}|^{2}$
$\displaystyle=$
$\displaystyle\frac{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]^{s-2}\left[I_{1}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]^{2}}{\int_{\mathbb{C}}d{h}\,e^{-{V}_{i}^{t}|h|^{2}}\left[I_{0}(2||T_{i}^{t}|+V_{i}^{t}h|)\right]^{s}}$
(119)

The integrals now depend on two real parameters, $V_{i}^{t}$ and
$|T_{i}^{t}|$, suitable to approximate with interpolation method( see Appendix
H). To continue, we need another approximation in which we take
$d\hat{x}_{k}/dT_{k}$ to be real. The justification is as follows. Since the
gradient direction of $\frac{d\hat{x}_{k}}{dT_{k}}$ along $\delta
T_{k}=\sqrt{\frac{\lambda}{N}}Y_{ki}\delta\hat{x}_{i}+O(1/N)$ is mostly
uncorrelated to $T_{k}$ due to large system size $N$, we can take the average
$\left<\frac{d\hat{x}_{k}}{dT_{k}}\right>_{\theta_{\delta
T_{k}}}=\mathbb{E}_{\theta_{\delta T_{k}}}\left[\lim_{\delta T_{k}\rightarrow
0}\frac{\delta\hat{x}_{k}}{\delta T_{k}}\right]$ over a uniformly distributed
angle $\theta_{\delta T_{k}}$ to be the value of
$\frac{d\hat{x}_{k}}{dT_{k}}$, which will be always real.

The state evolution is then written as follows. Consider in the general case
where the estimated $\hat{\lambda}$ is different from the true value of
$\lambda$,

$\displaystyle
T^{t}_{i}=\sum_{k}\sqrt{\frac{\hat{\lambda}}{N}}Y_{ik}{\hat{x}^{t}}_{{k}\rightarrow{}{ki}}$
$\displaystyle=$
$\displaystyle\sqrt{\lambda\hat{\lambda}}\left(\sum_{k}\frac{{\hat{x}^{t}}_{{k}\rightarrow{}{ki}}}{N}\right)x_{i}+\sqrt{\frac{\hat{\lambda}}{N}}W_{ik}{\hat{x}^{t}}_{{k}\rightarrow{}{ki}}$
(120) $\displaystyle=$
$\displaystyle\sqrt{\lambda\hat{\lambda}}m^{t}x_{i}+\sqrt{\hat{\lambda}q^{t}/2}z$
$\displaystyle V^{t}$ $\displaystyle=$ $\displaystyle\hat{\lambda}\Delta^{t}$
(121) $\displaystyle m^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[\hat{x}(T^{t},V^{t})\right]$ (122)
$\displaystyle q^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[|\hat{x}(T^{t},V^{t})|^{2}\right]$ (123)
$\displaystyle\Delta^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\left[\Delta(T^{t},V^{t},q^{t})\right]$ (124)

where $z$ is a complex variable distributed as
$z\sim\mathcal{N}(0,1)+i\mathcal{N}(0,1)$. Here, we have simplified the
expression with $\Delta^{t}=\frac{1}{N}\sum_{i}\Delta_{i}^{t}$. In deriving
$T_{i}^{t}$ we have referred to the expression before TAPyfication such that
we can use the variable ${\hat{x}^{t}}_{{k}\rightarrow{}{ki}}$, whose
correlation with the noise $W_{ik}$ we ignore [33].

##  Appendix G Free entropy computation

We compute the Bethe free entropy for the general case of $s\neq 1$, following
the recipe provided in [29][14]. Then setting $s=1$ we will recover the RS
free entropy. The replicated 1RSB free entropy (20) is simply the Bethe free
entropy in the $s-$replicated graphical model. To obtain a consistent
expression of the free entropy, we will resurrect some necessary terms that
were absorbed as normalisation factors for the posterior distribution and
messages. Particularly, the following form of the posterior measure is used

$P(x|Y)=\frac{1}{Z}\exp{\left[2\sqrt{\frac{\hat{\lambda}}{N}}\sum_{i<j}\sum_{a=1}^{s}\mathrm{Re}\left(Y_{ij}\overline{x}_{i}^{(a)}x_{j}^{(b)}\right)-\frac{\hat{\lambda}}{N}\sum_{i\leq
j}\sum_{a=1}^{s}|(x_{i}^{(a)})^{2}||(x_{j}^{(a)})^{2}|\right]}\prod_{i=1}^{N}\prod_{a=1}^{s}P_{X}(x_{i}^{(a)})$
(125)

where $P_{X}=\delta(|x_{i}|-1)/(2\pi)$ is the prior distribution, which for
simplicity we keep in implicit form. Notice also that since in practice the
prior $P_{X}$ forces the spins to be of norm one, the second term in the
exponential is just a constant. Exploiting the fact that our measure
factorizes according to a pairwise graphical model we write

$\displaystyle N\Phi_{\text{1RSB}}=\log Z=\sum_{i=1}^{N}\log
Z^{i}-\sum_{i,j:i<j}\log Z^{ij}$ (126) $\displaystyle Z^{i}=\int
d\vec{x}_{i}\,P_{X}\left(\vec{x}_{i}\right)\prod_{j}\int
d\vec{x}_{j}\exp\left[-s\frac{\hat{\lambda}}{N}+\sum_{a=1}^{s}2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right]{m}_{{j}\rightarrow{}{ij}}(\vec{x}_{j})$
(127) $\displaystyle Z^{ij}=\int
d\vec{x}_{i}d\vec{x}_{j}\,\exp\left[-s\frac{\hat{\lambda}}{N}+\sum_{a=1}^{s}2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}\overline{x}^{(a)}_{i}x^{(a)}_{j}\right)\right]{m}_{{j}\rightarrow{}{ji}}(\vec{x}_{j}){m}_{{i}\rightarrow{}{ij}}(\vec{x}_{i})$
(128)

First we find $Z^{i}$. We use the moment identities (96) to perform the inner
averages over $m_{j\to ij}$. When averaging we take a step very similar to (F)
but without absorbing the integration constant terms into the normalisation
factor.

$\displaystyle
Z^{i}\doteq\exp{\left[\sum_{j}s\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-{\Delta}_{{j}\rightarrow{}{ij}}-|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}\prod_{j\neq
i}\int d\vec{x}_{i}\,P_{X}(\vec{x}_{i})\times$
$\displaystyle\times\exp{\left[-s\frac{\hat{\lambda}}{N}+2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}{\hat{x}}_{{j}\rightarrow{}{ij}}\sum_{a=1}^{s}\overline{x}^{(a)}_{i}\right)+\frac{\hat{\lambda}}{N}\left|Y_{ij}^{2}\right|{\Delta}_{{j}\rightarrow{}{ij}}\left|\sum_{a=1}^{s}x_{i}^{(a)}\right|^{2}\right]},$
(129)

and under trick (103) we get,

$\displaystyle Z^{i}$ $\displaystyle\doteq$
$\displaystyle\exp{\left[-s\hat{\lambda}+\sum_{j}s\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-{\Delta}_{{j}\rightarrow{}{ij}}-|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}\frac{\frac{\hat{\lambda}}{N}\sum_{j=1}^{N}|Y_{ij}|^{2}{\Delta}_{{j}\rightarrow{}{ij}}}{\pi}\int_{\mathbb{C}}dh\,e^{-\frac{\hat{\lambda}}{N}\sum_{j}|Y_{ij}|^{2}{\Delta}_{{j}\rightarrow{}{ij}}|h|^{2}}\times$
(130) $\displaystyle\times\int
d\vec{x}_{i}\,P_{X}(\vec{x}_{i})\exp{\left\\{\sum_{j=1}^{N}2\mathrm{Re}\left[\left(\sqrt{\frac{\hat{\lambda}}{N}}Y_{ij}{\hat{x}}_{{j}\rightarrow{}{ij}}+\frac{\hat{\lambda}}{N}\left|Y_{ij}^{2}\right|{\Delta}_{{j}\rightarrow{}{ij}}h\right)\sum_{a=1}^{s}\overline{x}^{(a)}_{i}\right]\right\\}}$
$\displaystyle\doteq$
$\displaystyle\exp{\left[-s\hat{\lambda}+\sum_{j=1}^{N}s\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-{\Delta}_{{j}\rightarrow{}{ij}}-|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}\frac{V_{i}}{\pi}\int_{\mathbb{C}}dhe^{-V_{i}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s},$
(131)

where in the last step we took the unit norm prior into account. The
expression for $Z^{ij}$ can then be simplified: we proceed by expanding the
exponential, performing the averages and then re expressing the result in
exponential form. With that we have the final set of equations to compute the
Bethe free entropy,

$\displaystyle Z^{i}$ $\displaystyle=$
$\displaystyle\exp{\left[-s\hat{\lambda}+\sum_{j}s\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-{\Delta}_{{j}\rightarrow{}{ij}}-|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}\frac{V_{i}}{\pi}\int_{\mathbb{C}}dhe^{-V_{i}|h|^{2}}\left[I_{0}(2|T_{i}+V_{i}h|)\right]^{s}$
(132) $\displaystyle Z^{ij}$ $\displaystyle=$
$\displaystyle\exp\left[-s\frac{\hat{\lambda}}{N}+2s\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}{\overline{\hat{x}}}_{{i}\rightarrow{}{ij}}{\hat{x}}_{{j}\rightarrow{}{ij}}\right)\right.+$
(133)
$\displaystyle\left.+s\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}\left(1+(s-1)\left({\Delta}_{{i}\rightarrow{}{ij}}+|{\hat{x}}_{{i}\rightarrow{}{ij}}|^{2}\right)\left({\Delta}_{{j}\rightarrow{}{ij}}+|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}\right)-s|{\hat{x}}_{{i}\rightarrow{}{ij}}|^{2}|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}\right)\right]$

In the special case of $s=1$, we start computation from (127) and (128),
bearing in mind that there is only one replica, therefore cross terms of the
form $\sum_{a\neq b}x^{(a)}x^{(b)}$ disappear. Consequently, terms involving
$\Delta$ also disappear. The expression for $Z^{(i)}$ and $Z^{ij}$ are then
both greatly simplified:

$\displaystyle Z^{i}|_{s=1}$ $\displaystyle=$
$\displaystyle\exp{\left[-\hat{\lambda}+\sum_{j}\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}I_{0}\left(2|T_{i}|\right)$
(134) $\displaystyle Z^{ij}|_{s=1}$ $\displaystyle=$
$\displaystyle\exp{\left[-\frac{\hat{\lambda}}{N}+2\sqrt{\frac{\hat{\lambda}}{N}}\mathrm{Re}\left(Y_{ij}{\hat{\overline{x}}}_{{i}\rightarrow{}{ij}}{\hat{x}}_{{j}\rightarrow{}{ij}}\right)+\frac{\hat{\lambda}}{N}|Y_{ij}|^{2}(1-|{\hat{x}}_{{i}\rightarrow{}{ij}}|^{2}|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2})\right]}.$
(135)

From the expressions for $Z^{i}$ and $Z^{(ij)}$, using (126) we have the
single instance free entropy is
$\Phi_{\text{1RSB}}(Y,\hat{\lambda},\hat{\lambda},s)$. To compute the free
entropy we just plug in the values of $\\{\hat{x}_{i\to ij}\\}$ and
$\\{\Delta_{i\to ij}\\}$ to which the ASP algorithm converges. In the
$N\to\infty$ limit we expect
$\Phi_{\text{1RSB}}(Y,\hat{\lambda},\hat{\lambda},s)$ to concentrate around
its mean value
$\Phi_{\text{1RSB}}(\lambda,\hat{\lambda},s)=\mathbb{E}_{Y}[\Phi_{\text{1RSB}}(Y,\hat{\lambda},\hat{\lambda},s)]$.
We prove this concentration for some of the terms appearing in
$\Phi_{\text{1RSB}}$:

$\displaystyle\frac{1}{N}\sum_{j=1}^{N}|Y_{ij}|^{2}=\frac{\lambda}{N}+\frac{1}{N}\sum_{j=1}^{N}|W_{ij}|^{2}+\frac{\sqrt{\lambda}}{N^{3/2}}\operatorname{Re}\left(\sum_{j=1}^{N}W_{ij}\right)\rightsquigarrow
1$ (136)
$\displaystyle\frac{1}{N}\sum_{j}|Y_{ij}|^{2}|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}=\frac{\lambda}{N}q+\frac{1}{N}\sum_{j=1}^{N}|W_{ij}|^{2}|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}+2\frac{\sqrt{\lambda}}{N^{3/2}}\sum_{j=1}^{N}\operatorname{Re}(W_{ij})|{\hat{x}}_{{j}\rightarrow{}{ij}}|^{2}\rightsquigarrow\frac{1}{N}\sum_{j=1}^{N}|\hat{x}_{j}|^{2}=q$
(137)
$\displaystyle\frac{1}{N}\sqrt{\frac{\hat{\lambda}}{N}}\sum_{i<j}\mathrm{Re}\left(Y_{ij}{\overline{\hat{x}}}_{{i}\rightarrow{}{ij}}{\hat{x}}_{{j}\rightarrow{}{ij}}\right)=\frac{\sqrt{\lambda\hat{\lambda}}}{2N^{2}}\sum_{i\neq
j}\operatorname{Re}\left({\overline{\hat{x}}}_{{i}\rightarrow{}{ij}}{\hat{x}}_{{j}\rightarrow{}{ij}}\right)+\frac{\hat{\lambda}}{N^{3/2}}\sum_{i<j}\operatorname{Re}\left(W_{ij}{\overline{\hat{x}}}_{{i}\rightarrow{}{ij}}{\hat{x}}_{{j}\rightarrow{}{ij}}\right)\rightsquigarrow$
(138)
$\displaystyle\rightsquigarrow\frac{\sqrt{\lambda\hat{\lambda}}}{2N^{2}}\operatorname{Re}\left(\sum_{i,j}\overline{\hat{x}_{i}}\hat{x}_{j}\right)=\frac{1}{2}\sqrt{\lambda\hat{\lambda}}m^{2}.$
(139)

Here $\rightsquigarrow$ indicated convergence in probability with respect to
$W$. We also recall that $({\hat{x}}_{{j}\rightarrow{}{ij}},\Delta_{i\to
ij})\to(\hat{x}_{j},\Delta_{i})$. Moreover we assume that the correlation
between $W_{ij}$ and ${\hat{x}}_{{j}\rightarrow{}{ij}}$ is asymptotically
canceled by the Onsager term, so we can consider them independent. Similarly,
we can express other terms using $m$, $q$ and $\Delta$. For the terms that
involve $T_{i}$ and $V_{i}$ (or $h_{i}$ in the case of AMP), we can use their
values obtained by state evolution. In the end the averaged replicated free
entropy reads

$\displaystyle\Phi_{\mathrm{\text{1RSB}}}(s,\lambda,\hat{\lambda})$
$\displaystyle=$
$\displaystyle-s\sqrt{\lambda\hat{\lambda}}m^{2}+\frac{s\lambda}{2}\left(sq^{2}-2(\Delta+q)-(s-1)(\Delta+q)^{2}\right)+$
$\displaystyle+\mathbb{E}_{z}\left[\log\left(\frac{\lambda\Delta}{\pi}\int_{\mathbb{C}}dh\exp(-\lambda\Delta|h|^{2})[I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\lambda
q/2}z+\lambda\Delta h|)]^{s}\right)\right]$

The states selected by $s$ have free entropy

$\displaystyle f^{*}(s,\lambda,\hat{\lambda})$ $\displaystyle=$
$\displaystyle\frac{\partial\Phi_{\mathrm{\text{1RSB}}}(s,\lambda,\hat{\lambda})}{\partial
s}=-\sqrt{\lambda\hat{\lambda}}m^{2}+s\lambda
q^{2}-\lambda(q+\Delta)-\frac{(2s-1)\lambda}{2}(\Delta+q)^{2}$
$\displaystyle+$
$\displaystyle\mathbb{E}_{z}\left[\frac{\int_{\mathbb{C}}dh\exp(-\lambda\Delta|h|^{2})\log[I_{0}(2|T+\lambda\Delta
h|)][I_{0}(2|T+\lambda\Delta
h|)]^{s}}{\int_{\mathbb{C}}dh\exp(-\lambda\Delta|h|^{2})[I_{0}(2|T+\lambda\Delta
h|)]^{s}}\right]$

and correspondingly from (134) we get the following expression for AMP,

$\displaystyle\Phi_{\mathrm{RS}}(\lambda,\hat{\lambda})$ $\displaystyle=$
$\displaystyle-\sqrt{\lambda\hat{\lambda}}m^{2}+\frac{\lambda}{2}\left(q^{2}-2q\right)+\mathbb{E}_{z}\left[\log\left(I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\lambda
q/2}z|)\right)\right],$ (143)

where $\Phi_{\mathrm{RS}}$ denotes the Bethe free entropy when $s=1$.
Interestingly also setting $\Delta=0$ in (G) one recovers
$\Phi_{\mathrm{RS}}$, independently of the value of $s$.

###  G.1 Recovering state evolution Equations for AMP and ASP

In this section we show that the fixed points of the SE equations are
stationary points of the replicated free entropy. To verify this, we check the
first derivatives of the Bethe free entropy with respect to parameters of the
system. For AMP, the parameters are $m$ and $q$,

$\displaystyle\frac{\partial\Phi_{\mathrm{RS}}}{\partial m}$ $\displaystyle=$
$\displaystyle-2\sqrt{\lambda\hat{\lambda}}m+\mathbb{E}_{z}\left[\frac{\partial_{m}\oint
dx\,\exp{\left[2\mathrm{Re}(h\overline{x})\right]}}{2\pi I_{0}(2|h|)}\right]$
(144) $\displaystyle=$
$\displaystyle-2\sqrt{\lambda\hat{\lambda}}m+2\sqrt{\lambda\hat{\lambda}}\mathrm{Re}\\{\mathbb{E}_{z}\left[\frac{\oint
dx\overline{x}\,\exp{[2\mathrm{Re}(h\overline{x})]}}{2\pi
I_{0}(2|h|)}\right]\\}$ $\displaystyle=$
$\displaystyle-2\sqrt{\lambda\hat{\lambda}}m+2\sqrt{\lambda\hat{\lambda}}\mathrm{Re}\left\\{\mathbb{E}_{z}\left[\,\overline{\hat{x}}\,\right]\right\\}$
$\displaystyle\frac{\partial\Phi_{\mathrm{RS}}}{\partial q}$ $\displaystyle=$
$\displaystyle\lambda(q-1)+\mathbb{E}_{z}\left[\frac{\partial_{q}\oint
dx\exp{[2\mathrm{Re}(h\overline{x})]}}{2\pi I_{0}(2|h|)}\right]$ (145)
$\displaystyle=$ $\displaystyle\lambda(q-1)+\frac{\lambda}{\sqrt{2\lambda
q}}\mathrm{Re}\left\\{\mathbb{E}_{z}\left[\frac{\oint
dx\,z\overline{x}\exp{[2\mathrm{Re}(h\overline{x})]}}{2\pi
I_{0}(2|h|)}\right]\right\\}$ $\displaystyle=$
$\displaystyle\lambda(q-1)+\lambda\mathrm{Re}\left\\{\mathbb{E}_{z}\left[\frac{\oint
dx\,|x|^{2}\exp{[2\mathrm{Re}(h\overline{x})]}}{2\pi I_{0}(2|h|)}-\frac{|\oint
dx\,\overline{x}\exp{[2\mathrm{Re}(h\overline{x})]}|^{2}}{[2\pi
I_{0}(2|h|)]^{2}}\right]\right\\}$ $\displaystyle=$
$\displaystyle\lambda(q-1)+\lambda\mathbb{E}_{z}\left[1-|\hat{x}|^{2}\right]$

These two derivatives vanish under the SE equations given in eq.(11). When
deriving $\Phi_{RS}$ with respect to $q$, we have used Stein’s lemma and taken
derivatives with respect to $z_{1}$ and $z_{2}$. Note that the quantity $h$
here is consistent with the expression in eq.(13).

For ASP, similarly, we have the derivatives of $\Phi_{\mathrm{Bethe}}$ with
respect to $m$ and $q$, computed in a similar manner,

$\displaystyle\frac{\partial\Phi_{\mathrm{Bethe}}}{\partial m}$
$\displaystyle=$
$\displaystyle-2s\sqrt{\lambda\hat{\lambda}}m+\mathbb{E}_{z}\left[\frac{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\partial_{m}\left[I_{0}(2|T+Vh|)\right]^{s}}{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\left[I_{0}(2|T+Vh|)\right]^{s}}\right]$
(146) $\displaystyle=$
$\displaystyle-2s\sqrt{\lambda\hat{\lambda}}m+\mathbb{E}_{z}\left[\frac{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}I_{1}(2|T+Vh|)\left[I_{0}(2|T+Vh|)\right]^{s-1}2s\sqrt{\lambda\hat{\lambda}}\mathrm{Re}\left((T+Vh)\overline{x}_{0}\right)}{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\left[I_{0}(2|T+Vh|)\right]^{s}}\right]$
$\displaystyle=$
$\displaystyle-2s\sqrt{\lambda\hat{\lambda}}m+2s\sqrt{\lambda\hat{\lambda}}\mathrm{Re}\\{\mathbb{E}_{z}[\hat{x}]\\}$
$\displaystyle\frac{\partial\Phi_{\mathrm{Bethe}}}{\partial q}$
$\displaystyle=$ $\displaystyle s^{2}\lambda
q-s\lambda-s\lambda(s-1)(\Delta+q)+$ (147)
$\displaystyle+s\sqrt{\frac{\lambda}{2q}}\mathbb{E}_{z}\left[\frac{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\left[I_{0}(2|T+Vh|)\right]^{s-1}\int
dx\,\mathrm{Re}[z\overline{x}]\exp{2\mathrm{Re}[(T+Vh)\overline{x}]}}{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\left[I_{0}(2|T+Vh|)\right]^{s}}\right]$
$\displaystyle=$ $\displaystyle s^{2}\lambda
q-s\lambda-s\lambda(s-1)(\Delta+q)+$
$\displaystyle+s\lambda\mathrm{Re}\\{\mathbb{E}_{z}\left[(s-1)\frac{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}[I_{0}(2|T+Vh|]^{s-2}[I_{1}(2|T+Vh|)]^{2}}{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}\left[I_{0}(2|T+Vh|)\right]^{s}}+\right.$
$\displaystyle\left.+1-s\left(\frac{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}[I_{0}(2|T+Vh|)]^{s-1}I_{1}(2|T+Vh|)}{\int_{\mathbb{C}}dh\,e^{-V|h|^{2}}[I_{0}(2|T+Vh|)]^{s}}\right)^{2}\right]\\}$
$\displaystyle=$ $\displaystyle s^{2}\lambda
q-s\lambda-s\lambda(s-1)(\Delta+q)+s\lambda\mathrm{Re}\\{\mathbb{E}_{z}\left[(s-1)(\Delta+q)+1-sq\right].$

In addition, we also have the derivative with respect to $\Delta$, which
measures the overlap between different replica. To make the computation
simpler, we first perform a change of variable to the integral that appears in
$\Phi_{\mathrm{Bethe}}$,

$\displaystyle\frac{\lambda\Delta}{\pi}\int_{\mathbb{C}}dh\,\exp{(-\lambda\Delta|h|^{2})}\left[I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\lambda
q/2}z+\lambda\Delta h|)\right]^{s}$ (148) $\displaystyle=$
$\displaystyle\frac{\lambda}{\pi}\int_{\mathbb{C}}dh\,\exp{(-\lambda|h|^{2})}\left[I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\lambda
q/2}z+\lambda\sqrt{\Delta}h|)\right]^{s}.$

The integral can be seen as the expectation value of the function
$\left[I_{0}(2|\sqrt{\lambda\hat{\lambda}}m+\sqrt{\lambda
q/2}z+\lambda\sqrt{\Delta}h|)\right]^{s}$ under two Gaussian distributions of
the real and imaginary part of $h$, with mean $\mu=0$ and variance
$\sigma^{2}=\frac{1}{2\lambda}$.

$\displaystyle\frac{\partial\Phi_{\mathrm{Bethe}}}{\partial\Delta}$
$\displaystyle=$
$\displaystyle-s\lambda-s\lambda(s-1)(\Delta+q)+\mathbb{E}_{z}\left[\frac{\partial_{\Delta}\mathbb{E}_{h}\left\\{\left[\int
dx\,\exp{\left(2\mathrm{Re}[(T+Vh)\overline{x}]\right)}\right]^{s}\right\\}}{\mathbb{E}_{h}\left\\{\left[\int
dx\,\exp{\left(2\mathrm{Re}[(T+Vh)\overline{x}]\right)}\right]^{s}\right\\}}\right]$
(149) $\displaystyle=$ $\displaystyle-s\lambda-s\lambda(s-1)(\Delta+q)+$
$\displaystyle+\mathbb{E}_{z}\left[\frac{\frac{s\lambda}{\sqrt{\Delta}}\mathbb{E}_{h}\left\\{I_{0}(2|T+Vh|)^{s-1}\int
dx\,\mathrm{Re}[h\overline{x}]\exp{(2\mathrm{Re}[(T+Vh)\overline{x}])}\right\\}}{\mathbb{E}_{h}\left\\{\left[I_{0}(2|T+Vh|)\right]^{s}\right\\}}\right]$
$\displaystyle=$ $\displaystyle-s\lambda-s\lambda(s-1)(\Delta+q)+$
$\displaystyle+\mathbb{E}_{z}\left[\frac{\mathbb{E}_{h}\left\\{s\lambda(s-1)\left[I_{0}(2|T+Vh|)\right]^{s-2}\left[I_{1}(2|T+Vh|)\right]^{2}+s\lambda\left[I_{0}(2|T+Vh|)\right]^{s}\right\\}}{\mathbb{E}_{h}\left\\{\left[I_{0}(2|t+Vh|)\right]^{s}\right\\}}\right]$
$\displaystyle=$
$\displaystyle-s\lambda-s\lambda(s-1)(\Delta+q)+\mathbb{E}_{z}\left[s\lambda(s-1)(\Delta+q)+s\lambda\right]$

where we again used Stein’s formula on the real and imaginary parts of $h$.

At the equilibrium point, the first derivatives should vanish and we retrieve
the state evolution equations

$\displaystyle m=\mathbb{E}_{z}\left[{\hat{\overline{x}}}\right]$ (150)
$\displaystyle q=\mathbb{E}_{z}\left[|\hat{x}|^{2}\right]$ (151)
$\displaystyle\Delta=\mathbb{E}_{z}\left[\Delta\right]$ (152)

##  Appendix H Numerical details

###  H.1 Interpolating integral functions

We encounter a lot of integral functions in the updating steps, especially in
the ASP case where the integral dimension is larger compared to AMP.
Therefore, we use numerical interpolation of integrals Eq.(118,119) to speed
up both the ASP algorithm and state evolution. After try and error, the
interpolation (linear) on grid $(s,\log_{10}(T),\log_{10}(V))$ that we are
working with is:

  * •

$\log_{10}(T)$ in range $[-7,2]$ with 800 evenly spaced points,

  * •

$\log_{10}(V)$ in range $[-7,2]$ with 800 evenly spaced points,

  * •

$s$ in range $[0,1]$ with any grid configuration that one needs, and the grid
extensively used to generate plots shown in the manuscript is chosen in the
range $[0,1]$ with evenly spaced points.

Such interpolation grid was shown to give SE results with precision on the
order of ${10}^{-5}\text{\,}$. The most important parameters here are the
resolutions of $T$ and $V$, and also their range that should cover all the
parameter space of $T$ and $V$ during any iteration run. Such a precision is
the best we can do with reasonable file size and computation speed. The
achieved precision is good enough for retrieving most of the quantities we
study in this manuscript, but generally not enough to resolve complexity
functions at the level of ${10}^{-8}\text{\,}$.

###  H.2 Resolving complexity functions

Because the complexity functions encountered in the manuscript can be as small
as on the order of ${10}^{-8}\text{\,}$, the interpolation setting mentioned
in H.1 is not enough to resolve them. To obtain enough numerical precision to
resolve complexity functions, we use the following strategies:

All the derivative operations are done analytically so only numerical
integration is performed to ensure enough precision. The evaluations of
Eq.(G,G) are not done through interpolation, since they only need to be
evaluated once, and it’s not possible to achieve reasonable precision beyond
${10}^{-5}\text{\,}$ with interpolation methods. Instead, they are numerically
integrated directly using quadrature methods, with precision beyond
${10}^{-8}\text{\,}$. To obtain free energy fixed points with enough
precision, we first run SE with interpolated function until it converges at
the error level of ${10}^{-5}\text{\,}$ due to its fast computation speed.
Then we run SE with direct integration using quadrature method for 100 more
steps (much slower but can finish in a reasonable amount of time, targeting a
precision of fixed points on the level of ${10}^{-8}\text{\,}$), reducing the
error further and approach closer to the true fixed points. The comparison
with the AMP result indicates that this SE run method using a combination of
interpolated functions, and numerical integration of Eq.(G,G,118,119) using
quadrature methods, can generally achieved ${10}^{-8}\text{\,}$ numerical
precision on the results shown in Fig. 7 and 8, correctly revealing the
complexity functions on the order of ${10}^{-8}\text{\,}$.