A Lyapunov Analysis of Momentum Methods in Optimization

Introduction

Momentum is a powerful heuristic for accelerating the convergence of optimization methods. One can intuitively “add momentum” to a method by adding to the current step a weighted version of the previous step, encouraging the method to move along search directions that had been previously seen to be fruitful. Such methods were first studied formally by Polyak , and have been employed in many practical optimization solvers. As an example, since the 1980s, momentum methods have been popular in neural networks as a way to accelerate the backpropagation algorithm. The conventional intuition is that momentum allows local search to avoid “long ravines” and “sharp curvatures” in the sublevel sets of cost functions .

Polyak motivated momentum methods by an analogy to a “heavy ball” moving in a potential well defined by the cost function. However, Polyak’s physical intuition was difficult to make rigorous mathematically. For quadratic costs, Polyak was able to provide an eigenvalue argument that showed that his Heavy Ball Method required no more iterations than the method of conjugate gradients .Indeed, when applied to positive-definite quadratic cost functions, Polyak’s Heavy Ball Method is equivalent to Chebyshev’s Iterative Method . Despite its intuitive elegance, however, Polyak’s eigenvalue analysis does not apply globally for general convex cost functions. In fact, Lessard et al. derived a simple one-dimensional counterexample where the standard Heavy Ball Method does not converge .

In order to make momentum methods rigorous, a different approach was required. In celebrated work, Nesterov devised a general scheme to accelerate convex optimization methods, achieving optimal running times under oracle models in convex programming . To achieve such general applicability, Nesterov’s proof techniques abandoned the physical intuition of Polyak ; in lieu of differential equations and Lyapunov functions, Nesterov devised the method of estimate sequences to verify the correctness of these momentum-based methods. Researchers have struggled to understand the foundations and scope of the estimate sequence methodology since Nesterov’s initial papers. The associated proof techniques are often viewed as an “algebraic trick.”

To overcome the lack of fundamental understanding of the estimate sequence technique, several authors have recently proposed schemes to achieve acceleration without appealing to it . One promising general approach to the analysis of acceleration has been to analyze the continuous-time limit of accelerated methods , or to derive these limiting ODEs directly via an underlying Lagrangian , and to prove that the ODEs are stable via a Lyapunov function argument. However, these methods stop short of providing principles for deriving a discrete-time optimization algorithm from a continuous-time ODE. There are many ways to discretize ODEs, but not all of them give rise to convergent methods or to acceleration. Indeed, for unconstrained optimization on Euclidean spaces in the setting where the objective is strongly convex, Polyak’s Heavy Ball method and Nesterov’s accelerated gradient descent have the same continuous-time limit. One recent line of attack on the discretization problem is via the use of a time-varying Hamiltonian and symplectic integrators . In this paper, we present a different approach, one based on a fuller development of Lyapunov theory. In particular, we present Lyapunov functions for both the continuous and discrete settings, and we show how to move between these Lyapunov functions. Our Lyapunov functions are time-varying and they thus allow us to establish rates of convergence. They allow us to dispense with estimate sequences altogether, in favor of a dynamical-systems perspective that encompasses both continuous time and discrete time.

A Dynamical View of Momentum Methods

We are concerned with the following class of constrained optimization problems:

which is nonnegative since $h$ is convex. The Euclidean setting is obtained when $h(x)=\frac{1}{2}\|x\|^{2}$.

1 The Bregman Lagrangian

Wibisono, Wilson and Jordan recently introduced the following function on curves,

We introduce a second function on curves,

using the same definitions and scaling conditions. The Lagrangian (5) places a different damping on the kinetic energy than in the original Bregman Lagrangian (2).

Under the same scaling condition (3a), the Euler-Lagrange equation for the second Bregman Lagrangian (5) reduces to:

We provide a proof of Proposition 1 in Appendix A.1. In what follows, we pay close attention to the special case of the dynamics in (6) where $h$ is Euclidean and the damping $\beta_{t}=\gamma t$ is linear:

When $\gamma=\sqrt{\mu}$, we can discretize the dynamics in (7) to obtain accelerated gradient descent in the setting where $f$ is $\mu$-strongly convex.

2 Lyapunov function for the Euler-Lagrange equation

The existence of such a Lyapunov function guarantees that the dynamical system converges: if the function is positive yet strictly decreasing along all trajectories, then the dynamical system must eventually approach a region where ${\mathcal{E}}(X)$ is minimal. If this region coincides with the stationary points of the dynamics, then all trajectories must converge to a stationary point. We now discuss the derivation of time-dependent Lyapunov functions for dynamical systems with bounded level sets. The Lyapunov functions will imply convergence rates for dynamics (2) and (6).

Assume $f$ is convex, $h$ is strictly convex, and the second ideal scaling condition (3b) holds. The Euler-Lagrange equation (4) satisfies

when $x=x^{\ast}$. If the ideal scaling holds with equality, $\dot{\beta}_{t}=e^{\alpha_{t}}$, the solutions satisfy (8) for $\forall x\in\mathcal{X}$. Thus,

A similar proposition holds for the second family of dynamics (5) under the additional assumption that $f$ is $\mu$-uniformly convex with respect to $h$:

When $h(x)=\frac{1}{2}\|x\|^{2}$ is the Euclidean distance, (10) is equivalent to the standard assumption that $f$ is $\mu$-strongly convex. Another special family is obtained when $h(x)=\frac{1}{p}\|x\|^{p}$, which, as pointed out by Nesterov [20, Lemma 4], yields a Bregman divergence that is $\sigma$-uniformly convex with respect to the $p$-th power of the norm:

where $\sigma=2^{-p+2}$. Therefore, if $f$ is uniformly convex with respect to the Bregman divergence generated by the $p$-th power of the norm, it is also uniformly convex with respect to the $p$-th power of the norm itself. We are now ready to state the main proposition for the continuous-time dynamics.

Assume $f$ is $\mu$-uniformly convex with respect to $h$ (10), $h$ is strictly convex, and the second ideal scaling condition (3b) holds. Using dynamics (6), we have the following inequality:

for $x=x^{\ast}$. If the ideal scaling holds with equality, $\dot{\beta}_{t}=e^{\alpha_{t}}$, the inequality holds for $\forall x\in\mathcal{X}$. In sum, we can conclude that

The proof of both results, which can be found in Appendix A.2, uses the fundamental theorem of calculus and basic properties of dynamics (6). Taking $x=x^{\ast}$ and writing the Lyapunov property $\mathcal{E}_{t}\leq\mathcal{E}_{0}$ explicitly,

for (12), allows us to infer a $O(e^{-\beta_{t}})$ convergence rate for the function value for both families of dynamics (4) and (6).

So far, we have introduced two families of dynamics (4) and (6) and illustrated how to derive Lyapunov functions for these dynamics which certify a convergence rate to the minimum of an objective function $f$ under suitable smoothness conditions on $f$ and $h$. Next, we will discuss how various discretizations of dynamics (4) and (6) produce algorithms which are useful for convex optimization. A similar discretization of the Lyapunov functions (9) and (12) will provide us with tools we can use to analyze these algorithms. We defer discussion of additional mathematical properties of the dynamics that we introduce—such as existence and uniqueness—to Appendix C.4.

Discretization Analysis

In this section, we illustrate how to map from continuous-time dynamics to discrete-time sequences. We assume throughout this section that the second ideal scaling (3b) holds with equality, $\dot{\beta}_{t}=e^{\alpha_{t}}$.

The implicit Euler method, on the other hand, evaluates the vector field at the future point

An advantage of the explicit Euler method is that it is easier to implement in practice. The implicit Euler method has greater stability and convergence properties but requires solving an expensive implicit equation. We evaluate what happens when we apply these discretization techniques to both families of dynamics (4) and (6). To do so, we write these dynamics as systems of first-order equations. The implicit and explicit Euler method can be combined in four separate ways to obtain algorithms we can analyze; for both families, we provide results on several combinations of the explicit and implicit methods, focusing on the family that gives rise to accelerated methods.

1 Methods arising from the first Euler-Lagrange equation

We apply the implicit and explicit Euler schemes to dynamics (4), written as the following system of first-order equations:

Wibisono, Wilson and Jordan showed that the polynomial family $\beta_{t}=p\log t$ is the continuous-time limit of a family of accelerated disrete-time methods , Here, we consider any parameter $\beta_{t}$ whose time derivative $\frac{d}{dt}e^{\beta_{t}}=(A_{k+1}-A_{k})/\delta$ can be well-approximated by a discrete-time sequence $(A_{i})_{i=1}^{k}$. The advantage of choosing an arbitrary time scaling $\delta$ is that it leads to a broad family of algorithms. To illustrate this, make the approximations $Z_{t}=z_{k}$, $X_{t}=x_{k}$, $\frac{d}{dt}\nabla h(Z_{t})=\frac{\nabla h(z_{k+1})-\nabla h(z_{k})}{\delta}$, $\dot{X}_{t}=\frac{d}{dt}X_{t}=\frac{x_{k+1}-x_{k}}{\delta}$, and denote $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k}}:=\frac{\alpha_{k}}{A_{k}}$, so that $\frac{e^{\beta_{t}}}{\frac{d}{dt}e^{\beta_{t}}}=\delta/\tau_{k}$. With these approximations, we explore various combinations of the explicit and implicit discretizations.

Written as an algorithm, the implicit Euler method applied to (15a) and (15b) has the following update equations:

We now state our main proposition for the discrete-time dynamics.

Using the discrete-time Lyapunov function,

the bound $\frac{E_{k+1}-E_{k}}{\delta}\leq 0$ holds for algorithm (16).

In particular, this allows us to conclude a general $O(1/A_{k})$ convergence rate for the implicit method (16).

The implicit scheme (16), with the aforementioned discrete-time approximations, satisfies the following variational inequalities:

Using these identities, we have the following derivation:

The inequality on the last line follows from the convexity of $f$ and the strict convexity of $h$. ∎

We study families of algorithms which give rise to a family of accelerated methods. These methods can be thought of variations of the explicit Euler scheme applied to (15a) and the implicit Euler scheme applied to (15b).Here we make the identification $\tau_{k}=A_{k+1}-A_{k}/A_{k+1}:=\alpha_{k}/A_{k+1}$. The first family of methods can be written as the following general sequence:

where $\mathcal{G}:\mathcal{X}\rightarrow\mathcal{X}$ is an arbitrary map whose domain is the previous state, $x=(x_{k+1},z_{k+1},y_{k})$. The second family can be written:

where $\mathcal{G}:\mathcal{X}\rightarrow\mathcal{X}$ is an arbitrary map whose domain is the previous state, $x=(x_{k+1},z_{k},y_{k})$. When $\mathcal{G}(x)=x_{k+1}$ for either algorithm, we recover a classical explicit discretization applied to (15a) and implicit discretization applied to (15b). We will show that the additional sequence $y_{k}$ allows us to obtain better error bounds in our Lyapunov analysis. Indeed, we will show that accelerated gradient descent , accelerated higher-order methods , accelerated universal methods , accelerated proximal methods all involve particular choices for the map $\mathcal{G}$ and for the smoothness assumptions on $f$ and $h$. Furthermore, we demonstrate how the analyses contained in all of these papers implicitly show the following discrete-time Lyapunov function,

is decreasing for each iteration $k$. To show this, we begin with the following proposition.

Assume that the distance-generating function $h$ is $\sigma$-uniformly convex with respect to the $p$-th power of the norm $(p\geq 2)$ (11) and the objective function $f$ is convex. Using only the updates (19a) and (19b), and using the Lyapunov function (21), we have the following bound:

The error bounds in (23) were obtained using no smoothness assumption on $f$ and $h$; they also hold when full gradients of $f$ are replaced with elements in the subgradient of $f$. The proof of this proposition can be found in Appendix B.1. The bounds in Proposition 5 were obtained without using the arbitrary update $y_{k+1}=\mathcal{G}(x)$. In particular, accelerated methods are obtained by picking a map $\mathcal{G}$ that results in a better bound on the error than the straightforward discretization $y_{k+1}=x_{k+1}$. We immediately see that any algorithm for which the map $\mathcal{G}$ satisfies the progress condition $f(y_{k+1})-f(x_{k+1})\propto-\|\nabla f(x_{k+1})\|^{\frac{p}{p-1}}$ or $\langle\nabla f(y_{k+1}),y_{k+1}-x_{k+1}\rangle\propto-\|\nabla f(y_{k+1})\|^{\frac{p}{p-1}}$ will have a $O(1/\epsilon\sigma k^{p})$ convergence rate. We now show how this general analysis applied concretely to each of the aforementioned five methods.

The quasi-monotone subgradient method, which uses the map

for both algorithms (19) and (20), was introduced by Nesterov in 2015.Under this map, assuming the strong convexity of $h$ (which implies $p=2$), we can write the error (23) as

If we assume all the (sub)gradients of $f$ are upper bounded in norm, then maximizing $\sum_{i=1}^{k}\varepsilon_{i}/A_{i}$ results in an $O(1/\sqrt{k})$ convergence rate. This matches the lower bound for (sub)gradient methods designed for Lipschitz-convex functions.The same convergence bound can be shown to hold for the (sub)gradient method under this smoothness class, when one assesses convergence for the average/minimum iterate .

In 1983, Nesterov introduced accelerated gradient decent, which uses the following family of operators $\mathcal{G}\equiv\mathcal{G}_{\epsilon}$, parameterized by a scaling constant $\epsilon>0$:

Nesterov assumed the use of full gradients $\nabla f$ which are $(1/\epsilon)$-smooth; thus, the gradient map is scaled according to the Lipschitz parameter.

Assume $h$ is $\sigma$-strongly convex and $f$ is $(1/\epsilon)$-smooth. Using the gradient update, $y_{k+1}=\mathcal{G}_{\epsilon}(x_{k+1})$, for updates (19c) and (20b), where $\mathcal{G}_{\epsilon}$ is defined in (25), the error for algorithm (19) can be written as follows:

The optimality condition for the gradient update (25) is

The bound (26a) follows from smoothness of the objective function $f$,

For the second bound (26b), we use the $(1/\epsilon)$-smoothness of the gradient,

substituting (27) into (28), squaring both sides, and expanding the square on the left-hand side, yields the desired bound:

The error bounds we have just obtained depend explicitly on the scaling $\epsilon$. This restricts our choice of sequences $A_{k}$; they must satisfy the following inequality:

for the error to be bounded. Choosing $A_{k}$ to be a polynomial in $k$ of degree two, with leading coefficients $\epsilon\sigma$, optimizes the bound (29); from this we can conclude $f(y_{k})-f(x^{\ast})\leq O(1/\epsilon\sigma k^{2})$, which matches the lower bound for algorithms which only use full gradients of the objective function. Furthermore, if we take the discretization step to scale according to the smoothness as $\delta=\sqrt{\epsilon}$, then both $\|x_{k}-y_{k}\|=O(\sqrt{\epsilon})$ and $\varepsilon_{k}=O(\sqrt{\epsilon})$; therefore, as $\sqrt{\epsilon}\rightarrow 0$, we recover the dynamics (15) and the statement $\dot{\mathcal{E}}_{t}\leq 0$ for Lyapunov function (4) in the limit.

Typically, practitioners care about the setting where we have Hölder-continuous gradients ($p=2$) or Hölder-continuous Hessians ($p=3$), since methods which use higher-order information are often too computationally expensive. In the case $p\geq 3$, the gradient update

Assume $f$ has Hölder-continuous higher-order gradients. Using the map $y_{k+1}=\mathcal{G}_{\epsilon,p,\nu,N}(x_{k+1})$, defined by (31), in update (20b) yields the following progress condition:

Lemma 7 demonstrates that if the Taylor approximation is regularized according to the smoothness of the function, the progress condition scales as a function of the smoothness in a particularly nice way. Using this inequality, we can simplify the error (23b) in algorithm (20) to the following,

We end by mentioning that in the special case $p=2$, Nesterov showed that a slightly modified gradient map,

has the following property when applied to functions with Hölder-continuous gradients.

That is, if we take a gradient descent step with increased regularization and assume $h$ is $\sigma$-strongly convex, the error for algorithm (19) when $f$ is $(\epsilon,\nu)$-Hölder-continuous can be written as,

2 Methods arising from the second Euler-Lagrange equation

We apply the implicit and explicit Euler schemes to the dynamics (6) written as the following system of equations:

As in the previous setting, we consider any parameter $\beta_{t}$ whose time derivative $\frac{d}{dt}e^{\beta_{t}}=(A_{k+1}-A_{k})/\delta$ can be well-approximated by a discrete-time sequence $(A_{i})_{i=1}^{k}$. In addition, we make the discrete-time approximations $\frac{d}{dt}\nabla h(Z_{t})=\frac{\nabla h(z_{k+1})-\nabla h(z_{k})}{\delta}$ and $\frac{d}{dt}\dot{X}_{t}=\frac{x_{k+1}-x_{k}}{\delta}$, and denote $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k}}$. We have the following proposition.

Written as an algorithm, the implicit Euler scheme applied to (35a) and (35b) results in the following updates:

Using the following discrete-time Lyapunov function:

we obtain the bound $E_{k+1}-E_{k}\leq 0$ for algorithm (16). This allows us to conclude a general $O(1/A_{k})$ convergence rate for the implicit scheme (16).

The algorithm that follows from the implicit discretization of the dynamics (36) satisfies the variational conditions

where $\tau_{k}=\frac{\alpha_{k}}{A_{k}}$. Using these variational inequalities, we have the following argument:

The inequality uses the Bregman three-point identity (60) and $\mu$-uniform convexity of $f$ with respect to $h$ (10). ∎

We now focus on analyzing the accelerated gradient family, which can be viewed as a discretization that contains easier subproblems.

