Only tails matter: Average-Case Universal- ity and Robustness in the Convex Regime
Year: 2021

Abstract:
Recent works have studied the average convergence properties of first-order optimization methods on distributions of quadratic problems. The averagecase framework allows a more fine-grained and representative analysis of convergence than usual worst-case results, in exchange for a more precise hypothesis over the data generating process, namely assuming knowledge of the expected spectral distribution (e.s.d) of the random matrix associated with the problem. This work shows that a problem’s asymptotic average complexity is determined by the concentration of eigenvalues near the edges of the e.s.d. We argue that having a priori information on this concentration is a more grounded assumption than complete knowledge of the e.s.d. basing our analysis on the approximate concentration is effectively a middle ground between the coarseness of the worst-case scenario convergence and the restrictive previous average-case analysis. We introduce the Generalized Chebyshev method, asymptotically optimal under a hypothesis on this concentration, and globally optimal when the e.s.d. follows a Beta distribution. We compare its performance to classical optimization algorithms, such as Gradient Descent or Nesterov’s scheme, and we show that, asymptotically, Nesterov’s method is universally nearly optimal in the average-case.

heading: 1 Introduction
text: The analysis of the average complexity of algorithms has a long story in computer science. Average-case complexity, for instance, drives much of the decisions made in cryptography (Bogdanov & Trevisan, 2006).
Despite their relevance, average-case analyses are difficult to extend to other algorithms, partly because of the intrinsic issue of defining a typical distribution over problem instances. Recently though, Pedregosa & Scieur (2020) derived a framework to systemically evaluate the complexity of first-order methods when applied on distributions of quadratic minimization problems. This is done by relating the average-case convergence rate to the expected spectral distribution (e.s.d) of the objective function’s Hessian, which is a well-studied object on random matrix theory. Having access to this object in practice is a much stronger hypothesis when compared to the worst-case analysis that relies only on the values of the edges of this distribution.
Paquette et al. (2020) extended the average-case framework by introducing a noisy generative model for the problems. They further derived the average complexity of the Nesterov Accelerated Method (Nesterov, 2003) on a particular distribution. They showed the strong concentration of the metrics around a limiting value as dimensions go to infinity.
Scieur & Pedregosa (2020) showed that for a strongly convex problem with eigenvalues supported on a contiguous interval, the optimal average-case complexity converges asymptotically to the one given by the Polyak Heavy Ball method (Polyak, 1964) in the worst-case.

heading: 1.1 Current limitations of the average-case analysis
text: When analyzing the state of the art of average-case methods on quadratics problems, we observe significant limitations that we address in this paper. First, little is known
about the convergence rate on convex problems. Also, optimal average-case algorithms require an exact estimation of the e.s.d to guarantee an optimal convergence rate, their convergence rate under inexact e.s.d. is not known. Finally, the non-smooth is also discussed in (Pedregosa & Scieur, 2020), but with little details.
Convex problems. The minimization of non-strongly convex problems is drastically slower than their strongly convex counterpart, as Gradient Descent presents worst-case convergence in Θ( 1t ) and Nesterov is Θ( 1 t2 ). In the strongly convex case, both the worst-case and averagecase are asymptotically equal. However, little is known on optimal average-case rates for convex problems, as well as the average-case complexity of classical methods such as gradient descent or Nesterov’s method, see (Paquette et al., 2020).
Exact estimation of the e.s.d. In (Pedregosa & Scieur, 2020), the theoretical study of optimal algorithms in the average-case requires an exact estimation of the e.s.d. of the problem class. Such estimation may be hard, nor impossible to obtain in practical scenarios. Despite showing good performance when the e.s.d. is estimated with empirical quantities, there are no theoretical guarantees on the performance of the method when the e.s.d. is poorly estimated. There is therefore a need to analyze the algorithm’s performance under different notions of uncertainty on the spectrum. This allows a practitioner to choose the best algorithm for a practical problem, even with imperfect a priori information.
Non-smooth. Pedregosa & Scieur (2020) briefly introduce average-case optimal rates on non-smooth problems, when the e.s.d. is the Laguerre distribution e−λ. In this paper, we extend the analysis to the generalized Laguerre distribution λαe−λ, α > −1.

heading: 1.2 Contributions
text: Our main contribution is a fine-grained analysis of the average-case complexity on convex quadratic problems: we show that a problem’s complexity depends on the concentration of the eigenvalues of e.s.d. around the edges of their support. From this perspective, we propose a family of optimal algorithms in the average-case, analyze their robustness, and finally exhibit a universality result for Nesterov’s method. More precisely,
• (Optimal algorithms). In Section 3, we propose the Generalized Chebyshev Method (GCM, Algorithm 1), a family of algorithms whose parameters depend on the concentration of the e.s.d. around the edges of their support. If the parameters of the GCM method are set properly, the algorithm converges at an optimal average-case rate (Theorem 3 for smooth problems, Theorem 6 for non-smooth problems), a rate that we show is faster than worst-case optimal methods like Nesterov acceleration. We show these rates to be representative of the practical performance of the algorithms in Fig. 6, and retrieve the classical worst-case rates as limits of the average-case (see Table 1).
• (Robustness). Developing an optimal algorithm requires the knowledge of the exact e.s.d. However, in practical scenarios, we only have access to an approximation of the e.s.d. In Theorem 2 in Section 4 we analyze the rate of GCM in the presence of such a mismatch. We also analyze the optimal average-case rates of distributions representing the smooth convex, non-smooth convex, and strongly convex settings and compare them with the worst-case rates (Table 1).
• (Universality). Finally, in Theorem 4, we analyze the asymptotic average-case convergence rate of Nesterov’s method. We show that its convergence rate is nearly optimal (up to a logarithmic factor) under some natural assumptions over the data, namely a concentration of eigenvalues around 0 similar to the Marchenko-Pastur measure. This contributes to the theoretical understanding of the numerical efficiency of Nesterov’s acceleration.

