Understanding the Acceleration Phenomenon via High-Resolution Differential Equations

Bin Shi, Simon S. Du, Michael I. Jordan, Weijie J. Su

Introduction

Machine learning has become one of the major application areas for optimization algorithms during the past decade. While there have been many kinds of applications, to a wide variety of problems, the most prominent applications have involved large-scale problems in which the objective function is the sum over terms associated with individual data, such that stochastic gradients can be computed cheaply, while gradients are much more expensive and the computation (and/or storage) of Hessians is often infeasible. In this setting, simple first-order gradient descent algorithms have become dominant, and the effort to make these algorithms applicable to a broad range of machine learning problems has triggered a flurry of new research in optimization, both methodological and theoretical.

We will be considering unconstrained minimization problems,

where ff is a smooth convex function. Perhaps the simplest first-order method for solving this problem is gradient descent. Taking a fixed step size ss, gradient descent is implemented as the recursive rule

where α>0\alpha>0 is the momentum coefficient. While the heavy-ball method provably attains a faster rate of local convergence than gradient descent near a minimum of ff, it does not come with global guarantees. Indeed, [LRP16] demonstrate that even for strongly convex functions the method can fail to converge for some choices of the step size.[Pol64] considers s=4/(L+μ)2s=4/(\sqrt{L}+\sqrt{\mu})^{2} and α=(1−μs)2\alpha=(1-\sqrt{\mu s})^{2}. This momentum coefficient is basically the same as the choice α=1−μs1+μs\alpha=\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}} (adopted starting from Section 1.1) if ss is small.

The next major development in first-order methodology was due to Nesterov, who discovered a class of accelerated gradient methods that have a faster global convergence rate than gradient descent [Nes83, Nes13]. For a μ\mu-strongly convex objective ff with LL-Lipschitz gradients, Nesterov’s accelerated gradient method (NAG-SC) involves the following pair of update equations:

starting from x0x_{0} and x1=x0−2s∇f(x0)1+μsx_{1}=x_{0}-\frac{2s\nabla f(x_{0})}{1+\sqrt{\mu s}}. Like the heavy-ball method, NAG-SC blends gradient and momentum contributions into its update direction, but defines a specific momentum coefficient 1−μs1+μs\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}}. Nesterov also developed the estimate sequence technique to prove that NAG-SC achieves an accelerated linear convergence rate:

if the step size satisfies 0<s≤1/L0<s\leq 1/L. Moreover, for a (weakly) convex objective ff with LL-Lipschitz gradients, Nesterov defined a related accelerated gradient method (NAG-C), that takes the following form:

for any step size s≤1/Ls\leq 1/L. Under an oracle model of optimization complexity, the convergence rates achieved by NAG-SC and NAG-C are optimal for smooth strongly convex functions and smooth convex functions, respectively [NY83].

Throughout the present paper, we let α=1−μs1+μs\alpha=\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}} and x1=x0−2s∇f(x0)1+μsx_{1}=x_{0}-\frac{2s\nabla f(x_{0})}{1+\sqrt{\mu s}} to define a specific implementation of the heavy-ball method in (1.2). This choice of the momentum coefficient and the second initial point renders the heavy-ball method and NAG-SC identical except for the last (small) term in (1.4). Despite their close resemblance, however, the two methods are in fact fundamentally different, with contrasting convergence results (see, for example, [Bub15]). Notably, the former algorithm in general only achieves local acceleration, while the latter achieves acceleration method for all initial values of the iterate [LRP16]. As a numerical illustration, Figure 1 presents the trajectories that arise from the two methods when minimizing an ill-conditioned convex quadratic function. We see that the heavy-ball method exhibits pronounced oscillations throughout the iterations, whereas NAG-SC is monotone in the function value once the iteration counter exceeds 5050.

This striking difference between the two methods can only be attributed to the last term in (1.4):

which we refer to henceforth as the gradient correctionThe gradient correction for NAG-C is kk+3⋅s(∇f(xk)−∇f(xk−1))\frac{k}{k+3}\cdot s(\nabla f(x_{k})-\nabla f(x_{k-1})), as seen from the single-variable form of NAG-C: xk+1=xk+kk+3(xk−xk−1)−s∇f(xk)−kk+3⋅s(∇f(xk)−∇f(xk−1))x_{k+1}=x_{k}+\frac{k}{k+3}(x_{k}-x_{k-1})-s\nabla f(x_{k})-\frac{k}{k+3}\cdot s(\nabla f(x_{k})-\nabla f(x_{k-1})).. This term corrects the update direction in NAG-SC by contrasting the gradients at consecutive iterates. Although an essential ingredient in NAG-SC, the effect of the gradient correction is unclear from the vantage point of the estimate-sequence technique used in Nesterov’s proof. Accordingly, while the estimate-sequence technique delivers a proof of acceleration for NAG-SC, it does not explain why the absence of the gradient correction prevents the heavy-ball method from achieving acceleration for strongly convex functions.

A recent line of research has taken a different point of view on the theoretical analysis of acceleration, formulating the problem in continuous time and obtaining algorithms via discretization [SBC14, KBB15, WWJ16]). This can be done by taking continuous-time limits of existing algorithms to obtain ordinary differential equations (ODEs) that can be analyzed using the rich toolbox associated with ODEs, including Lyapunov functionsOne can think of the Lyapunov function as a generalization of the idea of the energy of a system. Then the method studies stability by looking at the rate of change of this measure of energy.. For instance, [SBC16] shows that

with initial conditions X(0)=x0X(0)=x_{0} and X˙(0)=0\dot{X}(0)=0, is the exact limit of NAG-C (1.5) by taking the step size s→0s\rightarrow 0. Alternatively, the starting point may be a Lagrangian or Hamiltonian framework [WWJ16]. In either case, the continuous-time perspective not only provides analytical power and intuition, but it also provides design tools for new accelerated algorithms.

Unfortunately, existing continuous-time formulations of acceleration stop short of differentiating between the heavy-ball method and NAG-SC. In particular, these two methods have the same limiting ODE (see, for example, [WRJ16]):

and, as a consequence, this ODE does not provide any insight into the stronger convergence results for NAG-SC as compared to the heavy-ball method. As will be shown in Section 2, this is because the gradient correction 1−μs1+μss(∇f(xk)−∇f(xk−1))=O(s1.5)\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}}s\left(\nabla f(x_{k})-\nabla f(x_{k-1})\right)=O(s^{1.5}) is an order-of-magnitude smaller than the other terms in (1.4) if s=o(1)s=o(1). Consequently, the gradient correction is not reflected in the low-resolution ODE (1.9) associated with NAG-SC, which is derived by simply taking s→0s\rightarrow 0 in both (1.2) and (1.4).

2 Overview of Contributions

Just as there is not a singled preferred way to discretize a differential equation, there is not a single preferred way to take a continuous-time limit of a difference equation. Inspired by dimensional-analysis strategies widely used in fluid mechanics in which physical phenomena are investigated at multiple scales via the inclusion of various orders of perturbations [Ped13], we propose to incorporate O(s)O(\sqrt{s}) terms into the limiting process for obtaining an ODE, including the (Hessian-driven) gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} in (1.7). This will yield high-resolution ODEs that differentiate between the NAG methods and the heavy-ball method.

We list the high-resolution ODEs that we derive in the paper hereWe note that the form of the initial conditions is fixed for each ODE throughout the paper. For example, while x0x_{0} is arbitrary, X(0)X(0) and X˙(0)\dot{X}(0) must always be equal to x0x_{0} and −2sf(x0)/(1+μs)-2\sqrt{s}f(x_{0})/(1+\sqrt{\mu s}) respectively in the high-resolution ODE of the heavy-ball method. This is in accordance with the choice of α=1−μs1+μs\alpha=\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}} and x1=x0−2s∇f(x0)1+μsx_{1}=x_{0}-\frac{2s\nabla f(x_{0})}{1+\sqrt{\mu s}}.:

The high-resolution ODE for the heavy-ball method (1.2):

with X(0)=x0X(0)=x_{0} and X˙(0)=−2s∇f(x0)1+μs\dot{X}(0)=-\frac{2\sqrt{s}\nabla f(x_{0})}{1+\sqrt{\mu s}}.

The high-resolution ODE for NAG-SC (1.3):

with X(0)=x0X(0)=x_{0} and X˙(0)=−2s∇f(x0)1+μs\dot{X}(0)=-\frac{2\sqrt{s}\nabla f(x_{0})}{1+\sqrt{\mu s}}.

for t≥3s/2t\geq 3\sqrt{s}/2, with X(3s/2)=x0X(3\sqrt{s}/2)=x_{0} and X˙(3s/2)=−s∇f(x0)\dot{X}(3\sqrt{s}/2)=-\sqrt{s}\nabla f(x_{0}).

High-resolution ODEs are more accurate continuous-time counterparts for the corresponding discrete algorithms than low-resolution ODEs, thus allowing for a better characterization of the accelerated methods. This is illustrated in Figure 2, which presents trajectories and convergence of the discrete methods, and the low- and high-resolution ODEs. For both NAGs, the high-resolution ODEs are in much better agreement with the discrete methods than the low-resolution ODEsNote that for the heavy-ball method, the trajectories of the high-resolution ODE and the low-resolution ODE are almost identical.. Moreover, for NAG-SC, its high-resolution ODE captures the non-oscillation pattern while the low-resolution ODE does not.

The three new ODEs include O(s)O(\sqrt{s}) terms that are not present in the corresponding low-resolution ODEs (compare, for example, (1.12) and (1.8)). Note also that if we let s→0s\rightarrow 0, each high-resolution ODE reduces to its low-resolution counterpart. Thus, the difference between the heavy-ball method and NAG-SC is reflected only in their high-resolution ODEs: the gradient correction (1.7) of NAG-SC is preserved only in its high-resolution ODE in the form s∇2f(X(t))X˙(t)\sqrt{s}\nabla^{2}f(X(t))\dot{X}(t). This term, which we refer to as the (Hessian-driven) gradient correction, is connected with the discrete gradient correction by the approximate identity:

for small ss, with the identification t=kst=k\sqrt{s}. The gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} in NAG-C arises in the same fashionHenceforth, the dependence of XX on tt is suppressed when clear from the context.. Interestingly, although both NAGs are first-order methods, their gradient corrections brings in second-order information from the objective function.

Despite being small, the gradient correction has a fundamental effect on the behavior of both NAGs, and this effect is revealed by inspection of the high-resolution ODEs. We provide two illustrations of this.

Effect of the gradient correction in acceleration. Viewing the coefficient of X˙\dot{X} as a damping ratio, the ratio 2μ+s∇2f(X)2\sqrt{\mu}+\sqrt{s}\nabla^{2}f(X) of X˙\dot{X} in the high-resolution ODE (1.11) of NAG-SC is adaptive to the position XX, in contrast to the fixed damping ratio 2μ2\sqrt{\mu} in the ODE (1.10) for the heavy-ball method. To appreciate the effect of this adaptivity, imagine that the velocity X˙\dot{X} is highly correlated with an eigenvector of ∇2f(X)\nabla^{2}f(X) with a large eigenvalue, such that the large friction (2μ+s∇2f(X))X˙(2\sqrt{\mu}+\sqrt{s}\nabla^{2}f(X))\dot{X} effectively “decelerates” along the trajectory of the ODE (1.11) of NAG-SC. This feature of NAG-SC is appealing as taking a cautious step in the presence of high curvature generally helps avoid oscillations. Figure 1 and the left plot of Figure 2 confirm the superiority of NAG-SC over the heavy-ball method in this respect.