We study a family of algorithms which can be thought of as slight variations of the implicit Euler scheme applied to (35a) and the explicit Euler scheme applied to (35b)

where $x=(x_{k},z_{k+1},y_{k})$ is the previous state and $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k+1}}$. Note that when $\mathcal{G}(x)=x_{k}$, we recover classical discretizations. The additional sequence $y_{k+1}=\mathcal{G}(x)$, however, allows us to obtain better error bounds using the Lyapunov analysis. To analyze the general algorithm (39), we use the following Lyapunov function:

We begin with the following proposition, which provides an initial error bound for algorithm (39) using the general update (39c).

Assume the objective function $f$ is $\mu$-uniformly convex with respect to $h$ (10) and $h$ is $\sigma$-strongly convex. In addition, assume $f$ is $(1/\epsilon)$-smooth. Using the sequences (39a) and (39b), the following bound holds:

where the error term has the following form:

When $h$ is Euclidean, the error simplifies to the following form

We present a proof of Proposition 10 in Appendix B.4. The result for accelerated gradient descent can be summed up in the following corollary, which is a consequence of Propositions 6 and 10.

for update (39c) results in an error which scales as

The parameter choice $\tau_{k}\leq\sqrt{\mu\epsilon}=1/\sqrt{\kappa}$ ensures the error is non-positive. With this choice, we obtain a linear $O(e^{-\sqrt{\mu\epsilon}k})=O(e^{-k/\sqrt{\kappa}})$ convergence rate. Again, if we take the discretization step to scale according to the smoothness as $\delta=\sqrt{\epsilon}$, then both $\|x_{k}-y_{k}\|=O(\sqrt{\epsilon})$ and $\varepsilon_{k}=O(\sqrt{\epsilon})$, so we recover the dynamics (7) and the continuous Lyapunov argument $\dot{\mathcal{E}}_{t}\leq 0$ in the limit $\sqrt{\epsilon}\rightarrow 0$.

2.2 Quasi-monotone method

We end this section by studying a family of algorithms which can be thought of as a variation of the implicit Euler scheme applied to (35b) and (35b),

where $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k}}:=\frac{\alpha_{k}}{A_{k}}$. In discretization (42a), the state $z_{k+1}$ has been replaced by the state $z_{k}$. When $h$ is Euclidean, we can write (42b) as the following update:

Assume $f$ is $\mu$-strongly convex with respect to $h$ and $h$ is $\sigma$-strongly convex. The following error bound:

can be shown for algorithm (42) using Lyapunov function (37), where the error scales as

No smoothness assumptions on $f$ and $h$ are needed to show this bound, and we can replace all the gradients with subgradients. If we assume that all the subgradients of $f$ are upper bounded in norm, then optimizing this bound results in an $f(x_{k})-f(x^{\ast})\leq O(1/k)$ convergence rate for the function value, which is optimal for subgradient methods designed for strongly convex functions.In particular, this rate is achieved by taking $\tau_{k}=\frac{2}{k+2}$.

3 Frank-Wolfe algorithms

In this section we describe how Frank-Wolfe algorithms can, in a sense, be considered as discrete-time mappings of dynamics which satisfy the conditions,

These dynamics are not guaranteed to exist; however, they are remarkably similar to the dynamics (4), where instead of using the Bregman divergence to ensure nonnegativity of the variational inequality $0\leq\dot{\beta}_{t}e^{\beta_{t}}\langle\nabla f(X_{t}),x-Z_{t}\rangle$, we simply assume (44b) holds on the domain $\mathcal{X}$. We summarize the usefulness of dynamics (44) in the following proposition.

Assume $f$ is convex and the ideal scaling (3b) holds. The following function:

is a Lyapunov function for the dynamics which satisfies (44). We can therefore conclude an $O(e^{-\beta_{t}})$ convergence rate of dynamics (44) to the minimizer of the function.

The proof of this Proposition is in Appendix B.6. Here, we will analyze two Frank-Wolfe algorithms that arise from dynamics (44). Applying the backward-Euler scheme to (44a) and (44b), with the same approximations, $\frac{d}{dt}X_{t}=\frac{x_{k+1}-x_{k}}{\delta}$, $\frac{d}{dt}e^{\beta_{t}}=\frac{A_{k+1}-A_{k}}{\delta}$, and denoting $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k+1}}$, we obtain the variational conditions for the following algorithm:

Update (46a) requires the assumptions that $\mathcal{X}$ be convex and compact; under this assumption, (46a) satisfies

consistent with (44b). The following proposition describes how a discretization of (45) can be used to analyze the behavior of algorithm (46).

Assume $f$ is convex and $\mathcal{X}$ is convex and compact. If f is $(1/\epsilon)$-smooth, using the Lyapunov function,

where the error for algorithm (46) scales as

If instead we assume $f$ has $(\epsilon,\nu)$-Hölder-continuous gradients (30), the error in algorithm (46) now scales as

Taking $x=x^{\ast}$ we infer the convergence rates $O(1/\epsilon k)$ and $O(1/\epsilon k^{\nu})$, respectively. We provide a proof of Proposition 14 in Appendix B.7.

Equivalence to Estimate Sequences

In this section, we connect our Lyapunov framework directly to estimate sequences. We derive continuous-time estimate sequences directly from our Lyapunov function and demonstrate how these two techniques are equivalent.

We provide a brief review of the technique of estimate sequences . We begin with the following definition.

[18, 2.2.1] A pair of sequences $\{\phi_{k}(x)\}_{k=1}^{\infty}$ and $\{A_{k}\}_{k=0}^{\infty}$ $A_{k}\geq 1$ is called an estimate sequence of function $f(x)$ if

The following lemma, due to Nesterov, explains why estimate sequences are useful.

[18, 2.2.1] If for some sequence $\{x_{k}\}_{k\geq 0}$ we have

then $f(x_{k})-f(x^{\ast})\leq A_{k}^{-1}[\phi_{0}(x^{\ast})-f(x^{\ast})]$.

Rearranging gives the desired inequality. ∎

