Stefano Sarao Mannelli

Florent Krzakala

Pierfrancesco Urbani

Lenka Zdeborová

Passed & Spurious: analysing descent algorithms and local minima in spiked matrix-tensor model

In this work we analyse quantitatively the interplay between the loss landscape and performance of descent algorithms in a prototypical inference problem, the spiked matrix-tensor model. We study a loss function that is the negative log-likelihood of the model. We analyse the number of local minima at a fixed distance from the signal/spike with the Kac-Rice formula, and locate trivialization of the landscape at large signal-to-noise ratios. We evaluate in a closed form the performance of a gradient flow algorithm using integro-differential PDEs as developed in physics of disordered systems for the Langevin dynamics. We analyze the performance of an approximate message passing algorithm estimating the maximum likelihood configuration via its state evolution. We conclude by comparing the above results: while we observe a drastic slow down of the gradient flow dynamics even in the region where the landscape is trivial, both the analyzed algorithms are shown to perform well even in the part of the region of parameters where spurious local minima are present.

Introduction

A central question in computational sciences is the algorithmic feasibility of optimization in high-dimensional non-convex landscapes. This question is particularly important in learning and inference problems where the value of the optimized function is not the ultimate criterium for quality of the result, instead the generalization error or the closeness to a ground-truth signal is more relevant.

Recent years brought a popular line of research into this question where various works show for a variety of systems that spurious local minima are not present in certain regimes of parameters and conclude that consequently optimization algorithms shall succeed, without the aim of being exhaustive these include [START_REF] Kawaguchi | Deep learning without poor local minima[END_REF][START_REF] Soudry | No bad local minima: Data independent training error guarantees for multilayer neural networks[END_REF][START_REF] Ge | Matrix completion has no spurious local minimum[END_REF][START_REF] Daniel | Topology and geometry of half-rectified network optimization[END_REF][START_REF] Bhojanapalli | Global optimality of local search for low rank matrix recovery[END_REF][START_REF] Park | Non-square matrix sensing without spurious local minima via the Burer-Monteiro approach[END_REF][START_REF] Simon S Du | Gradient descent learns one-hidden-layer cnn: Don't be afraid of spurious local minima[END_REF][START_REF] Ge | No spurious local minima in nonconvex low rank problems: A unified geometric analysis[END_REF][START_REF] Simon S Du | Gradient descent learns one-hidden-layer cnn: Don't be afraid of spurious local minima[END_REF][START_REF] Lu | Depth creates no bad local minima[END_REF]. The spuriosity of a minima is in some works defined by their distance from the global minima, in other works as local minimizers that lead to bad generalization or bad accuracy in reconstruction of the ground truth signal. These two notions are not always equivalent, and certainly the later is more relevant and will be used in the present work.

Many of the existing works stop at the statement that absence of spurious local minimizers leads to algorithmic feasibility and the presence of such spurious local minima leads to algorithmic difficulty, at least as far as gradient-descent-based algorithms are concerned. At the same time, even gradient-descent-based algorithms may be able to perform well even when spurious local minima are present. This is because the basins of attraction of the spurious minimas may be small and the dynamics might be able to avoid them. In the other direction, even if spurious local minima are absent, algorithms might take long time to find a minimizer for entropic reasons that in high-dimensional problems may play a crucial role.

Main results: In this work we provide a case-study of a high-dimensional inference problems -the spiked matrix-tensor model -for which we are able to describe and quantitatively compare the following:

• With the use of the Kac-Rice formula [START_REF] Yan V Fyodorov | Complexity of random energy landscapes, glass transition, and absolute value of the spectral determinant of random matrices[END_REF][START_REF] Ben Arous | The landscape of the spiked tensor model[END_REF] we compute the expected number of local minimizers of the associated likelihood at a given correlation with the ground truth signal.

• With the use of the recently introduced Langevin-state-evolution [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF], that is a generalization of an approach well known in physics for analysis of the Langevin dynamics [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 long-range spin-glass model[END_REF], we give a closed-form description for the behaviour of the gradient flow (GF) algorithm in the limit of large system sizes.

• We analyze the state evolution of the maximum-likelihood version of the approximate message passing algorithm (ML-AMP).

We show that the above two algorithms (GF and ML-AMP) achieve the same error in the regime where they succeed. That same value of the error is also deduced from the position of all the minima strongly correlated with the signal as obtained from the Kac-Rice approach (precise statement below). We quantify the region of parameters in which the two above algorithms succeed and show that, up to the degree of accurancy of the extrapolation of the GF performance, the two lines are not the same. We also show that the algorithmic performance is not driven by the absence of spurious local minima. These results are summarized in Fig. 1. These result show, that in order to obtain a complete picture for settings beyond the present model, the precise interplay between absence of spurious local minima and algorithmic performance remains to be further investigated.

Problem definition

In this paper we consider the spiked matrix-tensor model as studied in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF]. This is a statistical inference problem where the ground truth signal x * ∈ R N is sampled uniformly on the N -1-dimensional sphere, S N -1 ( √ N ). We (the full cyan line), Eq. ( 32), the gradient flow and the ML-AMP algorithm, respectively, converge close to the ground truth signal in time linear in the input size. While the results for Kac-Rice and ML-AMP are given in a closed form, the ones for GF are obtained by extrapolating a convergence time obtained by numerical solution of integro-differential equations that describe large size behaviour of the GF. We note that all the three lines ∆ triv 2 , ∆ GF 2 , and ∆ ML-AMP 2 converge to 1 as ∆ p → ∞, consistently with the spiked matrix model. These three lines, related to minimization of the landscape, and their mutual positions, are the main result of this paper. The colors in the background, separated by the black dashed-dotted lines, show for comparison the phase diagram for the Bayes-optimal inference, related to the ability to approximate the marginals of the corresponding posterior probability distribution, and are taken from [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF]. In the red region obtaining a positive correlation with the signal in information-theoretically impossible. In the green region it is possible to obtain optimal correlation with the signal using the Bayes-optimal AMP (BO-AMP). And in the orange the region the BO-AMP is not able to reach the Bayes-optimal performance. then obtain two types of observations about the signal, a symmetric matrix Y , and an order p symmetric tensor T , that given the signal x * are obtained as

