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

I Introduction

Characterizing rough multi-dimensional energy landscapes is a challenging task that is central in many different fields from physics to computer science, high-dimensional 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 moore ; kurchan ; crisantisommers ; giardinacavagnaparisi , see cavagnapedestrian 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 cavagnapedestrian . 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 footnote1 . Remarkably, it was also realized that pure systems can behave as disordered ones, as first found in long-range spin models bernasconi ; 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 fyodorov the Kac-Rice method has emerged as the mathematical framework suited to do that braydean ; auffingerbenaouscerny ; auffingerbenaous ; fyodorovnadal ; eliran . 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 tengyuma ; MontanariBenArous . 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) reviewkrzakalazdeborova , as well as in biology and in physics, in cases where a specific ground state competes with many random ones (e.g. protein folding remprotein and random pinning glass transition randompinning ). Another important and quite distinct research direction consists in studying the number of equilibria in non-conservative dynamical systems that arise in neuroscience sompolinsky1 and theoretical ecology bunin . 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 touboul ; fyodorovmay . 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 eliran ; in general, the two computations lead to different results even at leading order. Physics methods based on replica theory and super-symmetry provided guidance and results in specific cases, but as we shall discuss in the following they suffer important limitations Monasson ; Annibale ; CLR1 ; CLR2 ; Rizzo1 ; Rizzo2 ; Aspelmeier0 . 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 sherrington1 . 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 cavagnapedestrian ; crisom92 , 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 spiked-tensor model, which recently attracted a lot of attention montanari ; krzakala ; bandera ; Chen ; MontanariBenArous . The thermodynamics of the system we focus on, that we henceforth call generalized spiked-tensor, has been originally introduced in Ref. sherrington1, 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 $\sqrt{N}$, i.e. a given configuration is a vector ${\bf s}$ of $N$ components $\{s_{1},s_{2},\dots,s_{N}\}$ such that $\sum_{i}^{N}s_{i}^{2}=N$. The $N$-dimensional vector ${\bf v_{0}}$ points towards a specific direction, say ${\bf v_{0}}=\{1,1,\dots,1\}$ without loss of generality (we have imposed on ${\bf v_{0}}$ the same normalisation condition as ${\bf 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 cavagnapedestrian with random coupling $J_{i_{1},i_{2},\dots,i_{p}}$ normally distributed with zero mean and variance $\langle J^{2}\rangle=p!/2N^{p-1}$ . The second term represents an energetic gain when the system’s configuration ${\bf s}$ is aligned with ${\bf v_{0}}$. We generically describe this energetic gain by a function $f_{k}(x)$ of the scalar product $x=\sum_{i}^{N}s_{i}v_{0}^{i}/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 in $x=1$, is zero for $x=0$ and is monotonously increasing in $$. It also has a symmetric or an antisymmetric continuation for$x\leq 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. sherrington1, . When$p=k$, the model corresponds to the so called spiked-tensor model which has been the focus of several recent studies in the computer science literature montanari ; krzakala ; bandera ; Chen ; MontanariBenArous . In particular, for this limiting case the calculation of the average number of stationary points has been very recently performed in MontanariBenArous .

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 configurations 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 ${\bf v_{0}}$, and the less so when the overlap with ${\bf 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 ${\bf 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 FyodorovTopologyTrivialization ; fyodledou ; fyodalone and dubbed topological trivialization. The change in the global minima structure, which is directly accessible to a thermodynamic study, was already reported in Ref. sherrington1, . 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^{\prime}(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 crisantisommers ; cavagnagarrahan ; fyodledou . 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 $\overline{q}\in[\overline{q}_{m}(0),\overline{q}_{M}(0)]$, where $\overline{q}_{m}(0)=-\overline{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 $\overline{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 $\overline{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 ${\bf 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_{2\rm{ND}}$, at which a second-order phase transition takes place at the bottom of the landscape, as already found in Ref. sherrington1, . By increasing $r$ above $r_{2\rm{ND}}$ the deepest minima continuously detach from the equator, see the second sphere in Fig. 2 (due to the symmetry $x\rightarrow-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 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 treatmentsherrington1 . 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 $$and have vanishing first and second derivatives in$x=0$. For simplicity, we shall consider even and odd functions under$x\rightarrow-x$when$k$is odd and even respectively. The simplest example of such a function is$f_{k}(x)=x^{k}/k$with$k\geq 3$, first introduced in Ref. sherrington1, . With this choice and taking$p=k$,$H_{p,k}$corresponds to the spiked-tensor model recently investigated in Refs. montanari, ; krzakala, ; bandera, ; Chen, ; MontanariBenArous, . 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^{\prime\prime}(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$.

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_{1\rm{ST}}$, is reached. When $r$ reaches the value $r_{1\rm{ST}}$ 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.

III.3.2 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_{1\rm{ST}}$ and are the ones of the original (unperturbed) $p$-spin model. When $r$ reaches the value $r_{1\rm{ST}}$ 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 sherrington1 ). 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_{1\rm{ST}}$ could take place after that the band has divided and, depending whether $r_{1\rm{ST}}$ 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.

III.4 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 thermodynamic 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.

IV Structure of global minima by the replica method

Using the replica method, one can only partially characterize the energy landscape and its critical points. The aim of this section is to show and recall what kind of information can be gained in this way. The comparison with the Kac-Rice analysis is presented in Sec. VI. Previous studies can be found in Ref. sherrington1, (see also Ref. CrisantiLeuzzi2013, ). Our motivation and perspective on the equilibrium results are different from those, since we focus on where, and to which extent, the recovery of a signal is thermodynamically favored against the noise dispersion. The starting point of the thermodynamic analysis is the evaluation of the free energy $f(\beta,r)$, obtained by computing the $n$-times replicated partition function $\langle Z^{n}\rangle$:

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.

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 $\overline{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 $\overline{q}$. We denote this function $E(\overline{q})$. As it is expected by the $\overline{q}\rightarrow-\overline{q}$ symmetry of the original $p$-spin problem, $E(\overline{q})=E_{GS}+\frac{C_{p}}{2}\overline{q}^{2}+O(\overline{q}^{4})$ for small $\overline{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(\overline{q})-r\frac{\overline{q}^{k}}{k}$.

Case I: If $k=1$ then, no matter how small is $r$, the ground state is at $\overline{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 $\overline{q}=0$ for $r<r_{2{\rm ND}}=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\geq 3$ then a discontinuous transition is bound to take place: for $r<r_{1{\rm ST}}$ the ground state is at $\overline{q}=0$, whereas for $r>r_{1{\rm ST}}$ 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 $\overline{q}$, one can show that at fixed $p$ the solution always remains 1RSB until a given value of $\overline{q}_{c}$ is reached where the 1RSB-RS transition takes place. Moreover the replica structure is the same for identical values of $\overline{q}$. Only the way in which $\overline{q}$ changes as a function of $r$ depends on the value of $k$. In consequence, when the ground state is at $\overline{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 $\overline{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(\overline{q})$ close to the first order transition: one at $\overline{q}=0$ and one at $\overline{q}>0$. The high-overlap minimum can become 1RSB before or after the discontinuous transition depending on the value of $r_{c}$ and $r_{1{\rm ST}}$. The replica analysis that we present below, see also Refs. sherrington1, ; CrisantiLeuzzi2013, , 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_{1{\rm ST}}\leq r\leq r_{c}$.

We present below the replica computation. The following section will be also useful to fully understand the comparison with the Kac-Rice method discussed in Sec. VI. [Readers not interested in replica theory can skip Sec.IV.2 and jump directly to Sec.V]

IV.2 Replica Solution

The standard replica computation cavagnapedestrian for $f(\beta,r)$ leads to the following result

and $Q_{\alpha\beta}$ is an $(n+1)$x$(n+1)$ matrix ($\alpha,\beta\in[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\in[1,n]$ on the remaining $n$x$n$ 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}=\overline{q}$, and the usual RS or $1$RSB 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\neq b$. In the second case the replicas are classified according to $n/m$ different blocks, $Q_{a,b}=q_{1}$ for $a\neq 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 expression of the $1$RSB action $S_{\rm 1RSB}$ in the $N\rightarrow\infty$ and $n\rightarrow 0$ limit is reported in Appendix IX.1 for generic values of $\beta$. When $\beta\rightarrow\infty$ the action reads

where the parameters $\beta(1-q_{1})$, $q_{0}$, $\beta m$, and $\overline{q}$ have to be determined by the following saddle point equations:

with $f^{\prime}_{k}(x)=x^{k-1}$ being the first derivative of $f_{k}$. For each value of $r$, the value of $\overline{q}$ obtained solving the saddle point equations gives the latitude of the deepest minima of the landscape, while the function $-S_{\rm 1RSB}/\beta$ evaluated at the saddle point parameters gives their energy density.

IV.2.2 RS-111RSB transition for the high overlap phase

We now discuss the limit where the ground state energy ceases to be obtained by a $1$RSB 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 crisom92 for small $q_{1}-q_{0}$, and by keeping the lowest order non-zero terms. This gives four equations, see Appendix IX.1. Applying them to the $1$RSB solution with high $\overline{q}$, we get that the critical point occurs at

At this point the high $\overline{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 $\overline{q}$ solution with the solution with $\overline{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.

IV.3 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 $1$RSB $\overline{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-$\overline{q}$ solution which undergoes a continuous transition between a $1$RSB phase and a RS phase at $r_{c}$. The high-$\overline{q}$ solution ($\overline{q}\neq 0$) contains at least partial information about the signal, the amount of this information being represented by the overlap $\overline{q}$. This solution is at first metastable compared to the $\overline{q}=0$ state, but it becomes stable at higher $r$. This occurs through a first order transition at $r_{1{\rm ST}}$. If $r_{\rm 1ST}>r_{c}$ the first order-transition marks a thermodynamic discontinuity between a $1$RSB state (at $\overline{q}=0$) and a RS state (the high $\overline{q}$ one). This scenario is generally found for $k\geq p$ as shown in the phase diagram in Fig. 5. If $r_{\rm 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 $\overline{q}=0$ and high-$\overline{q}$, $1$RSB, states. A continuos transition between the $1$RSB and the RS phase will follow within the high-$\overline{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, $2$nd order-like, transition between the $\overline{q}=0$ state and the high-$\overline{q}$ state before the last one becomes RS at $r_{c}$. As explained before this can be rationalised thinking that the $1$RSB action is quadratic in $\overline{q}$. As such a term $f_{k}$ with higher power of $\overline{q}$ ($k>2$) cannot affect the local stability of the $\overline{q}=0$ state. When $k=2$ instead, $f_{k}$ can counterbalance the quadratic contribution of the $1$RSB action leading to the instability of the $\overline{q}=0$ solution at high enough $r_{2{\rm ND}}$. The 1RSB-RS transition happens for a strictly larger value of $r$ since it takes place for a finite value of $\overline{q}$. Case $k=1$. Finally the case $k=1$ has been extensively studied years ago crisom92 , it corresponds to the $p$-spin spherical model in an external magnetic field. In this case there are no competing $1$RSB states at all. The linear field immediately shifts $\overline{q}$ of the $1$RSB phase away from zero until the continuous transition at $r_{c}$ brings the $1$RSB phase into the RS solution. In order to show concrete examples, we report in Table 1 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\geq 3$ case whether the discontinuous transition to the high $\overline{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.1. We then compare this action to the one of the 1RSB phase with $\overline{q}=0$: $S_{1RSB}^{L}$. There are two possible cases:

$S^{Hc}_{RS}<S^{L}_{1RSB}$. 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 $1$RSB. 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}/\beta=S^{L}_{1RSB}/\beta$ at zero temperature. To the right (respectively left) of the line lie models in which the $1$RSB-RS transition takes place after (respectively before) the discontinuous transition to the high overlap phase. Whenever the $1$RSB-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 $\mathcal{N}_{N}(\epsilon,\overline{q})$ of local minima (or, more generally, of stationary points) of the energy functional, having a given energy density $\epsilon$ and a fixed overlap ${\bf s}\cdot{\bf v_{0}}=N\overline{q}$ with the special direction ${\bf v_{0}}$. The number $\mathcal{N}_{N}(\epsilon,\overline{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 $\mathcal{N}_{N}(\epsilon,\overline{q})$, that is, we aim at computing the quenched complexity $\Sigma_{p,k}(\epsilon,\overline{q};r)$ defined as

As anticipated, we perform the calculation making use of the Kac-Rice formula Kac ; AdlerRandomFields . This formalism has been recently exploited to characterize the topological properties of random landscapes associated to the pure and mixed $p$-spin models auffingerbenaouscerny ; auffingerbenaous , to the spiked-tensor model MontanariBenArous , as well as to count the equilibria of dynamical systems modeling large ecosystems fyodorovmay ; FyodorovTopologyTrivialization and neural networks touboul . 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 $\mathcal{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 crisantisommers ; eliran . 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 $\mathcal{N}_{N}$. The replicated version of the Kac-Rice formula allows us to obtain (to leading order in $N$) the moments $\langle\mathcal{N}_{N}^{n}\rangle$, for integers values of $n$. As we show in the following, the expression for $\langle\mathcal{N}_{N}^{n}\rangle$ that we obtain involves $n$ critical points ${\bf s}^{a}$, $a=1,\cdots,n$, each with energy density $\epsilon$ and overlap $N\overline{q}$ with the North Pole. Introducing their mutual overlaps ${\bf s}^{a}\cdot{\bf s}^{b}=Nq_{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\langle\mathcal{N}_{N}^{n}(\epsilon,\overline{q})\rangle/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}\equiv q$ for $a\neq b$, and take the $n\rightarrow 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.

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:

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. MontanariBenArous . The annealed complexity is an upper bound to the quenched one. As we argue in Sec. V.9, 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\to 0$ has to be performed, by analytically continuing (9). The result is

while $q_{\text{SP}}=q_{\text{SP}}(\epsilon,\overline{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 $\overline{q},\epsilon$. A substantial simplification comes from a general identity that we derive in Sec. V.8 and which relates, for fixed $p$ and $\overline{q}$, the complexities $\Sigma_{p,k}(\epsilon,\overline{q};r)$ for different values of $k$:

for $f_{k}\equiv f_{k}(\overline{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.

V.2 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 ${\bm{\sigma}}={\bf s}/\sqrt{N}$ and ${\bf w_{0}}={\bf v_{0}}/\sqrt{N}$ having unit norm, and we define the rescaled energy functional

with $h_{ps}[{\bm{\sigma}}]\equiv-\sum_{\langle i_{1},i_{2},\dots,i_{p}\rangle}J^{\prime}_{\bf i}\sigma_{i_{1}}\sigma_{i_{2}}\dots\sigma_{i_{p}}$ denoting the $p$-spin energy functional with rescaled coupling satisfying $\langle(J^{\prime}_{\bf i})^{2}\rangle=p!$. We count the stationary point ${\bm{\sigma}}$ of this functional satisfying $h\left[{\bm{\sigma}}\right]=\sqrt{2N}\,\epsilon$ and ${\bm{\sigma}}\cdot{\bf w_{0}}=\overline{q}$, which are in one-to-one correspondence with the stationary points of $H_{p,k}(r)$ with energy density $\epsilon$ and ${\bf s}\cdot{\bf v_{0}}=N\overline{q}$. The Kac-Rice formula incorporates the spherical constraint, as it counts the number of stationary points of the functional $h[{\bm{\sigma}}]$ restricted to the unit sphere; such points ${\bm{\sigma}}$ nullify the surface gradient of (16), which is a vector ${\bf g}[{\bm{\sigma}}]$ lying on the tangent plane to the sphere at the point ${\bm{\sigma}}$. Similarly, their stability is governed by the Hessian on the sphere, which we denote with $\mathcal{H}[{\bm{\sigma}}]$ (see Eq. (19) for a precise definition of this matrix). Given $n$ replicas ${\bm{\sigma}}^{a}$, $a=1,\cdots,n$, we introduce the shorthand notation ${\bf g}^{a}\equiv{\bf g}[{\bm{\sigma}}^{a}]$, $\mathcal{H}^{a}\equiv\mathcal{H}[{\bm{\sigma}}^{a}]$, $h^{a}\equiv h[{\bm{\sigma}}^{a}]$, and denote with $p_{\vec{{\bm{\sigma}}}}$ the joint density function of the $(N-1)n$ gradients components $g_{\alpha}^{a}$ and of the $n$ functionals $h^{a}$, induced by the distribution of the couplings $J^{\prime}_{\bf i}$ in (16). With this notation, the replicated Kac-Rice formula reads:

In (17) the integral is over $n$ replicas ${\bm{\sigma}}^{a}$ constrained to be in the unit sphere, at overlap $\overline{q}$ with the vector ${\bf w_{0}}$. The function $p_{\vec{{\bm{\sigma}}}}({\bf 0},\epsilon)$ is the joint density of gradients and energies evaluated at $g_{\alpha}^{a}=0$ and $h^{a}=\sqrt{2N}\,\epsilon$ for any $a=1,\cdots,n$. The expectation value (18) is over the joint distribution of the Hessians $\mathcal{H}^{a}$, conditioned on each ${\bm{\sigma}}^{a}$ being a stationary point with rescaled energy $\sqrt{2N}\epsilon$, and overlap $\overline{q}$ with ${\bf w_{0}}$. The computation of the moments (17) requires to determine, for each configuration of the replicas ${\bm{\sigma}}^{a}$, the joint distribution of the variables $\mathcal{H}^{a}_{\alpha\beta}$, $g_{\alpha}^{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 ${\bm{\sigma}}$, it can be shown (see the discussion below, and Refs. auffingerbenaouscerny ; fyodorov ) that (i) the gradient ${\bf g}\left[{\bm{\sigma}}\right]$ is statistically independent from $h\left[{\bm{\sigma}}\right]$ and from the Hessian, and (ii) the distributions depend on ${\bm{\sigma}}$ only through its overlap $\overline{q}$ with the special direction ${\bf w_{0}}$ (in absence of the signal, the distribution turns out to be independent on $\bm{\sigma}$). These features make the computation of the annealed complexity feasible; in particular, (ii) is crucial, as it allows to integrate out the variable ${\bm{\sigma}}$ and get an expression for $\langle\mathcal{N}_{N}\rangle$ which depends only on few parameters. Moreover, it suggests that the distributions of the random vector ${\bf g}\left[{\bm{\sigma}}\right]$ and of the random matrix $\mathcal{H}\left[{\bm{\sigma}}\right]$ satisfy some rotation invariant symmetry, hinting at the connection with the random matrix theory of invariant ensembles auffingerbenaouscerny ; fyodorov . 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 $\overline{q}$ and few additional order parameters, that are the overlaps $q_{ab}={\bm{\sigma}}^{a}\cdot{\bm{\sigma}}^{b}$ between the different replicas footnote2 . 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 ${\bm{\sigma}}^{a}$: this is evaluated with standard methods in Sec. V.4. 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.5. 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.6.1 and Sec. V.6.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.6.3. We discuss the stability of the stationary points, which is encoded in the statistics of the spectrum of the Hessians, in Sec. V.7. 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}\equiv q$ for $a\neq b$.

V.3 Structure of covariances and order parameters

As a first step, we analyze the structure of the correlations between the random variables $\mathcal{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 $\overline{q}$ and on the overlaps $q_{ab}={\bm{\sigma}}^{a}\cdot{\bm{\sigma}}^{b}$. To uncover this structure, we consider the gradients ${\bm{\nabla}}h^{a}\equiv{\bm{\nabla}}h[{{\bm{\sigma}}^{a}}]$ and Hessian ${\bm{\nabla}}^{2}h^{a}\equiv{\bm{\nabla}}^{2}h[{\bm{\sigma}}^{a}]$ of the functional (16) extended to the whole $N$-dimensional space footnote3 , and determine the covariances between their components along arbitrary directions in the $N$-dimensional space, given by some $N$-dimensional unit vectors ${\bf e}_{i}$. From here, the correlations of the components ${\bf g}^{a}$ and $\mathcal{H}^{a}$ are easily determined setting ${\bf e}_{i}\to{\bf e}_{\alpha}^{a}$, where $\left\{{\bf e}_{\alpha}^{a}\right\}_{\alpha=1}^{N-1}$ is an arbitrarily chosen basis of the tangent plane at ${\bm{\sigma}}^{a}$. This follows from the fact that ${\bf g}^{a}$ is an $(N-1)$-dimensional vector with components ${g}^{a}_{\alpha}={\bm{\nabla}}h^{a}\cdot{\bf e}_{\alpha}^{a}$, which is obtained from ${\bm{\nabla}}h^{a}$ by simply projecting it onto the tangent plane. Similarly, $\mathcal{H}^{a}$ is an $(N-1)\times(N-1)$ matrix with components

as it follows from imposing the spherical constraint with a Lagrange multiplier footnote4 . For arbitrary ${\bf e}_{i}$, taking the derivative of (16) and computing the expectation value we find:

Finally, the correlations between Hessians and gradients read:

Consider first the case of a single replica: choosing ${\bf e}_{i}\to{\bf e}_{\alpha}\left[{\bm{\sigma}}\right]$ to be vectors in the tangent plane, using (19) and ${\bf e}_{\alpha}\left[{\bm{\sigma}}\right]\cdot{\bm{\sigma}}=0$ one sees that ${\bf g}\left[{\bm{\sigma}}\right]$ is uncorrelated from $\mathcal{H}\left[{\bm{\sigma}}\right]$ and $h\left[{\bm{\sigma}}\right]$; 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 ${\bf e}_{\alpha}^{a}$ in the tangent plane at ${\bm{\sigma}}^{a}$ and the other replicas ${\bm{\sigma}}^{b}$. However, the correlations of the components along directions that are orthogonal to ${\bf w}_{0}$ and to all the ${\bm{\sigma}}^{a}$ hugely simplify. To exploit this, it is convenient to separate the $N$-dimensional space embedding the sphere into the $(n+1)$-dimensional subspace $S$ spanned by the vectors ${\bf w}_{0}$ and $\left\{{\bm{\sigma}}^{a}\right\}_{a=1}^{n}$, and its orthogonal complement $S^{\perp}$. The reference frame of the embedding space, which we denote with $\left\{{\bf x}_{i}\right\}_{i=1}^{N}$, can be chosen in such a way that the last $(n+1)$ vectors ${\bf x}_{N-n}\cdots,{\bf x}_{N}$ are a linear combination of ${\bf w_{0}}$ and of all the ${\bm{\sigma}}^{a}$, forming an orthonormal basis of $S$, while the remaining $N-n-1$ vectors ${\bf x}_{1},\cdots,{\bf x}_{N-n-1}$ generate $S^{\perp}$. Similarly, the basis vectors in the tangent planes ${\bf e}_{\alpha}^{a}$ can be chosen so that the last $n$ vectors ${\bf e}_{N-n}^{a},\cdots,{\bf e}^{a}_{N-1}$, together with the normal direction ${\bm{\sigma}}^{a}$, are a basis for $S$, while the remaining ${\bf e}_{\alpha}^{a}$ with $\alpha<N-n$ generate $S^{\perp}$. In particular, these can be chosen equal for any $a$, as ${\bf e}^{a}_{\alpha}=\delta_{\alpha,i}\,{\bf 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 ${\bf e}^{a}_{\alpha}$, 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.3. Regardless of these choices, Eq. (17) can be rewritten in terms of the overlaps alone, as:

where $\mathcal{E}$ and $p_{\hat{Q}}$ are the expectation value and the joint distribution in Eq.(17), now expressed as a function of the overlap matrix $\hat{Q}$ with components

is an entropic contribution. We determine the leading order term in $N$ of each of the three contributions in (24) for $q_{ab}\equiv q$, and subsequently perform the integral with the saddle point method. To simplify the calculation, we choose the bases ${\bf x}_{i}$ and ${\bf e}_{\alpha}^{a}$ so that only one vector has a non-zero overlap with the special direction ${\bf w_{0}}$: this can be done setting ${\bf x}_{N}={\bf w_{0}}$ (hence the name North Pole), and choosing ${\bf e}_{N-1}^{a}$ to be the projection of ${\bf w_{0}}$ on the tangent plane of ${\bm{\sigma}}^{a}$, ${\bf e}_{N-1}^{a}=\left({\bf w_{0}}-\overline{q}{\bm{\sigma}}^{a}\right)/\sqrt{1-\overline{q}^{2}}$.

V.4 The phase space factor: V​(Q^,q¯)𝑉^𝑄¯𝑞V\left(\hat{Q},\overline{q}\right)

The term $V(\hat{Q},\overline{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 $\hat{\Lambda}$ and $\hat{Q}$ are $n\times n$ matrices in replica space with elements $\Lambda_{ab}=(1+\delta_{ab})\lambda_{ab}$ and $Q_{ab}=(1-\delta_{ab})q_{ab}+\delta_{ab}$, and $v(\hat{Q},\overline{q};{\bm{\mu}},\hat{\Lambda},\vec{\bm{\sigma}})=\sum_{a\leq b}i\lambda_{ab}\left({\bm{\sigma}}^{a}\cdot{\bm{\sigma}}^{b}-q_{ab}\right)+\sum_{a}i\mu_{a}\left({\bm{\sigma}}^{a}\cdot{\bf w}_{0}-\overline{q}\right)$. Performing the Gaussian integrals over the variables $\sigma^{a}_{i}$ and $\mu_{a}$ we get:

This contribution is dominated by $q=\overline{q}^{2}$, which corresponds to configurations in which the replicas are almost independent with each others, correlated only through the constraint on $\overline{q}$ (indeed, it corresponds to replicas having zero mutual overlap in the portion of phase space that is orthogonal to the special direction ${\bf 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}$.

We now determine the joint distribution $p_{\hat{Q}}({\bf 0},\epsilon)$ of the $(N-1)n+n$ components $(g^{a}_{\alpha},h^{a})$. This can be obtained from the joint distribution of the gradient components ${\bm{\nabla}}h^{a}_{i}={\bm{\nabla}}h[{\bm{\sigma}}^{a}]\cdot{\bf x}_{i}$ in the enlarged, $N$-dimensional space, whose covariances read (see Eq. 21):

and averages $\langle{{\bm{\nabla}}h^{a}_{i}}\rangle=-\delta_{iN}\,\sqrt{2N}\,rf^{\prime}_{k}\left(\overline{q}\right)$. The joint density of the ${\bm{\nabla}}h^{a}_{i}$ is thus:

and where $[\hat{C}^{-1}]^{ab}$ is the $ab$ block (in replica space) of the inverse covariance matrix, of dimension $N\times N$. Due to our choice of the reference frame ${\bf x}_{i}$, each $\hat{C}^{ab}$ is block-diagonal, $\hat{C}^{ab}=\text{diag}(\hat{A}^{ab},{\hat{B}}^{ab})$, where $\hat{A}^{ab}$ is an $(N-n-1)\times(N-n-1)$ block with components $\hat{A}_{ij}^{ab}=\delta_{ij}[p\delta_{ab}+p(1-\delta_{ab})q^{p-1}]$ giving the covariances between the gradients components in $S^{\perp}$, while ${\hat{B}^{ab}}$ is an $(n+1)\times(n+1)$ block whose elements are the covariances of the gradients components in $S$, which are explicit functions of $q$ and of $\sigma_{i}^{a}$. To leading order in $N$ this smaller block can be neglected for the computation of the normalization, and one gets:

Note that in (32) the matrix $\hat{C}^{-1}$ is contracted with the vectors ${\bm{\sigma}}^{a}$, so that the quadratic form depends only on the overlaps $q,\overline{q}$. The exponent (32) can be explicitly computed noticing that the vectors ${\bm{\xi}}_{1}$ and ${\bm{\xi}}_{2}$, together with the vector ${\bm{\xi}}_{3}=\left(\sum_{a\neq 1}{\bm{\sigma}}^{a},\cdots,\sum_{a\neq n}{\bm{\sigma}}^{a}\right)$, form a closed set under the action of the matrix $\hat{C}^{-1}$: the inversion of the correlation matrix can be performed in the restricted subspace spanned by these three vectors, and the matrix elements of $\hat{C}^{-1}$ within this subspace suffice to get (32). We refer to the Appendix IX.2 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\to 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) corresponds to the Gaussian weight of the energy functional $h\left[{\bm{\sigma}}^{a}\right]$, while the second accounts for the non-zero average of the last component of the vector ${\bf g}^{a}$ (here we used that ${\bf e}_{N-1}^{a}=\left({\bf x}_{N}-\overline{q}{\bm{\sigma}}^{a}\right)/\sqrt{1-\overline{q}^{2}}$). To leading order in $n$, setting $Q_{p,k}^{(n)}\equiv n\,Q_{p,k}+O(n^{2})$, we obtain

This term is dominated by an energy dependent value of $q=q(\epsilon,\overline{q})$. As we argue in the following section, to leading order in $N$ the expectation value $\mathcal{E}$ turns out to be independent on $q$, so that (36) is the term responsible for shifting footnote5 the saddle-point solution away from the value $q=\overline{q}^{2}$ maximizing the phase space term (27).

V.6 The expectation value of the product of determinants ℰℰ\mathcal{E}

The expectation value $\mathcal{E}$ in (24) is over the joint distribution of the Hessian matrices $\mathcal{H}^{a}$, conditioned on a particular value of the gradients ${\bf g}^{a}$ and field $h^{a}$. Using the identities (19) and (31) we get that the Hessians can be written as:

Consider first a single matrix $\mathcal{M}^{a}$: before conditioning to the values of the gradients and energy functionals, the distribution of each $\mathcal{M}^{a}$ is the one of a GOE matrix, with independent entries with variance $\langle[\mathcal{M}^{a}_{\alpha\beta}]^{2}\rangle=p(p-1)(1+\delta_{\alpha\beta})$, see Eq.(22) (this follows from the fact that the vectors ${\bf e}_{i}^{a}$ in the tangent plane are orthogonal to ${\bm{\sigma}}^{a}$). This distribution is modified by the conditioning, as the entries of $\mathcal{M}^{a}$ are correlated to the gradients and energies of all the other replicas. To determine this effect, we partition each $(N-1)\times(N-1)$ matrix $\mathcal{M}^{a}$ into blocks,

V.6.2 Factorization of the expectation value of the determinants

V.6.3 The GOE computation

Given these observations, and given the symmetry between replicas, the expectation in (40) reduces to

where $I(y)=I(-y)=\int d\mu{\sqrt{2-\mu^{2}}}\log|\mu-y|/\pi$ is given in (V.1), and $\beta(\epsilon,\overline{q})$ in (10). The contribution of the determinants in (24) thus reads:

V.7 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:

where $N$ is here the dimension of the matrix and ${\bm{\sigma}}$ 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 $\overline{q},\epsilon$ at fixed realization of the field, and subsequently average of the random field itself. This leads to

where $\chi(\overline{q},\epsilon)=\delta({\bm{\sigma}}\cdot{\bf w_{0}}-\overline{q})\delta(h[{\bm{\sigma}}]-\sqrt{2N}\epsilon)$ 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\to 0}x^{n-1}$. Exploiting the replicated Kac-Rice formula, we get

Proceeding as before, and using that the resolvent of the Hessian at a stationary point ${\bm{\sigma}}^{1}$ is a function only of the eigenvalue density $\rho^{1}(\lambda)$, we find that (48) can be evaluated with the same saddle-point calculation discussed above, and:

V.8 Mapping between complexity at different k𝑘k

Before presenting the results, we derive the mapping (15) relating the complexity for different values of $k$. Suppose that ${\bm{\sigma}}$ 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}\left[{\bm{\sigma}}\right]$), with overlap $\overline{q}$ and with energy density $\epsilon=\epsilon_{k}$. Then, the point ${\bm{\sigma}}$ is also a stationary point of the functional (16) with $k=1$, provided that $r=r_{1}^{\text{eff}}$ is chosen so that:

In this case, ${\bm{\sigma}}$ has overlap $\overline{q}$ with ${\bf w_{0}}$, and has energy density:

Indeed, for ${\bm{\sigma}}$ to be a stationary point at a fixed $k$, it must hold ${\bm{\nabla}}h_{k,r_{k}}\cdot{\bf e}_{\alpha}\left[{\bm{\sigma}}\right]=0$, which implies:

where we exploited our choice of bases ${\bf x}_{N}={\bf w_{0}}$ and ${\bf e}_{N-1}\left[{\bm{\sigma}}\right]=\left({\bf x}_{N}-\overline{q}{\bm{\sigma}}\right)/\sqrt{1-\overline{q}^{2}}$. Moreover, $h_{\text{ps}}\left[{\bm{\sigma}}\right]=\epsilon_{k}+r_{k}f_{k}(\overline{q})$. This in turn implies that ${\bm{\nabla}}h_{1,r_{1}}\left[{\bm{\sigma}}\right]\cdot{\bf e}_{\alpha}\left[{\bm{\sigma}}\right]=0$ for $\alpha<N-1$, while ${\bm{\nabla}}h_{1,r_{1}}\left[{\bm{\sigma}}\right]\cdot{\bf e}_{N-1}\left[{\bm{\sigma}}\right]=(r_{k}f^{\prime}_{k}(\overline{q})-r_{1}f^{\prime}_{1}(\overline{q}))\sqrt{1-\overline{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.