# Only Tails Matter: Average-Case Universality and
Robustness in the Convex Regime

Leonardo Cunha  Gauthier Gidel  Fabian Pedregosa  Damien Scieur  Courtney
Paquette

###### Abstract

The recently developed average-case analysis of optimization methods allows a
more fine-grained and representative convergence analysis than usual worst-
case results. In exchange, this analysis requires a more precise hypothesis
over the data generating process, namely assuming knowledge of the expected
spectral distribution (ESD) of the random matrix associated with the problem.
This work shows that the concentration of eigenvalues near the edges of the
ESD determines a problem’s asymptotic average complexity. This a priori
information on this concentration is a more grounded assumption than complete
knowledge of the ESD. This approximate concentration is effectively a middle
ground between the coarseness of the worst-case scenario convergence and the
restrictive previous average-case analysis. We also introduce the Generalized
Chebyshev method, asymptotically optimal under a hypothesis on this
concentration and globally optimal when the ESD follows a Beta distribution.
We compare its performance to classical optimization algorithms, such as
gradient descent or Nesterov’s scheme, and we show that, in the average-case
context, Nesterov’s method is universally nearly optimal asymptotically.

Machine Learning, ICML

##  1 Introduction

The analysis of the average complexity of algorithms has a long tradition 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 to distributions of quadratic minimization
problems. This derivation is done by relating the average-case convergence
rate to the expected spectral distribution (ESD) of the objective function’s
Hessian, which is a well-studied object in random matrix theory. In practice,
however, the knowledge of the ESD is a much stronger requirement than the
worst-case analysis, which relies only on knowledge of the distribution’s
support.

Paquette et al. (2022) extended the average-case framework by introducing a
noisy generative model for the problems. Among other results, they derived the
average complexity of the Nesterov Accelerated Method (Nesterov, 2003) on a
particular distribution. Moreover, they showed the concentration of
convergence metrics in the infinite sample and dimensional limit.

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.

Figure 1:  Representation of different spectra with different concentrations
of eigenvalues around the support edges. These concentrations determine the
average-case rates for non-strongly problems.

###  1.1 Current limitations of the average-case analysis

This paper addresses some of the limitations in previous works on average-case
analysis. First, little is known about the convergence rate on convex but not
strongly-convex problems. Also, optimal average-case algorithms require an
exact estimation of the ESD to guarantee an optimal convergence rate;
otherwise, their convergence rate under inexact ESD is unknown. Finally, the
non-smooth has not yet been discussed in its full generality.

Convex problems. In the smooth and strongly convex case, the optimal worst-
case and average-case convergence rates are asymptotically equal (Scieur &
Pedregosa, 2020). However, little is known about optimal average-case rates
for non-strongly convex problems, as well as the average-case complexity of
classical methods such as gradient descent or Nesterov’s method, see (Paquette
et al., 2022).

Exact estimation of the ESD. Pedregosa & Scieur (2020) developed average-case
optimal algorithms that required an exact estimation of the ESD of the problem
class. Such an estimation may be complex, if not impossible, to obtain in
practice. Despite showing good performance when the ESD is estimated with
empirical quantities, there are no theoretical guarantees of the method’s
performance when the ESD is poorly estimated. Therefore, there is a need to
analyze the algorithm’s performance under different notions of uncertainty on
the spectrum, allowing a practitioner to choose the best algorithm for a
practical problem even with imperfect a priori information.

Non-smooth objectives. In this paper we provide an average-case analysis for
the generalized Laguerre distribution
$\lambda^{\alpha}e^{-\lambda},\,\alpha>-1$. This is a distribution with
unbounded support, i.e., when the largest eigenvalue (also known as smoothness
constant) is bounded. This extends the results of Pedregosa & Scieur (2020)
for the Laguerre distribution $e^{-\lambda}$.

###  1.2 Contributions

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 ESD around the edges of
their support. Using this analysis, we also derive a family of average-case
optimal algorithms, analyze their robustness, and finally exhibit a
universality result for Nesterov’s method.

  * •

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 ESD around the edges of their support. With the proper
set of parameters, the GCM method converges at an optimal average-case rate
(Theorem 4.2 for smooth problems, Theorem 4.6 for non-smooth problems). These
rates are faster than optimal worst-case methods like Nesterov acceleration,
and we recover the classical worst-case rates as limits of the average-case
(see Table 2). Fig. 7 shows our theoretical analysis to match the practical
performance of the algorithms.

  * •

Robustness. Developing an optimal algorithm requires knowledge of the exact
ESD. However, in practical scenarios, we only have access to an approximation
of the ESD. In Theorem 4.1 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 2).

  * •

Universality. Finally, in Theorem 4.4, we analyze the asymptotic average-case
convergence rate of Nesterov’s method. We show that its convergence rate is
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 distribution. This observation contributes to the theoretical
understanding of the empirical effectiveness of Nesterov’s acceleration.

Regime | Worst-case | Average-Case
---|---|---
Strongly conv. | $\left(1-\Theta(\nicefrac{{1}}{{\sqrt{\kappa}}})\right)^{t}$ | $\left(1-\Theta(\nicefrac{{1}}{{\sqrt{\kappa}}})\right)^{t}$
Smooth conv. | $\nicefrac{{1}}{{t^{2}}}$ | $\nicefrac{{1}}{{t^{2\xi+4}}}$
Convex | $\nicefrac{{1}}{{\sqrt{t}}}$ | $\nicefrac{{1}}{{t^{\alpha+2}}}$
Figure 2: Comparison between function value worst-case and average-case
convergence. $\kappa$ is the condition number in the smooth strongly convex
case. In the smooth convex case $\xi>-1$ is the concentration of eigenvalues
around $0$ and in the non-smooth case we consider
$\mathop{}\\!\mathrm{d}\mu\propto\lambda^{\alpha}e^{-\lambda}d\lambda$.

##  2 Average-Case Analysis

This section recalls the average-case analysis framework for random quadratic
problems. The main result is Theorem 2.5, which relates the expected error to
the expected spectral 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 as follows:

######  Problem 2.1.

Let ${\boldsymbol{H}}\in{\mathbb{R}}^{d\times d}$ be a random symmetric
positive-definite matrix independent of
${\boldsymbol{x}}^{\star}\in{\mathbb{R}}^{d}$, which is a random vector and a
solution to the problem. We define the random quadratic minimization problem
as

$\vphantom{\sum_{0}^{i}}\min_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}\Big{\\{}f({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\\!\mfrac{1}{2}({\boldsymbol{x}}\\!-\\!{\boldsymbol{x}}^{\star})^{\top}\\!{\boldsymbol{H}}({\boldsymbol{x}}\\!-\\!{\boldsymbol{x}}^{\star})\Big{\\}}\,.$
(OPT)

We are interested on minimizing the expected errors
${\mathbb{E}}\|f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})\|$, the
expected function-value gap, and ${\mathbb{E}}\|\nabla
f({\boldsymbol{x}}_{t})\|^{2}$, the expected gradient norm, where
${\boldsymbol{x}}_{t}$ is the $t$-th update of a first-order method starting
from ${\boldsymbol{x}}_{0}$ and ${\mathbb{E}}$ is the expectation over the
random variables ${\boldsymbol{H}},{\boldsymbol{x}}_{0}$ and
${\boldsymbol{x}}^{\star}$.

Note that the expectations we consider throughout the paper are over problem
instances, 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
${\boldsymbol{x}}_{t}$ as

${\boldsymbol{x}}_{t}\in{\boldsymbol{x}}_{0}+\operatorname*{\mathbf{span}}\\{\nabla
f({\boldsymbol{x}}_{0}),\ldots,\nabla f({\boldsymbol{x}}_{t-1})\\}\,,$ (1)

that is, ${\boldsymbol{x}}_{t}$ 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 Methods to Polynomials.

There is an intimate link between first-order methods and polynomials that
simplifies the analysis of quadratic objectives. The following proposition
shows that, with this link, we can assign to each optimization method a
polynomial that determines its convergence. We will denote $P_{t}(\lambda)$
and $P_{t}({\boldsymbol{H}})$ for the polynomials taking arguments in
${\mathbb{R}}$ and ${\mathbb{R}}^{n\times n}$ respectively with the same
coefficients. Following (Fischer, 1996), we will say a polynomial $P_{t}$ is
residual if $P_{t}(0)=1$.

######  Proposition 2.2.

(Hestenes et al., 1952) Let ${\boldsymbol{x}}_{t}$ be generated by a first-
order method. Then there exists a residual polynomial $P_{t}$ of degree $t$,
that verifies

$\vphantom{\sum^{n}}{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}=P_{t}({\boldsymbol{H}})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})\leavevmode\nobreak\
.$ (2)

######  Remark 2.3.

If the first-order method is further a momentum method, i.e.