Y ij = x * i x * j √ N + ξ ij , (1)
T i 1 ,...,ip = (p -1)! N (p-1)/2 x * i 1 . . . x * ip + ξ i 1 ,...,ip (2)
for 1 ≤ i < j ≤ N and 1 ≤ i 1 < • • • < i p ≤ N , using and symmetries to obtain the other non-diagonal components. Here ξ ij and ξ i 1 ,...,ip are for each i < j and each i 1 < • • • < i p independent Gaussian random numbers of zero mean and variance ∆ 2 and ∆ p , respectively. The goal in this spiked matrix-tensor inference problem is to estimate the signal x * from the knowledge of the matrix Y and tensor T . If only the matrix was present, this model reduces to well known model of low-rank perturbation of a random symmetric matrix, closely related to the spiked covariance model [START_REF] Johnstone | On the distribution of the largest eigenvalue in principal components analysis[END_REF]. If on the contrary only the tensor is observed then the above model reduces to the spiked tensor model as introduced in [START_REF] Richard | A statistical model for tensor PCA[END_REF] and studies in a range of subsequent papers.

In this paper we study the matrix-tensor model where the two observations are combined. Our motivation is similar to the one exposed in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF], that is, we aim to access a regime in which it is algorithmically tractable to obtain good performance with corresponding message passing algorithms yet it is challenging (e.g. leading to non-convex optimization) with sampling or gradient descent based algorithms, this happens when both ∆ 2 = Θ(1) and ∆ p = Θ(1), while N → ∞ [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF].

In this paper we focus on algorithms that aim to find the maximum likelihood estimator. The negative log-likelihood (Hamiltonian in the physics language, or loss function in the machine learning language) of the spiked matrix-tensor reads

L = i<j 1 2∆ 2 Y ij - x i x j √ N 2 + i 1 <•••<ip 1 2∆ p T i 1 ...ip - (p -1)! N (p-1)/2 x i 1 . . . x ip 2 , (3)
where x ∈ S N -1 ( √ N ) is constrained to the sphere. In a high-dimensional, N → ∞, noisy regime the maximum-likelihood estimator is not always optimal as it provides in general larger error than the Bayes-optimal estimator computing the marginals of the posterior, studied in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF]. At the same time the log-likelihood (3) can be seen as a loss function, that is non-convex and high-dimensional. The tractability and properties of such minimization problems are the most questioned in machine learning these days, and are worth detailed investigation in the present model. 19), is shown for different values of parameter ∆ 2 at fixed ∆ p = 4.0 in the case p = 4. As ∆ 2 is decreased (the signal to noise ratio increases) the complexity allows to identify three main scenarios in the topology of the loss landscape. In the first case (a) only a wide band of non-negative complexity around the point of zero correlation is present, in the second case (b) minima with non-trivial correlation with the signal appear but the band around m = 0 is still present, finally (c) the signal dominates over the noise and only minima with non-trivial correlation are present. The transition from case (b) to case (c), i.e. when the support of Σ(m) ≥ 0 becomes two discontinuous points, as the bulk close to m = 0 becomes negative, is called the landscape trivialization. The ∆ 2 at which this occurs is denoted ∆ triv 2 and depicted in dashed purple in Fig. 1.

The first goal of this paper is to characterize the structure of local minima of the loss function (equivalently local maxima of the log-likelihood) eq. ( 3) as a function of the noise parameters ∆ 2 and ∆ p . We compute the average number of local minimizers x having a given correlation with the ground truth signal m = lim N →∞ x•x * /N . This leads to a so-called complexity function Σ(m) defined as the logarithm of the expected number of local minima at correlation m with the ground truth.

A typical example of this function, resulting from our analysis, is depicted in Fig. 2 for p = 4, ∆ p = 4.0, and several values of ∆ 2 . We see from the figure that at large ∆ 2 local minima appear only in a narrow range of values of m close to zero, as ∆ 2 decreases the support of Σ(m) ≥ 0 widens. At yet smaller values of ∆ 2 the support Σ(m) ≥ 0 becomes disconnected so that it is supported on an interval of value close to m = 0 and on two (one negative, one positive) isolated points. For yet smaller ∆ 2 the complexity for values of m close to zero becomes negative, signalling what we call a trivialization of the landscape, where all remaining local minima are (in the leading order in N ) as correlated with the ground truth as the global minima. The support of Σ(m) ≥ 0 in the trivialized regions consists of two separated points. We call the value of ∆ 2 at which the trivialization happens ∆ triv 2 . In the phase diagram of Fig. 1 the trivialization of the energy landscape happens above the purple dashed line.