Notice that this definition is not constructive. Finding sequences which satisfy these conditions is a non-trivial task. The next proposition, formalized by Baes in as an extension of Nesterov’s Lemma 2.2.2 , provides guidance for constructing estimate sequences. This construction is used in , and is, to the best of our knowledge, the only known formal way to construct an estimate sequence. We will see below that this particular class of estimate sequences can be turned into our Lyapunov functions with a few algebraic manipulations (and vice versa).

Define recursively $A_{0}=1$, $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k+1}}:=\frac{\alpha_{k}}{A_{k}}$, and

for all $k\geq 0$. Then $\left(\{\phi_{k}\}_{k\geq 0},\{A_{k}\}_{k\geq 0}\right)$ is an estimate sequence.

From (51) and (53), we observe that the following invariant:

where $\varepsilon_{k}\geq 0,\forall k$. Rearranging, we have the following bound:

Notice that an argument analogous to that of Lemma 15 holds:

Rearranging, we obtain the desired bound,

In Table 1 “linear” is defined as $f_{i}(x)=f(x_{i})+\langle\nabla f(x_{i}),x-x_{i}\rangle,$ and “quadratic” is defined as $f_{i}(x)=f(x_{i})+\langle\nabla f(x_{i}),x-x_{i}\rangle+\frac{\mu}{2}\|x-x_{i}\|^{2}.$ The estimate-sequence argument is inductive; one must know the three sequences $\{\varepsilon_{k},A_{k},\phi_{k}(x)\}$ a priori in order to check the invariants hold. This aspect of the estimate-sequence technique has made it hard to discern its structure and scope.

2 Equivalence to Lyapunov functions

We now demonstrate an equivalence between these two frameworks. The continuous-time view shows that the errors in both the Lyapunov function and estimate sequences are due to discretization errors. We demonstrate how this works for accelerated methods, and defer the proofs for the other algorithms discussed earlier in the paper to Appendix C.

The discrete-time estimate sequence (53) for accelerated gradient descent can be written:

Multiplying through by $A_{k+1}$, we have the following argument, which follows directly from our definitions:

The last inequality follows from definition (52). Rearranging, we obtain the inequality $E_{k+1}\leq E_{k}$ for our Lyapunov function (21). Going the other direction, from our Lyapunov analysis we can derive the following bound:

Rearranging, we obtain the estimate sequence (50), with $A_{0}=1$:

Writing $\mathcal{E}_{t}\leq\mathcal{E}_{0}$, one can simply rearrange terms to extract an estimate sequence:

Comparing this to (55), matching terms allows us to extract the continuous-time estimate sequence $\{\phi_{t}(x),e^{\beta_{t}}\}$, where $\phi_{t}(x)=f(X_{t})+e^{-\beta_{t}}D_{h}(x,Z_{t})$.

Further Observations

The dynamical perspective can be extended to the derivation and analysis of a range of other methods. In this section, we provide sketches of some of these analyses, providing a detailed treatment in Appendix D.

Methods for minimizing the composite of two convex functions, $\varphi(x)=f(x)+\psi(x)$, were introduced by Nesterov and studied by Beck and Teboulle , Tseng and several others. In Appendix D.1, we present a dynamical perspective on these methods and show how to recover their convergence theory via the Lyapunov functions presented in this paper.

Discussion

The main contributions in this paper are twofold: We have presented a unified analysis of a wide variety of algorithms using three Lyapunov functions–(21), (40) and (47), and we have demonstrated the equivalence between Lyapunov functions and estimate sequences, under the formalization of the latter due to Baes . More generally, we have provided a dynamical-systems perspective that builds on Polyak’s early intuitions, and elucidates connections between discrete-time algorithms and continuous-time, dissipative second-order dynamics. We believe that the dynamical perspective renders the design and analysis of accelerated algorithms for optimization particularly transparent, and we also note in passing that Lyapunov analyses for non-accelerated gradient-based methods, such as mirror descent and natural gradient descent, can be readily derived from analyses of gradient-flow dynamics.

We close with a brief discussion of some possible directions for future work. First, we remark that requiring a continuous-time Lyapunov function to remain a Lyapunov function in discrete time places significant constraints on which ODE solvers can be used. In this paper, we show that we can derive new algorithms using a restricted set of ODE techniques (several of which are nonstandard) but it remains to be seen if other methods can be applied in this setting. Techniques such as the midpoint method and Runge Kutta provide more accurate solutions of ODEs than Euler methods . Is it possible to analyze such techniques as optimization methods? We expect that these methods do not achieve better asymptotic convergence rates, but may inherit additional favorable properties. Determining the advantages of such schemes could provide more robust optimization techniques in certain scenarios. In a similar vein, it would be of interest to analyze the symplectic integrators studied by within our Lyapunov framework.

Several restart schemes have been suggested for the strongly convex setting based on the momentum dynamics (4). In many settings, while the Lipschitz parameter can be estimated using backtracking line-search, the strong convexity parameter is often hard—if not impossible—to estimate . Therefore, many authors have developed heuristics to empirically speed up the convergence rate of the ODE (or discrete-time algorithm), based on model misspecification. In particular, both Su, Boyd, and Candes and Krichene, Bayen and Bartlett develop restart schemes designed for the strongly convex setting based on the momentum dynamics (4). Our analysis suggests that restart schemes based on the dynamics (6) might lead to better results.

Earlier work by Drori and Teboulle , Kim and Fessler , Taylor et al , and Lessard et al have shown that optimization algorithms can be analyzed by solving convex programming problems. In particular, Lessard et al show that Lyapunov-like potential functions called integral quadratic constraints can be found by solving a constant-sized semidefinite programming problem. It would be interesting to see if these results can be adapted to directly search for Lyapunov functions like those studied in this paper. This would provide a method to automate the analysis of new techniques, possibly moving beyond momentum methods to novel families of optimization techniques.

Acknowledgements

We would like to give special thanks to Andre Wibisono as well as Orianna Demassi and Stephen Tu for the many helpful discussions involving this paper. ACW was supported by an NSF Graduate Research Fellowship. This work was supported in part by the Army Research Office under grant number W911NF-17-1-0304 and by the Mathematical Data Science program of the Office of Naval Research.

