Stefano Sarao Mannelli

Giulio Biroli

Chiara Cammarota

Florent Krzakala

Pierfrancesco Urbani

Lenka Zdeborová

Marvels and Pitfalls of the Langevin Algorithm in Noisy High-dimensional Inference

Gradient-descent-based algorithms and their stochastic versions have widespread applications in machine learning and statistical inference. In this work we perform an analytic study of the performances of one of them, the Langevin algorithm, in the context of noisy high-dimensional inference. We employ the Langevin algorithm to sample the posterior probability measure for the spiked matrix-tensor model. The typical behaviour of this algorithm is described by a system of integro-differential equations that we call the Langevin state evolution, whose solution is compared with the one of the state evolution of approximate message passing (AMP). Our results show that, remarkably, the algorithmic threshold of the Langevin algorithm is sub-optimal with respect to the one given by AMP. We conjecture this phenomenon to be due to the residual glassiness present in that region of parameters. Finally we show how a landscape-annealing protocol, that uses the Langevin algorithm but violate the Bayes-optimality condition, can approach the performance of AMP.

Motivation

Algorithms based on noisy variants of gradients descent [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF][START_REF] Welling | Bayesian learning via stochastic gradient Langevin dynamics[END_REF] stand at the roots of many modern applications of data science, and are being used in a wide range of high-dimensional non-convex optimization problems. The widespread use of stochastic gradient descent in deep learning [START_REF] Lecun | Deep learning[END_REF] is certainly one of the most prominent examples. For such algorithms, the existing theoretical analysis mostly concentrate on convex functions, convex relaxations or on regimes where spurious local minima become irrelevant. For problems with complicated landscapes where, instead, useful convex relaxations are not known and spurious local minima cannot be ruled out, the theoretical understanding of the behaviour of gradient-descent-based algorithm remains poor and represents a major avenue of research.

The goal of this paper is to contribute to such an understanding in the context of statistical learning, and to transfer ideas and techniques developed for glassy dynamics [START_REF] Bouchaud | Out of equilibrium dynamics in spin-glasses and other glassy systems[END_REF] to the analysis of non-convex relaxations in high-dimensional inference. In statistical learning, the minimization of a cost function is not the goal per se, but rather a way to uncover an unknown structure in the data. One common way to model and analyze this situation is to generate data with a hidden structure, and to see if the structure can be recovered. This is easily set up as a teacher-student scenario [START_REF] Seung | Statistical mechanics of learning from examples[END_REF][START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF]: First a teacher generates latent variables and uses them as input of a prescribed model to generate a synthetic dataset. Then, the student observes the dataset and tries to infer the values of the latent variables. The analysis of this setting has been carried out rigorously in a wide range of teacher-student models for high-dimensional inference and learning tasks as diverse as planted clique [START_REF] Deshpande | Finding hidden cliques of size (N/e) in nearly linear time[END_REF], generalized linear models such as compressed sensing or phase retrieval [START_REF] Barbier | Phase transitions, optimal errors and optimality of message-passing in generalized linear models[END_REF], factorization of matrices and tensors [START_REF] Barbier | Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF] or simple models of neural networks [START_REF] Aubin | The committee machine: Computational to statistical gaps in learning a two-layers neural network[END_REF]. In these works, the information theoretically optimal performances -the one obtained by an ideal Bayes-optimal estimator, not limited in time and memory-have been computed.

The main question is, of course, how practical algorithms -operating in polynomial time with respect to the problem size-compare to these ideal performances. The last decade brought remarkable progress into our understanding of the performances achievable computationally. In particular, many algorithms based on message passing [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF], spectral methods [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], and semidefinite programs (SDP) [START_REF] Samuel | Bayesian estimation from few samples: community detection and related problems[END_REF] were analyzed. Depending on the signal-to-noise ratio, these algorithms were shown to be very efficient in many of those task. Interestingly, all these algorithm fail to reach good performance in the same region of the parameter space, and this striking observation has led to the identification of a well-defined hard phase. This is a regime of parameters in which the underlying statistical problem can be information-theoretically solved, but no efficient algorithms are known, rendering the problem essentially unsolvable for large instances. This stream of ideas is currently gaining momentum and impacting research in statistics, probability, and computer science.

The performance of the noisy-gradient descent algorithms -that are certainly the one currently most used in practice-remains an entirely open question. Do they allow to reach the same performances as message passing and SDPs? Can they enter the hard phase, do they stop to be efficient at the same moment as the other approaches, or are they worse? The ambition of the present paper is to address these questions by analyzing the performance of the Langevin algorithm in the high-dimensional limit of a spiked matrix-tensor model, defined in detail in the next section.

We argue that this spiked matrix-tensor problem is both generic and relevant. Similar models have played a fundamental role in statistics and random matrix theory [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF][START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF]. Tensor factorization is also an important topic in mathematics and machine learning [START_REF] Richard | A statistical model for tensor PCA[END_REF][START_REF] Samuel B Hopkins | Tensor principal component analysis via sum-ofsquare proofs[END_REF][START_REF] Ge | On the optimization landscape of tensor decompositions[END_REF][START_REF] Ben Arous | Algorithmic thresholds for tensor PCA[END_REF][START_REF] Ros | Complex energy landscapes in spiked-tensor and simple glassy models: ruggedness, arrangements of local minima and phase transitions[END_REF], and is a widely used in data analysis [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF]. At variance with the pure spiked tensor case [START_REF] Richard | A statistical model for tensor PCA[END_REF], this mixed matrix-tensor model has the advantage that algorithmic threshold appears at the same scale as the information-theoretic one, similarly to what is observed in simple models of neural networks [START_REF] Barbier | Phase transitions, optimal errors and optimality of message-passing in generalized linear models[END_REF][START_REF] Aubin | The committee machine: Computational to statistical gaps in learning a two-layers neural network[END_REF].

In the present paper, we analyse the behavior of a noisy gradient descent algorithm, also called Langevin algorithm. We explicitly compare its performance to the one of the Bayes optimal estimator and to the best known efficient algorithm so-far -the approximate message passing algorithm [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF]. In particular, contrary to what has been anticipated in [START_REF] Krzakala | Hiding quiet solutions in random constraint satisfaction problems[END_REF][START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF], but as surmised in [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF], we observe that the performance of the Langevin algorithm is hampered by the many spurious metastable states still present in the AMP-easy phase. We show that the Langevin algorithm can approach AMP performances if one uses a landscape-annealing protocol in which instead of reducing the temperature, as done in thermal annealing, one instead anneals the strength of the contribution of the tensor. A number of intriguing conclusions can be drawn by these results and are likely to be valid beyond the considered model.

Finally, the possibility to describe analytically the behavior of the Langevin algorithm in this model is enabled by the existence of the Crisanti-Horner-Sommers-Cugliandolo-Kurchan (CHSCK) equations in spin glass theory, describing the behavior of the Langevin dynamics in the so-called spherical p-spin model [START_REF] Crisanti | The spherical p-spin interaction spin-glass model[END_REF][START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF], where the method can be rigorously justified [START_REF] Ben Arous | Cugliandolo-Kurchan equations for dynamics of spin-glasses[END_REF]. These equations were a key development in the field of statistical physics of disordered systems that lead to detailed understanding and predictions about the slow dynamics of glasses [START_REF] Bouchaud | Out of equilibrium dynamics in spin-glasses and other glassy systems[END_REF]. In this paper, we bring these powerful methods and ideas into the realm of statistical learning.

The spiked matrix-tensor model

We now detail the spiked matrix-tensor problem: a teacher generates a N -dimensional vector x * by choosing each of its components independently from a normal Gaussian distribution of zero mean and unit variance. In the large N limit this is equivalent to have a flat distribution over the N -dimensional hypersphere S N -1 defined by |x * | 2 = N . The teacher then also generates a symmetric matrix Y ij and a symmetric order-p tensor T i 1 ,...,ip as

Y ij = 1 √ N x * i x * j + ξ ij ∀i < j , T i 1 ...ip = (p -1)! N (p-1)/2 x * i 1 . . . x * ip + ξ i 1 ...ip ∀i 1 < . . . < i p , (1)
where ξ ij and ξ i 1 ,...,ip are iid Gaussian components of a symmetric random matrix and tensor of zero mean and variance ∆ 2 and ∆ p , respectively; ξ ij and ξ i 1 ,...,ip correspond to noises corrupting the signal of the teacher. In the limit ∆ 2 → 0, and ∆ p → 0, the above model reduces to the canonical spiked Wigner model [START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF], and spiked tensor model [START_REF] Richard | A statistical model for tensor PCA[END_REF], respectively. The goal of the student is to infer the vector x * from the knowledge of the matrix Y , of the tensor T , of the values ∆ 2 and ∆ p , and the knowledge of the spherical prior. The scaling with N as specified in Eq. ( 1) is chosen in such a way that the information-theoretically best achievable error varies between perfectly reconstructed spike x * and random guess from the flat measure on S N -1 . Here, and in the rest of the paper we denote x ∈ S N -1 the N -dimensional vector, and x i with i = 1, . . . , N its components. Analogous matrix-tensor models, where next to a order-p tensor one observes a matrix created from the same spike are studied e.g. in [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF] in the context of topic modeling, or in [START_REF] Richard | A statistical model for tensor PCA[END_REF]. From the optimization-theory point of view, this model is highly non-trivial being high-dimensional and non-convex. For the purpose of the present paper this model is chosen because it involves three key ingredients: (a) It is in the class of models for which the Langevin algorithm can be analyzed exactly in large N limit. (b) The different phase transitions, both algorithmic and information theoretic, discussed hereafter, all happen at ∆ 2 = O(1), ∆ p = O(1). (c) The AMP algorithm is in this model conjectured to be optimal among polynomial algorithms. It is this second and third ingredient that are not present in the pure spiked tensor model [START_REF] Richard | A statistical model for tensor PCA[END_REF], making it unsuitable for our present study. We note that the Langevin algorithm was recently analyzed for the pure spiked tensor model in [START_REF] Ben Arous | Algorithmic thresholds for tensor PCA[END_REF] in a regime where the noise variance is very small ∆ ∼ N -p/2 , but we also note that in that model algorithms such as tensor unfolding and semidefinite programming work better, roughly up to ∆ ∼ N -p/4 [START_REF] Richard | A statistical model for tensor PCA[END_REF][START_REF] Samuel B Hopkins | Tensor principal component analysis via sum-ofsquare proofs[END_REF],

3 Bayes-optimal estimation and message-passing

In this section we present the performance of the Bayes-optimal estimator and of the approximate message passing algorithm. This theory is based on a straightforward adaptation of analogous results known for the pure spiked matrix model [START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF][START_REF] Barbier | Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula[END_REF] and for the pure spiked tensor model [START_REF] Richard | A statistical model for tensor PCA[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF].

The Bayes-optimal estimator x is defined as the one that among all estimators minimizes the mean-squared error (MSE) with the spike x * . Starting from the posterior probability distribution

P (x|Y, T ) = 1 Z(Y, T )   N i=1 e -x 2 i /2   i<j e -1 2∆ 2 Y ij - x i x j √ N 2 i 1 <•••<ip e -1 2∆p T i 1 ...ip - √ (p-1)! N (p-1)/2 x i 1 ...x ip 2 , (2)
the Bayes-optimal estimator reads xi = E P (x|Y,T ) (x i ) .

To simplify notation, and to make contact with the energy landscape and the statistical physics notations, it is convenient to introduce the energy cost function, or Hamiltonian, as

H(x) = H 2 + H p = - 1 ∆ 2 √ N i<j Y ij x i x j - (p -1)! ∆ p N (p-1)/2 i 1 <•••<ip T i 1 ...ip x i 1 . . . x ip (4)
so that keeping in mind that for N → ∞ the spherical constraint is satisfied |x| 2 = N , the posterior is written as

P (x|Y, T ) = exp[-H(x)]/ Z(Y, T )
, where Z is the normalizing partition function.

With the use of the replica theory and its recent proofs from [START_REF] Barbier | Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula[END_REF][START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF] one can establish rigorously that the mean squared error achieved by the Bayes-optimal estimator is given as MMSE = 1 -m * where m * ∈ R is the global maximizer of the so-called free entropy of the problem

Φ RS (m) = 1 2 log(1 -m) + m 2 + 1 + m 2 4∆ 2 + 1 + m p 2p∆ p . ( 5
)
This expression is derived, and proven, in the appendix Sec. B.1. We now turn to the approximate message-passing (AMP) [START_REF] Richard | A statistical model for tensor PCA[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF], that is the best known so far for this problem. AMP is an iterative algorithm inspired from the work of Thouless-Anderson and Palmer in statistical physics [START_REF] David J Thouless | Solution of'solvable model of a spin glass[END_REF]. We explicit its form in the Sec. B.1. Most remarkably performance of AMP can be evaluated by tracking its evolution with the iteration time and it is given in terms of the (possibly local) maximum of the above free entropy that is reached as a fixed point of the following iterative process In the easy (green) region AMP achieves the optimal error smaller than random pick from the prior. In the impossible region (red) the optimal error is as bad as random pick from the prior, and AMP achieves it as well. In the hard region (orange) the optimal error is low, but AMP does not find an estimator better than random pick from the prior. In the case of Langevin algorithm the performance is strictly worse than that for AMP in the sense that the hard region increases up to line depicted in green dots. The blue dashed-dotted line delimits the region of existence of stable 1RSB states.

m t+1 = 1 - 1 1 + m t /∆ 2 + (m t ) p-1 /∆ p (6)
with initial condition m t=0 = with 0 < 1. Eq. ( 6) is called the State Evolution of AMP and its validity is proven for closely related models in [START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF]. We denote the corresponding fixed point m AMP and the corresponding estimation error MSE AMP = 1 -m AMP .

The phase diagram presented in Fig. 1 summarizes this theory for the spiked 2 + 3-spin model. It is deduced by investigating the local maxima of the scalar function [START_REF] Seung | Statistical mechanics of learning from examples[END_REF]. Notably we observe that the phase diagram in terms of ∆ 2 and ∆ p splits into three phases • Easy in green for ∆ 2 < 1 and any ∆ p : The fixed point of the state evolution [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF] is the global maximizer of the free entropy [START_REF] Seung | Statistical mechanics of learning from examples[END_REF], and m * = m AMP > 0.

• Hard in orange for ∆ 2 > 1 and low ∆ p < ∆ IT p (∆ 2 ): The fixed point of the state evolution ( 6) is not the global maximizer of the free entropy [START_REF] Seung | Statistical mechanics of learning from examples[END_REF], and m * > m AMP = 0.

• Impossible in red for ∆ 2 > 1 and high ∆ p > ∆ IT p (∆ 2 ): The fixed point of the state evolution ( 6) is the global maximizer of the free entropy [START_REF] Seung | Statistical mechanics of learning from examples[END_REF], and m * = m AMP = 0.

For the 2 + p-spin model with p > 3 the phase diagram is slightly richer and is presented in the appendix.

Langevin Algorithm and its Analysis