If we can translate this argument to the discrete case we can understand why NAG-SC achieves acceleration globally for strongly convex functions but the heavy-ball method does not. We will be able to make this translation by leveraging the high-resolution ODEs to construct discrete-time Lyapunov functions that allow maximal step sizes to be characterized for the NAG-SC and the heavy-ball method. The detailed analyses is given in Section 3.

Effect of gradient correction in gradient norm minimization. We will also show how to exploit the high-resolution ODE of NAG-C to construct a continuous-time Lyapunov function to analyze convergence in the setting of a smooth convex objective with LL-Lipschitz gradients. Interestingly, the time derivative of the Lyapunov function is not only negative, but it is smaller than −O(st2∥∇f(X)∥2)-O(\sqrt{s}t^{2}\|\nabla f(X)\|^{2}). This bound arises from the gradient correction and, indeed, it cannot be obtained from the Lyapunov function studied in the low-resolution case by [SBC16]. This finer characterization in the high-resolution case allows us to establish a new phenomenon:

That is, we discover that NAG-C achieves an inverse cubic rate for minimizing the squared gradient norm. By comparison, from (1.6) and the LL-Lipschitz continuity of ∇f\nabla f we can only show that ∥∇f(xk)∥2≤O(L2/k2)\left\|\nabla f(x_{k})\right\|^{2}\leq O\left(L^{2}/k^{2}\right). See Section 4 for further elaboration on this cubic rate for NAG-C.

As we will see, the high-resolution ODEs are based on a phase-space representation that provides a systematic framework for translating from continuous-time Lyapunov functions to discrete-time Lyapunov functions. In sharp contrast, the process for obtaining a discrete-time Lyapunov function for low-resolution ODEs presented by [SBC16] relies on “algebraic tricks” (see, for example, Theorem 6 of [SBC16]). On a related note, a Hessian-driven damping term also appears in ODEs for modeling Newton’s method [AABR02, AMR12, APR16].

3 Related Work

There is a long history of using ODEs to analyze optimization methods [HM12, Sch00, Fio05]. Recently, the work of [SBC14, SBC16] has sparked a renewed interest in leveraging continuous dynamical systems to understand and design first-order methods and to provide more intuitive proofs for the discrete methods. Below is a rather incomplete review of recent work that uses continuous-time dynamical systems to study accelerated methods.

In the work of [WWJ16, WRJ16, BJW18], Lagrangian and Hamiltonian frameworks are used to generate a large class of continuous-time ODEs for a unified treatment of accelerated gradient-based methods. Indeed, [WWJ16] extend NAG-C to non-Euclidean settings, mirror descent and accelerated higher-order gradient methods, all from a single “Bregman Lagrangian.” In [WRJ16], the connection between ODEs and discrete algorithms is further strengthened by establishing an equivalence between the estimate sequence technique and Lyapunov function techniques, allowing for a principled analysis of the discretization of continuous-time ODEs. Recent papers have considered symplectic [BJW18] and Runge–Kutta [ZMSJ18] schemes for discretization of the low-resolution ODEs.

An ODE-based analysis of mirror descent has been pursued in another line of work by [KBB15, KBB16, KB17], delivering new connections between acceleration and constrained optimization, averaging and stochastic mirror descent.

In addition to the perspective of continuous-time dynamical systems, there has also been work on the acceleration from a control-theoretic point of view [LRP16, HL17, FRMP18] and from a geometric point of view [BLS15, CML17]. See also [OC15, FB15, GL16, DO17, LMH18, DFR18] for a number of other recent contributions to the study of the acceleration phenomenon.

4 Organization and Notation

The remainder of the paper is organized as follows. In Section 2, we briefly introduce our high-resolution ODE-based analysis framework. This framework is used in Section 3 to study the heavy-ball method and NAG-SC for smooth strongly convex functions. In Section 4, we turn our focus to NAG-C for a general smooth convex objective. In Section 5 we derive some extensions of NAG-C. We conclude the paper in Section 6 with a list of future research directions. Most technical proofs are deferred to the Appendix.

The High-Resolution ODE Framework

This section introduces a high-resolution ODE framework for analyzing gradient-based methods, with NAG-SC being a guiding example. Given a (discrete) optimization algorithm, the first step in this framework is to derive a high-resolution ODE using dimensional analysis, the next step is to construct a continuous-time Lyapunov function to analyze properties of the ODE, the third step is to derive a discrete-time Lyapunov function from its continuous counterpart and the last step is to translate properties of the ODE into that of the original algorithm. The overall framework is illustrated in Figure 3.

Our focus is on the single-variable form (1.4) of NAG-SC. For any nonnegative integer kk, let tk=kst_{k}=k\sqrt{s} and assume xk=X(tk)x_{k}=X(t_{k}) for some sufficiently smooth curve X(t)X(t). Performing a Taylor expansion in powers of s\sqrt{s}, we get

We now use a Taylor expansion for the gradient correction, which gives

Multiplying both sides of (1.4) by 1+μs1−μs⋅1s\frac{1+\sqrt{\mu s}}{1-\sqrt{\mu s}}\cdot\frac{1}{s} and rearranging the equality, we can rewrite NAG-SC as

Next, plugging (2.1) and (2.2) into (2.3), we haveNote that we use the approximation xk+1+xk−1−2xks=X¨(tk)+O(s)\frac{x_{k+1}+x_{k-1}-2x_{k}}{s}=\ddot{X}(t_{k})+O(s), whereas [SBC16] relies on the low-accuracy Taylor expansion xk+1+xk−1−2xks=X¨(tk)+o(1)\frac{x_{k+1}+x_{k-1}-2x_{k}}{s}=\ddot{X}(t_{k})+o(1) in the derivation of the low-resolution ODE of NAG-C. We illustrate this derivation of the three low-resolution ODEs in Appendix A.2; they can be compared to the high-resolution ODEs that we derive here.

Multiplying both sides of the last display by 1−μs1-\sqrt{\mu s}, we obtain the following high-resolution ODE of NAG-SC:

where we ignore any O(s)O(s) terms but retain the O(s)O(\sqrt{s}) terms (note that (1−μs)s=s+O(s)(1-\sqrt{\mu s})\sqrt{s}=\sqrt{s}+O(s)).

Our analysis is inspired by dimensional analysis [Ped13], a strategy widely used in physics to construct a series of differential equations that involve increasingly high-order terms corresponding to small perturbations. In more detail, taking a small ss, one first derives a differential equation that consists only of O(1)O(1) terms, then derives a differential equation consisting of both O(1)O(1) and O(s)O(\sqrt{s}), and next, one proceeds to obtain a differential equation consisting of O(1),O(s)O(1),O(\sqrt{s}) and O(s)O(s) terms. High-order terms in powers of s\sqrt{s} are introduced sequentially until the main characteristics of the original algorithms have been extracted from the resulting approximating differential equation. Thus, we aim to understand Nesterov acceleration by incorporating O(s)O(\sqrt{s}) terms into the ODE, including the (Hessian-driven) gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} which results from the (discrete) gradient correction (1.7) in the single-variable form (1.4) of NAG-SC. We also show (see Appendix A.1 for the detailed derivation) that this O(s)O(\sqrt{s}) term appears in the high-resolution ODE of NAG-C, but is not found in the high-resolution ODE of the heavy-ball method.

In fact, Proposititon 2.1 holds for T=∞T=\infty because both the discrete iterates and the ODE trajectories converge to the unique minimizer when the objective is stongly convex.

The proofs of these propositions are given in Appendix A.3.1 and Appendix A.3.2.

Step 2: Analyzing ODEs Using Lyapunov Functions

With these high-resolution ODEs in place, the next step is to construct Lyapunov functions for analyzing the dynamics of the corresponding ODEs, as is done in previous work [SBC16, WRJ16, LRP16]. For NAG-SC, we consider the Lyapunov function

The first and second terms (1+μs)(f(X)−f(x⋆))(1+\sqrt{\mu s})\left(f(X)-f(x^{\star})\right) and 14∥X˙∥2\frac{1}{4}\|\dot{X}\|^{2} can be regarded, respectively, as the potential energy and kinetic energy, and the last term is a mix. For the mixed term, it is interesting to note that the time derivative of X˙+2μ(X−x⋆)+s∇f(X)\dot{X}+2\sqrt{\mu}(X-x^{\star})+\sqrt{s}\nabla f(X) equals −(1+μs)∇f(X)-(1+\sqrt{\mu s})\nabla f(X).

The differentiability of E(t)\mathcal{E}(t) will allow us to investigate properties of the ODE (1.11) in a principled manner. For example, we will show that E(t)\mathcal{E}(t) decreases exponentially along the trajectories of (1.11), recovering the accelerated linear convergence rate of NAG-SC. Furthermore, a comparison between the Lyapunov function of NAG-SC and that of the heavy-ball method will explain why the gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} yields acceleration in the former case. This is discussed in Section 3.1.

Step 3: Constructing Discrete Lyapunov Functions

Our framework make it possible to translate continuous Lyapunov functions into discrete Lyapunov functions via a phase-space representation (see, for example, [Arn13]). We illustrate the procedure in the case of NAG-SC. The first step is formulate explicit position and velocity updates:

where the velocity variable vkv_{k} is defined as:

The initial velocity is v0=−2s1+μs∇f(x0)v_{0}=-\frac{2\sqrt{s}}{1+\sqrt{\mu s}}\nabla f(x_{0}). Interestingly, this phase-space representation has the flavor of symplectic discretization, in the sense that the update for xk−xk−1x_{k}-x_{k-1} is explicit (it only depends on the last iterate vk−1v_{k-1}) while the update for vk−vk−1v_{k}-v_{k-1} is implicit (it depends on the current iterates xkx_{k} and vkv_{k})Although this suggestion is a heuristic one, it is also possible to rigorously derive a symplectic integrator of the high-resolution ODE of NAG-SC; this integrator has the form: xk−xk−1=svk−1\displaystyle x_{k}-x_{k-1}=\sqrt{s}v_{k-1} vk−vk−1=−2μsvk−s∇2f(xk)vk−(1+μs)s∇f(xk).\displaystyle v_{k}-v_{k-1}=-2\sqrt{\mu s}v_{k}-s\nabla^{2}f(x_{k})v_{k}-(1+\sqrt{\mu s})\sqrt{s}\nabla f(x_{k}). .

The representation (2.5) suggests translating the continuous-time Lyapunov function (2.4) into a discrete-time Lyapunov function of the following form:

by replacing continuous terms (e.g., X˙\dot{X}) by their discrete counterparts (e.g., vkv_{k}). Akin to the continuous (2.4), here I\mathbf{I}, II\mathbf{II}, and III\mathbf{III} correspond to potential energy, kinetic energy, and mixed energy, respectively, from a mechanical perspective. To better appreciate this translation, note that the factor 1+μs1−μs\frac{1+\sqrt{\mu s}}{1-\sqrt{\mu s}} in I\mathbf{I} results from the term 1+μs1−μss∇f(xk)\frac{1+\sqrt{\mu s}}{1-\sqrt{\mu s}}\sqrt{s}\nabla f(x_{k}) in (2.5). Likewise, 2μ1−μs\frac{2\sqrt{\mu}}{1-\sqrt{\mu s}} in III\mathbf{III} is from the term 2μs1−μsvk\frac{2\sqrt{\mu s}}{1-\sqrt{\mu s}}v_{k} in (2.5). The need for the final (small) negative term is technical; we discuss it in Section 3.2.