References

Appendix A Dynamics

We compute the Euler-Lagrange equation for the second Bregman Lagrangian (5). Denote $z=x+e^{-\alpha_{t}}\dot{x}$. The partial derivatives of the Bregman Lagrangian can be written,

We also compute the time derivative of the momentum $p=\frac{\partial\mathcal{L}}{\partial v}(X_{t},\dot{X}_{t},t)$,

The terms involving $\frac{d}{dt}\nabla h(X)$ cancel and the terms involving the momentum will simplify under the scaling condition (3a) when computing the Euler-Lagrange equation $\frac{\partial\mathcal{L}}{\partial x}(X_{t},\dot{X}_{t},t)=\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial v}(X_{t},\dot{X}_{t},t)$. Compactly, the Euler-Lagrange equation can be written

It is interesting to compare with the partial derivatives of the first Bregman Lagrangian (2),

as well as the derivative of the momentum,

For Lagrangian (2), not only do the terms involving $\frac{d}{dt}\nabla h(X)$ cancel when computing the Euler-Lagrange equation, but the ideal scaling will also force the terms involving the momentum to cancel as well.

A.2 Deriving the Lyapunov functions

We demonstrate how to derive the Lyapunov function (21) for the momentum dynamics (4); this derivation is similar in spirit to the Lyapunov analysis of mirror descent by Nemirovski and Yudin. Denote $Z_{t}=X_{t}+e^{-\alpha_{t}}\dot{X}_{t}$. We have:

Using this identity, we obtain the following argument:

Here (58a) uses the momentum dynamics (15b) and (15a). The inequality (58b) follows from the convexity of $f$. If $\dot{\beta}_{t}=e^{\alpha_{t}}$, simply by rearranging terms and taking $x=x^{\ast}$, we have shown that the function (9) has nonpositive derivative for all $t$ and is hence a Lyapunov function for the family of momentum dynamics (4). If $\dot{\beta}_{t}\leq e^{\alpha_{t}}$, the Lyapunov function is only decreasing for $x=x^{\ast}$.

A.2.2 Proof of Proposition 3

We demonstrate how to derive the Lyapunov function (12) for the momentum dynamics (6). Using the same identity (57), we have the following initial,

will now be useful. Proceeding from the last line, we have

The first inequality follows from the $\mu$-uniform convexity of $f$ with respect to $h$. The second inequality follows from nonnegativity of the Bregman divergence, and the ideal scaling condition (3b), where we must take $x=x^{\ast}$ if $\dot{\beta}_{t}\leq e^{\alpha_{t}}$.

Appendix B Algorithms derived from dynamics (4)

We show the initial bounds (23a) and (23b). We begin with algorithm (19):

The first inequality follows from the $\sigma$-uniform convexity of $h$ with respect to the $p$-th power of the norm and the last inequality follows from the Fenchel Young inequality. If we continue with our argument, and plug in the identity (23a), it simply remains to use our second update (19a):

From here, we can conclude $E_{k+1}-E_{k}\leq\varepsilon_{k}$ using the convexity of $f$.

We now show the bound (23b) for algorithm (20) using a similar argument.

The first inequality follows from the uniform convexity of $h$ and the second uses the Fenchel Young inequality and definition (23b). Using the second update (20a), we obtain our initial error bound:

The last line can be upper bounded by the error $\varepsilon_{k+1}$ using convexity of $f$.

B.2 Proof of Proposition 7

A similar progress bound was proved in Wibisono, Wilson and Jordan [34, Lem 3.2]. Note that $y=\mathcal{G}(x)$ satisfies the optimality condition

Furthermore, since $\nabla^{p-1}f$ is Hölder-continuous (30), we have the following error bound on the $(p-2)$-nd order Taylor expansion of $\nabla f$,

Substituting (62) to (63) and writing $r=\|y-x\|$, we obtain

Now the argument proceeds as in . Squaring both sides, expanding, and rearranging the terms, we get the inequality

Substituting the values $A=\frac{\epsilon}{2N}\|\nabla f(y)\|_{*}^{2}$ and $B=\frac{1}{2N\epsilon}(N^{2}-1)$ from (65), we obtain

B.3 Proof of Universal Gradient Method

We present a convergence rate for higher-order gradient method $y_{k+1}=\mathcal{G}_{\epsilon,p,\nu,N}(x_{k+1})$ where $\mathcal{G}$ is given by (31) and $f$ has $(\epsilon,\nu)$-Hölder-continuous gradients (30). The proof is inspired by the proof of the rescaled gradient flow $\dot{X}_{t}=-\nabla f(X_{t})/\|\nabla f(X_{t})\|_{\ast}^{\frac{p-2}{p-1}},$ outlined in [34, Appendix G], which is the continuous-time limit of the algorithm. Using the Lyapunov function

the following argument can be made using the convexity of $f$ and the dynamics:

B.4 Proof of Proposition 10

We show the initial error bound (41). To do so, we define the Lyapunov function,

The first inequality uses the $\sigma$-strong convexity of $h$ and the Fenchel-Young inequality. The second inequality uses the $\mu$-strong convexity of $f$ with respect to $h$. The third inequality uses the strong convexity of $f$ and $\sigma$-strong convexity of $h$. The following line uses the Bregman three point identity (60) and the subsequent inequality uses the strong convexity of $f$. The last line follows from the smoothness of $f$. Now we turn to the case where $h$ is Euclidean (so $\sigma=1$):

In the second line we have expanded the square. The last line uses the update (39a).

B.5 Proof of Proposition 12

We show the convergence bound for the quasi-monotone method (42). We have,

The first inequality from the strong convexity of $h$ as well as Hölder’s inequality. The second inequality from the uniform convexity of $f$ with respect to $h$ and convexity of $f$. The last line follows from the Bregman three-point identity (60) and non-negativity of the Bregman divergence. Taking $\tau_{k}=\frac{A_{k+1}-A_{k}}{A_{k}}$ gives the desired error bound.

