Acceleration Methods

Chapter 1 Introduction

Optimization methods are a core component of the modern numerical toolkit. In many cases, iterative algorithms for solving convex optimization problems have reached a level of efficiency and reliability comparable to that of advanced linear algebra routines. This is largely true for medium scale-problems where interior point methods reign supreme, but less so for large-scale problems where the complexity of first-order methods is not as well understood and efficiency remains a concern.

The situation has improved markedly in recent years, driven in particular by the emergence of a number of applications in statistics, machine learning, and signal processing. Building on Nesterov’s path-breaking algorithm from the 80’s, several accelerated methods and numerical schemes have been developed that both improve the efficiency of optimization algorithms and refine their complexity bounds. Our objective in this monograph is to cover these recent developments using a few master templates.

The methods described in this manuscript can be arranged in roughly two categories. The first, stemming from the work of Nest83, produces variants of the gradient method with accelerated worst-case convergence rates that are provably optimal under classical regularity assumptions. The second uses outer iteration (a.k.a. nested) schemes to speed up convergence. In this second setting, accelerated schemes run both an inner loop and an outer loop, with the inner iterations being solved by classical optimization methods, and the outer loop containing the acceleration mechanism.

Ever since the original algorithm by Nest83, the acceleration phenomenon was regarded as somewhat of a mystery. While accelerated gradient methods can be seen as iteratively building a model for the function and using it to guide gradient computations, the argument is essentially algebraic and is simply an effective exploitation of regularity assumptions. This approach of collecting inequalities induced by regularity assumptions and cleverly chaining them to prove convergence was also used in e.g., [Beck09], to produce an optimal proximal gradient method. There too, however, the proof yielded little evidence as to why the method is actually faster.

Fortunately, we are now better equipped to push the proof mechanisms much further. Recent advances in the programmatic design of optimization algorithms allow us to design and analyze algorithms by following a more principled approach. In particular, the performance estimation approach, pioneered by Dror14, can be used to design optimal methods from scratch, selecting algorithmic parameters to optimize worst-case performance guarantees [Dror14, kim2016optimized]. Primal dual optimality conditions on the design problem then provide a blueprint for the accelerated algorithm structure and for its convergence proof.

Using this framework, acceleration is no longer a mystery: it is the main objective in the design of the algorithm. We recover the usual “soup of regularity inequalities” that forms the template of classical convergence proofs, but the optimality conditions of the design problem explicitly produce a method that optimizes the convergence guarantee. In this monograph, we cover accelerated first-order methods using this systematic template and describe a number of convergence proofs for classical variants of the accelerated gradient method, such as those of Nesterov (Nest83, Nest03a), Beck09, tseng2008accelerated as well as more recent ones [kim2016optimized].

The second category of acceleration techniques that we cover in this monograph is composed of outer iteration schemes, in which classical optimization algorithms are used as a black-box in the inner loop and acceleration is produced by an argument in the outer loop. We describe three acceleration results of this type.

The first scheme is based on nonlinear acceleration techniques. Based on arguments dating back to [Aitk27, Wynn56, Ande87], these techniques use a weighted average of iterates to extrapolate a better candidate solution than the last iterate. We begin by describing the Chebyshev method for solving quadratic problems, which interestingly qualifies both as a gradient method and as an outer iteration scheme. It takes its name from the use of Chebyshev polynomial coefficients to approximately minimize the gradient at the extrapolated solution. The argument can be extended to non-quadratic optimization problems provided the extrapolation procedure is regularized.

The second scheme, due to [guler1992new, monteiro2013accelerated, lin2015universal] relies on a conceptual accelerated proximal point algorithm, and uses classical iterative methods to approximate the proximal point in an inner loop. In particular, this framework produces accelerated gradient methods (in the same sense as Nesterov’s acceleration) when the approximate proximal points are computed using linearly converging gradient-based optimization methods, taking advantage of the fact that the inner problems are always strongly convex.

Finally, we describe restart schemes. These techniques exploit regularity properties called Hölderian error bounds, which extend strong convexity properties near the optimum and hold almost generically, to improve the convergence rates of most first-order methods. The parameters of the Hölderian error bounds are usually unknown, but the restart schemes are robust: that is, they are adaptive to the Hölderian parameters and their empirical performance is excellent on problems with reasonable precision targets.

We present a few convergence acceleration techniques that are particularly relevant in the context of (first-order) convex optimization. Our summary includes our own points of view on the topic and is focused on techniques that have received substantial attention since the early 2000’s, although some of the underlying ideas are much older. We do not pretend to be exhaustive, and we are aware that valuable references might not appear below.

The sections can be read nearly independently. However, we believe the insights of some sections can benefit the understanding of others. In particular, Chebyshev acceleration (Section 2) and nonlinear acceleration (Section 3) are clearly complementary readings. Similarly, Chebyshev acceleration (Section 2) and Nesterov acceleration (Section 4), Nesterov acceleration (Section 4) and proximal acceleration (Section 5), as well as Nesterov acceleration (Section 4) and restart schemes (Section 6) certainly belong together.

This monograph is not meant to be a general-purpose manuscript on convex optimization, for which we refer the reader to the now classical references [boyd2004convexopt, bonnans2006numerical, nocedal2006numerical]. Other directly related references are provided in the text.

We assume the reader to have a working knowledge of base linear algebra and convex analysis (such as of subdifferentials), as we do not detail the corresponding technical details while building on them. Classical references on the latter include [Rock70, rockafellar2009variational, hiriart2013convex].

Chapter 2 Chebyshev Acceleration

While “Chebyshev polynomials are everywhere dense in numerical analysis,” we would like to argue here that Chebyshev polynomials also provide one of the most direct and intuitive explanations for acceleration arguments in first-order methods. That is, one can form linear combinations of past gradients for optimizing a worst-case guarantee on the distance to an optimal solution. In quadratic optimization, these linear combinations emerge from a Chebyshev minimization problem, whose solution can also be computed iteratively, thereby yielding an algorithm called the Chebyshev method [Nemic84]. The Chebyshev method traces its roots to at least Flan50, who credit Tuckey and Grosch. Its recurrence matches asymptotically the one of the heavy-ball method and is detailed below.

In this section, we demonstrate basic acceleration results on quadratic minimization problems. In such problems, optimal points are the solutions of a linear system, and the basic gradient method can be seen as a simple iterative solver for this linear system. In this context, acceleration methods can be obtained using a classical argument involving Chebyshev polynomials.

Analyzing this simple scenario is useful in two ways. First, recursive formulations of the Chebyshev argument yield a basic algorithmic template for designing accelerated methods and provide first approach to their structures, such as the presence of a momentum term. Second, the arguments are robust to perturbations of the quadratic function $f$ and hence apply in more generic contexts. This property enables acceleration in a wider range of applications, which we cover later in the Section 3 and Section 4.

For now, consider the following unconstrained quadratic minimization problem

and calling $x_{\star}$ the optimum of problem (2.1) (satisfying $\textbf{H}x_{\star}=b$) yields

This means the iterates of gradient descent $x_{k}-x_{\star}$ can be computed from $x_{0}-x_{\star}$ via $x_{k}-x_{\star}=P_{k}^{\text{Grad}}(\textbf{H})(x_{0}-x_{\star})$ using the matrix polynomial

Suppose we set the step size $\gamma$ to ensure

Because the matrix H is symmetric and hence diagonalizable in an orthogonal basis, given $\gamma>0$, we obtain

To get the best possible worst-case convergence rate, we now minimize this quantity in $\gamma$ by solving

The optimal step size is obtained when both terms in the max are equal, reaching:

Denoting by $\kappa\triangleq\tfrac{L}{\mu}\geq 1$ the condition number of the function $f$, the bound in (2.4) finally becomes

which is a worst-case guarantee for the gradient method when minimizing smooth strongly convex quadratic functions.

2 Optimal Methods and Minimax Polynomials

In Equation (2.4) above, we saw that the worst-case convergence rate of the gradient method on quadratic functions can be controlled by the spectral norm of a matrix polynomial. Figure 2.1 plots the polynomial $P_{k}^{\text{Grad}}$ for several degrees $k$. We can extend this reasoning further to produce methods with accelerated worst-case convergence guarantees.

The bounds derived above for gradient descent can be extended to a broader class of first-order methods for quadratic optimization. We consider first-order algorithms in which each iterate belongs to the span of previous gradients, i.e.

and show that the iterates can be written using matrix polynomials as in (2.3) above.

for all $k=0,1,\ldots$, if and only if the errors $\{x_{k}-x_{\star}\}_{k=0,1,\ldots}$ can be written as

for all $k=0,1\ldots$, for some sequence of polynomials $\{P_{k}\}_{k=0,1,\ldots}$ with $P_{k}$ of degree at most ${k}$ and $P_{k}(0)=1$.

Since $\nabla f(x)$ is the gradient of a quadratic function, it reads

for any $x_{\star}$ satisfying $\textbf{H}x_{\star}=b$, where H is symmetric. We have

We now show recursively that $x_{k}-x_{\star}=P_{k}(\textbf{H})(x_{0}-x_{\star})$, where $P_{k}$ is a residual polynomial of degree at most $k$. Our assumption about the iterates (2.9) implies that, for some sequence of coefficients $\{\alpha_{i}^{(k+1)}\}_{i=0,\ldots,k}$,

Assuming recursively that (2.10) holds for all indices $i\leq k$,

Then, by writing $P_{k+1}(x)=1+x\sum_{i=0}^{k}\alpha_{i}^{(k+1)}P_{i}(x)$, we have

with $P_{k+1}(0)=1$ and $\deg(P_{k+1})\leq k+1$. Since the proof is a sequence of equalities, the equivalence readily follows.

Given a class $\mathcal{M}$ of problem matrices H, Proposition 2.2.1 provides a way to design algorithms. Indeed, we can extract a first-order method from a sequence of polynomials $\{P_{k}\}_{k=0,\ldots,N}$. We can therefore use tools from approximation theory to find optimal polynomials and extract corresponding methods from them. Given a matrix class $\mathcal{M}$, this involves minimizing the worst-case convergence bound over $\textbf{H}\in\mathcal{M}$ by solving

where ${P}_{k}$ is the set of polynomials of degree at most $k$. The polynomial $P^{*}_{k}$ is an optimal polynomial for $\mathcal{M}$ and yields a (worst-case) optimal algorithm for the class $\mathcal{M}$. In terms of the notation used in the proof of Proposition 2.2.1, we are by construction looking for coefficients $\{\alpha_{i}^{(j)}\}_{i}$ depending only on the problem class $\mathcal{M}$, but not on a specific instance of H.