${\boldsymbol{x}}_{t+1}={\boldsymbol{x}}_{t}+h_{t}\nabla
f({\boldsymbol{x}}_{t})+m_{t}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}_{t-1}),$

we can determine the polynomials by the recurrence $P_{0}=1$ and

$P_{t+1}(\lambda)=P_{t}(\lambda)+h_{t}\lambda
P_{t}(\lambda)+m_{t}(P_{t}(\lambda)-P_{t-1}(\lambda)).$

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 2.4 (Expected spectral distribution (ESD)).

. Let ${\boldsymbol{H}}$ be a random matrix with eigenvalues
$\\{\lambda_{1},\ldots,\lambda_{d}\\}$. The empirical spectral distribution of
${\boldsymbol{H}}$, called ${\mu}_{{\boldsymbol{H}}}$, is the probability
measure

$\mu_{{\boldsymbol{H}}}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{d}{\textstyle{\sum_{i=1}^{d}}}\delta_{\lambda_{i}},$
(3)

where $\delta_{\lambda_{i}}$ is the Dirac delta, a distribution equal to zero
everywhere except at $\lambda_{i}$ and whose integral over the entire real
line is equal to one.

Since ${\boldsymbol{H}}$ is random, the empirical spectral distribution
$\mu_{\boldsymbol{H}}$ is a random variable in the space of measures. Its
expectation over ${\boldsymbol{H}}$ is called the expected spectral
distribution (ESD) and we denote it

$\mu\stackrel{{\scriptstyle\text{def}}}{{=}}{\mathbb{E}}_{{\boldsymbol{H}}}[\mu_{{\boldsymbol{H}}}]\,.$
(4)

The following theorem links the ESD to the average-case convergence of a
first-order method when ${\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}$ and
${\boldsymbol{H}}$ are independent.

######  Theorem 2.5.

Let ${\boldsymbol{x}}_{t}$ be generated by a first-order method with
associated residual polynomial $P_{t}$, $\mu$ be the ESD, and
$\mathbb{E}[({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}]=R^{2}\textbf{I}$
for some constant $R$. Then we have the following identities for different
convergence metrics:

$\displaystyle\mathbb{E}[\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}]={R^{2}}\int{P_{t}^{2}(\lambda)\mathop{}\\!\mathrm{d}\mu(\lambda)}\,,$
(5)
$\displaystyle\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]=\frac{R^{2}}{2}\int
P_{t}^{2}(\lambda)\lambda\mathop{}\\!\mathrm{d}\mu(\lambda)\,,$ (6)
$\displaystyle\mathbb{E}[\|\nabla f({\boldsymbol{x}}_{t})\|^{2}_{2}]=R^{2}\int
P_{t}^{2}(\lambda)\lambda^{2}\mathop{}\\!\mathrm{d}\mu(\lambda)\,.$ (7)

This theorem states that polynomials are a powerful abstraction, allowing us
to write all our convergence metrics within the same framework. Therefore, in
the rest of the paper, we will refer directly to the polynomials associated
with a given method. For simplicity, we set $R^{2}=1$.

This framework is linked to the field of orthogonal polynomials by the
following proposition. We construct an optimal method w.r.t. a given
distribution through a family of orthogonal polynomials.

######  Proposition 2.6.

Let $\nu$ be a distribution with continuous support and $P_{t}^{l}$ be defined
as

$P_{t}^{l}\stackrel{{\scriptstyle\text{def}}}{{=}}{\arg\min}_{P_{t}(0)=1}\int
P_{t}^{2}(\lambda)\lambda^{l}d\nu(\lambda)\,.$ (8)

Then $(P_{t}^{l})_{t=0,1,\ldots}$ is the family of residual orthogonal
polynomials w.r.t. to $\lambda^{l+1}d\nu$.

We refer to objective $l$ as the one associated to the added $\lambda^{l}$
term, i.e. the function-value is objective $l=1$. This proposition further
implies that the optimal first-order method is a momentum method because
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,

$P_{t+1}(\lambda)=(a_{t}+b_{t}\lambda)P_{t}(\lambda)+(1-a_{t})P_{t-1}(\lambda)\,.$
(9)

Following Remark 2.3, the optimal method is derived from this recurrence as

${\boldsymbol{x}}_{t+1}={\boldsymbol{x}}_{t}+(a_{t}-1)({\boldsymbol{x}}_{t}-{\boldsymbol{x}}_{t-1})+b_{t}\nabla
f({\boldsymbol{x}}_{t})\,.$ (10)

##  3 Methods

Writing the rates in terms of the expected spectral distribution ties the
average-case framework to the field of random matrix theory. A classical
result of this field is that the same spectral distribution arises from
different input distributions in the data. For example, the Marchenko-Pastur
distribution arises in the large sample and dimension limit when we take the
Gram matrix of a matrix where the entries are generated i.i.d. from any
distribution with mean zero, variance $\sigma^{2}$ and with a bounded fourth
moment.

######  Definition 3.1.

The Marchenko-Pastur distribution associated with the parameter $r$ and with
scale $\sigma^{2}$ is given by

$\mathop{}\\!\mathrm{d}\mu_{MP}(\lambda)=\frac{1}{2\pi\sigma^{2}}\frac{\sqrt{(\lambda^{+}-\lambda)(\lambda-\lambda^{-})}}{r\lambda}\mathop{}\\!\mathrm{d}\lambda\,,$
(11)

with $\lambda^{+}=\sigma^{2}(1+\sqrt{r})^{2}$,
$\lambda^{-}=\sigma^{2}\max(0,(1-\sqrt{r})^{2})$.