Step 4: Analyzing Algorithms Using Discrete Lyapunov Functions

The last step is to map properties of high-resolution ODEs to corresponding properties of optimization methods. This step closely mimics Step 2 except that now the object is a discrete algorithm and the tool is a discrete Lyapunov function such as (2.6). Given that Step 2 has been performed, this translation is conceptually straightforward, albeit often calculation-intensive. For example, using the discrete Lyapunov function (2.6), we will recover the optimal linear rate of NAG-SC and gain insights into the fundamental effect of the gradient correction in accelerating NAG-SC. In addition, NAG-C is shown to minimize the squared gradient norm at an inverse cubic rate by a simple analysis of the decreasing rate of its discrete Lyapunov function.

Gradient Correction for Acceleration

Throughout this section, the strategy is to analyze the two methods in parallel, thereby highlighting the differences between the two methods. In particular, the comparison will demonstrate the vital role of the gradient correction, namely 1−μs1+μs⋅s(∇f(xk)−∇f(xk−1))\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}}\cdot s\left(\nabla f(x_{k})-\nabla f(x_{k-1})\right) in the discrete case and s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} in the ODE case, in making NAG-SC an accelerated method.

The following theorem characterizes the convergence rate of the high-resolution ODE corresponding to NAG-SC.

The proof of Theorem 1 is based on analyzing the Lyapunov function E(t)\mathcal{E}(t) for the high-resolution ODE of NAG-SC. Recall that E(t)\mathcal{E}(t) defined in (2.4) is

The next lemma states the key property we need from this Lyapunov function

The proof of this theorem relies on Lemma 3.1 through the inequality E˙(t)≤−μ4E(t)\dot{\mathcal{E}}(t)\leq-\frac{\sqrt{\mu}}{4}\mathcal{E}(t). The term s2(∥∇f(X)∥2+X˙⊤∇2f(X)X˙)≥0\frac{\sqrt{s}}{2}(\left\|\nabla f(X)\right\|^{2}+\dot{X}^{\top}\nabla^{2}f(X)\dot{X})\geq 0 plays no role at the moment, but Section 3.2 will shed light on its profound effect in the discretization of the high-resolution ODE of NAG-SC.

Lemma 3.1 implies E˙(t)≤−μ4E(t)\dot{\mathcal{E}}(t)\leq-\frac{\sqrt{\mu}}{4}\mathcal{E}(t), which amounts to

Recognizing the initial conditions X(0)=x0X(0)=x_{0} and X˙(0)=−2s∇f(x0)1+μs\dot{X}(0)=-\frac{2\sqrt{s}\nabla f(x_{0})}{1+\sqrt{\mu s}}, we write (3.2) as

Since f∈Sμ,L2f\in\mathcal{S}_{\mu,L}^{2}, we have that ∥∇f(x0)∥≤L∥x0−x⋆∥\|\nabla f(x_{0})\|\leq L\|x_{0}-x^{\star}\| and f(x0)−f(x⋆)≤L∥x0−x⋆∥2/2f(x_{0})-f(x^{\star})\leq L\|x_{0}-x^{\star}\|^{2}/2. Together with the Cauchy–Schwarz inequality, the two inequalities yield

Furthermore, a bit of analysis reveals that

since μs≤μ/L≤1\mu s\leq\mu/L\leq 1, and this step completes the proof of Theorem 1. ∎

We now consider the heavy-ball method (1.2). Recall that the momentum coefficient α\alpha is set to 1−μs1+μs\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}}. The following theorem characterizes the rate of convergence of this method.

As in the case of NAG-SC, the proof of Theorem 2 is based on a Lyapunov function:

which is the same as the Lyapunov function (2.4) for NAG-SC except for the lack of the s∇f(X)\sqrt{s}\nabla f(X) term. In particular, (2.4) and (3.3) are identical if s=0s=0. The following lemma considers the decay rate of (3.3).

in the high-resolution ODE of the heavy-ball method and using the LL-smoothness of ∇f\nabla f, Lemma 3.2 yields

if the step size s≤1/Ls\leq 1/L. Finally, since 0<μs≤μ/L≤10<\mu s\leq\mu/L\leq 1, the coefficient satisfies

The proofs of Lemma 3.1 and Lemma 3.2 share similar ideas. In view of this, we present only the proof of the former here, deferring the proof of Lemma 3.2 to Appendix B.1.

Along trajectories of (1.11) the Lyapunov function (2.4) satisfies

Furthermore, ⟨∇f(X),X−x⋆⟩\left\langle\nabla f(X),X-x^{\star}\right\rangle is greater than or equal to both f(X)−f(x⋆)+μ2∥X−x⋆∥2f(X)-f(x^{\star})+\frac{\mu}{2}\|X-x^{\star}\|^{2} and μ∥X−x⋆∥2\mu\|X-x^{\star}\|^{2} due to the μ\mu-strong convexity of ff. This yields

which together with (3.4) suggests that the time derivative of this Lyapunov function can be bounded as

Next, the Cauchy–Schwarz inequality yields

Combining (3.5) and (3.6) completes the proof of the theorem.

The only inequality in (3.4) is due to the term s2(∥∇f(X)∥2+X˙⊤∇2f(X)X˙)\frac{\sqrt{s}}{2}(\left\|\nabla f(X)\right\|^{2}+\dot{X}^{\top}\nabla^{2}f(X)\dot{X}), which is discussed right after the statement of Lemma 3.1. This term results from the gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} in the NAG-SC ODE. For comparison, this term does not appear in Lemma 3.2 in the case of the heavy-ball method as its ODE does not include the gradient correction and, accordingly, its Lyapunov function (3.3) is free of the s∇f(X)\sqrt{s}\nabla f(X) term.

2 The Discrete Case

In brief, the theorem states that log⁡(f(xk)−f(x⋆))≤−O(kμ/L)\log(f(x_{k})-f(x^{\star}))\leq-O(k\sqrt{\mu/L}), which matches the optimal rate for minimizing smooth strongly convex functions using only first-order information [Nes13]. More precisely, [Nes13] shows that f(xk)−f(x⋆)=O((1−μ/L)k)f(x_{k})-f(x^{\star})=O((1-\sqrt{\mu/L})^{k}) by taking s=1/Ls=1/L in NAG-SC. Although this optimal rate of NAG-SC is well known in the litetature, this is the first Lyapunov-function-based proof of this result.

As indicated in Section 2, the proof of Theorem 3 rests on the discrete Lyapunov function (2.6):

Recall that this functional is derived by writing NAG-SC in the phase-space representation (2.5). Analogous to Lemma 3.1, the following lemma gives an upper bound on the difference E(k+1)−E(k)\mathcal{E}(k+1)-\mathcal{E}(k).

The form of the inequality ensured by Lemma 3.4 is consistent with that of Lemma 3.1. Alternatively, it can be written as E(k+1)≤11+μs6E(k)\mathcal{E}(k+1)\leq\frac{1}{1+\frac{\sqrt{\mu s}}{6}}\mathcal{E}(k). With Lemma 3.4 in place, we give the proof of Theorem 3.

Next, we inductively apply Lemma 3.4, yielding

Recognizing the initial velocity v0=−2s∇f(x0)1+μsv_{0}=-\frac{2\sqrt{s}\nabla f(x_{0})}{1+\sqrt{\mu s}} in NAG-SC, one can show that

Taking s=1/(4L)s=1/(4L) in (3.9), it follows from (3.7) and (3.8) that

Here the constant factor Cμ/LC_{\mu/L} is a short-hand for

which is less than five by making use of the fact that μ/L≤1\mu/L\leq 1. This completes the proof.

We now turn to the heavy-ball method (1.2). Recall that α=1−μs1+μs\alpha=\frac{1-\sqrt{\mu s}}{1+\sqrt{\mu s}} and x1=x0−2s∇f(x0)1+μsx_{1}=x_{0}-\frac{2s\nabla f(x_{0})}{1+\sqrt{\mu s}}.

In addition to allowing us to complete the proof of Theorem 4, Lemma 3.5 will shed light on why the heavy-ball method needs a more conservative step size. To state this lemma, we consider the discrete Lyapunov function defined as

which is derived by discretizing the continuous Lyapunov function (3.3) using the phase-space representation of the heavy-ball method:

The proof of Lemma 3.5 can be found in Appendix B.3. To apply this lemma to prove Theorem 4, we need to ensure

A sufficient and necessary condition for (3.13) is

This is because ∥∇f(xk+1)∥2≤2L(f(xk+1)−f(x⋆))\left\|\nabla f(x_{k+1})\right\|^{2}\leq 2L\left(f(x_{k+1})-f(x^{\star})\right), which can be further reduced to an equality (for example, f(x)=L2∥x∥2f(x)=\frac{L}{2}\|x\|^{2}). Thus, the step size ss must obey

In particular, the choice of s=μ16L2s=\frac{\mu}{16L^{2}} fulfills (3.14) and, as a consequence, Lemma 3.5 implies

The remainder of the proof of Theorem 4 is similar to that of Theorem 3 and is therefore omitted. As an aside, [Pol64] uses s=4/(L+μ)2s=4/(\sqrt{L}+\sqrt{\mu})^{2} for local accelerated convergence of the heavy-ball method. This choice of step size is larger than our step size s=μ16L2s=\frac{\mu}{16L^{2}}, which yields a non-accelerated but global convergence rate.

The term s2(1+μs1−μs)2∥∇f(xk+1)∥2\frac{s}{2}\left(\frac{1+\sqrt{\mu s}}{1-\sqrt{\mu s}}\right)^{2}\left\|\nabla f(x_{k+1})\right\|^{2} in (3.12) that arises from finite differencing of (3.10) is a (small) term of order O(s)O(s) and, as a consequence, this term is not reflected in Lemma 3.2. In relating to the case of NAG-SC, one would be tempted to ask why this term does not appear in Lemma 3.4. In fact, a similar term can be found in E(k+1)−E(k)\mathcal{E}(k+1)-\mathcal{E}(k) by taking a closer look at the proof of Lemma 3.4. However, this term is canceled out by the discrete version of the quadratic term s2(∥∇f(X)∥2+X˙⊤∇2f(X)X˙)\frac{\sqrt{s}}{2}(\left\|\nabla f(X)\right\|^{2}+\dot{X}^{\top}\nabla^{2}f(X)\dot{X}) in Lemma 3.1 and is, therefore, not present in the statement of Lemma 3.4. Note that this quadratic term results from the gradient correction (see Remark 3.3). In light of the above, the gradient correction is the key ingredient that allows for a larger step size in NAG-SC, which is necessary for achieving acceleration.

For completeness, we finish Section 3.2 by proving Lemma 3.4.

Using the Cauchy–Schwarz inequality, we haveSee the definition of III\mathbf{III} in (2.6).

Next, as shown in Appendix B.2, the inequality

holds for s≤1/(2L)s\leq 1/(2L). Comparing the coefficients of the same terms in (3.15) for E(k+1)\mathcal{E}(k+1) and (3.16), we conclude that the first difference of the discrete Lyapunov function (2.6) must satisfy