heading: 2 Average-Case Analysis
text: In this section, we recall the average-case analysis framework for random quadratic problems. The main result is Theorem 1, which relates the expected error to the expected spectral
Regime Worst-case Average-Case
Strongly conv. (1 − Θ(1/√κ))t (1 − Θ(1/√κ))t
Smooth conv. 1/t2 1/t2ξ+4 Convex 1/√t 1/tα+2
Table 1: Comparison between function value worstcase and average-case convergence. κ is the condition number in the smooth strongly convex case. In the smooth convex case ξ > −1 is the concentration of eigenvalues around 0 (see Assumption 1) and in the non-smooth case we consider dµ ∝ λαe−λ
distribution and the residual polynomial. The one-to-one correspondence between the residual polynomials and first-order methods applied to quadratics will allow us to pose the problem of finding an optimal method as the best approximation problem in the space of polynomials.
We define a random quadratic problem: Problem 1. Let H ∈ Rd×d be a random symmetric positive-definite matrix independent to x? ∈ Rd, a random vector that is the solution to the problem. We define the random quadratic minimization problem as
min x∈Rd
{ f(x) := 1
2 (x−x?)>H(x−x?)
} . (OPT)
We are interested on minimizing the expected errors E‖f(xt)− f(x?)‖, the expected functionvalue gap, and E||∇f(xt)||2, the expected gradient norm, where xt is the t-th update of a first-order method starting from x0 and E is the expectation over the random variables H,x0 and x?.
The expectation we consider is over the problem and not over any randomness of the algorithm.
In this paper, we consider the class of first-order methods (F.O.M’s) to minimize (OPT). Methods in this class construct the iterates xt as
xt ∈ x0 + span{∇f(x0), . . . ,∇f(xt−1)} . (1) That is, xt belongs to the span of previous gradients. This class of algorithms includes for instance gradient descent and momentum, but not quasi-Newton methods since the preconditioner could allow the iterates to go outside of the span. Furthermore, we will only consider oblivious methods, that is, methods in which the coefficients of the update are known in advance and do not depend on previous updates. This leaves out some methods such as conjugate gradient or methods with line-search.
From First-Order Method to Polynomials. There is an intimate link between firstorder methods and polynomials that simplifies the analysis of quadratic objectives. The next proposition shows that, with this link, we can assign to each optimization method a polynomial that determines its convergence. Following Fischer (1996), we will say a polynomial Pt is residual if Pt(0) = 1. Proposition 1. (Hestenes et al., 1952) Let xt be generated by a first-order method. Then there exists a residual polynomial Pt of degree t, that verifies
xt − x? = Pt(H)(x0 − x?) . (2)
Remark 1. If the first-order method is further a momentum method, i.e.
xt+1 = xt + ht∇f(xt) +mt(xt − xt−1). We can determine the polynomials by the recurrence P0 = 1 and
Pt+1(λ) = Pt(λ) + htλPt(λ) +mt(Pt(λ)− Pt−1(λ)). We note that while most popular F.O.M’s can be posed as a momentum method, the Nesterov method cannot.
A convenient way to collect statistics on the spectrum of a matrix is through its empirical spectral distribution. Definition 1 (Expected spectral distribution (e.s.d)). . Let H be a random matrix with eigenvalues {λ1, . . . , λd}. The empirical spectral distribution of H, called µH , is the probability measure
µH := 1
d
∑d i=1δλi , (3)
where δλi is the Dirac delta, a distribution equal to zero everywhere except at λi and whose integral over the entire real line is equal to one.
Since H is random, the empirical spectral distribution µH is a random variable in the space of measures. Its expectation over H is called the expected spectral distribution and we denote it
µ := EH [µH ] . (4)
We can link the e.s.d. of H to the convergence of a first-order method on the distribution of H. In the following we will consider x0 − x? and H to be independent, with x0 − x? sampled isotropically. Theorem 1. Let xt be generated by a first-order method associated to the polynomial Pt, the measure µ the e.s.d. of H, and E[(x0 − x?)(x0 − x?)T ] = R2I for some constant R. Then we can write the convergence metrics at time step t as
E[‖xt − x?‖2] = R2 ∫ P 2t (λ)dµ(λ), E[f(xt)− f(x?)] = R2 ∫ P 2t (λ)λdµ(λ)
and E[||∇f(xt)||22] = R2 ∫ P 2t (λ)λ 2dµ(λ).
(5)
This shows that polynomials are a powerful abstraction as they allow us to write all of our convergence metrics within the same framework . For simplicity, we set R2 = 1 and we will refer directly to the polynomials associated to a given method. We will refer to objective l as the one associated to the added λl term, i.e. the function-value is objective l = 1.
This framework is linked to the field of orthogonal polynomials by the next proposition. We construct an optimal method w.r.t. a given distribution through a family of orthogonal polynomials associated to it. Proposition 2 ((Pedregosa & Scieur, 2020)). Let P lt be defined as
P lt := arg minPt(0)=1 ∫ P 2t (λ)λ ldν(λ). (6)
Then (P lt ) is the family of residual orthogonal polynomials w.r.t. to λl+1dν.
This theorem further implies that the optimal first-order method is a momentum method as Favard’s theorem Marcellán & Álvarez-Nodarse (2001) tells us the orthogonal polynomials w.r.t. a given distribution are related through a three term recurrence,
Pt+1(λ) = atPt(λ) + btλPt(λ) + (1− at)Pt−1(λ). (7) Following Remark 1, the optimal method is derived from this recurrence as
xt+1 = xt + (at − 1)(xt − xt−1) + bt∇f(xt) . (8)