We use the Kac-Rice formula to determine the complexity Σ(m) [START_REF] Robert | Random fields and geometry[END_REF][START_REF] Yan V Fyodorov | High-dimensional random fields and random matrix theory[END_REF]. Given an arbitrary continuous function, the Kac counting formula allows to compute the number of points where the function crosses a given value. The number of minima can be characterized using Kac's formula on the gradient of the loss (3), counting how many time the gradient crosses the zero value, under the condition of having a positive definite Hessian in order to count only local minima and not saddles. Since the spiked matrix-tensor model is characterized by a random landscape, due to the noise ξ ij and ξ i 1 ,...,ip , we will consider the expected number of minima obtaining the Kac-Rice formula [START_REF] Robert | Random fields and geometry[END_REF][START_REF] Yan V Fyodorov | High-dimensional random fields and random matrix theory[END_REF].

For mathematical convenience we will consider the rescaled configurations σ = x/ √ N ∈ S N -1 (1), and rescaled signal

σ * = x * / √ N . Call φ G,F 2 ,
Fp the joint probability density of the gradient G of the loss, and of the F 2 and F p the contributions of the matrix and tensor to the loss, respectively. Given the value of the two contributions to the loss F 2 = 2 N and F p = p N , and the correlation between the configuration and ground truth m ∈ [-1, +1] that we impose using a Dirac's delta, the averaged number of minimizers is

N (m, 2 , p ; ∆ 2 , ∆ p ) = e Σ∆ 2 ,∆p (m, 2 , p) = = S N -1 E[det H|G = 0, F 2 = N 2 , F p = N p , H 0]φ G,F 2 ,Fp (σ, 0, 2 , p ) δ(m -σ • σ * ) dσ . (4)
Rewrite the loss Eq. ( 3) neglecting terms that are constant with respect to the configuration and thus do not contribute to the complexity

L = N (p -1)! ∆ p i 1 <•••<ip ξ i 1 ...ip σ i 1 . . . σ ip + N (p -1)! ∆ p i 1 <•••<ip σ * i 1 σ i 1 . . . σ * ip σ ip + √ N ∆ 2 i<j ξ ij σ i σ j + N ∆ 2 i<j σ * i σ i σ * j σ j . (5)
In the following we will use small letters f 2 , f p , g, h to characterize losses, gradient and Hessian constrained on the sphere and capital letters for the same quantities unconstrained. Define I d the d-dimensional identity matrix.

The following lemma characterizes φ G,F 2 ,Fp .

Lemma 1. Given the loss function Eq. ( 5) and a configuration x such that the correlation and the signal is m, then there exists a reference frame such that the joint probability distribution of 1) is given by