3 A Numerical Stability Perspective on Acceleration

As shown in Section 3.2, the gradient correction is the fundamental cause of the difference in convergence rates between the heavy-ball method and NAG-SC. This section aims to further elucidate this distinction from the viewpoint of numerical stability. A numerical scheme is said to be stable if, roughly speaking, this scheme does not magnify errors in the input data. Accordingly, we address the question of what values of the step size ss are allowed for solving the high-resolution ODEs (1.10) and (1.11) in a stable fashion. While various discretization schemes on low-resolution ODEs have been explored in [WWJ16, WRJ16, ZMSJ18], we limit our attention to the forward Euler scheme to simplify the discussion (see [SB13] for an exposition on discretization schemes).

For the heavy-ball method, the forward Euler scheme applied to (1.10) is

Using the approximation ∇f(X(t−s)+ϵ)≈∇f(X(t−s))+∇2f(X(t−s))ϵ\nabla f(X(t-\sqrt{s})+\epsilon)\approx\nabla f(X(t-\sqrt{s}))+\nabla^{2}f(X(t-\sqrt{s}))\epsilon for a small perturbation ϵ\epsilon, we get the characteristic equation of (3.17):

where I\bm{I} denotes the n×nn\times n identity matrix. The numerical stability of (3.17) requires the roots of the characteristic equation to be no larger than one in absolute value. Therefore, a necessary condition for the stability is thatThe notation A⪯BA\preceq B indicates that B−AB-A is positive semidefinite for symmetric matrices AA and BB.

By the LL-smoothness of ff, the largest singular value of ∇2f(X(t−s))\nabla^{2}f(X(t-\sqrt{s})) can be as large as LL. Therefore, (3.18) is guaranteed in the worst case analysis only if

Next, we turn to the high-resolution ODE (1.11) of NAG-SC, for which the forward Euler scheme reads

which, as earlier, suggests that the numerical stability condition of (3.20) is

This inequality is ensured by setting the step size

As constraints on the step sizes, both (3.19) and (3.21) are in agreement with the discussion in Section 3.2, albeit from a different perspective. In short, a comparison between (3.17) and (3.20) reveals that the Hessian s∇2f(X(t−s))\sqrt{s}\nabla^{2}f(X(t-\sqrt{s})) makes the forward Euler scheme for the NAG-SC ODE numerically stable with a larger step size, namely s=O(1/L)s=O(1/L). This is yet another reflection of the vital importance of the gradient correction in yielding acceleration for NAG-SC.

Gradient Correction for Gradient Norm Minimization

In this section, we extend the use of the high-resolution ODE framework to NAG-C (1.5) in the setting of minimizing an LL-smooth convex function ff. The main result is an improved rate of NAG-SC for minimizing the squared gradient norm. Indeed, we show that NAG-C achieves the O(L2/k3)O(L^{2}/k^{3}) rate of convergence for minimizing ∥∇f(xk)∥2\|\nabla f(x_{k})\|^{2}. To the best of our knowledge, this is the sharpest known bound for this problem using NAG-C without any modification. Moreover, we will show that the gradient correction in NAG-C is responsible for this rate and, as it is therefore unsurprising that this inverse cubic rate was not perceived within the low-resolution ODE frameworks such as that of [SBC16]. In Section 4.3, we propose a new accelerated method with the same rate O(L2/k3)O(L^{2}/k^{3}) and briefly discuss the benefit of the phase-space representation in simplifying technical proofs.

By taking the step size s=1/Ls=1/L, this theorem shows that

where the infimum operator is necessary as the squared gradient norm is generally not decreasing in tt. In contrast, directly combining the convergence rate of the function value (see Corollary 4.2) and inequality ∥∇f(X)∥2≤2L(f(X)−f(x⋆))\|\nabla f(X)\|^{2}\leq 2L(f(X)-f(x^{\star})) only gives a O(L/t2)O(L/t^{2}) rate for squared gradient norm minimization.

The proof of the theorem is based on the continuous Lyapunov function

which reduces to the continuous Lyapunov function in [SBC16] when setting s=0s=0.

The decreasing rate of E(t)\mathcal{E}(t) as specified in the lemma is sufficient for the proof of Theorem 5. First, note that Lemma 4.1 readily gives

where the last step is due to the fact E(t)≥0\mathcal{E}(t)\geq 0. Thus, it follows that

Recognizing the initial conditions of the ODE (1.12), we get

This bound reduces to the one claimed by Theorem 5 by only keeping the first term s(t3−t03)/3\sqrt{s}(t^{3}-t_{0}^{3})/3 in the denominator.

The gradient correction s∇2f(X)X˙\sqrt{s}\nabla^{2}f(X)\dot{X} in the high-resolution ODE (1.12) plays a pivotal role in Lemma 4.1 and is, thus, key to Theorem 5. As will be seen in the proof of the lemma, the factor ∥∇f(X)∥2\left\|\nabla f(X)\right\|^{2} in (4.2) results from the term ts∇f(X)t\sqrt{s}\nabla f(X) in the Lyapunov function (4.1), which arises from the gradient correction in the ODE (1.12). In light of this, the low-resolution ODE (1.8) of NAG-C cannot yield a result similar to Lemma 4.1 and; furthermore, we conjecture that the O(L/t3)O(\sqrt{L}/t^{3}) rate does applies to this ODE. Section 4.2 will discuss this point further in the discrete case.

In passing, it is worth pointing out that the analysis above applies to the case of s=0s=0. In this case, we have t0=0t_{0}=0, and (4.4) turns out to be

This result is similar to that of the low-resolution ODE in [SBC16]To see this, recall that [SBC16] shows that f(X(t))−f(x⋆)≤2∥x0−x⋆∥2t2f(X(t))-f(x^{\star})\leq\frac{2\|x_{0}-x^{\star}\|^{2}}{t^{2}}, where X=X(t)X=X(t) is the solution to (4.4) with s=0s=0. Using the LL-smoothness of ff, we get ∥∇f(X(t))∥2≤2L(f(X(t))−f(x⋆))≤4L∥x0−x⋆∥2t2\|\nabla f(X(t))\|^{2}\leq 2L(f(X(t))-f(x^{\star}))\leq\frac{4L\|x_{0}-x^{\star}\|^{2}}{t^{2}}..

This section is concluded with the proof of Lemma 4.1.

The time derivative of the Lyapunov function (4.1) obeys

Note that Lemma 4.1 shows E(t)\mathcal{E}(t) is a decreasing function, from which we get

by recognizing the initial conditions of the high-resolution ODE (1.12). This gives the following corollary.

Under the same assumptions as in Theorem 5, for any t>t0t>t_{0}, we have

2 The Discrete Case

for all k≥0k\geq 0. In additional, we have

Taking s=1/(3L)s=1/(3L), Theorem 6 shows that NAG-C minimizes the squared gradient norm at the rate O(L2/k3)O(L^{2}/k^{3}). This theoretical prediction is in agreement with two numerical examples illustrated in Figure 4. To our knowledge, the bound O(L2/k3)O(L^{2}/k^{3}) is sharper than any existing bounds in the literature for NAG-C for squared gradient norm minimization. In fact, the convergence result f(xk)−f(x⋆)=O(L/k2)f(x_{k})-f(x^{\star})=O(L/k^{2}) for NAG-C and the LL-smoothness of the objective immediately give ∥∇f(xk)∥2≤O(L2/k2)\|\nabla f(x_{k})\|^{2}\leq O(L^{2}/k^{2}). This well-known but loose bound can be improved by using a recent result from [AP16], which shows that a slightly modified version NAG-C satisfies f(xk)−f(x⋆)=o(L/k2)f(x_{k})-f(x^{\star})=o(L/k^{2}) (see Section 5.2 for more discussion of this improved rate). This reveals

which, however, remains looser than that of Theorem 6. In addition, the rate o(L2/k2)o(L^{2}/k^{2}) is not valid for k≤n/2k\leq n/2 and, as such, the bound o(L2/k2)o(L^{2}/k^{2}) on the squared gradient norm is dimension-dependent [AP16]. For completeness, the rate O(L2/k3)O(L^{2}/k^{3}) can be achieved by introducing an additional sequence of iterates and a more aggressive step size policy in a variant of NAG-C [GL16]. In stark contrast, our result shows that no adjustments are needed for NAG-C to yield an accelerated convergence rate for minimizing the gradient norm.

An Ω(L2/k4)\Omega(L^{2}/k^{4}) lower bound has been established by [Nes12] as the optimal convergence rate for minimizing ∥∇f∥2\|\nabla f\|^{2} with access to only first-order information. (For completeness, Appendix C.3 presents an exposition of this fundamental barrier.) In the same paper, a regularization technique is used in conjunction with NAG-SC to obtain a matching upper bound (up to a logarithmic factor). This method, however, takes as input the distance between the initial point and the minimizer, which is not practical in general [KF18].

Returning to Theorem 6, we present a proof of this theorem using a Lyapunov function argument. By way of comparison, we remark that Nesterov’s estimate sequence technique is unlikely to be useful for characterizing the convergence of the gradient norm as this technique is essentially based on local quadratic approximations. The phase-space representation of NAG-C (1.5) takes the following form:

for any initial position x0x_{0} and the initial velocity v0=−s∇f(x0)v_{0}=-\sqrt{s}\nabla f(x_{0}). This representation allows us to discretize the continuous Lyapunov function (4.1) into

The following lemma characterizes the dynamics of this Lyapunov function.

Under the assumptions of Theorem 6, we have

for k≥2k\geq 2. To show this, note that it suffices to guarantee

which is self-evident since s≤1/(3L)s\leq 1/(3L) by assumption.

Next, by a telescoping-sum argument, Lemma 4.3 leads to the following inequalities for k≥4k\geq 4:

where the second inequality is due to (4.7). To further simplify the bound, observe that

for k≥4k\geq 4. Plugging this inequality into (4.9) yields

for s≤1/(3L)s\leq 1/(3L). As a consequence of this, (4.10) gives

For completeness, Appendix C.1 proves, via a brute-force calculation, that ∥∇f(x0)∥2,∥∇f(x1)∥2,∥∇f(x2)∥2\left\|\nabla f(x_{0})\right\|^{2},\left\|\nabla f(x_{1})\right\|^{2},\left\|\nabla f(x_{2})\right\|^{2}, and ∥∇f(x3)∥2\left\|\nabla f(x_{3})\right\|^{2} are all bounded above by the right-hand side of (4.11). This completes the proof of the first inequality claimed by Theorem 6.

For the second claim in Theorem 6, the definition of the Lyapunov function and its decreasing property ensured by (4.7) implies

for all k≥2k\geq 2. Appendix C.1 establishes that f(x0)−f(x⋆)f(x_{0})-f(x^{\star}) and f(x1)−f(x⋆)f(x_{1})-f(x^{\star}) are bounded by the right-hand side of (4.12). This completes the proof.

The difference of the Lyapunov function (4.6) satisfies

which follows from the phase-space representation (4.5). Rearranging the identity for E(k+1)−E(k)\mathcal{E}(k+1)-\mathcal{E}(k), we get

The next step is to recognize that the convexity and the LL-smoothness of ff gives

Plugging these two inequalities into (4.14), we have

where the second inequality uses the fact that ⟨∇f(xk+1),xk+1−x⋆⟩≥0\left\langle\nabla f(x_{k+1}),x_{k+1}-x^{\star}\right\rangle\geq 0.