We now turn to the core of the paper and the analysis of the Langevin algorithm. In statistics, the most commonly used way to compute the Bayes-optimal estimator (3) is to attempt to sample the posterior distribution (2) and use several independent samples to compute the expectation in [START_REF] Lecun | Deep learning[END_REF]. In order to do that one needs to set up a stochastic dynamics on x that has a stationary measure at long times given by the posterior measure [START_REF] Welling | Bayesian learning via stochastic gradient Langevin dynamics[END_REF]. The Langevin algorithm is one of the possibilities (others include notably Monte Carlo Markov chain). The common bottleneck is that the time needed to achieve stationarity can be in general exponential in the system size. In which case the algorithm is practically useless. However, this is not always the case and there are regions in parameter space where one can expect that the relaxation to the posterior measure happens on finite timescales. Therefore it is crucial to understand where this happens and what are the associated relaxation timescales.

The Langevin algorithm on the hypersphere with Hamiltonian given by Eq. ( 4) reads

ẋi (t) = -µ(t)x i (t) - ∂H ∂x i + η i (t) , (7)
where η i (t) is a zero mean noise term, with η i (t)η j (t ) = 2δ ij δ(tt ) where the average • is with respect to the realizations of the noise. The Lagrange multiplier µ(t) is chosen in such a way that the dynamics remains on the hypersphere. In the large N -limit one finds µ(t) = 1 -2H 2 (t) -pH p (t) where the H 2 (t) is the 1st term from (4) evaluated at x(t), and H p (t) is the value of the 2nd term from (4).

The presented spiked matrix-tensor model falls into the particular class of spherical 2+p-spin glasses [START_REF] Crisanti | Spherical 2+ p spin-glass model: An exactly solvable model for glass to spin-glass transition[END_REF][START_REF] Crisanti | Spherical 2+ p spin-glass model: An analytically solvable model with a glass-to-glass transition[END_REF] for which the performance of the Langevin algorithm can be tracked exactly in the large-N limit via a set of integro-partial differential equations [START_REF] Crisanti | The spherical p-spin interaction spin-glass model[END_REF][START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF], beforehand dubbed CHSCK. We call this generalised version of the CHSCK equations Langevin State Evolution (LSE) equations in analogy with the state evolution of AMP.

In order to write the LSE equations, we defined three dynamical correlation functions

C N (t, t ) ≡ 1 N N i=1 x i (t)x i (t ) , (8)
C N (t) ≡ 1 N N i=1 x i (t)x * i , (9)
R N (t, t ) ≡ 1 N N i=1 ∂x i (t)/∂h i (t )| h i =0 , (10)
where h i is a pointwise external field applied at time t to the Hamiltonian as H + i h i x i . We note that the correlation functions defined above depend on the realization of the thermal history (i.e. of the noise η(t)) and on the disorder (here the matrix Y and tensor T ). However, in the large-N limit they all concentrate around their averages. We thus define C(t, t ) = lim N →∞ E Y,T C N (t, t ) η and analogously for C(t) and R(t, t ). Standard field theoretical methods [START_REF] Cecil Martin | Statistical dynamics of classical systems[END_REF] or dynamical cavity method arguments [START_REF] Mézard | Spin glass theory and beyond[END_REF] can then be used to obtain a closed set of integro-differential equations for the averaged dynamical correlation functions, describing the average global evolution of the system under the Langevin algorithm. The resulting LSE equations are (see the appendix Sec. C for a complete derivation)

∂ ∂t C(t, t ) = 2R(t , t) -µ(t)C(t, t ) + Q (C(t))C(t ) + t 0 dt R(t, t )Q (C(t, t ))C(t , t )+ t 0 dt R(t , t )Q (C(t, t )) , ∂ ∂t R(t, t ) = δ(t -t ) -µ(t)R(t, t ) + t t dt R(t, t )Q (C(t, t ))R(t , t ) , ∂ ∂t C(t) = -µ(t)C(t) + Q (C(t)) + t 0 dt R(t, t )C(t )Q(C(t, t )) , (11)
where we have defined Q(x) = x 2 /(2∆ 2 ) + x p /(p∆ p ). The Lagrange multiplier, µ(t), is fixed by the spherical constraint, through the condition C(t, t) = 1 ∀t. Furthermore causality implies that R(t, t ) = 0 if t < t . Finally the Ito convention on the stochastic equation [START_REF] Deshpande | Finding hidden cliques of size (N/e) in nearly linear time[END_REF] gives ∀t lim t →t -R(t, t ) = 1.

Behavior of the Langevin algorithm

In order to assess the perfomances of the Langevin algorithm and compare it with AMP, we notice that the correlation function C(t) is directly related to accuracy of the algorithm. We solve the differential equations [START_REF] Aubin | The committee machine: Computational to statistical gaps in learning a two-layers neural network[END_REF] numerically along the lines of [START_REF] Kim | The dynamics of the spherical p-spin model: From microscopic to asymptotic[END_REF][START_REF] Berthier | Spontaneous and induced dynamic fluctuations in glass formers. I. General results and dependence on ensemble and dynamics[END_REF], for a detailed procedure see the appendix Sec. D, codes available at [START_REF] Sarao Mannelli | Langevin state evolution integrators[END_REF]. In Fig. 2 we plot the correlation with the spike C(t) as a function of the running time t for p = 3, fixed ∆ 2 = 0.7 and several values of ∆ p . In the inset of the plot we compare it to the same quantity obtained from the state evolution of the AMP algorithm. For the Langevin algorithm in Fig. 2 we see a pattern that is striking. One would expect that as the noise ∆ p decreases the inference problem is getting easier, the correlation with the signal is larger and is reached sooner in the iteration. This is, after all, exactly what we observe for the AMP algorithm in the inset of Fig. 2. Also for the Langevin algorithm the plateau reached for large times t becomes higher (better accuracy) as the noise ∆ p is reduced. Furthermore the height of the plateau coincides with that reached by AMP, thus testifying the algorithm reached equilibrium. However, contrary to AMP, the relaxation time for the Langevin algorithm increases dramatically when diminishing ∆ p (notice the log scale on x-axes of Fig. 2, as compared to the linear scale of the inset). We will define τ as the time it takes for the correlation to reach a value C plateau /2. We then plot the value of this equilibration time in the insets of Fig. 3 as a function of the noise ∆ p or ∆ 2 having fixed ∆ 2 (upper panel) or ∆ p (lower panel) respectively. The data are clearly consistent with a divergence of τ at a certain finite value of ∆ * p and ∆ * 2 . We fit the data by a power law fit τ

(∆) = 1 ∆ -1 ∆ * -γ
and obtain in the particular case of fixed ∆ 2 = 0.7 a fit with γ = 1.75 and ∆ * p = 0.58, whereas for fixed ∆ p = 1.0 we obtain γ = 1.42 and ∆ * 2 = 0.86. However, we are not able to strictly prove that the divergence of the relaxation time truly occurs, but at least our results imply that for ∆ p < ∆ * p and ∆ 2 > ∆ * 2 the Langevin algorithm [START_REF] Deshpande | Finding hidden cliques of size (N/e) in nearly linear time[END_REF] is not a practical solver for the spiked mixed matrix-tensor problem. We will call the region ∆ p < ∆ * p and ∆ 2 < 1 where the AMP algorithm works optimally without problems yet Langevin algorithm does not, the Langevin-hard region. Both ∆ * p and ∆ * 2 are then plotted in Fig. 1 with green points (circles and stars) and consistently delimit the Langevin-hard region that extends considerably into the region where the AMP algorithm works optimally in a small number or iterations. Our main conclusion is thus that the Langevin algorithm designed to sample the posterior measure works efficiently in a considerably smaller region of parameters than the AMP as quantified in Fig. 1.

Fig. 4 presents another way to depict the observed data, the correlation C(t) reached after time t is plotted as a function of the tensor noise variance ∆ p . The results of AMP are depicted with dotted lines and, as one would expect, decrease monotonically as the noise ∆ p increases. The equilibrium value (black dashed) is reached within few dozens of iterations. On the contrary, the correlation reached by the Langevin algorithm after time t is non-monotonic and close to zero for small values of noise ∆ p signalling again a diverging relaxation time when ∆ p is decreased. : Fit using a power law of the relaxation times of the Langevin algorithm, the point of divergence marks the limit of the Langevin easy region. These fits have been performed both, for fixed ∆ 2 = 0.7 (blue circles and yellow crosses) and for fixed ∆ p = 1 (pink circles). The circles are obtained with numerical solution of LSE that uses the dynamical grid while crosses are obtained using a fixed-grid (details in Sec. D.1).

Glassy nature of the Langevin-hard phase

We notice that for ∆ p → 0 and any finite ∆ 2 , and after a suitable rescaling of time and temperature, the LSE equations coincide with the ones of the pure spiked-tensor (p = 3) at very low temperature and slightly perturbed by additional terms. The glassy nature of the landscape of such model [START_REF] Ros | Complex energy landscapes in spiked-tensor and simple glassy models: ruggedness, arrangements of local minima and phase transitions[END_REF][START_REF] Ben Arous | Algorithmic thresholds for tensor PCA[END_REF] qualitatively justifies why the Langevin algorithm remains trapped in one of the spurious minima, hence leading to a hard-Langevin phase for ∆ p → 0 and any finite ∆ 2 (in particular ∆ 2 < 1).

In order to obtain a quantitative estimate of the extent of the glassy region, we repeated the analysis of [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF] in the present model (details in the appendix Sec. E). In particular, we compute the entropy of metastable states (finite temperature extensions of spurious minima), a.k.a. the complexity, using the one-step replica symmetry breaking (1RSB) assumption for the structure of these states. We then analyze the stability of the 1RSB solution with respect to further levels of RSB and draw the region of paramaters ∆ 2 , ∆ p for which states of positive complexity exist and are stable. This happens to be below the blue dashed-dotted line depicted in Fig. 1.

We observe that indeed in the regime where the 1RSB solution indicates the existence of an exponential number of spurious metastable states, the Langevin algorithm does not work. Moreover, the curve delimiting existence of the 1RSB stable states has the same trend as the boundary of the Langevin-hard regime. Yet the Langevin-hard phase extends to larger values of ∆ p , see Fig. 1. This quantitative discrepancy can be due to either the fact that we did not take into account effects of full-step replica symmetry breaking [START_REF] Rizzo | Replica-symmetry-breaking transitions and off-equilibrium dynamics[END_REF] or rather because the relation between the static replica calculation and the behavior of the Langevin algorithm is more complex than anticipated in the literature [START_REF] Rizzo | Replica-symmetry-breaking transitions and off-equilibrium dynamics[END_REF]. We let this point for investigation in future work.

Better performance by annealing the landscape

A generic strategy to avoid metastable states is by simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] where the thermal noise is given by η i (t)η j (t ) = 2δ ij T (t)δ(tt ) and the temperature T (t) is time dependent and slowly decreased during the dynamics starting from very large values towards the target temperature, T = 1 in our case. We tried this algorithm and found that it does not succeed to recover the signal.

As discussed above, it is the tensor part of the Hamiltonian that induces glassiness. Therefore, we consider a protocol in which contribution of the tensor is weak at initial stages of the dynamics and increases gradually into its Bayes-optimal strength. Specifically, we focus on a time-dependent Hamiltonian, H 2 (Y ) + 1 Tp(t) H p (T ), and change dynamically T p (t) according to T p (t) = 1 + C ∆p e -t/τann , C being a constant. The target value of T p (t → ∞) = 1 is the Bayes-optimal one and τ ann is the characteristic timescale of the annealing procedure. In Fig. 5 we show that if τ ann is sufficiently large, the Langevin hard phase is destroyed and this procedure approaches the performances as AMP. More specifically, before the H p (T ) component becomes relevant, i.e. for times that are τ ann log C ∆p , the modified Langevin algorithm tends to the equilibrium solution of the pure spiked matrix model. The optimal AMP reconstruction is subsequently allowed when the tensor component gains its full weight, which occurs later for slower annealing. Conversely if τ ann is too small, the Langevin algorithm does not have the time for escaping the glassy region of the tensor component before its own contribution to the dynamics becomes relevant and irremediably hampers the final reconstruction of the signal.

It is interesting to underline that the above finding is somewhat paradoxical from the point of view of Bayesian inference. In the present setting we know perfectly the model that generated the data and all its parameters, yet we see that for the Langevin algorithm it is computationally advantageous to mismatch the parameter T p in order to reach faster convergence to equilibrium. This is particularly striking given that for AMP it has been proven in [START_REF] Deshpande | Finding hidden cliques of size (N/e) in nearly linear time[END_REF] that mismatching the parameters can never improve the performance.

Perspectives

In this work we have investigated the performances of the Langevin algorithm considered as a tool to sample the posterior measure in the spiked matrix-tensor model. We have shown that the Langevin algorithm fails to find the signal in part of the AMP-easy region. Our analysis is based on the Langevin State Evolution equations that describe the evolution of the algorithm in the large size limit.

In this work we managed to find the landscape-annealing protocol under which the Langevin algorithm is able to match the performance of AMP by relying on knowledge about generative model. It would be an interesting direction for future work to investigate whether the performance of the Langevin algorithm can be improved with some model-blind manner.

While we studied here the spiked matrix-tensor model, we expect that our findings are universal because they are due to the glassiness of the hard phase and therefore should apply to any local sampling dynamics, e.g. to Monte Carlo Markov chains. An interesting extension of this work would investigate algorithms closer to stochastic gradient descent and models closer to current neural network architectures. Figure 5: Behaviour of the correlation with the signal in Langevin algorithm where the tensor-related temperature T p is annealed as T p (t) = 1 + C ∆p e -t/τann with C = 100 in the Langevin-hard regime with ∆ 2 = 0.70, ∆ p = 0.10. We show the behavior for several rates τ (solid lines) and compare to the quenched Langevin algorithm (dashed line close to zero), and to the value reached by AMP for the full model (upper dotted line) and for the pure matrix model (lower dotted line). We see that unless the annealing is very fast, it reaches the AMP value. When the annealing is very slow it takes time to reach the AMP value.

Appendix A Definition of the spiked matrix-tensor model

We consider a teacher-student setting in which the teacher constructs a matrix and a tensor from a randomly sampled signal and the student is asked to recover the signal from the observation of the matrix and tensor provided by the teacher [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF].

The signal, x * is an N -dimensional vector whose entries are real i.i.d. random variables sampled from the normal distribution (i.e. the prior is P X ∼ N (0, 1)). The teacher generates from the signal a symmetric matrix and a symmetric tensor of order p. Those two objects are then transmitted through two noisy channels with variances ∆ 2 and ∆ p , so that at the end one has two noisy observations given by

Y ij = x * i x * j √ N + ξ ij , (12)
T i 1 ,...,ip = (p -1)! N (p-1)/2 x * i 1 . . . x * ip + ξ i 1 ,...,ip , (13)
where, for i < j and i 1 < • • • < i p , ξ ij and ξ i 1 ,...,ip are i.i.d. random variables distributed according to ξ ij ∼ N (0, ∆ 2 ) and ξ i 1 ,...,ip ∼ N (0, ∆ p ). The ξ ij and ξ i 1 ,...,ip are symmetric random matrix and tensor, respectively. Given Y ij and T i 1 ,...,ip the inference task is to reconstruct the signal x * . In order to solve this problem we consider the Bayesian approach. This starts from the assumption that both the matrix and tensor have been produced from a process of the same kind of the one described by Eq. [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]. Furthermore we assume to know the statistical properties of the channel, namely the two variances ∆ 2 and ∆ p , and the prior on x. Given this, the posterior probability distribution over the signal is obtained through the Bayes formula