f 2 , f p ∈ R, g ∈ R N -1 and h ∈ R (N -1)×(N -
f k N ∼ 1 k∆ k m k + 1 √ k∆ k 1 √ N Z k ; (6) g N ∼ 1 ∆ p m p-1 + 1 ∆ 2 m 1 -m 2 e 1 - 1 ∆ p + 1 ∆ 2 1 √ N Z ; (7) h N ∼ p -1 ∆ p m p-2 + 1 ∆ 2 (1 -m 2 )e 1 e T 1 + p -1 ∆ p + 1 ∆ 2 N -1 N W -(pf p + 2f 2 )I N -1 ; (8)
with Z k standard Gaussians and k ∈ {2, p}, Z ∼ N (0, I N -1 ) a standard multivariate Gaussian and W ∼ GOE(N -1) a random matrix from the Gaussian orthogonal ensemble.

Proof sketch. Starting from Eq. ( 5), split the contributions of the matrix and tensor in F 2 and F p , two Gaussian variables and impose the spherical constrain with a Lagrange multiplier µ.

f 2 (σ) + f p (σ) = F 2 (σ) + F p (σ) - µ 2 i σ 2 i -1 , (9)
g i (σ) = G i (σ) -µσ i , (10)
h ij (σ) = H ij (σ) -µ . (11)
The expression for µ in a critical point can be derived as follows. Given g i (σ) ≡ 0, multiply Eq. ( 10) by σ i , sum over the indices and we obtain:

µ = i G i (σ)σ i = 2f 2 (σ) + pf p (σ)
. We now restrict our study to the unconstrained random variables and substitute µ. Since the quantities f 2 , f p , g, h, µ are linear functionals of Gaussians they will be distributed as Gaussian random variables and therefore can be characterized by computing expected values and covariances. Starting from the losses coming from the matrix and the tensor in Eq. ( 5), F 2 (σ) and F p (σ), respectively, consider the moments with respect to the realization of the noise, ξ i 1 ...ip , ξ ij . For k ∈ {2, p} the first moment leads to

E[F k (σ)] = N k∆ k (σ • σ * ) k + O(1) . (12)
Let's consider the second moment but having two different configurations σ and τ ,

E F k (σ)F k (τ ) = N k∆ k (σ • τ ) k + O(1) . ( 13
)
Using standard results for derivatives of Gaussians (see e.g. [START_REF] Robert | Random fields and geometry[END_REF] Eq. 5.5.4) we can obtain means and covariances of the random variables taking derivatives with respect to σ and τ . Then set τ = σ, imposing the spherical constrain and using σ • σ * = m. The last step is the definition of a convenient reference frame {e j } j=1,...,N . Align the configuration along the last coordinate e N = σ and the signal with a combination of the first and last coordinates σ * = √ 1 -m 2 e 1 +me N . Finally, project on the sphere by discarding the last coordinate.

We can now rewrite the determinant of the conditioned Hessian by grouping the multiplicative factor in front of the GOE in Eq. ( 8)

det h = p -1 ∆ p + 1 ∆ 2 N -1 2 N N -1 -N -1 2 det W -t N I N -1 + θ N e 1 e T 1 ( 14
)
with t N and θ N given by

t N → t = 2 p p + 2 2 p-1 ∆p + 1 ∆ 2 , ( 15
)
θ N → θ = p-1 ∆p m p-2 + 1 ∆ 2 p-1 ∆p + 1 ∆ 2 (1 -m 2 ) ( 16
)
in the large N -limit. Therefore the Hessian behaves like a GOE shifted by t with a rank one perturbation of strength θ. This exact same problem has already been studied in [START_REF] Ben Arous | The landscape of the spiked tensor model[END_REF] and we can thus deduce the expression for the complexity as

Σ∆ 2 ,∆p (m, 2 , p ) = 1 2 log p-1 ∆p + 1 ∆ 2 1 ∆p + 1 ∆ 2 + 1 2 log(1 -m 2 ) - 1 2 m p-1 ∆p + m ∆ 2 2 1 ∆p + 1 ∆ 2 (1 -m 2 ) -2p∆ p p - m p 2p∆ p 2 -4∆ 2 2 - m 2 4∆ 2 2 + Φ(t) -L(θ, t), (17)
with

Φ(t) = t 2 4 + 1 |t|>2   log t 2 4 -1 + |t| 2 - |t| 4 t 2 -4   L(θ, t) =                1 4 t θ+ 1 θ y 2 -4dy - θ 2 t -θ + 1 θ + t 2 -θ + 1 θ 2 8 θ > 1, 2 ≤ t < θ 2 + 1 θ ∞ t < 2 0 otherwise.
We note at this point that for the case of the pure spiked tensor model ∆ 2 → ∞ the above expression reduces exactly to the complexity derived in [START_REF] Ben Arous | The landscape of the spiked tensor model[END_REF]. The following theorem states that to the leading order Eq. ( 17) represents the complexity of our problem.

Theorem 1. Given ∆ 2 and ∆ p , for any

( 2 , p ) ∈ R 2 and m ∈ [-1, +1] it holds lim N →∞ 1 N log EN (m, 2 , p ; ∆ 2 , ∆ p ) = = Σ∆ 2 ,∆p (m, p , 2 ) (18)
Proof sketch. The proof comes immediately from [START_REF] Ben Arous | The landscape of the spiked tensor model[END_REF] Thm. 2, see also Sec. 4.1.

The quantity that we are interested in is the projection of Eq. ( 17) to the maximal values of 2 and p :

Σ(m) = max 2 , p Σ∆ 2 ,∆p (m, 2 , p ). (19)
Eq. ( 19) allows to understand if at a given correlation with the signal, there are regions with an exponential expected number of minima, see Fig. 2. Thus it allows to locate parameters where the landscapes is trivial. We computed the expected number of minima, i.e. the so-called annealed average. The annealed average might be dominated by rare samples, and in general provides only an upper bound for typical samples. The quenched complexity, i.e. the average of the logarithm of the number of minima, is more involved. The quenched calculation was done in the case of a the spiked tensor model [START_REF] Ros | Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions[END_REF]. It is interesting to notice that in [START_REF] Ros | Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions[END_REF] the authors found that the annealed complexity does not differ from the quenched complexity for m = 0. This combined with analogous preliminary results for the spiked matrix-tensor model, suggest that considering the quenched complexity would not change the conclusions of this paper presented in the phase diagrams Fig. 1.

Gradient flow analysis

In this section we analyze the performance of the gradient flow descent in the loss function (3)

d dt x i (t) = -µ(t)x i (t) - δL δx i (t) , (20)
where the Lagrange parameter µ(t) is set in a way to ensure the spherical constraint x ∈ S N -1 ( √ N ). Our aim is to understand the final correlation between the ground truth signal and the configuration reached by the gradient flow in large but finite time, while N → ∞.

The gradient flow [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF] can be seen as a zero-temperature limit of the Langevin algorithm where 24) characterizes the evolution of the correlation of the gradient flow with the ground truth signal, evaluated for several values of ∆ 2 , at ∆ p = 4.0 starting from m(0) = 10 -10 . The dynamics displays a fast increase of the convergence time as ∆ 2 increases. At large times, the plateau we observe has the same value of correlation m as the minima best correlated with the signal, as predicted via Kac-Rice approach.

d dt x i (t) = -µ(t)x i (t) - δL δx i (t) -η i (t) , (21)
with η i (t) being the Langevin noise with zero mean and covariance η i (t)η j (t ) = 2T δ ij δ(t -t ), where T has the physical meaning of temperature, the notation . . . stands for the average over the noises ξ ij and ξ i 1 ,...,ip . As we take the limit T → 0, the noise becomes peaked around zero, effectively recovering the gradient flow.