$r$ in the Gram matrix case, is determined as $n/d$. The Marchenko-Pastur
distribution $\mu_{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
$\mathop{}\\!\mathrm{d}\mu_{MP}\propto\lambda^{-1/2}\sqrt{\lambda^{+}-\lambda}\mathop{}\\!\mathrm{d}\lambda$.

Pedregosa & Scieur (2020) first derived the optimal method w.r.t. $\mu_{MP}$,
and Paquette et al. (2022) derived Nesterov’s rates under the distribution. In
this paper we take a more general view and consider Beta distribution.

######  Definition 3.2.

The (generalized) Beta distribution with parameters $\tau,\xi$ and scale $L$
are given by the (non-normalized) probability density function

$\mathop{}\\!\mathrm{d}\mu_{\tau,\xi}(\lambda)\propto\lambda^{\xi}(L-\lambda)^{\tau}\mathop{}\\!\mathrm{d}\lambda\,,$
(12)

where $\mathop{}\\!\mathrm{d}\lambda$ denotes the standard Lebesgues measure
on $\mathbf{R}$. This family of distributions generalizes the Marchenko-Pastur
distribution, and both have similar concentrations near $0$ when
$\xi\approx-1/2$.

The optimal method w.r.t. $\mu_{\tau,\xi}$ and objective $l$ is associated to
a shifted Jacobi polynomial $\tilde{P}_{t}^{\alpha,\beta}$ with
$\beta=\xi+l+1,\alpha=\tau$. This is a direct implication of the definition of
the Jacobi polynomials as the orthogonal family w.r.t. the weights
$(1-\lambda)^{\alpha}(1+\lambda)^{\beta}$. When $\alpha=\beta=-1/2$, we
retrieve the Chebyshev Method (Flanders & Shortley, 1950). It is also related
to the “$\nu$-method” of Brakhage (1987), who considers the conjugate gradient
algorithm applied to a matrix with a Gamma spectral distribution
($\beta=-1/2$).

We name the following method the Generalized Chebyshev Method (GCM).

Algorithm 1 Generalized Chebyshev Method GCM($\alpha,\beta$)

Inputs: Initial vector ${\boldsymbol{x}}_{0}$, distribution parameters
$\alpha,\beta$, $L=$largest eigenvalue of ${\boldsymbol{H}}$ .
${\boldsymbol{x}}_{-1}\leftarrow{\boldsymbol{x}}_{0},\delta_{0}\leftarrow 0$

for $t=1,\ldots,T$ do

$a_{t}\leftarrow-\frac{2\left(\beta^{2}+\alpha\beta+(2t+1)(\alpha+\beta)+2t^{2}+2\right)(2t+\alpha+\beta+1)}{2(t+1)(t+\alpha+\beta+1)(2t+\alpha+\beta)}$

$b_{t}\leftarrow\frac{(2t+\alpha+\beta+1)(2t+\alpha+\beta+2)}{L(t+1)(t+\alpha+\beta+1)}$

$\gamma_{t}\leftarrow-\frac{(t+\alpha)(t+\beta)(2t+\alpha+\beta+2)}{(t+1)(t+\alpha+\beta+1)(2t+\alpha+\beta)}$

$\delta_{t}\leftarrow\frac{1}{a_{t}+\gamma_{t}\delta_{t-1}}$

${\boldsymbol{x}}_{t}\leftarrow{\boldsymbol{x}}_{t-1}+(\delta_{t}a_{t}-1)({\boldsymbol{x}}_{t-1}-{\boldsymbol{x}}_{t-2})+\delta_{t}b_{t}\nabla
f({\boldsymbol{x}}_{t-1})$

end for

We consider the Nesterov’s method used in (Paquette et al., 2022), which is
defined by the iterations:

$\displaystyle{\boldsymbol{x}}_{t+1}$
$\displaystyle={\boldsymbol{y}}_{t}-\frac{1}{L}\nabla
f({\boldsymbol{y}}_{t}),$ (13) $\displaystyle{\boldsymbol{y}}_{t+1}$
$\displaystyle={\boldsymbol{x}}_{t+1}+\frac{t}{t+3}({\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}_{t})\,.$
(14)

We also consider the Laguerre method (Pedregosa & Scieur, 2020), which is
optimal w.r.t.
$\mathop{}\\!\mathrm{d}\mu(\lambda)=\frac{\lambda^{\alpha}e^{-\lambda}}{\Gamma(\alpha+1)}\mathop{}\\!\mathrm{d}\lambda$,
taking $\alpha$ as a parameter. This method is proposed to optimize non-smooth
functions.

Both these methods are generalizations of the ones in (Pedregosa & Scieur,
2020). We show that Algorithm 1 corresponded to polynomials
$\tilde{P}_{t}^{\alpha,\beta}$ and derive the Laguerre method in Appendix B.

######  Remark 3.3.

The Generalized Chebyshev takes the largest eigenvalue $L$ as a parameter, but
the rates we will show are robust to an overestimation of $L$.

##  4 Robust Average-Case Rates

Throughout the rest of the paper, we make the following assumption about the
spectral distributions:

######  Assumption.

Let $\nu_{\tau,\xi}$ be a distribution supported in $(0,L]$ s.t.
$\nu_{\tau,\xi}^{\prime}(x)>0$ for $x\in[0,L]$,
$d\nu_{\tau,\xi}=\Theta(\lambda^{\xi})$ near $0$ and
$d\nu_{\tau,\xi}=\Theta((L-\lambda)^{\tau})$ near $L$.

We characterize our distributions of interest only in terms of the behavior at
the edges. This assumption is sufficient to determine the asymptotic
convergence of algorithms. This (mild) assumption excludes only distributions
that decay exponentially near their edges, in which case they effectively
behave as strongly convex problems. Moreover, this assumption allows
considering a broader class of problems than previous work that assumed
complete knowledge of the spectrum of the Hessian.

The $\xi$ parameter measures how close we are to the worst-case scenario as it
approaches $-1$. Samples in finite dimension of distributions with high values
of $\xi$ will work as strongly convex functions in practice, as samples of low
eigenvalues are rarer and rarer with increasing $\xi$.

We show that the coefficients $\xi,\,\tau$ determine 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 for all
distributions with the same concentrations.

######  Theorem 4.1 (GCM average-case rates).

The Generalized Chebyshev Method with parameters $(\alpha,\beta)$ applied to a
problem with expected spectral distribution $\nu_{\tau,\xi}$ has average-case
rates ${\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]}\sim
L\cdot C^{\alpha,\beta}_{1,\nu}t^{e_{1}}$ and $\mathbb{E}[\|\nabla
f({\boldsymbol{x}}_{t})\|^{2}_{2}]\sim L^{2}\cdot
C^{\alpha,\beta}_{2,\nu}t^{e_{2}}$, with the exponents

$\displaystyle e_{1}=$
$\displaystyle\left\\{\begin{array}[]{ll}-1-2\beta\hskip
8.5359pt\mbox{if}\hskip
8.5359pt\alpha<\tau+\tfrac{1}{2}\;\;\land\;\;\beta<\xi+\tfrac{3}{2},&\\\
-2(\xi+2)\log t\hskip 8.5359pt\mbox{if}\hskip
8.5359pt\alpha=\tau+\tfrac{1}{2}\;\;\land\;\;\beta=\xi+\tfrac{3}{2},&\\\
2(\max\\{\alpha-\beta-\tau,-\xi-1\\}-1)\hskip
8.5359pt\mbox{otherwise},&\end{array}\right.$ $\displaystyle e_{2}=$
$\displaystyle\left\\{\begin{array}[]{ll}-1-2\beta\hskip
8.5359pt\mbox{if}\hskip
8.5359pt\alpha<\tau+\tfrac{1}{2}\;\;\land\;\;\beta<\xi+\tfrac{5}{2},&\\\
-2(\xi+3)\log t\hskip 8.5359pt\mbox{if}\hskip
8.5359pt\alpha=\tau+\tfrac{1}{2}\;\;\land\;\;\beta=\xi+\tfrac{5}{2},&\\\
2(\max\\{\alpha-\beta-\tau,-\xi-2\\}-1)\hskip
8.5359pt\mbox{otherwise}\,,&\end{array}\right.$

where $\land$ denotes the logical “ _and_ ” operator and
$C^{\alpha,\beta}_{.,\nu}$ are distribution dependent constants.

###### Proof Sketch.

We first compute the convergence rates assuming that the ESD is a generalized
Beta distribution with parameters $\tau,\xi$. For this we use Theorem 2.5 and
Lemma C.1. For a general $\nu_{\tau,\xi}$, with the aid of Lemma 64, we then
show that the mass away from the edges is negligible, i.e.,

$\int_{\epsilon}^{L-\epsilon}P_{t}^{\alpha,\beta}(\lambda)^{2}\lambda^{l}d\nu_{\tau,\xi}\,,$
(15)

is $O(t^{-1-2\beta})=O(t^{e_{.}}),\forall\epsilon>0$. As the mass near the
edges must be similar to the one for the Beta weights, the result follows. ∎

Theorem 4.1, illustrated by Fig. 2, shows that overestimating $\beta$ and
underestimating $\alpha$ will still leave us with the optimal asymptotic
rates. So a good rule of thumb for calibrating the algorithm is to use high
$\beta$ and low $\alpha$.

Theorem 4.2 shows that a proper choice of $\alpha,\beta$ can indeed make the
Jacobi polynomial asymptotically optimal w.r.t. to any $\nu_{\tau,\xi}$.

######  Theorem 4.2 (Optimal Rates).

Consider the distribution $\nu_{\tau,\xi}$. The optimal asymptotic average-
case rates for
$\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]$ and
$\mathbb{E}[\|\nabla f({\boldsymbol{x}}_{t})\|^{2}_{2}]$ are attained by the
GCM with parameters $(\tau,\xi+2)$ and $(\tau,\xi+3)$, respectively, and read

$\displaystyle\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]=\Theta(t^{-2(\xi+2)})\,,$
(16) $\displaystyle\mathbb{E}[\|\nabla
f({\boldsymbol{x}}_{t})\|^{2}_{2}]=\Theta(t^{-2(\xi+3)})\,.$ (17)

###### Proof Sketch.

When $\nu_{\tau,\xi}$ is a Beta distribution, the result follows from Theorem
4.1.

For a general $\nu_{\tau,\xi}$, we consider the optimal method for this
expected spectral distribution. We show that its performance on the
corresponding Beta weight $\mu_{\tau,\xi}$ is similar to that of
$\nu_{\tau,\xi}$, as the integrals concentrate on the edges. It follows that
the optimal performance for $\nu_{\tau,\xi}$ is asymptotically the same as for
$\mu_{\tau,\xi}$. ∎

Figure 3:  The figure illustrates the robustness of the Generalized Chebyshev Method with parameters $(\alpha,\beta)$ for a _fixed problem_ corresponding to the Marchenko-Pastur distribution ${(\tau=\tfrac{1}{2},\xi=-\tfrac{1}{2})}$. The color represents the exponent $a$ of the average-case rate $O(t^{a})$ of the method for different values of $\alpha$ and $\beta$. The white star represents the optimal tuning and the blue area is the set of parameters for which the method converges. Note that a large region guarantees the same optimal asymptotic rate. Method |  Parameters $(\tau,\xi)$
---|---
$(\frac{1}{2},\frac{1}{2})$ | $(\frac{1}{2},-\frac{1}{2})$
GCM$(\alpha=\frac{1}{2},\beta=\frac{5}{2})$ | $t^{-5}$ | $t^{-3}$
GCM$(\alpha=\frac{1}{2},\beta=\frac{3}{2})$ | $t^{-4}$ | $t^{-3}$
Nesterov | $t^{-4}$ | $t^{-3}\log t$
Gradient descent | $t^{\frac{-5}{2}}$ | $t^{\frac{-3}{2}}$
Table 1:  The table compares the asymptotic average-case rates for the
function-value for different methods and pairs $(\tau,\xi)$.

For the function value ($l=1$), we find rates that approach $t^{-2}$ as
$\xi\rightarrow-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 $\beta$ up to $1/2$ below the
optimal choice and infinitely above.