B.6 Proof of Proposition 13

We show that (47) is a Lyapunov function for dynamics (44). The argument is simple:

B.7 Proof of Proposition 14

If we take $\nu=1$ bound (49) implies (48); therefore we simply show the bound (49). To that end,

The first inequality follows from the Hölder continuity and convexity of $f$. The rest simply follows from plugging in our identities.

Appendix C Estimate Sequences

The discrete-time estimate sequence (53) for quasi-monotone subgradient method can be written:

Multiplying through by $A_{k+1}$, we have

Rearranging, we obtain our Lyapunov argument $E_{k+1}\leq E_{k}+\varepsilon_{k+1}$ for (21):

Going the other direction, from our Lyapunov analysis we can derive the following bound:

Rearranging, we obtain our estimate sequence (50) ($A_{0}=1$) with an additional error term:

C.2 Frank-Wolfe

The discrete-time estimate sequence (53) for conditional gradient method can be written:

Multiplying through by $A_{k+1}$, we have

Rearranging, we obtain our Lyapunov argument $E_{k+1}-E_{k}\leq\varepsilon_{k+1}$ for (47) :

Going the other direction, from our Lyapunov analysis we can derive the following bound:

Rearranging, we obtain our estimate sequence (50) ($A_{0}=1$) with an additional error term:

Given that the Lyapunov function property allows us to write

we can extract $\{f(X_{t}),e^{\beta_{t}}\}$ as the continuous-time estimate sequence for Frank-Wolfe.

C.3 Accelerated gradient descent (strong convexity)

The discrete-time estimate sequence (53) for accelerated gradient descent can be written:

Summing over the right-hand side, we obtain the estimate sequence (50):

Since the Lyapunov function property allows us to write

we can extract $\{f(X_{t})+\frac{\mu}{2}\|x-Z_{t}\|^{2},e^{\beta_{t}}\}$ as the continuous-time estimate sequence for accelerated gradient descent in the strongly convex setting.

C.4 Existence and uniqueness

In this section, we show existence and uniqueness of solutions for the differential equations (6), when $h$ is Euclidean. To do so, we write the dynamics as the following system of equations

Inverting the first of these relations, we get

Computing the time-dilated Euler-Lagrange equation, we get

for the first equation, as well as the identity

which is the Euler-Lagrange equation for the sped-up curve, where the ideal scaling holds with equality. Finally, we mention that we can deduce the existence/uniqueness of solution for the proximal dynamics (74) and (80) from the existence/uniqueness of solution for dynamics (4) and (6), given the difference between these dynamics is that (74) (80) have an extra Lipschitz-continuous vector field. Thus, the Cauchy-Lipschitz theorem can be readily applied to the proximal dynamics and the same arguments can be made regarding time-dilation.

Appendix D Additional Observations

In 2009, Beck and Teboulle introduced FISTA, which is a method for minimizing the composite of two convex functions

Define $f=\varphi+\psi$ and assume $\varphi$ and $\psi$ are convex. Under the ideal scaling condition (3b), Lyapunov function (9) can be used to show that solutions to dynamics

satisfy $f(X_{t})-f(x^{\ast})\leq O(e^{-\beta_{t}})$.

The first line follows from the Bregman identity (57). The second line plugs in the dynamics (80a) and (80b). The third lines follows from (58). The fourth and fifth lines follow from convexity. The sixth line plugs in the dynamics (80b) and the last line follows from application of the chain rule. ∎

Next, to show results for dynamics when subgradients of the function are used, we adopt the setting of Su, Boyd and Candes [30, p.35]. First, we define the subgradient through the following lemma.

This guarantees the existence of a directional derivative. Now we establish the following theorem (similar to [30, Thm 24]):

Given the sum of two convex functions $f(x)=\varphi(x)+\psi(x)$ with directional subgradient $G_{\psi}(x,v)$, assume that the second-order ODE

admits a solution $X_{t}$ on $[0,\alpha)$ for some $\alpha>0$. Then for any $0<t<\alpha$, we have $f(X_{t})-f(x)\leq O(1/t^{p})$.

We follow the framework of Su, Boyd and Candes [30, pg. 36]. It suffices to establish that our Lyapunov function is monotonically decreasing. Although $\mathcal{E}_{t}$ may not be differentiable, we can study $\mathcal{E}(t+\Delta t)-\mathcal{E}(t))/\Delta t$ for small $\Delta t>0$. For the first term, note that

where the second line follows since we assume $f$ is locally Lipschitz. The $o(\Delta t)$ does not affect the function in the limit:

The second term, $D_{h}(x,X_{t}+\frac{t}{p}\dot{X}_{t})$, is differentiable, with derivative $-\left\langle\frac{d}{dt}\nabla h(Z_{t}),x-Z_{t}\right\rangle$. Hence,

The last two inequalities follows from the convexity of $f=\varphi+\psi$. In the last inequality, we have used the identity $Z_{t}-X_{t}=\frac{t}{p}\dot{X}_{t}$ in the term $pt^{p-1}\langle G_{\psi}(X_{t},Z_{t}-X_{t}),Z_{t}-X_{t}\rangle$. Combining everything we have shown

which along with the continuity of $\mathcal{E}_{t}$, ensures $\mathcal{E}_{t}$ is a non-increasing of time. ∎

Now we will discretize the dynamics (74). We assume the ideal scaling (3b) holds with equality. Using the same identifications $\dot{X}_{t}=\frac{x_{k+1}-x_{k}}{\delta}$, $\frac{d}{dt}\nabla h(Z_{t})=\frac{\nabla h(z_{k+1})-\nabla h(z_{k})}{\delta}$ and $\frac{d}{dt}e^{\beta_{t}}=\frac{A_{k+1}-A_{k}}{\delta}$ , we apply the implicit-Euler scheme to (74b) and the explicit-Euler scheme to (74a). Doing so, we obtain a proximal mirror descent update,

and the sequence (19a), respectively. We write the algorithm as