To further bound E(k+1)−E(k)\mathcal{E}(k+1)-\mathcal{E}(k), making use of (4.13) with k+1k+1 in place of kk, we get

In passing, we remark that the gradient correction sheds light on the superiority of the high-resolution ODE over its low-resolution counterpart, just as in Section 3. Indeed, the absence of the gradient correction in the low-resolution ODE leads to the lack of the term (k+1)s∇f(xk)(k+1)s\nabla f(x_{k}) in the Lyapunov function (see Section 4 of [SBC16]), as opposed to the high-resolution Lyapunov function (4.6). Accordingly, it is unlikely to carry over the bound E(k+1)−E(k)≤−O(s2k2∥∇f(xk+1)∥2)\mathcal{E}(k+1)-\mathcal{E}(k)\leq-O(s^{2}k^{2}\|\nabla f(x_{k+1})\|^{2}) of Lemma 4.3 to the low-resolution case and, consequently, the low-resolution ODE approach pioneered by [SBC16] is insufficient to obtain the O(L2/k3)O(L^{2}/k^{3}) rate for squared gradient norm minimization.

3 A Modified NAG-C without a Phase-Space Representation

This section proposes a new accelerated method that also achieves the O(L2/k3)O(L^{2}/k^{3}) rate for minimizing the squared gradient norm. This method takes the following form:

starting with x0x_{0} and y0=x0y_{0}=x_{0}. As shown by the following theorem, this new method has the same convergence rates as NAG-C.

We refer readers to Appendix C.2 for the proof of Theorem 7, which is, as earlier, based on a Lyapunov function. However, since both f(xk)f(x_{k}) and f(yk)f(y_{k}) appear in the iteration, (4.15) does not admit a phase-space representation. As a consequence, the construction of the Lyapunov function is complex; we arrived at it via trial and error. Our initial aim was to seek possible improved rates of the original NAG-C without using the phase-space representation, but the enormous challenges arising in this process motivated us to (1) modify NAG-C to the current (4.15), and (2) to adopt the phase-space representation. Employing the phase-space representation yields a simple proof of the O(L2/k3)O(L^{2}/k^{3}) rate for the original NAG-C and this technique turned out to be useful for other accelerated methods.

Extensions

Motivated by the high-resolution ODE (1.12) of NAG-C, this section considers a family of generalized high-resolution ODEs that take the form

for t≥αs/2t\geq\alpha\sqrt{s}/2, with initial conditions X(αs/2)=x0X(\alpha\sqrt{s}/2)=x_{0} and X˙(αs/2)=−s∇f(x0)\dot{X}(\alpha\sqrt{s}/2)=-\sqrt{s}\nabla f(x_{0}). As demonstrated in [SBC16, ACR17, VJFC18], the low-resolution counterpart (that is, set s=0s=0) of (5.1) achieves acceleration if and only if α≥3\alpha\geq 3. Accordingly, we focus on the case where the friction parameter α≥3\alpha\geq 3 and the gradient correction parameter β>0\beta>0. An investigation of the case of α<3\alpha<3 is left for future work.

By discretizing the ODE (5.1), we obtain a family of new accelerated methods for minimizing smooth convex functions:

starting with x0=y0x_{0}=y_{0}. The second line of the iteration is equivalent to

In Section 5.1, we study the convergence rates of this family of generalized NAC-C algorithms along the lines of Section 4. To further our understanding of (5.2), Section 5.2 shows that this method in the super-critical regime (that is, α>3\alpha>3) converges to the optimum actually faster than O(1/(sk2))O(1/(sk^{2})). As earlier, the proofs of all the results follow the high-resolution ODE framework introduced in Section 2. Proofs are deferred to Appendix D. Finally, we note that Section 6 briefly sketches the extensions along this direction for NAG-SC.

The theorem below characterizes the convergence rates of the generalized NAG-C (5.2).

for all k≥0k\geq 0. The constants cα,βc_{\alpha,\beta} and Cα,βC_{\alpha,\beta} only depend on α\alpha and β\beta.

The proof of Theorem 8 is given in Appendix D.1 for α=3\alpha=3 and Appendix D.2 for α>3\alpha>3. This theorem shows that the generalized NAG-C achieves the same rates as the original NAG-C in both squared gradient norm and function value minimization. The constraint β>12\beta>\frac{1}{2} reveals that further leveraging of the gradient correction does not hurt acceleration, but perhaps not the other way around (note that NAG-C in its original form corresponds to β=1\beta=1). It is an open question whether this constraint is a technical artifact or is fundamental to acceleration.

2 Faster Convergence in Super-Critical Regime

We turn to the case in which α>3\alpha>3, where we show that the generalized NAG-C in this regime attains a faster rate for minimizing the function value. The following proposition provides a technical inequality that motivates the derivation of the improved rate.

where the constants cα,β′c^{\prime}_{\alpha,\beta} and Cα,β′C^{\prime}_{\alpha,\beta} only depend on α\alpha and β\beta.

In relating to Theorem 8, one can show that Proposition 5.1 in fact implies (5.3) in Theorem 8. To see this, note that for k≥1k\geq 1, one has

where the second inequality follows from Proposition 5.1.

Proposition 5.1 can be thought of as a generalization of Theorem 6 of [SBC16]. In particular, this result implies an intriguing and important message. To see this, first note that, by taking s=O(1/L)s=O(1/L), Proposition 5.1 gives

which would not be valid if f(xk)−f(x⋆)≥cL∥x0−x⋆∥2/k2f(x_{k})-f(x^{\star})\geq cL\left\|x_{0}-x^{\star}\right\|^{2}/k^{2} for a constant c>0c>0. Thus, it is tempting to suggest that there might exist a faster convergence rate in the sense that

This faster rate is indeed achievable as we show next, though there are examples where (5.4) and f(xk)−f(x⋆)=O(L∥x0−x⋆∥2/k2)f(x_{k})-f(x^{\star})=O(L\left\|x_{0}-x^{\star}\right\|^{2}/k^{2}) are both satisfied but (5.5) does not hold (a counterexample is given in Appendx D.3).

Under the same assumptions as in Proposition 5.1, taking the step size s=cα,β′/Ls=c^{\prime}_{\alpha,\beta}/L, the iterates {xk}k=0∞\{x_{k}\}_{k=0}^{\infty} generated by the generalized NAG-C (5.2) starting from any x0≠x⋆x_{0}\neq x^{\star} satisfy

Figures 5 and 6 present several numerical studies concerning the prediction of Theorem 9. For a fixed dimension nn, the convergence in Theorem 9 is uniform over functions in F1=∪L>0FL1\mathcal{F}^{1}=\cup_{L>0}\mathcal{F}_{L}^{1} and, consequently, is independent of the Lipschitz constant LL and the initial point x0x_{0}. In addition to following the high-resolution ODE framework, the proof of this theorem reposes on the finiteness of the series in Proposition 5.1. See Appendix D.2 and Appendix D.4 for the full proofs of the proposition and the theorem, respectively.

In the literature, [AP16, May17, ACPR18] use low-resolution ODEs to establish the faster rate o(1/k2)o(1/k^{2}) for the generalized NAG-C (5.2) in the special case of β=1\beta=1. In contrast, our proof of Theorem 9 is more general and applies to a broader class of methods.

In passing, we make the observation that Proposition 5.1 reveals that

which would not hold if min⁡0≤i≤k∥∇f(xi)∥2≥c∥x0−x⋆∥2/(s2k3)\min_{0\leq i\leq k}\|\nabla f(x_{i})\|^{2}\geq c\|x_{0}-x^{\star}\|^{2}/(s^{2}k^{3}) for all kk and a constant c>0c>0. In view of the above, it might be true that the rate of the generalized NAG-C for minimizing the squared gradient norm can be improved to

We leave the confirmation or disconfirmation of this asymptotic result for future research.

Discussion

In this paper, we have proposed high-resolution ODEs for modeling three first-order optimization methods—the heavy-ball method, NAG-SC, and NAG-C. These new ODEs are more faithful surrogates for the corresponding discrete optimization methods than existing ODEs in the literature, thus serving as a more effective tool for understanding, analyzing, and generalizing first-order methods. Using this tool, we identified a term that we refer to as “gradient correction” in NAG-SC and in its high-resolution ODE, and we demonstrate its critical effect in making NAG-SC an accelerated method, as compared to the heavy-ball method. We also showed via the high-resolution ODE of NAG-C that this method minimizes the squared norm of the gradient at a faster rate than expected for smooth convex functions, and again the gradient correction is the key to this rate. Finally, the analysis of this tool suggested a new family of accelerated methods with the same optimal convergence rates as NAG-C.

The aforementioned results are obtained using the high-resolution ODEs in conjunction with a new framework for translating findings concerning the amenable ODEs into those of the less “user-friendly” discrete methods. This framework encodes an optimization property under investigation to a continuous-time Lyapunov function for an ODE and a discrete-time Lyapunov function for the discrete method. As an appealing feature of this framework, the transformation from the continuous Lyapunov function to its discrete version is through a phase-space representation. This representation links continuous objects such as position and velocity variables to their discrete counterparts in a faithful manner, permitting a transparent analysis of the three discrete methods that we studied.

There are a number of avenues open for future research using the high-resolution ODE framework. First, the discussion of Section 5 can carry over to the heavy-ball method and NAG-SC, which correspond to the high-resolution ODE

with β=0\beta=0 and β=1\beta=1, respectively. This ODE with a general 0<β<10<\beta<1 corresponds to a new algorithm that can be thought of as an interpolation between the two methods. It is of interest to investigate the convergence properties of this class of algorithms. Second, we recognize that new optimization algorithms are obtained in [WWJ16, WRJ16] by using different discretization schemes on low-resolution ODE. Hence, a direction of interest is to apply the techniques therein to our high-resolution ODEs and to explore possible appealing properties of the new methods. Third, the technique of dimensional analysis, which we have used to derive high-resolution ODEs, can be further used to incorporate even higher-order powers of s\sqrt{s} into the ODEs. This might lead to further fine-grained findings concerning the discrete methods.

More broadly, we wish to remark on possible extensions of the high-resolution ODE framework beyond smooth convex optimization in the Euclidean setting. In the non-Euclidean case, it would be interesting to derive a high-resolution ODE for mirror descent [KBB15, WWJ16]. This framework might also admit extensions to non-smooth optimization and stochastic optimization, where the ODEs are replaced, respectively, by differential inclusions [ORX+16, VJFC18] and stochastic differential equations [KB17, HLLL17, LTE17, LS17, XWG18, HMC+18, GGZ18]. Finally, recognizing that the high-resolution ODEs are well-defined for non-convex functions, we believe that this framework will provide more accurate characterization of local behaviors of first-order algorithms near saddle points [JGN+17, DJL+17, HLS17]. On a related note, given the centrality of the problem of finding an approximate stationary point in the non-convex setting [CDHS17a, CDHS17b, AZ18], it is worth using the high-resolution ODE framework to explore possible applications of the faster rate for minimizing the squared gradient norm that we have uncovered.

B. S. is indebted to Xiaoping Yuan for teaching him the modern theory of ordinary differential equations and would like to thank Rui Xin Huang for teaching him how to leverage intuitions from physics to understand differential equations. We would like to thank Nicolas Flammarion for suggesting references. This work was supported in part by the NSF via grant CCF-1763314 and Army Research Office via grant W911NF-17-1-0304.