P (X|Y, T ) = P (Y, T |X)P (X) P (Y, T ) , (14)
where

P (Y, T |X) = i<j P Y Y ij x i x j √ N i 1 <•••<ip P T T i 1 ...ip (p -1)! N (p-1)/2 x i 1 . . . x ip = ∝ i<j e -1 2∆ 2 Y ij - x i x j √ N 2 i 1 <•••<ip e -1 2∆p T i 1 ...ip - √ (p-1)! N (p-1)/2 x i 1 ...x ip 2 . (15)
Therefore we have

P (X|Y, T ) = 1 Z(Y, T ) i e -1 2 x 2 i i<j e -1 2∆ 2 Y ij - x i x j √ N 2 i 1 <•••<ip e -1 2∆p T i 1 ...ip - √ (p-1)! N (p-1)/2 x i 1 ...x ip 2 , (16)
where Z(Y, T ) is a normalization constant. Plugging Eqs. [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] into Eq. ( 16) allows to rewrite the posterior measure in the form of a Boltzmann distribution of the mixed 2 + p-spin Hamiltonian [START_REF] Crisanti | Spherical 2+ p spin-glass model: An exactly solvable model for glass to spin-glass transition[END_REF][START_REF] Crisanti | Spherical 2+ p spin-glass model: An analytically solvable model with a glass-to-glass transition[END_REF][START_REF] Crisanti | Exactly solvable spin-glass models with ferromagnetic couplings: The spherical multi-p-spin model in a self-induced field[END_REF]]

H = - 1 ∆ 2 √ N i<j ξ ij x i x j - (p -1)! ∆ p N p-1 2 i 1 <•••<ip ξ i 1 ...ip x i 1 . . . x ip - N 2∆ 2   1 N i x i x * i   2 + - N p∆ p   1 N i x i x * i   p - 1 2 N i=1 x 2 i + const. ( 17
)
so that

P (X|Y, T ) = 1 Z(Y, T ) e -H . (18)
In the following we will refer to Z(Y, T ) as the partition function. We note here that in the large N limit, using a Gaussian prior on the variables x i is equivalent to consider a flat measure over the N -dimensional hypersphere N i=1 x 2 i = N . This choice will be used when we will describe the Langevin algorithm and in this case the last term in the Hamiltonian will become an irrelevant constant.

B Approximate Message Passing, state evolution and phase diagrams

Approximate Message Passing (AMP) is a powerful iterative algorithms to compute the local magnetizations x i given the observed matrix and tensor. It is rooted in the cavity method of statistical physics of disordered systems [START_REF] David J Thouless | Solution of'solvable model of a spin glass[END_REF][START_REF] Mézard | Spin glass theory and beyond[END_REF] and it has been recently developed in the context of statistical inference [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF], where in the Bayes optimal case it has been conjectured to be optimal among all local iterative algorithms. Among the properties that make AMP extremely useful is the fact that its performances can be analyzed in the thermodynamic limit. Indeed in such limit, its dynamical evolution is described by the so called State Evolution (SE) equations [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF]. In this section we derive the AMP equations and their SE description for the spiked matrix-tensor model and solve them to obtain the phase diagram of the model as a function of the variances ∆ 2 and ∆ p of the two noisy channels.

B.1 Approximate Message Passing and Bethe free entropy

AMP can be obtained as a relaxed Gaussian closure of the Belief Propagation (BP) algorithm. The derivation that we present follows the same lines of [START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF][START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF]. The posterior probability can be represented as a factor graph where all the variables are represented by circles and are linked to squares representing the interactions [START_REF] Mézard | Information, physics, and computation[END_REF].

x 1 x 2 x N P X P X P X P Y P Y P Y P T . . . m12 → 1 (x 2 ) m 2 → 1 2 ( x 2 ) t 1 2 N → 1 ( x 2 ) t 2 → 1 2 N ( x 2 )
Figure 6: The factor graph representation of the posterior measure of the matrix-tensor factorization model. The variable nodes represented with white circles are the components of the signal while black squares are factor nodes that denote interactions between the variable nodes that appear in the interaction terms of the Boltzmann distribution in Eqs. [START_REF] Richard | A statistical model for tensor PCA[END_REF][START_REF] Samuel B Hopkins | Tensor principal component analysis via sum-ofsquare proofs[END_REF]. There are three types of factor nodes: P X is the prior that depends on a single variable, P Y that is the probability of observing a matrix element Y ij given the values of the variables x i and x j , and finally P T that is the probability of observing a tensor element T i 1 ,...,ip . The posterior, apart from the normalization factor, is simply given by the product of all the factor nodes.

This representation is very convenient to write down the BP equations. In the BP algorithm we iteratively update until convergence a set of variables, which are beliefs of the (cavity) magnetization of the nodes. The intuitive underlying reasoning behind how BP works is the following. Given the current state of the variable nodes, take a factor node and exclude one node among its neighbors. The remaining neighbors through the factor node express a belief on the state of the excluded node. This belief is mathematically described by a probability distribution called message, mt ij→i (x i ) and tt ii 2 ...ip→i (x i ) depending on which factor node is selected. At the same time, another belief on the state of the excluded node is given by the rest of the network but the factor node previously taken into account, m i→ij (x i ) and t i→ii 2 ...ip (x i ) respectively. All these messages travel in the factor graph carrying partial information on the real magnetization of the single nodes, and they are iterated until convergence1 . The iterative scheme is described by the following equations

mt ij→i (x i ) ∝ dx j m t j→ij (x j )P Y Y ij x i x j √ N , (19)
m t+1 i→ij (x i ) ∝ P X (x i ) l =j mt il→i (x i ) i 2 <•••<ip tt ii 2 ...ip→i (x i ), (20)
tt ii 2 ...ip→i (x i ) ∝ l=2...p dx l t t i l →ii 2 ...ip (x l ) P T T ii 2 ...ip (p -1)! N (p-1)/2 x i x i 2 . . . x ip , (21)
t t+1 i→ii 2 ...ip (x i ) ∝ P X (x i ) l mt il→i (x i ) k 2 <•••<kp =i 2 ...ip tt ik 2 ...kp→i (x i ) (22)
and we have omitted the normalization constants that guarantee that the messages are probability distributions.

When the messages have converged to a fixed point, the estimation of the local magnetizations can be obtained through the computation of the real marginal probability distribution of the variables given by

µ i (x i ) =   j( =i) dx j   P (X|Y, T ) = P X (x i ) l mt il→i (x i ) i 2 <•••<ip tt ii 2 ...ip→i (x i ) . (23)
We note that the computational cost to produce an iteration of BP scales as O(N p ). Furthermore Eqs. (19 -22) are iterative equations for continuous functions and therefore are extremely hard to solve when dealing with continuous variables. The advantage of AMP is to reduce drastically the computational complexity of the algorithm by closing the equations on a Gaussian ansatz for the messages. This is justified in the present context since the factor graph is fully connected and therefore each iteration step of the algorithm involves sums of a large number of independent random variables that give rise to Gaussian distributions. Gaussian random variables are characterized by their mean and covariance that are readily obtained for N 1 expanding the factor nodes for small

ω ij = x i x j / √ N and ω i 1 ...ip = (p -1)!x 1 . . . x p /N p-1
2 . Once the BP equations are relaxed on Gaussian messages, the final step to obtain the AMP algorithm is the so-called TAPyfication procedure [START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF][START_REF] David J Thouless | Solution of'solvable model of a spin glass[END_REF], which exploits the fact that the procedure of removing one node or one factor produces only a weak perturbation to the real marginals and therefore can be described in terms of the real marginals of the variable nodes themselves. By applying this scheme we obtain the AMP equations, which are described by a set of auxiliary variables A (k) and B (k) i and by the mean x i and variance σ i = x 2 i of the marginals of variable nodes. The AMP iterative equations are

B (2),t i = 1 ∆ 2 √ N k Y ki xt k - 1 ∆ 2   1 N k σ t k   xt-1 i ; (24)
A (2),t = 1 ∆ 2 N k xt k 2 ; (25)
B (p),t i = (p -1)! ∆ p N (p-1)/2 k 2 ...kp T ik 2 ...kp xt k 2 . . . xt kp - p -1 ∆ p      1 N k σ t k     1 N k xt k 2   p-2    xt-1 i ; (26)
A (p),t = 1 ∆ p   1 N k xt k 2   p-1 ; ( 27
) xt+1 i = f (A (2) + A (p) , B (2)
i + B (p) i ) ; (28)
σ t+1 i = ∂ ∂B f (A, B) A=A (2) +A (p) ,B=B (2)
i +B

(p) i , (29)
f (A, B) ≡ dx 1 Z(A, B) xP X (x)e Bx-1 2 Ax 2 = B 1 + A . ( 30
)
It can be shown that these equations can be obtained as saddle point equations from the so called Bethe free entropy defined as Φ Bethe = log Z Bethe (Y, T )/N where Z Bethe is the Bethe approximation to the partition function which is defined as the normalization of the posterior measure. The expression of the Bethe free entropy per variable can be computed in a standard way (see [START_REF] Mézard | Information, physics, and computation[END_REF]) and it is given by

Φ Bethe = 1 N    i log Z i + i≤j log Z ij + i 1 ≤•••≤ip log Z i 1 ...ip - i(ij) log Z i,ij - i(ii 2 ...ip) log Z i(ii 2 ...ip)    , (31)
where

Z i = dx i P X (x i ) j mij→i (x i ) (i 2 ...ip) tii 2 ...ip→i (x i ), Z ij = j( =i) dx j m j→ij (x j ) i<j e -1 2∆ 2 Y ij - x i x j √ N 2 , Z i 1 ...ip = p l=1 dx i l t i l →i 1 ...ip (x i l ) i 1 <•••<ip e -1 2∆p T i 1 ...ip - √ (p-1)! N (p-1)/2 x i 1 ...x ip 2 , Z i(ij) = dx i m i→ij (x) mij→i (x i ), Z i(ii 2 ...ip) = dx i t i→ii 2 ...ip (x) tii 2 ...ip→i (x i )
are a set of normalization factors. Using the Gaussian approximation for the messages and employing the same TAPyification procedure used to get the AMP equations we obtain the Bethe free entropy density as

Φ Bethe = 1 N i log Z(A (p) + A (2) , B (p) i + B (2) i ) + p -1 p 1 N i -B (p) i xi + A (p) i x2 i + σ i 2 + + p -1 2p∆ p i x2 i N p-1 i σ i N + 1 2N i -B (2)
i xi + A (2) i x2 i + σ i 2 + 1 4∆ 2 i x2 i N i σ i N , (32)
where we used the variables defined in eqs. [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF][START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF][START_REF] Crisanti | The spherical p-spin interaction spin-glass model[END_REF][START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF] for sake of compactness and Z(A, B) is defined as

Z(A, B) = dxP X (x)e Bx-Ax 2 2 = 1 √ A + 1 e B 2 2(A+1) . (33)

B.2 Averaged free entropy and its proof

Eq. ( 32) represents the Bethe free entropy for a single realization of the factor nodes in the large size limit.

Here we wish to discuss the actual, exact, value of this free entropy, that is:

f N (Y, T ) = log Z(Y, T ) N
This is a random variable, since it depends a priori on the planted signal and the noise in the tensor and matrices. However one expects that, since free entropy is an intensive quantity, we expect from the statistical physics intuition that it should be self averaging and concentrate around its mean value in the large N limit [START_REF] Mézard | Spin glass theory and beyond[END_REF]. In fact, this is easily proven. First, since the spherical model has a rotational symmetry, one may assume the planted assignment could be any vector on the hyper-sphere, and we might as well suppose it is the uniform one x * i = 1∀ i: the true source of fluctuation comes from the noise Y and T . These can be controlled by noticing that the free entropy is a Lipshitz function of the Gaussian random variable Y and T . Indeed:

∂ Y ij f N (Y, T ) = 1 ∆ 2 N √ N x i x j
so that the free energy f N is Lipschitz with respect to Y with constant

L = 1 ∆ 2 N √ N i<j x i x j 2 ≤ 1 ∆ 2 N √ N 1 2 i,j x i x j 2 = 1 ∆ 2 N √ N 1 2 i,j x i xi x j xj
where x represent a copy (or replica) of the system. In this case

L ≤ 1 ∆ 2 N √ N N 2 i x i xi N 2 = q 2 ∆ 2 √ N
where q is the overlap between the two replica x and x, that is bounded by one on the sphere, so

L ≤ 1 ∆ 2 √
N . Therefore, by Gaussian concentration of Lipschitz functions (the Tsirelson-Ibragimov-Sudakov inequality [START_REF] Boucheron | Concentration inequalities[END_REF]), we have for some constant K:

Pr |f n -E Y f n | ≥ t ≤ 2e -N t 2 /K (34)
and it particular any fluctuation larger than O(1/ √ N ) is (exponentially) rare. A similar computaton shows that f N also concentrates with respect to the tensor T . This shows that in the large size limit, we can consider the averaged free entropy:

F N ≡ 1 N E [log Z N ]
With our (non-rigorous) statistical physics tools, this can be obtained by averaging Eq. ( 32) over the disorder, see for instance [START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF], and this yields an expression for the free energy called the replica symmetric (RS) formula:

Φ RS = lim N →∞ E Y,T log Z(Y, T ) N . (35)
We now state precisely the form of Φ RS and prove the validity of Eq. ( 35). The RS free entropy for any prior distribution P X reads as

Φ RS ≡ max m ΦRS (m) where ΦRS (m) = E W,x *     log   Z   m ∆ 2 + m p-1 ∆ p , m ∆ 2 + m p-1 ∆ p x * + m ∆ 2 + m p-1 ∆ p W          - 1 4∆ 2 m 2 - p -1 2p∆ p m p , ( 36
)
where W is a Gaussian random variable of zero mean and unit variance and x * is a random variables taken from the prior P X . We remind that the function Z(A, B) is defined via Eq. [START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF].

For Gaussian prior P X , which is the one of interest here, we obtain

ΦRS (m) = - 1 2 log m ∆ 2 + m p-1 ∆ p + 1 + 1 2 m ∆ 2 + m p-1 ∆ p - 1 4∆ 2 m 2 - p -1 2p∆ p m p . ( 37
)
The expression given in the main text is slightly different but can be obtained as follow. First notice that the extremization condition for ΦRS (m) reads

m = 1 - 1 1 + m ∆ 2 + m p-1 ∆p (38)
and by plugging this expression in Eq. ( 37) we recover the more compact expression Φ RS (m) showed in the main text:

Φ RS (m) = 1 2 log (1 -m) + m 2 + m 2 4∆ 2 + m p 2p∆ p . ( 39
)
The two expressions Φ RS (m) and ΦRS (m) are thus equal for each value of m that satisfy Eq. ( 38). The parameter m can be interpreted as the average correlation between the true and the estimated signal