The performance of the Langevin algorithm was characterized recently in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] using equations developed in physics of disordered systems [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 long-range spin-glass model[END_REF]. In [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] this characterization was given for an arbitrary temperature T and compared to the landscape of the Bayes-optimal estimator [START_REF] Antenucci | On the glassy nature of the hard phase in inference problems[END_REF]. Here we hence summarize and use the results of [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] corresponding to the limit T → 0.

The Langevin dynamics with generic temperature is in the large size limit, N → ∞, characterized by a set of PDEs for the self-correlation C(t, t ) = lim N →∞

1 N x i (t)x i (t ) , the response function R(t, t ) = lim N →∞ 1 N δx i (t)
δη i (t ) , and the correlation with the signal m(t) = lim N →∞ 1 N

x i (t)x * i . Ref. [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] established that as the gradient flow evolves these quantities satisfy eqs. (74)-(76) in that paper. Taking the zero-temperature limit in those equations we obtain

∂ ∂t C(t, t ) = -μ(t)C(t, t ) + Q (m(t))m(t ) + t 0 dt R(t, t )Q (C(t, t ))C(t , t ) + t 0 dt R(t , t )Q (C(t, t )) , (22)
∂ ∂t R(t, t ) = -μ(t)R(t, t ) + t t dt R(t, t )Q (C(t, t ))R(t , t ) , (23)
∂ ∂t m(t) = -μ(t)m(t) + Q (m(t)) + t 0 dt R(t, t )m(t )Q(C(t, t )) , (24)
with Q(f ) = f p /(p∆ p ) + f 2 /(2∆ 2 ) and μ(t) = lim T →0 T µ(t) the rescaled spherical constraint. Boundary conditions for the equations are C(t, t) = 1 ∀t, R(t, t ) = 0 for all t < t and lim t →t -R(t, t ) = 1 ∀t. An additional equation for μ(t) is obtained by fixing C(t, t) = 1 in Eq. ( 22). In the context of disordered systems those equations have been established rigorously for a related case of the matrix-tensor model without the spike [START_REF] Ben Arous | Cugliandolo-Kurchan equations for dynamics of spin-glasses[END_REF].

Performance of the gradient flow

Eqs. [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF][START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF] are integrated numerically showing the large-size-limit performance of the gradient flow algorithm.

Example of this evolution is given in Fig. 3 for p = 3, ∆ p = 4. The code will be made available and linked to this paper. For consistency we confirm numerically that at large times the gradient flow reaches values of the correlation that correspond exactly to the value of the correlation of the minima correlated to the signal as obtained in the Kac-Rice approach.

As the variance ∆ 2 increases the time it takes to the gradient flow to acquire good correlation with the signal increases. We define the convergence time t c as the time it takes to reach 1/2 of the final plateau. The dependence of t c on ∆ 2 is consistent with a power law divergence at ∆ GF 2 . This is illustrated in Fig. 4 where we plot the convergence time as a function of ∆ 2 and show the power-law fit in the inset. The points ∆ GF 2 are collected and plotted in Fig. 1, dotted blue line.

From Fig. 4 we see that the gradient flow algorithm undergoes a considerable slow-down even in the region where the landscape is trivial, i.e. does not have spurious local minimizers. At the same time divergence of the convergence time happens only well inside the phase where spurious local minimizers do exist.

Maximum-likelihood approximate message passing

Approximate Message Passing (AMP) is a popular iterative algorithm [START_REF] David L Donoho | Message-passing algorithms for compressed sensing[END_REF] with a key advantage of being analyzable via its state evolution [START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF]. The maximum-likehood AMP (ML-AMP) algorithm studied in this paper is a generalization of AMP for the pure spiked tensor model from [START_REF] Richard | A statistical model for tensor PCA[END_REF] to the spiked matrix-tensor model. We will show that its fixed points correspond to stationary points of the loss function [START_REF] Soudry | No bad local minima: Data independent training error guarantees for multilayer neural networks[END_REF]. This should be contrasted with the Bayes-optimal AMP (BO-AMO) that was studied in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] and aims to approximate the marginals of the corresponding posterior probability distribution. The ML-AMP instead aims to estimate the maximum-likelihood solution, x. In information theory the BO-AMP would correspond to the sum-product algorithm, while the present one to the max-sum algorithm. In statistical physics language the BO-AMP corresponds to temperature one, while the present one to zero temperature. In the appendix, Sec. C, we provide a schematic derivation of the ML-AMP as a zero-temperature limit of the BO-AMP, using a scheme similar to [START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF]. 19). The vertical purple dashed line at 1/∆ triv 2 = 2.57 corresponds to the trivialization of the landscape, beyond which only local minima well correlated with the signal remain. We note that all these approaches agree on the value of the MSE up to numerical precision. For the sake of comparison we show (the dashed-dotted grey line) also the minimal-MSE achieved in the Bayes-optimal setting. The ML-AMP algorithm reads

B t i = (p -1)! N (p-1)/2 k 2 <•••<kp T ik 2 ...kp ∆ p xt k 2 . . . xt kp + 1 √ N k Y ik ∆ 2 xt k -r t xt-1 i , (25)
xt+1 i = B t i 1 √ N ||B t || 2 , ( 26
) σt+1 = 1 1 √ N ||B t || 2 (27) with || • • • || 2
2 the 2 -norm and r t the Onsager reaction term

r t = 1 ∆ 2 1 N k σt k + p -1 ∆ p 1 N k σt k   1 N k xt k xt-1 k   p-2 . ( 28
)

ML-AMP & stationary points of the loss

Using an argument similar to Prop. 5.1 in [START_REF] Montanari | Graphical models concepts in compressed sensing[END_REF] we can show that a fixed points found by ML-AMP corresponds to finding a stationary point of the loss Eq. ( 3) with a ridge regularizer.

Property 1. Given (x * , σ * ) a fixed point of ML-AMP, then x * satisfies the stationary condition of the loss.

Proof sketch. Let us denote B * , r * the fixed point of Eqs. ( 25) and (28). From Eq. ( 26) and Eq. ( 25) we have

1 √ N ||B * || 2 + r * x * = 1 √ N k Y ik ∆ 2 x * i + (p -1)! N (p-1)/2 k 2 <•••<kp T ik 2 ...kp ∆ p x * k 2 . . . x * kp (29)
which is exactly solution of the derivative of Eq. (3) with respect to x i when the spherical constraint is enforced by a Lagrange multiplier µ

