Landscape Complexity for the Empirical Risk of Generalized Linear Models

We present a method to obtain the average and the typical value of the number of critical points of the empirical risk landscape for generalized linear estimation problems and variants. This represents a substantial extension of previous applications of the Kac-Rice method since it allows to analyze the critical points of high dimensional non-Gaussian random functions. We obtain a rigorous explicit variational formula for the annealed complexity, which is the logarithm of the average number of critical points at fixed value of the empirical risk. This result is simplified, and extended, using the non-rigorous Kac-Rice replicated method from theoretical physics. In this way we find an explicit variational formula for the quenched complexity, which is generally different from its annealed counterpart, and allows to obtain the number of critical points for typical instances up to exponential accuracy.

1 Introduction and main results

Introduction
Characterizing the landscape of the empirical risk is a key issue in several contexts. Many current machine learning problems are both non-convex and high-dimensional. In these cases, the analysis of optimization algorithms, such as gradient descent and its stochastic variants, represents a very hard challenge. In recent years, there has been a series of works that developed a landscape-based approach to tackle this challenge.
The key idea is to study the statistical properties of the empirical risk landscape, and to use these findings to obtain results on the performance of algorithms. Without the aim of being exhaustive this research avenue includes analysis of the landscape of neural networks, matrix completion, tensor factorization and tensor principal component analysis [Fyodorov, 2004, Fyodorov and Nadal, 2012, Kawaguchi, 2016, Soudry and Carmon, 2016, Ge et al., 2016, Freeman and Bruna, 2019, Bhojanapalli et al., 2016, Park et al., 2017, Du et al., 2018, Ge and Ma, 2017, Lu and Kawaguchi, 2017, Ling et al., 2019, Ben Arous et al., 2019, Ros et al., 2019, Sarao Mannelli et al., 2019b,a, Ben Arous et al., 2020, Biroli et al., 2019. The majority of these works identifies the region of parameters where the landscape is "easy", i.e. it focuses on the regime where there shouldn't be any bad local minima and it proves that indeed there are none. However, gradient descent and other landscape-based algorithms are often observed to work even very far from the region described above where the landscape can be proved mathematically to be "easy". A possible reason is that the bounds obtained rigorously are not tight enough. Another, more interesting, is that the landscape is "hard", i.e. spurious minima are present, but their basins of attraction are small and the dynamics is able to avoid them [Sarao Mannelli et al., 2019b,a].
Here we develop a general method that allows to study and focus directly on the "hard" regime, where the empirical risk displays a huge number of bad minima. Our aim is to obtain explicit formulas for the

Main results
We consider two classes of high-dimensional random functions. The first one is a kind of energy that arises in a simple model of neural networks (the perceptron, cf. Engel and Van den Broeck [2001]) and in mean-field glass models [Franz and Parisi, 2016]: where φ : R → R is a smooth activation function (the hypotheses on φ are precised later), x ∈ S n−1 , the unit sphere in n dimensions, and ξ µ are i.i.d. random variables generated from the standard Gaussian distribution in R n . The second class of functions we will consider are related to the loss functions of generalized linear models (GLMs) Nelder and Wedderburn [1972], McCullagh [1984]. In generalized linear estimation an observer has to infer a hidden vector x ∈ S n−1 from the observation of the mdimensional output vector Y = {φ(ξ µ · x )} m µ=1 . In this sense, the GLMs generalize the usual linear regression by allowing the output function to be non-linear 1 . We consider here random GLMs, meaning that the data (or measurement) matrix ξ is taken random, with an i.i.d. standard Gaussian distribution, and we assume that the function φ and the data matrix ξ are given to the observer. This naturally leads to the mean square loss L 2 : GLMs (and their random versions) arise in many different areas of statistics, such as e.g. compressed sensing, phase retrieval, logistic regression, or in random artificial neural networks; we refer to Barbier et al. [2019] for a review of its numerous applications.
Here, we are interested in the statistics of the number of critical points, or more precisely the complexity of the associated empirical risk (2). For any open intervals B ⊆ R + and Q ⊆ (−1, 1), we consider the (random) number Crit n,L 2 (B, Q) of critical points of the function L 2 with loss value in B and overlap with the signal q ≡ x · x in Q: Here grad is the Riemannian gradient on S n−1 . For L 1 we define the similar quantity Crit n,L 1 (B), dropping the notion of overlap. The average E Crit n,L 2 (B, Q) is the quantity that can be analyzed rigorously. Its logarithm divided by n is called the annealed complexity. However, since the random variable Crit n,L 2 (B, Q) is in general strongly fluctuating and scales exponentially with n, its typical value is different from the mean and can be obtained by taking the exponential of E ln Crit n,L 2 (B, Q). This last quantity (divided by n) is called the quenched complexity. It is in general different from the annealed one, with very few exceptions [Subag, 2017, Crisanti andSommers, 1995].
Our main results consist in explicit formulas for the annealed and quenched complexities for L 1 and L 2 . The formula for the annealed case is obtained by a rigorous Kac-Rice method, whereas the one for the quenched complexity is obtained by theoretical physics methods combining the Kac-Rice method with replica theory. We consider the limit n, m → ∞ with m/n → α > 1, a setting called in the statistical physics literature the thermodynamic limit. The condition α > 1 is essential, as can be seen for instance in eq. (1) : if m < n, for each realization of {ξ µ }, the function L 1 has an infinite number of critical points in the set of unit-norm x orthogonal to all the {ξ µ }, and counting the critical points in this case is meaningless (or one has to quotient the space to take care of the degeneracy). Our results hold for many classical activation functions φ, such as e.g. the hyperbolic tangent, the arctangent, the sigmoid, or a smoothed and leaky version of the ReLU activation function 1 . Henceforth we shall denote C + the strict upper-half complex plane, and M(R k ) the set of probability measures on R k . For two probability measures µ and ν we define the relative entropy H(µ|ν) ≡ ln(dµ/dν)dµ if µ is absolutely continuous with respect to ν, and +∞ otherwise. Finally, µ G is a generic notation for the standard Gaussian measure on any R k .
We can now present our main results.
Theorem 1 (The annealed complexity of L 1 ). Let B ⊆ R a non-empty open interval and denote M φ (B) the set of probability measures ν on R such that ν(dt)φ(t) ∈ B. Given: • Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and y ∈ R m a vector with components taken i.i.d. from a probability measure ν. Let D (ν) the diagonal matrix of size m with elements D (ν) µ = φ (y µ ). We define µ α,φ [ν] as the asymptotic spectral measure of zD (ν) z /m.
The precise hypotheses on the activation function φ are precised in Section 2. one has 1 : A note on free probability Interestingly, the measure µ α,φ [ν] can be interpreted as the free multiplicative convolution of the Marchenko-Pastur law (at ratio α) and the asymptotic spectral distribution of the matrix D (ν) , cf. e.g. Voiculescu [1987], Anderson et al.
[2010] 2 . We describe in Section 3 how to explicitly compute the density of µ α,φ [ν], or its linear spectral statistics (as e.g. κ α,φ (ν, C)), via the computation of its Stieljtes transform. We turn to our second annealed result: Theorem 2 (The annealed complexity of L 2 ). Let B ⊆ R + and Q ⊆ (−1, 1) two non-empty open intervals. For q ∈ (−1, 1) we denote M φ (B, q) the set of probability measures ν on R 2 such that Given: • Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and Y ∈ R m×2 with components taken i.i.d. from a probability measure ν on R 2 . Let D (ν,q) the diagonal matrix of size m with elements D (ν,q) µ = f q (Y µ ). We define µ α,φ [q, ν] as the asymptotic spectral measure of zD (ν,q) z /m.
then one has: The proof of Theorem 1 is presented in Section 2. The proof of Theorem 2 is a straightforward generalization, and sketched in Appendix C.1. The variational problems in Theorems 1 and 2 are challenging, as they imply an optimization on a set of measures, and they involve transforms of this measure that are very hard to access numerically. In Section 3 we present a drastic simplification: a heuristic calculation that allows one to reduce the supremum over the probability measure ν to a much more straightforward optimization over a small number of parameters.
As we have already stressed, the annealed complexity, although interesting in itself, is generically not representative of the landscape corresponding to a given typical instance of the empirical risk. In 1 A fully rigorous statement would imply a lower and an upper bound given by a supremum over the adherence and the interior of M φ (B). For reasons of lightness and clarity of the presentation we write it in the simpler presented form.
2 Free multiplication is usually defined for positively-supported measures, however one can generalize it here by explicitly separating the positive and negative parts of φ (we can show freeness of the resulting two random matrices). order to obtain the value of the quenched complexity we use the replicated Kac-Rice method, which is an extension to non-Gaussian functions of the one developed in Ros et al. [2019]. Although the replica method is non-rigorous, it is considered an exact method in theoretical physics and it has been proven to give correct results for spin-glasses and inference problems [Talagrand, 2006, Barbier et al., 2019. We have obtained an explicit formula 1 for the quenched complexity of L 1 and L 2 at fixed values of the empirical risk, and overlap with the solution (in the L 2 case).
For L 1 , using the notations of Theorem 1 we have: Dξ ln I(ξ) .
Here, ξ = (ξ q , ξ a , ξ c , ξ c ) and Dξ is the standard Gaussian probability measure on R 4 . The extr denotes extremization with respect to all variables (A,Â, a,â, C,Ĉ,ĉ, {g(λ)}). We denoted In this formula, the notation extr denotes that one should set the partial derivatives with respect to the involved variables to zero. This notation arises from the replica calculation, which mixes saddle-point computations with Lagrange multipliers associated to certain constraints, and the precise meaning of this extremization (as a supremum or infimum) would have to be clarified by a more rigorous method. On a numerical point of view, one would have to solve the associated saddle-point equations, so that this precise meaning is not crucial for applications. We can state a very similar result for L 2 : Result 2 (Quenched complexity of L 2 ). Let B ⊆ R, Q ⊆ (−1, 1) two open intervals and: • For m ∈ (−1, 1), M φ (B, m) is the space of probability measures ν on R 2 that satisfy: Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and Y ∈ R m×2 with components taken i.i.d. from a probability measure ν on R 2 . Let D (ν) the diagonal matrix of size m with elements D as the asymptotic spectral measure of zD (ν) z /m.
1 Here we used a replica symmetric structure, which is correct in many cases, and a very good approximation in others were replica symmetry has to be broken [Mézard et al., 1987].

Conclusion, outcomes and perspectives
We have obtained analytical results for the annealed and quenched complexities of statistical models with non-Gaussian loss functions arising in generalized linear estimation and simple models of glasses and neural networks. Our method is versatile and can be easily extended to other cases. We describe in Appendix C.2 three other inference models to which it applies: binary linear classification, a mixture of two Gaussians, and a simple model of unsupervised learning.
As a sanity check of our results, we have analytically verified by explicit solution that for a linear activation function, the annealed complexities of L 1 is null. It is again a tedious but straightforward computation to check that the annealed complexity of L 1 with a quadratic activation φ(x) = x 2 is also null, as the number of critical points in this case is linear with n. Note that for L 2 , even the case of a linear activation is non trivial, as shown in the very recent analysis of Fyodorov and Tublin [2020].
Our results allow for a complete characterization of the empirical loss landscapes of generalized linear models. The main issue ahead is determining for which class of functions φ and in which regimes (e.g. values of α), the annealed and quenched complexities become positive, i.e. when the associated landscape is rough. This will allow to study the connection between landscape properties and dynamics induced by local algorithms. In particular, it will shed light on the relationship between the roughness of the empirical loss landscape and the existence of "hard" phases in the learning of generalized linear models [Barbier et al., 2019]. It will also provide an interesting benchmark for obtaining the algorithmic thresholds of gradient descent (and variants) only through the knowledge of the landscape properties [Sarao Mannelli et al., 2019a. Based on ongoing works, we can for instance conjecture the existence of a rough landscape for small enough α in phase retrieval [Lucibello et al., 2019] and retarded learning [Engel and Van den Broeck, 2001].
Addressing these questions requires additional work, which is beyond the scope of this paper. We have explained in Section 3 how to make tractable the variational problem associated to the annealed and quenched complexity. Its analysis for specific models is an ongoing direction of research and will be presented elsewhere.
Another important extension of our results consists in counting the critical points of a fixed index (i.e. with a fixed number of negative directions in the spectrum of the Hessian). This would provide additional interesting information, in particular it would allow to differentiate local minima from the other critical points of the landscape, as was done for spin glass models in Auffinger et al. [2013]. Such a counting would require to understand the large deviations properties of the eigenvalues of the Hessian arising for generalized linear models, a random matrix problem that is hard but hopefully tractable by building on recent developments [Maida, 2007, Ben Arous et al., 2019. This is an ongoing work that we are pursuing, with already promising results.
As final note, it is an open problem to generalize our methods to neural network models with many nodes and hidden layers; the random matrix analysis of the Hessian in this case is a particularly exciting challenge.

Proof of Theorem 1 for the annealed complexity
In this section we prove Theorem 1. The technique leverages the Kac-Rice formula and Sanov's theorem on the large deviations of the empirical measure of i.i.d. variables. First we precise our hypotheses on φ, that we will take in the following set of "well-behaved" activation functions: Definition 1. φ : R → R is "well-behaved" if it is of class C 3 and if, for y ∼ N (0, 1), the random variable a = φ (y) admits a continuous probability density in a neighborhood of a = 0.

The Kac-Rice formula
The first step is to apply the Kac-Rice formula to the random function L 1 : Lemma 3 (Kac-Rice formula). For any x ∈ S n−1 , denote gradL 1 (x) and HessL 1 (x) the (Riemannian) gradient and Hessian of L 1 at the point x. Then gradL 1 (x) has a well defined density (on the tangent space T x S n−1 R n−1 ) in a neighborhood of zero, that we denote ϕ grad L 1 (x) . Denote µ S the usual surface measure on S n−1 . One has: The proof of this lemma uses necessary conditions for a random function to be a.s. Morse 1 , that are stated in Azaïs and Wschebor [2009]. The details are given in Appendix A.1.

The complexity at finite n
In this section we state the result of the Kac-Rice method. For y ∈ R m , let Λ(y) ∈ R m×m : in which we (abusively) denote φ (y) ≡ (φ (y µ )) m µ=1 , and D(y) ∈ R m×m the diagonal matrix with elements D(y) µ = D(y µ ) = n m φ (y µ ). The main result of this section is: Lemma 4 (Complexity at finite n).
in which C n is exponentially trivial, meaning lim n→∞ (1/n) ln C n = 0. The variable y ∈ R m follows N (0, I m ), and z ∈ R (n−1)×m has i.i.d. standard Gaussian matrix elements, independent of y. H Λ n (y) is a square matrix of size (n − 1) with the following distribution : The rest of the section is devoted to the proof of Lemma 4. We start from the result of Lemma 3. The following proposition specifies the joint distribution of (L 1 (x), gradL 1 (x), HessL 1 (x)): Proposition 5 (Distribution of the gradient and Hessian). Let x ∈ S n−1 . Then (L 1 (x), grad L 1 (x), Hess L 1 (x)) follows the following joint distribution: ∼ N (0, I n−1 ), and all {y µ , z ν } are independent. We identified in these equations the tangent spaces T x S n−1 with R n−1 .
Proof. Denote P ⊥ x the orthogonal projection on {x} ⊥ . For a smooth function f : S n−1 → R, ∇f and ∇ 2 f are its Euclidean gradient and Hessian. The Riemannian structure on S n−1 induces the gradient and Hessian of f as grad f (x) = P ⊥ x ∇f and Hess f ( Applying these formulas yields: Letting y µ ≡ ξ µ · x and z µ ≡ P ⊥ x ξ µ (identified to an element of R n−1 ) yields the result.
The joint distribution of eq. (10) is invariant with respect to x, thus we can chose x to be the North pole x = e n = (δ i,n ) n i=1 . With ω n ≡ 2π n/2 /Γ(n/2) the volume of S n−1 , we obtain from Lemma 3: Removing the e n indication and conditioning on the distribution of y in eq. (10), we reach: Once conditioned on y, eq. (10b) describes a Gaussian density so we can directly compute: in which ln C n = O n (n) (using Stirling's formula). The conditioning of the Hessian by grad L 1 = 0 at fixed y reduces to a linear conditioning on z. One thus obtains by classical Gaussian conditioning: in which H Λ n (y) is defined by eq. (9). This ends the proof of Lemma 4.

Large deviations
This section is devoted to the end of the proof of Theorem 1. Denote ν m y ≡ 1 m m µ=1 δ yµ the empirical distribution of y. We take the notations of Theorem 1 and Lemma 4. We first state an important lemma on the concentration of E z | det H Λ n (y)| 1 : Lemma 6. There exists η > 0 such that for all t > 0: The proof of Lemma 6 is detailed in Appendix A.2. Note that we expect this result to actually be valid up to η = 1, as the large deviations of the spectral distribution of random matrices is typically on the n 2 scale [Ben Arous andGuionnet, 1997, Hiai andPetz, 1998]. The following moment condition, proven in Appendix A.3, will be important: Lemma 7. For every γ ∈ (1, α) we have: Then we can conclude using Sanov's large deviation principle [Sanov, 1958, Dembo andZeitouni, 1998], and Varadhan's lemma. Using our previous lemmas, we obtain the statement of Theorem 1: The proof of eq. (18) is detailed in Appendix A.4.
3 Towards a numerical solution to the variational problem is given by the unique solution in C + to the implicit equation (as shown for instance in Silverstein and Bai [1995]): For any µ ∈ M(R) and t ∈ R we define the logarithmic potential as and an arbitrary t ∈ R, see Theorem 1. For clarity, we will write µ for µ α,φ [ν] for the remainder of this section. Let us define for any z ∈ C + , is well defined and holomorphic on C + . Moreover, from the Chapter II of Faraut [2014], we know that From this it is clear that in order to get the logarithmic potential in the bulk, we need to be able to evaluate G(z) for z ∈ C + . Define, for z, g ∈ C + 2 : At any fixed z, F (z, g) is an holomorphic function of g on C + . Its Wirtinger derivative is: Thus g(z) (the Stieltjes transform of µ, cf. eq. (19)) is the only g ∈ C + such that ∂F ∂g (z, g) = 0. Moreover, by definition g(z) is an holomorphic function on C + with values in C + . We can thus apply the usual composition of derivatives and obtain: Furthermore we know dG dz = −g(z). Computing the remaining constant by investigating the limit Re[z] → ∞, we reach that G(z) = F (z, g(z)) for every z ∈ C + . We thus have the crucial relation: This allows for an efficient numerical derivation of the logarithmic potential of µ α,φ [ν], as we will see in more details below.

Heuristic derivation of the simplified fixed point equations corresponding to Theorem 1
We present here an heuristic derivation of scalar fixed point equations for the numerical resolution of Theorem 1. This technique could be easily extended to Theorem 2 as well as the quenched calculations presented afterwards, but we restrict to this simpler case for the sake of the presentation.

Expressing
From eq. (23) we know that for every t ∈ R: For every t ∈ R, g(t + i ) is the only solution in C + to the partial derivative of the previous equation: So heuristically, we can write that for a small enough : with g = g r + ig i (in practice one iterates over g r and g i successively).

Heuristic solution to Theorem 1
We start from the result of Theorem 1. For a function f , we write E[f (X)] ≡ ν(dt)f (t). We introduce Lagrange multipliers to fix the conditions E[φ(X)] ∈ B, and we fix the values of E[φ (X) 2 ] and E[Xφ (X)]. We obtain: Note that now the supremum over ν is unconstrained over the set M(R) of probability distributions. We now make use of eq. (25) to write, with 2K(α) ≡ −1 + ln α − 2α ln α and a small > 0: For any scalar function F , the maximum sup ν [E[F (X)] − αH(ν|µ G )] is attained in ν * with density proportional to e −x 2 /2+F (x)/α , which is called the Gibbs measure in statistical physics. This gives (D is the standard Gaussian measure on R): Plugging this into our previous equation for the annealed complexity yields: This can be further simplified, as the extrema over A, t are trivially solved and give the value of λ 2 = Re[g] and λ 1 = (2A) −1 . Thus we obtain: Let us now denote then the fixed point equations of eq. (30) can be written as: These equations are to be iterated over λ 0 , λ 1 , g, and l (while enforcing the constraint l ∈ B). From experience, the best procedure is to start from the solution of the unconstrained problem (without any constraint on the loss value), before smoothly following the solution while adding the constraint. In the case of L 2 (x), one would follow a similar procedure.

The quenched complexity and the replica method
In this section we detail the principle of the quenched calculation that gives rise to Results 1-2. For the sake of the presentation we restrict to Result 1, while Result 2 will be discussed in Appendix C.1. We therefore focus on the function L 1 of eq. (1). As the very basis of this calculation is non-rigorous we present this calculation in a fashion closer to theoretical physics standards, differently from Section 2 in which we present rigorous results on the annealed complexity. Some technicalities will be postponed to Appendix B.

The replica trick and the p-th moment
The replica method is a heuristic tool of theoretical physics that allows to compute the quenched values of observables in the thermodynamic (i.e. n → ∞) limit from the knowledge of their integer moments, under some assumptions. It is based on the non-rigorous identity (note that it involves an inversion of limits), for a strictly positive function f of a n-dimensional random vector x: Mézard et al. [1987] gives a comprehensive introduction to the replica method and its (many) Here, ϕ {grad L(x a )} p a=1 (0) represents the joint density of the p gradients, taken at 0. Note that the nonlinearity L 1 (x) only depends on the parameters y a µ ≡ ξ µ · x a , so we will often write L 1 (y) ≡ L 1 (x). Proceeding as in the annealed case, we can rewrite the expectations by conditioning over {y a } p a=1 : The gradient and Hessian at x a live in the tangent plane to the sphere at x a , identified with R n−1 . Note that the {y a µ } are Gaussian variables with zero mean and covariance E[y a µ y b ν ] = δ µν q ab , with q ab ≡ x a · x b the "overlap" between replicas a and b. We introduce the variables {q ab } via delta functions in eq. (34): Since we fixed the {q ab }, the distribution of the {y a } is fixed, as well as the joint distribution of the loss, gradients and Hessians, as we will explicit in the following. As the number of overlap variables is p(p − 1)/2 = O n (1), we will perform a saddle-point over the variables {q ab } in the thermodynamic limit. The replica-symmetric assumption (see Mézard et al. [1987]) is a crucial hypothesis that can be made in the framework of the replica method. It amounts to assume that, once the saddle-point is performed, the extremizing {q ab } are "symmetric" over the different replicas of the system. Concretely, we assume that the variables {q ab } satisfy q aa = 1, q ab = q for a = b. Assuming this structure of the overlap matrix allows to extend the expression of the moments to arbitrary non-integer p, and then to take the p → 0 + limit as needed in eq. (32). We used a replica symmetric structure, which is correct in many cases, and a very good approximation in others were replica symmetry has to be broken.

The phase volume factor
Let us first compute the phase space factor in eq. (35). More precisely, the term: in which we denoted q aa = 1. As we detail in Appendix B.1 we reach, when p → 0 + and n → ∞:

The joint density of the gradients
We will now compute the joint density of the gradients at {x a }, conditioned on the values of {y a }. The calculation is an extension of Sections V.C and V.E of Ros et al. [2019]. We consider two vectors x a and x b of overlap q ab = q. It is easy to see that E[grad L(x a )|{y b } p b=1 ] = 0 from eq. (10b), so we will focus on the covariance matrix E[grad L(x a )grad L(x b ) |{y c } p c=1 ]. After some calculations detailed in Appendix B.2 we get the gradient density at leading exponential order: in which the auxiliary functions (z p (q), f p (q), f 0 p (q)) are defined in eq. (86).

Factorization of the mean product of determinants
The argument of this section is very close to Section V.F of Ros et al. [2019]. We consider the term: We make two important remarks, which are straightforward transpositions of the arguments of Ros et al.
[2019] to our problem, and we refer to this work for more extensive physical justifications.
• The conditioning over the gradients being zero, as in the annealed calculation, only gives a finiterank change to the Hessians Hess L 1 (x a ) and thus does not modify the limit at the scale e Θ(n) . At this scale, the statistics of the p matrices {Hess L 1 (x a )} p a=1 are identical.
• The spectral measure of Hess L 1 (x a ) concentrates at a rate at least n 1+ for a small enough > 0 (we expect that the actual rate is n 2 ), so that at the order e Θ(n) the expectation value factorizes and we can assume all the Hessians to be independent.
Before stating the consequences of such remarks, we give some definitions: • µ G,q is the Gaussian probability measure on R p with zero mean and covariance E[X a X b ] = (1 − q)δ ab + q. Note that {y µ } m µ=1 are i.i.d. variables distributed according to µ G,q .
• We define ν y as the empirical measure of (y 1 , · · · , y m ), that is ν y ≡ 1 m µ δ y µ . For every a, we denote ν a y its marginal distribution: ν a y (dλ a ) ≡ b( =a) ν y (dλ). Then ν a y is also the empirical distribution of (y a µ ) m µ=1 . Our remarks show that we can use the results of the annealed calculation, and we have here by factorization of the expectation of the determinants, at leading exponential order: 4.2 Decoupling the replicas and the p → 0 + limit We can then apply Sanov's theorem to the empirical measure ν y ∈ M(R p ). Recall that we have constraints on this measure by the density of the gradient and the fixation of the energy level. More precisely, we denote M (p) φ (q, B) the set of probability measures on R p that satisfy the following: Recall that the functions (z p (q), f 0 p (q), f p (q)) are defined in eq. (86). Leveraging from the results of eqs. (36), (37) and (39), we obtain from Sanov's theorem and Varadhan's lemma: Recall that ν a is the marginal distribution of ν for the variable λ a . We can then decouple the replicas under an assumption on the measure ν that amounts for replica symmetry. We stress that this replica symmetric assumption in the Kac-Rice calculation actually corresponds to a 1-step replica symmetry breaking (1RSB) structure of the zero-temperature Gibbs measure, that is an exponential number of single-point metastable states that all have the same two-point overlap. While possibly not exact, this assumption should already yield a good approximation to the landscape, and could be analytically checked by studying the stability of the replica-symmetric ansatz within replica theory. This allows to take subsequently the p → 0 + limit, and after some simplifications, we reach from eq. (41)  A Technical steps of the proof of Theorem 1

A.1 Proof of Lemma 3
We will apply the Kac-Rice machinery in the form of the remark made in Paragraph 6.1.4 of Azäis and Wschebor [2008]. We recall it as a theorem: Theorem 8 (Azais-Wschebor). Let k, d ∈ N . Let Z : U → R d be a random field, in which U is an open subset of R d . Assume that for every t ∈ U , we can write Z(t) = H [Y (t)], such that: (i) {Y (t), t ∈ U } is a Gaussian random field with values in R k , C 1 paths, and such that for every t ∈ U , the distribution of Y (t) is non-degenerate.
(ii) H : R k → R d is a C 1 function.
(iii) For all t ∈ U , Z(t) has a density ϕ Z(t) (x), which is a continuous function of (t, x) ∈ U × R d .
Define, for every compact set B ⊆ U , N (Z, B) to be the (finite) number of zeros of Z in B. Then: We wish to apply this theorem to the gradient gradL 1 (x). Verifying its hypotheses will end the proof of Lemma 3. We denote ξ ∈ R n×m the matrix {ξ iµ } = {(ξ µ ) i }, ∇L 1 the Euclidean gradient of L 1 , and P ⊥ x the orthogonal projection on T x S n−1 . Since gradL 1 (x) = P ⊥ x ∇L 1 (x) we have: We will apply Theorem 8 with d = n − 1 and k = m × n. The Gaussian random field Y (x) ∈ R n×m is defined as Y (x) ≡ P ⊥ x ξ 1 · · · P ⊥ x ξ m ξ 1 · x · · · ξ m · x . Since Y (x) is just ξ written in an orthonormal basis of R n whose first vector is x, its distributions is non-degenerate. H : R n×m → R n−1 is defined as: Since φ is C 2 , H is C 1 . This verifies (i) and (ii). We turn our attention to verifying (iii). One can write the distribution of the gradient of eq. (43) as grad L 1 (x) ∼ N (0, I n−1 ), and all {y µ , z ν } are independent. Since the distribution of grad L 1 (x) does not depend on x, it is enough to check that its density exists and is a continuous function. To do so, we will show that its characteristic functionφ gradL 1 (x) ∈ L 1 (R n−1 ). We denoteφ a the characteristic function of the random variable a ≡ φ (y), and one obtains: Since α > 1, if [qE zφa (qz)] = O q→∞ (1) we can conclude that φ gradL 1 (x) 1 < ∞. And: 2q 2 e iaz = 1 q E e − q 2 a 2 2 , by Fubini's theorem. Therefore q E zφa (qz) → q→∞ ϕ a (0) by continuity of ϕ a around a = 0 (Definition 1), so φ gradL 1 (x) 1 < ∞. Thus, gradL 1 (x) admits the following probability density: which is a continuous function of u, sinceφ gradL 1 (x) ∈ L 1 (R n−1 ). This shows (iii). In order to show (iv), we will use Proposition 6.5 of Azäis and Wschebor [2008], that we recall here: Lemma 9 (Azais-Wschebor). Let d ∈ N * , and U a compact subset of R d . Consider Z : U → R d a random field, such that (a): The paths of Z are of class C 2 , and (b): There exists C > 0 such that for all t ∈ U and all u in a neighborhood of 0, the density ϕ Z(t) of Z verifies ϕ Z(t) (u) ≤ C. Then P [∃t ∈ U s.t. Z(t) = 0 and det Z (t) = 0] = 0.
Since φ is assumed to be of class C 3 , hypothesis (a) is verified for Z = grad L 1 . Notice then that we can fix C > 0 such that E z∼N (0,1)φa (qz) ≤ C 1+q for all q ≥ 0. Starting from eq. (45): with C c , D n constants depending only on n, using that m ≥ n (α > 1). This shows (b), so by Lemma 9, hypothesis (iv) of Theorem 8 follows. This ends the proof of Lemma 3.

A.2 Proof of Lemma 6
The proof is done in several parts, of which some are inspired by arguments of Silverstein [1995], Silverstein and Bai [1995], Silverstein and Choi [1995], Bai and Silverstein [2010].

A.2.1 Technicalities on the Hessian
We begin by a quick lemma on Λ(y), defined in eq. (8).
Lemma 10 (Low-rank perturbation). Since the distributions of z and y are independent, by rotation invariance we can assume that Λ(y) is a diagonal matrix with elements Λ µ (y). There exists a constant, denoted ||D|| ∞ , such that for all n, y, |D(y)| ≤ ||D|| ∞ . Then we have: (ii) Let Z ∈ R (n−1)×m be i.i.d. variables with zero mean and unit variance. We denote µ Proof. Recall that |D(y)| = (n/m)|φ (y)|. Since m/n → α > 1 and φ is bounded, |D(y)| is bounded (uniformly over n, y) by a constant that we denote ||D|| ∞ . Note that sup 1≤µ≤m |Λ µ (y)| = sup ||u||=1 u Λ(y)u. Using eq. (8) and denoting v(y) ≡ φ (y)/||φ (y)||, we reach in which we used the uniform boundedness of |D(y µ )|, and the Cauchy-Schwarz inequality. This proves (i). We note that 1 n zΛ(y)z and 1 n zD(y)z differ by a rank-2 matrix. (ii) is thus an immediate application of the following lemma (from a course of C. Bordenave): Lemma 11. Let n ≥ 1, and A, B two symmetric matrices of size n, such that the rank of A − B is r. Denote F A (resp. F B ) the c.d.f. of the empirical spectral distribution of A (resp. B). Then This ends the proof of Lemma 10.
Lemma 11. We note λ 1 (A) ≥ · · · ≥ λ n (A) the eigenvalues of A (and similarly for B). Recall weak Weyl's interlacing inequalities [Weyl, 1912]: for every 1 ≤ i ≤ n, λ i+r (A) ≤ λ i (B) ≤ λ i−r (A) (we use the convention λ 1−i = +∞ and λ n+i = −∞ for i ≥ 1). Let t ∈ R, and i, j be the smallest indices such that λ i (A) ≤ t and λ j (B) < t. By the interlacing inequalities, |i − j| ≤ r. And We have some control of the boundedness of the Hessian, summarized in two subsequent lemmas: Lemma 12 (Moment bound). For all γ > 0, one has: Proof. We begin by bounding the extremal eigenvalues λ min , λ max of H Λ n (y) (denoted H for lightness): We used (i) of Lemma 10. Note that this bound is independent of y. In the same way we can bound λ min , and we reach: Using this identity, we have the bound: It is then a very classical result of random matrix theory that the largest eigenvalue of a Wishart matrix zz /n satisfies a large deviation principle in the scale n. This is stated e.g. in Theorem 2.4 of Biroli and Guionnet [2020],which gives moreover the behavior of the rate function I(x). We recall some of its properties: • I(x) : [s max , +∞) → R + is continuous and increasing.
Combining this inequality with eq. (46) ends the proof of Lemma 12.
(ii) The support of µ α,φ [ν m y ] is included in (−C, C) uniformly over y and n.
(iii) For all y ∈ R m , µ α,φ [ν m y ] has a well-defined and continuous density outside x = 0.
Points (ii) and (iii) of Lemma 13 are consequences of Theorem 1.1 of Silverstein and Choi [1995], while item (i) follows from the boundedness of Λ(y) by Lemma 10, and the one of xφ (x).

A.2.2 The cut-off and the logarithmic potential
For any > 0, define ln : x ∈ R + → ln (max(x, )), then x → ln |x| is a −1 -Lipschitz function on R. Let δ ∈ (0, 1). In this section, we show that a cut-off n = n −δ on the smallest eigenvalues does not perturb the logarithmic potential at the thermodynamical scale. We rely on the following result, proven in an ongoing work of Ben Arous, Bourgade and McKenna. For any δ ∈ (0, 1), there exists η > 0 such that for all t > 0: Physically, this makes explicit that, with large probability, there should not be enough eigenvalues of H Λ n (y) around zero so that they contribute macroscopically to the logarithmic potential. This is a consequence the natural fluctuations and repulsion of the eigenvalues of H N n (y). Denote {λ i } n−1 i=1 the (sorted) eigenvalues of H Λ n (y). We can now state: Lemma 14. There exists η > 0 such that for all K > 0: Proof. We consider η given by eq. (47a). Let t > 0. We denote A (n) We have for all y and t > 0 (Ā (n) t being the complementary event to A (n) t ): So that (using ln n (x) ≥ ln(x) for all x > 0): We know ln n |x| ≥ −δ ln(n). By Lemma 12, for all γ > 0: lim sup n→∞ sup y∈R m 1 n ln E z e γ n−1 i=1 ln n |λ i | < +∞.
Fixing γ > 1 and using Hölder's inequality, there exists therefore C > 0 such that for all K > 0 and t ∈ (0, K): in which we used the Markov inequality in (a) and eq. (47a) in (b).
In the following, for simplicity we will often abusively denote ln n | det H Λ n (y)| ≡ e n−1 i=1 ln n |λ i | and ln n E| det H Λ n (y)| ≡ ln E e n i=1 ln n |λ i | .

A.2.3 Concentration of the logarithmic potential with a cut-off
We show here that discarding the eigenvalues of the Hessian that are close to 0 using a cut-off n ≡ n −δ , we have concentration of the logarithmic potential.

Let us recall the Hoffman-Wielandt inequality Hoffman and Wielandt [2003]:
Lemma 16 (Hoffman-Wielandt inequality for the L 2 norm). Let k ∈ N , and A, B ∈ S k (R) be two symmetric matrices with respective eigenvalues λ 1 (A) ≤ · · · ≤ λ k (A) and λ 1 (B) ≤ · · · ≤ λ k (B). Then In particular if A and B are positive matrices one has Tr[AB] ≤ i λ i (A)λ i (B). We use this in eq. (50) along with the n δ -Lipschizity of f n : in which we used Lemma 10. We denote A the event It is a classical concentration result (cf e.g. Chapter 3.1 of Vershynin [2018]) that there exists c > 0 such that: In particular, this implies Let us now show that it suffices to prove the bound of Proposition 15 assuming that A does not occur. Indeed, n −2(1−δ) ln P z [A] ≤ −cn 2δ for a constant c > 0, and: We used Markov's inequality in (a). The inequality (b) can be obtained by very similar arguments than the one used to prove Lemma 12: by rescaling z by z , it is easy to see that the event A will not change the scaling of the large deviations of the largest eigenvalue of the Wishart matrix zz /n, so that the bound of Lemma 12 will also apply when conditioning by the event A. All in all, Using the law of total expectation, we can safely ignore the occurrence of A when showing Proposition 15. Assuming that A is not occurring yields: Recall the Lipschitz concentration of independent variables with laws satisfying the logarithmic Sobolev inequality with a uniform constant c (see for instance Anderson et al. [2010] for a proof and an introduction to the logarithmic Sobolev inequalities): Lemma 17 (Herbst). Let n ∈ N and P be a probability distribution on R n satisfying the Logarithmic Sobolev Inequality with constant c > 0. Let G be a Lipschitz function on R n with Lipschitz constant G L . Then for all t > 0, P[|G − EG| ≥ t] ≤ 2 exp[−t 2 /(2c G 2 L )]. It is easy to check that the Gaussian standard law of z, conditioned by the (extremely probable) event A, satisfies the Logarithmic Sobolev Inequality with constant c = 1 + O n (1). Applying Lemma 17 alongside eq. (54) finishes the proof.

A.2.4 The logarithmic potential of the asymptotic measure
In this part, we relate the expected logarithmic potential to the logarithmic potential of the measure µ α,φ [ν m y ], cf Theorem 1. Proposition 18 (Concentration on κ α,φ ). There exists η > 0 such that for all t > 0: Proof. The proof goes in two parts. First, we show that there exists η 1 > 0 such that 1 : We will then conclude by showing that there exists η 2 > 0 such that for all t > 0: We begin by eq. (57). We take η 2 given by eq. (47b). We have Since ln (x) ≥ ln(x) , we reach from this: and using eq. (47b), we reach eq. (57). Let us show eq. (56). Its proof is based on the following lemma, a consequence of the analysis of Silverstein and Bai [1995], Bai and Silverstein [2010]: Lemma 19. Denote g n (z) the Stieltjes transform of 1 n zΛ(y)z , and g α,φ [ν m y ](z) the one of µ α,φ [ν m y ], for z ∈ C + . Then there exists η ∈ (0, 1) such that for all z ∈ C + : The proof of this lemma is postponed to Appendix A.5. Let us fix η given by this lemma. As stated for instance in Theorem 2.4.4 of Anderson et al. [2010], a consequence of the Stieltjes inversion theorem is that for every Borel set E ⊆ R, we have in which µ n is the empirical spectral distribution of 1 n zΛ(y)z . Fix η 1 < η. We have, uniformly over y: Let us fix C > 0 given by item (i) of Lemma 13. Moreover, we can also bound t φ (ν m y ) by ||xφ || ∞ . This gives that for n large enough the quantity of eq. (60) is bounded (uniformly over y) by: Using eq. (59) and since η 1 < η, this shows eq. (56).

A.2.5 Conclusion of the proof
Let us conclude the proof of Lemma 6 from all our results of Section A.2. We fix δ > 0 such that δ < 1/2. Note that Proposition 15 is a uniform result on y, much stronger than what we required. Since it shows that the concentration (as a function of z) of the log-determinant is in a scale n 1+ with > 0, it implies that there exists η > 0 such that for all t > 0: Combining this identity with Lemma 14 and Proposition 18, we reach the conclusion of Lemma 6.
A.4 Proof of eq. (18) Let t > 0, and fix η > 0 given by Lemma 6. We define E A n is related to the complexity by Lemma 4, and by Varadhan's lemma (and Lemma 7), we have: The factor α in front of the relative entropy arises as we consider the empirical distribution of m = αn i.i.d. variables. For all t > 0, we have by definition of A n , B n : Using Hölder's inequality and Lemma 7, there exists γ > 1 and a constant C > 0 such that: • Assume that lim B n = −∞ and lim sup A n > −∞. Let us fix a lower-bounded sub-sequence A ϕ(n) of A n , so that lim[A ϕ(n) − B ϕ(n) ] = +∞. However, by eq. (67) and Lemma 6, we have n ] → −∞. So we showed that lim B n = −∞ ⇒ lim A n = −∞, which shows eq. (18) in this case.
• Let us now assume that lim B n > −∞. Using the left inequality of eq. (67) and Lemma 6, we reach in the same way that lim inf[A n − B n ] ≥ −t, which implies that lim inf A n > −∞. Thus we can use the right inequality of eq. (67) to show similarly that lim sup[A n − B n ] ≤ t. Taking the t → 0 limit finishes the proof of eq. (18).

B Technical aspects of the quenched calculation B.1 The phase volume factor
Introducing the Fourier transform of the delta functions, we reach at leading exponential order in n: The replica symmetric assumption can be made on the variablesq that achieve this supremum :q aa =q 0 andq ab = −q for a = b. This leads to detq = (q 0 +q) p−1 (q 0 − (p − 1)q), and after taking the p → 0 + limit, we reach: The diverging term − 1 2 log n will be canceled out by the joint density of the gradients as we will see later. The solution of the supremum is easy to carry out, and we finally reach eq. (36).

B.2 The joint density of the gradients
We denote S = Span ({x a } p a=1 ). Following Ros et al. [2019], for every 1 ≤ a ≤ p we can construct an orthonormal basis of S, denoted (e a b ) 1≤b≤p for which x a is the first vector, that is e a a = x a . This basis is convenient, since {x a } ⊥ ∩ S = Span ({e a b } b( =a) ). We can also chose an arbitrary orthonormal basis (e p+1 , · · · , e n ) of S ⊥ . With this choice of basis, we can see that the gradient grad L(x a ) is identified with the vector in R n−1 with components: Recall that ∇L(x a ) = 1 m µ ξ µ φ (y a µ ). Let us make a few remarks: • For every a, the basis (e a b ) p b=1 is only a function of the values of the overlaps {q ab }.
• We consider the joint density of the gradients conditioned by the value of {y a }. In particular, this means that for every a = b, ∇L(x a ) · e a b is fixed by the values of {y c } p c=1 and the overlaps q ab . In particular, the first (p − 1) components of eq. (84) are deterministic, thus their density will yield delta functions that are constraints on {y a } and {q ab }.
• The last n − p components of eq. (84) are (at fixed {y a }) zero mean Gaussian random variables with covariance given by E gradL(x a ) i gradL( . Their joint density taken at 0 is thus at leading exponential order in the n → ∞ limit: Given these remarks and eq. (85), in order to complete the calculation of the joint gradient density we need to compute the quantities (∇ L(x a ) · e a b ) for every a = b as a function of {y a µ } and {q ab }. In order to simplify the calculation, we will already make use of the replica-symmetric assumption on q, that is we assume q aa = 1 and q ab = q for a = b. Let us now describe a possible construction for the basis (e a b ) p b=1 . We introduce three auxiliary quantities that are functions of q and p: With these definitions, we can consider: It is straightforward to check from eq. (87) that we have for all a, b, c that e a b · e a c = δ bc . We can now see that the delta term of the joint density of the gradients taken at 0 is: The product of eq. (85) and eq. (88) gives eq. (37).
B.3 Decoupling the replicas and taking the p → 0 + limit B.3.1 Replica symmetry and decoupling In order to apply the replica method, we need to be able to take the p → 0 limit, by analytically extending eq. (41) to all p > 0. The main idea is that we expect replica symmetry to influence the measure ν that appears in eq. (41). More precisely, we expect that for all permutation π ∈ S p , we have ν(dλ 1 , · · · , dλ p ) = ν(dλ π(1) , · · · , dλ π(p) ). We make in eq. (41) the substitution:  (89) In this last expression, the replica symmetric assumption leads us in particular to assume that µ a = µ for all a. In order to make the remaining calculation tractable we will also need to fix some linear statistics of ν via Lagrange multipliers: • For every a ≤ b, we will fix the linear statistics ν(dλ) φ (λ a )φ (λ b ) = A ab , with Lagrange multipliersÂ ab . Note that by replica symmetry, we actually assume that A ab = a for a = b and A aa = A (and samely for the Lagrange multipliers).
• For all a, b we will fix the linear statistics ν(dλ) φ (λ a )λ b = B ab , with Lagrange multipliersB ab . By replica symmetry, we assume that B aa = B and B ab = b.
Combining these remarks, we reach that the ν-dependent term of eq. (41) is equal to: In the Kac-Rice formula, one has now to integrate over the overlap q ≡ x · x as well. Moreover, we condition over the joint values of a µ ≡ ξ µ · x and b µ ≡ (1 − q 2 ) −1/2 [(ξ µ · x ) − qa µ ], rather than just ξ µ · x (as we did for L 1 ). Note that (a µ , b µ ) follows a joint standard Gaussian distribution.Using these definitions, we can obtain the counterpart of Lemma 4 for L 2 : E Crit n,L2 (B, Q) = C n Q dq e n(1+ln α+ln(1−q 2 )) 2 E a,b δ(P n (a, b))1 L2(a,b)∈B e −nEn(a,b) E z |det H n (a, b)| , in which C n is exponentially trivial, and we defined: Here z ∈ R (n−2)×m is an i.i.d. standard Gaussian matrix. T The condition P q (a, b) = 0 arises from the conditioning on the nullity of the gradient in the linear subspace spanned by x , and E n (a, b) from the density of the gradient in the subspace orthogonal to {x, x }. A crucial feature of this equation is that the joint distribution of (L 2 (x), gradL 2 (x), HessL 2 (x)) only depends on x via the overlap q = x · x with the "true" solution. Once conditioned over the values of q, it thus becomes clear why the calculations made for L 1 will generalize here. As in Section 2.3, one can then show the concentration of the empirical logarithmic potential on the functional κ α,φ (q, ν m a,b ), in which ν m a,b ∈ M(R 2 ) is now the empirical measure of {a µ , b µ } m µ=1 . We obtain the counterpart of Lemma 6: there exists η > 0 such that for all t > 0: lim n→∞ 1 n 1+η ln P 1 n ln E z |det H n (a, b)| − κ α,φ (q, ν m a,b ) ≥ t = −∞.
Thanks to this result, we perform then a saddle point on the overlap q and the empirical measure ν ∈ M(R 2 ), using Sanov's theorem and Varadhan's lemma. This yields the result of Theorem 2. As a final note, there exists similar results to the one presented in Section 3 that allow to compute the density (and the logarithmic potential) of µ α,φ [q, ν], via the computation of its Stieltjes transform.

C.2 Generalizations to other models
Our calculations, both annealed and quenched, can be generalized straightforwardly to many other loss functions and models. As is clear for instance in the annealed computation of Section 2, the key features that must be present are a Gaussian distribution for the data, and a loss function L(x) that only depends on the data samples ξ µ via their projection over a few vectors (as x for L 1 (x) and x, x for L 2 (x)). We give thereafter three examples of models, that can be found in Engel and Van den Broeck [2001], Mei et al. [2018], and for which the calculations can be performed.
Model 2 (Mixture of two Gaussians). Consider n, m ≥ 1 such that m/n → α > 1. We are given m samples y µ ∈ R n , generated as y µ i.i.d.
∼ 2 a=1 p a N (θ 0 a , I n ). The proportions p 1 , p 2 are known, and we wish to recover θ 0 1 and θ 0 2 by minimizing the maximum-likelihood estimator: Model 3 (Simple unsupervised learning model). Consider n, m ≥ 1 such that m/n → α > 1. Let φ : R → R a smooth activation function, V : R → R + a "potential", and x 0 ∈ S n−1 a fixed vector. We assume that we are given i.i.d. data samples {ξ µ } m µ=1 ∈ R n distributed such that their projection on x 0 has a probability density P (ξ µ · x = h) ∝ e − 1 2 h 2 −V (h) , and the other coordinates of ξ µ are i.i.d. standard Gaussian variables. We wish to recover the vector x 0 by minimizing: For each of these three models one can replicate the annealed and quenched calculations of Sections 2 and 4, under suitable technical hypotheses.
A note on non-spherical priors Finally, it is clear from the calculation of Appendix B (particularly Section B.1) that we can also generalize these techniques (at least heuristically) to non-spherical prior distributions on the vectors x. The most natural hypothesis that allows the computation to be generalized is that the prior distribution takes the decoupled form P (dx) = i P (dx i ).