heading: 3 Methods
text: Being able to write the rates in terms of the expected spectral distribution ties the average-case framework to the field of random matrix theory. Indeed, because of results from this field, certain e.s.d’s are considered more natural than others. Indeed, it can be shown that the same distribution arises when we take the gram matrix of random centered i.i.d. features with variance σ2: the Marchenko Pastur distribution. Definition 2 (MP distribution). The Marchenko Pastur distribution associated with the parameter r and with scale σ2 is given by
dµMP (λ) = 1
2πσ2
√ (λ+ − λ)(λ− λ−)
rλ , (9)
with λ+ = σ2(1 + √ r)2, λ− = σ2 max(0, (1−√r)2).
The Marchenko Pastur distribution µMP can be considered a natural first model for e.s.d’s as it arises universally from matrices with i.i.d. entries,under mild low moment assumptions, there is no specific distribution of the matrix to be considered. It can be seen as a model for the white-noise in the data. When r = 1, i.e. n = d, we have dµMP ∝ λ−1/2 √ λ+ − λ.
Pedregosa & Scieur (2020) first derived the optimal method w.r.t. µMP , and Paquette et al. (2020) derived Nesterov’s rates under the distribution. As we are concerned with being robust, a natural step is to consider the Beta weights. Definition 3. The (generalized) Beta weights with parameters τ, ξ and scale L are given by the (non-normalized) pdf
dµ(λ) = λξ(L− λ)τ . (10)
This family of distribution generalizes the MP distribution, and both have similar concentrations near 0 when ξ ≈ −1/2. The optimal method w.r.t. µ and objective l is associated to a shifted Jacobi polynomial P̃α,βt with β = ξ + l + 1, α = τ . When α = β = −1/2, we retrieve the Chebyshev Method (Flanders & Shortley, 1950). As such, we name our proposed methods the Generalized Chebyshev Method (GCM).
Algorithm 1: GCM(α, β) Inputs: Initial vector x0, function f , smoothness parameter estimate L x−1 ← 0, δ0 ← 0 for t = 1, . . . , T do
at ← − 2(β2+αβ+(2t+1)(α+β)+2t2+2)(2t+α+β+1)
2(t+1)(t+α+β+1)(2t+α+β)
bt ← (2t+α+β+1)(2t+α+β+2)L(t+1)(t+α+β+1) γt ← − (t+α)(t+β)(2t+α+β+2)(t+1)(t+α+β+1)(2t+α+β) δt ← 1at+γtδt−1 xt ← xt−1 + (δtat − 1)(xt−1 − xt−2) + δtbt∇f(xt−1)
We’ll consider the Nesterov’s method used in Paquette et al. (2020), which is defined by the iterations:
xt+1 = yt − 1
L ∇f(yt) (11)
yt+1 = xt+1 + t
t+ 3 (xt+1 − xt) (12)
We also consider the Laguerre method, which is optimal w.r.t. dµ(x) = x αe−x
Γ(α+1) , taking α as a parameter. This method is proposed to optimize non-smooth functions. Both these methods are generalizations of ones that have been proposed in Pedregosa &
Scieur (2020). We show that Algorithm 1 corresponded to polynomials P̃α,βt and derive the Laguerre method in appendix B.
Remark 2. The Generalized Chebyshev takes the largest eigenvalue L as a parameter, but the rates we will show are robust to an overestimation of L.