0 = -µx i + 1 √ N k Y ik ∆ 2 x i + (p -1)! N (p-1)/2 k 2 <•••<kp T ik 2 ...kp ∆ p x k 2 . . . x kp .
Moreover ML-AMP by construction preserves the spherical constrain at every time iteration, thus the fixed point value of the Lagrange multiplier in (29) is the one leading to the spherical constraint.

State evolution

The evolution of ML-AMP can be tracked through a set of equations called state evolution (SE). The state evolution can be characterized via an order parameter:

m t = 1 N i xt i x *
i , the correlation of the ML-AMPestimator with the ground truth signal at time t. According to the SE, as proven in [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], this parameter evolves in the large N limit as

m t+1 = m t ∆ 2 + (m t ) p-1 ∆p 1 ∆ 2 + 1 ∆p + m t ∆ 2 + (m t ) p-1 ∆p 2 , (30)
and the mean square error correspondingly

MSE t = 2(1 -m t ). (31)
A derivation of this state evolution is presented in the appendix, Sec. C. Analysis of the simple scalar SE, Eq. (30), allows to identify the error reached by the ML-AMP algorithm. We first observe that m = 0 is always a fixed point. For the performance of ML-AMP is the stability of this fixed point that determines whether the ML-AMP will be able to find a positive correlation with the signal or not. Analyzing Eq. ( 30) we obtain that the m = 0 is a stable fixed point for ∆ 2 > ∆ ML-AMP 2 where

∆ ML-AMP 2 (∆ p ) = -∆ p + ∆ 2 p + 4∆ p 2 . ( 32
)
Consequently for ∆ 2 > ∆ ML-AMP 2 the ML-AMP algorithm converges to m = 0, i.e. zero correlation with the signal. The line ∆ ML-AMP 2 is the line plotted in Fig. 1. For p = 3 and p = 4, we obtain that for ∆ 2 < ∆ ML-AMP 2 the ML-AMP algorithm converges to a positive m * > 0 correlation with the signal, depicted in Fig. 5. In Fig. 5 we also observe that this correlation agrees (up to the numerical precision) with the position of the maximum having largest value of m in the complexity function Σ(m), this is also depicted in the figure. The trivialization of the landscape occurs at ∆ triv 2 < ∆ ML-AMP 2 , thus showing that for ∆ triv 2 < ∆ < ∆ ML-AMP 2 the ML-AMP algorithm is able to ignore a good portion of the spurious local minima and to converge to the local minima best correlated with the signal.

In Fig. 5 we also compared to the MSE obtained by the Bayes-optimal AMP that provably minimizes the MSE in the case depicted in the figure [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF]. We see that the gap between the Bayes-optimal error and the one reached by the loss minimization approaches goes rapidly to zero as ∆ 2 increases.

Discussion

We analyzed the behavior of two descent algorithms in optimizing a rough high-dimensional loss landscape of the spiked matrix-tensor model. We used the Kac-Rice formula to count the average number of minima of the loss function having a given correlation with the signal. Analyzing the resulting formula we defined and located where the energy landscape becomes trivial in the sence that spurious local minima disappear. We analyzed the performance of gradient flow via integro-differential state-evolution-like equations. We numerically solved the equations and extrapolated the divergence of their convergence-time. We delimited a region of parameters for which the gradient flow is able to avoid the spurious minima and obtain a good correlation with the signal in time linear in the input size. We also analyzed the maximum-likelihood AMP algorithm, located the region of parameters in which this algorithm works, which is larger than the (numerically extrapolated) region for which the gradient flow works. We found that in cases when both the algorithms converge to an informative minima, the corresponding error is the same in both and also corresponds to the position of all the minima well correlated with the signal in the Kac-Rice approach. The relation between existence or absence of spurious local minima in the loss landscapes of a generic optimization problems and the actual performance of optimization algorithm is yet to be understood. Our analysis of the spiked matrix-tensor model brings a case-study where we were able to specify this relation quantitatively.

A Kac-Rice formula 19) at m = 0, the dashed line considers the whole non-informative band.

In the cases in which the order of the tensor p is odd we encounter an interesting phenomenon due to the different symmetries of the two types of observation. The matrix is symmetric by inverting the sign of the signal, x → -x, while the tensor is not symmetric for odd p. This creates an asymmetry in the complexity, Fig. 6 (to be compared with Fig. 2) and causes a shift toward lower correlations of the band characterizing the non-informative minima. Therefor observing when the complexity at m = 0 becomes negative does not guarantee that the non-informative minima disappeared. To do so, one must check that the whole non-informative band disappears. This should be contrasted with the case of even p where a maximum of the complexity Σ(m) is always at m = 0. These two definitions of the threshold have little, but not negligible, difference, see Fig. 7.

Observe that as ∆ p increases the peak of the complexity decreases, since the loss Eq. ( 3) tends to the simple matrix-factorization problem where the landscape is characterized by two isolated minima. This implies that the two definitions become indistinguishable for large ∆ p . In the main text we use the definition taking into account the maximum (even when it is not strictly at m = 0) because gives a more accurate characterization of the trivialization threshold. The dynamics of the gradient flow shows a dependence on the initial conditions, because formally zero correlation is a (unstable) fixed point of the GF state evolution. In practice we observe for both GF and ML-AMP that instability of the fixed point is sufficient for good performance of the algorithm. However, this makes the definition of the convergence time depend of the initial condition. We observed from our numerical solution of the GF state evolution equations that the initial condition add a factor a log m(0) to the convergence times. Thus by fitting this term and rescaling the convergence times, the different initializations collapse into a single curve, see inset of Fig. 8. Finally, the collapsed points were used to extrapolate the critical line as shown in the main text, Fig. 4.