where we have similarly substituted the state $x_{k}$ with a sequence $y_{k}$, and added the update $y_{k+1}=\mathcal{G}(x)$. We summarize how the initial bound scales for algorithm (76) in the following proposition.

Assume $h$ is strongly convex, $\varphi$ is $(1/\epsilon)$-smooth and $\psi$ is simple but not necessarily smooth. Using the Lyapunov function (21), the following initial bound

can be shown for algorithm (76), where the error scales as

Tseng [32, Algorithm 1] showed that the map

can be used to simplify the error to the following,

Notice that the condition necessary for the error to be non-positive is the same as the condition for accelerated gradient descent (29). Using the same polynomial, we can conclude an $O(1/\epsilon\sigma k^{2})$ convergence rate.

We begin with the observation that the update (77) and the convexity of $\psi$ allow us to show the inequality

Thus we can conclude $A_{k+1}\psi(y_{k+1})-A_{k}\psi(y_{k})\leq\alpha_{k}\psi(z_{k+1})$. With this, the standard Lyapunov analysis follows:

The first inequality uses the identity (79). The second inequality follows from the convexity of $\psi$. The last line uses the $\frac{1}{\epsilon}$-smoothness of $\varphi$. It simply remains to use the $\sigma$-strong convexity of $h$ and the identities (19a) and $x_{k+1}-y_{k+1}=\tau_{k}(z_{k+1}-z_{k})$. Continuing from the last line, and using these properties, we have

The last line follows from the convexity of $\varphi$. ∎

D.1.2 Strongly convex functions

We study the problem of minimizing the composite objective $f=\varphi+\psi$ in the setting where $\varphi$ is $(1/\epsilon)$-smooth and $\mu$-strongly convex and $\psi$ is simple but not smooth. Like the setting where $f$ is weakly convex, we begin with the following proposition concerning dynamics that are relevant for this setting.

Define $f=\varphi+\psi$ and assume $\varphi$ is $\mu$-strongly convex with respect to $h$ and $\psi$ is convex. Under the ideal scaling condition (3b), Lyapunov function (12) can be used to show that solutions to dynamics,

satisfy $f(X_{t})-f(x)\leq O(e^{-\beta_{t}})$.

The second line comes from plugging in dynamics (80b). The third line uses the Bregman three-point identity (60). We continue by using the strong convexity assumption:

The fourth line follows the strong convexity of $\varphi$ and convexity of $\psi$. The fifth line (second inequality) uses the convexity of $\psi$ once again. The third inequality plugs in the definition of $Z_{t}-X_{t}$ and the second-last inequality follows from the chain rule and the ideal scaling condition (3b). ∎

Assume $h$ is Euclidean and the ideal scaling (3b) holds with equality $\dot{\beta}_{t}=e^{\alpha_{t}}$. To discretize the dynamics (80b), we split the vector field (80b) into two components, $v_{1}(x,z,t)=\dot{\beta}_{t}(X_{t}-Z_{t}-(1/\mu)\nabla\varphi(X_{t}))$, and $v_{2}(x,z,t)=-\dot{\beta}_{t}/\mu\nabla\psi(Z_{t})$ and apply the explicit Euler scheme to $v_{2}(x,z,t)$ and the implicit Euler scheme to $v_{1}(x,z,t)$, with the same identification $\dot{\beta}_{t}=\tau_{k}/\delta$ for both vector fields.While using the same identification of $\dot{\beta}_{t}$ for both vector fields is problematic—since one is being evaluated forward in time and the other backward in time—the error bounds only scale sensibly in the setting where $\dot{\beta}_{t}=\gamma\leq\sqrt{\mu}$ is a constant. This results in the proximal update

We summarize how the initial bound changes with this modified update in the following proposition.

Assume $h$ is Euclidean, $\varphi$ is strongly convex, $\varphi$ is $(1/\epsilon)$-smooth, and $\psi$ is convex and simple. Using the Lyapunov function (40), we have

The condition necessary for the error to be non-positive, $\tau_{k}\leq\sqrt{\epsilon\mu}=1/\sqrt{\kappa},$ results in a $O(e^{-k/\sqrt{\kappa}})$ convergence rate. This matches the lower bound for the class of $(1/\epsilon)$-smooth and $\mu$-strongly convex functions. As in continuous time, this analysis also allows for the use of subgradients of $\psi$.

The first inequality follows from the strong convexity and $(1/\epsilon)$-smoothness of $\varphi$ and (79), from which we can conclude $\varphi(y_{k+1})-\varphi(y_{k})\leq-\tau_{k}\varphi(y_{k})-\tau_{k}\varphi(z_{k+1})$. The second inequality follows from the convexity of $\psi$. The third inequality uses the strong convexity of $f$. Next, we use identity (60) and the smoothness of $\varphi$ to simplify the bound as follows:

D.2 Stochastic methods

Assume $h$ is $\sigma$-strongly convex and $f$ is convex. For algorithm (19), where stochastic gradients are used instead of full gradients and $\mathcal{G}(x)=x_{k+1}$ , we can show the following error bound:

for Lyapunov function (17), where the error scales as

For algorithm (42), where stochastic gradients are used instead of full gradients, we can show the following error bound:

for Lyapunov function (40), where the error scales as

The proof of this claim follows from the proof of Proposition 5 and 12, where we simply take $\nabla f$ to be stochastic. Maximizing over this sequence gives a $O(1/\sqrt{k})$ for the first algorithm and $O(1/k)$ for the second. This convergence rate is optimal and matches the rate of SGD. Notice, however, that the convergence rate is for the entire sequence of iterates, unlike SGD.

Having introduced the dynamics (6), it is clear that the following stochastic dynamics

where the inequality follows from the proof of proposition 3 which can be found in Appendix A.2.2. That is, we can conclude

In particular, choosing $\beta_{t}=2\log t+\log(1/2)$, we obtain a $O(1/t^{2})+O(1/t)$ convergence rate. We can compare this upper bound to the bound (85) with the identifications $\dot{\beta}_{t}=\tau_{k}$ and $e^{\beta_{t}}=A_{k}$.