3 The Chebyshev Method

In the case where $\mathcal{M}$ is the set of positive definite matrices with a bounded spectrum, namely

the optimal polynomial can be found by solving

Polynomials that solve (2.12) are derived from Chebyshev polynomials of the first kind in approximation theory and can be formed explicitly to produce an optimal algorithm called the Chebyshev method. This section describes this method and provides its corresponding worst-case convergence guarantees.

We now explicitly introduce the Chebyshev polynomials. A more complete treatment of these polynomials is available in, e.g., mason2002chebyshev. Chebyshev polynomials of the first kind are defined recursively as follows

There exists a compact explicit solution for Chebyshev polynomials that involves trigonometric functions:

It is possible to show that Chebyshev polynomials satisfy the minimax property

where a monic polynomial is a polynomial whose coefficient associated with the highest power is equal to one. From this minimax definition, that defines the minimal polynomial over $$, we can transform it by shifting it to the interval$[\mu,L]$, then rescaling it to obtain a polynomial such that$P(0)=1$. More precisely, using a simple linear mapping from$[\mu,\,L]$to$[-1,\,1]$,

where we have enforced the normalization constraint $C^{[\mu,L]}_{k}(0)=1$. Under these transformations, the shifted Chebyshev polynomials keep some minimax property, and can be shown to be solutions to (2.12).

3.2 Chebyshev Algorithm

The following recursion follows from (2.13) together with (2.15) and a few simplifications:

where $\delta_{1}=\frac{L-\mu}{L+\mu}$ and

The sequence $\delta_{k}$ ensures $C_{k}^{[\mu,L]}(0)=1$. We present this recursion for simplicity, but one should note that it might present some numerical stability issues in practice. There exist numerically more stable first-order methods based on Chebyshev polynomials, see, e.g., [gutknecht2002chebyshev, Algorithm 1]. We plot $C_{1}^{[\mu,L]}(x)$, $C_{3}^{[\mu,L]}(x)$ and $C_{5}^{[\mu,L]}(x)$ in Figure 2.2 for illustration.

Let us quickly detail how to arrive to (2.3.2) using (2.13) and (2.15). We first expand \mathcal{T}_{k}\big{(}t^{[\mu,L]}(x)\big{)} using (2.13),

Since C_{k}^{[\mu,L]}(x)=\frac{\mathcal{T}_{k}\big{(}t^{[\mu,L]}(x)\big{)}}{\mathcal{T}_{k}\big{(}t^{[\mu,L]}(0)\big{)}}, we substitute \mathcal{T}_{k}\big{(}t^{[\mu,L]}(x)\big{)} by \mathcal{T}_{k}\big{(}t^{[\mu,L]}(0))\cdot C_{k}^{[\mu,L]}(x) in the equation above and obtain

For obtaining a simple recursion on $\delta_{k}$, note that $C_{k}^{[\mu,L]}(0)=1$ for all $k\geq 0$ by construction. It follows that

3.3 Chebyshev and Polyak’s Heavy-Ball Methods

We now present the resulting algorithm, called Chebyshev semi-iterative method [Golu61]. We define iterates using $C_{k}^{[\mu,L]}(x)$ as follows,

we can simplify away $x_{\star}$ to get the following recursion:

which describes iterates of the Chebyshev method. We summarize it as Algorithm 2. By construction, the Chebyshev method is a worst-case optimal first-order method for minimizing quadratics whose spectrum lies in $[\mu,L]$. Surprisingly, its iteration structure is simple and somewhat intuitive: it involves a gradient step with variable step size ${4\delta_{k}}/({L-\mu})$, combined with a variable momentum term.

Perhaps more surprisingly, the Chebyshev method has a stationary regime that is even simpler. Indeed, when $k\rightarrow\infty$, the coefficients of the recursion from Algorithm 2 converge to the ones of Polyak’s heavy-ball method,

To see this, it suffices to compute the limit of $\delta_{k}$, written as $\delta_{\infty}$, by solving

We obtain Polyak’s heavy-ball method by replacing $\delta_{k}$ with $\delta_{\infty}$ in Algorithm 2.

3.4 Worst-case Convergence Bounds

The shifted Chebyshev polynomials are solutions to (2.12). Therefore, using the same trick as for gradient descent, we can obtain the following worst-case bound for the Chebyshev method

The maximum value is determined by evaluating the polynomial at one of the extremities of the interval [mason2002chebyshev, Chapter 2] (see also Figure 2.2), i.e,

We obtain the worst-case convergence guarantee of the Chebyshev method, as stated in the following theorem.

After plugging this result into the cosh, we get

It may be difficult to compare the convergence rate of the Chebyshev method with that of gradient descent, due to its more complex expression. However, by neglecting the denominator term $\xi^{-N}$, we obtain the following upper bound:

Note that the convergence rate of Polyak’s heavy-ball method matches (up to a multiplicative factor) that of Chebyshev’s method asymptotically, which is better than that of gradient descent in (2.7), which reads

We summarize this result in the following corollary, which compares the number of iterations required to reach a target accuracy $\epsilon$.

$N\geq\frac{L}{2\mu}~{}\log\left(\frac{\|x_{0}-x_{\star}\|_{2}}{\epsilon}\right)$ iterations of gradient descent (Algorithm 1 with (2.6)), or

$N\geq\sqrt{\frac{L}{2\mu}}~{}\log\left(\frac{\|x_{0}-x_{\star}\|_{2}}{\epsilon}\right)$ iterations of Chebyshev’s method (Algorithm 2),

For gradient descent, a sufficient condition on the number of iterations required to reach an accuracy of $\epsilon$ reads

Using the bound $\log\left(\frac{\frac{1}{x}+1}{\frac{1}{x}-1}\right)>2x$, the above condition can be simplified to the following stronger condition on $N$

This gives the desired result for gradient descent. With the same approach, we also get the result for the Chebyshev algorithm.

This corollary shows that the Chebyshev method can be $\sqrt{\kappa}$ faster than gradient descent. This translate to a speedup factor of $100$ in problems with a (reasonable) condition number of $10^{4}$, which is very significant.

When the dimension $d$ of the ambient space is sufficiently large, and without further assumptions on the spectrum of H, the worst-case guarantee on Chebyshev’s method is essentially unimprovable. Informally, given a budget $N$, a problem class $\mathcal{M}$ and some $R>0$, the best worst-case guarantee on the distance to optimality $\|x_{N}-x_{\star}\|_{2}$ that can be achieved by a first-order method is given by

which corresponds to the worst-case performance of the best performing method on any problem of the class. More precisely, for any first-order method satisfying $x_{k}\in x_{0}+\mbox{{span}}\left\{\nabla f(x_{0}),\,\nabla f(x_{1}),\,\ldots,\,\nabla f(x_{k-1})\right\}$ for all $k\geq 1$ (the “span assumption”) applied on the quadratic problem (2.1), it holds that:

and therefore, the conjugate gradient-like method

In other words, the $\max$ term in (2.22) searches for the “most difficult” quadratic function, while the $\min$ term represents the best first-order method for a specific quadratic function. Of course, the method (2.21) is much more powerful than the Chebyshev method, since it is optimal for any specific function. However, it is possible to show that despite being more powerful, this optimal algorithm has the same worst-case performance as that of Chebyshev’s method when the dimension $d$ of the ambient space is large enough. The lower bound result is summarized by the next theorem.

with $\xi=\frac{\sqrt{\frac{L}{\mu}}+1}{\sqrt{\frac{L}{\mu}}-1}$.

More details on this topic can be found in [nemirovskinotes1995, Section 12.3].

4 Notes and References

The Chebyshev method presented in this section is worst-case optimal for the class of quadratic functions with Hessians H satisfying $\mu\mathbf{I}\preceq\textbf{H}\preceq L\mathbf{I}$. More detailed discussions and developments on the topic of Chebyshev polynomials, for quadratic minimization, are provided in [Nemi83, nemirovsky1992information, Nest03a], as well as in the lecture notes [nemirovskinotes1995, Chapter 10]. Those references include the treatment of the case where the smallest eigenvalue is $\mu=0$. Finally, one should note that the optimal convergence bounds achieved by the Chebyshev method requires knowledge of the problem class parameters, $\mu$ and $L$, which might or might not be an issue, depending on the problem at hand.

Probably the most celebrated method for unconstrained quadratic optimization problems is the conjugate gradient (CG) method. Its origin is usually attributed to stiefel1952methods, straeter1971extension. As for the setup of this section, it turns out that CG methods are instance-optimal, in the sense that they are the best performing first-order methods on every particular problem instance in the range of unconstrained quadratic minimization problems (in particular, the CG variant presented in (2.21) achieves the lower bound from Theorem 2.3.5). The classical CG produces iterates $\{x_{k}\}_{k\geq 0}$ such that

which admits efficient formulation; see, e.g., [nocedal2006numerical]. Another variant of CG is often referred to as MINRES, which produces iterates $\{x_{k}\}_{k\geq 0}$ in the form

Its generalization GMRES [saad1986gmres] is popular for solving linear systems of the form $\textbf{H}x=b$ when H is not required to be either symmetric or invertible.

Yet another alternative for dealing with quadratic minimization is to resort on Anderson-type acceleration schemes. As for conjugate gradient methods, those schemes do not readily extend beyond quadratic minimization with the same nice theoretical guarantees. This is the topic of the next section.

Beyond quadratic optimization, properties of Chebyshev polynomials is the focus of [mason2002chebyshev]. The use of Chebyshev polynomials in the context of solving linear systems is covered at length in [fischerpolynomial]. In particular, the theory of [fischerpolynomial] can be instantiated for the convex quadratic minimization in average-case analyses, where Chebyshev polynomials (along with their heavy-ball limits) also naturally appear [pedregosa2020acceleration, lacotte2020optimal, scieur2020universal].

Chapter 3 Nonlinear Acceleration

In this section, we see that the main argument used in the Chebyshev method can be adapted beyond quadratic problems. The extension that we present here, called nonlinear acceleration, follows a pattern that is known in numerical analysis as vector extrapolation methods: it seeks to accelerate the convergence of sequences by extrapolation using nonlinear averages. Different such strategies are known under various names, starting with Aitken’s $\Delta^{2}$ [Aitk27], Wynn’s epsilon algorithm [Wynn56], and Anderson acceleration [anderson1965iterative]; a survey of these techniques can be found in [Sidi86]. The vector extrapolation techniques, generic by nature, can be applied to optimization, as explained in what follows.