heading: 4 Robust Average-Case Rates
text: We will state our assumption over the spectral distributions. It effectively allows us to parametrize all of our distributions of interest in a way that characterizes the asymptotic convergence, diving them into equivalence classes. Assumption 1. We will write ντ,ξ for a continuous distribution supported in (0, L] s.t. ν′τ,ξ(x) > 0 for x ∈ [0, L], dντ,ξ = Θ(λξ) near 0 and dντ,ξ = Θ((L− λ)τ ) near L.
Assumption 1 is quite nonrestrictive, in that, the spectral distribution of the Hessian for any smooth convex problem can be identified with some τ, ξ in this class. It is a milder assumption than (1) assuming complete knowledge of the spectrum of the Hessian or (2) the specific distribution on the entries of your data. Moreover the assumption encompasses the frequently used MP (e.g. Martin & Mahoney (2021); Pennington & Bahri (2017)) and Uniform distributions We note there’s no need to consider eigenvalues situated at 0 as they do not contribute to the optimization process.
The ξ works as a measure of how close we are to the worst-case scenario, as it approaches −1. Samples in finite dimension of distributions with high values of ξ will work as strongly convex functions in practice.
We show that ντ,ξ indeed behaves like an equivalence class when considering the asymptotics of the convergence of the methods: only the concentrations near the edge matter. We do this by singling out from each of these classes the beta distributions for which we can compute the rates, then show the rates to be the same inside ντ,ξ. Theorem 2 (GCM average-case rates). A Generalized Chebyshev Method with parameters (α, β) applied to a problem with e.s.d. as in Assumption 1 has average-case rates
E[f(xt)− f(x?)] ∼ L · Cα,β1,ν  t −1−2β if α < τ + 1/2 and β < ξ + 3/2 t−2(ξ+2) log t if α = τ + 1/2 and β = ξ + 3/2 t2(max{α−β−τ,−ξ−1}−1) if α > τ + 1/2 or β > ξ + 3/2 ,
(13)
E[||∇f(xt)||22] ∼ L2 · Cα,β2,ν  t −1−2β if α < τ + 1/2 and β < ξ + 5/2 t−2(ξ+3) log t if α = τ + 1/2 and β = ξ + 5/2 t2(max{α−β−τ,−ξ−2}−1) if α > τ + 1/2 or β > ξ + 5/2 ,
(14)
where Cα,βν is a distribution dependent constant.
Theorem 2, which is illustrated by fig. 2 shows that overestimating β, and underestimating α will still leave us with the optimal asymptotic rates, so a good rule of thumb for calibrating the algorithm is to use high β and low α.
Theorem 3 shows that a proper choice of α, β can indeed make the Jacobi polynomial asymptotically optimal w.r.t. to any ντ,ξ. Theorem 3 (Optimal Rates). Let ν follow Assumption 1. The optimal asymptotic averagecase rates for E[f(xt)− f(x?)] and E[||∇f(xt)||22] are attained by the GCM with parameters (τ, ξ + 2) and (τ, ξ + 3), respectively, and read
E[f(xt)− f(x?)] = Θ(t−2(ξ+2)), E[||∇f(xt)||22] = Θ(t−2(ξ+3)).
For the function value (l = 1), we find rates that approach t−2 as ξ → −1, showing the worst-case as a limit (over the considered distribution) on the average-case.
We remark that the above theorems imply that, at least asymptotically, the GCM is robust for a suboptimal choice of parameter β up to 1/2 below the optimal choice and infinitely above.
For completeness, we also derive worst-case rates for the GCM: Proposition 3 (GCM worst-case rates). Let f be a convex, L-smooth quadratic function. Then, for the Generalized Chebyshev Method with parameters (α, β), we have worst-case rates
f(xt)− f(x?) ≤ C1L  t 2(α−β) if α > β − 1 t−1−2β , if α ≤ β − 1 β ≤ 12 t−2, if α ≤ β − 1 β ≥ 12 , (15)
||∇f(xt)− f(x?)|| ≤ C2L2  t 2(α−β) if α > β − 2 t−1−2β , if α ≤ β − 2 β ≤ 3/2 t−4, if α ≤ β − 2 β ≥ 3/2 . (16)
For a reasonable choice of α, β, i.e. β ≥ 12 , α ≤ β − 1. the function value achieves the theoretical lower bound of t−2.
We now analyze the convergence of the Nesterov method. Nesterov (2003) has shown that it matches up to a constant factor a lower bound on the worst-case complexity of non strongly convex problems. A natural question is if this performance would translate to good average-case rates. To do so, we will extend Paquette et al. (2020) proof for the Nesterov method rates under the MP distribution. Theorem 4 (Nesterov average-case rates). Let ν as in Assumption 1. Then for the Nesterov method, we have average-case rates
E[f(xt)− f(x?)] ∼ C ′1,ν { t−2(ξ+2) if ξ < −1/2 t−3 log t if ξ = −1/2 t−(ξ+7/2) if ξ > −1/2 , E[||∇f(xt)||22] ∼ C ′2,νt−(ξ+9/2). (17)
The difference between the asymptotic average-case rates of Nesterov and the optimal ones are tξ+l−1/2, when ξ + l > 1/2, log t when ξ + l = 1/2 and 0 otherwise. This shows that Nesterov is almost optimal when the concentrations near 0 are relatively high, i.e. low ξ. Theorem 5 (Gradient Descent average-case rates). Let ν as in Assumption 1. Then for gradient descent
E[f(xt)− f(x?)] = Θ(t−(ξ+2)), E[||∇f(xt)||22] = Θ(t−(ξ+3)). (18)
Observe for the function value that the rate for Nesterov is t−2 and the rate for Gradient Descent is t−1 when ξ → −1. Lastly, we consider the optimal rates for a Gamma distribution. Theorem 6 (Laguerre method rates). Let α > −1 and µα, α > −1 be a Gamma distribution, i.e. dµα(x) = x αe−x
Γ(α+1) . The optimal rates are given by the Laguerre method of appropriate tuning and
E[f(xt)− f(x?)] = Θ(t−(α+2)). (19)
Note that this result does not have the same universality as the others because of the non-compacity of the distribution’s support. These rates are contrasted to the worst-case lower bound on the optimization of non-smooth functions by first-order methods, which gives
f(xk)− f(x?) ≥ C√ t .
These rates are not found when α → −1, indicating that the worst-case is especially pessimistic in this scenario. Remark 3. All of the expected rates we state are almost deterministic on the high dimensional setting as per the concentration results shown in Paquette et al. (2020)

heading: 5 Experiments
text: We simulate the e.s.d’s in two ways. The Marchenko Pastur distribution, which we sample by taking H = XXT where X has i.i.d. gaussian samples. This enables us to simulate (τ, ξ) values of (1/2,−1/2). Other values of (τ, ξ) are simulated by sampling Λ ∈ Rd from the corresponding Beta distribution and taking H = U diag(Λ)UT , where U is an independently sampled orthonormal matrix.
We let x? = 0 and sample x0 from a centered Gaussian distribution, the dynamics are the same as in the general case. In all experiments we use the problem’s instance largest eigenvalue to calibrate each method (e.g. Gradient Descent’s stepsize is 1/L). Our theoretical rates in Theorem 5 and Theorem 4 respectively for the Nesterov method and Gradient Descent are precise under the approximate range −1 < ξ < 0 as we show in Figure 6. Distributions with higher ξ need many samples otherwise they behave as strongly convex functions.
The same is not true for the Generalized Chebyshev Method. If β < β? or ξ is low the empirical findings diverge from the theoretical. We believe this is due to numerical instability
under these regimes as the metrics also have much larger variance than in the other regimes. We’ve not been able though to pinpoint the exact source of this supposed instability. This is shown in appendix D.
The GCM with β > β? performs corresponding to the theory, and it’s non-asymptotically very close to the performance of β?. High values of β also perform very well on non-synthetic data, suggesting in practice we should use these values.

heading: 6 Conclusion and further work
text: In this paper, we’ve established that the assymptotic convergence of first order methods on quadratic problems in the convex regime depend on the concentration of the Hessian’s eigenvalue near the edges of the spectrum’s support. We further contributed to the theoretic understanding of the Nesterov’s method performance and established the contrast between the worst-case and average-case in the main regimes considered in Optimization.

heading: A Proofs of Section 2
text: Theorem 1. Let xt be generated by a first-order method associated to the polynomial Pt, the measure µ the e.s.d. of H, and E[(x0 − x?)(x0 − x?)T ] = R2I for some constant R. Then we can write the convergence metrics at time step t as
E[‖xt − x?‖2] = R2 ∫ P 2t (λ)dµ(λ), E[f(xt)− f(x?)] = R2 ∫ P 2t (λ)λdµ(λ)
and E[||∇f(xt)||22] = R2 ∫ P 2t (λ)λ 2dµ(λ).
(5)
Proof. We remark that by the definition of the expected spectral distribution µ of H, we have for continuous g
EH [g(tr(H))] = ∫ g(λ) dµ(λ) (20)
We know that xt − x? = Pt(H)(x0 − x?). We can write ||xt − x?||2 in terms of a trace and use the independence of H and x0 − x? to connect it to the e.s.d.:
E||xt − x?||2 = E[tr((x0 − x?)TPt(H)2(x0 − x?))] (21) = EH,x0−x? [tr(Pt(H)2(x0 − x?)(x0 − x?)T ] (22) = EH [ Pt(H) 2Ex0−x? [(x0 − x?)(x0 − x?)T ]) ] (23)
= R2EH [Pt(tr(H))2] = R2 ∫ Pt(λ) 2 dµ(λ) (24)
For the gradient and function value the reasoning is the same by noticing that
E[f(xt)− f(x?)] = E[tr((x0 − x?)TPt(H)HPt(H)(x0 − x?))] (25) = EH [(λPt)(tr(H))2], (26)
where λPt is also a polynomial. As ∇f(xt) = H(xt − x?).
E||∇f(xt))||2 = E[tr((x0 − x?)TPt(H)H2Pt(H)(x0 − x?))] (27) = EH [(λ2Pt)(tr(H))2] (28)
Proposition 2 ((Pedregosa & Scieur, 2020)). Let P lt be defined as
P lt := arg minPt(0)=1 ∫ P 2t (λ)λ ldν(λ). (6)
Then (P lt ) is the family of residual orthogonal polynomials w.r.t. to λl+1dν.
Proof. We differentiate the expression for the metrics w.r.t. to the coefficients of the polynomials:
d
dak
(∫ λlP 2t (λ)dµ(λ) ) = ∫ λl d
dak
( t∑
k=0
akλ kPt(λ) ) dµ(λ) =
= 2 · (∫ λl+kPt(λ)dµ(λ) ) = 0
This means that Pt(λ) is orthogonal to any polynomial of degree t− 1 w.r.t to the intern product 〈., .〉λl+1dµ