m = 1 N N i=1 x * i xi . ( 40
)
The average minimal mean squared error (MMSE) can be obtained from the maximizer m of the average Bethe free entropy as

MMSE ≡ 1 N N i=1 (x * i -xi ) 2 = 1 -m * , where m * = argmax ΦRS (m) . ( 41
)
where the overbar stands for the average over the signal x * and the noise of the two Gaussian channels. The validity of Eq. ( 36) can be proven rigorously for every prior having a bounded second moment. The proof we shall present is a straightforward generalization of the one presented in [START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF] for the pure tensor case, and in [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF] for the matrix case, and it is based on two main ingredients. The first one is the Guerra interpolation method applied on the Nishimori line [START_REF] Babu | Exact solution of the gauge symmetric p-spin glass model on a complete graph[END_REF][START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF][START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF], in which we construct an interpolating Hamiltonian that depends on a parameter t ∈ [0; 1] that is used to move from the original Hamiltonian of Eq. ( 17), to the one corresponding to a scalar denoising problem whose free entropy is given by the first term in Eq. [START_REF] Cecil Martin | Statistical dynamics of classical systems[END_REF]. The second ingredient is the Aizenman-Sims-Starr method [START_REF] Michael Aizenman | Extended variational principle for the sherringtonkirkpatrick spin-glass model[END_REF] which is the mathematical version of the cavity method (note that other techniques could also be employed to obtain the same results, see [START_REF] Barbier | Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula[END_REF][START_REF] Barbier | The adaptive interpolation method: A simple scheme to prove replica formulas in bayesian inference[END_REF][START_REF] Krzakala | Estimation in the spiked wigner model: A short proof of the replica formula[END_REF][START_REF] Mourrat | Hamilton-Jacobi equations for mean-field disordered systems[END_REF]). The theorem we want to prove is: Theorem 1 (Replica-Symmetric formula for the free energy). Let P X be a probability distribution over R, with finite second moment Σ X . Then, for all ∆ 2 > 0 and ∆ p > 0

F N ≡ 1 N E [log Z N ] ----→ N →∞ sup m≥0 ΦRS (m) ≡ Φ RS (∆ 2 , ∆ p ) . (42)
For almost every ∆ 2 > 0 and ∆ p > 0, ΦRS admits a unique maximizer m over R + × R + and

T -MMSE N ----→ N →∞ Σ p X -(m * ) p , M -MMSE N ----→ N →∞ Σ 2 X -(m * ) 2 .
Here, we have defined the tensor-MMSE T-MMSE N by the error in reconstructing the tensor:

T-MMSE N (∆ 2 , ∆ p ) = inf θ    p! N p i 1 <•••<ip x 0 i 1 . . . x 0 ip -θ(Y ) i 1 ...ip 2    ,
and the matrix-MMSE M-MMSE N by the error in reconstructing the matrix:

M-MMSE N (∆ 2 , ∆ p ) = inf θ    2 N 2 i<j x 0 i x 0 j -θ(Y ) i,j 2    ,
where in both cases the infimum is taken over all measurable functions θ of the observations Y . The result concerning the MMSE is a simple application of the I-MMSE theorem [START_REF] Dongning Guo | Mutual information and minimum mean-square error in gaussian channels[END_REF], that relates the derivative of the free energy with respect to the noise variances and the MMSE. The details of the arguments are the same than in the matrix (p = 2) case ( [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF], corollary 17) and the tensor one ([10], theorem 2). Indeed, as discussed in [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF], these M-MMSE and T-MMSE results implies the vector MMSE result of Eq. ( 41) when p is odd, and thus in particular for the p = 3 case discussed in the main text.

Sketch of proof

In this section we give a detailed sketch of the proof theorem 1. Following the techniques used in many recent works [START_REF] Babu | Exact solution of the gauge symmetric p-spin glass model on a complete graph[END_REF][START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF][START_REF] Barbier | Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula[END_REF][START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF][START_REF] Barbier | The adaptive interpolation method: A simple scheme to prove replica formulas in bayesian inference[END_REF][START_REF] Krzakala | Estimation in the spiked wigner model: A short proof of the replica formula[END_REF][START_REF] Barbier | Phase transitions, optimal errors and optimality of message-passing in generalized linear models[END_REF], we shall make few technical remarks:

• We will consider only priors with bounded support, supp(P X ) = S ⊂ [-K; K]. This allows to switch integrals and derivatives without worries. This condition can then be relaxed to unbounded distributions with bounded second moment using the same techniques as the ones that we are going to present, and the proof is therefore valid in this case. This is detailed for instance in [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF] sec. 6.2.2.

• Another key ingredient is the introduction of a small perturbation in the model that takes the form of a small amount of side information. This kind of techniques are frequently used in statistical physics, where a small "magnetic field" forces the Gibbs measure to be in a single pure state [START_REF] Georgii | Gibbs measures and phase transitions[END_REF]. It has also been used in the context of coding theory [START_REF] Macris | Griffith-Kelly-Sherman correlation inequalities: A useful tool in the theory of error correcting codes[END_REF] for the same reason. In the context of Bayesian inference, we follow the generic scheme proposed by Montanari in [START_REF] Montanari | Estimating random variables from random sparse observations[END_REF] (see also [START_REF] Coja-Oghlan | Information-theoretic thresholds from the cavity method[END_REF]) and add a small additional source of information that allows the system to be in a single pure state so that the overlap concentrates on a single value. This source depends on Bernoulli random variables

L i i.i.d. ∼ Bern( ), i ∈ [N ]; if L i = 1
, the channel, call it A, transmits the correct information. We can then consider the posterior of this new problem, P (X|A, Y, T ), and focus on the associated free energy density F N, defined as the expected value of the average of the logarithm of normalization constant divided by the number of spins. Then we can immediately prove that for all N ≥ 1 and , ∈ [0; 1] it follows:

F N, -F N, ≤ K 2p ∆p + K 4 ∆ 2 | -|.
This allows (see for instance [START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF]) to obtain the concentration of the posterior distribution around the replica parameter (q = 1 N x (1) • x (2) )

E x (1) • x (2) N -q 2 N →∞ -→ 0 ; (43)
E x * • x N -q 2 N →∞ -→ 0 , (44)
where x, x (1) , x (2) are sampled from the posterior distribution and the averages • and E[•] are respectively the average over the posterior measure and the remaining random variables.

• Finally, a fundamental property of inference problems which is a direct consequence of the Bayes theorem2 , is the so-called Nishimori symmetry [START_REF] Nishimori | Statistical physics of spin glasses and information processing: an introduction[END_REF][START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF]: Let (X, Y ) be a couple of random variables on a polish space. Let k ≥ 1 and let X (1) , . . . , X (k) be k i. 1) , . . . , X (k-1) , X) .