This section focuses on the convex minimization problem:

We aim at adapting some of the ideas behind Chebyshev’s acceleration (see Section 2) to a broader class of convex minimization problems beyond quadratic minimization. These adaptations stems from a local quadratic approximation of the objective:

where $\textbf{H}=\nabla^{2}f(x_{\star})\in{\mbox{\bf S}}_{d}$ (the set of symmetric $d\times d$ matrices) is the Hessian of $f$ at $x_{\star}$, which we assume to satisfy $\mu\mathbf{I}\preceq\textbf{H}\preceq L\mathbf{I}$ for some $0<\mu<L$. Of course, neglecting the second-order term in (3.2) allows recovering a quadratic minimization problem for which one could apply Chebyshev’s method as is.

so that $\alpha_{k}$ does not depend on the particular problem instance $(\textbf{H},x_{\star})$ and on the initialization $x_{0}$, but only on the problem class described by $\mu$ and $L$ (we note that in Section 2, (3.3) was expressed in terms of optimizing a polynomial (2.12)). A natural alternative to Chebyshev’s method consists in choosing those weights adaptively. That is, depending on the particular instance of the problem at hand. For doing so, we have to choose another way to measure performance (because minimizing $\|x_{k}-x_{\star}\|_{2}$ would require knowledge of $x_{\star}$); one such possibility is to rely on function values or gradient norms. One could then rely on conjugate gradient-type methods which are very attractive for unconstrained quadratic minimization (see, e.g., discussions in Section 2.4). In this section, we consider the case where a first-order optimization method provided us with a sequence of pairs $\{(x_{i},\nabla f(x_{i})\}_{i=0,\ldots,k}$ satisfying (for $i=1,\ldots,k$)

and we study methods producing approximations of $x_{\star}$ as linear combinations $\sum_{i=0}^{k}c_{i}x_{i}$ of the previous iterates $\{x_{i}\}_{i=0,\ldots,k}$. For choosing the corresponding weights, we minimize the norm of the gradient at the approximated point. In unconstrained convex quadratic minimization problems, this approach is closely related to the so-called MINRES [paige1975solution] and GMRES [saad1986gmres] methods (conjugate gradient-type methods minimizing gradient norms; see discussions by walker2011anderson and Section 2.4). That is, when $f$ is quadratic, we choose the weights $\{c_{i}\}_{i=0,\ldots,k}$ by solving

Whereas the new approximation is a linear combination of previous iterates $\{x_{i}\}_{i=0,\ldots,k}$, the coefficients $\{c_{i}\}_{i=0,\ldots,k}$ depend nonlinearly on both $\nabla f$ and on $\{x_{i}\}_{i=0,\ldots,k}$. This technique is known under a few different names including those of Anderson acceleration and minimal polynomial extrapolation (see discussions and references in Section 3.6 for more details). In this section, we refer to all these methods as “nonlinear acceleration” techniques.

2 Nonlinear Acceleration for Quadratic Minimization

In this section, we present the simplest form of nonlinear acceleration, which is often referred to as the offline nonlinear acceleration mechanism. We start with the main arguments underlying the technique, and present a few variants later in this section.

The core idea of the mechanism is to use a sequence of iterates $\{x_{i}\}_{i=0,\ldots,k}$ provided by a first-order method for solving (3.1). On this basis, we generate a new approximation of a solution to (3.1) as a linear combination of past iterates, in the form $x_{\text{extr}}=\sum_{i=0}^{k}c_{i}x_{i}$. The point $x_{\text{extr}}$ is commonly referred to as the extrapolation and can be chosen in different ways. In classical nonlinear acceleration mechanisms, it is chosen for making $\|\nabla f(x_{\text{extr}})\|_{2}$ small, as in (3.5). In general, solving (3.5) is just as costly as solving (3.1), but the mechanism turns out to have an efficient formulation when minimizing quadratic functions of the form

In this case, it is possible to find an explicit formula for (3.5). Indeed, the gradient of $f$ is then a linear function, and because the coefficients $\{c_{i}\}_{i=0,\ldots,k}$ sum to one, the gradient of the linear combination is equal to a linear combination of gradients:

It follows that (3.5) reduces to a simple quadratic program involving gradients of the past iterates. It can be formulated as

For convenience, we use the following more compact form in the sequel

where $\textbf{G}=[\nabla f(x_{0}),\ldots,\nabla f(x_{k})]$ is the matrix formed by concatenating past gradients. This quadratic subproblem requires solving a small linear system of $(k+1)$ equations. When $\textbf{G}^{T}\textbf{G}$ is invertible, an explicit solution is provided by

This mechanism is summarized in Algorithm 3.

In this section, we quantify the accuracy of nonlinear acceleration. In particular, we show that it is instance-optimal and achieves the same worst-case convergence rate as that Chebyshev’s method (see Theorem 2.3.1 and Theorem 2.3.5) in the worst-case, as soon as the sequence $\{x_{i}\}_{i=0,\ldots,k}$ is generated by a reasonable first-order method. Before going into the analysis, we introduce a few technical ingredients specifying what is a reasonable first-order method. In short, we require that $\{x_{i}\}_{i=0,\ldots,k}$ is obtained from a “nondegenerate” first-order method. That is, we assume that the method uses $\nabla f(x_{i})$ non-trivially for generating $x_{i+1}$ (for all $i=0,\ldots,k-1$).

The proposition below shows that when the sequence $\{x_{i}\}_{i=0,1,\ldots}$ is generated by a nondegenerate first-order method, the gradient of any iterate $x_{i}$ can be written using a polynomial of degree exactly $i$.

We proceed by induction. First, we have $P_{0}=1$ ($P_{0}$ has degree and $P_{0}(0)=1$) and hence

thereby trivially reaching the desired conclusion for $i=0$.

We proceed with the induction hypothesis, assuming that the desired result holds for $x_{i}$. By Definition 3.2.1 (nondegenerate first-order method), and because $f$ is a quadratic function (3.6), we have

Thanks to the induction hypothesis, we have $\nabla f(x_{j})=P_{j}(\textbf{H})\nabla f(x_{0})$ with $P_{j}(0)=1$ and $\deg(P_{j})=j$ for $j\leq i$. For showing that $P_{i+1}$ satisfies the desired claim, we start by expressing $\nabla f(x_{i+1})$ in terms of a polynomial:

It is relatively straightforward to verify $P_{i+1}(0)=1$ using the previous expression:

Finally, a minor reorganization of the expression of $P_{i+1}$ allows writing

Nondegeneracy of the first-order method implies that there exists some $\alpha_{i}^{(i)}\neq 0$ such that the previous expression holds. Finally it follows from $\deg(P_{i})=i$ (induction hypothesis) that $\deg(\textbf{H}P_{i}(\textbf{H}))=i+1$, thereby reaching the desired claim.

Equipped with previous technical ingredients, one can show that Algorithm 3 is “instance-optimal” when applied to a nondegenerate first-order method. This means that the nonlinear acceleration algorithm finds the best polynomial given a specific quadratic function $f$—in opposition to the Chebyshev method that finds the best polynomial for a class of functions (Theorem 2.3.1). In other words, nonlinear acceleration adaptively looks for the best combination of previous iterates given the information stored in the previous gradients, while Chebyshev’s method uses the same worst-case optimal polynomial in all cases. Moreover, Algorithm 3 does not require knowledge of the smoothness or strong convexity parameters.

where $x_{\text{extr}}$ is obtained from Algorithm 3 $\{(x_{i},\nabla f(x_{i}))\}_{i=0,\ldots,k}$, and $\mathcal{P}_{k}$ is the set of polynomials of degree at most $k$.

Since the coefficients $c_{i}$ sum to one, we have the following equalities:

We now use the definition of a first-order method, which yields iterates $x_{i}$ such that

From Section 2.3, Proposition 2.2.1, it follows that $x_{i}$ can be written as

It also holds that its gradient can be written using the same polynomial:

By substituting this expression in the objective of (3.8), we obtain

Since the iterates $\{x_{i}\}_{i=0,\ldots,k}$ are generated by a nondegenerate first-order method, all polynomials $P_{i}$ have differents degrees. Therefore, the polynomials $\{P_{i}\}_{i=0,\ldots,k}$ are linearly independent, and hence $\{P_{i}\}_{i=0,\ldots,k}$ is a basis for the space $\mathcal{P}_{k}$. Finally, because $P_{i}(0)=1$ and $c^{T}\mathbf{1}=1$, we can rephrase the objective (3.13) as

This theorem shows that the worst-case convergence rate is essentially controlled by the optimal value of a minimization problem. The minimum in (3.11) can be bounded using a Chebyshev argument similar to the main argument used in Section 2.

where $x_{\text{extr}}$ is the output of Algorithm 3 applied to $\{(x_{i},\,\nabla f(x_{i}))\}_{i=0,\ldots,k}$.

We use the shifted Chebyshev polynomial from Theorem 2.3.1 as a feasible solution of the minimization problem (3.11).

As for the Chebyshev method, it is possible to show that the convergence rate cannot be improved as it matches that of the corresponding lower bound, which can be obtained by adapting Theorem 2.3.5 to gradient norms; see e.g., [nemirovskinotes1995, Proposition 12.3.2].

When $\{(x_{i},\,\nabla f(x_{i}))\}_{i=0,\ldots,k}$ is generated by a nondegenerate first-order method and when $k$ is large enough, nonlinear acceleration eventually converges exactly to the minimizer of the quadratic function $f$ (this follows easily from analogies with conjugate gradient-type methods; see, e.g.,Section 2.4). More formally, for all $k\geq d$ it holds that

where $x_{\text{extr}}$ is the output of Algorithm 3 applied to $\{(x_{i},\,\nabla f(x_{i}))\}_{i=0,\ldots,k}$. This is a natural consequence of the fact that $x_{\text{extr}}$ is the best point in $x_{0}+\mbox{{span}}\{\nabla f(x_{0}),\ldots,\nabla f(x_{k-1})\}$, as provided by (3.8).

2.2 Computational Complexity

The computational complexity of Algorithm 3 is $O(dk^{2}+k^{3})$, where $k$ is the length of the sequence of gradients and $d$ is the dimension of the ambient space. The first term originates from the matrix-matrix multiplication $\textbf{G}^{T}\textbf{G}$ in step 1, and the second one comes from solving a $(k+1)\times(k+1)$ matrix in step 2. The length $k$ being often much smaller than $d$, the resulting complexity is typically $O(dk^{2})$.

When using the nonlinear acceleration method in parallel with a first-order method generating a growing sequence of iterates $\{x_{i}\}_{i=0,\ldots,k}$, an extrapolation step can be computed each time a new iterate is produced. we can reduce the per iteration complexity up to $O(dk)$ by computing the matrix $\textbf{G}^{T}\textbf{G}$ and the coefficients $c_{\star}$ using low-rank updates [sidi1991efficient].

Because the iteration complexity of nonlinear acceleration grows with the length of the input sequence, it might become costly to compute an extrapolation. It is therefore common to use nonlinear acceleration only with the last few iterates produced by the first-order method. More formally, if we impose a maximum memory of $m$ pairs, we compute the extrapolation as follow when $k\geq m$. More formally, if we impose a maximum memory of $m$ pairs, one can use Algorithm 3 with input $\{(x_{k-m+i},\,\nabla f(x_{k-m+i}))\}_{i=1,\ldots,m}$.

2.3 Online Nonlinear Acceleration

So far, we have seen a post-processing procedure that generates an extrapolated point $x_{\text{extr}}$ from a sequence of pairs $\{x_{i},\nabla f(x_{i})\}_{i=0,\ldots,k}$. If this sequence is generated by a nondegenerate first-order method, then Corollary 3.2.6 shows that the gradient of the extrapolated point $x_{\text{extr}}$ converges to zero at an optimal worst-case convergence rate, without any hyper-parameters. Perhaps surprisingly, Algorithm 3 is itself not a nondegenerate first-order method, and can therefore not be used recursively as is.

In what follows, we introduce a mixing parameter. This parameter transforms the nonlinear acceleration mechanism to a nondegenerate first-order method without hurting its worst-case performance. This enables using nonlinear acceleration recursively for generating the whole sequence $\{x_{i}\}_{i=1,2,\ldots}$. This technique is often referred to as online nonlinear acceleration.

The idea underlying the mixing parameter is fairly simple: instead of combining previous iterates, we combine gradient steps as follows:

as provided by Algorithm 4. Furthermore, the use of an appropriate step size $h$ can even slightly improve the worst-case convergence speed of Algorithm 3.

Intuitively, this mixing between iterates and gradients emulates a gradient step on the extrapolated point from Algorithm 3, that we call here $x_{\text{extr}}^{\text{offline}}=\sum_{i=0}^{k}c_{i}x_{i}$,

where the second equality comes from the fact that $f$ is a quadratic function, and that $\sum_{i=0}^{k}c_{i}=1$.

This mixing parameter requires tuning one hyper-parameter $h$, which can be chosen in various ways. Proposition 3.2.9 shows that the mixing strategy slightly improves the performance of nonlinear acceleration if $h$ is set properly.

where $x_{\text{extr}}^{\text{offline}}$ is obtained from Algorithm 3 applied to $\{(x_{i},\,\nabla f(x_{i}))\}_{i=0,\ldots,k}$, $x_{\text{extr}}^{\text{mixing}}$ is the extrapolation with mixing from (3.15), and

Moreover, the factor $C_{h}$ is guaranteed to be smaller than one if $h\in(0,\,\frac{2}{L})$, and takes its minimal value at $h=\frac{2}{L+\mu}$.

Because $f$ is the quadratic function (3.6), the gradient of $x_{\text{extr}}^{\text{mixing}}$ reads

and the desired result follows from $\mu\mathbf{I}\preceq\textbf{H}\preceq L\mathbf{I}$.

As previously underlined, the mixing parameter transforms the nonlinear acceleration method into a nondegenerate first-order method. We can therefore use it recursively. The online variant of the nonlinear acceleration technique, with limited memory, is provided in Algorithm 5, using Algorithm 4 as a subroutine. One should note that when $m=\infty$ (no memory restriction), the worst-case performance of offline version of nonlinear acceleration with mixing (Algorithm 4) is also valid for its online variant (Algorithm 5). It follows from Proposition 3.2.9 that the worst-case performance of Algorithm 4 and Algorithm 5 is no worse than that of offline version of nonlinear acceleration (Algorithm 3), provided by Corollary 3.2.6.

In the next section, we see that nonlinear acceleration technique might suffer from serious instability issues when applied beyond quadratic minimization. Perhaps luckily, a simple regularization technique allows stabilizing the procedure beyond quadratics.

3 Regularized Nonlinear Acceleration Beyond Quadratics

where $\textbf{E}\triangleq[e_{0},\,\ldots,\,e_{k}]$ is the noise matrix. In this bound, the perturbation impacts the solution proportionally to its norm and to the conditioning of $\textbf{G}+\textbf{E}$.

Unfortunately, even for small perturbations, the condition number of $\textbf{G}+\textbf{E}$ and the norm of the vector $c$ are usually huge. In fact, G has a Krylov matrix structure, which is notoriously poorly conditioned [Tyrt94]. Thereby, even a very small perturbation E might have a significant impact on performance. For illustrating this, let us briefly illustrate the link between G and Krylov matrices: consider using gradient descent with step size $\gamma$ on a quadratic function; the iterates follow the rule

Thereby, a matrix G formed by these expressions has the form

which shows that G is in fact a Krylov matrix—by definition, a Krylov matrix $K$ associated with a matrix $A$ and vector $v$ is defined as $K=[v,\,Av,\,A^{2}v,\ldots]$.

In Figure 3.1, we show the norm of $c$ and the condition number of the matrix $\textbf{G}^{T}\textbf{G}$ when it is formed from iterates of gradient descent, accelerated gradient descent (see Section 4), and nonlinear acceleration (in the online setting, see Algorithm 5) for minimizing some randomly generated quadratic function. Figure 3.1(a) shows that even after 3 iterations, the system can already be considered singular (i.e., the condition number exceeds $10^{16}$).

For stabilizing the method, it is common to regularize the linear system. The resulting algorithm is often referred to as regularized nonlinear acceleration (RNA) [scieur2016regularized]. The following section is devoted to some theoretical properties of this method.

Regularized nonlinear acceleration (RNA) consists of using Algorithm 3 with a regularization, thereby rendering the method less sensitive to noise. In short, the base operation underlying RNA is to solve

As in the quadratic case, Algorithm 6 and Algorithm 7 could be used as subroutines in the online nonlinear acceleration method (Algorithm 5), thereby forming the regularized version of the online acceleration algorithm.

3.2 Perturbed Linear Gradients

In this section, we consider the problem of minimizing a twice continuously differentiable convex function $f$, as in (3.1), beyond quadratic problems. For doing that, we introduce perturbed linear gradients. As before, we consider iterates $\{x_{i}\}_{i=0,1,\ldots}$ originating from a first-order method satisfying

However, $\nabla f(x)$ is now no longer the gradient of a quadratic function. Instead, gradients of $f$ can be written as a sum of the gradients of a quadratic function with a perturbation term $e(x)$, as follows:

Indeed, it follows from twice continuous differentiability of $f$ that

Therefore, $f(x)$ can be approximated by the quadratic function (3.18) with $\textbf{H}=\nabla^{2}f(x_{\star})$. Similarly, its gradient reads

where $e(x)$ is the first-order Taylor remainder of the gradient. Thus, minimizing a non-quadratic function is equivalent to minimizing a perturbed quadratic one with a second-order error on its gradient:

where $\textbf{H}=\nabla^{2}f(x_{\star})$.

3.3 Convergence Bound

Using a perturbation argument, it is possible to derive a convergence guarantee for RNA. We state here a simplified version of [scieur2018online, Theorem 3.2], which describes how regularization balances acceleration and stability in Algorithm 7. We discuss convergence rates in greater detail in what follows.

where $x_{\text{extr}}$ is the output of Algorithm 7 applied to $\{(x_{i},\,\nabla f(x_{i}))\}_{i=0,\ldots,k}$ with parameters $h,\,\lambda>0$, and $V_{k}^{[\mu,L]}(\lambda)$ is a constant that corresponds to the maximum value on the interval $[\mu,L]$ of the regularized Chebyshev polynomial, i.e,

where $\|P\|_{2}$ is the norm of the vector of coefficients of the polynomial $P$.

This theorem states that regularization helps stabilizing the algorithm while slowing down the convergence rate. The regularized Chebyshev polynomial is somehow a mid-point between the classical shifted Chebyshev polynomial $C_{k}^{[\mu,L]}$ (from (2.15)) and the polynomial whose coefficients are defined by $\mathbf{1}/k$ (the polynomial that averages the iterates $\{x_{i}\}_{i=0,\ldots,k}$). By construction, its maximum value is always larger than that of the Chebyshev polynomial, but the norm of its coefficients is smaller. Unfortunately, there is as yet no known explicit expression of the regularized Chebyshev polynomial. To the best of our knowledge, its value can nevertheless be computed numerically [Barr20].

This mid-point between Chebyshev coefficients and the simple averaging of iterates is also natural in the context of noisy iterations. When the noise is negligible, a small regularization parameter combines the iterates $\{x_{i}\}_{i=0,\ldots,k}$ using nearly the classical Chebyshev weights. When the noise is more substantial, a larger regularization parameter brings the vector of coefficients $c$ closer to the average $\mathbf{1}/k$, thereby improving the “stability term” while rendering the “acceleration” less effective.

3.4 Asymptotic Convergence Rate

We briefly discuss the behavior of RNA when the initial point $x_{0}$ approaches the solution $x_{\star}$. In particular, the next proposition shows that if the perturbation magnitude $\epsilon$ decreases faster than $\|\textbf{H}(x_{0}-x_{\star})\|_{2}$, the parameter $\lambda$ can be adjusted to ensure an asymptotic convergence rate comparable to that of the Chebyshev method on quadratic problems (see Section 2).

Informally, the theorem exploits the fact that as $x_{0}$ approaches $x_{\star}$, $f$ gets closer to its quadratic approximation around $x_{\star}$. Thereby, an appropriate tuning of RNA allows matching (asymptotically) the convergence rate of nonlinear acceleration on quadratics (see Theorem 3.2.4).

and if we set $\lambda\propto\|x_{0}-x_{\star}\|_{2}^{s}$ (proportional to $\|x_{0}-x_{\star}\|_{2}^{s}$), where $0<s<2(\alpha-1)$, then it holds that

where $x_{\text{extr}}$ is the output of Algorithm 7 applied to the sequence $\{(x_{i},\nabla f(x_{i}))\}_{i=0,\ldots,k}$ with parameters $h,\,\lambda>0$.

To simplify the notation, set $R\triangleq\|\textbf{H}(x_{0}-x_{\star})\|_{2}$. We start from the result of Theorem 3.3.1 and divide both sides by $R$:

Since $\lambda\propto R^{s}$ and $\epsilon=O(R^{\alpha})$,

When $x_{0}\rightarrow x_{\star}$, we have $R\rightarrow 0$ and

Finally, in (3.20) the regularization parameter $\lambda\|\textbf{G}^{T}\textbf{G}\|_{2}=O(R^{2})\rightarrow 0$. Since the (non-regularized) shifted Chebyshev polynomial $C_{k}^{[\mu,L]}=C_{k}^{[\mu,L],0}$ is a feasible solution of (3.20), we have the following bounds:

As $\lambda\rightarrow 0$ we have that the upper bound on $(V_{k}^{[\mu,L]}(\lambda)$ converges to the maximum value of the regular (shifted) Chebyshev polynomial, thereby reaching the desired claim.

In short, the previous theorem states that the asymptotic convergence rate matches the rate of Chebyshev’s method as soon as $\lambda\propto\|x_{0}-x_{\star}\|_{2}^{s}$ is decreasing (condition $s>0$), but not too quickly compared to the perturbation magnitude $\epsilon=O(\|x_{0}-x_{\star}\|_{2}^{\alpha})$ (condition $s<2(\alpha-1)$), which is achievable only when $\alpha>1$. This condition is met, for instance, when accelerating twice continuously differentiable functions with gradient descent: the error decreases as $O(\|x_{0}-x_{\star}\|_{2}^{2})$, see (3.19), and therefore $\alpha=2>1$.

4 Extensions

The previous sections presented the nonlinear acceleration mechanism for unconstrained convex quadratic minimization. It also contained an analysis of its regularized version when applied beyond quadratics. In this section, we briefly cover two natural extensions: (i) the application of nonlinear acceleration to iterates that are corrupted by a stochastic noise, and (ii) the application of nonlinear acceleration to constrained/composite convex optimization problems when a projection/proximal operator is used.

to $e_{i}$ being the sum of a stochastic noise with a Taylor remainder. This is typically the case when applying RNA to stochastic gradient descent (SGD) and related methods. In this case, Theorem 3.3.1 holds in expectation under standard assumptions [Scie17], such as a bounded variance of $e_{i}$. However, the asymptotic convergence result from Proposition 3.3.2 may not be achieved. Indeed, in this setting, Proposition 3.3.2 also holds in expectation under the condition

Unfortunately, an asymptotic acceleration is not always possible. For instance, when trying to accelerate the fixed step SGD, we have $\epsilon=O(1)$ (i.e., $\alpha=0$) and the asymptotic guarantee does not apply. This is probably not a surprise as this SGD does not converge to the optimum, hence there is no apparent reason for any sequence extrapolation technique to work at all. Fortunately, Algorithm 6 does usually work for “variance reduced” first-order methods [Scie17], such as SAG [schmidt2017minimizing], SAGA [defazio2014saga], or SVRG [johnson2013accelerating].

It is common to apply first-order methods to composite convex minimization problems of the form:

where $f$ is a smooth strongly convex function (this class of functions is used intensively in Section 4; see Definition 4.1.1) and $h$ is a closed, proper, and convex function (i.e., $h$ has a closed, non-empty, and convex epigraph) and whose proximal operator is available:

for some step size $\gamma>0$. Problem (3.21) can then be approached iteratively via the proximal gradient method:

We omit most of the details on proximal algorithms; see Section 4 and Section 5 for more details and references. For instance, when $h$ is the indicator function of a non-empty closed convex set $\mathcal{C}$, the proximal operator corresponds to an orthogonal projection onto $\mathcal{C}$ and the proximal gradient method reduces to the projected gradient method.

Unfortunately, a naive use of nonlinear acceleration on the iterates $\{x_{i}\}_{i=0,\ldots,k}$ does not immediately work in this context, for several reasons. In particular, it is not possible to ensure that the extrapolated point $x_{\text{extr}}$ belongs to $\mathop{\bf dom}(h)$ (or to the set $\mathcal{C}$ when $h$ is an indicator function for $\mathcal{C}$). Moreover, due to the use of the proximal operator, the iterates $\{x_{i}\}_{i=0,\ldots,k}$ do not necessarily satisfy the span assumption (3.4).

Recently, mai2019anderson adapted the Anderson Acceleration method to handle a large class of constrained and non-smooth composite problems. The main idea is as follows: instead of accelerating the sequence $\{x_{i}\}_{i=0,\ldots,k}$ generated by the proximal gradient method (3.23), we accelerate an alternate sequence $\{z_{i}\}_{i=0,\ldots,k}$ which satisfies

This sequence $\{z_{i}\}_{i=0,\ldots,k}$ corresponds to the sequence generated by (3.23) with the ordering of the gradient and proximal steps being swapped.

This trick allows obtaining convergence bounds for nonlinear acceleration in the proximal setup under very few changes in the algorithm. In particular, mai2019anderson show that using Algorithm 8 in the presence of a proximal operator does not change the convergence analysis—using Clarke’s generalized Jacobian [clarke1990optimization], semi-smoothness [mifflin1977semismooth, qi1993nonsmooth] and assuming that the function $h$ is twice epi-differentiable and that $h$ is twice-differentiable around the solution $x_{\star}$. We refer the reader to [rockafellar2009variational, Section 13] for a comprehensive treatment of epi-differentiability.

5 Globalization Strategies and Speeding-up Heuristics

As for many standard optimization methods, such as quasi-Newton methods, RNA only has local convergence guarantees beyond quadratics. Therefore, it is common to embed the mechanism with some globalization strategies, a.k.a. safeguards. Those strategies ensure not to deteriorate too much the performance of the initial first-order method in cases where RNA is used beyond its guaranteed range of applications. Those strategies can also be seen as speeding-up heuristics.

It is in general not guaranteed that the extrapolated point $x_{\text{extr}}$ is better than any iterate of the sequence $\{x_{i}\}_{i=0,\ldots,k}$ produced by the original method. This situation might for example occur when extrapolating with a bad mixing or regularization parameter, or simply when the error terms are too large. One classical way of limiting the impact of such problems is by checking some descent condition. For instance, one might consider “accepting” $x_{\text{extr}}$ only if it is better than previous iterates $\{x_{i}\}_{i=0,\ldots,k}$:

Nonlinear acceleration requires the selection of a mixing parameter, which might be difficult to tune in practice. One common trick is to choose it via a line-search strategy. That is, defining:

one can choose $h$ by approximately solving $\mathop{\rm argmin}_{h}f\left(x_{\text{extr}}(h)\right)$.

6 Notes and References

Nonlinear acceleration techniques have been studied extensively during recent decades, and excellent reviews can be found in [smith1987extrapolation, jbilou1991some, brezinski1991extrapolation, jbilou1995analysis, jbilou2000vector, brezinski2001convergence, brezinski2019genesis]. The first usage of an acceleration technique for fixed point iteration can be traced back to [gekeler1972solution, brezinski1971algorithme, brezinski1970application].

There are numerous independent works leading to methods similar to those described here. The most classical, and probably the most similar, is Anderson acceleration [anderson1965iterative], which corresponds exactly to the online mode of nonlinear acceleration (without regularization). Despite it being an old algorithm, there has been a recent uptake of interest in the convergence analysis [walker2011anderson, toth2015convergence] of Anderson acceleration thanks to its good empirical performance, and strong connection with quasi-Newton methods [fang2009two].

Other versions of nonlinear acceleration use different arguments but behave similarly. For instance, minimal polynomial extrapolation (MPE), which uses the properties of the minimal polynomial of a matrix [cabay1976polynomial]; reduced rank extrapolation (RRE); and the Mesina method [mevsina1977convergence, eddy1979extrapolating] are also variants of Anderson acceleration. The properties and equivalences of these approaches whave been studied extensively during the past decades [sidi1988extrapolation, ford1988recursive, sidi1991efficient, jbilou1991some, sidi1998upper, sidi2008vector, sidi2017minimal, sidi2017vector, brezinski2018shanks, brezinski2020shanks]. Unfortunately, these methods do not extend well to nonlinear functions, especially due to conditioning problems [sidi1986convergence, sidi1988convergence, scieur2016regularized]. Recent works have nevertheless proven the convergence of such methods, provided that good conditioning of the linear system [sidi2019convergence] can be ensured.

There are also other classes of nonlinear acceleration algorithms, based on existing algorithms, for accelerating the convergence of scalar sequences [brezinski1975generalisations]. For instance, the topological epsilon vector algorithm (TEA) extends the idea of the scalar $\varepsilon$-algorithm of [Wynn56] to vectors.

Chapter 4 Nesterov Acceleration

This section presents a systematic interpretation of the acceleration of the gradient method stemming from Nesterov’s original work [Nest83]. The early parts of the section are devoted to the gradient method and the “optimized gradient method,” due to Dror14 and kim2016optimized. The motivations and ideas underlying the latter are intuitive and very similar to those behind the introduction of Chebyshev methods for optimizing quadratic functions (see Section 2). Furthermore, the optimized gradient method has a relatively simple format and proof and can be used as an inspiration for developing numerous variants with wider ranges of applications, including Nesterov’s early accelerated gradient methods [Nest83, Nest13] and the fast iterative shrinkage-thresholding algorithm [Beck09, FISTA]. Although some parts of this section are more technical, we believe all the ideas can be reasonably well understood even when skipping, or skimming through the algebraic proofs. The section and the proofs are organized so that each time an additional ingredient (strong convexity, constraints, etc.) is included, its inclusion only requires a few additional ingredients compared to the previous (simpler) proofs of the base versions of the method.

We start with the theory and interpretation of acceleration in a simple setting: smooth unconstrained convex minimization in a Euclidean space. All subsequent developments follow from the same template, namely a linear combination of regularity inequalities, with additional ingredients being added one by one. The next part is devoted to methods that take advantage of strong convexity by using the same ideas and algorithmic structures. On the way, we provide a few different (equivalent) templates for the algorithms, since in more advanced settings, those templates do not generalize in the same way. We then recap and discuss a few practical extensions for handling constrained problems, nonsmooth regularization terms, unknown problem parameters/line-searches, and non-Euclidean geometries. Finally, we briefly discuss a popular ordinary differential equation (ODE)-based interpretation of Nesterov’s method. Techniques for obtaining the worst-case analyses presented throughout this text are presented in Appendix C, and notebooks for simpler reproduction of the proofs are provided in Section 4.9.

In the first part of this section, we consider smooth unconstrained convex minimization problems. This type of problems is a direct extension of unconstrained convex quadratic minimization problems where the quadratic function has eigenvalues bounded above by some constant. More precisely, we consider the simple unconstrained differentiable convex minimization problem

where $f$ is convex with an $L$-Lipschitz gradient (we call such functions convex and $L$-smooth, see Definition 4.1.1 below), and we assume throughout that there exists a minimizer $x_{\star}$. The goal of the methods presented below is to find a candidate solution $x$ satisfying $f(x)-f_{\star}\leq\epsilon$ for some $\epsilon>0$. Depending on the target application, other quality measures, such as guarantees on $\|\nabla f(x)\|_{2}$ or $\|x-x_{\star}\|_{2}$, might be preferred. We refer to Section 4.9 “Changing the performance measure”, for discussions on this topic.

We start with the analysis of gradient descent and then show that its iteration complexity can be significantly improved using an acceleration technique proposed by Nest83.

After presenting the theory for the smooth convex case we see how it goes in the smooth strongly convex one. This class of problems extends to that of unconstrained convex quadratic minimization problems where the quadratic function has eigenvalues respectively bounded above and below by some constants $L$ and $\mu$.

Furthermore, we denote by $q=\tfrac{\mu}{L}$ the inverse condition number (that is, $q=\tfrac{1}{\kappa}$, with $\kappa$ is the usual condition number as used e.g., in Section 2) of functions in the class $\mathcal{F}_{\mu,L}$.

Figure 4.1 provides an illustration of the global quadratic upper approximation (with curvature $L$) on $f(\cdot)$ due to smoothness and of the global quadratic lower approximation (with curvature $\mu$) on $f(\cdot)$ due to strong convexity.

A number of inequalities can be written to characterize functions in $\mathcal{F}_{\mu,L}$: see, for example, Nest03a. When analyzing methods for minimizing functions in this class, it is crucial to have the right inequalities at our disposal, as worst-case analyses essentially boil down to appropriately combining such inequalities. We provide the most important inequalities along with their interpretations and proofs in Appendix A. In this section, we only use three. First, we use the quadratic upper and lower bounds arising from the definition of smooth strongly convex functions, that is, (4.2) and (4.3). For some analyses however, we need an additional inequality, provided by the following theorem. This inequality is often referred to as an interpolation (or extension) inequality. Its proof is relatively simple: it only consists of requiring all quadratic lower bounds from (4.3) to be below all quadratic upper bounds from (4.2) (details in Appendix A.1). It can be shown that worst-case analyses of all first-order methods for minimizing smooth strongly convex functions can be performed using only this inequality for some specific values of $x$ and $y$ (details in Appendix C).

As discussed later in Section 4.3.3, this inequality has some flaws. Therefore, we only use (4.2) and (4.3) whenever possible.

Before continuing to the next section, we mention that both smoothness and strong convexity are strong assumptions. More generic assumptions are discussed in Section 6 to obtain improved rates under weaker assumptions.

2 Gradient Method and Potential Functions

In this section, we analyze gradient descent using the concept of potential functions. The resulting proofs are technically simple, although they might not seem to provide any direct intuition about the method at hand. We use the same ideas to analyze a few improvements on gradient descent before providing interpretations underlying this mechanism.

The simplest and probably most natural method for minimizing differentiable functions is gradient descent. It is often attributed to cauchy1847methode and consists of iterating

where $\gamma_{k}$ is some step size. There are many different techniques for picking $\gamma_{k}$, the simplest of which is to set $\gamma_{k}=1/L$, assuming $L$ is known—otherwise, line-search techniques are typically used; see Section 4.7. Our present objective is to bound the number of iterations required by gradient descent to obtain an approximate minimizer $x_{k}$ of $f$ that satisfies $f(x_{k})-f_{\star}\leq\epsilon$.

2.2 A Simple Proof Mechanism: Potential Functions

Potential (or Lyapunov/energy) functions are classical tools for proving convergence rates in the first-order literature, and a nice recent review of this topic is given by bansal2019potential. For gradient descent, the idea consists in recursively using a simple inequality (proof below),

as a potential and use $\phi_{k+1}\leq\phi_{k}$ as the building block for the worst-case analysis. Once such a potential inequality $\phi_{k+1}\leq\phi_{k}$ is established, a worst-case guarantee can easily be deduced through a recursive argument, yielding

and hence, $f(x_{N})-f_{\star}\leq\tfrac{L}{2N}\|x_{0}-x_{\star}\|^{2}_{2}$. We also conclude that the worst-case accuracy of gradient descent is $O(N^{-1})$ or equivalently, that its iteration complexity is $O(\epsilon^{-1})$. Therefore, the main inequality to be proved for this worst-case analysis to work is the potential inequality $\phi_{k+1}\leq\phi_{k}$. In other words, the analysis of $N$ iterations of gradient descent is reduced to the analysis of a single iteration, using an appropriate potential. This kind of approach was already used for example by Nest83, and many different variants of the potential function can be used to prove convergence of gradient descent and related methods in similar ways.

with $x_{k+1}=x_{k}-\tfrac{1}{L}\nabla f(x_{k})$ and $A_{k+1}=1+A_{k}$.

The proof consists of performing a weighted sum of the following inequalities:

convexity of $f$ between $x_{k}$ and $x_{\star}$, with weight $\lambda_{1}=A_{k+1}-A_{k}$:

smoothness of $f$ between $x_{k}$ and $x_{k+1}$ with weight $\lambda_{2}=A_{k+1}$:

The last inequality is often referred to as the descent lemma since substituting $x_{k+1}$ allows to obtain $f(x_{k+1})\leq f(x_{k})-\frac{1}{2L}\|\nabla f(x_{k})\|^{2}_{2}$.

The weighted sum forms a valid inequality:

Using $x_{k+1}=x_{k}-\tfrac{1}{L}\nabla f(x_{k})$, this inequality can be rewritten (by completing the squares or simply extending both expressions and verifying that they match on a term-by-term basis) as follows:

which can be reorganized and simplified to

where the last inequality follows from picking $A_{k+1}=A_{k}+1$ and neglecting the last residual term $-\tfrac{A_{k}}{2L}\|\nabla f(x_{k})\|^{2}_{2}$ (which is nonpositive) on the right-hand side.

A convergence rate for gradient descent can be obtained directly as a consequence of Theorem 4.2.1, following the reasoning of (4.5), and the worst-case guarantee corresponds to $f(x_{N})-f_{\star}=O(A_{N}^{-1})=O(N^{-1})$. We detail this in the next corollary.

Following the reasoning of (4.5), we recursively use Theorem 4.2.1, starting with $A_{0}=0$. That is, we define

and recursively use the inequality $\phi_{k+1}\leq\phi_{k}$ from Theorem 4.2.1, with $A_{k+1}=A_{k}+1$ and $A_{0}=0$; hence, $A_{k}=k$. We thus obtain

2.3 How Conservative is this Worst-case Guarantee?

Before moving to other methods, we show that the worst-case rate $O(N^{-1})$ of gradient descent is attained on very simple problems, motivating the search for alternate methods with better guarantees. This rate is observed on, e.g., all functions that are nearly linear over large regions. One such common function is the Huber loss (with $x_{\star}=0$, arbitrarily):

with $a_{\tau}=L\tau$ and $b_{\tau}=-\frac{L}{2}\tau^{2}$ to ensure its continuity and differentiability. On this function, as long as the iterates of gradient descent satisfy $|x_{k}|\geq\tau$, they behave as if the function were linear, and the gradient is constant. It is therefore relatively easy to explicitly compute all iterates. In particular, by picking $\tau=\frac{|x_{0}|}{2N+1}$, we get $f(x_{N})-f_{\star}=\frac{L\|x_{0}-x_{\star}\|^{2}_{2}}{2(2N+1)}$ and reach the $O(N^{-1})$ worst-case bound; see Dror14. Therefore, it appears that the worst-case bound from Corollary 4.2.3 for gradient descent can only be improved in terms of the constants, but the rate itself is the best possible one for this simple method; see, for example, [Dror14, drori2014contributions] for the corresponding tight expressions.

In the next section, we show that similar reasoning based on potential functions produces methods with improved worst-case convergence rate $O(N^{-2})$, compared to the $O(N^{-1})$ of vanilla gradient descent.

3 Optimized Gradient Method

Given that the complexity bound for gradient descent cannot be improved, it is reasonable to look for alternate, hopefully better, methods. In this section, we show that accelerated methods can be designed by optimizing their worst-case performance. To do so, we start with a reasonably broad family of candidate first-order methods described by

Of course, methods in this form are impractical since they require keeping track of all previous gradients. Neglecting this potential problem for now, one possibility for choosing the step size $\{h_{i,j}\}$ is to solve a minimax problem:

In other words, we are looking for the best possible worst-case ratio among methods of the form (4.6). Of course, different target notions of accuracy could be considered instead of $\left.(f(y_{N})-f_{\star})\middle/\|y_{0}-x_{\star}\|^{2}_{2}\right.$, but we proceed with this notion for now.

It turns out that (4.7) has a clean solution, obtained by kim2016optimized, based on clever reformulations and relaxations of (4.7) developed by Dror14 (some details are provided in Section 4.9). Furthermore, this method has “factorized” forms that do not require keeping track of previous gradients. The optimized gradient method (OGM) is parameterized by a sequence $\{\theta_{k,N}\}_{k}$ that is constructed recursively starting from $\theta_{-1,N}=0$ (or equivalently $\theta_{0,N}=1$), using

We also mention that optimized gradient methods can be stated in various equivalent formats, we provide two variants in Algorithm 9 and Algorithm 10 (a rigorous equivalence statement is provided in Appendix B.1.1). While the shape of Algorithm 10 is more common in accelerated methods, the equivalent formulation provided in Algorithm 9 allows for slightly more direct proofs.

Direct approaches to (4.7) are rather technical—see details in [Dror14, kim2016optimized]. However, showing that the OGM is indeed optimal on the class of smooth convex functions can be accomplished indirectly by providing an upper bound on its worst-case complexity guarantees and by showing that no first-order method can have a better worst-case guarantee on this class of problems. We detail a fully explicit worst-case guarantee for OGM in the next section. It consists in showing that

is a potential function for the optimized gradient method when $k<N$ (Theorem 4.3.1, below). For $k=N$, we need a minor adjustment (Lemma 4.3.3, below) to obtain a bound on $f(y_{N})-f_{\star}$ and not in terms of $f(y_{N})-f_{\star}-\tfrac{1}{2L}\lVert\nabla f(y_{N})\rVert^{2}_{2}$, which appears in the potential.

As in the case of gradient descent, the proof relies on potential functions. Following the recursive argument from (4.5), the convergence guarantee is driven by the convergence speed of $\theta_{k,N}^{-2}$ towards . We note that when $k<N-1$,

and therefore, $\theta_{k,N}\geq\tfrac{k}{2}+1$. We also directly obtain

and hence, $\theta_{N,N}^{-2}=O(N^{-2})$. Before providing the proof, we mention that it heavily relies on inequality (4.4) with $\mu=0$. This inequality is key for formulating (4.7) in a tractable way.

The main point now is to prove that (4.9) is indeed a potential for the optimized gradient method. We emphasize again that our main motivation for proving this is to show that the OGM provides a good template algorithm for acceleration (i.e., a method involving two or three sequences) and that the corresponding potential functions can also be used as a template for the analysis of more advanced methods.

Note that the potential structure does not seem immediately intuitive: it was actually found using computer-assisted proof design techniques; see Section 4.9 “On obtaining the proofs in this section” and Appendix C for further references. In particular, the following theorem can be found in taylor19bach.

when $y_{k}$ and $z_{k+1}$ are obtained from Algorithm 9.

Recall that the algorithm can be written as

The proof consists of performing a weighted sum of the following inequalities.

Smoothness and convexity of $f$ between $y_{k-1}$ and $y_{k}$ with weight $\lambda_{1}={2\theta_{k-1,N}^{2}}$:

Smoothness and convexity of $f$ between $x_{\star}$ and $y_{k}$ with weight $\lambda_{2}={2\theta_{k,N}}$:

Since the weights are nonnegative, the weighted sum produces a valid inequality:

which (either by completing the squares or simply by extending both expressions and verifying that they match on a term-by-term basis) can be reformulated as

The desired conclusion follows from picking $\theta_{k,N}\geq\theta_{k-1,N}$ satisfying

and hence the choice (4.8), thus reaching

A final technical fix is required now. To show that the optimized gradient method is an optimal solution to (4.7), we need an upper bound on the function values, rather than on the function values minus a squared gradient norm. This discrepancy is handled by the following technical lemma.

where $y_{N}$ is obtained from Algorithm 9.

The proof consists of performing a weighted sum of the following inequalities.

Smoothness and convexity of $f$ between $y_{N-1}$ and $y_{N}$ with weight $\lambda_{1}={2\theta_{N-1,N}^{2}}$:

Smoothness and convexity of $f$ between $x_{\star}$ and $y_{N}$ with weight $\lambda_{2}=\theta_{N,N}$:

Since the weights are nonnegative, the weighted sum produces a valid inequality:

The conclusion follows from choosing $\theta_{N,N}\geq\theta_{N-1,N}$ such that

By combining Theorem 4.3.1 and the technical Lemma 4.3.3, we get the final worst-case performance bound of the OGM on function values, detailed in the corollary below.

using Theorem 4.3.1 and technical Lemma 4.3.3. We obtain the last bound by using $\theta_{N,N}\geq(N+1)/\sqrt{2}$; see (4.11).

In the following section, we mostly use potential functions, relying directly on the function value $f(x_{k})$ instead of $f(y_{k})$ for practical reasons discussed below. Note that using the descent lemma (i.e., the inequality $f(x_{k+1})\leq f(y_{k})-\tfrac{1}{2L}\|\nabla f(y_{k})\|^{2}_{2}$) directly on the potential function allows us to obtain a bound on $f(x_{k+1})$ for the OGM. This result can be found in [kim2017convergence, Theorem 3.1] without the “potential function” mechanism.

3.2 Optimality of Optimized Gradient Methods

A nice, commonly used guide for designing optimal methods consists of constructing problems that are difficult for all methods within a certain class. This strategy results in lower complexity bounds, and it is often deployed via the concept of minimax risk (of a class of problems and a class of methods)—see, e.g., guzman2015lower—which corresponds to the worst-case performance of the best method within the prescribed class. In this section, we briefly discuss such results in the context of smooth convex minimization, on the particular class of black-box first-order methods. The term black-box is used to emphasize that the method has no prior knowledge of $f$ (beyond the class of functions to which $f$ belongs, so methods are allowed to use $L$) and that it can only obtain information about $f$ by a evaluating its gradient/function value through an oracle.

Of particular interest to us, drori2017exact established that the worst-case performance achieved by the optimized gradient method (see Corollary 4.3.5) on the class of smooth convex functions cannot in general be improved by any black-box first-order method.

In the previous sections, we showed that $\theta_{N,N}^{2}\geq\tfrac{(N+1)^{2}}{2}$. It is also relatively easy to establish that

thereby obtaining $\theta_{k,N}\leq k+1$ (because $\theta_{0,N}=1$) as well as

We conclude (through Theorem 4.3.7) that the lower bound has the form

While Drori’s approach to obtaining this lower bound is rather technical (a slightly simplified and weaker version of this result can be found in [drori2021exact, Corollary 5]), there are simpler approaches that allow us to show that the rate $\Omega(N^{-2})$ (that is, neglecting the tight constants) cannot in general be beaten in black-box smooth convex minimization. For one such example, we refer to [Nest03a, Theorem 2.1.6]. In a closely related line of work, [nemirovsky1991optimality] established similar exact bounds in the context of solving linear systems of equations and for minimizing convex quadratic functions (see also Section 2.3.4). For convex quadratic problems whose Hessian has bounded eigenvalues between and $L$, these lower bounds are attained by the Chebyshev (see Section 2) and by conjugate gradient methods [nemirovsky1991optimality, nemirovsky1992information].

Perhaps surprisingly, the conjugate gradient method also achieves the lower complexity bound of smooth convex minimization provided by Theorem 4.3.7. Furthermore, the proof follows essentially the same structure as that for the OGM. In particular, it relies on the same potential function (see Appendix B.2).

3.3 Optimized Gradient Method: Summary

Before going further, we quickly summarize what we have learned from the optimized gradient method. First of all, the optimized gradient method can be seen as a counterpart of the Chebyshev method for minimizing quadratics, applied to smooth convex minimization. It is an optimal method in the sense that it has the smallest possible worst-case ratio $\tfrac{f(y_{N})-f_{\star}}{\|y_{0}-x_{\star}\|^{2}_{2}}$ over the class $f\in\mathcal{F}_{0,L}$ among all black-box first-order methods, given a fixed computational budget of $N$ gradient evaluations. Furthermore, although this method has a few drawbacks (we mention a few below), it can be seen as a template for designing other accelerated methods using the same algorithmic and proof structures. We extensively use variants of this template below. In other words, most variants of accelerated gradient methods rely on the same two (or three) sequence structures, and on similar potential functions. Such variants usually rely on slight variations in the choice of the parameters used throughout the iterative process, typically involving less aggressive step size strategies (i.e., smaller values for $\{h_{i,j}\}$ in (4.6)).

Second, the OGM is not a very practical method as such: it is fined-tuned for unconstrained smooth convex minimization and does not readily extend to other situations, such as situations involving constraints, for which (4.4) does not hold in general; see the discussions in [drori2018properties] and Remark A.1.6.

On the other hand, we see in what follows that it is relatively easy to design other methods that follow the same template and achieve the same $O(N^{-2})$ rate, while resolving the issues of the OGM listed above. Such methods use slightly less aggressive step size strategies, at the cost of being slightly suboptimal for (4.7), i.e., they have slightly worse worst-case guarantees. In this vein, we start by discussing the original accelerated gradient method, proposed by Nest83.

4 Nesterov’s Acceleration

Motivated by the format of the optimized gradient method, we detail a potential-based proof for Nesterov’s method. We then quickly review the concept of estimate sequences and show that they provide an interpretation of potential functions as increasingly good models of the function to be minimized. Finally, we extend these results to strongly convex minimization.

In this section, we follow the algorithmic template provided by the optimized gradient method. In this spirit, we start by discussing the first accelerated method in its simplest form (Algorithm 11) as well as its potential function, originally proposed by Nest83, but the presentation here is different.

Our goal is to derive the simplest algebraic proof for this scheme. We follow the algorithmic template of the optimized gradient method (which is further motivated in Section 4.6.1). Once a potential is chosen, the proofs are quite straightforward as simple combinations of inequalities and basic algebra. Our choice of potential function is not immediately obvious but allows for simple extensions afterwards. Other choices are possible, for example, incorporating $f(y_{k})$ as (in the OGM) or additional terms such as $\|\nabla f(x_{k})\|^{2}_{2}$. We pick a potential function similar to that used for gradient descent and that of the optimized gradient method, which is written

where one iteration of the algorithm has the following form, reminiscent of the OGM:

Our goal is to select algorithmic parameters $\{(\tau_{k},\alpha_{k},\gamma_{k})\}_{k}$ so as to greedily make $A_{k+1}$ as large as possible as a function of $A_{k}$ since the convergence rate of the method is controlled by the inverse of the growth rate of $A_{k}$, i.e., $f(x_{N})-f_{\star}=O(A_{N}^{-1})$.

In practice, we can pick $A_{k+1}=A_{k}+\frac{1}{2}(1+\sqrt{4A_{k}+1})$ by choosing $\tau_{k}=1-A_{k}/A_{k+1}$, $\alpha_{k}=\frac{1}{L}$, and $\gamma_{k}=(A_{k+1}-A_{k})/{L}$ (see Algorithm 11), and the proof is then quite compact.

Before continuing to the proof of the potential inequality, we show that $A_{k}^{-1}=O(k^{-2})$. Indeed, we have,

where the last inequality follows from a recursive application of the previous one, along with $A_{0}=0$.

with $A_{k+1}=A_{k}+\tfrac{1+\sqrt{4A_{k}+1}}{2}$.

The proof consists of a weighted sum of the following inequalities.

Convexity of $f$ between $x_{\star}$ and $y_{k}$ with weight $\lambda_{1}=A_{k+1}-A_{k}$:

Convexity of $f$ between $x_{k}$ and $y_{k}$ with weight $\lambda_{2}=A_{k}$:

Smoothness of $f$ between $y_{k}$ and $x_{k+1}$ (a.k.a., descent lemma) with weight $\lambda_{3}=A_{k+1}$:

We therefore arrive at the following valid inequality

For the sake of simplicity, we do not substitute $A_{k+1}$ by its expression until the last stage of the reformulation. Substituting $y_{k}$, $x_{k+1}$, and $z_{k+1}$ by their expressions in (4.13) along with $\tau_{k}=1-A_{k}/A_{k+1}$, $\alpha_{k}=\frac{1}{L}$, and $\gamma_{k}=\frac{A_{k+1}-A_{k}}{L}$, basic algebra shows that the previous inequality can be reorganized as

The claim follows from selecting $A_{k+1}\geq A_{k}$ such that $A_{k+1}-(A_{k}-A_{k+1})^{2}=0$, thereby reaching

The final worst-case guarantee is obtained by using the same chaining argument as in (4.5), combined with an upper bound on $A_{N}$.

Following the argument of (4.5), we recursively use Theorem 4.4.1 with $A_{0}=0$:

where we used $A_{N}\geq N^{2}/4$ from (4.14) to reach the last inequality.

Before moving on, we emphasize that the rate of $O(N^{-2})$ matches that of lower bounds (see, e.g., Theorem 4.3.7) up to absolute constants.

Finally, note that Nesterov’s method is often written in a slightly different format, similar to that of Algorithm 10. The alternate formulation omits the third sequence $z_{k}$ and is provided in Algorithm 12. It is preferred in many references on the topic due to its simplicity. A third equivalent variant is provided in Algorithm 13; this variant turns out to be useful when generalizing the method beyond Euclidean spaces. The equivalence statements between Algorithm 11, Algorithm 12, and Algorithm 13 are relatively simple and are provided in Appendix B.1.2. Many references tend to favor one of these formulations, and we want to point out that they are equivalent in the base problem setup of unconstrained smooth convex minimization. Although the expression of the different formats in terms of the same external sequence $\{A_{k}\}_{k}$ does not always correspond to their simplest forms (i.e., alternate parameterizations might be simpler, particularly in the strongly convex case which follows), we proceed with this sequence to avoid introducing too many variations on the same theme.

4.2 Estimate Sequence Interpretation

We now relate the potential function approach to estimate sequences. That is, we relate acceleration to first-order methods maintaining a model of the function throughout the iterative procedure. This approach was originally developed in [Nest03a, Section 2.2], and it has since been used in numerous works to obtain accelerated first-order methods in various settings (see discussions in Section 4.9). We present a slightly modified version, related to those of [baes2009estimate, wilson2016lyapunov], which simplifies our comparisons with the previous material.

(ii) $A_{k}\rightarrow\infty$ as $k\rightarrow\infty$. If in addition, an estimate sequence satisfies (iii) for all $k\geq 0$, there exists some $x_{k}$ such that $f(x_{k})\leq\varphi_{k}(x_{\star})$, then we can guarantee that $f(x_{k})-f_{\star}=O(A_{k}^{-1})$.

To develop such models and the corresponding methods, three sequences of points are commonly used: (a) minimizers of our models $\varphi_{k}$ that correspond to iterates $z_{k}$ of the corresponding method; (b) a sequence $y_{k}$ of points, whose first-order information is used to update the model of the function; and (c) the iterates $x_{k}$, corresponding to the best possible $f(x_{k})$ that we can form. (The iteratives often do not correspond to the minimum of the model, $\varphi_{k}$, which is not necessarily an upper bound on the function.)

Regarding (iii), this condition ensures that the models $\varphi_{k}$ remain upper bounds on the optimal value $f_{\star}$. That is, it ensures that $f_{\star}\leq\varphi_{k}(x_{\star})$ (since $f_{\star}\leq f(x_{k})$) and hence that $\varphi_{k}(x_{\star})-f_{\star}\geq 0$. From previous bullet point, this ensures that the modeling error of $f_{\star}$ goes to asymptotically as $k$ increases. More formally, conditions (ii) and (iii) allow us to construct proofs similar to potential functions and to obtain convergence rates. That is, under (iii), we get that

and that therefore $f(x_{k})-f_{\star}\leq O(A_{k}^{-1})$. The convergence rate is thereby dictated by the rate of $A_{k}^{-1}$, which goes to by (ii).

Now, the game consists of picking appropriate sequences $\{(A_{k},\varphi_{k})\}$ that correspond to simple algorithms. We thus translate our potential function results in terms of estimate sequences.

One can observe that potential functions and estimate sequences are closely related. First, in both cases, the convergence speed is dictated by that of a scalar sequence $A_{k}^{-1}$. In fact, there is one subtle but important difference between the two approaches: whereas $\varphi_{k}(x)$ should be an increasingly good approximation of $f$ for all $x$ in the context of estimate sequences, potential functions require a model to be an increasingly good approximation of only $f_{\star}$, which is less restrictive. Hence, estimate sequences are more general but may not effectively handle situations in which the analysis actually requires having a weaker model that holds only on $f_{\star}$, and not of $f(x)$, for all $x$. We make this discussion more concrete via three examples, namely gradient descent, Nesterov’s method, and the optimized gradient method.

Gradient descent: the potential inequality from Theorem 4.2.1 actually holds for all $x$, and not only $x_{\star}$, as the proof does not exploit the optimality of $x_{\star}$. That is, it is proved that:

and $A_{k}=A_{0}+k$ (with $A_{0}>0$) is an estimate sequence for gradient descent.

and $A_{k}=A_{k-1}+\tfrac{1+\sqrt{4A_{k-1}+1}}{2}$ (with $A_{0}>0$) is an estimate sequence for Nesterov’s method.

Optimized gradient method: the potential inequality from Theorem 4.3.1 exploits the fact that $x_{\star}$ is an optimal point. Indeed, the proof relies on

which is an instance of Equation (4.4) exploiting $\nabla f(x_{\star})=0$. This does not mean that there is no estimate sequence-type model of the function as the algorithm proceeds, but the potential does not directly correspond to one. Alternatively, one can interpret

as an increasingly good model of $f_{\star}$ (i.e., it is an increasingly good approximation of $f(x)$ for all $x$ such that $\nabla f(x)=0$).

A similar conclusion holds for the conjugate gradient method (CG), from Appendix B.2. We are not aware of any estimate sequence that can be used to prove that CG reaches the lower bound from Theorem 4.3.7.

These discussions can be extended to the strongly convex setting, which we now address.

5 Acceleration under Strong Convexity

As in the smooth convex case, the smooth strongly convex case can be studied through potential functions. There are many ways to prove convergence rates for this setting, but we only consider one that allows us to recover the $\mu=0$ case as its limit such that the results are well-defined even in degenerate cases. The next proof is essentially the same as that for the smooth convex case in Theorem 4.2.1, and the same inequalities are used, with strong convexity instead of convexity. The potential is only slightly modified, thereby allowing $A_{k}$ to have a geometric growth rate:

For notational convenience, we use $q=\tfrac{\mu}{L}$ to denote the inverse condition ratio. This quantity plays a key role in the geometric convergence of first-order methods in the presence of strong convexity.

with $x_{k+1}=x_{k}-\tfrac{1}{L}\nabla f(x_{k})$, $A_{k+1}=(1+A_{k})/(1-q)$, and $q=\tfrac{\mu}{L}$.

The proof consists of performing a weighted sum of the following inequalities.

Strong convexity of $f$ between $x_{k}$ and $x_{\star}$, with weight $\lambda_{1}=A_{k+1}-A_{k}$:

Smoothness of $f$ between $x_{k}$ and $x_{k+1}$ with weight $\lambda_{2}=A_{k+1}$

This weighted sum yields a valid inequality:

Using $x_{k+1}=x_{k}-\tfrac{1}{L}\nabla f(x_{k})$, this inequality can be rewritten exactly as

The desired inequality follows from $A_{k+1}=(1+A_{k})/(1-q)$ and the sign of $A_{k}$, making one of the last two terms nonpositive and the other equal to zero, thus reaching

From this theorem, we observe that adding strong convexity to the problem allows $A_{k}$ to follow a geometric rate given by $(1-q)^{-1}$ (where we again denote by $q=\tfrac{\mu}{L}$ the inverse condition number). The corresponding iteration complexity of gradient descent to find an approximate solution $f(x_{k})-f_{\star}\leq\epsilon$ for smooth strongly convex minimization is therefore $O(\tfrac{L}{\mu}\log\tfrac{1}{\epsilon})$. This rate is essentially tight, as can be verified on quadratic functions (see, e.g., Section 2), and it follows from the following corollary whose result can be translated to iteration complexity using the same arguments as in Corollary 2.3.3.

with the inverse condition number $q=\tfrac{\mu}{L}$.

Following the reasoning of (4.5), we recursively use Theorem 4.2.1 starting with $A_{0}=0$; that is,

and we notice that the recurrence equation $A_{k+1}=(A_{k}+1)/(1-q)$ has the solution $A_{k}=((1-q)^{-k}-1)/q$. The final bound is obtained by using $f(x_{N})-f_{\star}\leq\tfrac{L\|x_{0}-x_{\star}\|^{2}_{2}}{2A_{N}}$ again.

Note that as $\mu\rightarrow 0$, the result of Corollary 4.5.3 tends to that of Corollary 4.2.3.

As in the smooth convex case, one can derive lower complexity bounds for smooth strongly convex optimization. Using the lower bounds from smooth strongly convex quadratic minimization (for which Chebyshev’s methods have optimal iteration complexity), one can conclude that no black-box first-order method can behave better than $f(x_{k})-f_{\star}=O(\rho^{k})$ with $\rho=\tfrac{(1-\sqrt{q})^{2}}{(1+\sqrt{q})^{2}}$ (see Section 2). In other words, lower complexity bounds from the quadratic optimization setting have the form $f(x_{k})-f_{\star}=\Omega(\rho^{k})$. We refer the reader to Nest03a, nemirovsky1992information for more details.

For smooth strongly convex problems beyond quadratics, this lower bound can be improved to $f(x_{k})-f_{\star}=\Omega((1-\sqrt{q})^{2k})$ as provided in [drori2021exact, Corollary 4]. In this context, we see that Nesterov’s acceleration satisfies

That is, it has an $O(\sqrt{\tfrac{L}{\mu}}\log\tfrac{1}{\epsilon})$ iteration complexity (using similar simplifications as those of Corollary 2.3.3), reaching the lower complexity bound up to a constant factor. As for the optimized gradient method provided in Section 4.3, an optimal method for the smooth strongly convex case is detailed in Section 4.6.1, and it can be shown to match exactly the corresponding worst-case lower complexity bound.

5.2 Acceleration for Smooth Strongly Convex Objectives

To adapt our proofs of convergence of accelerated methods to the strongly convex case, we need to make a small adjustment to the shape of the previous accelerated method