heading: B GCM and Laguerre method derivation
text: We will first state two lemmas that allow us to construct the optimal polynomials. With them in hand the procedure is trivial.
Lemma 7. Let (P̃t) be a family polynomials following
P̃t(λ) = (αt + βtλ)P̃t−1)λ) + γtP̃t−2(λ),
with P̃0 a constant polynomial and P̃t 6= 0,∀t. Then Pt(λ) = (at + btλ)Pt−1(λ) + (1− at)Pt−2(λ) (29)
is the recurrence for Pt(λ) = P̃t(λ)/P̃t(0). With:
at = δtαt (30) bt = δtβt (31) δt = (αt + γtδt−1) (δ0 = 0) (32)
The proof of this is presented in Pedregosa & Scieur (2020). Further, we know how to compute the recurrence for the polynomials of a shifted distribution:
Lemma 8. Let (P̃t) be a family polynomials orthogonal w.r.t following
P̃t(λ) = (αt + βtλ)P̃t−1(λ) + γtP̃t−2(λ), (33)
and define polynomials Pt s.t. : Pt(m(λ)) = P̃t(λ),
with m(λ) = aλ+ b a non singular affine transform. Then Pt follows a recurrence like in eq. (33), with:
α′t = αt + bβt (34) β′t = aβt (35) γ′t = γt (36)
The lemma is self-evident by considering eq. (33) with argument m−1(λ)
These results are enough to get the recurrence relation for the residual polynomial w.r.t xβ(L − x)α. We begin by the standard Jacobi polynomials, which are orthogonal w.r.t (1− x)α(1 + x)β and follow a recurrence according to αt, βt, γt below Szegö (1975):
αt = (2n+ α+ β)(2n+ α+ β − 1)
2n(n+ α+ β) (37)
βt = (α2 − β2)(2n+ α+ β − 1)
2n(n+ α+ β)(2n+ α+ β − 2) (38)
γt = −2(n+ α− 1)(n+ β − 1)(2n+ α+ β)
2n(n+ α+ β)(2n+ α+ β − 2) (39)
We then shift the distribution according to η(x), and then transform to the residual ones. We slightly simplify these computations and use remark 1 to get Algorithm 1.
We know (Szegö, 1975) that the Laguerre polynomials Lαt , with usual normalization, follow the recurrence
Lαt (λ) =
( 2t+ α− 1
t − 1 t λ
) Lαt−1(λ) +
t+ α− 1 t Lαt−2(λ) (40)
As we don’t have to shift the domains, we have only to apply lemma 7 to get the Laguerre method.Further, we can get a explicit expression for δt = tt+α , simplifying the expression.
Algorithm 2: Laguerre(α) Inputs: Initial vector x0, function f x−1 ← 0 for t = 1, . . . , T do
xt ← xt−1 + t−1t+α (xt−1 − xt−2)− 1t+α∇f(xt−1)

