Complex energy landscapes in spiked-tensor and simple glassy models: ruggedness, arrangements of local minima and phase transitions

We study rough high-dimensional landscapes in which an increasingly stronger preference for a given configuration emerges. Such energy landscapes arise in glass physics and inference. In particular we focus on random Gaussian functions, and on the spiked-tensor model and generalizations. We thoroughly analyze the statistical properties of the corresponding landscapes and characterize the associated geometrical phase transitions. In order to perform our study, we develop a framework based on the Kac-Rice method that allows to compute the complexity of the landscape, i.e. the logarithm of the typical number of stationary points and their Hessian. This approach generalizes the one used to compute rigorously the annealed complexity of mean-field glass models. We discuss its advantages with respect to previous frameworks, in particular the thermodynamical replica method which is shown to lead to partially incorrect predictions.


I. INTRODUCTION
Characterizing rough multi-dimensional energy landscapes is a challenging task that is central in many different fields from physics to computer science, highdimensional statistics, machine learning and biology. In a nutshell this problem consists in analyzing the statistical properties of functions defined on very high dimensional spaces. Relevant information that one wants to obtain is for instance the number of minima at a given energy, and more generally of the critical points, and the spectral properties of their corresponding Hessian. This issue is crucial to understand the dynamics within these landscapes, in particular gradient descent which has many physical and practical applications. Depending on the context, the landscape can correspond to the energy of a physical system, to the loss-function of a machine learning algorithm, to the cost function of an optimization problem or to the fitness function of a biological system. Pioneering works on this subject were done in physics, in the context of mean field spin-glasses, starting from the 80s [1][2][3][4], see [5] for a review. One of the essential results, besides the explicit computations in several models, was the understanding that the statistical properties of rough energy landscapes are the ones characteristic of two different physical systems: spin-glasses and glasses [5]. The origin of this universality lies in replica theory: the properties of the landscape are actually encoded in the type of mean-field solution obtained by the replica method, respectively full replica symmetry breaking and one step replica symmetry breaking [6]. Remarkably, it was also realized that pure systems can behave as disordered ones, as first found in long-range spin models [7] ; accordingly energy functions of several complex systems are qualitatively similar to random functions. In mathematics, in particular in probability theory, there has been a recent and growing research activity aimed at developing rigorous analysis of rough energy landscapes. Starting from the seminal work [8] the Kac-Rice method has emerged as the mathematical framework suited to do that [9][10][11][12][13]. It allowed to put on a firmer basis previous results obtained in the physics literature, and it highlighted important relationships with random matrix theory. Moreover, it has been recently exploited to analyze landscape properties of machine learning and inference models [14,15]. The recent results and questions concerning the statistical properties of rough landscapes make clear that what found for mean-field glassy systems represents only a facet of a much more general challenge. There are several different directions in which further investigations are timely and interesting. One of them is the characterization of landscapes in current problems central in machine learning and high-dimensional statistics, such as the analysis of rough energy landscapes and associated phase transitions when an increasingly stronger preference for a given configuration arises. This problem is central in data science (the signal versus noise problem) [16], as well as in biology and in physics, in cases where a arXiv:1804.02686v2 [cond-mat.dis-nn] 24 Apr 2018 specific ground state competes with many random ones (e.g. protein folding [17] and random pinning glass transition [18]). Another important and quite distinct research direction consists in studying the number of equilibria in non-conservative dynamical systems that arise in neuroscience [19] and theoretical ecology [20]. In this case, forces do not derive from a potential, hence there is no landscape to start with, but nevertheless information about the number of equilibria and their stability can be obtained by methods similar to the one used for the conservative case [21,22]. From the methodological point of view, the main open crucial issue is developing the Kac-Rice method to compute the typical number of critical points, related to the average of the logarithm of the number of critical points (called quenched entropy). Computing the logarithm of the average (called annealed entropy), as done until now, is correct in a few cases only [13]; in general, the two computations lead to different results even at leading order. Physics methods based on replica theory and supersymmetry provided guidance and results in specific cases, but as we shall discuss in the following they suffer important limitations [23][24][25][26][27][28][29]. Our work has a double valence. One is conceptual: we present a general analysis of the properties and the phase transitions occurring in rough energy landscapes whenever an increasingly stronger preference for a given configuration arises, an interesting and timely issue as discussed above. The other is methodological: we develop the sought generalization of the Kac-Rice method to compute the typical number of critical points and the corresponding quenched entropy, a theoretical framework expected to have multiple applications in several fields. Overall, our work opens the way to throughout analysis of the statistical properties of rough landscapes in topical problems relevant in several different fields, from physics to machine learning and biology. We focus on the p-spin spherical model and add to its Hamiltonian a term favoring all configurations that are close to a given one [30]. This choice is natural from different points of view. First, the system without the additional extra term has already proven to be an instrumental paradigm for rough energy landscapes [5,31], so it is a natural starting point to study the effect of a preferred configuration on a random landscape. Second, it is directly relevant for very recent problems studied in the computer science literature; in fact a particular realization of it corresponds to the so called spikedtensor model, which recently attracted a lot of attention [15,[32][33][34][35]. The thermodynamics of the system we focus on, that we henceforth call generalized spiked-tensor, has been originally introduced in Ref. 30 to study the effect of a ferromagnetic coupling on a p-spin spherical model. Here we investigate in detail its energy landscape. Depending on the functional form of the additional term, we generically find different scenarios and different types of energy landscape (or geometric) phase transitions. Although this model is certainly extremely simplified, we think that the lessons that can be learnt from its analysis provide instrumental guidelines and extend to more realistic cases. Moreover, because of its relation with the spiked-tensor model, our results are directly relevant to current issues investigated in high-dimensional statistics and inference. As stressed above, one of the main outcome of our work is the construction of a general Kac-Rice method which allows one to analyze cases in which the so-called quenched entropy does not coincide with the annealed one, as it happens for the model we consider. Since we use replicas in a rather innocuous way-we remain at the replica symmetric level-transforming it from a theoretical physics technique to a fully rigorous one should be within reach in a not too distant future. In the following two sections we present a summary of the main results. In Section IV we discuss the zero temperature thermodynamics of the model by means of the replica method. In Section V we present the new Kac-Rice method for the quenched complexity, and we compare its findings with the ones obtained with the replica method in Section VI. After reviewing the implications of these findings for the special case of the spiked-tensor model in Section VII, we present our conclusions in Section VIII.

II. DEFINITION OF THE MODEL
We consider the Hamiltonian or energy functional: where the first sum is over all distinct p-uples and the subindices run from 1 to N . The configuration space of the model is the sphere of radius √ N , i.e. a given configuration is a vector s of N components {s 1 , s 2 , . . . , s N } such that N i s 2 i = N . The N -dimensional vector v 0 points towards a specific direction, say v 0 = {1, 1, . . . , 1} without loss of generality (we have imposed on v 0 the same normalisation condition as s). In the following we are going to refer to this preferential direction of the model as the North Pole. The first term of H p,k is the Hamiltonian of the standard spherical p-spin model [5] with random coupling J i1,i2,...,ip normally distributed with zero mean and variance J 2 = p!/2N p−1 . The second term represents an energetic gain when the system's configuration s is aligned with v 0 . We generically describe this energetic gain by a function f k (x) of the scalar product x = N i s i v i 0 /N . Our aim is to use f k (x) as a template of a smooth function defined on the N dimensional sphere, with a deep minimum in a specific direction. We found that the main relevant features of f k (x) are its derivatives in x = 0: the sub-index k indicates what is the first non zero derivative in x = 0. We assume that the function f k (x) reaches its highest value FIG. 1. Evolution of the energy landscape in case I. In this drawing we illustrate the evolution of the energy landscape due to the increase of r in Case I (k = 1). The red strip denotes the region on the sphere where minima lie in an exponential number. The continuous yellow line corresponds to the parallel where the deepest minima are located. The dashed yellow line corresponds to the parallel where the most numerous minima are located. At rc the energy landscape has a transition: for r < rc it is rough and full of minima, for r > rc it is smooth and only contains one minimum (represented by the yellow dot in the figure).

FIG. 2.
Evolution of the energy landscape in case II. In this drawing we illustrate the evolution of the energy landscape due to the increase of r in Case II (k = 2). The red strip denotes the region on the sphere where minima lie in an exponential number. The continuous (dashed) yellow line corresponds to the parallel where the deepest (most numerous) minima are located. The energy landscape has several transitions. At r2ND the deepest minima are no longer on the equator and move toward the poles. Afterwards the band containing the exponential number of minima fractures in three parts, one around the equator and two symmetric ones closer to the poles. At rc the bands closer to the pole implode and are replaced by two isolated global minima (the one on the south hemisphere is not visible since it is on the back of the sphere) but the band at the equator persists. Finally, for even larger values of r, the landscape becomes completely smooth with only two symmetric minima. in x = 1, is zero for x = 0 and is monotonously increasing in [0,1]. It also has a symmetric or an antisymmetric continuation for x ≤ 0 depending whether k is even or odd respectively. For concreteness, we shall often refer to the case f k (x) = x k /k, which was first introduced in Ref. 30. When p = k, the model corresponds to the so called spikedtensor model which has been the focus of several recent studies in the computer science literature [15,[32][33][34][35]. In particular, for this limiting case the calculation of the average number of stationary points has been very recently performed in [15].

III. SUMMARY OF RESULTS
The energy function H p,k contains two terms. The first is a random Gaussian function, whereas the second one is deterministic. These contributions are competing: the random fluctuations encoded in the former lead to an exponential (in N ) number of critical points. On the sphere in very high dimensions the majority of the con-figurations are orthogonal to the North Pole, thus it is on the equator that we expect the deepest minima created by the first term alone. Since there are exponentially less configurations in the direction of v 0 , and the less so when the overlap with v 0 is higher, the random fluctuations alone lead to minima of higher energy on parallels closer to the north pole. On the other hand, the deterministic term energetically favors configurations aligned with v 0 . In consequence, depending on the relative strength of the two terms, that can be tuned by changing the value of r, and on the form of the function f k (x) the resulting rough energy landscape changes shape and the low lying energy minima change position and nature from many to a single one. As we shall see, all that corresponds to phase transitions in the geometry of the landscape. Topological phase transitions, occurring when the landscape changes from being complex to simple, have been recently studied in [47,49,50] and dubbed topological trivialization. The change in the global minima structure, which is directly accessible to a thermodynamic study, was already reported in Ref. 30. In the following we present our main results on the evolution with r of the full energy landscape. For the sake of the presentation we group the different scenarios in three classes, associated with the behavior of the global minima as a function of r.
This case corresponds to functions f k (x) which are monotonically increasing and such that f (0) > 0. The simplest example, f k (x) = x, corresponds to the p-spin spherical model in an external magnetic field (with r playing the role of the field), and has been studied in [3,48,49]. In agreement with those analyses, we find that the energy landscape evolves as illustrated in Fig. 1. First, at r = 0, there are an exponential number of minima located around the equator, i.e. for q ∈ [q m (0), q M (0)], where q m (0) = −q M (0). This corresponds to the first sphere on the left, in which the presence of minima is indicated with a red strip. The deepest minima, not exponentially numerous, are at q = 0 and correspond to the continuous yellow line. The most numerous states, which are also the marginally stable ones since the density of states of their Hessian is a Wigner semicircle with left edge touching zero, are also at q = 0 for r = 0 (and they are of course at higher energy). By increasing r the strip containing all the minima moves toward the north pole, see the second sphere from the left in Fig. 1. The deepest ones are on a parallel closer to the north pole as soon as r > 0. The most numerous ones, always marginally stable, are now on a different parallel with smaller latitude, as it can be expected on general grounds since in order to have a lot of minima it is better to avoid too large latitudes at which less configurations are available (they are represented by a yellow dashed line in the figure).
By increasing r the landscape becomes smoother due to a larger deterministic term and, accordingly, the number of minima and the strip where they are located shrink until reaching a value r c above which only one minimum remains. This corresponds to a phase transition of the landscape, which is associated to recovering a replica symmetric solution for the global minimum within the replica method, and hence also to a phase transition in the thermodynamics (related to the structure of the global minima). For r > r c there is only one minimum in the energy landscape. In this case the random contribution due to the first term in the Hamiltonian is no longer strong enough to create a rugged landscape but still deforms it sufficiently to move the global minimum at a finite overlap with v 0 . This corresponds to the rightmost sphere in Fig. 1. As we shall see in the following, a much richer energy landscape evolution is found for k > 1. In these cases the behavior of the global minima is only a facet of a more general complex organization in configurations space. This regime corresponds to functions f k (x) which have vanishing derivative in x = 0 but finite second derivative and are monotonically increasing from x = 0 to x = 1. In order to simplify the discussion we consider the symmetric case in which f k (−x) = f k (x). The simplest example of such a function is f k (x) = x 2 /2. With this choice, H p,k corresponds to a p-spin spherical model with an extra ferromagnetic interaction among spins (r plays the role of the coupling). The evolution of the energy landscape is now different from Case I and it is illustrated in Fig. 2. The starting point at r = 0 is the same. However, by increasing r the strip containing all the minima widens and the deepest ones and the most numerous ones (always marginally stable) remain stuck on the equator. Actually they are exactly the same ones found for r = 0 since f k (x) has no effect on the equation that determines the critical points on the equator (this is due to the vanishing of the first derivative in x = 0). This situation persists until r = r 2ND , at which a secondorder phase transition takes place at the bottom of the landscape, as already found in Ref. 30. By increasing r above r 2ND the deepest minima continuously detach from the equator, see the second sphere in Fig. 2 (due to the symmetry x → −x they are located both in the north and south hemispheres). The behavior for larger r is different from case I: there is first a transition in the structure of the energy landscape in which the strip separates in three bands, two closer to the north and south poles respectively, to which the deepest minima belong, and one around the equator where the most numerous ones are located, see the middle sphere in Fig. 2. At r = r c there is another transition at which the two bands closer to the north and south poles FIG. 3. Evolution of the energy landscape in case III (option A). In this drawing we illustrate one of the possible evolutions of the energy landscape due to the increase of r in Case III. The red strips denote the regions on the sphere where minima lie in an exponential number. The continuous yellow lines corresponds to the parallel where the global minimum is located. At rc an isolated local minimum appears. The dotted yellow line denotes that it is not yet the global one. At r1ST > rc the deepest minimum is no longer on the equator and switches discontinuously to the isolated one close to the north-pole. For larger values of r the global minimum approaches the north pole and the band around the equator shrinks but does not disappear for any finite r. The most numerous states, denoted by a dashed line, are always located on the equator.
FIG. 4. Evolution of the energy landscape in case III (option B). In this drawing we illustrate one of the possible evolutions of the energy landscape due to the increase of r in Case III. The red strips denote the regions on the sphere where minima lie in an exponential number. At r1ST the deepest minimum is no longer on the equator and switches discontinuously to one at higher latitude. For larger values of r the band of local minima divides in two: one around the equator and one around the global minimum. For r = rc > r1ST the latter band shrinks to zero and the global minimum becomes isolated. The band around the equator shrinks but does not disappear for any finite r. The most numerous states, denoted by a dashed line, are always located on the equator.
containing an exponential number of minima shrink to zero and are replaced by an isolated global minimum per hemisphere (fourth sphere from the left in Fig. 2). This corresponds to recovering the RS solution in the thermodynamic treatment [30]. Finally, at even larger r all minima around the equator disappear and a final transition toward a fully smooth landscape characterized by only two minima takes place. This corresponds to the rightmost sphere in Fig. 2. The most numerous minima remain always at the equator for any value of r until this final transition at which they disappear. However, they change nature when increasing r: at the beginning all the eigenvalues of their Hessian are distributed along a Wigner semi-circle whose left edge touches zero (so-called threshold states), whereas at large values of r they are all distributed along a Wigner semi-circle whose support is strictly positive except for one eigenvalue, corresponding to an eigenvector oriented toward the north pole, which pops out from the semi-circle and is located exactly in zero. Thus, in both cases they are marginally stable but in a very different way. In conclusion, in the k = 2 case in which the strength of the deterministic part is weaker in particular around the equator, the spurious local minima created by the random fluctuations are more stable. This results in a different evolution of the landscape, that before becoming fully smooth is characterized by isolated islands of ruggedness around the equator and close to the global minima.
This regime corresponds to functions f k (x) which are monotonically increasing in [0, 1] and have vanishing first and second derivatives in x = 0. For simplicity, we shall consider even and odd functions under x → −x when k is odd and even respectively. The simplest example of such a function is f k (x) = x k /k with k ≥ 3, first introduced in Ref. 30. With this choice and taking p = k, H p,k corresponds to the spiked-tensor model recently investigated in Refs. 15, 32-35. The particularity of Case III is that the critical points on the equator are not affected at all by the deterministic perturbation, not even their Hessian (contrary to case II) since f (0) = 0: they remain stable and unperturbed for any finite value of r. In consequence, there is always a strip of minima around the equator. We have found that different evolution are possible in Case III depending on p, k.

Option A
This is the case found for example for spiked-tensor models such as p = k = 3 and p = k = 4. For concreteness we focus on p = k = 3 (p = k = 4 is analogous but one has to take into account that f k (x) is even instead of being odd). A band of minima, growing with r, is found around the equator. At a value r c an isolated minimum detaches from the top of the band, and for larger values of r it moves to higher latitudes, while the rest of the band shrinks around the equator. The deepest minima are located on the equator and are the ones of the original (unperturbed) p-spin model until a value of r, that we call r 1ST , is reached. When r reaches the value r 1ST the global minimum switches from the equator to the single minimum outside the band and close to the north pole. Increasing r further the isolated global minimum approaches the north pole and the band around the equator shrinks but never disappears for any finite r. The most numerous states are on the equator and are the threshold states of the unperturbed p-spin model. The evolution of the energy landscape and its transitions are illustrated in Fig. 3.

Option B
This is the case found for example for p = 3 and k = 4. A band of minima, which first grows with r, is found around the equator. The deepest minima are located on the equator until r 1ST and are the ones of the original (unperturbed) p-spin model. When r reaches the value r 1ST the global minimum switches discontinuously from the equator to another minimum inside the band, at higher latitude. Increasing r further, the band divides in two: one closer to the equator and one around the global minimum. For r = r c , the band around the global minimum shrinks to zero (this corresponds to recovering the RS solution in the thermodynamics treatment [30]). For r > r c the global minimum is isolated. The remaining band around the equator shrinks but never disappears for any finite r. The most numerous states are on the equator and are the threshold states of the unperturbed p-spin model. The evolution of the energy landscape and its transitions are illustrated in Fig. 4.
Two other options are possible: the discontinuous transition at r 1ST could take place after that the band has divided and, depending whether r 1ST is larger or smaller than r c , it could take place when the global minimum is isolated (option C) or is still surrounded by many other local minima (option D). We did scan a few more (see Fig. 5), but not all possible values of p, k, nor analyzed all possible functions f k (x) to search for these two behaviors but this can be easily (even though painfully) done if specific interest in these intermediate cases arises.

D. Randomness versus deterministic contribution
A short conclusion of the results presented above is that the evolution of an energy landscape in which random fluctuations compete with a deterministic contribution favoring a single minimum depends on the behavior of the deterministic part on the portion of configuration space where the majority of minima created by randomness lie. If the deterministic part affects and deforms these minima then the evolution is quite simple: the number of minima decreases and they become more and more oriented toward the direction favoured by the deterministic part until a point at which only one isolated global minimum remains. A different behavior is instead found when the deterministic part does not deform the majority of minima created by randomness. In this case, the competition between random and deterministic contributions is resolved in two different ways: it deforms the landscape in the proximity of the configurations favoured by the deterministic part, which can even result in an island of ruggedness and many local minima, and it creates a very rugged landscape in the region where the deterministic part has no effect, where the majority of the configuration lie. As it can be easily guessed, this landscape structure can have crucial consequences on dynamical properties. We shall discuss these issues and, more generally the implications and consequences of our results in the Conclusion. In the following we present the methods we used, and our findings in more details. We first recall the ther-modynamic analysis of the model, focusing on the zero temperature limit, in order to discuss the behavior of the global minima of the landscape. Subsequently, we analyze the evolution of the full set of minima, encoded in the quenched complexity. The starting point of the thermodynamic analysis is the evaluation of the free energy f (β, r), obtained by computing the n-times replicated partition function Z n : and where the signal contribution to the Hamiltonian is represented by f k (x) = x k /k. To gain direct information on the energy landscape we focus on the zero temperature limit, when the equilibrium states dominating the partition function (2) coincide with the absolute minima (or minimum) of the energy landscape. This thermodynamic analysis gives then access to the equilibrium transitions, which occur whenever these global minima detach from the equator and move at higher latitudes in the sphere, becoming correlated to the signal. Moreover, it allows to determine whether the bottom of the energy landscape is simple, i.e. just one global minimum, or has a more complicated structure, encoded in the Replica Symmetry Breaking (RSB) formalism.
A. Energy at fixed overlap and the cases I,II,III We first describe the main results of the replica analysis, the computation is shown later. One important remark is that the signal affects the model's solution only through the value of the typical overlap q with the north-pole. To get the zero temperature solution, it is interesting then to focus on the intensive ground state energy of the original p-spin spherical model, i.e. without the function f k (x), for configurations constrained to have a fixed overlap q. We denote this function E(q). As it is expected by the q → −q symmetry of the original p-spin problem, E(q) = E GS + Cp 2 q 2 +O(q 4 ) for small q, where E GS is the intensive ground state energy of the p-spin spherical model, and C p happens to be a positive constant. This result already allows us to show the existence of the three regimes discussed in the previous section because now we can obtain and study the ground state energy of our model as E(q) − r q k k . • Case I: If k = 1 then, no matter how small is r, the ground state is at q > 0 and increases when r is augmented. This is the first scenario described in the previous section.
• Case II: If k = 2 then the ground state is at q = 0 for r < r 2ND = C p and becomes continuously different from zero by increasing r above C p . This corresponds to a second-order like transition and to the second scenario discussed before.
• Case III: if k ≥ 3 then a discontinuous transition is bound to take place: for r < r 1ST the ground state is at q = 0, whereas for r > r 1ST it jumps to a finite value. This corresponds to a first-order like transition and to the third scenario discussed before.
An insight on the changes in the structure of the bottom of the landscape can be obtained along the same lines. At r = 0 the replica solution is 1RSB. Proceeding as before, i.e. studying the p-spin spherical model at fixed q, one can show that at fixed p the solution always remains 1RSB until a given value of q c is reached where the 1RSB-RS transition takes place. Moreover the replica structure is the same for identical values of q.
Only the way in which q changes as a function of r depends on the value of k. In consequence, when the ground state is at q = 0, there is a 1RSB structure of the energy landscape close to the global minimum (roughly speaking the energy landscape is rough close to the bottom). When the ground state is at q > 0, for r larger than a critical value r c , the structure of the energy landscape close to the global minimum become RS (roughly speaking the energy landscape is convex close to the bottom). In case III there are two minima of E(q) close to the first order transition: one at q = 0 and one at q > 0. The high-overlap minimum can become 1RSB before or after the discontinuous transition depending on the value of r c and r 1ST . The replica analysis that we present below, see also Refs. 30 and 36, allows to find the models in which this happens, see Fig. 5. The yellow sheet identifies (on its right) models that display a regime in which a rough landscape around the high-overlap global minimum is present for r 1ST ≤ r ≤ r c .
We present below the replica computation. The following section will be also useful to fully understand the The standard replica computation [5] for f (β, r) leads to the following result and Q αβ is an (n + 1)x(n + 1) matrix (α, β ∈ [0, n]) composed by 1 on the diagonal, Q 0,a on the 2n entries of the first line and column, and a matrix Q a,b with a, b ∈ [1, n] on the remaining nxn block. The action has 3 parts: the energy of the p-spin part of the original Hamiltonian, the energy due to the added potential f k controlled by the parameter r, and the entropy of a N -dimensional spherical system with one special direction.
To proceed in the calculation, we use a RS ansatz on the entries Q 0,a , Q 0,a = q, and the usual RS or 1RSB ansatz for the Q a,b matrix (no additional breaking of replica symmetry is expected). The first case corresponds to Q a,b = q if a = b. In the second case the replicas are classified according to n/m different blocks, Q a,b = q 1 for a = b with a and b in the same block of size m, and Q a,b = q 0 when a and b belong to different blocks. The 1RSB ansatz contains the RS one: the second can be recovered by setting either m = 1 or q 1 = q 0 . We thus only focus on the first.

The 1RSB saddle point equations
The expression of the 1RSB action S 1RSB in the N → ∞ and n → 0 limit is reported in Appendix IX A for generic values of β. When β → ∞ the action reads where the parameters β(1 − q 1 ), q 0 , βm, and q have to be determined by the following saddle point equations: with f k (x) = x k−1 being the first derivative of f k . For each value of r, the value of q obtained solving the saddle point equations gives the latitude of the deepest minima of the landscape, while the function −S 1RSB /β evaluated at the saddle point parameters gives their energy density.

RS-1RSB transition for the high overlap phase
We now discuss the limit where the ground state energy ceases to be obtained by a 1RSB solution and is instead determined by a RS one. The transition between these two regimes signals a change in the structure at the bottom of the landscape from many low-lying minima to one single global minimum. We call r c the corresponding critical value of r at which this occurs. The first piece of information about this change of structure is obtained by expanding the four 1RSB saddle point equations [31] for small q 1 −q 0 , and by keeping the lowest order non-zero terms. This gives four equations, see Appendix IX A. Applying them to the 1RSB solution with high q, we get that the critical point occurs at when At this point the high q solution recovers a RS structure, i.e. it becomes a single minimum. Still one has to consider whether this solution is a global minimum of the energy landscape or only a local one. This piece of information is recovered by comparing the energy cost of the high q solution with the solution with q = 0, when this does still exist. A full account of all the possible models' solution obtained by using all the gathered information is presented in the next section.

C. Results
From the numerical study of all the T = 0 equations above we recover the three distinct scenarios accounted for in Sec. III. Case k = 3 and higher. For k > 2, we find a stable 1RSB q = 0 solution at every value of r. This solution is orthogonal to the signal and completely dominated by the noise represented by the p-spin part. Beside this solution, when r increases we find a second, high-q solution which undergoes a continuous transition between a 1RSB phase and a RS phase at r c . The high-q solution (q = 0) contains at least partial information about the signal, the amount of this information being represented by the overlap q. This solution is at first metastable compared to the q = 0 state, but it becomes stable at higher r. This occurs through a first order transition at r 1ST . If r 1ST > r c the first order-transition marks a thermodynamic discontinuity between a 1RSB state (at q = 0) and a RS state (the high q one). This scenario is generally found for k ≥ p as shown in the phase diagram in Fig. 5. If r 1ST < r c instead two transitions are observed when r increases. A first order transition will occur at lower r showing the exchange of stability between the q = 0 and high-q, 1RSB, states. A continuos transition between the 1RSB and the RS phase will follow within the high-q state at higher r. An intermediate complex phase, related to a rugged landscape, already containing partial information on the signal, or North Pole, then emerges in this case.
Case k = 2. The case k = 2 is qualitatively different from k = 3. The first order transition is replaced by a continuous, 2nd order-like, transition between the q = 0 state and the high-q state before the last one becomes RS at r c . As explained before this can be rationalised thinking that the 1RSB action is quadratic in q. As such a term f k with higher power of q (k > 2) cannot affect the local stability of the q = 0 state. When k = 2 instead, f k can counterbalance the quadratic contribution of the 1RSB action leading to the instability of the q = 0 solution at high enough r 2ND . The 1RSB-RS transition happens for a strictly larger value of r since it takes place for a finite value of q. Case k = 1. Finally the case k = 1 has been extensively studied years ago [31], it corresponds to the p-spin spherical model in an external magnetic field. In this case there are no competing 1RSB states at all. The linear field immediately shifts q of the 1RSB phase away from zero until the continuous transition at r c brings the 1RSB phase into the RS solution.
In order to show concrete examples, we report in Table I the different transition values for different p and k, in particular for the spike-tensor model p = k = 3. As shown above, in the k ≥ 3 case whether the discontinuous transition to the high q phase takes place before or after the 1RSB-RS transition depends on the model, i.e. on p and k. In order to find a general criterion we evaluate the action of the high overlap phase at the 1RSB-RS transition: see Eqs. (61,62,63,64) in Appendix IX A. We then compare this action to the one of the 1RSB phase with q = 0: S L 1RSB . There are two possible cases: In this case the high-overlap phase becomes energetically favorable after the 1RSB-RS transition takes place. Thus, the rough energy landscape around the north pole described by the 1RSB phase does not contain the global minima of the landscape, but only some local (metastable) ones.
• S Hc RS > S L 1RSB . In this case the high-overlap phase becomes energetically favorable before the 1RSB-RS transition takes place, hence there is a range of r where the stable high-overlap phase is 1RSB. This region extends up to r c .
In Fig. 5 we show a diagram in the p, k−space, with the black line representing the point where S Hc RS /β = S L 1RSB /β at zero temperature. To the right (respectively left) of the line lie models in which the 1RSB-RS transition takes place after (respectively before) the discontinuous transition to the high overlap phase. Whenever the 1RSB-RS transition takes place after the discontinuous transition, the complex phase with partial information about the signal contains ground state minima for a finite range of r. The height of the coloured sheet in Fig. 5 represents the range of r for which this holds, i.e. r c − r 1ST . Note that the range becomes larger and larger when p increases for fixed k, or k decreases at fixed p. The two complex phases we have discussed present a multi-minima structure that is worth studying to get insights on the possibility to recover the signal through different sampling dynamics. We perform this study in the following section, making use of the Kac-Rice formalism. A comparison with the results obtained by means of the replica formalism is postponed to Sec. VI.

V. LANDSCAPE ANALYSIS VIA REPLICATED KAC-RICE FORMULA
In this section, we present the analysis of the energy landscape of H p,k (r) performed through the replicated version of the Kac-Rice method.
Our aim is to determine the number N N ( , q) of local minima (or, more generally, of stationary points) of the energy functional, having a given energy density and a fixed overlap s · v 0 = N q with the special direction v 0 . The number N N ( , q) is a random variable that, when the random fluctuations dominate over the signal, scales exponentially with N . This occurs over a finite range of energies; among the exponentially-many local minima, the lowest-energy ones dominate the thermodynamics of the model (described in detail in the previous section), while the higher-energy ones are expected to play a relevant role when discussing the dynamical evolution on the energy landscape. We are interested in determining the exponential scaling of the typical value of N N ( , q), that is, we aim at computing the quenched complexity Σ p,k ( , q; r) defined as As anticipated, we perform the calculation making use of the Kac-Rice formula [51,52]. This formalism has been recently exploited to characterize the topological properties of random landscapes associated to the pure and mixed p-spin models [10,11], to the spiked-tensor model [15], as well as to count the equilibria of dynamical systems modeling large ecosystems [22,47] and neural networks [21]. In these contexts, results have been given for annealed complexity, which governs the exponential scaling of the average number of stationary points, or equilibria. This corresponds to averaging N N over the disorder realization before taking the logarithm, at variance with Eq. (6).
For the Hamiltonian H p,k (r) with r = 0, it is known that the quenched and annealed prescriptions give the same result for the complexity [3,13]. In presence of a signal, however, this equivalence does not hold (as we show below), so that the quenched calculation becomes necessary. We perform the latter by means of the replica trick, via the identity analytically continuing the expression for the higher moments of N N . The replicated version of the Kac-Rice formula allows us to obtain (to leading order in N ) the moments N n N , for integers values of n. As we show in the following, the expression for N n N that we obtain involves n critical points s a , a = 1, · · · , n, each with energy density and overlap N q with the North Pole. Introducing their mutual overlaps s a · s b = N q ab , we find that we can parametrize the moments as: where the integral can be computed with the saddle point approximation, optimizing over the order parameters q ab . In consequence the action evaluated at the saddled point directly gives ln N n N ( , q) /N up to vanishing corrections in the large N limit. To get the complexity, i.e the typical value of the number of critical points, we perform this calculation assuming replica symmetry, meaning that we set q ab ≡ q for a = b, and take the n → 0 limit; we expect this to give accurate results, in view of the fact that H p,k (r) does not exhibit full-RSB but only 1-RSB in the statics. Before entering into the details of the calculation, we collect the main resulting expressions in the following subsection.
A. The main results: quenched complexity and mapping between the three cases For arbitrary values of n and assuming replica symmetry, we find that the action in Eq. (8) is given by where I(y) is an even function of its argument, equal to: The dependence on the overlap q between the various replicas enters in the termQ (n) p,k , which reads: where we defined From this result, we can readily obtain the expression for the annealed complexity, which is obtained setting n = 1. In this case, the dependence on q drops (as it is natural to expect, since there is only one replica and thus no overlap with any other one), and the action reduces to: For p = k, this expression agrees with the results in Ref. [15]. The annealed complexity is an upper bound to the quenched one. As we argue in Sec. V I, it captures correctly the properties of the energy landscape whenever this is smooth and has only one isolated minimum, i.e., in the regime r > r c . Our result at fixed n provides all the integer moments of the number of critical points. To get the quenched complexity, the limit n → 0 has to be performed, by analytically continuing (9). The result is whereQ p,k is defined fromQ (n) p.k = nQ p,k + O(n 2 ) and reads while q SP = q SP ( , q) is the saddle point extremizing the function (14). The evaluation of the quenched complexity therefore requires to compute a saddle point on q for given values of the parameters q, . A substantial simplification comes from a general identity that we derive in Sec. V H and which relates, for fixed p and q, the complexities Σ p,k ( , q; r) for different values of k: for f k ≡ f k (q), meaning that all complexity curves for k > 1 can be derived from the ones at k = 1. This is convenient, as it allows us to solve the saddle point equations for q in one single case. We remark however that not all the properties of the landscape at k > 1 can be deduced from the case k = 1: in particular, the analysis of the stability of the stationary points (i.e., of the spectrum of their hessian) has to be performed separately for any k, as we discuss in more detail in the following subsections.

B. The replicated Kac-Rice formula
Here we present the replicated version of the Kac-Rice formula, and outline the main steps of the subsequent calculation. For convenience, we introduce the vectors σ = s/ √ N and w 0 = v 0 / √ N having unit norm, and we define the rescaled energy functional . . σ ip denoting the p-spin energy functional with rescaled coupling satisfying (J i ) 2 = p!. We count the stationary point σ of this functional satisfying h [σ] = √ 2N and σ·w 0 = q, which are in one-to-one correspondence with the stationary points of H p,k (r) with energy density and s · v 0 = N q. The Kac-Rice formula incorporates the spherical constraint, as it counts the number of stationary points of the functional h[σ] restricted to the unit sphere; such points σ nullify the surface gradient of (16), which is a vector g[σ] lying on the tangent plane to the sphere at the point σ. Similarly, their stability is governed by the Hessian on the sphere, which we denote with H[σ] (see Eq. (19) for a precise definition of this matrix). Given n replicas σ a , a = 1, · · · , n, we introduce the shorthand notation , and denote with p σ the joint density function of the (N − 1)n gradients components g a α and of the n functionals h a , induced by the distribution of the couplings J i in (16). With this notation, the replicated Kac-Rice formula reads: with In (17) the integral is over n replicas σ a constrained to be in the unit sphere, at overlap q with the vector w 0 . The function p σ (0, ) is the joint density of gradients and energies evaluated at g a α = 0 and h a = √ 2N for any a = 1, · · · , n. The expectation value (18) is over the joint distribution of the Hessians H a , conditioned on each σ a being a stationary point with rescaled energy √ 2N , and overlap q with w 0 . The computation of the moments (17) requires to determine, for each configuration of the replicas σ a , the joint distribution of the variables H a αβ , g a α and h a , which are all mutually correlated and whose distribution depends, in principle, on the coordinates of all the replicas. For the simplest case (n = 1) of a single replica σ, it can be shown (see the discussion below, and Refs. [8,10]) that (i) the gradient g [σ] is statistically independent from h [σ] and from the Hessian, and (ii) the distributions depend on σ only through its overlap q with the special direction w 0 (in absence of the signal, the distribution turns out to be independent on σ). These features make the computation of the annealed complexity feasible; in particular, (ii) is crucial, as it allows to integrate out the variable σ and get an expression for N N which depends only on few parameters. Moreover, it suggests that the distributions of the random vector g [σ] and of the random matrix H [σ] satisfy some rotation invariant symmetry, hinting at the connection with the random matrix theory of invariant ensembles [8,10]. When the number of replicas is larger than one, the situation is more involved, as the random variables associated to different replicas are non-trivially correlated. However, it remains true that their joint distribution can be parametrized in terms of q and few additional order parameters, that are the overlaps q ab = σ a · σ b between the different replicas [37]. In the following subsection, we discuss in more detail this structure, which allows us to re-express the moments (17) as an integral over the order parameters q ab of three terms scaling exponentially with N , see Eq. (24). The first term is a volume factor, emerging when integrating over the variables σ a : this is evaluated with standard methods in Sec. V D. The second terms is the joint distribution of gradients and energy fields; the difficulty in computing this term relies in the inversion of the correlation matrix of the gradients: we overcome it by realizing that it is sufficient to invert the projection of the matrix on a restricted portion of replica space, see Sec. V E. Finally, the third term is the conditional expectation value of the product of determinants. We find that the conditioned Hessians of the various replicas are coupled, weakly perturbed GOE matrices, such that their mutual correlations can be neglected when computing the expectation value to leading order in N (Sec. V F 1 and Sec. V F 2). As a consequence, we find that this term contributes with a factor that is independent on q ab , and which is governed by the properties of the GOE invariant ensembles, see Sec. V F 3. We discuss the stability of the stationary points, which is encoded in the statistics of the spectrum of the Hessians, in Sec. V G. The final result of the calculation is Eq. (13), where we remind that the integral over the order parameters has been performed within the saddle point approximation, assuming a replica-symmetric structure of the overlap matrix, q ab ≡ q for a = b.

C. Structure of covariances and order parameters
As a first step, we analyze the structure of the correlations between the random variables H a , g a and h a : since they are Gaussian, their statistics is fully determined by their averages and mutual covariances, which turn out to depend only on q and on the overlaps q ab = σ a · σ b . To uncover this structure, we consider the gradients ∇h a ≡ ∇h[σ a ] and Hessian ∇ 2 h a ≡ ∇ 2 h[σ a ] of the functional (16) extended to the whole N -dimensional space [38], and determine the covariances between their components along arbitrary directions in the Ndimensional space, given by some N -dimensional unit vectors e i . From here, the correlations of the components g a and H a are easily determined setting is an arbitrarily chosen basis of the tangent plane at σ a . This follows from the fact that g a is an (N − 1)-dimensional vector with components g a α = ∇h a · e a α , which is obtained from ∇h a by simply projecting it onto the tangent plane. Similarly, H a is an as it follows from imposing the spherical constraint with a Lagrange multiplier [39]. For arbitrary e i , taking the derivative of (16) and computing the expectation value we find: and: where the subscript "c" indicates that the correlation function is connected. The same computation for the second derivatives gives: and Finally, the correlations between Hessians and gradients read: Consider first the case of a single replica: choosing e i → e α [σ] to be vectors in the tangent plane, using (19) and e α [σ] · σ = 0 one sees that g [σ] is uncorrelated from H [σ] and h [σ]; moreover, irrespectively of the choice of the basis in the tangent plane, the components of the gradient are independent Gaussian variables with variance p, while the Hessian is a GOE matrix with variance p(p − 1), shifted by a random diagonal matrix. For more than one replica, correlations arise because of the non-zero overlaps between some directions e a α in the tangent plane at σ a and the other replicas σ b . However, the correlations of the components along directions that are orthogonal to w 0 and to all the σ a hugely simplify. To exploit this, it is convenient to separate the Ndimensional space embedding the sphere into the (n + 1)dimensional subspace S spanned by the vectors w 0 and {σ a } n a=1 , and its orthogonal complement S ⊥ . The reference frame of the embedding space, which we denote with , can be chosen in such a way that the last (n + 1) vectors x N −n · · · , x N are a linear combination of w 0 and of all the σ a , forming an orthonormal basis of S, while the remaining N −n−1 vectors x 1 , · · · , x N −n−1 generate S ⊥ . Similarly, the basis vectors in the tangent planes e a α can be chosen so that the last n vectors e a N −n , · · · , e a N −1 , together with the normal direction σ a , are a basis for S, while the remaining e a α with α < N − n generate S ⊥ . In particular, these can be chosen equal for any a, as e a α = δ α,i x i for i < N − n. With this choice, the covariances between the first N − n − 1 components of the gradients do not depend on the corresponding directions e a α , and depend trivially on the overlaps q ab . The covariances between the last components are instead more complicated functions of q ab , which depend explicitly on the choice of the basis in S. Optimal choices for the basis can be made, to simplify the calculations; we discuss an example in Appendix IX C. Regardless of these choices, Eq. (17) can be rewritten in terms of the overlaps alone, as: where E and pQ are the expectation value and the joint distribution in Eq. (17), now expressed as a function of the overlap matrixQ with components while (26) is an entropic contribution. We determine the leading order term in N of each of the three contributions in (24) for q ab ≡ q, and subsequently perform the integral with the saddle point method. To simplify the calculation, we choose the bases x i and e a α so that only one vector has a non-zero overlap with the special direction w 0 : this can be done setting x N = w 0 (hence the name North Pole), and choosing e a N −1 to be the projection of w 0 on the tangent plane of σ a , e a The term V (Q, q) is a phase space factor, which accounts for the multiplicity of configurations of replicas satisfying the constraints on the overlap. Its large-N limit can be obtained from the representation: whereΛ andQ are n × n matrices in replica space with elements Λ ab = (1 + δ ab )λ ab and Q ab = (1 − δ ab )q ab + δ ab , and v(Q, q; µ,Λ, σ) = a≤b iλ ab σ a · σ b − q ab + a iµ a (σ a · w 0 − q). Performing the Gaussian integrals over the variables σ a i and µ a we get: After rescaling −iΛ → NΛ , the remaining integral can be computed with a saddle point (which gives Λ =Q −1 ), leading to: where within the RS ansatz: To leading order in N n, we find: This contribution is dominated by q = q 2 , which corresponds to configurations in which the replicas are almost independent with each others, correlated only through the constraint on q (indeed, it corresponds to replicas having zero mutual overlap in the portion of phase space that is orthogonal to the special direction w 0 ). These configurations are the most numerous, and reproduce the phase space factor obtained in the annealed calculation (when n = 1), since in that case: However, they are disfavored by the other terms in (24), which depend non-trivially on q; the competition between these terms leads to a more complicated global saddle point solution q SP .
E. Joint density of the gradients and energies and pQ(0, ) We now determine the joint distribution pQ(0, ) of the (N − 1)n + n components (g a α , h a ). This can be obtained from the joint distribution of the gradient components ∇h a i = ∇h[σ a ]·x i in the enlarged, N -dimensional space, whose covariances read (see Eq. 21): and averages ∇h a i = −δ iN √ 2N rf k (q). The joint density of the ∇h a i is thus: and where [Ĉ −1 ] ab is the ab block (in replica space) of the inverse covariance matrix, of dimension N × N . Due to our choice of the reference frame giving the covariances between the gradients components in S ⊥ , whileB ab is an (n + 1) × (n + 1) block whose elements are the covariances of the gradients components in S, which are explicit functions of q and of σ a i . To leading order in N this smaller block can be neglected for the computation of the normalization, and one gets: To get the statistics of the components on the tangent planes, we consider the and it is thus related to the value of the functional at the point σ a . The vectors ∇h[σ a ] andg a are related by a unitary rotation: the joint density pQ(0, ) of (g a , h[σ a ]) evaluated at (0, √ 2N ) is easily obtained from (29) with a change of variables.
To determine (29), we introduce the N n-dimensional vectors ξ 1 = σ 1 , · · · , σ n and ξ 2 = (w 0 , · · · , w 0 ) (we remind that we chose x N = w 0 ), so that: and Note that in (32) the matrixĈ −1 is contracted with the vectors σ a , so that the quadratic form depends only on the overlaps q, q. The exponent (32) can be explicitly computed noticing that the vectors ξ 1 and ξ 2 , together with the vector ξ 3 = a =1 σ a , · · · , a =n σ a , form a closed set under the action of the matrixĈ −1 : the inversion of the correlation matrix can be performed in the restricted subspace spanned by these three vectors, and the matrix elements ofĈ −1 within this subspace suffice to get (32). We refer to the Appendix IX B for the details of this computation. As a result, we obtain where D(q) is given in (11). In the limit of a single replica n → 1, the quadratic form reduces to: which is consistent with (20), as it reflects the factorization of the distribution of the gradients and of the rescaled energy fields: the first term in (34)  p,k ≡ n Q p,k + O(n 2 ), we obtain and combining (29), (30) and (35) we get This term is dominated by an energy dependent value of q = q( , q). As we argue in the following section, to leading order in N the expectation value E turns out to be independent on q, so that (36) is the term responsible for shifting [40] the saddle-point solution away from the value q = q 2 maximizing the phase space term (27).
F. The expectation value of the product of determinants E The expectation value E in (24) is over the joint distribution of the Hessian matrices H a , conditioned on a particular value of the gradients g a and field h a . Using the identities (19) and (31) we get that the Hessians can be written as: Here αβ ≡ e a α · ∇ 2 h a ps · e a β , andg a N are random variables. When conditioning to h a = √ 2N , the random variableg a N becomes a deterministic function equal to √ 2N u( , q), see (33), so that the conditional law of H a can be easily obtained from the one of the matrices M a . In the following, we discuss the conditional law of M a . We denote withM a the random matrices obeying this law, and similarly forH a . As we show below, the exponential scaling of E is determined by the leading order term (in N ) of the density of states of the conditioned matricesH a / √ N . This is simply the density of states ofM a / √ N , shifted by the constant term √ 2u( , q): indeed, the term proportional to θ r,k (q) is a rank-1 perturbation that modifies the density of states only to lower order in N , and does not affect the result for E to exponential accuracy in N . Notice however that despite this term is irrelevant in the computation of the number of stationary points, it has to be taken into account when discussing their stability, see Sec. V G and Appendix IX D.

Conditional law of the Hessians
Consider first a single matrix M a : before conditioning to the values of the gradients and energy functionals, the distribution of each M a is the one of a GOE matrix, with independent entries with variance [M a αβ ] 2 = p(p− 1)(1 + δ αβ ), see Eq.(22) (this follows from the fact that the vectors e a i in the tangent plane are orthogonal to σ a ). This distribution is modified by the conditioning, as the entries of M a are correlated to the gradients and energies of all the other replicas. To determine this effect, we partition each (N −1)×(N −1) matrix M a into blocks, where the larger block M a 0 has dimension (N − n − 1) × (N − n − 1), and contains the components M a αβ along directions α and β that both belong to the subspace S ⊥ , M a 1 is n × n and contains the components with both α and β belonging to S, and M a 1/2 contains the remaining, mixed components. The same partitioning can be done for each gradient vector g b = g b 0 , g b 1 . As we argue in Appendix IX C, the block structure (38) is preserved when conditioning, meaning that correlations between components in different blocks are not induced. Moreover, from (23) it appears that the only components that are affected by the conditioning are the ones in the blocks M a 1/2 and M a 1 , which are correlated with the vector g b 0 , and with the vector g b 1 and the fields h b , respectively. In particular, the conditioning induces correlations between the components of M a 1/2 that belong to the same row, and between all the components in the smaller block M a 1 . Similarly, non-zero averages are induced for the components in the block M a 1 . As a result, M a can be written as where the blocks are independent with respect to each others, the largest one M a 0 has a GOE statistics, while the others have correlated entries. Such correlated entries, since their number is not O(N 2 ), do not impact the density of eigenvalues ofM a in the large-N limit, which is instead dominated by the larger block, and is thus a semicircle law, see Appendix IX C. Since, as we argue in the following subsection, the eigenvalues density is the only object needed to compute the expectation value in (24) to leading order in N , it follows that the correlations induced by the conditioning can be neglected to exponential accuracy.

Factorization of the expectation value of the determinants
As it follows from (22) associated to the rescaled matricesH a / √ N , as: where Z = a Dρ a exp [F N ({ρ a })] is a normalization. The functional F N ({ρ a }) couples the eigenvalues densities of the different matrices. The crucial point is that the leading order term in F N ({ρ a }) scales as N 2 , as it follows from the fact that the matricesH a are shifted GOE matrices deformed by finite-rank perturbations. When computing the functional integral (40) with a saddle point, the saddle point solutions ρ a sp are determined just by the minimization of F. As such, they coincide with the marginals of the joint distribution in the space of measures, since This implies that to leading order in N , the expectation value of the product of determinants is just the product of expectation values computed with the marginal distribution of the matricesM a , properly shifted according to (37). In turn, each expectation value can be computed by means of the eigenvalues density ofM a .

The GOE computation
Given these observations, and given the symmetry between replicas, the expectation in (40) reduces to where ρ sp (λ) is the eigenvalues density of the rescaled matricesH/ √ N , which coincides with the density ofM/ √ N − √ 2u( , q)1.
Given that the density ofM/ √ N is a semicircle law with support [−2 p(p − 1), 2 p(p − 1)], one finds From here it follows that: where I(y) = I(−y) = dµ 2 − µ 2 log |µ − y|/π is given in (V A), and β( , q) in (10). The contribution of the determinants in (24) thus reads: G. Threshold energy, isolated eigenvalue and complexity of the stable stationary points The results obtained so far suffice to derive the explicit expression for the quenched complexity, since combining everything we get: has to be evaluated, in the limit of large N , at the saddle point q SP ( , q) which maximizesQ p,k , as it is prescribed by the replica method.
To conclude this analysis, it is necessary to address the stability of the stationary points counted by the quenched complexity. Indeed, when presenting the results of the saddle point calculation in the following section, we shall focus only on those stationary points which are local minima, meaning that we restrict to values of the parameters q, for which the typical Hessian of the stationary points is positive-definite.
The density of eigenvalues of the Hessian at a typical stationary point is the one of a random matrix obeying the same law as the conditioned onesH, in the limit n → 0. As a matter of fact, the information on the eigenvalues density of the Hessian H[σ] is encoded in the resolvent function: where N is here the dimension of the matrix and σ is assumed to be a stationary point. The quantity (46) is a fluctuating variable even for fixed realization of the random field, as it changes from stationary point to stationary point. To capture its typical behavior, we first average it over all stationary points at given q, at fixed realization of the field, and subsequently average of the random field itself. This leads to where χ(q, ) = δ(σ · w 0 − q)δ(h[σ] − √ 2N ) enforces the constraints on the overlap with the signal and on the energy density. The above average can be performed by means of the replica trick, using the identity x −1 = lim n→0 x n−1 . Exploiting the replicated Kac-Rice formula, we get where now Proceeding as before, and using that the resolvent of the Hessian at a stationary point σ 1 is a function only of the eigenvalue density ρ 1 (λ), we find that (48) can be evaluated with the same saddle-point calculation discussed above, and: Here, as before, ρ sp (λ) denotes the eigenvalues density of a matrix distributed asH/ √ N , which depends on the parameters q, as well as on the mutual overlap between replicas, evaluated at its saddle point value q SP ( , q). Therefore, to discuss the stability one needs to characterize in detail this density, in the limit n → 0. Note that this computation is a quenched one, as it accounts for the fluctuations between the various stationary points at fixed q, . The annealed approximation would correspond to averaging separately the numerator and the denominator in (47) over the random field. This does not account for the correlations between the stationary points, and it is reproduced setting n = 1 in the above expression (instead of taking n → 0). To leading order in N , the density of states ofH/ √ N is dominated by the large GOE block, and it is therefore the shifted semicircle in Eq. (43): stationary points for which part of the support of the semicircle lies in the negative semi-axis have an extensive number ∼ O(N ) of unstable directions in phase space. Since the location of the support of (43) depends on the energy density , a threshold energy can be defined, as the energy at which the support touches zero: For fixed value of q and of the parameters p, k, r, stationary points having energy larger than (51) have an Hessian with extensively many negative eigenvalues, and therefore can not be considered as trapping minima of the energy landscape. The semicircle law accounts for the continuous part of the spectrum of the Hessian, and does not give information on possible isolated eigenvalues ofH/ √ N , that may not belong to the support of (43). Such eigenvalues correspond to isolated poles of the resolvent R(z); they contribute to the density of states with subleading terms of order N −1 , and are thus irrelevant when computing (45). However, if an isolated eigenvalue exists, it can become smaller than zero at values of < th (q), leading to the instability of the point σ along some direction in phase space. For the matricesH, isolated eigenvalues can be generated by the subset of entries that do not belong to the large GOE block (as they are distributed with a different average and variance, induced by the conditioning to the gradients), as well as by the rank-1 perturbation proportional to θ r,k (q) in (37). As a matter of fact, for large random matrices perturbed by low-rank operators, it is known that when the perturbation exceeds a critical value, the extreme eigenvalues detach from the boundary of the support of the density of states of the unperturbed matrix. This transition is akin to the one proved in [53] for the Wishart ensamble and known as the BBP transition, and it has been shown to occur under quite general conditions [54,55].
To inspect whether such a transition occurs in the case under consideration, we need to characterize in more detail the law of the matricesH/ √ N . This requires to explicitly compute the averages and covariances of all the entries ofH/ √ N with respect to some fixed basis in the subspace S, see the details in Appendix IX C. This in turn allows to derive general equation for the isolated poles of R(z) in the limit N → ∞, see Appendix IX D. As a result, we find that the isolated eigenvalues, if any, are given in terms of the roots of a third order equation having coefficients that are functions of the parameters q, and of the saddle point value q SP (q, ). Imposing the minimal eigenvalue to be zero gives the boundary of stability of the stationary points: for fixed r and q, it defines a critical stability energy st (q, r) such that the typical stationary points are local minima for < st (q, r), are saddles of finite index for st (q, r) < < th (q, r), and are saddles of extensive index otherwise. In particular, for the values of parameters that we inspect, we find that in the regime st (q, r) < < th (q, r) there is one single isolated eigenvalue that is negative, whose eigenvector has a finite projection on the direction of w 0 ; thus, this instability is an instability toward the signal. In summary, the complexity of the stable stationary points is therefore given by (13), endowed with the condition of the energy being smaller than the threshold energy (51) or, when the isolated eigenvalue exists, of the stability energy st (q, r).

H. Mapping between complexity at different k
Before presenting the results, we derive the mapping (15) relating the complexity for different values of k. Suppose that σ is a stationary point of the functional (16) for a fixed k and for a given value of r = r k (we now make explicit the dependence on k and r writing h k,r [σ]), with overlap q and with energy density = k . Then, the point σ is also a stationary point of the functional (16) with k = 1, provided that r = r eff 1 is chosen so that: In this case, σ has overlap q with w 0 , and has energy density: Indeed, for σ to be a stationary point at a fixed k, it must hold ∇h k,r k · e α [σ] = 0, which implies: where we exploited our choice of bases x N = w 0 and e N −1 [σ] = (x N − qσ) / 1 − q 2 . Moreover, h ps [σ] = k +r k f k (q). This in turn implies that ∇h 1,r1 [σ]·e α [σ] = 0 for α < N − 1, while ∇h 1,r1 [σ] · e N −1 [σ] = (r k f k (q) − r 1 f 1 (q)) 1 − q 2 . Thus, this is a stationary point for k = 1 if r 1 is chosen as (52). In this case, it is easy to check that its energy density equals (53). It follows from this that the knowledge of the curves (13) for k = 1 is sufficient to reconstruct the curves at any larger k, via the mapping (15). We however remind that the analysis of the instability of the stationary point induced by the isolated eigenvalue is strongly dependent on k, and thus has to be performed separately for any case.

I. The results of the Kac-Rice calculation
We are now ready to discuss concrete results. In the following we report the curves resulting from the computation of (13), focusing on the cases p = 3 and p = 4. For each of the values of k that we consider, we find the following general features: as long as r < r c (and, in most cases, also for r > r c ), there are values of q, for which the quenched complexity is positive and the typical Hessian of the stationary points is positive definite, indicating the presence of exponentially many local minima of the energy functional. In particular, at fixed latitude q this occurs over a finite range of energies * (q) ≤ ≤ th (q), with th (q) replaced by st (q) whenever the isolated eigenvalue exists. We find that Σ p,k ( , q) is monotone increasing in this energy range, implying that the most numerous stable stationary points at a given q are the ones at higher energy. At the other extreme of the support * (q), the quenched complexity vanishes, Σ p,k ( * (q), q) = 0. We denote with * (r) the absolute minimum of the energies over all q, and with q * (r) the corresponding latitude; these values coincide with the ones found solving the RSB equations in Sec. IV B. We use the notation q num (r) for the latitude where the largest number of stationary points is found, for any fixed r. At the transition point r c and at the latitude q c given in (12), the support of the positive part of the complexity shrinks to a single point * = th = c , where Σ p,k ( c , q c ) = 0. Moreover, the whole complexity curve at this latitude coincides with the annealed one, Eq. (12). The same remains true for larger r: the annealed complexity is exactly zero at values of q * , * which coincide with the solution of the RS limit of the saddle point equations in Sec. IV B, and which give the latitude and energy of an isolated minimum of the energy landscape. For k > 1 and for some values of r, beyond this isolated minimum there is a residual band containing exponentially many local minima, at smaller overlap q with the signal. In the following, we present in more detail the results for each of the cases presented qualitatively in Sec. III.

Case I
Instances of the complexity curves Σ p,k ( , q) in the case k = 1, f 1 (x) = x are given in Fig. 6, for p = 3 and fixed r < r c . The curves are obtained solving numerically the saddle point equations for q for each value of the parameters q, . For k = 1 we find that there is no isolated eigenvalue exiting the bulk of the semicircle: thus, for each q the maximal energy where stable stationary points are found is th (q), which is marked with the squares in Fig. 6. The curves show the following trend: below a minimum value q m , the complexity is positive only for the states which have energy above the threshold, and are therefore unstable. At q m , the equality * (q m , r) = th (q m , r) holds, meaning that at this latitude there are only marginally stable (and unstable) stationary points. For larger latitudes, as q increases the energy interval in which the complexity is positive (and the points are stable) gets wider and moves toward smaller energies, until the maximal width is reached at q num . At larger q, the energy interval start shrinking, and the minimal energies * (q) decrease until the absolute minimum is reached at q * ; for q > q * the trend is reversed and * (q) starts increasing, until it collapses to th (q) at q M . Analogous results are obtained for different values of r below r c , as well as for p = 4. In Fig. 7, we plot the bands q m (r) ≤ q ≤ q M (r) con- taining exponentially many local minima, as a function of r. These bands correspond to the red ones plotted pictorially in Fig. 1. For each of the q within the bands, the quenched complexity behaves as in Fig. 6. As r increases, the bands gets wider and subsequently shrink and collapse to q c at r = r c (p), corresponding to the black points in the figures. Here the minimum becomes unique, and it is marginally stable. This landscape phase transition at r c is signaled by the fact that the saddle point solution q SP converges to 1, meaning that all the replicas coincide, and that the quenched complexity be-comes equal to the annealed one. This corresponds to the recovery of the RS symmetry in the replica calculation of Sec. IV B.

Case II
According to the analysis of Sec. IV and of Ref.
[30], for k = 2, f 2 (x) = x 2 /2, the minima of the energy landscape undergo a second order transition at r = r 2ND < r c . The transition marks the boundary between two different behaviors of the complexity curves, see Fig. 8: for r < r 2ND , the energy interval containing exponentially many states is maximally large at the equator, where both the deepest and the most numerous states lie. For r 2ND < r < r c , instead, the most numerous states remain at the equator and have = th , but the deepest states move toward a higher overlap q * > 0 with the signal. At r 2ND , the states of minimal energy * (r) detach from the equator, moving toward larger latitudes. The features of the bottom of the landscape (that is, the spectrum of the minimal energies * (q, r), the thermodynamic energies * (r) and the value of r 2N D ) can all be obtained from the corresponding k = 1 curves * 1 (q, r 1 ) satisfying Σ p,1 ( * 1 (q, r 1 ), q; r 1 ) = 0, as we discuss in more detail in Appendix IX E. Consider now the other transition in the energy landscape, which occurs when the strip containing exponentially many stationary points splits into three different bands (in the following, we restrict to positive values of q: due to the symmetry, the landscape at negative overlap is specular to the one at positive overlap). The strip containing the stationary points for k = 2 can be identified exploiting again the mapping (15), with the caveat that the stationary points so determined are stable only in the sense of the threshold, and the analysis of the sign of the isolated eigenvalue has to be performed separately. We give the details of the mapping in Appendix IX E. The resulting bands are shown in Fig. 9, where one sees that they split at r ≈ 1.55 for p = 3, and r ≈ 2.05 for p = 4. For r larger than this splitting point, the band closer to the North Pole, which is the one containing the deepest minima, shrinks until it collapses to a single state at the RS transition, while the band enclosing the equator, which is the one containing the most numerous minima, shrinks to zero only asymptotically (this corresponds to the dashed lines in Fig. 9). This implies that, for any value of r, there is a strip of finite width and small overlap q with the signal, containing exponentially many stationary point with energy smaller than the threshold. To conclude the analysis of the landscape, it is necessary to investigate the possible instability of these points due to the presence of a negative, isolated eigenvalue. For k = 2, the isolated eigenvalue exists only for sufficiently large r, and it renders unstable, for each q for which it exists, the stationary points at higher energy st (q, r) < < th (q, r). We refer to the Appendix IX E for a more detailed analysis of this instability, and report here only its main consequences. First, if this instability is accounted for, we find that for r > r c the most numerous non-unstable points are still at q = 0, but are no longer marginally stable. Rather, they have an energy st (0, r) smaller than the threshold energy, and have one flat direction in their Hessian, corresponding to the isolated eigenvalue being zero. For general q, as r increases st (q, r) decreases, until it becomes smaller than the lower bound * (q, r), implying that all the points at the given latitude are unstable because of a single negative eigenvalue (See Fig. 16 in Appendix IX E). This happens first for the larger values of q belonging to the band: thus, the band of those stationary points gets narrower around the equator, from above. At a finite value of r (r ≈ 2.06 for p = 3 and r ≈ 3.21 for p = 4), also the last stationary points at the equator become unstable (this value of r can be computed within the annealed approximation, see the comments at the end of Appendix IX D). For larger r, there is a unique stable minimum, that is the minimum of the annealed complexity.

Case III
In this case, the transition at the bottom of the energy landscape is of first order. What distinguishes the two options presented in Sec. III is whether this thermodynamic transition occurs before or after the band of stationary points separates into two distinct strips, and the strip at larger overlap q undergoes the RS transition. The second case (Option B ) is realized, for instance, for k = 3, p = 4. In this case, the curves Σ p,k ( , q) behave in the following way: for small r, they are monotone decreasing for increasing q (they look like their k = 2 counterpart in Fig. 8 (a)), so that both the deepest and the most numerous states are at the equator. At a spinodal point r 1SP ≈ 2.01, a local minimum in * (q) appears at a latitude q * 2 (r) > 0, so that for r > r 1SP the curves are no longer monotone, see Fig. 10 (a). The absolute minima remain however at the equator, q * = 0. The latitude q * 2 of the second minimum increases with r, and its energy decreases; at the first order transition r 1ST , its energy become smaller than the energy of the minima at the equator (that is the ground states of the unperturbed p-spin model), and q * jumps discontinuously from zero to a finite value q * 2 (r 1ST ), see Fig. 10 (b). The value of r 1SP , the latitudes of the second minima q * 2 (r) and the corresponding energies can be obtained via the mapping from the curves at k = 1, as we discuss in Appendix IX E. The bands of latitudes corresponding to positive complexities below the threshold energy can also be obtained from k = 1, in a way analogous to the one discussed in the Appendix for k = 2. A major difference with respect to Case II concerns the effect of the isolated eigenvalue, since for k ≥ 3 the states at the equator are not destabilized by it (see the details in Appendix IX E). Thus, in this case the threshold states of the unperturbed p-spin model are the most numerous stable minima, for any r. This case is summarized in Fig. 11 (b). Finally, we consider the case k = 3, p = 3, which realizes Option A of Sec. III. In this case we find that r 1SP = r c . The curves Σ p,k ( , q) behave similarly to the ones in Fig. 8 (a) for any r < r c ≈ 2.45. As r approaches r c from below, the band of stationary points rapidly grows, and at r c it reaches its maximal width, incorporating q c (i.e., q M (r c ) = q c ). Exactly at this latitude q c , the saddle point q SP = 1 reaches one, and the quenched complexity becomes equal to the annealed one, having positive support for a single value of the energy density c . The curve of minimal energies * (q) has a minimum at q = 0, and it is flat at q c , where it intersects the threshold energy th (which for p = k is independent of q and r, and equals to the threshold of the unperturbed p-spin model). Therefore, the second minimum of * (q) appears exactly at r c , and at this point it coincides with the RS solution. At larger values of r, the minimum of the annealed complexity is isolated (it departs from the band containing all the other minima), and becomes energetically favorable at r 1ST ≈ 2.56. The band at small overlap shrinks asymptotically around the equator. Thus, in this case the band of minima is connected up to r c , and it splits exactly at the critical point, see Fig. 11 (b). The analysis of the isolated eigenvalue shows that for large enough r, the eigenvalue renders unstable the points at higher overlap in the strip enclosing the equator, but it does not affect the most numerous, marginally stable states at the equator, nor the minimum of the annealed complexity, which is stable for any r > r c .

VI. COMPARISON BETWEEN KAC-RICE AND REPLICA METHOD
As pointed out in the previous section, the information on the thermodynamics provided by the replica calculation is fully recovered from the Kac-Rice results, by analyzing the spectrum of minima * (q, r) satisfying Σ p,k ( * (q, r), q; r) = 0. As first pointed out in [23], the thermodynamical replica method can also be used to obtain information on the number of critical points. In this section, by comparing the predictions of the two calculation schemes concerning the configurational entropy, i.e. the complexity of the most numerous stationary points Σ c ( ), we show that the thermodynamical replica method is not able to reproduce the full Kac-Rice results and leads to partially incorrect predictions. This is an important point since, although the method could be probably amended, its present form which is often used for the purpose of computing the configurational entropy fails for the models we consider. We highlight below the two different reasons for failure. The replica formalism allows to sample local minima at energy higher than the equilibrium one by not imposing the saddle point condition on βm (the third among Eqs. 3), and using m as a parameter, which plays the role of an effective inverse temperature. By lowering m (βm in the T = 0 case) the remaining saddle point equations describe the macroscopic features of the most numerous local minima at higher energy. In particular, the expression contained in the disregarded saddle point equation (3) gives the corresponding intensive log multiplicity of these minima, i.e. their configurational entropy S c : The stability of the metastable minima, whose multiplicity is accounted for by (54), is checked by analyzing the stability with respect to fluctuations in the overlap matrix Q ab , see Appendix IX F. The entropy S c can then be compared with the Kac-Rice complexity of the most numerous stationary points. The latter is obtained, for each energy density , as the maximum of the curves Σ p,k ( , q) over those latitudes q that correspond to stationary points that are stable at the energy . We find that there are regimes in which the two calculations are not equivalent, with the replica calculation failing to identify part of the complexity curve resulting from Kac-Rice. As an illustrative case, we consider the parameters k = 2, p = 4. For r < r c , the stability of the stationary points at fixed latitude is determined only by the bulk of the eigenvalues density of the Hessian, since no isolated eigenvalue is present. Therefore, the constrained maximization of the Kac-Rice complexities reads: As long as r < r 2ND , at fixed the curves Σ( , q) are monotone decreasing in q, with a maximum at q = 0. In this case, Σ c ( ) coincides with the complexity of the stationary points at the equator, and the quantity (54) reproduces it. For r 2ND < r < r c , the two complexity curves coincide only in the lowest part of the energy domain (see Fig. 12 FIG. 12. Comparison between the configurational entropies Sc( ) and Σc( ) obtained with the replica and Kac-Rice calculation, respectively. The two curves coincide below ≈ −1.245 (blue point), which is the energy at which the replica solution with q = 0 becomes unstable. For higher energies, the curve computed with replicas is contributed by the states at q = 0 (which are a solution of the replica equations for every energy density), while the Kac-Rice curve is contributed by marginally stable states whose q does not satisfy the stationarity condition of the replica action. for a comparison between the curves obtained with the two methods for r = 0.9). More precisely, the curves coincide for the energies for which the maximum in (55) is attained inside the interval, at a q s satisfying < th (q s ). This means that the most numerous states at these energies have Hessian gapped away from zero, and are at latitudes satisfying ∂Σ p,k ( , q s )/∂q = 0. In this case, q s coincides with the value of q selected by the saddle point equations of the replica calculation (see Sec. IV B), and we recover S c ( ) = Σ c ( ). In the second part of the curve, instead, the maximum is attained at the boundary of the interval, at latitudes q s such that = th (q s ). This part of the curve Σ c ( ) is thus contributed by points that are marginally stable, and which do not fulfill the stationarity condition ∂Σ p,k ( , q)/∂q = 0. This piece of curve is not recovered by the replica scheme: rather, in this energy regime the replica solution corresponding to saddle point values q = 0 results in a different entropy curve (the black-dashed line in the Inset of Fig. 13) that has to be disregarded, being unstable with respect to the replicon criterion recalled in Appendix IX F. On the other hand, the dashed blue line in the figure corresponds to q = 0, which is always a solution of the saddle point equations in the replica calculation, that coincides with the Kac-Rice curve only for the highest energies. Thus, the replica result is inconsistent with the Kac-Rice one at intermediate energy densities. As r increases toward r c , the interval of energies in which the two curves coincide shrinks, so that the replica calculation allows to recover only a very small portion of the configurational entropy obtained via Kac-Rice, the part of curve contributed by strictly stable points. The situation outlined above highlights the first way in which the replica method can fail: the correct result is recovered only when the largest contribution to the configurational entropy at fixed energy is given by a q such that ∂Σ p,k ( , q)/∂q = 0. The reason is that the configurations taken into account by the replica method if ∂Σ p,k ( , q)/∂q = 0 do not correspond to true minima since they have a non-zero gradient in the north-pole direction. Let's now focus on the other way in which the usual replica method to compute the configurational entropy can fail. For r > r c (but smaller than the value of r at which the landscape becomes completely convex), for the smaller energies the curves Σ( , q) are monotonic in q, with a local minimum at q = 0 and a global maximum at a latitude q > 0 (see the Inset in Fig. 14), while at larger energies the minimum at the equator becomes the maximum of the curve. Therefore, at small energies Σ c ( ) is contributed by the points at the latitudes q > 0 of the maximum. The replica calculation reproduces nevertheless only the complexity at q = 0, even when there is a full spectrum of more numerous (stable) points at higher overlap with the signal (see Fig. 14 for the case r = 3). The reason for this is that precisely at the latitude where the complexity has a maximum, the isolated eigenvalue of the Hessian is exactly equal to zero. Thus, the lower-energy part of the curve Σ c ( ) obtained from Kac-Rice is contributed by stationary points that have the Hessian with a single zero mode, which are known to be not captured by the standard replica calculation [24][25][26]. The physical reason is that these stationary points do not correspond to the zero-temperature limit of stable states. In fact, as shown in [56], they correspond to minima characterized by finite barriers. This situation has been already found in computation of the multiplicity of TAP solutions in mixed models [24] and in models exhibiting a full replica symmetry breaking phase [25, 26].
FIG. 14. Comparison between the complexity of the most numerous states at fixed energy computed with the replica and Kac-Rice calculation, respectively. The two curves coincide above ≈ −1.255, where the curves are contributed by points at q = 0. At lower energy densities, the replica equations reproduce the complexity of the points at q = 0, while from the Kac-Rice calculation it emerges that there are more numerous stable points at higher overlap q > 0, among which the most numerous ones have an Hessian with a single zero mode. Their complexity is given by the red curve. Inset. Complexity as a function of the latitude q for a fixed value of energy = −1.266. The curve has a maximum at q ≈ 0.18, where the isolated eigenvalue of the Hessian is exactly zero. The points at larger latitudes are saddles of index one.

VII. ON THE SPIKED TENSOR CASE
As we have previously remarked, the Hamiltonian H p,k (r) in the case k = p is related to spiked tensor model [32], i.e., to the inference problem of detecting a low-rank, additive perturbation of a symmetric Gaussian tensor, which has attracted a lot of attention recently. In this section we specifically present our analysis on this system, focusing in the case k = p ≥ 3 for concrete results. Some of these observations are already stated in [32][33][34][35]. We also discuss which properties the annealed computation [15] cannot capture. The inference task in the spiked tensor problem consists in reconstructing the unknown vector v 0 from the observation of a random p-tensor with components where the random couplings J i1,...,ip , symmetric with respect to a permutation of the indices, correspond to the noise and v 0 (the signal) is generated at random from a spherical prior distribution. In particular, one is interested in identifying the strong detection threshold [34], i.e. the critical signal-to-noise ratio below which the spiked model is statistically indistinguishable from the un-spiked one with r = 0, and the detection threshold, above which an estimatorŝ of the signal having a finite overlap with v 0 in the limit N → ∞ exists (for a precise definition of statistical indistinguishability see [34]) . In the matrix case p = 2, the two thresholds coincide [43][44][45]. They are given by the signal-to-noise ratio at which the smallest eigenvalue of the matrix pops out from the semi-circle; the corresponding eigenvector is correlated with the signal. In the tensor case, rigorous bounds on both thresholds are given in [32,34,44], while the sharp threshold given by the Minimal Mean Squared Error estimator is determined in [33].
The connection with the analysis presented above emerges when considering the maximum-likelihood estimatorŝ ML of v 0 [32]. It is immediate to see that this corresponds to the vector that maximizes the injective norm of the tensor: where ·, · denotes the tensor product. This coincides, up to a global sign flip of the energy functional H p,k (r), with its absolute minimum; therefore, a reliable estimate of the signal by means ofŝ ML is possible whenever the global minimum of the landscape acquires a non-zero overlap with the special direction of the signal, i.e., whenever r ≥ r 1ST . The thermodynamic transition thus gives (in general) and upper bound to the detection threshold. On the other hand, the performance of algorithms [32] aiming at reconstructing the signal is expected to depend on the full structure of metastable states, encoded in the complexity. The analysis presented in the previous sections and in  lead to the following picture for the spikedtensor model (p = k ≥ 3): (i) The spinodal point r 1SP , where a high-overlap metastable minimum appears in the curve * (q) (or, equivalently, where a second solution appears for the replicas equations) is exactly equal to the point r c (p) where the trivialization of the portion of the landscape at high overlap with the signal occurs. Moreover, there is no splitting of the band of minima for r < r c . This implies that the portion of landscape close to the high-overlap minimum is not rugged for all rs such that the highoverlap minimum exists, and no intermediate phase with metastability at high-overlap with the signal is present.
(ii) The transition points r c and r 1ST are related, respectively, to the dynamical and statical transition temperatures (β d and β s ) of the pure spherical pspin model; more precisely: (iii) For r 1SP ≤ r ≤ r 1ST and for most energy densities , the complexity is non-monotonic in the overlap q; however, for any the most numerous minima are found to be orthogonal to the signal, at q = 0.
The equality r 1SP = r c is shown in Appendix IX E. It implies that for p = k the thermodynamic transition always occurs when the high-overlap part of the energy landscape is convex. This allows for the annealed Kac-Rice computation [15] to correctly capture the transition value r c although quenched and annealed complexity do not coincide for r < r c . For other models, such as p = 3, k = 4 for which the landscape is rugged close to global highoverlap minimum at r 1ST , this is no longer the case and the quenched computation is needed to also correctly describe the transition. The first identity in (58) can be read explicitly from Eq. 4. The second identity is naturally true for Bayes optimal estimates [33], as it holds in general along the Nishimori line (which corresponds to the line, in the (r, β) phase diagram, where β = 2r/p). The fact that the same detection threshold is found with the maximum likelihood estimator follows from the properties of the thermodynamic phase diagram [30], where the first-order transition line appears to be independent of temperature, thus implying that the same r 1ST found on the Nishimori line is recovered at T = 0 (note however that this is a peculiarity of the spherical case and does not hold in general [33]). The property (iii) implies that the quenched complexity Σ c ( ) is identical for the spiked and the un-spiked model for r SP ≤ r ≤ r 1ST , which is consistent with the strong detection threshold being at r 1ST . Given the structure of the energy landscape, we expect that for all r not diverging with N , physical dynamics starting from random initial conditions behaves as in the un-spiked model, i.e., the system remains stuck in the vicinity of the most numerous, marginally stable states that lie at the equator q = 0, thus being unable to recover any information on the signal. The approximate message passing algorithm is known to fail as well [33]. On the other hand, for r > r 1ST the dynamics with a warm start should converge to the global minimum of the energy landscapes over time scales of O(1), due to the smoothness of the landscape in its vicinity [46]. Polynomial-time algorithms are instead known to succeed for r scaling as N (p−2)/4 , see [32] and references therein. Finally, it is proven in [35] that for the Ising spiked tensor defined on the hypercube (s i = ±1), the strong detection and detection threshold coincide, being both equal to the threshold given by the minimal mean square error estimator [33]. The proof relies on the bound (for large N ) of the fluctuations of the free energy of the Ising p-spin model around its average value, in the high-T phase. A similar bound should hold for the spherical case, since the variance of the intensive free-energy is found to be of order 1/N by the replica method (the variance can be directly obtained using the RS approximation to compute the O(n 2 ) term of the replicated free energy [3]). In consequence, we expect that this argument can be extended to the spherical case, thus implying that both thresholds are given by the maximum-likelihood estimator.

VIII. DISCUSSION AND CONCLUSION
We have analyzed the evolution of an archetypical model of high-dimensional landscapes generated by an energy function in which random fluctuations compete with a deterministic contribution favoring a single minimum. For entropic reasons the overall majority of the minima created by the randomness lie in a region different from the one favored by the deterministic contribution. By increasing the strength of the deterministic contribution, and depending on the form of the latter, different behaviors and geometric phase transitions, that we have classified and thoroughly analyzed, can take place. As discussed in the introduction, our results provide guidelines for current problems in several different fields, and a full analysis of the energy landscape of the spiked-tensor model which recently attracted a lot of attention [15,[32][33][34][35]. In particular, our analysis is useful to understand how the dynamics governed by gradient descent (and stochastic versions of it) proceed in such landscapes. The region of bad and numerous local minima that we called the equator is a trap for the dynamics. Only in case of a sufficiently warm start, i.e. if the initial condition of the dynamics has a finite overlap with the special direction v 0 selected by the deterministic contribution, the system can end up close to v 0 , although not necessarily in the global minimum since many additional good local minima can be present. The other main contribution of our work is methodological. We have developed a framework based on the Kac-Rice method that allows to compute the quenched complexity, opening the way to full analysis of random landscapes in many different contexts. We have shown that it is superior to previous frameworks used in the literature. Indeed, the usual replica method fails in some cases, as demonstrated in this work, whereas the super-symmetry one is in comparison quite obscure. Instead, the Kac-Rice formalism we developed is free of ambiguities, straightforward although complex, and likely to be transformed in a rigorous formalism in a not too distant future.
[40] In absence of the signal (r = 0), the expression (36) becomes independent of q, and it is maximal at q = 0. When optimizing Eq. (24) over q as well, one finds that the complexity is dominated by q = q = 0. This is consisent with the fact that for the p-spin Hamiltonian, the quenched calculation for the complexity reproduces the annealed result, that is indeed recovered setting q = q = 0.
[41] Note nevertheless that they are correlated in blocks: the blocksM a γ correlate only with the correspondent blocks M b γ , for a = b and γ = 0, 1/2, 1.
[42] Let pjoint ({λ a i } | , 0) be the joint distribution of the (N − 1)n eigenvalues of the matrices H a conditioned on gradients and energies. The expectation value of the product over determinants equals with P ({ρ a }) ≡ e −F N ({ρ a }) /Z a measure over the space of functionals ρ a , given by [55] Edwards, S.F, Jones, R.

A. Details on the replica calculation
In this Appendix we provide some additional details on the replica analysis presented in Sec. IV. At finite β, the replicated action S evaluated within the 1RSB ansatz for the overlap matrix Q ab reads: The saddle point equations for the four parameters q 1 , q 0 , m and q equal to: In the zero temperature limit β → ∞, the variables of order one are β(1 − q 1 ) and βm. Performing this limit in the above equations one recovers the expressions given in the main text.
The expansion of the saddle point equations in q 1 − q 0 gives rise to four different equations; three of them are useful to get the variables q, q 0 = q 1 = q and m at the continuous transition between 1RSB structure to RS structure of the high overlap phase. The fourth equation fixes the position of the continuous transition line on the phase diagram T, r, i.e. for a given T it gives the corresponding r c (T ). At zero temperature, when expressed in terms of β(1 − q), βm, and q, the equations read as follows: and Their solution gives the generic expression for r c and q c reported in the main text, Eqs. 4 and 5.

B. Computation of the quadratic form Eq. (32)
In this Appendix we provide some details on the computation of the inverse correlation matrixĈ −1 in (32). As pointed out in the main text, for the purpose of computing the quadratic form (32) it suffices to invert C within the subspace spanned by the N n-dimensional vectors ξ 1 , ξ 2 and ξ 3 , which is closed under the action ofĈ. For convenience, we separate the matrixĈ into its diagonal and off-diagonal parts in replica space, The operator1 +ÔD −1 acts on the chosen vectors as follows: Given thatD −1 ξ 1 = p −1 ξ 1 andD −1 ξ 2 = ξ 2 − q(p − 1)p −1 ξ 1 , the quadratic form (32) can be straightforwardly rewritten in terms of matrix elements of the op- To invert this operator, we introduce an orthonormal basis of the subspace spanned by the vectors ξ i : with A = n(n − 1)(1 − q) 1 − nq 2 + (n − 1)q , and we write Y ij = v i ·Ŷ · v j . In this basis, using f k (q) = q k /k and u( ) = p + r(p/k − 1)q k , we obtain for (32): Using (66) and (67), we find that the operator 1 +ÔD −1 acts on the basis v i as follows: while the relevant matrix elements of its inverse are: and with D(q) given in (11). The result (V E) is recovered substituting these expressions into (68).

C. Conditional distribution of Hessians
In this Appendix, we analyze the structure of the (N − 1)n × (N − 1)n covariance matrix of the Hessians components H a αβ , conditioned to the gradients and energy fields of all the n replicas. We remind that, given (37), the Hessians can be written as where M a denotes the p-spin part of the Hessian. We compute the conditional law of M a , and denote with M a the random matrix obeying this law (and similarly forH a αβ ). Since the last two terms in (70) are deterministic, the covariance matrix ofH a is the same as the one ofM a , while the averages of the components are shifted by the deterministic terms. We show that, for each a, the matricesM a are perturbed GOE matrices; in particular, eachM a / √ N can be written as a sum of a stochastic matrix S a with zero average, and a deterministic matrix D a , The stochastic part S a has the block structure: m 11 m 12 · · ·· · · m 1M m 21 m 22 · · ·· · · m 2M · · · · · · · · ·· · · · · · · · · · · · · · ·· · · · · · m M 1 m M 2 · · ·· · · m M M n M +11 n M +12 · · ·· · ·n M +1M · · · · · · · · ·· · · · · · n N −11 n N −12 · · ·· · ·n N −1M n 1M +1 · · · n 1N −1 n 2M +1 · · · n 2N −1 · · · · · · · · · · · · · · · · · · n M M +1 · · · n M N −1 where (i) the larger diagonal block has size M × M = (N − n − 1) × (N − n − 1) and it is made of elements m αβ that are independent with variance σ 2 = p(p − 1), (ii) for a generic choice of basis in the subspace S, only the elements n αβ belonging to the same row are correlated, (iii) the smaller diagonal block has size n × n and its elements q αβ are all mutually correlated. The deterministic matrix D a is zero everywhere, except in the small n × n block. Form this structure it follows that the reduced density of eigenvalues ofM a / √ N is, to leading order in N , the one of a GOE matrix with σ 2 = p(p − 1), since the fraction of entries having a modified variance and average is vanishing in the large-N limit. This information suffices to perform the quenched calculation given in the main text, as only the density of eigenvalues is needed to compute the quenched complexity to leading order in N . We remark that the partitioning ofM a into blocks and the properties (i-iii) follow solely from the separation of the coordinates into the subspaces S ⊥ and S, and are independent on the choice of the basis in both S ⊥ and S. The covariances of the elements n αβ , q α,β instead depend on the choice of the basis in S: in particular, for a particular choice of the basis, that we shall discuss in the following, the elements n αβ are uncorrelated with eachothers, with variances that differ from the ones of the m αβ . As a first step, we discuss how the general structure (71) is recovered.

Block structure of the conditioned Hessian
To compute the averages and covariances of the compo-nentsM a αβ , we group all the independent components of the unconditioned matrices M a into an nN and M ≡ N − n − 1. The vectors M 0 , M 1 and M 1/2 and have dimension nM (M + 1)/2, n 2 (n + 1)/2 and n 2 M , respectively. They group the Hessians coordinates along directions that belong both to S ⊥ , or both to S, or one to each subspace, respectively. This decomposition reflects the one in (38), except that now we consider all the n replicas. Analogously, we define the nN -dimensional vectorg = (g 0 ,g 1 ), withg γ = (g 1 γ , · · · ,g n γ ), and: g a 0 = (g a 1 , · · · , g a M ), g a 1 = (g a M +1 , · · · , g a N −1 ,g a N ).
We recall thatg a N = ∇h a · σ a = p h a + √ 2N r (pf k (q) − f k (q)q), and thus for any fixed q, conditioning to h a = √ 2N is equivalent to conditioning tõ so that the conditional covariances read: (73) As claimed in the main text, the covariances of the largest blocks M a 0 are left untouched by the conditioning, so that the components of this block form a GOE matrix with variance σ 2 = p(p − 1). We now analyze the structure ofΣ M|g , to show that, for generic choices of the basis in the tangent planes, correlations are induced in the blocks M a 1/2 only between elements belonging to the same line. One has: where S ab is a block of size n × n, equal for every α, with components and with Q ab = δ ab + (1 − δ ab )q. Moreover α 01 α 11 · · · α 11 α 11 α 01 · · · · · · · · · · · · · · · · · · α 11 α 11 · · · α 01 where the identity matrices have dimension M × M , and α 0 = q + (n − 2)q p (1 − q p−1 )(q + (n − 1)q p ) , Doing the matrix multiplication, we find for γ, δ = M + 1, · · · , N − 1. Therefore, correlations arise only between elements H a αβ and H b αγ , where α is a direction in S ⊥ while γ, δ are directions is S. For arbitrary n the conditioning does not induce any non-zero average for these components, since such averages are proportional to the elements ofg a 0 , which are all set to zero. We now come to the n × n blocks M a 1 . From (73), one sees that the conditional covariance matrix of these components is in general a dense matrix, meaning that all components are correlated with each others. Furthermore, the conditioning induces non-zero averages for these components. We denote with µ 1 M|g the n 2 (n+1)/2dimensional vector whose components are the conditional averages M a γδ of the elements in M 1 . We also introduce the (n + 1)-dimensional vectors τ a = (σ a N −n , · · · , σ a N ) collecting the (n + 1) non-zero components of the σ a , and τ 0 = (0, · · · , 1), as well as the n(n + 1)-dimensional vectors χ 1 = τ 1 , · · · , τ n , χ 2 = τ 0 , · · · , τ 0 , and χ 3 = a =1 τ a , · · · , a =n τ a . With this notation, it holds: where the second term arises because of the signal, that induces non-zero averages to the components ofg. This vector can be determined with the same strategy exploited in Sec. V E. In fact, Eq. (78) can be reexpressed in terms of the inverse correlation matrixĈ −1 of the N -dimensional vectors ∇h a , or more precisely of its n(n + 1) × n(n + 1) block associated to the last n + 1 components for each replica. We introduce the n(n + 1) × n(n + 1) rotation matrix: where each blockR a is (n + 1) × (n + 1)-dimensional, with columns given by the non-zero components of the vectors e a M +1 , e a M +2 , · · · , σ a .

Explicit covariances in a given basis
As we have previously remarked, the structure of (77) and (81) is a sole consequence of the separation of the subspaces S and S ⊥ , and it is independent on the choice of the basis in S ⊥ (provided the same choice is made for each tangent plane). Therefore, the correlations and averages depend explicitly only on the last n basis vectors in each tangent plane (having α = M + 1, · · · , N − 1), which are the ones having non-zero projections on the subspace S. We now discuss two possible choices of these basis vectors that strongly simplify the covariances (77) of the single-replica matrixM a . The first possibility is to choose the basis in such a way that: (i) e a N −1 is the projection on the tangent plane at σ a of the special vector w 0 , (ii) e a N −2 is the projection on the tangent plane at σ a of the vector b =a σ b , made orthogonal to e a N −1 , (iii) the remaining n−2 basis vectors are of the form σ b1 + σ b2 + · · · σ b k − kσ b k+1 for non-repeating indices b i = a. For a = 1, this leads to It can be checked that these vectors, together with σ 1 , form an orthonormal basis of the subspace S. Analogous choices can be made for any replica a. Plugging these vectors into (77) with a = b, we find that for any γ = M +1, · · · , N −3 it holds c( =a) (e a γ ·σ c )(e a δ · σ c ) = δ γδ (1 − q), and c( =a) (e a γ · σ c ) = 0. This implies that the components M a αγ for α ≤ M and M + 1 ≤ γ ≤ N − 3 are uncorrelated with each others, and have a modified variance with respect to the one of the larger block, given by: with .
Thus, with this choice of basis in S, the pair of elements in M a 1/2 belonging to the same row and to the last two columns are correlated with each others, and other than that all elements are independent, with variances that depend on the column to which they belongs to. For what concerns the averages (81), we find that in this basis, for γ = M + 1, · · · , N − 3, it holds while with the functions µ ij = µ ij (n, , q, q; r) being a linear combination of the ζ i . For n → 1, one can check that the equality µ 22 (1, , q, q; r) = (1 − q)ζ 1 (1, , q, q; r) holds, while µ 11 (1, , q, q; r) = 0 = µ 12 (1, , q, q; r). Thus, in this limit there are (n − 1) → 0 columns having equal, non-zero average. Similarly, σ 22 → σ γ , while σ 2 11 = p(p − 1) and σ 12 = 0. Therefore, the GOE-structure of the unconditioned matrix is recovered, since the number of columns with different average and variance σ 2 γ is (n−1), and goes to zero. This is the way the annealed limit is recovered. This choice of basis in S is such that there is a unique vector in each tangent plane, e a N −1 , having a non-zero overlap with the direction of the signal w 0 . An alternative choice can be made, for instance to have the mutual in- that diagonalize the matrix (83). This is given by: with z N −2 = (n − 1)(1 − q)(1 + (n − 1)q) and z N −1 = (1 − nq 2 + (n − 1)q)(1 + (n − 1)q). In this new basis: meaning that the components of the last two columns ofM/ √ N are now independent with varianceσ 2 γ and p(p−1), respectively (the fact that the variance along the direction e N −1 is equal to the unconditioned one stems for the fact that e N −1 is orthogonal to any σ a ). Note that in the eigenstates basis, the covariances depend only on the overlap q and are independent on q: this is natural to expect, since the fluctuating part of the Hamiltonian is the p-spin part, that is blind to the special direction w 0 . The information on the signal is carried by the deterministic part, which in the rotated basis reads: 22 + rf k (q) a (2) 22 rf k (q) a 12 rf k (q) a 12 0 , with a (1) 22 2q − (1 − (n − 1)q)[q p + q(p(1 − q) + q)] (n−1)q 2p + q p (1 − p(1 − q)(1 + (n−1)q) − (n−1)q 2 ) − q 2 , a 12 p − 1 = 2A/n q p + (n − 1)q .
In summary, with this second choice of basis vectors in S we find that the decomposition (71) holds, with a deterministic matrix D a equal to 0 0 · · ·· · · 0 0 0 · · ·· · · 0 · · ·· · ·· · ·· · ·· · · · · ·· · ·· · ·· · ·· · · 0 0 · · ·· · · 0 0 0 · · ·· · · 0 · · ·· · ·· · ·· · ·· · · 0 0 · · ·· · · 0 0 0 · · ·· · · 0 0 0 · · ·· · · 0 0 · · ·· · · · · · 0 0 · · ·· · · · · · 0 · · ·· · ·· · · · · · · · · · · ·· · ·· · · · · · · · · 0 · · ·· · · · · · 0 µ γ 0 · · · · · · 0 0 µ γ · · · · · · 0 0 0 µ γ · · · 0 0 · · ·· · ·μ 22μ12 0 · · ·· · ·μ 12μ11 and a stochastic matrix S a having the block structure (72) with (i) m αβ gaussian iid with variance σ 2 = p(p−1), (ii) n ik gaussian iid with variances σ 2 γ ,σ 2 γ and σ 2 for k = M + 1, · · · , N − 3, k = N − 2 and k = N − 1, respectively, (iii) q ij Gaussian, in general mutually correlated. We have not determined their covariances since, as we argue in Appendix IX D, their expression is not necessary to characterize the lowest order corrections to the density of states of the matricesH a . Finally, we point out that with this second choice of basis, the term deriving from the rank-1 perturbation in (70) is no longer diagonal. For later convenience, we define the constants: D. Isolated eigenvalues of the conditioned Hessians In this Appendix, we derive the equations satisfied by the isolated eigenvalues of the conditioned Hessian matrices H, whenever they exist. We focus of the spectrum of the centered matrix A defined bỹ and having itself the block structure where the largest (N − 1 − n) × (N − 1 − n) block A 0 is a GOE with σ 2 = p(p − 1), and A 1/2 and A 1 have the statistics described in Appendix IX C. In particular, we choose the basis in the subspace S to be equal to the second one discussed in the Appendix. In the large-N limit, the bulk of the density of eigenvalues of A is controlled by the largest block A 0 , and is thus a centered semicircle. We aim at determining the poles of the resolvent of (92) that lie on the real axis outside the support of the semicircle, meaning that they are smaller than −2 p(p − 1). From the block structure of A it follows that the trace of (z − A) −1 has two contributions, one coming from the largest (N − 1 − n) × (N − 1 − n) block, and one given by the small n × n block. We focus on this second contribution, since the corresponding matrix elements lie in the subspace S, and have therefore a non-zero overlap with the signal w 0 . The poles of the part of the resolvent coming from this block correspond to isolated eigenvalues having an eigenvector with a nonzero component in the direction of the signal. The quantity to determine are thus the poles of Tr {1/N · D(z)} , where and where now the average is over the distribution of the entries of the matrix A. To compute these poles, we exploit the fact that in the large N limit: This can be shown setting D = D + δD and making use of the expansion: (95) Taking the average of the trace, we find that the corrections to the leading order term in (94) are given by: where the sum is over indices taking n distinct values. The fluctuating part δD of (93) is contributed by two independent terms: the first one is made by the fluctuating components q ij / √ N of the block A 1 , while the second term is made by the fluctuating part of around its mean value. Since the covariances of the q ij are O(1), the first term contributes to the sum (96) with O(1/N ). We now consider the contribution of the second term. For large N : is the resolvent of a GOE matrix with variance σ 2 , while, given the results of Appendix IX C, we have Then: where the constants are of O(1) in N . Thus, the behavior in N of this second contribution is controlled by the decay of the covariances of the matrix elements of the resolvent of a GOE matrix; since the latter go to zero with N as it can be readily checked in perturbation theory, it follows that this is a subleading correction to the leading term in (94). Therefore, in the large-N limit:

Isolated eigenvalue: quenched calculation
We now perform the quenched calculation of the isolated eigenvalue, which accounts for the correlations between minima at fixed, quenched realization of the random Gaussian field. This requires to determine the zeros of Π n (z) in the limit n → 0, which are solutions of the equation: Taking the square of this equation and rearranging the components we find the third order equation with coefficients Note that for k = 1, µ 11 = 0 and thus (113) reduces to a second order equation. For k > 1, the Eq.(113) has three solutions: for the values of parameters that we are considering, we find that at most one out of these three solutions is real and, for some values of parameters, exits the support of the semicircle. We denote this solution with z I . To determine its domain of existence, we require the consistence with the equation (112), i.e., we ask that, when evaluated at z = z I , the RHS of (112) has a sign that is opposite to the sign of the expression −σ 2 z I + σ 2 µ 22 +σ 2 γ µ 11 . This condition has to be imposed separately, since (113) is obtained from (112) by taking a square, and thus it is insensitive to the sign in front of the square root in the definition of the resolvent (99) (it is a generalization to the quenched case of the condition µ < −σ found in the annealed approximation). For those z I that meet this condition, the eigenvalue of the full HessianH/ √ N is given by z I − √ 2u( , q): imposing this expression to be zero gives the critical energy st (q, r) discussed in the main text. The quenched calculation of the isolated eigenvalue gives results that are in general quantitatively different (although qualitatively very similar) with respect to the annealed one. The annealed limit is exact only at the equator q = 0, since in this case q SP ( , 0) = 0 and thus the equation Π n (z) = 0 reduces to (105). In this case, one finds that no eigenvalue exists for k = 1 and k ≥ 3, since the denominator in (108) vanishes, while it exists for r ≥ r c for k = 2.

E. Kac-Rice calculation: additional results
This Appendix contains some additional results related to the content of Sec. V I: we discuss how the mapping (15) is exploited to derive the bands in Figs. 9 and 11, comment on some impliecations on the thermodynamical transitions, and provide some details on the isolated eigenvalue of the Hessian of the stationary points. For k = 1, the band containing the stable stationary points is delimited by the curves q m (r) and q M (r) plotted in Fig. (7). To determine the analogous curves for k = 2 (and fixed r), it is sufficient to consider the functions q → q m (r eff 1 ) and q → q M (r eff 1 ), with r eff 1 (r, q) = rq. Indeed, the latitudes q satisfying q m (r eff 1 ) ≤ q ≤ q M (r eff 1 ) are such that the complexity Σ p,2 ( , q) is positive over a finite energy interval. In Fig. 15, we give an example of this mapping for different values of r: for the smaller r, there is a connected strip containing exponentially many stationary points, which encloses the equator. For the intermediate r, the strip is instead separated into a larger band enclosing to the equator, and a thinner one at larger overlap q (the thinner band at larger overlap has its counterpart at negative overlap). The landscape phase transition in which the band splits into disconnected components occurs between these values of r. The largest r in Fig. 15 corresponds to r > r c ; in this case, the band of states enclosing the equator has shrunk but it is still finite, while the strip closer to the North Pole has collapsed to a single state. The bands at k = 3 can be obtained with an analogous procedure, using r eff 1 = rq 2 . In the same figure we show the case p = 3 and r = r c = √ 6, to illustrate that at the critical point there is a unique connected band containing the equator, with maximal latitude q M (r c ) = q c . We now re-examine the thermodynamic properties of the system at k > 1, and clarify how they can be deduced from the case k = 1. Consider the k = 1 curves * 1 (q, r 1 ) satisfying Σ p,1 ( * 1 (q, r 1 ), q; r 1 ) = 0. For k > 1 and fixed r, it follows from (15) that Σ p,k ( * k (q, r), q; r) = 0 with * k (q, r) = * 1 (q, rq k−1 ) + r(1 − 1/k)q k . This allows one to reconstruct the full spectrum on minima * k (q, r) from  FIG. 17. Comparison between the energies st(q), th (q) and * (q) for k = 3 = p. The orange points are obtained from the solution of the saddle point equations for k = 3, the black squares are obtained exploiting the mapping from k = 1. The intersection points between th and * correspond to the dashed lines in Fig. 11, while the intersection between st and * gives q M . The yellow strip identifies the energies of the stable stationary points. (a) Exactly at r = rc, all states at the latitudes q m < q < q c = 0.707 are stable. The isolated eigenvalue exists only for the states q = q c , which are at their threshold energy (the eigenvalue is attached to the lower edge of the semicircle, that touches zero for the points at this latitude). (b) For r ≈ 2.62 > rc, some of the stationary points in the interval 0.415 q < q M ≈ 0.49 are unstable because of the eigenvalue (the ones at higher energy). The stationary points in the interval q M < q < 0.6 are all unstable irrespective of their energies. For this value of r, the deepest stable minimum (not in the figure) is the isolated one at q * ≈ 0.88.