B Gradient Flow

C AMP C.1 From AMP to ML-AMP

In this section we consider the spiked-tensor model in a Bayesian way. We show how the Bayes-optimal AMP leads to the Maximum Likelihood AMP using a temperature-like parameter T . We will introduce the algorithm AMP for a generic T , and show that as T → 0 we recover ML-AMP as presented in the main text. The probability distribution we consider is

P (X|Y, T ) ∝ e -µ x 2 i<j e -1 2T ∆ 2 Y ij - x i x j √ N 2 i 1 <•••<ip e -1 2T ∆p T i 1 ...ip - √ (p-1)! N (p-1)/2 x i 1 ...x ip 2 . (33)
The scheme for deriving AMP estimating marginals of such a probability distribution can be found in [START_REF] Thibault Lesieur | Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications[END_REF][START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF] and consist in making a Gaussian assumption on the distribution of the messages in the belief propagation (BP) algorithm and neglecting the node-dependence in the messages. A final consideration to be used in order to derive the algorithm is that the spherical constrain can be imposed by setting 1 N i ( xi 2 + σ i ) = 1 at every iteration. The resulting AMP algorithm will iterate on the following equations:

B t i = (p -1)! N (p-1)/2 k 2 <•••<kp T ik 2 ...kp T ∆ p xt k 2 . . . xt kp + 1 √ N k Y ik T ∆ 2 xt k -r t xt-1 i ( 34
) xt+1 i = 2 B t i 1 + 1 + 4 N ||B t || 2 2 , ( 35
)
σ t+1 = 2 1 + 1 + 4 N ||B t || 2 2 . ( 36
)
with || • • • || 2 2 the 2 -norm and r t the Onsager reaction term

r t = 1 ∆ 2 T 2 1 N k σ t k + p -1 ∆ p T 2 1 N k σ t k   1 N k xt k xt-1 k   p-2 . (37)
In the limit T → 0 AMP defined by Eqs. (34-37) is equivalent to ML-AMP, Eqs. [START_REF] Montanari | Graphical models concepts in compressed sensing[END_REF][START_REF] Ricci-Tersenghi | Typology of phase transitions in bayesian inference problems[END_REF](27)(28). To see this we define the rescaled variables σt .

= σ t /T , Bt . = T B t and rt . = T r t . Taking the limit T → 0 the expression for xt+1 i Eq. ( 35) and the expression for σt+1 i Eq. (36) simplify as Eq. ( 26) and as Eq. ( 26) respectively. Dropping the tildes we obtain ML-AMP as presented in the main text.

C.2 State evolution

The generic T version of AMP has a slightly more complicated SE that depends of two order parameters: the already introduced m t = 1 N i xt i x * i and q t = 1 N i (x t i ) 2 the self-overlap of the estimator. The SE equations are:

m t+1 = 2 z t (T ) 1 + 1 + 4y t (T ) , (38)
q t+1 = 4 y t (T )

1 + 1 + 4y t (T )

2

(39) and MSE t = 1 -2m t + q t , (40) with y t (T ) = z t (T )

2 + 1 T 2 q t ∆ 2 + 1 T 2
(q t ) p-1 ∆p and z t (T ) = 1

T m t ∆ 2 + 1 T (m t ) p-1
∆p . Given 1 N ||x 0 || 2 2 = 0, in the limit T → 0 AMP SE Eqs. (38-39) simplify, to a single equation corresponding to ML-AMP SE Eq. (41). This is seen by taking the limit for Eq. (39) which gives q t = 1 ∀t > 0, implying MSE t = 2(1 -m t ). Then, using the result for q t , we show that Eq. (38) tends to Eq. (30). that admits a single solution in R + :

z = 2γ ∆ 4 2 3 -γ . (44)
An important point in the phase diagram is where the algorithmic and dynamical spinodals meet, this is called the tricritical point. Its value is obtained for different p, numerically (for p > 4) and analytically (for p = 4), and is reported in Table 1. The case p = 3 does not show any tricritical point for any finite ∆ p , the two lines eventually meet at ∆ p = ∞ when the spiked matrix problem is recovered.

For the cases p > 4 we observe additionally the zero temperature analog of what is called hybrid phase in AMP in Bayes-optimal regime [START_REF] Ricci-Tersenghi | Typology of phase transitions in bayesian inference problems[END_REF]. The hybrid phase is illustrated in Fig. 10. This phase is defined as a region where the ML-AMP algorithm initialized at random converges to a solution with positive correlation but that is less correlated then the solution achievable starting from the solution. In these cases Eq. (32) does not correspond to the algorithmic spinodal but it is just the stability of the non-informative solution.