heading: C Proofs of section 3
text: In the following we will consider shifted versions of the spectral distributions. This shift is written as an affine transform m(λ) : [0, L]→ [−1, 1] because most results in the theory of orthogonal polynomials are stated in terms of distributions supported in [−1, 1]. This can be seen as an additional layer of abstraction because the quantities evaluated with the shifted distributions and polynomials are proportional, i.e. if Pt(x) = P̃t(m(x)) and µ′(x) = µ̃′(m(x)): ∫
P 2t (x)µ ′(x) dx ∝ ∫ P̃ 2t (x)µ̃ ′(x) dx (41)
So all the asymptotics are the same. The Jacobi polynomials Pα,βt are those orthogonal w.r.t dµ(x) = (1− x)α(1 + x)β . Most works use the normalization P̃α,βt (−1) = (−1)t ( t+β t ) . We will write P̃α,βt for this normalization and P α,β t for the residual polynomials
Theorem 2 (GCM average-case rates). A Generalized Chebyshev Method with parameters (α, β) applied to a problem with e.s.d. as in Assumption 1 has average-case rates
E[f(xt)− f(x?)] ∼ L · Cα,β1,ν  t −1−2β if α < τ + 1/2 and β < ξ + 3/2 t−2(ξ+2) log t if α = τ + 1/2 and β = ξ + 3/2 t2(max{α−β−τ,−ξ−1}−1) if α > τ + 1/2 or β > ξ + 3/2 ,
(13)
E[||∇f(xt)||22] ∼ L2 · Cα,β2,ν  t −1−2β if α < τ + 1/2 and β < ξ + 5/2 t−2(ξ+3) log t if α = τ + 1/2 and β = ξ + 5/2 t2(max{α−β−τ,−ξ−2}−1) if α > τ + 1/2 or β > ξ + 5/2 ,
(14)
where Cα,βν is a distribution dependent constant.
Proof. We will prove that for any α and β, ξ, τ > −1, l > 0 and ν following Assumption 1, we have∫ Pα,βt (x) 2xldντ,ξ−l(x) ∼ LlCα,βν  t −1−2β if α < τ + 1/2 and β < ξ + 1/2 t−2(ξ+1) log t if α = τ + 1/2 and β = ξ + 1/2 t2(max{α−β−τ,−ξ}−1) if α > τ + 1/2 or β > ξ + 1/2
We will first show this result for the Beta weights, then show that distributions with the same concentration behave similarly. The normalization of P̃α,βt is s.t.[ Szegö (1975) (4.3.3)]:∫ 1 −1 P̃α,βt (x)(1−x)α(1+x)β dx = 2α+β+1 2n+ α+ β + 1 Γ(n+ α+ 1)Γ(n+ β + 1) Γ(n+ 1)Γ(n+ α+ β + 1) = Θ(t−1) (42)
Further, the residual polynomials are s.t. |Pα,βt | = Θ(t−β)|P̃α,βt |, from the definition of the classical normalization. We state the result (Exercise 91, Generalisation of 7.34.1) from Szegö (1975):
Lemma 9. We have ∫ 1 0 (1− x)τPα,βt (x)2dx ∼ Θ(hατ ) (43)
hατ :=  t 2(α−τ−1) if α > τ + 1/2 t−1 log t if α = τ + 1/2 t−1 if α < τ + 1/2
(44)
Noting that P̃α,βt (x) = (−1)tP̃ β,αt (−x), we can write:∫ 1 −1 P̃t(x) 2(1−x)τ (1+x)ξdx = Θ (∫ 1 0 (1− x)τ |P̃α,βt (x)|2dx ) +Θ (∫ 1 0 (1− x)ξ|P̃ β,αt (x)|2dx ) (45) We can then show our result for dντ,ξ−l(x) = xξ−l(L− x)α by carefully considering each of the cases on Lemma 9 and the maximum of each term in eq. 45, and an added t−2β from the different normalization. With this, we have the wanted result for the Beta weights It remains to show:∫ 1
0
P̃α,βt (x) 2 dντ,ξ(x) = Θ (∫ 1 0 (1− x)τ P̃α,βt (x)2dx )
(46)
And the rest follows from the same arguments. We do this with the help of this lemma shown in Van Assche (1995) relating to the weak convergence of the orthogonal polynomials:
Lemma 10. Let µ be a measure and (pt) it‘s family of orthonormal polynomials s.t. pt follow the recurrence:
xpt(x) = atpt+1(x) + btpt(x) + at−1pt−1(x)
and at, bt converge respectively to a, b. Then for any f continuous and bounded:∫ f(x)p2t (x)dµ(x)→ 1
π ∫ 1 −1 f(x)√ 1− x2 dx (47)
Let s.t. x ≥ 1− ⇒ |dντ,ξ −A(1− x)τ | ≤ B(1− x)τ (48)
We observe that for 0 < x < 1− , f(x) = dντ,ξ (1−x)α(1+x)β is bounded. We get from an application of 10, and the observation that P̃α,βt = Ntpα,βt , with Nt = Θ(t−1/2):∫ 1
0 (1− x)τ P̃α,βt (x)2 dx︸ ︷︷ ︸ Θ(hατ ) =
∫ 1− 0
(1− x)τ P̃α,βt (x)2 dx︸ ︷︷ ︸ Θ(t−1) +
∫ 1 1− (1− x)τ P̃α,βt (x)2 dx⇒ (49)
∫ 1 1−
(1− x)τ P̃α,βt (x)2 dx = Θ(hατ ) (50)∫ 1 0 P̃α,βt (x) 2 dντ,ξ(x) = ∫ 1− 0 P̃α,βt (x) 2f(x)(1− x)α(1 + x)β dx)︸ ︷︷ ︸
Θ(t−1)
+ Θ  ∫ 1
1− (1− x)τ P̃α,βt (x)2 dx︸ ︷︷ ︸
Θ(hατ )
 (51)