References

Appendix A Technical Details in Section 2

By only ignoring the O(s)O(s) term, we obtain the high-resolution ODE (1.10) for the heavy-ball method

For convenience, we slightly change the definition tk=ks+(3/2)st_{k}=k\sqrt{s}+(3/2)\sqrt{s} instead of tk=kst_{k}=k\sqrt{s}. Plugging (2.1) into (A.2), we have

Ignoring any O(s)O(s) terms, we obtain the high-resolution ODE (1.12) for NAG-C

A.2 Derivation of Low-Resolution ODEs

In this section, we derive low-resolution ODEs of accelerated gradient methods for comparison. The results presented here are well-known in the literature and the purpose is for ease of reading. In [SBC16], the second-order Taylor expansions at both xk−1x_{k-1} and xk+1x_{k+1} with the step size s\sqrt{s} are,

With the Taylor expansion (A.3), we obtain the gradient correction

From (A.3) and (A.4), we can derive the following low-resolution ODEs.

Recall the equivalent form (2.3) of NAG-SC (1.3) is

Plugging (A.3) and (A.4) into (2.3), we have

Hence, taking s→0s\rightarrow 0, we obtain the low-resolution ODE (1.9) of NAG-SC

Recall the equivalent form (A.1) of the heavy-ball method (1.2) is

Plugging (A.3) and (A.4) into (A.1), we have

Hence, taking s→0s\rightarrow 0, we obtain the low-resolution ODE (1.9) of the heavy-ball method

Notably, NAG-SC and the heavy-ball method share the same low-resolution ODE (1.9), which is almost consistent with (1.10). Thus the low-resolution ODE fails to capture the information from the “gradient correction” of NAG-SC.

Plugging (A.3) and (A.4) into (A.2), we have

Thus, by taking s→0s\rightarrow 0, we obtain the low-resolution ODE (1.8) of NAG-C

A.3 Solution Approximating Optimization Algorithms

To investigate the property about the high-resolution ODEs (1.10), (1.11) and (1.12), we need to state the relationship between them and their low-resolution corresponding ODEs. Here, we denote the solution to high-order ODE by Xs=Xs(t)X_{s}=X_{s}(t). Actually, the low-resolution ODE is the special case of high-resolution ODE with s=0s=0. Take NAG-SC for example

In other words, we consider a family of ODEs about the step size parameter ss.

To prove the global existence and uniqueness of solution to the high-resolution ODEs (1.10) and (1.11), we first emphasize a fact that if Xs=Xs(t)X_{s}=X_{s}(t) is the solution of (1.10) or (1.11), there exists some constant C1>0\mathcal{C}_{1}>0 such that

which is only according to the following Lyapunov function

of which the classical theory about global existence and uniqueness of solution is shown as below.

For the heavy-ball method, the phase-space representation of high-resolution ODE (1.10) is

For any (Xs,X˙s)⊤,(Ys,Y˙s)⊤∈MC1(X_{s},\dot{X}_{s})^{\top},(Y_{s},\dot{Y}_{s})^{\top}\in M_{\mathcal{C}_{1}}, we have

For NAG-SC, the phase-space representation of high-resolution ODE (1.11) is

For any (Xs,X˙s)⊤,(Ys,Y˙s)⊤∈MC1(X_{s},\dot{X}_{s})^{\top},(Y_{s},\dot{Y}_{s})^{\top}\in M_{\mathcal{C}_{1}}, we have

Based on the phase-space representation (A.8) and (A.10), together with the Lipschitz condition (• ‣ A.3.1) and (• ‣ A.3.1), Theorem 10 leads to the following Corollary.

Based on the Lyapunov function (A.6), the gradient norm is bounded along the solution of (1.10) or (1.11), that is,

Recall the low-resolution ODE (1.9), the phase-space representation is proposed as

Similarly, using a Lyapunov function argument, we can show that if X=X(t)X=X(t) is a solution of (1.9), we have

Simple calculation tells us that there exists some constant L1>0\mathcal{L}_{1}>0 such that

Now, we proceed to show the approximation.

Let the solution to high-resolution ODEs (1.10) and (1.11) as X=Xs(t)X=X_{s}(t) and that of (1.9) as X=X(t)X=X(t), then we have

In order to prove (A.16), we prove a stronger result as

Before we start to prove (A.17), we first describe the standard Gronwall-inequality as below.

Let m(t)m(t), t∈[0,T]t\in[0,T], be a nonnegative function satisfying the relation

The proof is only according to simple calculus, here we omit it.

For the heavy-ball method, the phase-space representations (A.8) and (A.13) tell us that

By the boundedness (A.12), (A.5) and (A.14) and the inequality (A.15), we have

For NAG-SC, the phase-space representations (A.10) and (A.13) tell us that

Similarly, by the boundedness (A.12), (A.5) and (A.14) and the inequality (A.15), we have

The two methods, heavy-ball method and NAG-SC, converge to their low-resolution ODE (1.9) in the sense that

This result has bee studied in [WRJ16] and the method for proof refer to [SBC16, Appendix 2]. Combined with Corollary A.1, Lemma A.2 and Lemma A.4, we complete the proof of Proposition 2.1.

A.3.2 Proof of Proposition 2.2

Similar as Appendix A.3.1, we first emphasize the fact that if Xs=Xs(t)X_{s}=X_{s}(t) is the solution of high-resolution ODE (1.12), there exists some constant C6\mathcal{C}_{6} such that

which is only according to the following Lyapunov function

of which the classical theory about global existence and uniqueness of solution is shown as below.

for all (x,t),(y,t)∈M×I(x,t),(y,t)\in M\times I. Then for any x0∈Mx_{0}\in M, the IVP (A.20) has a unique solution x(t)x(t) defined for all t∈It\in I.

The proof is consistent with Theorem 33 and Theorem 44 of Chapter 3.13.1 in [Per13] except the Lipschitz condition for the vector field

for any x,y∈Mx,y\in M. The readers can also refer to [GH13]. Similarly, the set

is a compact manifold satisfying Theorem 11 with m=2nm=2n.

For NAG-C, the phase-space representation of high-resolution ODE (1.11) is

For any (Xs,X˙s,t),(Ys,Y˙s,t)∈MC6×[(3/2)s,∞)(X_{s},\dot{X}_{s},t),(Y_{s},\dot{Y}_{s},t)\in M_{\mathcal{C}_{6}}\times\left[(3/2)\sqrt{s},\infty\right), we have

Based on the phase-space representation (A.21), together with (A.3.2), Theorem 11 leads the following Corollary.

Using a linear transformation t+(3/2)st+(3/2)\sqrt{s} instead of tt, we can rewrite high-resolution ODE (1.12) as

for t≥0t\geq 0, with initial Xs(0)=x0X_{s}(0)=x_{0} and X˙s(0)=−s∇f(x0)\dot{X}_{s}(0)=-\sqrt{s}\nabla f(x_{0}), of which the phase-space representation is

Here, we adopt the technique max⁡{δ,t}\max\{\delta,t\} instead of tt for any δ>0\delta>0 to overcome the singular point t=0t=0, which is used firstly in [SBC16]. Then (A.24) is replaced into

with the initial Xsδ(0)=x0X_{s}^{\delta}(0)=x_{0} and X˙sδ(0)=−s∇f(x0)\dot{X}_{s}^{\delta}(0)=-\sqrt{s}\nabla f(x_{0}). Recall the low-resolution ODE (1.8), with the above technique, the phase-space representation is proposed as

with the initial Xsδ(0)=x0X_{s}^{\delta}(0)=x_{0} and X˙sδ(0)=0\dot{X}_{s}^{\delta}(0)=0. Then according to (A.25) and (A.26), if we can prove for any δ>0\delta>0 and any t∈[0,T]t\in[0,T], the following equality holds

Then, we can obtain the desired result as

Similarly, using Lyapunov function argument, we can show that the solutions XsδX_{s}^{\delta} and XδX^{\delta} satisfy

for all t≥0t\geq 0. Now, we proceed to show the approximation.

Denote the solution to high-resolution ODE (1.12) as X=Xs(t)X=X_{s}(t) and that to (1.8) as X=X(t)X=X(t). We have

In order to prove (A.30), we prove a stronger result

The phase-space representation (A.25) and (A.26) tell us that

By the boundedness (A.27) and (A.28) and the Lipschitz inequality (A.3.2), we have

According to Lemma A.3, we obtain the result as (A.31)

NAG-C converges to its low-resolution ODE in the sense that

Combined with Corollary A.5, Lemma A.6 and Lemma A.7, we complete the proof of Proposition 2.2.

A.4 Closed-Form Solutions for Quadratic Functions

In this section, we propose the closed-form solutions to the three high-resolution ODEs for the quadratic objective function

The closed-form solution of (A.33) can be shown from the theory of ODE, as below.

When θ>μ\theta>\mu, that is, 4μ−4θ<04\mu-4\theta<0, the closed-form solution is the superimposition of two independent oscillation solutions

When θ=μ\theta=\mu, that is, 4μ−4θ=04\mu-4\theta=0, the closed-form solution is the superimposition of two independent non-oscillation solutions

Second, plugging the quadratic objective (A.32) into the high-resolution ODE (1.11) of NAG-SC, we have

The closed-form solutions to (A.34) are shown as below.

When s<4(θ−μ)θ2s<\frac{4(\theta-\mu)}{\theta^{2}}, that is, 4(μ−θ)+sθ2<04(\mu-\theta)+s\theta^{2}<0, the closed-form solution is the superimposition of two independent oscillation solutions

When s=4(θ−μ)θ2s=\frac{4(\theta-\mu)}{\theta^{2}}, that is, 4(μ−θ)+sθ2=04(\mu-\theta)+s\theta^{2}=0, the closed-form solution is the superimposition of two independent non-oscillation solutions

When s>4(θ−μ)θ2s>\frac{4(\theta-\mu)}{\theta^{2}}, that is, 4(μ−θ)+sθ2>04(\mu-\theta)+s\theta^{2}>0, the closed-form solution is also the superimposition of two independent non-oscillation solutions

Hence, when the step size satisfies s≥2s\geq 2, there is always no oscillation in the closed-form solution of (A.34).

Finally, plugging the quadratic objective (A.32) into the high-resolution ODE (1.10) of the heavy-ball method, we have

Since 4μ−4(1+μs)θ<04\mu-4(1+\sqrt{\mu s})\theta<0 is well established, the closed-form solution of (A.35) is the superimposition of two independent oscillation solutions

In summary, both the closed-form solutions to (A.33) and (A.35) are oscillated except the fragile condition θ=μ\theta=\mu and the speed of linear convergence is Θ(e−μt)\Theta\left(e^{-\sqrt{\mu}t}\right). However, the rate of convergence in the closed-form solution to the high-resolution ODE (A.34) is always faster than Θ(e−μt)\Theta\left(e^{-\sqrt{\mu}t}\right). Additionally, when the step size s≥2s\geq 2, there is always no oscillation in the closed-form solution of the high-resolution ODE (A.34).

A.4.2 Kummer’s Equation and Confluent Hypergeometric Function

the closed-form solution of which has been proposed in [SBC16]

where J1(⋅)J_{1}(\cdot) and Y1(⋅)Y_{1}(\cdot) are the Bessel function of the first kind and the second kind, respectively. According to the asymptotic property of Bessel functions,