Figure 1 : 2 (

 12 Figure 1: The figure summarizes the main results of this paper for the spiked matrix-tensor model with p = 3 (left) and p = 4 (right). As a function of the tensor-noise parameter ∆ p on the x-axes, we plot the values of 1/∆ 2 above which the following happens (from above): Above ∆ triv 2 (the dashed purple line) the landscape of the problem becomes trivial in the sence that all spurious local minima disappear. Above ∆ GF 2 (the dotted blue line) and ∆ ML-AMP 2

5 Figure 2 :

 52 Figure2: The complexity Σ(m), Eq.[START_REF] Ros | Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions[END_REF], is shown for different values of parameter ∆ 2 at fixed ∆ p = 4.0 in the case p = 4. As ∆ 2 is decreased (the signal to noise ratio increases) the complexity allows to identify three main scenarios in the topology of the loss landscape. In the first case (a) only a wide band of non-negative complexity around the point of zero correlation is present, in the second case (b) minima with non-trivial correlation with the signal appear but the band around m = 0 is still present, finally (c) the signal dominates over the noise and only minima with non-trivial correlation are present. The transition from case (b) to case (c), i.e. when the support of Σ(m) ≥ 0 becomes two discontinuous points, as the bulk close to m = 0 becomes negative, is called the landscape trivialization. The ∆ 2 at which this occurs is denoted ∆ triv

Figure 3 :

 3 Figure 3: Eq. (24) characterizes the evolution of the correlation of the gradient flow with the ground truth signal, evaluated for several values of ∆ 2 , at ∆ p = 4.0 starting from m(0) = 10 -10 . The dynamics displays a fast increase of the convergence time as ∆ 2 increases. At large times, the plateau we observe has the same value of correlation m as the minima best correlated with the signal, as predicted via Kac-Rice approach.

Figure 4 : 2 = 1 .

 421 Figure 4: The convergence time the gradient flow takes to find a configuration well correlated with the signal for ∆ p = 4.0, p = 3 as a function of ∆ 2 , starting from m(0) = 10 -10 . The points are fitted with a power law consistent with a divergence point 1/∆ GF 2 = 1.35 (vertical dotted line, log-log scale of the fit shown in the inset) while landscape trivialization occurs at 1/∆ triv 2 = 1.57 (vertical dashed line).

Figure 5 :

 5 Figure 5: We show the mean-squared error (MSE) as achieved by the analyzed algorithms, for p = 3, ∆ p = 1.0 as a function of the signal-to-noise (snr) ratio 1/∆ 2 . The full cyan line corresponds to the error reached by the ML-AMP algorithm, it jumps discontinuously at 1/∆ ML-AMP 2 = 1.62. The blue points is the error reached by the gradient flow in time t < 1000. The divergence of the convergence time is extrapolated to occur at 1/∆ GF 2 = 1.97, blue dotted vertical line. The purple dotted line represents the maximum having the largest m of the complexity function Σ(m), Eq. (19). The vertical purple dashed line at 1/∆ triv

Figure 6 : 2 Figure 7 :

 627 Figure 6: Analogously to Fig. 2, the figures show the complexity, Eq. (19), as a function of the correlation with the signal for different values of parameter ∆ 2 at fixed ∆ p = 4.0 in the case p = 3.

B. 1 Figure 8 : 2 = 1 .

 1821 Figure 8: The time corresponding to convergence close to the signal is shown for ∆ p = 4.0 in the case p = 3. Different shades of grey correspond to different initial conditions, from m(0) = 10 -10 (light grey) to m(0) = 10 -42(dark grey). The different initializations collapse to a single line when the time is rescaled by a log m(0) with a = 1.3, see inset. In the figure we fit only the case m(0) = 10 -10 with a power law and use the same parameters for all the other fits with a vertical translation. The divergence point extrapolated is 1/∆ GF 2 = 1.35 and is represented by the vertical dotted line, while the dashed line identifies the landscape trivialization predicted with the Kac-Rice formula, 1/∆ triv 2 = 1.57.

Figure 10 :

 10 Figure 10: Phase diagram as shown in Fig. 9 for the case p = 6. The difference between p = 3, 4 and p > 4 is that a new phase appears, called hybrid hard phase, where two fixed points of ML-AMP aligned with the signal are present and the convergence to one or the other depends on the initialization. The region is highlighted in the inset. In the phase diagram the grey dashed line represent the threshold above which the non-informative fixed point becomes unstable.

  p

Table 1 :

 1 Table of the values of tricritical points for p ≥ 4.

		∆ 2	∆ p
	4 2 3	0.667 4 3	1.333
	5	0.470	0.451
	6	0.384	0.305
	7	0.322	0.220
	8	0.279	0.172
	9	0.246	0.147
	10	0.220	0.121

Acknowledgments

We thank G. Ben Arous, G. Biroli, C. Cammarota, G. Folena, and V. Ros for precious discussions. We acknowledge funding from the ERC under the European Unions Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe; from the French National Research Agency (ANR) grant PAIL; and from "Investissements dAvenir" LabEx PALM (ANR-10-LABX-0039-PALM) (SaMURai and StatPhysDisSys). This research was supported in part by the National Science Foundation under Grant No. PHY-1748958.

 1with two additional lines. The dashed brown line is the limit predicted by Kac-Rice formula where the support of the Σ(m) ≥ 0 becomes disconnected (above the line). The full orange line is related to the ML-AMP algorithm, is called dynamical spinodal, below it the algorithm converenges to m = 0 even if initialized in the solution. In the insets we show the large ∆ p behaviour of the thresholds, where we can observe that the lines merge at infinity.

C.3 Derivation of spinodals

From SE Eq. ( 30) we can obtain analytical equations for the spinodals, the threshold of stability of the different ML-AMP fixed points. We have xt+1 = f SE (z t ) with

. We can now define either

As remarked in [START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference[END_REF], the spinodals are given by the following condition:

A trivial solution is given by z → 0 corresponding to stability of the non-informative solution m = 0, and gives the algorithmic spinodal for the cases p ∈ {3, 4}. This solution and has a very simple equation for every p:

√ γ giving Eq. (32), already presented in the main text. An interesting implication of Eq. ( 32) is that it is independent from the value of p, it is in some sense universal among the 2 + p-models.

The expression for the stability of the informative solution, dynamical spinodal, is less straightforward, but analytical progresses can be done in the cases p = 3 and p = 6 (using Cardano formula) and in the case p = 4 for which it is equivalent to a second order polynomial