For completeness, we derive worst-case rates for the GCM.

######  Proposition 4.3 (GCM worst-case rates).

Let $f$ be a convex, L-smooth quadratic function. Then, for the Generalized
Chebyshev Method with parameters $(\alpha,\beta)$, we have worst-case rates
$f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})\leq C_{1}Lt^{e_{3}}$ and
$\|\nabla f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})\|\leq
C_{2}L^{2}t^{e_{4}}$, with the exponents

$\displaystyle e_{3}=$
$\displaystyle\left\\{\begin{array}[]{ll}2(\alpha-\beta)\hskip
8.5359pt\mbox{if}\hskip 8.5359pt\alpha>\beta-1,&\\\ -1-2\beta\hskip
8.5359pt\mbox{if}\hskip 8.5359pt\alpha\leq\beta-1\land\beta\leq 1/2,&\\\
-2\hskip 8.5359pt\mbox{otherwise},&\end{array}\right.$ (21) $\displaystyle
e_{4}=$ $\displaystyle\left\\{\begin{array}[]{ll}2(\alpha-\beta)\hskip
8.5359pt\mbox{if}\hskip 8.5359pt\alpha>\beta-2,&\\\ -1-2\beta\hskip
8.5359pt\mbox{if}\hskip 8.5359pt\alpha\leq\beta-2\land\beta\leq 3/2,&\\\
-4\hskip 8.5359pt\mbox{otherwise}.&\end{array}\right.$ (25)

The logical and is noted $\land$.

For common choice of $\alpha,\beta$, i.e. $\beta\geq\frac{1}{2}$,
$\alpha\leq\beta-1$, the rate of decay of the function suboptimality 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 extend
the proof of Paquette et al. (2022, Lemma B.2) for the Nesterov method under
the Marchenko-Pastur distribution.

######  Theorem 4.4 (Nesterov average-case rates).

Consider the distribution $\nu_{\tau,\xi}$. Then for the Nesterov method, we
have average-case rates

$\displaystyle\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]\sim
C^{\prime}_{1,\nu}{\Bigg{\\{}}\begin{array}[]{ll}t^{-2(\xi+2)}&\mbox{if
}\xi<-1/2,\\\ t^{-3}\log t&\mbox{if }\xi=-1/2,\\\ t^{-(\xi+7/2)}&\mbox{if
}\xi>-1/2,\end{array}$ (29) $\displaystyle\mathbb{E}[\|\nabla
f({\boldsymbol{x}}_{t})\|^{2}_{2}]\sim C^{\prime}_{2,\nu}t^{-(\xi+9/2)}.$ (30)

The gap between the asymptotic average-case rates of Nesterov and the optimal
ones is of the order of $t^{\xi+l-1/2}$, when $\xi+l>1/2$, $\log t$ when
$\xi+l=1/2$ and $0$ otherwise. This result shows that Nesterov is almost
optimal when the concentrations near $0$ are relatively high, i.e., low $\xi$.

######  Theorem 4.5 (Gradient descent average-case rates).

Consider the distribution $\nu_{\tau,\xi}$. Then for gradient descent

$\displaystyle\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]=\Theta(t^{-(\xi+2)}),$
(31) $\displaystyle\mathbb{E}[\|\nabla
f({\boldsymbol{x}}_{t})\|^{2}_{2}]=\Theta(t^{-(\xi+3)}).$ (32)

###### Proof Sketch.

From the simple form of the polynomial associated with gradient descent, we
have a closed-form (Eq. (90)) expression for the performance of gradient
descent on the Beta weights. For a general $\nu_{\tau,\xi}$, the mass that is
$\epsilon$ away from $0$ is $O((1-\epsilon)^{2t})$; thus, the performance of
gradient descent is asymptotically the same as on the Beta weights. ∎

Observe for the function value that the rate for Nesterov is $t^{-2}$ and the
rate for gradient descent is $t^{-1}$ when $\xi\rightarrow-1$.

We now consider the optimal rates for a Gamma distribution.

######  Theorem 4.6 (Laguerre method rates).

Let $\alpha>-1$ and $\mu_{\alpha}$ be a Gamma distribution, i.e.
$\mathop{}\\!\mathrm{d}\mu_{\alpha}(\lambda)={\lambda^{\alpha}e^{-\lambda}}/{\Gamma(\alpha+1)}\mathop{}\\!\mathrm{d}\lambda$.
The optimal rates are given by the Laguerre method of appropriate tuning and

$\mathbb{E}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]=\Theta(t^{-(\alpha+2)})\,.$
(33)

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({\boldsymbol{x}}_{k})-f({\boldsymbol{x}}^{\star})\geq\frac{C}{\sqrt{t}}\,.$

These rates are not found when $\alpha\rightarrow-1$, indicating that the
worst-case is particularly pessimistic in this scenario.

######  Remark 4.7.

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., 2022)

##  5 Experiments

Figure 4: Rates for a synthetic problem, simulating the Marchenko-Pastur
distribution. CG stands for conjugate gradient and GD for gradient descent.
Shades are standard deviation for 8 different samplings of the distribution
and of the random initialization. Note that both tunings of the GCM achieve
performance in function value very close to the one of Conjugate Gradient,
which is optimal for every draw of the problem.

Figure 5: Function values for different methods where the data comes from
CIFAR-10 Inception features (left) and MNIST features (right). The properly
tuned GCM achieves remarkable performance under these non-synthetic spectra.

Figure 6: Empirical spectrum for the covariance matrix of the features. On the
left we considered the MNIST raw features and on the right we considered the
features of an inception net applied on the CIFAR10 dataset.

Figure 7:  Comparison between experiments run on synthetic Beta distribution
and theoretical asymptotic. Shades are standard deviation for 8 different
samplings of the distribution and of the random initialization. Regression
Assymptotic is the linear regression of the last 700 values, in logarithmic
scale . Y-axis is the function value.

We generate random quadratic problems in two different ways. First, by
sampling ${\boldsymbol{H}}={\boldsymbol{X}}{\boldsymbol{X}}^{\top}$ where the
entries of ${\boldsymbol{X}}$ are i.i.d. Gaussian. Since the ESD of
${\boldsymbol{H}}$ converges to the Marchenko-Pastur distribution as the
number of samples and features, go to infinity. This corresponds to a Beta
distribution with parameters $(\tau,\xi)=(1/2,-1/2)$.

The other process we use to generate random quadratic problems is by sampling
$\Lambda\in{\mathbb{R}}^{d}$ from the corresponding Beta distribution and
taking
${\boldsymbol{H}}={\boldsymbol{U}}\operatorname*{\mathbf{diag}}(\Lambda){\boldsymbol{U}}^{\top}$,
where ${\boldsymbol{U}}$ is an independently sampled orthonormal matrix.

We let ${\boldsymbol{x}}^{\star}=\boldsymbol{0}$ and sample
${\boldsymbol{x}}_{0}$ from a standard Gaussian distribution. In all
experiments, we use the problem’s instance largest eigenvalue to calibrate
each method (e.g., gradient descent’s step size is $1/L$).

Our theoretical rates in Theorem 4.5 and Theorem 4.4 respectively for the
Nesterov method and gradient descent are precise under the approximate range
$-1<\xi<0$ as we show in Figure 7. Distributions with higher $\xi$ need many
samples. Otherwise, they behave as strongly convex functions.

The Generalized Chebyshev Method shows large variance for the different
objectives $r_{t}$ in some settings. Consider $\beta^{\star}$ the optimal
$\beta$ value for a given $\xi$. If $\beta<\beta^{\star}$ or $\xi$ is low the
function value $f({\boldsymbol{x}}_{t})$, when seen as random variable w.r.t.
the random initialization is slow to converge with increasing dimension. This
is shown in appendix D.

The GCM with $\beta>\beta^{\star}$ performs corresponding to the theory, and
it is non-asymptotically very close to the performance of $\beta^{\star}$.
High values of $\beta$ also perform very well on non-synthetic data,
suggesting we should use these values in practice.

##  6 Conclusion and further work

In this paper, we’ve established that the asymptotic 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.

### Acknowledgements

Gauthier Gidel is supported by an IVADO grant and a Canada CIFAR AI Chair.