E f (Y, X (1) , . . . , X (k) ) = E f (Y, X (
While the consequences of this identity are important, the proof is rather simple: It is equivalent to sample the couple (X, Y ) according to its joint distribution or to sample first Y according to its marginal distribution and then to sample X conditionally to Y from its conditional distribution P (X = • | Y ). Thus the (k + 1)-tuple (Y, X (1) , . . . , X (k) ) is equal in law to (Y, X (1) , . . . , X (k-1) , X).

The proof of Theorem 1 is obtained by using the Guerra interpolation technique to prove a lower bound for the free entropy and then by applying the Aizenman-Sims-Star scheme to get a matching upper bound.

Lower bound: Guerra interpolation We now move to the core of the proof. The first part combines the Guerra interpolation method [START_REF] Guerra | The thermodynamic limit in mean field spin glass models[END_REF] developed for matrices in [START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF] and tensors in [START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF].

Consider the interpolating Hamiltonian depending of t ∈ [0, 1]

H N,t = - i<j √ t ∆ 2 √ N Y ij x i x j + t 2∆ 2 N (x i x j ) 2 + - i 1 <•••<ip   t(p -1)! ∆ p N p-1 2
T i 1 ...ip x i 1 . . . x ip + t(p -1)! 2∆ p N p-1 (x i 1 . . . x ip ) 2   + - j   √ 1 -t m p-1 ∆ p + m ∆ 2 W j x j + (1 -t) m p-1 ∆ p + m ∆ 2 x * j x j + 1 -t 2 m p-1 ∆ p + m ∆ 2 x 2 j   , (45)
where we have for t = 1 the regular Hamiltonian and for t = 0 the first term of Eq. ( 36) where W j are i.i.d. canonical Gaussian variables. More importantly, for all t ∈ [0, 1] we can show that the Hamiltonian above can be seen as the one emerging for an appropriate inference problem, so that the Nishimori property is kept valid for generic t ∈ [0, 1] [START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF].

Given the interpolating Hamiltonian we can write the corresponding Gibbs measure,

P (x|W, Y, T ) = 1 Z N,t P X (x)e H N,t (x) , (46)
and the interpolating free entropy

ψ N (t) . = 1 N E log Z N,t , (47)
whose boundaries are ψ N (1) = 1 N F N (our target) and ψ N (0) = 1 N ΦRS + 1 4∆ 2 m 2 + p-1 2p∆p m p . We then use the fundamental theorem of calculus to write

F N = ψ N (1) = ψ N (0) + 1 N E 1 0 - ∂ log Z N,t ∂t dt . =R . (48)
We work with the second term and use Stein's lemma which, given a well behaving function g, provides the useful relation for a canonical Gaussian variable Z:

E Z [Zg(Z)] = E Z [g (Z)]. This yields R = -E 1 0 1 Z N,t dx N ∂H N,t (x) ∂t P X (x)e H N,t (x) dt = -E 1 0 ∂H N,t (x) ∂t dt = -E 1 0 i<j 1 ∆ 2 N (x * i x i x * j x j ) + i 1 <•••<ip (p -1)! ∆ 2 N p-1 (x * i 1 x i 1 . . . x * ip x ip ) - i m 2∆ 2 + m p-1 2∆ p x * i x i dt = E 1 0   1 4∆ 2 x • x * N 2 -2m x • x * N + 1 2p∆ p x • x * N p -pm p-1 x • x * N   dt .
where we have used the Nishimori property to replace terms such as x 2 by xx * . At this point, we can write

R = E 1 0   1 4∆ 2 x • x * N 2 -2m x • x * N   dt + E 1 0   1 2p∆ p x • x * N p -pm p-1 x • x * N   dt = - m 2 4∆ 2 + 1 4∆ 2 E 1 0 1 4∆ 2 x • x * N -m 2 dt + 1 2p∆ p E 1 0 x • x * N p -pm p-1 x • x * N dt . (49)
The first integral is clearly positive. The second one, however, seems harder to estimate. We may, however, use a simple convexity argument on the function f (x) = x k . Indeed observe that ∀ a, b ≥ 0 and p ≥ 1: a ppb p-1 a ≥ (1 -p)b p . We would like to use this property but there is the subtlety that we need x • x * to be non-negative. To bypass this problem we can add again a small perturbation that forces x • x * to concentrate around a non-negative value, without affecting the "interpolating free entropy" ψ N (t) in the N → ∞ limit. This is, again, the argument used in [START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF] and originally in [START_REF] Babu | Exact solution of the gauge symmetric p-spin glass model on a complete graph[END_REF]. In this way we can write

R ≥ - m 2 4∆ 2 + E 1 0   1 4∆ 2 x • x * N 2 -2m x • x * N   dt + (1 -p)m p 4∆ 2 ≥ - m 2 4∆ 2 - (p -1)m p 4∆ 2 . ( 50
)
This concludes the proof and yields the lower bound:

F N ≥ ψ N (0) - 1 4∆ 2 m 2 - p -1 2p∆ p m p = 1 N ΦRS (m) , (51)
so that for all m ≥ 0 lim inf

N →∞ F N = lim inf N →∞ ψ N (1) = lim inf N →∞ ψ N (0) + 1 0 ψ N (t)dt ≥ ΦRS (m) .
Upper bound: Aizenman-Sims-Starr scheme. The matching upper bound is obtained using the Aizenman-Sims-Starr scheme [START_REF] Michael Aizenman | Extended variational principle for the sherringtonkirkpatrick spin-glass model[END_REF] . This is a particularly effective tool that has been already used for these problems, see for example [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF][START_REF] Coja-Oghlan | Information-theoretic thresholds from the cavity method[END_REF][START_REF] Thibault Lesieur | Statistical and computational phase transitions in spiked tensor estimation[END_REF]. The method goes as follows. Consider the original system with N variables, H N and add an new variable x 0 so that we get an Hamiltonian H N +1 . Define the Gibbs measures of the two systems, the first with N variables and the second with N + 1 variables, and consider the two relative free entropies.

Call

A N = E [log Z N +1 ] -E [log Z N ] their difference.
First, we notice that we have lim sup N F N ≤ lim sup N A N because

F N = E 1 N log Z N = 1 N E log Z N Z N -1 Z N -1 Z N -2 . . . Z 1 Z 0 = 1 N i A i ≤ sup i A i .
Moreover, we can separate the contribution of the additional variable in the Hamiltonian H N +1 so that H N +1 = HN + x 0 z(x) + x 2 0 s(x), with x = (x 1 , . . . , x N ), and

z(x) = 1 ∆ 2 (N + 1) N i=1 Z 0i x i + (p -1)! ∆ p (N + 1) (p-1)/2 1≤i 1 <•••<i p-1 ≤N Z 0i 1 ...i p-1 x i 1 . . . x i p-1 + + 1 ∆ 2 (N + 1) N i=1 x * 0 x * i x i + (p -1)! ∆ p (N + 1) p-1 1≤i 1 <•••<i p-1 ≤N x * 0 x * i 1 x i 1 . . . x * i p-1 x i p-1 s(x) = - 1 2∆ 2 (N + 1) N i=1 x 2 i - (p -1)! 2∆ p (N + 1) p-1 1≤i 1 <•••<i p-1 ≤N (x i 1 . . . x i p-1 ) 2
and H N +1 is the same expression as Eq. ( 17) where the N in the denominators are replaced by N + 1. We rewrite also H N (x) as a perturbation of HN :

H N (x) = HN (x) + y(x) + O(1) with y(x) = 1 √ ∆ 2 N i<j V ij x i x j + p -1 (p -1)! ∆ p N p/2 i 1 <•••<ip V i 1 ...ip x i 1 . . . x ip + + 1 N 2 i<j x * i x i x * j x j - 1 2 (x i x j ) 2 + (p -1)! p -1 N p i 1 <•••<ip x * i 1 x i 1 . . . x * ip x ip - 1 2 (x i 1 . . . x ip ) 2 ,
where the Zs and the V s are standard Gaussian random variables.

Finally we can observe the partition functions Z N can be interpreted as ensemble averages with respect to HN . Thus A N = E log P X (x 0 )e x 0 z(x)+x 2 0 s(x) dx 0 HN -E log e y(x)

HN

. Now, using the Nishimori property and the concentration of the overlap around a non-negative value -that we denote m(Y, T ) since it depends explicitly on the disorder-it yields (see [START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF], see section 4.3 for details) eq. ( 36) in the thermodynamic limit, with m(Y, T ) instead of m. From this, we can now obtain the upper bound that concludes the proof: The phase diagram identifies four regions: easy (green), impossible (red), and hard (orange). The lines correspond to different phase transitions namely the stability threshold (dashed black), the information theoretic threshold (solid red), the algorithmic threshold (solid cyan), and the dynamical threshold (dotted orange). The vertical cuts represent the section along which the magnetization is plotted in Fig. 9. On the right: Phase diagram of the spiked matrixtensor model for p = 4. The main difference with respect to case p = 3, is that the algorithmic spinodal (solid cyan) is strictly above the stability threshold (dashed black). The hybrid-hard phase appears between these two lines (combined green and orange color). The vertical cuts represent the section along which the magnetization is plotted in Fig. 10.

lim sup N F N ≤ lim sup N A N ≤ lim sup N E Y,T ΦRS [m(Y, T )] ≤ lim sup N sup m Φ RS (m) ≤ ΦRS . (52
The dynamical evolution of the AMP algorithm in the large N limit is described by the so-called State Evolution (SE) equations. The derivation of these equations can be straightforwardly done using the same techniques developed in [START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF]. They can be written in terms of two dynamical order parameters namely m t = i xt i x * i /N , which encodes for the alignment of the current estimation xt i of the components of the signal with the signal itself at time t and q t = i xt i xt i /N . Finally, using the Nishimori symmetry it can be shown that m t = q t at all times3 , see e.g. [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF], and therefore the evolution of the algorithm is characterized by a single order parameter m t whose dynamical evolution is given by

m t+1 = 1 - 1 1 + m t ∆ 2 + (m t ) p-1 ∆p . ( 53
)
If we initialize the configuration of the estimator x at random, the initial value of m will be equal to zero on average. However finite size fluctuations will produce by chance a small bias towards the signal and therefore it is more meaningful to consider the initialization to be m t=0 = being an arbitrarily small positive number. We will call m AMP the fixed point of Eq. ( 53) reached from this infinitesimal initialization. The mean-square-error (MSE) reached by AMP after convergence is then given by MSE AMP = 1 -m AMP . We underline that Eq. ( 53) can be proven rigorously following [START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF][START_REF] Richard | A statistical model for tensor PCA[END_REF]. Finally we note that the fixed point of the SE satisfies the very same Eq. ( 38) that gives the replica free entropy. In the rest of this section we will study the fixed points of Eq. ( 53). This will allow too determine the phase diagram of the spiked matrix-tensor model.

We start by observing that m = 0 is a fixed point of Eq. [START_REF] Georgii | Gibbs measures and phase transitions[END_REF]. However in order to understand wether it is a possible attractor of the AMP dynamics we need to understand its local stability. This can be obtained perturbatively by expanding Eq. ( 53) around m = 0

m t+1 = m t ∆ 2 + m t ∆ 2 2 - (m t ) p-1 ∆ p + O (m t ) 3 . (54)
It is clear that the non-informative fixed point m = 0 is stable as long as ∆ 2 > 1. We will call ∆ 2 = 1 the stability threshold.

When p = 3 the SE equations are particularly simple and the fixed points are written explicitly as

m 0 = 0 ; m ± = 1 2   1 - ∆ 3 ∆ 2 ± 1 + ∆ 3 ∆ 2 2 -4∆ 3   . (55)
In the regime where ∆ 2 > 1, m 0 and m + are stable while m -in unstable. When ∆ 2 becomes smaller than one, m + becomes the only non-negative stable solution and therefore ∆ 2 = 1 is also known as the algorithmic spinodal since it corresponds to the point where the AMP algorithm converges to the informative fixed point. The informative solution m + exists as long as ∆ 2 ≤ ∆ dyn 2 , where we have defined the dynamical spinodal by

∆ dyn 2 = ∆ 3 2 √ ∆ 3 -1 . ( 56
)
For a generic p we cannot determine the values of the informative fixed points explicitly but we can easily study Eq. ( 53) numerically to get the full phase diagram.

Furthermore we can obtain the spinodal transition lines as follows. The key observation is that the two spinodals are critical points of the equation ∆ p (m; ∆ 2 ) where ∆ 2 is fixed, or analogously ∆ 2 (m; ∆ p ) where ∆ p is fixed (to have a pictorial representation of the idea you can see Fig. 10). We call x = m/∆ 2 + m p-1 /∆ p , and

f SE (x) ≡ 1 -1 1+x , then ∆ p ≡ ∆ p (x; ∆ 2 ) = f SE (x) p-1 x -f SE (x) ∆ 2 . ( 57
)
Then the stationary points are implicitly defined by

0 = d log ∆ p dm = ∂ log ∆ p ∂x (1 + x) 2 ∝ (p -1) f SE (x) f SE (x) - 1 - f SE (x) ∆ 2 x - f SE (x) ∆ 2 = 2-p ∆ 2 + (1 + x)(p -x -2) x(1 + x) x + 1 -1 ∆p , giving x ± (∆ 2 ) = 1 2 p -3 ± (p -1) 2 - 4 ∆ 2 (p -2) . (58)
Finally ∆ p x ± (∆ 2 ); ∆ 2 describes the two spinodals. We can also derive the tri-critical point, when the two spinodals meet, which is given by the zero discriminant condition on eq. ( 58)

∆ tri p ; 1/∆ tri 2 =     4(p -2) p-3 p-1 p-1 (p -3) 2 ; (p -1) 2 4(p -2)     . (59)

B.4 Phase diagrams of spiked matrix-tensor model

In this section we present the phase diagrams for the spiked matrix-tensor model as a function of the two noise levels ∆ 2 and ∆ p and for several values of p. These phase diagrams are plotted in Figs. 7 and8. Generically we can have four regions:

• Easy phase (green), where the MSE obtained through AMP coincides with the MMSE which is better than random sampling of the prior.

• Impossible phase (red), where the MMSE and MSE of AMP coincide and are equal to 1 (meaning that m * = m AM P =0).

• Hard phase (orange), where the MMSE is smaller than the MSE obtained from AMP and m * > m AMP ≥ 0.

• Hybrid-hard phase [START_REF] Ricci-Tersenghi | Typology of phase transitions in bayesian inference problems[END_REF] (mix of green and orange), is a part of the hard phase where the AMP performance is strictly better than random sampling from the prior, but still the MSE obtained this way does not match the MMSE, i.e. m * > m AMP > 0. The hybrid-hard phase can be found for p ≥ 4.

All these phases are separated by the following transition lines:

• The stability threshold (dashed black line) at ∆ 2 = 1 for all p. This corresponds to the point where the uninformative fixed point m = 0 looses its local stability.

• The information theoretic threshold (solid red line). Here m * > 0 and the MMSE jumps to a value strictly smaller than one.

• The algorithmic threshold (solid cyan line). This is where the fixed point of AMP jumps to the MMSE< 1. For p = 3 this line coincides with a segment of the stability threshold while for p ≥ 4 it is strictly above.

• The dynamic threshold (dotted orange line). Here the most informative fixed point (the one with largest m AMP ) disappears.

In Figs. 9, and 10 we plot the evolution of the magnetization m, as found through the fixed points of the SE equation, for several fixed values of ∆ p and p = 3 and p = 4, respectively. The values of ∆ p are identified by the vertical cuts in the phase diagrams of Fig. 7. The situation is qualitatively similar to Fig. 9, the difference being only the presence of the hybrid-hard phase. We can observe that when the transition is discontinuous, figure from (a) to (e), for 1/∆ 2 > 1.0 the uninformative solution becomes unstable and continuously goes to a stable-informative solution which is not the optimal one.

C Langevin Algorithm and its state evolution

The main goal of our analysis is to compare AMP with the performance of the Langevin dynamics. The advantage of the spiked matrix-tensor model is that in this case the Langevin dynamics can be studied in the large N limit through integro-differential equations for the correlation function, C(t, t ) = lim N →∞ i x i (t)x i (t ) /N , the response function R(t, t ) = lim N →∞

1 N i d x i (t)
dη i (t ) and the magnetization C(t) = lim N →∞ i x i (t)x * i /N . To obtain these equations we use the techniques developed in the context of mean-field spin glass systems [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Cugliandolo | Course 7: Dynamics of glassy systems[END_REF]. We call η i (t) a time dependent noise and we indicate with • the average with respect to it. The noise is Gaussian and characterized by η i (t) = 0 for all t and i = 1, . . . N and η i (t)η j (t ) = 2δ ij δ(tt ). As before we will denote by E[. . . ] the average with respect to the realization of disorder that in this case goes back to the specific realization of the signal.

Before proceeding, it is useful to introduce a set of auxiliary variables that will help in the following. For k ∈ {2, p} we define

r k ≡ r k (t) = 2/ kT k (t)∆ k , f k (x) = x k /2 and m(t) . = 1 N i x i (t)x * i , and the random variable ξi 1 ...i k ≡ 1 ∆ k ξ i 1 ...i k ∼ N (0, 1/∆ k ).
The time dependence in T k , will be used in the smart annealing protocol that will be used to avoid part of the Langevin hard phase. We introduce a time dependent Hamiltonian

H(t) = - 1 T 2 (t) √ N i<j ξij x i (t)x j (t) - (p -1)! T p (t)N p-1 2 i 1 <•••<ip ξi 1 ...ip x i 1 (t) . . . x ip (t) -N r 2 (t)f 2 (m(t)) -N r p (t)f p (m(t)) ,
and the associated Langevin dynamics

ẋi (t) = -µ(t)x i (t) - ∂H ∂x i (t) -η i (t) = -µ(t)x i (t) - 1 T 2 (t) √ N j( =i) ξij x j (t)+ + r 2 (t)f 2 (m(t)) - (p -1)! T p (t)N p-1 2 (i,i 1 ,...,i p-1 )\i ξii 1 ...i p-1 x i 1 (t) • • • x i p-1 (t) + r p (t)f p (m(t)) -η i (t) , (60)
with µ a Langrange multiplier that enforces the spherical constraint

N i=1 x 2 i (t) = N . If T k (t) = 1
for all k = 2, p, the stationary equilibrium distribution for the Langevin dynamics is given by the posterior measure. Using Ito's lemma one finds 1 N

d dt i x 2 i (t) = 2 N i x i (t) ẋi (t) + 2 .
Since the spherical constraint imposes the left-hand-side to be zero, one obtains a condition on the right-handside. By plugging the expression ( 60) in it, one gets that in the large N limit

µ(t) = 1 -2H 2 (t) -pH p (t) (61)
where

H k = - (k -1)! T k (t)N k-1 2 i 1 <•••<i k ξi 1 ...i k x i 1 (t) . . . x i k (t) -N r k (t)f k (m(t)) k = 2, p (62)
are the parts of the Hamiltonian defined in Eq. (C) relative to the matrix (k = 2) and to the tensor (k = p). Note that we have not specified any initial condition for the variables x i (t = 0). Therefore, since we always employ the spherical constraint, the initial condition for the dynamics is a point on the N dimensional hypersphere |x| 2 = N extracted with the flat measure.

In order to analyze the Langevin dynamics in the large N limit, we will use the dynamical cavity method [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Castellani | Spin glass theory for pedestrians[END_REF][START_REF] Agoritsas | Out-of-equilibrium dynamical mean-field equations for the perceptron model[END_REF]. We will consider a system of N variables, with N 1, and add a new one. This new variable will be considered as a small perturbation to the original system but at the same time will be treated self consistently.

C.1 Dynamical Mean-Field Equations

In the following we will drop the time dependence for simplicity restoring it only when it is needed. Given the system with N variables i = 1 . . . N , we add a new one, say i = 0, and define m = 1

N +1 N i=0 x i x * i 1 N N i=0 x i x * i
(henceforth we use the symbol to denote two quantities that are equal up to terms that vanish in the large-N limit). The Langevin equation associated to the new variable is

ẋ0 = -µx 0 - 1 T 2 (t) √ N j( =0) ξ0j x j + r 2 f 2 ( m) - (p -1)! T p (t)N p-1 2 (0,i 1 ,...,i p-1 )\0 ξ0i 1 ...i p-1 x i 1 • • • x i p-1 + r p f p ( m) -η 0 , (63)
where we used that N N + 1 for N 1. We will consider the contribution of the new variable on the others in perturbation theory. In the dynamical equations for the variables i = 1, . . . , N we can isolate the variable i = 0 and write

ẋi = -µx i - 1 T 2 (t) √ N j( =i,0) ξij x j + r 2 f 2 (m) - (p -1)! T p (t)N p-1 2 (i,i 1 ,...,i p-1 )\i,0 ξii 1 ...i p-1 x i 1 • • • x i p-1 + r p f p (m) -η i + H i , (64)
with

H i (t) = r 2 f 2 (m) + r p f p (m) 1 N x 0 - 1 T 2 (t) √ N ξ0i x 0 - (p -1)! T p (t)N p-1 2 (i,0,i 1 ,...,i p-2 )\i,0 ξi0i 1 ...i p-2 x 0 x i 1 • • • x i p-2 . ( 65
)
Consider the unperturbed variables x 0 i = x i H i =0 . At leading order in N we can write

x i x 0 i + t to dt δx i (t) δH i (t ) H i =0 H i (t ) . ( 66
)
In the dynamical equation for the variable 0 we can identify a piece associated to the unperturbed variables x 0 i . This term can be thought of collectively as a stochastic term Ξ(t)

ẋ0 = -µx 0 . = Ξ(t) - 1 T 2 (t) √ N j( =0) ξ0j x 0 j - (p -1)! T p (t)N p-1 2 (0,i 1 ,...,i p-1 )\0 ξ0i 1 ...i p-1 x 0 i 1 • • • x 0 i p-1 -η 0 + + r 2 f 2 (m) + r p f p (m) + r 2 f 2 (m) + r p f p (m) 1 N x 0 - 1 T 2 (t) √ N j( =0) ξ0j t to dt δx j (t) δH j (t ) H j =0 H j (t )+ - (p -1)! T p (t)N p-1 2 (0,i 1 ,...,i p-1 )\0 ξ0i 1 ...i p-1 t to dt δx i 1 (t) δH i 1 (t ) H i 1 =0 H i 1 (t )x 0 i 2 • • • x 0 i p-1 + permutations . ( 67
)
Indeed Ξ(t) encodes the effect of a kind of bath made by of the unperturbed variables i = 1, . . . , N to the new one. We can show that at leading order in N , Ξ(t) is a Gaussian noise with zero mean and variance given by

E Ξ(t)Ξ(t ) = 2δ(t -t ) -E   1 T 2 (t)T 2 (t )N j( =0) l( =0) ξ0j ξ0l x 0 j (t)x 0 l (t )   + -E    (p -1)! T p (t)T p (t )N p-1 (0,i 1 ,...,i p-1 )\0 (0,j 1 ,...,j p-1 )\0 ξ0i 1 ...i p-1 ξ0j 1 ...j p-1 x 0 i 1 • • • x 0 i p-1 x 0 j 1 • • • x 0 j p-1   
and the second term can be simplified as

E (p -1)! T p (t)T p (t )N p-1 (0,i 1 ,...,i p-1 )\0 (0,j 1 ,...,j p-1 )\0 ξ0i 1 ...i p-1 ξ0j 1 ...j p-1 x 0 i 1 • • • x 0 i p-1 x 0 j 1 • • • x 0 j p-1 = (p -1)! N p-1 1 T p (t)T p (t )∆ p (0,i 1 ,...,i p-1 )\0 x 0 i 1 (t)x 0 i 1 (t ) • • • x 0 i p-1 (t)x 0 i p-1 (t ) = 1 T p (t)T p (t )∆ p C p-1 (t, t ) ,
where we used (i 1 ,...,i k ) = 1 k! 1≤i 1 ,...,i k ≤N , we neglected terms sub-leading in N , and we used the definition of the dynamical correlation function

C(t, t ) = 1 N N i=1
x i (t)x i (t ) .

Therefore we have

E Ξ(t) = 0 ; (68) E Ξ(t)Ξ(t ) = 2δ(t -t ) + 1 T 2 (t)T 2 (t ) C(t, t ) + 1 T p (t)T p (t )∆ p C p-1 (t, t ) . ( 69
)
Now we can focus of the deterministic term coming from the first order perturbation in eq. ( 67). Consider just the integral for the p-body term, the other will be given by setting p = 2

(p -1)! T p (t)N p-1 2 (0,i 1 ,...,i p-1 )\0 ξ0i 1 ...i p-1 t to dt δx i 1 (t) δH i 1 (t ) H i 1 =0 H i 1 (t )x 0 i 2 • • • x 0 i p-1 + permutations = (p -1)! T p (t)N p-1 (0,i 1 ,...,i p-1 )\0 ξ2 0i 1 ...i p-1 t to dt 1 T p (t ) δx i 1 (t) δH i 1 (t ) H i 1 =0 x 0 i 1 (t)x 0 i 1 (t ) • • • x 0 i p-2 (t)x 0 i p-2 (t )x 0 (t )+ + permutations -p(p -1) 1 2T p (t)∆ p t to dt 1 T p (t ) R(t, t )C p-2 (t, t )x 0 (t ) ( 70
)
where we have used the definition of the response function

R(t, t ) = 1 N N i=1 δx i (t) δH i (t ) .
Plugging eq. ( 70) into eq. ( 67) we obtain an effective dynamical equation for the new variable in terms of the correlation and response function of the system with N variables

ẋ0 (t) = -µ(t)x 0 (t) + Ξ(t) + r p f p (C(t)) + r 2 f 2 (C(t))+ + (p -1) 1 T p (t)∆ p t to dt 1 T p (t ) R(t, t )C p-2 (t, t )x 0 (t ) + 1 T 2 (t)∆ 2 t to dt 1 T p (t ) R(t, t )x 0 (t ) . (71)
In order to close Eq. ( 71) we need to give the recipe to compute the correlation and response function.

C.2 Integro-differential equations

In order to obtain the final equations for dynamical order parameters we will assume that the new variable x 0 is a typical one, namely it has the same statistical nature of all the others. Therefore we can assume that

C(t, t ) . = E x 0 (t)x 0 (t ) R(t, t ) . = E δx 0 (t) δΞ(t ) C(t) . = E x 0 (t)x * 0 . (72)
Eqs. (72) give a way to obtain the equation for all the correlation functions. Indeed we can consider Eq. ( 71), multiply it by x 0 (t ), or differentiate it with respect to an external field h 0 (t ), or multiply it it by x * 0 and we can average the results over the disorder and thermal noise. Using the following identity

E Ξ(t)x 0 (t ) = DΞ(t) Ξ(t)x 0 (t )e -d td tΞ( t)K -1 ( t, t)Ξ( t) = = -dt DΞ(t) x 0 (t ) δ δΞ(t ) e -d td tΞ( t)K -1 ( t, t)Ξ( t) K(t, t ) = = dt E δx 0 (t ) δΞ(t ) K(t, t ) = dt R(t , t )K(t, t ) = = 2R(t , t) + 1 T p (t)∆ p t to dt 1 T p (t ) R(t , t )C p-1 (t, t ) + 1 T 2 (t)∆ 2 t to dt 1 T p (t ) R(t , t )C(t, t ) (73)
we get the following Langevin State Evolution (LSE) equations

∂ ∂t C(t, t ) = E ẋ0 (t)x 0 (t ) = 2R(t , t) -µ(t)C(t, t ) + r p (t)f p (C(t))C(t ) + r 2 (t)f 2 (C(t))C(t )+ + (p -1) 1 T p (t)∆ p t to dt 1 T p (t ) R(t, t )C p-2 (t, t )C(t , t )+ + 1 T p (t)∆ p t to dt 1 T p (t ) R(t , t )C p-1 (t, t )+ + 1 T 2 (t)∆ 2 t to dt 1 T 2 (t ) R(t, t )C(t , t ) + 1 T 2 (t)∆ 2 t to dt 1 T p (t ) R(t , t )C(t, t ) ; (74)
∂ ∂t R(t, t ) = E δ ẋ0 (t) δΞ(t ) = = δ(t -t ) -µ(t)R(t, t ) + (p -1) 1 T p (t)∆ p t t dt 1 T p (t ) R(t, t )R(t , t )C p-2 (t, t )+ + 1 T 2 (t)∆ 2 t t dt 1 T 2 (t ) R(t, t )R(t , t ) ; (75)
∂ ∂t C(t) = E ẋ0 (t)x * 0 = = -µ(t)C(t) + r p (t)f p (C(t)) + r 2 (t)f 2 (C(t))+ + (p -1) 1 T p (t)∆ p t to dt 1 T p (t ) R(t, t )C p-2 (t, t )C(t ) + 1 T 2 (t)∆ 2 t to dt 1 T 2 (t ) R(t, t )C(t ) ; (76)
µ(t) = 1 + r p (t)f p (C(t))C(t) + r 2 (t)f 2 (C(t))C(t)+ + p 1 T p (t)∆ p t to dt 1 T p (t ) R(t, t )C p-1 (t, t ) + 2 1 T 2 (t)∆ 2 t to dt 1 T 2 (t ) R(t, t )C(t, t ) . (77)
Note that the last equation for µ(t) is obtained by imposing the spherical constraint C(t, t) = 1 ∀t using the fact that 0 = dC(t,t) dt

= ∂C(t,t ) ∂t t =t + ∂C(t ,t) ∂t t =t
. The boundary conditions of this equations are: C(t, t) = 1 the spherical constrain, R(t, t) = 0 which comes from causality in the Itô approach and R(t, t → t -) = 1. The initial condition for C(0) = C 0 is the overlap with the initial configuration with the true signal. If the initial configuration is random, C 0 = 0 but will have finite size fluctuations, as in the case of AMP. Therefore we can think that C 0 = being an arbitrary small positive number.

D Numerical solution of the LSE equations

The dynamical equations (74-75-76-77) were integrated numerically using two schemes:

• fixed time-grid: the derivatives were discretized and integrated according to their causal structure. This method is suited only for short times (up to 500 time units);

• dynamic time-grid: the step size is doubled after a given number of steps and the equations are solved self-consistently for every waiting-time. This is the approach proposed in [START_REF] Kim | The dynamics of the spherical p-spin model: From microscopic to asymptotic[END_REF] and described in Appendix C of [START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. II. Model calculations and comparison to numerical simulations[END_REF]. It allows integration up to very large times (up to 10 6 time units).

The results of these algorithms are concisely reported in the phase diagram shown in the main paper. In what follows we will present the algorithms and a series of investigations that we carried out to check their stability, we will explain the procedure followed to delimit the Langevin hard region, and we will discuss how we can enter into part of that region by choosing a proper annealing protocol. The codes are available online [START_REF] Sarao Mannelli | Langevin state evolution integrators[END_REF].

D.1 Fixed time-grid (2 + p)-spin

In this approach time-derivatives and integrals were discretized using ∂ ∂t f (t, t )

1 ∆t f (t + ∆t, t ) -f (t, t
) , and the trapezoidal rule for integration t 0 f (t)dt ∆t 2 t/∆t-1 l=0 f (l∆t) + f ((l + 1)∆t) . For instance we defined a function for computing the update in the the response function, eq. ( 75) as follows

R(t + ∆t, t ) = R(t, t ) -∆t µ(t)R(t, t ) + 1 2 ∆t 2 ∆ 2 t/∆t-1 l=t /∆t R(t, l∆t)R(l∆t, t ) + R(t, (l + 1)∆t)R((l + 1)∆t, t ) + + (p -1) ∆t 2 ∆ p t/∆t-1 l=t /∆t C p-2 (t, l∆t)R(t, l∆t)R(l∆t, t ) + C p-2 (t, (l + 1)∆t)R(t, (l + 1)∆t)R((l + 1)∆t, t ) .
Analogously we defined the other integrators. A simple causal integration scheme, being careful with the Itô prescription, gives the pseudo-code below. Figure 11: Representation of the initialization and the first two iterations for the evaluation of a two-times observable using the dynamic-grid algorithm. The empty circles represent slots allocated in memory but not associated to any specific value, while the full circles are memory slots already associated. For any two time function, it first allocates the memory (a), than it fills half of the grid by linear propagation (b). Still using linear propagation it fills the slots with tt 1 (c), and it sets the other values by imposing self-consistency (d). Finally it halves the grid (e), doubles time step and it allocates the memory (f). Then the algorithm loops following the same scheme as in (b-c-d-e).

C(0, 0) ← 1; R(0, 0) ← 0; C(0) ← C 0 ; for t ≤ t max do C(t + ∆t, t + ∆t) ← 1; R(t + ∆t, t + ∆t) ← 0; µ(t) ← compute mu(C, R, C, t); C(t + ∆t) ← compute mag(µ, C, R, C, t); for t ≤ t do C(t + ∆t, t ) ← compute C(µ, C, R, C, t); R(t + ∆t, t ) ← compute R(µ, C, R
The numerical scheme we are going to discuss is presented in the Bayes-Optimal case where T 2 (t) ≡ T p (t) ≡ 1. However the derivation that we propose can be easily generalized to the case where the T s assume different values, but are constants 4 . It is convenient to manipulate the equations to obtain an equivalent set of equations for the functions

C(t, t ), Q(t, t ) . = 1 -C(t, t ) - t t R(t, t )dt , C(t)
, where Q(t, t ) represents the deviation from Fluctuation Dissipation Theorem (FDT) at time t starting from time t . Indeed when the FDT theorem holds, it states that R(t, t ) = -∂ t C(t, t ).

We briefly anticipate the strategy that the algorithm uses to solve the equations. The algorithm discretizes the times into N t intervals, first starting from the boundary conditions, C(t, t) = 1, Q(t, t) = 0 and C = C 0 ∈ [0; 1], it fills the grid for small times (or small time differences τ = tt 1) using linear propagation. Given a time t and the initial guess for the Lagrange multiplier obtained by the linear propagator, the integrals are discretized and evaluated, then the results is used to update the value of the Lagrange multiplier. This procedure is repeated iteratively until convergence. Once that the first grid is filled, it follows a coarse-graining procedure where the sizes of the time intervals is doubled and only half of the information is retained. This procedure is repeated a fixed number of doubling of the original grid. The doubling scheme allows to explore exponentially long times at the cost of loosing part of the information, the direct consequence of this is the loss of stability for very large times (especially when the functions C(t, t ), R(t, t ), C(t) undergo fast changes at large times).

Dynamical equations in the algorithm. We recall the function f k (x) = x k 2 and its derivatives,

f k (x) = kx k-1 2 and f k (x) = k(k-1)x k-2 2
. For simplicity in the notation, we introduce also

f k (t, t ) . = f k C(t, t ) (∂ t + µ(t))C(t, t ) = 2R(t , t) + r 2 C(t )f 2 (C(t)) + r p C(t )f p (C(t))+ + 1 ∆ 2 t 0 dt f 2 (t, t )R(t , t ) + 1 ∆ 2 t 0 dt f 2 (t, t )R(t, t )C(t , t )+ + 2 p∆ p t 0 dt f p (t, t )R(t , t ) + 2 p∆ p t 0 dt f p (t, t )R(t, t )C(t , t ) , (∂ t + µ(t))R(t, t ) = δ(t -t ) + 1 ∆ 2 t t dt f 2 (t, t )R(t, t )R(t , t )+ + 2 p∆ p t t dt f p (t, t )R(t, t )R(t , t ) , (∂ t + µ(t))C(t) = r 2 f 2 (C(t)) + r p f p (C(t))+ + 1 ∆ 2 t 0 dt f 2 (t, t )R(t, t )C(t ) + 2 p∆ p t 0 dt f p (t, t )R(t, t )C(t ) , µ(t) = 1 + r 2 C(t)f 2 (C(t)) + r p C(t)f p (C(t))+ + 2 ∆ 2 t 0 dt f 2 (t, t )R(t, t ) + 2 ∆ p t 0 dt f p (t, t )R(t, t ) .
Following the lines of [START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. II. Model calculations and comparison to numerical simulations[END_REF], we introduce the FDT violation function, Q(t, t ), and after some manipulation the systems becomes

(∂ t +µ(t))C(t, t ) = C(t ) r 2 f 2 (C(t)) + r p f p (C(t)) + + 1 ∆ 2 t 0 dt f 2 (t, t ) ∂Q(t , t ) ∂t + f 2 (t, t ) ∂Q(t, t ) ∂t C(t , t ) + - t t dt f 2 (t, t ) ∂C(t , t ) ∂t -f 2 (t, t ) ∂Q(t, t ) ∂t C(t , t ) + f 2 (1)C(t, t ) -f 2 (t, 0)C(t , 0) + + 2 p∆ p t 0 dt f p (t, t ) ∂Q(t , t ) ∂t + f p (t, t ) ∂Q(t, t ) ∂t C(t , t ) + - t t dt f p (t, t ) ∂C(t , t ) ∂t -f p (t, t ) ∂Q(t, t ) ∂t C(t , t ) + f p (1)C(t, t ) -f p (t, 0)C(t , 0) , (78)
(∂ t +µ(t))Q(t, t ) = µ(t) -1 + 1 ∆ 2 - t t dt f 2 (t, t ) ∂Q(t , t ) ∂t + t t dt f 2 (t, t ) ∂Q(t, t ) ∂t [Q(t , t ) -1]+ + f 2 (1)[Q(t, t ) -1] + f 2 (t, 0)C(t , 0) - t 0 dt f 2 (t,
(∂ t +µ(t))C(t) = r 2 f 2 (C(t)) + r p f p (C(t)) + 1 ∆ 2 f 2 (1)C(t) -f 2 (t, 0)C(0) - t 0 dt f 2 (t, (79)
µ(t) = 1 + r 2 C(t)f 2 (C(t)) + r p C(t)f p (C(t)) + 2 ∆ 2 [f 2 (1) -f 2 (t, 0)] + 2 ∆ p [f p (1) -f p (t, 0)]+ + t 0 dt 2 ∆ 2 f 2 (t, t ) + 2 ∆ p f p (t, t ) ∂Q(t, t ) ∂t , (80)
further simplifications can be obtained introducing

µ (t) = µ(t) -2 ∆ 2 f 2 (1) -2 ∆p f p (1)
µ (t) = 1 + r 2 C(t)f 2 (C(t)) + r p C(t)f p (C(t)) - 2 ∆ 2 f 2 (t,
First order expansion coefficients. In the numerics we will initialize the grid by a linear propagation of the initial conditions. To determine the coefficients to use we can expand the functions up the second term for small values of τ (and in the last equation of t)

C(t + τ, t ) = C(t , t ) + C (1,0) (t , t )τ + 1 2 C (2,0) (t , t ) + O(τ 3 ) , Q(t + τ, t ) = Q(t , t ) + Q (1,0) (t , t )τ + 1 2 Q (2,0) (t , t ) + O(τ 3 ) , C(t) = C(0) + C (1) (0)τ + 1 2 C (2) (0) + O(t 3 ) . (86)
This gives the following coefficients:

C(t, t) = 1, C (1,0) (t, t) = -1, Q(t, t) = 0, Q (1,0) (t, t) = 0, C(0) = C 0 and C (1) (0) = r 2 f 2 (C 0 ) + r p f p (C 0 ) (1 -(C 0 ) 2 ) -C 0
, where C 0 is the initial value of the overlap with the signal.

Numerical integration and derivation. The set of equations derived above presents six types of integrals 2. Fill the grid (small τ ). Continue to propagate the values for small time differences τ = tt 1. In terms of the algorithm it means that we have some elements of the grid, N c of them, close to the diagonal that will be updated by linear propagation because the approximation of small τ is still valid. In our simulation the first ∆t is of the order 10 -7 and N c = 2.

I (1AB) ij = t i t j dt A(t i ,
3. Fill the grid (larger τ ). The rest of the values will be copied from the previous t (A t+∆t,t = A t,t ).

These values are the initial guess for solving the self-consistent equations (94-95-96) and (97), in this procedure the derivatives are updated using the 2nd order discretization.

4.

Half the grid and expand. The grid is decimated which means that each observable is contracted A i,j ← A 2i,2j and the derivate are updated as follows dA

(h) i,j ← 1 2 (dA (h) 2i,2j + dA (h) 2i-1,2j ), dA (v)
i,j ← 1 2 (dA

(v) 2i,2j + dA (v)
2i,2j-1 ). The new time step is now: ∆t ← 2∆t. 5. Start over from step 2.

D.3 Numerical checks on the dynamical algorithm

The dynamic-grid algorithm has been checked in a variety of ways. Cross-checking using the fixed-grid algorithm. For short times the dynamical equations were solved using the fixed-grid algorithm and compared with the outcome of the dynamic-grid algorithm, obtaining the same results, see Fig. 12. In the figure we used the fixed-grid with t max = 100 and the ∆t = 6.25 * 10 -3 .

Same magnetization in the easy region.

In the impossible and easy regions, the overlap with the signal of both AMP and dynamic-grid integration, converges to the same value. In Fig. 13 we show the overlap obtained with AMP, black dashed line, and the overlap achieved by the integration scheme at a given time. We can see that the overlap with the signal as obtained solving the LSE equations converges to the same value of the fixed point of AMP. Given a fixed ∆ 2 we can observe that the time to convergence increase very rapidly as we decrease ∆ p . We fitted this increase of the relaxation time to get the boundary of the Langevin hard region. Dynamical transition. The dynamical transition where the finite magnetization fixed point disappears can be regarded as a clustering or dynamical glass transition. Indeed coming from the impossible phase, going towards the hard phase, at the dynamical transition the free energy landscape changes and the unique ergodic paramagnetic minimum of the impossible phase gets clustered into an exponential number of metastable glassy states (see Sec. E). Correspondingly the relaxation time of the Langevin algorithm diverges. Fitting this divergence with a power law we obtain an alternative estimation of the dynamical line. In the right panel of Fig. 14 we plot with yellow points the dynamical transition line as extracted from the fit of the relaxation time of the Langevin algorithm extracted coming from the impossible phase and entering in the hard phase. (circles) and conversely, at fixed ∆ 2 and decreasing ∆ p (stars). The blue dashed line marks a region above which we do not observe anymore a stable positive 1RSB complexity. Finally we plot with orange and yellow dots the dynamical transition line as extracted from the relaxation time of the Langevin algorithm coming from respectively the hard and impossible phase. On the right: phase diagram of the spiked matrix-tensor model for p = 4 as presented in the right panel of Fig. 7 with the additional Langevin hard phase boundary. Also in this case we observe that Langevin hard phase extends in the AMP easy phase. Interestingly the Langevin hard phase here folds and closes near the tri-critical point. The blue dashed line marks a region above which we do not observe anymore a stable positive 1RSB complexity.

D.4 Extrapolation procedure

In order to determine the Langevin Hard region, given a fixed value of ∆ p (∆ 2 ), we measure the relaxation time that it takes to relax to equilibrium. On approaching the Langevin hard region, this relaxation time increases and we extrapolate the divergence to obtain the critical ∆ * p (∆ * 2 respectively) where the relaxation time diverges. The extrapolation is done assuming a power law divergence. Fig. 14 shows the results of this procedure for the cases 2 + 3 and 2 + 4.

Numerical checks on the extrapolation procedure. To test the quality of the fits we use a similar numerical procedure to locate the spinodal of the informative solution, which is given by the points where the informative solution ceases to exist. This spinodal must be the same for both the AMP and the Langevin algorithm [START_REF] Zdeborová | Statistical physics of inference: thresholds and algorithms[END_REF].

Since we aim at studying the spinodal of the informative solution, we initialize the LSE with C 0 = 1 and let it relax, measuring the time it takes to equilibrate at the value of C given by the informative fixed point of AMP. We do this fixing ∆ 2 and changing ∆ p . As we approach the critical ∆ p,dyn , the relaxation time will diverge and we can fit this divergence with a power law. The dynamic threshold extracted in this way is finally compared with the one obtained from AMP. In Fig. 15 we show how this scheme has been applied for ∆ 2 ∈ {1.01, 1.05, 1.10, 2.00}. As we get closer to the critical line ∆ 2 = 1 the relaxation time increases (and the height of the plateau decreases), making the fit harder. All in all, we observe a very good agreement between the points found with these extrapolation procedures and the prediction obtaine with AMP, as shown in Fig. 14.

D.5 Annealing protocol

In this section we show that using specific protocols we are able to enter in the Langevin hard region. A generic annealing scheme would lower the noises of both the channels simultaneously, which, as we tested and discussed in the main text, will not be able to avoid the Langevin hard region. Instead we can use the following protocol

T 2 ≡ 1 , T p ≡ T p (t) = 1 + C ∆ p e -t τann . (105)
The constant C allows to select at the initial time the desired effective ∆ p,eff = ∆ p + Ce -t τann far from (and much larger than) the original one. Instead τ ann chooses the speed of the annealing protocol. Fig. 16 shows that using this protocol we are able to enter in Langevin Hard region even with ∆ 2 close to the AMP threshold. To this purpose we initiated the effective ∆ p,eff close to 100 (i.e. C = 100), very far from the Langevin hard region, and we used different speeds for the annealing of ∆ p (different colors in the figures). In the figures we can observe that approaching the ∆ 2 = 1 we need slower and slower protocols (larger and larger τ ann ). The reason for this behavior is due to the fact that approaching ∆ 2 = 1 with ∆ p = 100 a longer time is required to gain a non trivial overlap with the solution. Evidence of this growing timescale at ∆ p = 100 is given in Fig. 17 where we show the relaxation time for magnetizing to the solution varying ∆ 2 . In particular, we can observe that at ∆ 2 = 0.70 the relaxation time is of the order of 100 time units.

For the protocol to be successful it is therefore crucial that the annealing time τ ann is large enough to give the possibility of magnetizing the solution before ∆ p,eff has significantly decreased towards ∆ p . According to this analysis, it is not surprising that in Fig. 16 for ∆ 2 = 0.90 the proposed protocol seems not to be successful. For this value of the parameter ∆ 2 , the time to find a solution even with ∆ p = 100 should be larger than 1000 time units, which is much larger than the used τ ann and anyway out of the time window of our numerical solution. However, with an annealing time large enough it would be in principle possible to recover exactly the same boundaries of the AMP easy region.

E Glassy nature of the Langevin hard phase: the replica approach

In this section we study the landscape of the spiked matrix-tensor problem following the approach of [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF]. We underline here that we are interested in studying the free energy landscape problem rather than the energy landscape since the former is the relevant quantity for finite temperatures (β = 1 in our case, as discussed in Sec. A). The results of [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF] suggest that the AMP-hard phase and part of the AMP-easy phase are glassy. Therefore we could expect that low magnetization glassy states trap the Langevin algorithm and forbid the relaxation to the equilibrium configurations that surrounds the signal. This may happen also in a region where AMP instead is perfectly fine in producing configurations strongly correlated with the signal. In order to check this hypothesis we compute the logarithm of the number of glassy states, called the complexity by using the replica method [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF][START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF]. The goal of this analysis is to trace an additional line in the phase diagram that delimits the region where stable one step replica symmetry breaking (1RSB) metastable states exist. We conjecture that this provides a physical lower bound to the Langevin hard phase in the (∆ p , 1/∆ 2 ) phase diagram.

E.1 Computation of the complexity through the replica method

The replica trick is based on the simple identity: E log x = lim n→0 ∂ ∂n Ex n . Using this observation we can compute the expected value of the free energy, Φ = -(log Z)/N , averaging the Z n and taking the limit n → 0. This is in general as difficult as the initial problem, however, if we consider only integer n and extrapolate to 0, the computation becomes much less involved due to the fact that for integer n the average Ex n can be sometimes performed analytically. Indeed in this case the replicated partition function Z n can be regarded as the partition function of n identical uncoupled systems or replicas. Averaging over the disorder we obtain a clean system of interacting replicas. The Hamiltonian of this system displays an emerging replica symmetry since it is left unchanged by a permutation of replicas. This symmetry can be spontaneously broken in certain disordered models where frustration is sufficiently strong [START_REF] Mézard | Spin glass theory and beyond[END_REF].

In mean field models characterized by fully connected factor graphs, the resulting Hamiltonian of interacting replicas depends on the configuration of the system only through a simple order parameter, the overlap Q between them, which is a n × n matrix that describes the similarities of the configurations of different replicas in phase space. Furthermore the Hamiltonian is proportional to N which means that in the thermodynamic limit N → ∞, the model can be solved using the saddle point method. In this case one needs to consider a simple ansatz for the saddle point structure of the matrix Q that allows to take the analytic continuation for n → 0. The solution to this problem comes from spin glass theory and general details can be found in [START_REF] Mézard | Spin glass theory and beyond[END_REF]. The saddle point solutions for Q can be classified according to the replica symmetry breaking level going from the replica symmetric solution where replica symmetry is not spontaneously broken to various degree of spontaneous replica symmetry breaking (including full-replica symmetry breaking). Here we will not review this subject but the interested reader can find details in [START_REF] Mézard | Spin glass theory and beyond[END_REF]. The model we are analyzing can be studied in full generality at any degree of RSB (see for example [START_REF] Crisanti | Spherical 2+ p spin-glass model: An exactly solvable model for glass to spin-glass transition[END_REF][START_REF] Crisanti | Spherical 2+ p spin-glass model: An analytically solvable model with a glass-to-glass transition[END_REF][START_REF] Crisanti | Exactly solvable spin-glass models with ferromagnetic couplings: The spherical multi-p-spin model in a self-induced field[END_REF] where the same models have been studied in absence of a signal). However here we will limit ourselves to consider saddle point solutions up to a 1RSB level.

The complexity of the landscape can be directly related to replica symmetry breaking. A replica symmetric solution implies an ergodic free energy landscape characterized by a single pure state. When replica symmetry is broken instead, a large number of pure states arises and the phase space gets clustered in a hierarchical way [START_REF] Mézard | Spin glass theory and beyond[END_REF]. Making a 1RSB approximation means to look for a situation in which the hierarchical organization contains just one level: the phase space gets clustered into an exponential number of pure states with no further internal structure.

If we assume a 1RSB glassy landscape, we can compute the complexity of metastable states using a recipe due to Monasson [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF] (see also [START_REF] Zamponi | Mean field theory of spin glasses[END_REF] for a pedagogical introduction). The argument goes as follows.

Let us consider system with x real replicas infinitesimally coupled. If the free energy landscape is clustered into an exponential number of metastable states, the replicated partition function, namely the partition function of the system of x real replicas, can be written as

Z x e N [Σ(f * )-xβf * ]
where f * is the internal free energy of the dominant metastable states that is determined by the saddle point condition dΣ df (f * ) = βx and β the inverse temperature. Note that since we are interested in the Bayes optimal case, this corresponds to set β = 1. The function Σ(f ) is the complexity of metastable states having internal entropy f . Therefore, using the free parameter x we can reconstruct the form of Σ(f ) from the replicated free energy. In order to compute the replicated free energy we need to apply the replica trick on the replicated system, log Z x = lim n→0 ∂ ∂n (Z x ) n . Calling the replicated free energy Φ = -1 N log Z x , we get the complexity as Σ = x ∂Φ ∂x -Φ. We can now specify the computation to our case where the partition function is the normalization of the posterior measure. With simple manipulations of the equations [START_REF] Castellani | Spin glass theory for pedestrians[END_REF], the partition function can be expressed as the integral over the overlap matrix

(Z x ) n = Z n x ∝ ab dQ ab e N nxS(Q) lim sup Q e N nxS(Q) ; ( 106
)
where the overlap Q is a (nx + 1) × (nx + 1) matrix

Q =      1 m • • • m m Q . . . m     
that contains a special row and column that encodes the overlap between different replicas with the signal and therefore the corresponding overlap is the magnetization m. The 1RSB structure for the matrix Q can be obtained by defining the following nx × nx matrices: the identity matrix 1 ij = δ ij , the full matrix J a positive complexity as plotted in Fig. 18. At this point the equilibrium states that dominate the posterior measure are the so called threshold states for which the complexity is maximal. For those states the eigenvalue λ II = 0 which confirms that these states are marginally stable [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF]. Decreasing ∆ 2 one crosses the information theoretic phase transition where the relevant metastable states that dominate the posterior measure have zero complexity. This corresponds to a freezing/condensation/Kauzmann transition. Below the information theoretic phase transition the thermodynamics of the posterior measure is dominated by the state containing the signal. However one can neglect the high magnetization solution of the 1RSB equations to get the properties of the metastable branch and computing the complexity of states that have zero overlap with the signal. The complexity curves as a function of the Parisi parameter x for decreasing values of ∆ 2 are plotted in Fig. 18 for fixed ∆ p = 0.5 and several ∆ 2 . The curves contain a stable 1RSB part and an unstable one where λ II is negative. The 1RSB line shown in Figs. 14 is obtained by looking at when the states with positive complexity and λ II = 0 disappear. This means that it gives the point where the 1RSB marginally stable states disappear and therefore it is expected to be a lower bound for the disappearance of glassiness in the phase diagram. The important outcome of this analysis is that for ∆ 2 < 1 but not sufficiently small, namely in part of the AMP-easy phase, the replica analysis predicts the existence of 1RSB marginally stable glassy states that may trap the Langevin algorithm from relaxing towards the signal [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF] and therefore supports the existence of the Langevin hard phase. Finally in Fig. 19 we plot the complexity as a function of the internal free energy of the metastable states for some values of ∆ 2 and ∆ p .

E.2 Breakdown of the fluctuation-dissipation theorem in the Langevin hard phase

When the Langevin algorithm is able to reach equilibrium, being it the signal or the paramagnetic state, it should satisfy the Fluctuation-Dissipation Theorem (FDT) according to which the response function is related to the correlation function through R(t, t ) = -∂C(t,t ) ∂t . Furthermore, time translational invariance (TTI) should arise implying that both correlation and response functions should be functions of only the time difference meaning that R(t, t ) = R(tt ) and C(t, t ) = C(tt ) 5 . When the dynamics is run in the glass phase, metastable states may forbid equilibration. In this case time translational invariance is never reached at long times 6 and the dynamics displays aging violating at the same time the FDT relation. The analysis of the asymptotic aging dynamics has been cracked by Cugliandolo and Kurchan in [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF][START_REF] Lf Cugliandolo | Weak ergodicity breaking in mean-field spin-glass models[END_REF] (see also [START_REF] Cugliandolo | Course 7: Dynamics of glassy systems[END_REF] for a pedagogical review) in the simplest spin glass model (see also [START_REF] Lf Cugliandolo | On the out-of-equilibrium relaxation of the Sherrington-Kirkpatrick model[END_REF] for a much more complex situation) where no signal is present. The outcome of this work is that when the dynamics started from a random initial conditions is run in the glass phase, it drives the system to surf on the threshold states. In the model analyzed in [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF] these states correspond to the 1RSB marginally stable glassy states that maximize the complexity. In this section we analyze the Cugliandolo-Kurchan scenario by contrasting the numerical solution of the dynamical equations with the replica analysis of the complexity. According to [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF], the long time Langevin dynamics7 can be characterized by two time regimes. For short times differences tt ∼ O(1) and t → ∞, the system obeys the FDT theorem and TTI; this regime can be understood as a first fast local equilibration in the nearest metastable state available. On a longer timescale tt → ∞ and t/t < ∞, the dynamics surfs on threshold states and FDT and TTI are both violated. In this time window, both the response and correlation functions become functions of λ = h(t)/h(t ) being h(t) an arbitrary reparametrization of the time variable 8 . By defining C(λ) = C(t, t ) and R(λ) = tR(t, t ) the Cugliandolo-Kurchan solution implies that in this aging regime the FDT relation can be generalized to

R (λ) = x C (λ) (112)
with x an effective FDT ratio that controls how much the FDT is violated. In the scenario of [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF], the value of x coincides with the 1RSB Parisi parameter that corresponds to threshold states computed within the replica

Figure 1 :

 1 Figure1: Phase diagram of the spiked 2 + 3-spin model (matrix plus order 3 tensor are observed). In the easy (green) region AMP achieves the optimal error smaller than random pick from the prior. In the impossible region (red) the optimal error is as bad as random pick from the prior, and AMP achieves it as well. In the hard region (orange) the optimal error is low, but AMP does not find an estimator better than random pick from the prior. In the case of Langevin algorithm the performance is strictly worse than that for AMP in the sense that the hard region increases up to line depicted in green dots. The blue dashed-dotted line delimits the region of existence of stable 1RSB states.

Figure 2 :

 2 Figure 2: Evolution of the correlation with the signal C(t) in the Langevin algorithm at fixed noise on the matrix (∆ 2 = 0.7) and different noises on the tensor (∆ p ). As we approach the transition, estimated to ∆ * p 0.58, the time required to jump to the solution diverges. Inset: the behavior of C(t) as a function of the iteration time for the AMP algorithm for the same values of ∆ p .

10 Figure 3

 103 Figure3: Fit using a power law of the relaxation times of the Langevin algorithm, the point of divergence marks the limit of the Langevin easy region. These fits have been performed both, for fixed ∆ 2 = 0.7 (blue circles and yellow crosses) and for fixed ∆ p = 1 (pink circles). The circles are obtained with numerical solution of LSE that uses the dynamical grid while crosses are obtained using a fixed-grid (details in Sec. D.1).

Figure 4 :

 4 Figure 4: Correlation with the signal of AMP and Langevin at the kth iteration (at time t) for fixed ∆ 2 = 0.7.

  , τ ann = 10 C = 100, τ ann = 25 C = 100, τ ann = 50 C = 100, τ ann = 75 C = 100, τ ann = 100 C = 100, τ ann =

  i.d. samples (given Y ) from the distribution P (X = • | Y ), independently of every other random variables. Let us denote • the expectation with respect to P (X = • | Y ) and E the expectation with respect to (X, Y ). Then, for all continuous bounded function f

)B. 3 Figure 7 :

 37 Figure7: On the left: Phase diagram of the spiked matrix-tensor model for p = 3. The phase diagram identifies four regions: easy (green), impossible (red), and hard (orange). The lines correspond to different phase transitions namely the stability threshold (dashed black), the information theoretic threshold (solid red), the algorithmic threshold (solid cyan), and the dynamical threshold (dotted orange). The vertical cuts represent the section along which the magnetization is plotted in Fig.9. On the right: Phase diagram of the spiked matrixtensor model for p = 4. The main difference with respect to case p = 3, is that the algorithmic spinodal (solid cyan) is strictly above the stability threshold (dashed black). The hybrid-hard phase appears between these two lines (combined green and orange color). The vertical cuts represent the section along which the magnetization is plotted in Fig.10.

Figure 8 :

 8 Figure 8: On the left: Phase diagram of the spiked matrix-tensor model for p = 5. On the right: Phase diagram of the spiked matrix-tensor model for p = 10. In both cases we observe qualitatively the same scenario found in the right panel of Fig. 7.

1 Figure 9 : 3 Figure 10 :

 19310 Figure 9: Fixed points of Eq. (53) as a function of ∆ 2 for p = 3 and several fixed values of ∆ p . The values of ∆ p correspond to the vertical cuts in the left panel of Fig.7. Solid lines are stable fixed point, dashed lines are unstable fixed points. The blue line represent informative fixed points with positive overlap with the signal while the orange line represent a uninformative fixed points with no overlap with the signal. Starting from high ∆ 2 an informative fixed point appears at the dynamical threshold (vertical dashed line) but is energetically disfavored until the information theoretic threshold (vertical dotted line) and finally it becomes the only stable solution crossing the algorithmic threshold (vertical dotted-dashed line). When the transition is continuous the three vertical threshold lines merge and we have a single second order phase transition, which here occurs at ∆ p ≥ 1.

Figure 12 :

 12 Figure12: Evolution of the correlation with the signal at fixed ∆ 2 for different ∆ p . The dotted red line overlapping with other lines, is the same quantity evaluated using the fixed grid algorithm up to time 100. We have started the LSE from an informative initial condition.

Figure 13 :

 13 Figure13: Correlation with the signal of AMP (dotted lines) and Langevin (solid lines) at kth iteration and t time respectively. The black dashed line is the asymptotic value predicted with AMP. In the easy region, provided enough running time, Langevin dynamics finds the same alignment as AMP. The figures show qualitatively the same behaviour for different values of ∆ 2 .

Figure 14 :

 14 Figure14: On the left: phase diagram of the spiked matrix-tensor model for p = 3 as presented in the left panel of Fig.7with the additional boundary of the Langevin hard phase (green circles and stars). The boundary of the Langevin hard phase has been obtained by fitting the relaxation time at fixed ∆ p and increasing ∆ 2 (circles) and conversely, at fixed ∆ 2 and decreasing ∆ p (stars). The blue dashed line marks a region above which we do not observe anymore a stable positive 1RSB complexity. Finally we plot with orange and yellow dots the dynamical transition line as extracted from the relaxation time of the Langevin algorithm coming from respectively the hard and impossible phase. On the right: phase diagram of the spiked matrix-tensor model for p = 4 as presented in the right panel of Fig.7with the additional Langevin hard phase boundary. Also in this case we observe that Langevin hard phase extends in the AMP easy phase. Interestingly the Langevin hard phase here folds and closes near the tri-critical point. The blue dashed line marks a region above which we do not observe anymore a stable positive 1RSB complexity.

Figure 15 :

 15 Figure 15: Relaxation time obtained from the LSE starting from an informative initial condition C 0 = 1. The four cases refers to the 2 + 3 model and are fitted with a power law and the relaxation time appears to diverge very close to point predicted by AMP (the AMP prediction is given in the captions).

Figure 17 :

 17 Figure 17: Relaxation times of Langevin dynamics at ∆ p = 100 without using protocols.

  nx,ij = 1, andJ (1) nx = diag(J (0) x , . . . , J(0)x ) a block diagonal matrix

Figure 18 :

 18 Figure 18: Complexity with as a function of the Parisi parameter x for p = 3 on the line ∆ p = 0.5. The solid line characterizes the stable part of the complexity while the dashed line the unstable one.

Figure 19 :

 19 Figure19: The stable part of the 1RSB complexity as a function of the free energy for p = 3 and ∆ p = 0.5.

Note that the pointwise convergence of the algorithm depends on the situations.

And the fact that we are in the Bayes optimal setting where we know the statistical properties of the signal, namely the prior, and the statistical properties of the channels, namely ∆2 and ∆p.

Note that AMP satisfies the Nishimori property at all times while this condition is violated on the run by the Langevin dynamics. In that case the Nishimori symmetry is recovered only when equilibrium is reached and therefore it is violated when the Langevin algorithm gets trapped in the glass phase, see below.

Therefore we do not employ this algorithm to solve the LSE equations in the smart annealing protocol for which instead we use the fixed time-grid algorithm.

All one time quantities are constant in equilibrium.

It is supposed to be reached only on exponential timescales in the system size.

But still for times that are not exponentially large in the system size N .

The function h(t) must be a monotonously increasing function. The asymptotic reparametrization invariance is a key property of the dynamical equations[START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF].

Acknowledgments

We thank G. Folena and G. Ben Arous for precious discussions. We thank K. Miyazaki for sharing his code for the numerical integration of CHSCK equations. We acknowledge funding from the ERC under the European Union's Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe; from the French National Research Agency (ANR) grant PAIL; from "Investissements d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM) (SaMURai and StatPhysDisSys); and from the Simons Foundation (#454935, Giulio Biroli).

The integrals can be easily discretized In particular the 6 integrals become

(87)

where the superscript (v) and (h) represent the vertical (t ) and horizontal (t) derivatives in the discretized times, see Fig. 11 for an intuitive understanding. We also discretized the derivative using the last two time steps

Given the time indices i and j, we will define and evaluate the following quantities

plus the respective vertical and horizontal derivatives.

Calling

p∆p M p ii , the original dynamical equations are integrated as follow

In the systems we used İ to characterize the integrals where we remove from the sum the term present in the left-hand side (e.g. for C ij eq. 94). Using Simpson's integration formula we define the increments

and we determine µ as

with δµ initially set to 0.

Algorithm: Here we describe the main steps of the algorithm, pictorially represented Fig. 11. Discretize the time (t, t ) in N t (even) intervals, the results shown use N t = 1024.

1. Initialization. Fill the first N t /2 times by linear propagation of the value obtained from the perturbative analysis where the diagonal blocks J (0)

x have size x × x and are matrices full of 1. In this case the 1RSB ansatz for Q reads Q = (1 -q M )1 nx + (q M -q m )J (1) nx + q m J (0) nx . Using this ansatz we can compute S(Q) that is given by

We can observe that starting from this expression we can derive the RS free energy, eq. ( 36), taking q M = q m or equivalently in the limit x → 1. From eq. ( 107) we obtain the saddle point equations

The low magnetization solution to these equations gives the complexity of the metastable branch of the posterior measure which is given by

The free parameter x allows us to tune the free energy of the states of which we compute the complexity. Thus we can characterize the part of the phase diagram where an exponential number of states is present.

To complete the 1RSB analysis we must compute the stability of the 1RSB saddle point solution for Q. In particular it is important to compute the so called replicon eigenvalues that characterizes the instability of this solution towards further replica symmetry breaking. Following [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | Exactly solvable spin-glass models with ferromagnetic couplings: The spherical multi-p-spin model in a self-induced field[END_REF] we have two replicon eigenvalues given by

We can analyze what happens to the landscape when we fix ∆ p < 1 and we start from a large value of ∆ 2 < ∆ 2,dyn (∆ p ) and we decrease ∆ 2 . In this case for sufficiently high ∆ 2 and large enough ∆ p the system is in a paramagnetic phase and no glassy states are present. At the dynamical transition line instead we find The vertical dotted line is the point where we observe a kink, which we denote by C = qEA and should be equal to the saddle point value of q M as extracted from the 1RSB threshold states in the replica computation [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF]: qEA = 0.633 and q M = 0.638. For C smaller than q EA the FDT is violated and is replaced by a generalized version as in Eq. ( 112). We can obtain the value of the FDT ratio from a fit of the slope of the asymptotic curves for C < qEA . We obtain x = 0.397 which should be compared with the Parisi parameter that corresponds to 1RSB marginally stable states obtained from the replica computation that is x = 0.408. Right panel: parametric plot of the integrated response as a function of the correlation for p = 3 and ∆ 2 = 1.4 and ∆ p = 0.2. In this case the value of the FDT ratio extracted from fitting the data is x = 0.397 to be compared with the value of the Parisi parameter for the 1RSB threshold states that is x = 0.408. At the same time data gives qEA = 0.633 while the replica computation gives q M = 0.638. approach. In order to test this picture we follow Cugliandolo and Kurchan [START_REF] Leticia | Aging and effective temperatures in the low temperature modecoupling equations[END_REF] and we plot the integrated response F(t, t ) = -t t R(t, t )dt as a function of C(t, t ) in a parametric way. This is done in Fig. 20.

If FDT holds at all timescales, one should see a straight line with slope -1. Instead what we see in the Langevin hard phase is that for large values of t the curves approach asymptotically for t 1 two straight lines. For high values of C, meaning for short time differences tt ∼ O(1), the slope of the straight line is -1 which means that F = 1 -C as implied by the short time FDT relation. On longer timescales FDT is violated, confirming the glassiness of the Langevin hard phase. By doing a linear fit we can use the data plotted in Fig. 20 to estimate the FDT ratio x appearing in Eq. ( 112). This can be compared with the Parisi parameter x for which we have marginally stable 1RSB states. We find an overall very good agreement (data coming from the fit is reported in the caption of Fig. 20). The small discrepancy between the two values of x can be either due to the numerical accuracy in solving the dynamical equations as well as the possibility that the 1RSB threshold is not exactly the one that characterizes the long time dynamics. Further investigations are needed to clarify this point. Finally, according to [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a longrange spin-glass model[END_REF] the value of C at which the two straight line cross should coincide with the value of q M computed for the threshold states within the 1RSB solution. Again we find a very good agreement.