Theorem 3 (Optimal Rates). Let ν follow Assumption 1. The optimal asymptotic averagecase rates for E[f(xt)− f(x?)] and E[||∇f(xt)||22] are attained by the GCM with parameters (τ, ξ + 2) and (τ, ξ + 3), respectively, and read
E[f(xt)− f(x?)] = Θ(t−2(ξ+2)), E[||∇f(xt)||22] = Θ(t−2(ξ+3)).
Proof. We will prove that for τ, ξ > −1 If α = τ and β = ξ + l + 1 (i.e., are optimal), the rate of convergence reads
min Pt(0)=1
∫ P 2t (λ)λ ldν(λ) = Θ (∫ l 0 P̃α,βt (λ) 2(L− λ)τλξ+l dλ ) = Θ(t−2(ξ+l+1)) (52)
Showing the second equality is easy by considering theorem 2, and that is further the minimum asymptotic rate for the Beta distribution µτ,ξ. By setting pνt and P νt = pνt pνt (0)
the optimal orthonormal and residual and polynomials w.r.t. ν we show that P νt must have the same rate on µτ,ξ as it does on ν, thus the optimal rate of ν cannot be lower than the optimal rate of µτ,ξ. Indeed, setting 1, 2 as in eq. 48:∫ 1
1− 2 P νt (x) 2dν(x) =
( Θ ∫ 1 1− 2 P νt (x) 2dµτ,ξ(x) ) (53)∫ −1+ 1
−1 P νt (x) 2dν(x) = Θ (∫ −1+ 1 −1 P νt (x) 2dµτ,ξ(x) ) (54)∫ 1− 2
−1+ 1 P νt (x) 2dν(x) = Θ (∫ 1− 2 −1+ 1 P νt (x) 2dµτ,ξ(x) ) = Θ ( 1 pνt (−1)2 )
(55)
(56) Where the first two equations come from the fact that ν = Θ(µτ,ξ) near −1 and 1 and the third from lemma 10. This effectively upper bounds the rates on ν because the rates of P νt on µτ,ξ can’t be lower than −2(ξ + l + 1). Proposition 3 (GCM worst-case rates). Let f be a convex, L-smooth quadratic function. Then, for the Generalized Chebyshev Method with parameters (α, β), we have worst-case rates
f(xt)− f(x?) ≤ C1L  t 2(α−β) if α > β − 1 t−1−2β , if α ≤ β − 1 β ≤ 12 t−2, if α ≤ β − 1 β ≥ 12 , (15)
||∇f(xt)− f(x?)|| ≤ C2L2  t 2(α−β) if α > β − 2 t−1−2β , if α ≤ β − 2 β ≤ 3/2 t−4, if α ≤ β − 2 β ≥ 3/2 . (16)
Proof. rates] We will prove that: supx∈[0,L] xlP α,β t (x) 2 = O(Lltv(α,β,l)). Where:
v(α, β, l) =  2(α− β) if α > β − l −1− 2β, if α ≤ β − l β ≤ l − 12 −2l, if α ≤ β − l β ≥ l − 12
(57)
From Szegö (1975), Theorem 7.32.2, if θ < π2 :
P̃α,βt (cos θ) =  O(t −1/2) if α < − 12
O(tα) if α ≥ − 12 , 0 ≤ θ ≤ ct−1 θ−α−1/2O(t−1/2) if α ≥ − 12 , θ > ct−1
(58)
We observe that, from the symmetry of the Jacobi polynomials:
sup x∈[0,L]
xlPα,βt (x) 2 = Θ
( max { sup x∈[0,1] Pα,βt (x) 2, sup x∈[0,1] (1− x)lP β,αt (x)2 })
(59)
The (1 − x)l term, corresponds to (2 sin( θ2 ))l in the variable θ, which is O(θ2l). The rest follows from carefully considering the expressions given by eq. 58.
Theorem 4 (Nesterov average-case rates). Let ν as in Assumption 1. Then for the Nesterov method, we have average-case rates
E[f(xt)− f(x?)] ∼ C ′1,ν { t−2(ξ+2) if ξ < −1/2 t−3 log t if ξ = −1/2 t−(ξ+7/2) if ξ > −1/2 , E[||∇f(xt)||22] ∼ C ′2,νt−(ξ+9/2). (17)
Proof. We will prove:∫ 1 0 PNest (λ) 2λl dντ,ξ−l ∼ C ′ν { t−2(ξ+1) if 0 < ξ < 1/2 t−3 log t if ξ = 1/2 t−(ξ+5/2) if ξ > 1/2
(60)
Paquette et al. (2020) has shown that the nesterov polynomials Pt are asymptotically, in t:
Pt(λ) ∼ 2J1(t
√ αλ)
t √ αλ e−αλt/2 (61)
In the sense that:∫ 1 0 ul [ ˜P 2t (u)− 4J21 (t √ u) t2u e−ut ] 4 dµMP (u) = O(t −(l+25/12)) (62)
The arguments can be easily used to show that such an integral is O(t−(α+l+31/12)) when evaluated wrt a general dµ s.t µ′ = Θ(λα) near 0. We can thus consider our integral of interest substituting PNest by it’s Bessel asymptotic and dividing it into three regions, i.e. [0, 1] = [0, t ] ∪ [ t , √t ] ∪ [ √ t , 1] corresponding to two different regimes for the Bessel function. The first region will give us the asymptotic and the others we will bound. We consider first, for some > 0:∫ √
t
t
uξ 4J21 (t
√ u)
t2u e−ut du (63)
We note the asymptotic for J21 :
J21 ( √ tv) ∼ 1
π √ tv
(1 + cos(2 √ tv + 2γ)) (64)
Doing the change of variable v = tu, and identifying the upper limit of the interval, which is t1/2 to ∞: ∫ √
t
t
uξ 4J21 (t
√ u)
t2u e−ut du = Θ
( t−2−ξ ∫ ∞ vξ−1J21 ( √ tv)e−v dv ) (65)
= Θ ( t−2−ξ ∫ ∞ vξ−1 1
π √ tv e−v dv
) (66)
= Θ t− 52−ξ ∫ ∞ vξ− 3 2 1
π √ tv e−v dv︸ ︷︷ ︸
Γ(ξ− 12 , )  (67) Where the cosinus term goes to 0 from the Riemann-Lebesgue lemma and Γ is the incomplete Gamma function. The term corresponding to the interval [ t−1/2, 1] is exponentially small. Indeed, because of the exponential e−ut it is O(e− √ t). This shows that the integral concentrates in a region that is closer and closer to 0 and that only the behaviour of the distribution near 0 matters. We have for the [0, t ] region, doing the change of variables v = t
2u:∫ t
0
uξ 4J21 (t
√ u)
t2u e−ut du = Θ
( t−2(ξ+1) ∫ t 0 vξ J21 ( √ v) v e− v t dv ) (68)
And the e −v t is Θ(1). We have the following Bessel asymptotics:
J21 ( √ v)
v ∼ 1 4 , v → 0 (69)
J21 ( √ v)
v = O(v−3/2), v →∞ (70)
So we divide this integral aswell:
t−2(ξ+1) ∫ t
1
vξ J21 ( √ v)
v e− v t dv = Θ
( t−2(ξ+1) ∫ t vξ− 3 2 dv ) = Θ ( Iξ(t)t −ξ− 52 )
(71)
t−2(ξ+1) ∫ 1
0
vξ J21 ( √ v)
v e− v t dv = Θ
( t−2(ξ+1) ∫ 0 1vξ dv ) = Θ ( t−2(ξ+1) ) (72)
Where Iξ(t) = log t if ξ = 12 and 1 otherwise. The nesterov rate is then Iξ(t)t−ξ− 5 2 if ξ ≥ 12 and t−2(ξ+1) if 0 < ξ < 12
Theorem 5 (Gradient Descent average-case rates). Let ν as in Assumption 1. Then for gradient descent
E[f(xt)− f(x?)] = Θ(t−(ξ+2)), E[||∇f(xt)||22] = Θ(t−(ξ+3)). (18)
Proof. Considering that PGDt (λ) = (1− λL )t we will prove :∫ 1 0 (1− λ)2tλl dντ,ξ−l = Θ(t−(ξ+l+1) (73)
We know, for the Beta weights, that:∫ 1 0 (1− λ)2t+τλξ+l dλ = Γ(l + ξ + 1)Γ(2t+ τ + 1) Γ(2t+ l + ξ + τ + 2) = Θ(t−(ξ+l+1)) (74)
We can indentify this asymptotic to the interval ∫
0 for any because:∫ 1
(1− λ)2t+τλξ+l dλ = O((1− )2t) (75)
Then: ∫ 1
(1− λ)2tλl dντ,ξ−l = O((1− )2t) (76)∫ 0 (1− λ)2tλl dντ,ξ−l = Θ (∫ 0 (1− λ)2t+τλξ+l dλ ) = Θ(t−(ξ+l+1)) (77)
Theorem 6 (Laguerre method rates). Let α > −1 and µα, α > −1 be a Gamma distribution, i.e. dµα(x) = x αe−x
Γ(α+1) . The optimal rates are given by the Laguerre method of appropriate tuning and
E[f(xt)− f(x?)] = Θ(t−(α+2)). (19)
Proof. Let Lαt be the Laguerre polynomials with the usual normalization Szegö (1975):∫ Lαt (x) 2dµα(x) = L α t (0) = ( n+ α
n
) (78)
We further now [Szegö (1975) (5.1.13)]]: t∑
k=0
Lαt (x) = L α+1 t (x) (79)
Thus, letting Pαt be the residual laguerre polynomial, we consider: E[f(xt)− f(x?)] = ∫ Pα+2t (λ) 2dµα+1(λ) = ( t+ α+ 2
t
)−2 ∫ Lα+2t dµα+1(λ)
=
( t+ α+ 2
t
)−2 t∑ k=0 [∫ Lα+1k (λ)dµα+1(λ) ]
=
( t+ α+ 2
t
)−2 t∑ k=0 ( k + α+ 1 k ) = ( t+ α+ 2 t )−2( t+ α+ 2 t )
=
( t+ α+ 2
t
)−1 = Θ(t−(α+2))
(80)

heading: D Additional Experiments
text: We note that in the regimes where the empirical average performance doesn’t match the theoretical one, we can still find samples of problems who do match. This and the much larger variance on the function-value, this discrepancy is due to numerical unstability in these regimes.

Review summaries:
-- EDIT: I have updated my scores in response to clarifications --  The problem of optimizing a convex quadratic function via first order methods is considered. This is a well-understood problem from the worst case point of view, and its complexity will depend on the largest and smallest eigenvalues of the associated Hessian matrix. However, a nice recent average-case analysis by Pedrogosa and Scieur gives the following result: if the spectral density of the Hessian converges to a "nice" probability measure (such as the Manchenko-Pastur law), then first-order methods that are tailored to this density may converge faster. In fact, such methods can be obtained from the three-term recurrence for the orthogonal polynomials of the measure.  The present manuscript follows up on the Pedregosa/Scieur, but departs from in two specific ways. Firstly, it considers the case where the limiting spectral law is a beta distribution, and derives the corresponding iteration. Secondly, the authors analyze the performance of this iteration under "mispecification": that is, they allow for more general spectral measures than the beta, and only specify how these measures behave near the extremes. In spite of this, they are able to derive bounds of the same order as in the beta case.   A few other theoretical results are presented: worst-case rates for the method derived from the beta distribution; an analysis of Nesterov's method under conditions on the tail of the spectral measure, and a result on Laguerre spectral measures. These results generalize other theorems obtained in the Pedregosa/Scieur paper. The theory is complemented by a small numerical study.  

The paper considers the problem of average convergence rate of first order methods on a given ensemble of quadratic problems. The authors propose the Generalized Chebyshev Method (GCM) and show that it is optimal when the e.s.d. is beta distribution. They also show that so long as we know the behavior of e.s.d. near the edges of its support, GCM still achieves the optimal rate. The authors finally consider the Nestrov method and derive its asymptotic average-case convergence rate.