## References

  * Bogdanov & Trevisan (2006) Bogdanov, A. and Trevisan, L.  Average-case complexity.  _arXiv preprint cs/0606037_ , 2006.
  * Brakhage (1987) Brakhage, H.  On ill-posed problems and the method of conjugate gradients.  In _Inverse and ill-posed Problems_. Elsevier, 1987.
  * Fischer (1996) Fischer, B.  _Polynomial based iteration methods for symmetric linear systems_.  SIAM, 1996.
  * Flanders & Shortley (1950) Flanders, D. A. and Shortley, G.  Numerical determination of fundamental modes.  _Journal of Applied Physics_ , 1950.
  * Hestenes et al. (1952) Hestenes, M. R., Stiefel, E., et al.  _Methods of conjugate gradients for solving linear systems_.  NBS Washington, DC, 1952.
  * Marcellán & Álvarez-Nodarse (2001) Marcellán, F. and Álvarez-Nodarse, R.  On the “Favard theorem” and its extensions.  _Journal of computational and applied mathematics_ , 127(1-2):231–254, 2001.
  * Nesterov (2003) Nesterov, Y.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media, 2003.
  * Paquette et al. (2022) Paquette, C., van Merriënboer, B., Paquette, E., and Pedregosa, F.  Halting time is predictable for large models: A universality property and average-case analysis.  _Foundations of Computational Mathematics_ , 2022.
  * Pedregosa & Scieur (2020) Pedregosa, F. and Scieur, D.  Acceleration through spectral density estimation.  In _International Conference on Machine Learning_. PMLR, 2020.
  * Polyak (1964) Polyak, B.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 04, 1964\.
  * Scieur & Pedregosa (2020) Scieur, D. and Pedregosa, F.  Universal average-case optimality of Polyak momentum.  _International Conference on Machine Learning_ , 2020.
  * Szegö (1975) Szegö, G.  Orthogonal polynomials.  In _American Mathematical Society Colloquium Publications_ , 1975\.
  * Van Assche (1995) Van Assche, W.  Weak convergence of orthogonal polynomials.  _Indagationes Mathematicae_ , 1995.

##  Appendix A Proofs of Section 4

See 2.5

###### Proof.

We remark that by the definition of the expected spectral distribution $\mu$
of ${\boldsymbol{H}}$, we have for a polynomial $g$

${\mathbb{E}}_{H}[\operatorname{tr}(g({\boldsymbol{H}}))]=\int
g(\lambda)\mathop{}\\!\mathrm{d}\mu(\lambda).$ (34)

We know that
${\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}=P_{t}({\boldsymbol{H}})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})$.
We can write $\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$ in terms
of a trace and use the independence of ${\boldsymbol{H}}$ and
${\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}$ to connect it to the ESD:

$\displaystyle\mathbb{E}\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle={\mathbb{E}}[\operatorname{tr}(({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}P_{t}({\boldsymbol{H}})^{2}({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}))],$
(35)
$\displaystyle={\mathbb{E}}_{{\boldsymbol{H}},{\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}}[\operatorname{tr}(P_{t}({\boldsymbol{H}})^{2}({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}]$
(36)
$\displaystyle={\mathbb{E}}_{\boldsymbol{H}}\left[\operatorname{tr}(P_{t}({\boldsymbol{H}})^{2}{\mathbb{E}}_{{\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}}[({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}])\right)],$
(37)
$\displaystyle=R^{2}{\mathbb{E}}_{\boldsymbol{H}}[P_{t}(\operatorname{tr}({\boldsymbol{H}}))^{2}]=R^{2}\int
P_{t}(\lambda)^{2}\mathop{}\\!\mathrm{d}\mu(\lambda)\,.$ (38)

For the gradient and function value the reasoning is the same by noticing that

$\displaystyle{\mathbb{E}}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]$
$\displaystyle=\frac{R^{2}}{2}{\mathbb{E}}[\operatorname{tr}(({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}P_{t}({\boldsymbol{H}}){\boldsymbol{H}}P_{t}({\boldsymbol{H}})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}))]$
(39) $\displaystyle={\mathbb{E}}_{\boldsymbol{H}}[\operatorname{tr}((\lambda
P_{t})({\boldsymbol{H}})^{2})]=\frac{R^{2}}{2}\int\lambda
P_{t}(\lambda)^{2}\mathop{}\\!\mathrm{d}\mu(\lambda),$ (40)

where $\lambda P_{t}$ is also a polynomial. As $\nabla
f({\boldsymbol{x}}_{t})={\boldsymbol{H}}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})$.

$\displaystyle{\mathbb{E}}\|\nabla f({\boldsymbol{x}}_{t}))\|^{2}$
$\displaystyle={\mathbb{E}}[\operatorname{tr}(({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star})^{\top}P_{t}({\boldsymbol{H}}){\boldsymbol{H}}^{2}P_{t}({\boldsymbol{H}})({\boldsymbol{x}}_{0}-{\boldsymbol{x}}^{\star}))]$
(41)
$\displaystyle={\mathbb{E}}_{\boldsymbol{H}}[\operatorname{tr}((\lambda^{2}P_{t})({\boldsymbol{H}})^{2})]=R^{2}\int\lambda^{2}P_{t}(\lambda)^{2}\mathop{}\\!\mathrm{d}\mu(\lambda)\,.$
(42)

∎

See 2.6

###### Proof.

This is a direct application of the following lemma w.r.t. to the distribution
$\nu$ defined as
$d\nu(\lambda)=\lambda^{l}\mathop{}\\!\mathrm{d}\mu(\lambda)$.

######  Lemma A.1.

(Fischer, 1996) The residual polynomial of degree $t$ that minimizes $\int
P_{t}^{2}(\lambda)d\nu(\lambda)$ is given by the degree $t$ residual
orthogonal polynomial with respect to the weight function $\lambda
d\nu(\lambda)$.

Let $P_{t}(\lambda;(a_{k}))$ be the polynomial of degree $k$ defined by the
coefficients $a_{k}$ and
$F((a_{k}))=\int\lambda^{l}P_{t}^{2}(\lambda;(a_{k}))\mathop{}\\!\mathrm{d}\mu(\lambda)$.
The restriction $P_{t}(0)=1$ is expressed as $g(a_{k})=a_{0}-1=0$. We then
consider the minimization problem:

$\displaystyle\min_{a_{k}}F(a_{k})$ s.t. $\displaystyle g(a_{k})=0$

We observe that the problem is convex in the variables $a_{k}$ and consider
the lagrangian $\mathcal{L}(a_{k};\alpha)=F(a_{k})-\alpha g(a_{k})$. By the
first-order optimality conditions the derivatives w.r.t. the coefficients
$a_{k},k>0$ of the polynomials should be $0$:

$\displaystyle\frac{d}{da_{k}}\left(\int\lambda^{l}P_{t}^{2}(\lambda)\mathop{}\\!\mathrm{d}\mu(\lambda)\right)$
$\displaystyle=\int\lambda^{l}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}a_{k}}\left(\sum_{k=0}^{\top}a_{k}\lambda^{k}P_{t}(\lambda)\right)\mathop{}\\!\mathrm{d}\mu(\lambda),$
$\displaystyle=2\cdot\left(\int\lambda^{l+k}P_{t}(\lambda)\mathop{}\\!\mathrm{d}\mu(\lambda)\right)=0.$

This means that $P_{t}(\lambda)$ is orthogonal to any polynomial of degree
$t-1$ w.r.t to the intern product
$\langle.,.\rangle_{\lambda^{l+1}\mathop{}\\!\mathrm{d}\mu}$

∎

##  Appendix B GCM and Laguerre method derivation

We will first state two lemmas that allow us to construct the optimal
polynomials.

######  Lemma B.1.

Let $(\tilde{P}_{t})$ be a family polynomials with recurrence

$\tilde{P}_{t}(\lambda)=(\alpha_{t}+\beta_{t}\lambda)\tilde{P}_{t-1}(\lambda)+\gamma_{t}\tilde{P}_{t-2}(\lambda),$

with $\tilde{P}_{0}$ a constant polynomial and $\tilde{P}_{t}(0)\neq 0$ for
all $t$. Then,

$P_{t}(\lambda)=(a_{t}+b_{t}\lambda)P_{t-1}(\lambda)+(1-a_{t})P_{t-2}(\lambda)\,,$
(43)

is the recurrence for
$P_{t}(\lambda)=\tilde{P}_{t}(\lambda)/\tilde{P}_{t}(0)$. With:

$\displaystyle a_{t}$ $\displaystyle=\delta_{t}\alpha_{t},$ (44)
$\displaystyle b_{t}$ $\displaystyle=\delta_{t}\beta_{t},$ (45)
$\displaystyle\delta_{t}$
$\displaystyle=(\alpha_{t}+\gamma_{t}\delta_{t-1})\hskip
14.22636pt(\delta_{0}=0).$ (46)

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 B.2.

Let $(\tilde{P}_{t})$ be a family polynomials following

$\tilde{P}_{t}(\lambda)=(\alpha_{t}+\beta_{t}\lambda)\tilde{P}_{t-1}(\lambda)+\gamma_{t}\tilde{P}_{t-2}(\lambda),$
(47)

and define polynomials $P_{t}$ such that :

$P_{t}(m(\lambda))=\tilde{P}_{t}(\lambda),$

with $m(\lambda)=a\lambda+b$ a non-singular affine transform. Then $P_{t}$
follows a recurrence like in eq. (47), with:

$\displaystyle\alpha_{t}^{\prime}$ $\displaystyle=\alpha_{t}+b\beta_{t},$ (48)
$\displaystyle\beta^{\prime}_{t}$ $\displaystyle=a\beta_{t}\,,$ (49)
$\displaystyle\gamma^{\prime}_{t}$ $\displaystyle=\gamma_{t}\,.$ (50)

###### Proof.

The result follows from considering Eq. (47) with argument $m^{-1}(\lambda)$.
∎

These results are sufficient to obtain the recurrence relation for the
residual polynomial w.r.t $x^{\beta}(L-x)^{\alpha}$. We begin by the standard
Jacobi polynomials, which are orthogonal w.r.t $(1-x)^{\alpha}(1+x)^{\beta}$
and follow a recurrence according to $\alpha_{t},\beta_{t},\gamma_{t}$ below
(Szegö, 1975):

$\displaystyle\alpha_{t}$
$\displaystyle=\frac{(2n+\alpha+\beta)(2n+\alpha+\beta-1)}{2n(n+\alpha+\beta)},$
(51) $\displaystyle\beta_{t}$
$\displaystyle=\frac{(\alpha^{2}-\beta^{2})(2n+\alpha+\beta-1)}{2n(n+\alpha+\beta)(2n+\alpha+\beta-2)},$
(52) $\displaystyle\gamma_{t}$
$\displaystyle=\frac{-2(n+\alpha-1)(n+\beta-1)(2n+\alpha+\beta)}{2n(n+\alpha+\beta)(2n+\alpha+\beta-2)}\,.$
(53)

We then shift the distribution according to $\eta(x)$, and then transform to
the residual ones. We slightly simplify these computations and use remark 2.3
to get Algorithm 1.

We know (Szegö, 1975) that the Laguerre polynomials $L_{t}^{\alpha}$, with
usual normalization, follow the recurrence

$L_{t}^{\alpha}(\lambda)=\left(\frac{2t+\alpha-1}{t}-\frac{1}{t}\lambda\right)L_{t-1}^{\alpha}(\lambda)+\frac{t+\alpha-1}{t}L_{t-2}^{\alpha}(\lambda)\,.$
(54)

As we don’t have to shift the domains, we have only to apply lemma B.1 to get
the Laguerre method. Further, we can get a explicit expression for
$\delta_{t}=\frac{t}{t+\alpha}$ by simplifying the expression.

Algorithm 2 Laguerre($\alpha$)

Inputs: Initial vector ${\boldsymbol{x}}_{0}$, function $f$, distributional
parameter $\alpha$ ${\boldsymbol{x}}_{-1}\leftarrow\textbf{0}$

for $t=1,\ldots,T$ do

${\boldsymbol{x}}_{t}\leftarrow{\boldsymbol{x}}_{t-1}+\frac{t-1}{t+\alpha}({\boldsymbol{x}}_{t-1}-{\boldsymbol{x}}_{t-2})-\frac{1}{t+\alpha}\nabla
f({\boldsymbol{x}}_{t-1})$

end for

##  Appendix C Proofs of section 4

In the following we will consider shifted versions of the spectral
distributions. This shift is the affine function
$m(\lambda):[0,L]\rightarrow[-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 $P_{t}(x)=\tilde{P}_{t}(m(x))$ and
$\mu^{\prime}(x)=\tilde{\mu}^{\prime}(m(x))$, then,

$\int
P_{t}^{2}(x)\mu^{\prime}(x)\mathop{}\\!\mathrm{d}x\propto\int\tilde{P}_{t}^{2}(x)\tilde{\mu}^{\prime}(x)\mathop{}\\!\mathrm{d}x,$
(55)

so all the asymptotics are the same. The Jacobi polynomials
$P^{\alpha,\beta}_{t}$ are those orthogonal w.r.t
$\mathop{}\\!\mathrm{d}\mu(x)=(1-x)^{\alpha}(1+x)^{\beta}\mathop{}\\!\mathrm{d}x$.
Most works use the normalization
$\tilde{P}_{t}^{\alpha,\beta}(-1)=(-1)^{t}\binom{t+\beta}{t}$. We will write
$\tilde{P}^{\alpha,\beta}_{t}$ for this normalization and
$P^{\alpha,\beta}_{t}$ for the residual polynomials See 4.1

###### Proof.

We will prove that for any $\alpha$ and $\beta$, $\xi,\tau>-1$, $l>0$ and
distributions $\nu_{\tau,\xi-l}$, we have

$\int P_{t}^{\alpha,\beta}(x)^{2}x^{l}d\nu_{\tau,\xi-l}(x)\sim
L^{l}C^{\alpha,\beta}_{\nu}\left\\{\begin{array}[]{ll}t^{-1-2\beta}&\mbox{if
}\alpha<\tau+1/2\text{ and }\beta<\xi+1/2,\\\ t^{-2(\xi+1)}\log t&\mbox{if
}\alpha=\tau+1/2\text{ and }\beta=\xi+1/2,\\\
t^{2(\max\\{\alpha-\beta-\tau,-\xi\\}-1)}&\mbox{if }\alpha>\tau+1/2\text{ or
}\beta>\xi+1/2.\end{array}\right.$

We will first show this result for the beta weights, then show that
distributions with the same concentration behave similarly.
The normalization of $\tilde{P}^{\alpha,\beta}_{t}$ is s.t. (Szegö, 1975)
(4.3.3):

$\int_{-1}^{1}\tilde{P}_{t}^{\alpha,\beta}(x)(1-x)^{\alpha}(1+x)^{\beta}\mathop{}\\!\mathrm{d}x=\frac{2^{\alpha+\beta+1}}{2t+\alpha+\beta+1}\frac{\Gamma(t+\alpha+1)\Gamma(t+\beta+1)}{\Gamma(t+1)\Gamma(t+\alpha+\beta+1)}=\Theta(t^{-1})\,.$
(56)

Further, the residual polynomials are s.t.
$|P^{\alpha,\beta}_{t}|=\Theta(t^{-\beta})|\tilde{P}^{\alpha,\beta}_{t}|$,
from the definition of the classical normalization.
We state the result (Exercise 91, Generalisation of 7.34.1) from (Szegö,
1975):

######  Lemma C.1 (Szegö (1975)).

We have,

$\displaystyle\int_{0}^{1}(1-x)^{\tau}P_{t}^{\alpha,\beta}(x)^{2}dx\sim\Theta(h_{\tau}^{\alpha}),$
(57) $\displaystyle
h_{\tau}^{\alpha}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\\{\begin{array}[]{ll}t^{2(\alpha-\tau-1)}&\mbox{if
}\alpha>\tau+1/2,\\\ t^{-1}\log t&\mbox{if }\alpha=\tau+1/2,\\\
t^{-1}&\mbox{if }\alpha<\tau+1/2\,.\end{array}\right.$ (61)

Noting that
$\tilde{P}^{\alpha,\beta}_{t}(x)=(-1)^{\top}\tilde{P}_{t}^{\beta,\alpha}(-x)$,
we can write:

$\int_{-1}^{1}\tilde{P}_{t}(x)^{2}(1-x)^{\tau}(1+x)^{\xi}dx=\Theta\left(\int_{0}^{1}(1-x)^{\tau}|\tilde{P}_{t}^{\alpha,\beta}(x)|^{2}dx\right)+\Theta\left(\int_{0}^{1}(1-x)^{\xi}|\tilde{P}_{t}^{\beta,\alpha}(x)|^{2}dx\right).$
(62)

We can then show our result for
$\mathop{}\\!\mathrm{d}\nu_{\tau,\xi-l}(x)=x^{\xi-l}(L-x)^{\alpha}\mathop{}\\!\mathrm{d}x$
by carefully considering each of the cases on Lemma C.1 and the maximum of
each term in eq. 62, and an added $t^{-2\beta}$ from the different
normalization. With this, we have the wanted result for the Beta weights
It remains to show:

$\int_{0}^{1}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi}(x)=\Theta\left(\int_{0}^{1}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}dx\right).$
(63)

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 C.2 (Van Assche (1995)).

Let $\mu$ be a measure and $(p_{t})$ a family of associated orthonormal
polynomials such that $p_{t}$ follow the recurrence:

$xp_{t}(x)=a_{t}p_{t+1}(x)+b_{t}p_{t}(x)+a_{t-1}p_{t-1}(x),$

and $a_{t},b_{t}$ converge respectively to $a,b$. Then for any $f$ continuous
and bounded we have:

$\int
f(x)p_{t}^{2}(x)\mathop{}\\!\mathrm{d}\mu(x)\rightarrow\frac{1}{\pi}\int_{-1}^{1}\frac{f(x)}{\sqrt{1-x^{2}}}dx\,.$
(64)

Let $\epsilon$ such that

$x\geq
1-\epsilon\Rightarrow|\mathop{}\\!\mathrm{d}\nu_{\tau,\xi}-A(1-x)^{\tau}|\leq
B(1-x)^{\tau}.$ (65)

We observe that for $0<x<1-\epsilon$,
$f(x)=\frac{d\nu_{\tau,\xi}}{(1-x)^{\alpha}(1+x)^{\beta}}$ is bounded.
We get from an application of 64, and the observation that
$\tilde{P}_{t}^{\alpha,\beta}=\mathcal{N}_{t}p_{t}^{\alpha,\beta}$, with
constants $\mathcal{N}_{t}=\Theta(t^{-1/2})$:

$\displaystyle\underbrace{\int_{0}^{1}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}x}_{\Theta(h_{\tau}^{\alpha})}$
$\displaystyle=\underbrace{\int_{0}^{1-\epsilon}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}x}_{\Theta(t^{-1})}+\int_{1-\epsilon}^{1}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}x\Rightarrow$
(66)
$\displaystyle\int_{1-\epsilon}^{1}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}x=\Theta(h_{\tau}^{\alpha}),$
(67)
$\begin{split}\int_{0}^{1}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi}(x)&=\underbrace{\int_{0}^{1-\epsilon}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}f(x)(1-x)^{\alpha}(1+x)^{\beta}\mathop{}\\!\mathrm{d}x)}_{\Theta(t^{-1})}\\\
&+\Theta\left(\underbrace{\int_{1-\epsilon}^{1}(1-x)^{\tau}\tilde{P}_{t}^{\alpha,\beta}(x)^{2}\mathop{}\\!\mathrm{d}x}_{\Theta(h_{\tau}^{\alpha})}\right).\end{split}$
(68)