Now, we plug the quadratic objective (A.32) into the high-resolution ODE (1.12) of NAG-C and obtain

For convenience, we define two new parameters as

which actually corresponds to the Kummer’s equation. According to the closed-form solution to Kummer’s equation, the high-resolution ODE (A.36) for quadratic function can be solved analytically as

where M(⋅,⋅,⋅)M(\cdot,\cdot,\cdot) and U(⋅,⋅,⋅)U(\cdot,\cdot,\cdot) are the confluent hypergeometric functions of the first kind and the second kind. The integral expressions of M(⋅,⋅,⋅)M(\cdot,\cdot,\cdot) and U(⋅,⋅,⋅)U(\cdot,\cdot,\cdot) are given as

Since the possible value of arg⁡(ξt)\arg(\xi t) either or π/2\pi/2, we have

Apparently, from the asymptotic estimate of (A.38), we have

When s<4/θs<4/\theta, that is, sθ2−4θ<0s\theta^{2}-4\theta<0, the closed-form solution (A.37) is estimated as

Hence, when the step size satisfies s<4/Ls<4/L, the above upper bound always holds.

When s≥4/θs\geq 4/\theta, that is, sθ2−4θ≥0s\theta^{2}-4\theta\geq 0, the closed-form solution (A.37) is estimated as

Appendix B Technical Details in Section 3

the Lyapunov function (3.3) can be estimated as

Along the solution to the high-resolution ODE (1.10), the time derivative of the Lyapunov function (3.3) is

the time derivative of the Lyapunov function can be estimated as

B.2 Completing the Proof of Lemma 3.4

when the step size satisfies s≤1/(2L)≤1/Ls\leq 1/(2L)\leq 1/L, we have

B.2.2 Derivation of (B.2.1)

Now, we show the derivation of (B.2.1). Recall the discrete Lyapunov function (2.6),

For convenience, we calculate the difference between E(k)\mathcal{E}(k) and E(k+1)\mathcal{E}(k+1) by the three parts, I\mathbf{I}, II\mathbf{II} and III\mathbf{III} respectively.

For the part I\mathbf{I}, potential, with the convexity, we have

For the part II\mathbf{II}, kinetic energy, with the phase representation of NAG-SC (2.5), we have

For the part III\mathbf{III}, mixed energy, with the phase representation of NAG-SC (2.5), we have

Both II2\mathbf{II}_{2} and III3\mathbf{III}_{3} above are the discrete correspondence of the terms −s2∥∇f(X(t))∥2-\frac{\sqrt{s}}{2}\left\|\nabla f(X(t))\right\|^{2} and −s2X˙(t)⊤∇2f(X(t))X˙(t)-\frac{\sqrt{s}}{2}\dot{X}(t)^{\top}\nabla^{2}f(X(t))\dot{X}(t) in (3.1). The impact can be found in the calculation. Now, we calculate the difference of discrete Lyapunov function (2.6) at kk-th iteration by the simple operation

Now, the term, (1/2)I1+II5+II6+III4(1/2)\mathbf{I}_{1}+\mathbf{II}_{5}+\mathbf{II}_{6}+\mathbf{III}_{4}, can be calculated as

With phase representation of NAG-SC (2.5), we have

For convenience, we note the term IV=(1/2)I1+III2+III3\mathbf{IV}=(1/2)\mathbf{I}_{1}+\mathbf{III}_{2}+\mathbf{III}_{3}. Then, with phase representation of NAG-SC (2.5), the difference of Lyapunov function (2.6) is

Now, we can find the impact of additional term in the Lyapunov function (2.6). In other words, the II4+IV1\mathbf{II}_{4}+\mathbf{IV}_{1} term added the additional term is a perfect square, as below

Merging all the similar items, II4+IV1+additional  term\mathbf{II}_{4}+\mathbf{IV}_{1}+\mathbf{additional\;term}, I2+II3\mathbf{I}_{2}+\mathbf{II}_{3}, we have

Now, we obtain that the difference of Lyapunov function (2.6) is

B.3 Proof of Lemma 3.5

With the phase representation of the heavy-ball method (3.11) and Cauchy-Schwarz inequality, we have

The discrete Lyapunov function (3.10) can be estimated as

For convenience, we also split the discrete Lyapunov function (3.10) into three parts and mark them as below

where the three parts I\mathbf{I}, II\mathbf{II} and III\mathbf{III} are corresponding to potential, kinetic energy and mixed energy in classical mechanics, respectively.

For the part II\mathbf{II}, kinetic energy, with the phase representation of the heavy-ball method (3.11), we have

For the part III\mathbf{III}, mixed energy, with the phase representation of the heavy-ball method (3.11), we have

Now, we calculate the difference of discrete Lyapunov function (2.6) at the kk-th iteration by the simple operation as

With the phase representation of the heavy-ball method (3.11), we have

Now, the difference of discrete Lyapunov function (3.10) can be rewritten as

Comparing the coefficient of the estimate of Lyapunov function (B.3), we have

Appendix C Technical Details in Section 4

When k=2k=2, the iterate (xk,yk)(x_{k},y_{k}) is

When k=3k=3, the iterate (xk,yk)(x_{k},y_{k}) is

Taking s≤1/(3L)s\leq 1/(3L) and using (C.4), (C.5) and (C.6), we have

Taking s≤1/(3L)s\leq 1/(3L), (C.7) tells us that

C.1.4 Estimate for Lyapunov function ℰ​(2)ℰ2\mathcal{E}(2) and ℰ​(3)ℰ3\mathcal{E}(3)

With the phase-space representation form (4.5), we have

According to (4.6), the Lyapunov function E(2)\mathcal{E}(2) can be written as

With (C.8) and Cauchy-Schwarz inequality, we have

Finally, with (C.4)-(C.5), Cauchy-Schwarz inequality tells

By Lemma 4.3, when the step size s≤1/(3L)s\leq 1/(3L), (C.9) tells us

C.2 Proof of Theorem 7

Let wk=(1/2)[(k+2)xk−kyk+(k−1)s∇f(yk)]w_{k}=(1/2)\left[(k+2)x_{k}-ky_{k}+(k-1)s\nabla f(y_{k})\right] for convenience. Using the dynamics of {(xk,yk)}k=0∞\{(x_{k},y_{k})\}_{k=0}^{\infty} generated by the modified NAG-C (4.15), we have

Hence, the difference between ∥wk+1−x⋆∥2\left\|w_{k+1}-x^{\star}\right\|^{2} and ∥wk−x⋆∥2\left\|w_{k}-x^{\star}\right\|^{2} is

Hence, the difference between E(k+1)\mathcal{E}(k+1) and E(k)\mathcal{E}(k) in (C.11) is

Similarly, when s≤1/Ls\leq 1/L, for k=0k=0, we have

for k=1k=1, following the modified NAG-C (4.15), we obtain (x1,y1)(x_{1},y_{1}) as

C.3 Nesterov’s Lower Bound

Appendix D Technical Details in Section 5

Before starting to prove Theorem 8, we first look back our high-resolution ODE framework in Section 2.

Step 11, the generalized high-resolution ODE has been given in (5.1).

Step 22, the continuous Lyapunov function is constructed as

Following this Lyapunov function (D.1), we can definitely obtain similar results as Theorem 5 and Corollary 4.2. The detailed calculation, about the estimate of the optimal constant β\beta and how the constant β\beta influence the initial point, is left for readers.

Step 33, before constructing discrete Lyapunov functions, we show the phase-space representation (5.2) as

Now, we show how to construct the discrete Lyapunov function and analyze the algorithms (5.2) with α=3\alpha=3 in order to prove Theorem 8.

When β<1\beta<1, we know that the function

decreases monotonically. Hence we can construct the discrete Lyapunov function as

which is slightly different from the discrete Lyapunov function (4.6) for NAG-C. When β→1\beta\rightarrow 1, the discrete Lyapunov function (D.3) approximate to (4.6) as k→∞k\rightarrow\infty.

With the phase-space representation (D.2) for α=3\alpha=3, we can obtain

The difference of the discrete Lyapunov function (D.3) of the kk-th iteration is

the difference of the discrete Lyapunov function (D.3) can be estimated as

Utilizing the phase-space representation (D.2) again, we calculate the difference of the discrete Lyapunov function (D.3) as

To guarantee that the Lyapunov function E(k)\mathcal{E}(k) is decreasing, a sufficient condition is

Simple calculation tells us that (D.5) can be rewritten as

Apparently, when β→1\beta\rightarrow 1, the step size satisfies

which is consistent with (4.8). Now, we turn to discuss the parameter 0≤β<10\leq\beta<1 case by case.

When the parameter β≤1/2\beta\leq 1/2, the sufficient condition (D.5) for the Lyapunov function E(k)\mathcal{E}(k) decreasing cannot be satisfied for sufficiently large kk.

When the parameter 1/2<β<11/2<\beta<1, since the function h(k)=1Lβ2(2β−1+β−3k+1)h(k)=\frac{1}{L\beta^{2}}\left(2\beta-1+\frac{\beta-3}{k+1}\right) increases monotonically for k≥0k\geq 0, there exists k3,β=⌊4−3β2β−1⌋+1k_{3,\beta}=\left\lfloor\frac{4-3\beta}{2\beta-1}\right\rfloor+1 such that the step size

works for any k≥k3,βk\geq k_{3,\beta} (k3,β→2k_{3,\beta}\rightarrow 2 with β→1\beta\rightarrow 1). Then, the difference of the discrete Lyapunov function (D.3) can be estimated as

Here, the proof is actually complete. Without loss of generality, we briefly show the expression is consistent with Theorem 8 and omit the proofs for the following facts. When k≥k3,β+1k\geq k_{3,\beta}+1, there exists some constant C3,β0>0\mathfrak{C}^{0}_{3,\beta}>0 such that

For k≤k3,βk\leq k_{3,\beta}, using mathematic induction, there also exists some constant C3,β1>0\mathfrak{C}^{1}_{3,\beta}>0 such that for s=O(1/L)s=O(1/L), we have

D.1.2 Case: β≥1𝛽1\beta\geq 1

When β≥1\beta\geq 1, we know that the function

decreases monotonically. Hence we can construct the discrete Lyapunov function as

which for β=1\beta=1 is consistent with the discrete Lyapunov function (4.6) for NAG-C.

the difference of the discrete Lyapunov function (D.7) of the kk-th iteration is

the difference of the discrete Lyapunov function (D.7) can be estimated as

Utilize the phase-space representation (D.2) again, we calculate the difference of the discrete Lyapunov function (D.7) as

Consistently, we can obtain the sufficient condition for the Lyapunov function E(k)\mathcal{E}(k) decreasing (D.5) and the sufficient condition for step size (D.6).

Now, we turn to discuss the parameter β≥1\beta\geq 1 case by case.

When the parameter β≥3\beta\geq 3, since the function h(k)=1Lβ2(2β−1+β−3k+1)h(k)=\frac{1}{L\beta^{2}}\left(2\beta-1+\frac{\beta-3}{k+1}\right) decreases monotonically for k≥0k\geq 0, then the condition of the step size

holds for (D.5), where ϵ>0\epsilon>0 is a real number. Hence, when k≥k3,β+1k\geq k_{3,\beta}+1, where

the difference of the discrete Lyapunov function (D.7) can be estimated as

When the parameter 1≤β<31\leq\beta<3, since the function h(k)=1Lβ2(2β−1+β−3k+1)h(k)=\frac{1}{L\beta^{2}}\left(2\beta-1+\frac{\beta-3}{k+1}\right) increases monotonically for k≥0k\geq 0, there exists k3,β=max⁡{0,⌊β−3⌋+1,⌊4−3β2β−1⌋+1}k_{3,\beta}=\max\left\{0,\left\lfloor\beta-3\right\rfloor+1,\left\lfloor\frac{4-3\beta}{2\beta-1}\right\rfloor+1\right\} such that the step size

works for any k≥k3,βk\geq k_{3,\beta}. When β=1\beta=1, the step size satisfies

which is consistent with (4.8) and k3,β=2k_{3,\beta}=2. Then, the difference of the discrete Lyapunov function (D.3) can be estimated as

By simple calculation, we complete the proof.

D.2 Proof of Theorem 8: Case α>3𝛼3\alpha>3

Before starting to prove Theorem 8: Case α>3\alpha>3, we first also look back our high-resolution ODE framework in Section 2.

Step 11, the generalized high-resolution ODE has been given in (5.1).

Step 22, the continuous Lyapunov function is constructed as

which is consistent with (D.1) for α→3\alpha\rightarrow 3. Following this Lyapunov function (D.8), we can obtain

for any t>t0=max⁡{s(α/2−β)(α−2)/(α−3),s(α/2)}t>t_{0}=\max\left\{\sqrt{s}(\alpha/2-\beta)(\alpha-2)/(\alpha-3),\sqrt{s}(\alpha/2)\right\}. The two inequalities of (D.9) for the convergence rate of function value is stronger than Corollary 4.2. The detailed calculation, about the estimate of the optimal constant β\beta and how the constant β\beta influences the initial point, is left for readers.

Step 33, before constructing discrete Lyapunov functions, we look back the phase-space representation (D.2)

The discrete functional is constructed as

When β=1\beta=1, with α→3\alpha\rightarrow 3, the discrete Lyapunov function E(k)\mathcal{E}(k) degenerates to (4.6).

Now, we procced to Step 44 to analyze the algorithms (5.2) with α>3\alpha>3 in order to prove Theorem 5.1. The simple transformation of (D.2) for α>3\alpha>3 is

Thus, the difference of the Lyapunov function (D.10) on the kk-th iteration is

the difference of the discrete Lyapunov function (D.10) can be estimated as

Utilizing the phase-space representation (D.2) again, we calculate the difference of the discrete Lyapunov function (D.10) as

To guarantee the Lyapunov function E(k)\mathcal{E}(k) decreasing, a sufficient condition is

With the inequality (D.12), the step size can be estimated as

When the parameter β>1/2\beta>1/2 and α<β\alpha<\beta, since the function h(k)=2β−1Lβ2−α−β(k+1)Lβ2h(k)=\frac{2\beta-1}{L\beta^{2}}-\frac{\alpha-\beta}{(k+1)L\beta^{2}} decreases monotonically for k≥0k\geq 0, thus the step size

holds for (D.12), where ϵ>0\epsilon>0 is a real number. Hence, when k≥kα,β+1k\geq k_{\alpha,\beta}+1, where

the difference of the discrete Lyapunov function (D.10) can be estimated as

When the parameter β>1/2\beta>1/2 and α≥β\alpha\geq\beta, since the function h(k)=2β−1Lβ2−α−β(k+1)Lβ2h(k)=\frac{2\beta-1}{L\beta^{2}}-\frac{\alpha-\beta}{(k+1)L\beta^{2}} increases monotonically for k≥0k\geq 0, there exists

which is consistent with (4.8). Then, the difference of the discrete Lyapunov function (D.10) can be estimated as

D.3 A Simple Counterexample

The simple counterexample is constructed as

Hence, Proposition 5.1 cannot guarantee the faster convergence rate.

Here, we still turn back to our high-resolution ODE framework in Section 2. The generalized high-resolution ODE has been still shown in (5.1). A more general Lyapunov function is constructed as

where 2<ν≤α−12<\nu\leq\alpha-1. When ν=α−1\nu=\alpha-1, the Lyapunov function (D.13) degenerates to (D.8). Furthermore, when ν=α−1→2\nu=\alpha-1\rightarrow 2, the Lyapunov function (D.13) degenerates to (D.1). Finally, when 2=ν=α−12=\nu=\alpha-1 and β=1\beta=1, the Lyapunov function (D.13) is consistent with (4.1). We assume that initial time is

Based on the Lyapunov function (D.13), we have the following results.

for all t≥tα,β,νt\geq t_{\alpha,\beta,\nu}, where the positive constant Cα,β,ν2\mathfrak{C}^{2}_{\alpha,\beta,\nu} and the integer tα,β,νt_{\alpha,\beta,\nu} depend only on α\alpha, β\beta and ν\nu. In other words, the equivalent expression of (D.14) is

Now, we start to show the proof. Since X=X(t)X=X(t) is the solution of the ODE (5.1) with α>3\alpha>3 and β>0\beta>0, when t>tα,β,νt>t_{\alpha,\beta,\nu}, the time derivative of Lyapunov function (D.13) is

the time derivative of Lyapunov function (D.4.1) can be estimated as

With the Lyapunov function Eν(t)≥0\mathcal{E}_{\nu}(t)\geq 0 and the technique for integral, for any t>t0t>t_{0} we have

where δ<t−t0\delta<t-t_{0}. Thus, we can obtain the following Lemma.

Under the same assumption of Theorem 12, the following limits exist

With (D.4.1) and Lemma D.1, the following Lemma holds.

Under the same assumption of Theorem 12, the following limit exists

Under the same assumption of Theorem 12, the following limits exist

Taking ν≠ν′∈[2,γ−1]\nu\neq\nu^{\prime}\in[2,\gamma-1], we have

With Lemma D.1 and (D.9), the following limit exists

Define a new function about time variable tt:

If we can prove the existence of the limit π(t)\pi(t) with t→∞t\rightarrow\infty, we can guarantee lim⁡t→∞∥X(t)−x⋆∥\lim\limits_{t\rightarrow\infty}\left\|X(t)-x^{\star}\right\| exists with Lemma D.2. We observe the following equality

With (D.16) and Lemma D.2, we obtain that the following limit exists

that is, there exists some constant C3\mathfrak{C}^{3} such that the following equality holds,

For any ϵ>0\epsilon>0, there exists t0>0t_{0}>0 such that when t≥t0t\geq t_{0}, we have

Finally, we finish the proof for Theorem 12.

When t>tα,β,νt>t_{\alpha,\beta,\nu}, we expand the Lyapunov function (D.13) as

With Lemma D.1 and Lemma D.3, we obtain the first equation of (D.14). Furthermore, Cauchy-Scharwz inequality tells that

With Lemma D.1, we obtain the second equation of (D.14). With basic calculation, we complete the proof. ∎

D.4.2 Proof of Theorem 9

Similarly, under the assumption of Theorem 9, if we can show a discrete version of (D.14), that is, there exists some constant Cα,β,ν4>0\mathfrak{C}^{4}_{\alpha,\beta,\nu}>0 and cα,β,ν>0\mathfrak{c}_{\alpha,\beta,\nu}>0 such that when the step size satisfies 0<s≤cα,β,ν/L0<s\leq\mathfrak{c}_{\alpha,\beta,\nu}/L, the following relationship holds

Thus, we obtain the sharper convergence rate as

Now we show the derivation of the inequality (D.17). The discrete Lyapunov function is constructed as

where 2≤ν<α−12\leq\nu<\alpha-1 and parts I\mathbf{I}, II\mathbf{II} and III\mathbf{III} are potential, Euclidean distance and mixed energy respectively. Apparently, when ν=α−1\nu=\alpha-1, the discrete Lyapunov function (D.18) is consistent with (D.10). When β=1\beta=1 and ν=α−1→2\nu=\alpha-1\rightarrow 2, the discrete Lyapunov function (D.18) degenerates to (4.6), Now, we turn to estimate the difference of Lyapunov function (D.18).

For the part I\mathbf{I}, potential, we have

where the last inequality follows k+α+2>k+α+1>k+2k+\alpha+2>k+\alpha+1>k+2.

For the part II\mathbf{II}, Euclidean distance, we have

For the part III\mathbf{III}, mixed energy, with the simple transformation (D.11) for α>3\alpha>3

When k≥n(α−1−ν)β−(α+1−β)k\geq n(\alpha-1-\nu)\beta-(\alpha+1-\beta), we have

With the monotonicity of the following function about kk

we know there exists some constant cα,β,ν\mathfrak{c}_{\alpha,\beta,\nu} and k1,α,β,νk_{1,\alpha,\beta,\nu} such that the step size satisfies 0<s≤cα,β,ν/L0<s\leq\mathfrak{c}_{\alpha,\beta,\nu}/L. When k≥k1,α,β,νk\geq k_{1,\alpha,\beta,\nu}, the following inequality holds

we know that there exists k2,α,β,νk_{2,\alpha,\beta,\nu} such that when k≥k2,α,β,νk\geq k_{2,\alpha,\beta,\nu},

Let kα,β,ν=max⁡{k1,α,β,ν,k2,α,β,ν}+1k_{\alpha,\beta,\nu}=\max\{k_{1,\alpha,\beta,\nu},k_{2,\alpha,\beta,\nu}\}+1. Summing up all the estimates above, when β>1/2\beta>1/2, the difference of discrete Lyapunov function, for any k≥kα,β,νk\geq k_{\alpha,\beta,\nu},

Under the same assumption of Theorem 9, the following limit exists

and the summation of the following series exist

Under the same assumption of Theorem 9, the following limits exist

Taking ν≠ν′∈(2,γ−1]\nu\neq\nu^{\prime}\in(2,\gamma-1], we have

With Lemma D.4, the following limit exists

If we can show the existence of the limit π(k)\pi(k) with k→∞k\rightarrow\infty, we can guarantee lim⁡k→∞∥xk+1−x⋆∥\lim\limits_{k\rightarrow\infty}\left\|x_{k+1}-x^{\star}\right\| exists with Lemma D.4. We observe the following equality

Lemma D.4 and (D.19) tell us there exists some constant C5\mathfrak{C}^{5} such that

that is, taking a simple translation π′(k)=π(k)−C5/(γ−1)\pi^{\prime}(k)=\pi(k)-\mathfrak{C}^{5}/(\gamma-1), we have

Since E(k)\mathcal{E}(k) decreases for k≥kα,β,νk\geq k_{\alpha,\beta,\nu}, thus, ∥xk−x⋆∥2\left\|x_{k}-x^{\star}\right\|^{2} is bounded. With Lemma D.4, we obtain that π(k)\pi(k) is bounded, that is, π′(k)\pi^{\prime}(k) is bounded. Then we have

that is, for any ϵ>0\epsilon>0, there exists k0′>0k^{\prime}_{0}>0 such that

With arbitrary ϵ>0\epsilon>0, we complete the proof of Lemma D.5. ∎

When k≥kα,β,νk\geq k_{\alpha,\beta,\nu}, we expand the discrete Lyapunov function (D.18) as

With Lemma D.4 and Lemma D.5, we obtain the first equation of (D.17). Additionally, we have

With Lemma D.4, we obtain the second equation of (D.17). ∎