∎

See 4.2

###### Proof.

We will prove that for $\tau,\xi>-1$ If $\alpha=\tau$ and $\beta=\xi+l+1$
(i.e., are optimal), the rate of convergence reads

$\min_{P_{t}(0)=1}\int
P_{t}^{2}(\lambda)\lambda^{l}d\nu(\lambda)=\Theta\left(\int_{0}^{l}\tilde{P}_{t}^{\alpha,\beta}(\lambda)^{2}(L-\lambda)^{\tau}\lambda^{\xi+l}\mathop{}\\!\mathrm{d}\lambda\right)=\Theta(t^{-2(\xi+l+1)}).$
(69)

Showing the second equality is easy by considering theorem 4.1, and that is
further the minimum asymptotic rate for the Beta distribution
$\mu_{\tau,\xi}$.
By setting $p_{t}^{\nu}$ and $P^{\nu}_{t}=\frac{p_{t}^{\nu}}{p^{\nu}_{t}(0)}$
the optimal orthonormal and residual and polynomials w.r.t. $\nu$ we show that
$P^{\nu}_{t}$ must have the same rate on $\mu_{\tau,\xi}$ as it does on $\nu$,
thus the optimal rate of $\nu$ cannot be lower than the optimal rate of
$\mu_{\tau,\xi}$. Indeed, setting $\epsilon_{1},\epsilon_{2}$ as in eq. 65:

$\displaystyle\int_{1-\epsilon_{2}}^{1}P_{t}^{\nu}(x)^{2}d\nu(x)$
$\displaystyle=\Theta\left(\int_{1-\epsilon_{2}}^{1}P_{t}^{\nu}(x)^{2}\mathop{}\\!\mathrm{d}\mu_{\tau,\xi}(x)\right),$
(70) $\displaystyle\int_{-1}^{-1+\epsilon_{1}}P_{t}^{\nu}(x)^{2}d\nu(x)$
$\displaystyle=\Theta\left(\int_{-1}^{-1+\epsilon_{1}}P_{t}^{\nu}(x)^{2}\mathop{}\\!\mathrm{d}\mu_{\tau,\xi}(x)\right),$
(71)
$\displaystyle\int_{-1+\epsilon_{1}}^{1-\epsilon_{2}}P_{t}^{\nu}(x)^{2}d\nu(x)$
$\displaystyle=\Theta\left(\int_{-1+\epsilon_{1}}^{1-\epsilon_{2}}P_{t}^{\nu}(x)^{2}\mathop{}\\!\mathrm{d}\mu_{\tau,\xi}(x)\right)=\Theta\left(\frac{1}{p_{t}^{\nu}(-1)^{2}}\right),$
(72)

where the first two equations come from the fact that
$\nu=\Theta(\mu_{\tau,\xi})$ near $-1$ and $1$ and the third from lemma 64.
This effectively lower bounds the rates on $\nu$ because the rates of
$P_{t}^{\nu}$ on $\mu_{\tau,\xi}$ can’t be lower than $-2(\xi+l+1)$. ∎

See 4.3

###### Proof.

We will prove that:
$\sup_{x\in[0,L]}x^{l}P_{t}^{\alpha,\beta}(x)^{2}=O(L^{l}t^{v(\alpha,\beta,l)})$,
where:

$v(\alpha,\beta,l)=\left\\{\begin{array}[]{cc}2(\alpha-\beta)&\text{if}\hskip
14.22636pt\alpha>\beta-l\\\ -1-2\beta&\text{if}\hskip
14.22636pt\alpha\leq\beta-l\hskip 14.22636pt\beta\leq l-\frac{1}{2},\\\
-2l,&\text{if}\hskip 14.22636pt\alpha\leq\beta-l\hskip 14.22636pt\beta\geq
l-\frac{1}{2}.\end{array}\right.$ (73)

From (Szegö, 1975), Theorem 7.32.2, if $\theta<\frac{\pi}{2}$:

$\tilde{P}_{t}^{\alpha,\beta}(\cos\theta)=\left\\{\begin{array}[]{cc}O(t^{-1/2})&\text{if}\hskip
14.22636pt\alpha<-\frac{1}{2},\\\ O(t^{\alpha})&\text{if}\hskip
14.22636pt\alpha\geq-\frac{1}{2}\land 0\leq\theta\leq ct^{-1},\\\
\theta^{-\alpha-1/2}O(t^{-1/2})&\text{if}\hskip
14.22636pt\alpha\geq-\frac{1}{2}\land\theta>ct^{-1}.\end{array}\right.$ (74)

We observe that, from the symmetry of the Jacobi polynomials:

$\sup_{x\in[0,L]}x^{l}P_{t}^{\alpha,\beta}(x)^{2}=\Theta\left(\max\left\\{\sup_{x\in[0,1]}x\tilde{P}_{t}^{\alpha,\beta}(x)^{2},\sup_{x\in[0,1]}(1-x)^{l}\tilde{P}_{t}^{\beta,\alpha}(x)^{2}\right\\}\right).$
(75)

The $(1-x)^{l}$ term, corresponds to $(2\sin(\frac{\theta}{2}))^{2l}$ in the
variable $\theta$, which is $O(\theta^{2l})$. The rest follows from carefully
considering the expressions given by eq. 74. ∎

See 4.4

###### Proof.

We will prove:

$\int_{0}^{1}P_{t}^{\text{Nes}}(\lambda)^{2}\lambda^{l}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi-l}\sim
C^{\prime}_{\nu}\Big{\\{}\begin{array}[]{ll}t^{-2(\xi+1)}&\mbox{if
}0<\xi<1/2,\\\ t^{-3}\log t&\mbox{if }\xi=1/2,\\\ t^{-(\xi+5/2)}&\mbox{if
}\xi>1/2.\end{array}$ (76)

(Paquette et al., 2022) has shown that the nesterov polynomials $P_{t}$ are
asymptotically, in $t$:

$P_{t}(\lambda)\sim\frac{2J_{1}(t\sqrt{\alpha\lambda})}{t\sqrt{\alpha\lambda}}e^{-\alpha\lambda
t/2},$ (77)

in the sense that:

$\int_{0}^{1}u^{l}\left[\tilde{P_{t}^{2}(u)}-\frac{4J_{1}^{2}(t\sqrt{u})}{t^{2}u}e^{-ut}\right]\mathop{}\\!\mathrm{d}\mu_{MP}(u)=O(t^{-(l+25/12))}).$
(78)

The arguments can be easily used to show that such an integral is
$O(t^{-(\alpha+l+31/12)})$ when evaluated w.r.t. a general
$\mathop{}\\!\mathrm{d}\mu$ s.t $\mu^{\prime}=\Theta(\lambda^{\alpha})$ near
$0$.
We can thus consider our integral of interest substituting
$P_{t}^{\text{Nes}}$ by it’s Bessel asymptotic and dividing it into three
regions, i.e.
$[0,1]=[0,\frac{\epsilon}{t}]\cup[\frac{\epsilon}{t},\frac{\epsilon}{\sqrt{t}}]\cup[\frac{\epsilon}{\sqrt{t}},1]$
corresponding to two different regimes for the Bessel function. The first
region will give us the asymptotic and we will bound the others.
We consider first, the "middle region",i .e. for some $\epsilon>0$:

$\int_{\frac{\epsilon}{t}}^{\frac{\epsilon}{\sqrt{t}}}u^{\xi}\frac{4J_{1}^{2}(t\sqrt{u})}{t^{2}u}e^{-ut}\mathop{}\\!\mathrm{d}u.$
(79)

We note the asymptotic for $J_{1}^{2}$:

$J_{1}^{2}(\sqrt{tv})\sim\frac{1}{\pi\sqrt{tv}}(1+\cos(2\sqrt{tv}+2\gamma)).$
(80)

Doing the change of variable $v=tu$, and identifying the upper limit of the
interval, which is $\epsilon t^{1/2}$, to $\infty$:

$\displaystyle\int_{\frac{\epsilon}{t}}^{\frac{\epsilon}{\sqrt{t}}}u^{\xi}\frac{4J_{1}^{2}(t\sqrt{u})}{t^{2}u}e^{-ut}\mathop{}\\!\mathrm{d}u$
$\displaystyle=\Theta\left(t^{-2-\xi}\int_{\epsilon}^{\infty}v^{\xi-1}J_{1}^{2}(\sqrt{tv})e^{-v}\mathop{}\\!\mathrm{d}v\right),$
(81)
$\displaystyle=\Theta\left(t^{-2-\xi}\int_{\epsilon}^{\infty}v^{\xi-1}\frac{1}{\pi\sqrt{tv}}e^{-v}\mathop{}\\!\mathrm{d}v\right),$
(82)
$\displaystyle=\Theta\left(t^{-\frac{5}{2}-\xi}\underbrace{\int_{\epsilon}^{\infty}v^{\xi-\frac{3}{2}}\frac{1}{\pi\sqrt{tv}}e^{-v}\mathop{}\\!\mathrm{d}v}_{\Gamma(\xi-\frac{1}{2},\epsilon)}\right).$
(83)

Where, from the Riemann-Lebesgue, the cosine term goes to $0$ lemma and
$\Gamma$ is the incomplete Gamma function.
The mass from the rightmost region ,i.e, the term corresponding to the
interval $[\epsilon t^{-1/2},1]$ is exponentially small. Indeed, because of
the exponential $e^{-ut}$ it is $O(e^{-\epsilon\sqrt{t}})$. Lastly, we have
for the $[0,\frac{\epsilon}{t}]$ region, doing the change of variables
$v=t^{2}u$:

$\int_{0}^{\frac{\epsilon}{t}}u^{\xi}\frac{4J_{1}^{2}(t\sqrt{u})}{t^{2}u}e^{-ut}\mathop{}\\!\mathrm{d}u=\Theta\left(t^{-2(\xi+1)}\int_{0}^{t\epsilon}v^{\xi}\frac{J_{1}^{2}(\sqrt{v})}{v}e^{-\frac{v}{t}}\mathop{}\\!\mathrm{d}v\right).$
(84)

And the $e^{\frac{-v}{t}}$ term is $\Theta(1)$. We have the following Bessel
asymptotics:

$\displaystyle\frac{J_{1}^{2}(\sqrt{v})}{v}$
$\displaystyle\sim\frac{1}{4},\hskip 56.9055ptv\rightarrow 0,$ (85)
$\displaystyle\frac{J_{1}^{2}(\sqrt{v})}{v}$ $\displaystyle=O(v^{-3/2}),\hskip
28.45274ptv\rightarrow\infty,$ (86)

so we divide this integral aswell:

$\displaystyle
t^{-2(\xi+1)}\int_{1}^{t\epsilon}v^{\xi}\frac{J_{1}^{2}(\sqrt{v})}{v}e^{-\frac{v}{t}}\mathop{}\\!\mathrm{d}v,$
$\displaystyle=\Theta\left(t^{-2(\xi+1)}\int_{\epsilon}^{t\epsilon}v^{\xi-\frac{3}{2}}\mathop{}\\!\mathrm{d}v\right)=\Theta\left(I_{\xi}(t)t^{-\xi-\frac{5}{2}}\right)$
(87) $\displaystyle
t^{-2(\xi+1)}\int_{0}^{1}v^{\xi}\frac{J_{1}^{2}(\sqrt{v})}{v}e^{-\frac{v}{t}}\mathop{}\\!\mathrm{d}v$
$\displaystyle=\Theta\left(t^{-2(\xi+1)}\int_{0}\epsilon^{1}v^{\xi}\mathop{}\\!\mathrm{d}v\right)=\Theta\left(t^{-2(\xi+1)}\right),$
(88)

where $I_{\xi}(t)=\log t$ if $\xi=\frac{1}{2}$ and $1$ otherwise.
The Nesterov rate is then $I_{\xi}(t)t^{-\xi-\frac{5}{2}}$ if
$\xi\geq\frac{1}{2}$ and $t^{-2(\xi+1)}$ if $0<\xi<\frac{1}{2}$ ∎

See 4.5

###### Proof.

Considering that $P_{t}^{\text{GD}}(\lambda)=(1-\frac{\lambda}{L})^{\top}$ we
will prove :

$\int_{0}^{1}(1-\lambda)^{2t}\lambda^{l}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi-l}=\Theta(t^{-(\xi+l+1)}).$
(89)

We know, for the Beta weights, that:

$\int_{0}^{1}(1-\lambda)^{2t+\tau}\lambda^{\xi+l}\mathop{}\\!\mathrm{d}\lambda=\frac{\Gamma(l+\xi+1)\Gamma(2t+\tau+1)}{\Gamma(2t+l+\xi+\tau+2)}=\Theta(t^{-(\xi+l+1)}).$
(90)

We can identify this asymptotic to the interval $\int_{0}^{\epsilon}$ for any
$\epsilon$ because:

$\int_{\epsilon}^{1}(1-\lambda)^{2t+\tau}\lambda^{\xi+l}\mathop{}\\!\mathrm{d}\lambda=\mathcal{O}((1-\epsilon)^{2t}),$
(91)

then:

$\displaystyle\int_{\epsilon}^{1}(1-\lambda)^{2t}\lambda^{l}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi-l}$
$\displaystyle=\mathcal{O}((1-\epsilon)^{2t}),$ (92)
$\displaystyle\int_{0}^{\epsilon}(1-\lambda)^{2t}\lambda^{l}\mathop{}\\!\mathrm{d}\nu_{\tau,\xi-l}$
$\displaystyle=\Theta\left(\int_{0}^{\epsilon}(1-\lambda)^{2t+\tau}\lambda^{\xi+l}\mathop{}\\!\mathrm{d}\lambda\right)=\Theta(t^{-(\xi+l+1)}).$
(93)

∎

See 4.6

###### Proof.

Let $L_{t}^{\alpha}$ be the Laguerre polynomials with the usual normalization
(Szegö, 1975):

$\int
L_{t}^{\alpha}(x)^{2}\mathop{}\\!\mathrm{d}\mu_{\alpha}(x)=L_{t}^{\alpha}(0)=\binom{n+\alpha}{n}$
(94)

Further [(Szegö, 1975) (5.1.13)]]:

$\sum_{k=0}^{T}L_{t}^{\alpha}(x)=L_{t}^{\alpha+1}(x).$ (95)

Thus, letting $P_{t}^{\alpha}$ be the residual Laguerre polynomial, we
consider:

$\begin{split}{\mathbb{E}}[f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}^{\star})]&=\int
P_{t}^{\alpha+2}(\lambda)^{2}\mathop{}\\!\mathrm{d}\mu_{\alpha+1}(\lambda)=\binom{t+\alpha+2}{t}^{-2}\int
L_{t}^{\alpha+2}\mathop{}\\!\mathrm{d}\mu_{\alpha+1}(\lambda),\\\
&=\binom{t+\alpha+2}{t}^{-2}\sum_{k=0}^{t}\left[\int
L_{k}^{\alpha+1}(\lambda)\mathop{}\\!\mathrm{d}\mu_{\alpha+1}(\lambda)\right],\\\
&=\binom{t+\alpha+2}{t}^{-2}\sum_{k=0}^{\top}\binom{k+\alpha+1}{k}=\binom{t+\alpha+2}{t}^{-2}\binom{t+\alpha+2}{t},\\\
&=\binom{t+\alpha+2}{t}^{-1}=\Theta(t^{-(\alpha+2)}).\end{split}$ (96)

∎

##  Appendix D Additional Experiments

Figure 8: Function value for $\beta<\beta^{\star}$ for different dimensions of
the problem

We observed that in the regimes $\beta<\beta^{\star}$ and $\xi\approx-1$ the
objectives exhibit much higher variance than outside these regimes, where they
concentrate tightly around the expected value. Despite this the objectives
still converge with increasing dimension.