The connections between Lyapunov functions for some optimization algorithms and differential equations

Introduction

This paper studies Lyapunov functions for differential equations with damping, their discretizations, and optimization algorithms.

which is of course the result of applying Euler’s rule, with step-size $\alpha_{k}>0$, to the gradient system

The value of $f$ decreases along solutions $x(t)$ of this system and, correspondingly, it may be hoped that, for GD, $f(x_{k+1})\leq f(x_{k})$ for sufficiently small $\alpha_{k}$. In fact, that is the case for $\alpha_{k}<2/L$ if $f$ is $L$-smooth, i.e. if $\nabla f(x)$ is $L$-Lipschitz continuous. In this paper we are mainly interested in problems where $f$ belongs the set $\mathcal{F}_{m,L}$ of $m$-strongly convex and $L$-smooth functions, a class that plays an important role in optimization . For $f$ in this class and the constant step-size $\alpha=2/(m+L)$, GD has a bound [19, Theorem 2.1.15]

where $x^{\star}$ is the (unique) minimizer of $f$ and $\kappa=L/m\geq 1$ is the condition number of $f$.

The $1-\mathcal{O}(1/\kappa)$ rate of decay in $f$ in the preceding bound is unsatisfactory because in many applications of interest one has $\kappa\gg 1$. It is possible to improve on GD by resorting to accelerated algorithms with rates $1-\mathcal{O}(1/\sqrt{\kappa})$; for instance, for the method

introduced by Nesterov, it may be shown [19, Theorem 2.2.3] that, if $y_{0}=x_{0}$,

The factor $1-\sqrt{1/\kappa}$ here is close to the optimal possible factor $(1-\sqrt{1/\kappa})^{2}/(1+\sqrt{1/\kappa})^{2}$ one can achieve for minimization algorithms when $f\in\mathcal{F}_{m,L}$ [19, Theorem 2.1.13]. The algorithm (1.2) is also related to ODEs, because it may be seen as a discretization of of the Polyak damped oscillator equation

whose solutions $x(t)$ approach $x^{\star}$ as $t\rightarrow\infty$ if $f$ is $m$-strongly convex [32, Proposition 3].

In recent years, there has been a revived interest, beginning with , in the connections between differential equations and optimization algorithms (see also ). In particular, there has been several papers (see e.g. ) that proposed accelerated algorithms, both in Euclidean and non- Euclidean geometry, based on discretizations of second order dissipative ODEs. The structure of these ODEs and the fact that they can been viewed as describing Hamiltonian systems with dissipation, led to a number of research works that tried to construct or explain optimization algorithms using concepts such as shadowing , symplecticity , discrete gradients , and backward error analysis .

A common feature of the analysis presented in many of the papers mentioned above was the construction of a discrete Lyapunov function that was used in order to deduce the convergence rate of the underlying algorithm. In a general analysis of optimization methods based on the derivation of Lyapunov functions that mimic ODE Lyapunov functions was carried out; that paper presents a Lyapunov function for (1.4). A Lyapunov function for (1.2) may be seen in , where it was also used to study stochastic versions of the algorithm. The paper , among other contributions, constructs a Lyapunov function for a one-parameter family of optimization algorithms that includes (1.2) as a particular case. Outside the field of optimization, Lyapunov functions are important in establishing ergodicity of random dynamical systems , as well as ergodicity of Markov Chain Monte Carlo algorithms, see for example . The construction of Lyapunov functions for optimization algorithms from the perspective of control theory was the subject of study in . The authors extend the work in and derive Linear Matrix Inequalities (LMIs) that guarantee the existence of suitable Lyapunov functions that may be used to establish the convergence rate of the algorithm under study. In addition, develops an LMI framework to construct Lyapunov functions for systems of ODEs. Typically, the LMIs that appear in this context have been solved numerically in the literature.

For $f\in\mathcal{F}_{m,L}$, we use the LMI framework from to derive analytically Lyapunov functions for a two-parameter family of Nesterov optimization methods (see (3.1) below); this family includes the one-parameter family of algorithms in . In this way we find, as a function of the two parameters in (3.1), a convergence rate for the methods in the family. It turns out that the best convergence rate is achieved when the parameters are chosen as in (1.2). The relation between the Lyapunov function constructed in the present work and its counterpart in is discussed in Remark 3.5.

By taking an appropriate limit of the parameters as in e.g. the optimization algorithms in the family may be seen as discretizations of second-order ODEs of the form

where $\bar{b}>0$ is a friction parameter. We obtain analytically Lyapunov functions for (1.5) and determine, as a function of $\bar{b}$, a convergence rate of $f$ to $f(x^{\star})$ along solutions $x(t)$. We prove that the value $\bar{b}=2$ in the Polyak ODE (1.4) yields the optimal convergence rate if $f$ is $m$-strongly convex. Additionally we show that if one is to take explicitly into account the value of $L$ into this calculation, the optimal value of $\bar{b}$ becomes strictly larger than $2$ and yields slightly better convergence rates.

We show that, in the limit where the optimization algorithms approximate the ODEs, the discrete Lyapunov functions converge to the ODE Lyapunov function. Using this correspondence we show, by means of the Heavy Ball method and other examples, that typically, optimization algorithms that are discretizations of (1.5) do not possess discrete Lyapunov functions that mimic the Lyapunov function of the differential equation in item 2 above and lead to acceleration. This emphasizes the well-known fact that, when designing optimization methods, it is not sufficient to ensure that the algorithm may be seen as a consistent discretization of a well-behaved ODE. Unfortunately, discretizations do not necessarily inherit the good long-time properties of the differential equation, as seen for example in the case of discretization of gradient flows , and Hamiltonian problems .

The rest of the paper is organized as follows. In Section 2 we briefly review the approach in that provides a basis for our constructions. In Section 3 we find analytically Lyapunov functions/rates of convergence for a two-parameter family of optimization methods that contains (1.2) as a particular case. Section 4 analyzes the ODE (1.5) and Section 5 studies the connection between the discrete and continuous Lyapunov functions. The Heavy Ball method and other methods that do not possess suitable Lyapunov functions are discussed in Section 6. Finally, we present in the appendix the calculations that allows us to deduce that while the choice $\bar{b}=2$ in (1.5) is optimal if $f$ is only assumed to be $m$-strongly convex, slightly better rates of convergence may be achieved for $f\in\mathcal{F}_{m,L}$ by taking $\bar{b}>2$.

Preliminaries

We will now briefly describe the framework introduced in for the construction of Lyapunov functions of optimization methods and differential equations. The presentation here is adapted from the material in to suit our specific needs.

The following material is limited to results needed to study strongly convex optimization. However the LMI approach in also works in convex optimization.

Optimization algorithms can often be represented as linear dynamical systems interacting with one or more static nonlinearities (see ). In this paper we will consider first-order algorithms that have the following state-space representation

As example, consider algorithms of the well-known form ()

in the optimization context $u^{\star}=0$, and $y^{\star}=x^{\star}$ is the minimizer sought.

To study the convergence rate of optimization algorithms, considers functions of the form

where $a_{0}>0$ and $P$ is positive semi-definite (denoted by $P\succeq 0$). If along the trajectories of (2.1)

we can conclude that $\rho^{-2k}a_{0}(f(x_{k})-f(x^{\star}))\leq V_{k}(\xi_{k})\leq V_{0}(\xi_{0})$ or

If $\rho<1$, we have found a convergence rate for $f(x_{k})$ towards the optimal value $f(x^{\star})$. The following theorem defines an LMI that, when $f\in\mathcal{F}_{m,L}$, guarantees that the property (2.4) holds and therefore (2.3) provides a Lyapunov function for the system .

Then, for $f\in\mathcal{F}_{m,L}$, the sequence $\{x_{k}\}$ satisfies

2 Continuous-time systems

We also consider continuous-time dynamical systems in state space form (throughout the paper we often use a bar over symbols related to ODEs)

in our context $u^{\star}=0$ and $y^{\star}=x^{\star}$. We can replicate the convergence analysis of the discrete case using now functions of the form

where $\lambda>0$. If $\bar{P}\succeq 0$ and, along solutions, $(d/dt)\bar{V}(\xi(t))\leq 0$, then we have $\bar{V}(\xi(t))\leq\bar{V}(\xi(0))$ which in turns implies

The following theorem similarly to the discrete time case, formulates an LMI that guarantees the existence of such a Lyapunov function.

Suppose that, for (2.6), there exist $\lambda>0$, $\bar{P}\succeq 0$, and $\sigma\geq 0$ that satisfy

Then the following inequality holds for $f\in\mathcal{F}_{m,L}$, $t\geq 0$,

A Lyapunov function for Nesterov’s optimization algorithm

We study the optimization method (cf. (2.2))

$k=0,1,\dots$, with parameters $\alpha>0$ and $\beta$. As noted before, the choice $\beta=0$ gives GD and $\beta\neq 0$ corresponds to Nesterov’s accelerated algorithm.

and the divided difference, $k=0,1,\dots$,

the recursion (3.1) may be rewritten ($k=0,1,\dots$)

For future reference, it is useful to observe that, from a dimensional analysis point of view, $m$, $L$ and $1/\alpha$ have the dimensions of the quotient $f/\|x\|^{2}$. Therefore $\delta$ is a non-dimensional version of $\sqrt{\alpha}$. The parameter $\beta$ is non-dimensional. The divided difference (3.2) shares the dimensions of $x$.

For $\beta=0$ (gradient descent), the first equation in (3.3) is a reformulation of the second: it would be more natural to use the simpler state $\xi_{k}=x_{k}$.

The matrix $A$ in (3.4) is a Kronecker product of a $2\times 2$ matrix and $I_{d}$,

the factor $I_{d}$ originates from the dimensionality of the decision variable $x$ and the $2\times 2$ factor is independent of $d$ and arises from the optimization algorithm. The matrices $B$, $C$ and $E$ have a similar Kronecker product structure. It is then natural to consider symmetric matrices $P$ of the form

and then $T$ will also have a Kronecker product structure

where the $t_{ij}$ are explicitly given by the following complicated expressions obtained from (3.4) and the recipes for $M^{(0)}$, $M^{(1)}$ and $M^{(2)}$ in Theorem 2.2:

Our task is to find $\rho\in[0,1)$, $p_{11}$, $p_{12}$, and $p_{22}$ that lead to $\widehat{T}\preceq 0$ and $\widehat{P}\succeq 0$ (which imply $T\preceq 0$ and $P\succeq 0$ ). The algebra becomes simpler if we represent $\beta$ and $\rho^{2}$ as:

Note that we are interested in $r\in(0,1/\delta]$ so as to get $\rho^{2}\in[0,1)$. We proceed in steps as follows.

First step. Impose the condition $t_{23}=0$. This leads to

Second step. Impose the condition $t_{13}=0$. This results in

Third step. Impose the condition ${\rm det}(\widehat{P})=p_{11}p_{22}-p_{12}^{2}=0$. Using (3.9) and (3.10), we have a linear equation for $p_{22}$ with solution

We now take this value to (3.9) and (3.10) and get

a matrix that is positive semi-definite (but not positive definite).

Fourth step. Impose $t_{33}\leq 0$. After using (3.11) in the expression for $t_{33}$ in (3.7), this condition is seen to be equivalent to $\alpha^{2}L-\alpha\leq 0$ or

(for $\alpha=1/L$, $t_{33}$ actually vanishes). In what follows we assume that this bound on $\alpha$ holds; note that then $\delta=\sqrt{m\alpha}\leq\sqrt{m/L}<1$.

Fifth step. We impose $t_{22}\leq 0$. This may be written as $(p_{22}-m/2)r\delta\leq 0$, which leads to $p_{22}\leq m/2$. From (3.11)

which sets a lower limit $\rho^{2}\leq 1-\delta$ for the rate of convergence. For $r^{2}<1$, $t_{22}<0$.

Sixth step. Impose $t_{11}t_{22}-t_{12}^{2}=0$. From (3.11) and (3.7), some algebra yields

Since $\delta<1$ and, after step five, $r\in(0,1]$, we must have $\Xi=0$. For fixed $\delta\in(0,1)$, the condition $\Xi_{\delta}=0$ establishes a relation between the values of $r$ and $b$ or, in other words, the rate of convergence $\rho^{2}$ and the parameter $\beta$ in (3.1). In order to study this relation, we now make a digression and describe, for fixed $\delta\in(0,1)$, the algebraic curve of equation $\Xi_{\delta}(r,b)=0$ in the real plane $(r,b)$; in this description we allow arbitrary real values of $r$ and $b$ (even though in our problem $r\in(0,1]$).

The formula for the roots of a quadratic equation yields

We now return to the construction of $T$. Recall that for our purposes, we need $r>0$ (so as to have $\rho<1$); this requirement holds for $b\in(b_{\rm min},b_{\rm max})$, where

are the intersections of the curve $\Xi_{\delta}=0$ with the vertical axis. As $\delta\downarrow 0$,

The limits on $b$ just found are equivalent to

For the maximum value $r=1$ found in step five above, the formula (3.13) gives the double root $b=2/(1+\delta)$ or $\beta=(1-\delta)/(1+\delta)$. Values $r\in(0,1)$ correspond to two different choices of $b\in(b_{\rm min},b_{\rm max})$.

We are now ready to present the following result.

Consider the minimization algorithm (3.1) (or (3.3)) with parameters subject to

Set $\delta=\sqrt{m\alpha}$ and let $r>0$ be the value determined by $\Xi_{\delta}(r,b)=0$ (see (3.12)), set $\rho^{2}=1-r\delta<1$ and define the positive semi-definite matrix $P$ by (3.5) and (3.11). Then the matrix $T$ in (3.6)–(3.7) is negative semi-definite.

As a result, for any $x_{-1}$, $x_{0}$, the sequence

decreases monotonically, which, in particular, implies

Using Theorem 2.2, we only have to prove that $\widehat{T}\preceq 0$. The second, first and fourth steps of our construction respectively ensure that $t_{13}=t_{23}=0$ and $t_{33}\leq 0$ and therefore we are left with the task of checking that the $2\times 2$ matrix $\widehat{T}^{12}$ obtained by suppressing the last row and last column of $\widehat{T}$ is $\preceq 0$. If $r<1$, we know from step five that $t_{22}<0$ and from step six that the determinant of $\widehat{T}^{12}$ vanishes and therefore $\widehat{T}^{12}\preceq 0$. For $r=1$, $t_{22}=0$, but again $\widehat{T}^{12}\preceq 0$, because in this case $t_{11}=-(m/2)\delta(1-\delta)^{3}/(1+\delta)<0$. ∎

For fixed $\alpha\leq 1/L$, as noted above, $\rho^{2}$ is minimized by the choice

When $\alpha$ is allowed to vary in the interval $(0,1/L]$, increasing $\alpha$ results in an improvement of $\rho^{2}$, so that the best rate $\rho^{2}=1-\sqrt{m/L}=1-\sqrt{1/\kappa}$ is obtained by setting $\alpha=1/L$ and then (3.1) coincides with (1.2). The parameter values $\alpha=1/L$, $\beta=(1-\sqrt{1/\kappa})/(1+\sqrt{1/\kappa})$ in (1.2) are of course the “standard” choice for Nesterov’s algorithm (see e.g. [15, Proposition 12]). For this choice of parameters and $x_{-1}=x_{0}$, the bound in Theorem 3.3 exactly coincides (including the value of $C$) with that in (1.3), which is derived in [19, Theorem 2.2.3] without using Lyapunov functions. Numerical experiments in show that for $\kappa^{-1}=m/L$ small the rate of convergence $\rho^{2}=1-\sqrt{1/\kappa}$ is essentially the best that the algorithm achieves.

The theorem may also be applied to the GD algorithm with $\beta=0$ and $b=1/\delta$, even though (see Remark 3.2) in this case the preceding treatment is unnatural. One finds $r=\delta$, so that the decay per step in $f(x_{k})-f(x_{\star})$ provided by Theorem 3.3 is $\rho^{2}=1-\delta^{2}=1-m\alpha$, for $\alpha\leq 1/L$. When $\alpha=2/(m+L)$, the decay per step guaranteed by Theorem 3.3 is $\rho^{2}=\frac{1-1/\kappa}{1+1/\kappa}$; this is worse than the bound in (1.1) valid for the same value of $\alpha$.

The decay rate $\rho^{2}$ provided by the theorem is a non-dimensional quantity that only depends on the non-dimensional variables $b$ and $\delta$. The bound $\alpha\leq 1/L$ may be rewritten in the non-dimensional form as $\delta^{2}\leq m/L=1/\kappa$. These facts guarantee that the theorem is equivariant with respect to changes in scale of $f$ and $x$. The Lyapunov function in (3.16) has the dimensions of $f$ because, according to (3.11), $P$ has the dimensions of $m$, i.e. those of $f/\|x\|^{2}$.

For the particular choice of $\alpha$ and $\beta$ leading to (1.2), the Lyapunov function in the theorem above was derived in by means of an alternative technique (see Remark 5.2). In a Lyapunov function that contains the gradient $\nabla f(x)$ is constructed analytically for the situation where the learning rate $\alpha$ in (3.1) is a free parameter and the momentum parameter is fixed as $\beta=(1-\sqrt{m\alpha})/(1+\sqrt{m\alpha})$ (i.e. at the value that according to the analysis above optimizes $\rho^{2}$). The analysis in requires (see Lemma 3.4 in that reference) $\alpha\leq 1/(4L)$, while here $\alpha\leq 1/L$. In addition for $\alpha=1/(4L)$, [28, Theorem 3] proves a rate $1/(1+(1/12)\sqrt{m/L})$ which, while establishing acceleration, compares unfavourably with the value $1-(1/2)\sqrt{m/L}$ provided by Theorem 3.3.

2 Optimality

The path leading to Theorem 3.3 has a degree of arbitrariness and it may be asked whether, by following an alternative construction, it is possible to determine the parameters $\rho$, $p_{11}$, $p_{12}$, $p_{22}$ and in such a way that $\widehat{T}\preceq 0$, $\widehat{P}\succeq 0$ and the value of $\rho$ is larger than the value provided in Theorem 3.3. We conclude this section by presenting a result in this direction. We fix the parameters in the algorithm at the standard choices i.e. $\alpha=1/L$, $\beta=(1-\delta)/(1+\delta)$, $\delta=\sqrt{m/L}$, and denote by $\rho^{\star}=\sqrt{1-\delta}$, $p^{\star}_{11}=(m/2)(1-\delta)^{2}$, $p^{\star}_{12}=(m/2)(1-\delta)$, $p^{\star}_{22}=m/2$ the values yielded by Theorem 3.3. In the space of the decision variables $\rho$, $p_{11}$, $p_{22}$, $p_{33}$ we pose the convex optimization problem of minimizing $\rho$ subject to the constraints $\widehat{T}\preceq 0$, $\widehat{P}\succeq 0$. We then have the following result that shows that the rate provided in Theorem 3.3 cannot be improved with an alternative choice of $\widehat{P}$.

With the notation just described, the unique solution of the minimization problem is $(\rho^{\star},p_{11}^{\star},p_{12}^{\star},p_{22}^{\star})$.

We use the notation $\sigma=\rho^{2}$, $\sigma^{\star}=(\rho^{\star})^{2}$ and write $\sigma=\sigma^{\star}+\widetilde{\sigma}$, $p_{11}=p_{11}^{\star}+\widetilde{p}_{11}$, $p_{12}=p_{12}^{\star}+\widetilde{p}_{12}$, $p_{22}=p_{22}^{\star}+\widetilde{p}_{22}$. Since the minimization problem is convex, it is sufficient to show that $\rho^{\star}$, $p_{11}^{\star}$, $p_{12}^{\star}$, $p_{22}^{\star}$ provide a local minimum, i.e. that if the increments $\widetilde{\sigma}\leq 0$, $\widetilde{p}_{11}$, $\widetilde{p}_{12}$, $\widetilde{p}_{22}$ are of sufficiently small magnitude and $(\sigma,p_{11},p_{12},p_{22})$ is feasible, then $\sigma=\sigma^{\star}$, $p_{11}=p_{11}^{\star}$, $p_{12}=p_{12}^{\star}$, $p_{22}=p_{22}^{\star}$.

We study three requirements that feasibility imposes on $\widetilde{\sigma}$, $\widetilde{p}_{11}$, $\widetilde{p}_{12}$, $\widetilde{p}_{22}$.

(1) First, the constraint $\widehat{P}\succeq 0$ implies that $p_{11}p_{22}-p_{12}^{2}\geq 0$ or

Because we are carrying a local study, we replace the constraint by its linearization

or, after using the known values of the symbols with a star,

(2) Then, the constraint $\widehat{T}\preceq 0$ implies $t_{22}t_{33}-t_{23}^{2}\geq 0$ or, using (3.7),

This time the leading terms in the right hand-side are quadratic in the increments and we discard the cubic terms to get:

By completing the square in the quadratic form, this may be equivalently rewritten as

(3) Finally $\widehat{T}\preceq 0$ requires $t_{22}\leq 0$ or $\widetilde{p}_{22}(\delta-\widetilde{\sigma})\leq 0$; discarding the quadratic term, we get

The proof concludes by applying the lemma below. ∎

If the increments $\widetilde{\sigma}\leq 0$, $\widetilde{p}_{11}$, $\widetilde{p}_{12}$, $\widetilde{p}_{22}$ satisfy the constraints (3.17)–(3.20), then $\widetilde{\sigma}=0$, $\widetilde{p}_{11}=0$, $\widetilde{p}_{12}=0$, $\widetilde{p}_{22}=0$.

We combine this inequality with (3.17) to get

Since the three quantities being added in the first bracket in (3.19) are now known to be $\leq 0$, it is enough to consider hereafter the worst case $\widetilde{\sigma}=0$.

Since $\delta\widetilde{p}_{12}+\delta^{2}\widetilde{p}_{22}\leq 0$, we must have

which implies (see (3.20), (3.22), (3.23))

By combining this inequality and (3.18) (with $\widetilde{\sigma}=0$), we obtain a relation

that shows that $\widetilde{p}_{12}=0$. Then comparing (3.17), (3.20) and (3.23), we conclude that $\widetilde{p}_{11}=\widetilde{p}_{22}=0$, which in turn concludes the proof. ∎

The differential equation

which, if $x_{k}$ is seen as an approximation to $x(kh)$, provides a consistent discretization of the differential equation (1.5). An example is provided by the choice $\beta=(1-\delta)/(1+\delta)=(1-\sqrt{m}h)/(1+\sqrt{m}h)$, where $\bar{b}=2$ and (1.5) is the equation (1.4) used by Polyak.

In general, this two-step discretization is, not a linear multistep formula. Note:

$\nabla f$ is evaluated at $y_{k}$, a linear combination of $x_{k}$ and $x_{k-1}$. In this regard, (3.1) is similar to the one-leg methods introduced by Dahlquist in his study of the long-time properties of multistep methods applied to nonlinear differential equations (see e.g. )

The unconventional factor $(1-\beta_{h})/(\sqrt{m}h)$ that converges to $\bar{b}$ as $h\downarrow 0$. From the point of view of discretization methods for ODEs having $\bar{b}$ instead of this factor, or equivalently having $\beta=1-\bar{b}\sqrt{m}h$, would be more natural. But note that, when $\beta=(1-\sqrt{m}h)/(1+\sqrt{m}h)$, the algorithm (3.1) becomes GD for $h=1/\sqrt{L}$ and $\kappa=1$; the choice $\beta=1-\bar{b}\sqrt{m}h$ does not share this favourable property.

and rewrite (1.5) as a first-order system

In a dimensional analysis as in Remarks 3.1 and 3.4, $h$ has the same units as $t$. It is then a dimensional time-step, to be compablue with the non-dimensional $\delta$. The units of $v$ are those of $x$. Of course, the divided difference (3.2) is a discrete version of $v=\dot{x}/\sqrt{m}$.

Now according to Theorem 2.3, in order to find a Lyapunov function of the form (2.7) it is sufficient to find a matrix $\bar{P}\succeq 0$ and parameters $\lambda>0$, $\sigma\geq 0$ such that the matrix $\bar{T}$ in (2.8) is negative semi-definite. Similarly to the discrete case, we will simplify the subsequent analysis by considering the case $\sigma=0$. (The case $\sigma>0$ is studied in the Appendix.) The Lipschitz constant $L$ only enters $T$ in Theorem 2.3 through $\bar{M}^{(3)}$; under the assumption $\sigma=0$, $\bar{T}$ is independent of $L$. This has an important implication: the analysis in this section applies to $f$ strongly $m$-convex but not necessarily $L$-smooth.

where the $\bar{t}_{ij}$ have the following expressions:

We now determine $\lambda$ and $\widehat{\bar{P}}$. The algebra is simplified if we set $\lambda=\sqrt{m}\>{\bar{r}}$.

First step. Since $\bar{t}_{33}=0$, the requirement $\widehat{\bar{T}}\preceq 0$ implies $\bar{t}_{13}=0$ and $\bar{t}_{23}=0$ and accordingly

Second step. We choose $\bar{p}_{22}$ to ensure ${\rm det}(\widehat{\bar{P}})=\bar{p}_{11}\bar{p}_{22}-\bar{p}_{12}^{2}=0$. This yields

a matrix that is positive-semidefinite (but not positive definite).

Third step. Since, $\widehat{\bar{T}}\preceq 0$ implies $\bar{t}_{22}\leq 0$, we may write $0\geq\bar{p}_{22}-m/2=(m/2)({\bar{r}}^{2}-1)$, and therefore we have

this imposes a bound $\lambda\leq\sqrt{m}$ on the convergence rate.

Fourth step. We impose the condition $\bar{t}_{11}\bar{t}_{22}-{\bar{t}}_{12}^{2}=0$. This results in an equation $\bar{\Xi}=0$,

that relates ${\bar{r}}$ (or equivalently the rate $\lambda$) and the parameter $\bar{b}$ in the differential equation (1.5).

We observe that the polynomial $\bar{\Xi}$ is the limit as $\delta\downarrow 0$ of the polynomial $\Xi_{\delta}$ in (3.12) (except of course for the symbols used to denote the variables: $r$ and $b$ for $\Xi_{\delta}$ and ${\bar{r}}$ and $\bar{b}$ for $\bar{\Xi}$). As a consequence, the discontinuous line in Figure 1, presented there as a limit of curves $\Xi_{\delta}=0$, also describes the curve $\bar{\Xi}=0$ (again after renaming the variables).

The curve of equation $\bar{\Xi}({\bar{r}},\bar{b})=0$ in the $({\bar{r}},\bar{b})$ plane is invariant with respect to the symmetry $({\bar{r}},\bar{b})\mapsto(-{\bar{r}},-\bar{b})$ (this is a consequence of the fact that changing $\bar{b}$ into $-\bar{b}$ in the differential equation is equivalent to reversing the sign of independent variable $t$).The curves $\Xi_{\delta}(r,b)=0$, $\delta>0$ do not possess any symmetry because in the discrete algorithm (3.1), $x_{k+1}$ and $x_{k-1}$ do nor play a symmetric role (or in the terminology of differential equation integrators we are not dealing with time-symmetric algorithms). The formula for the roots of a quadratic equation gives

From here one may prove that to each real $\bar{b}$ there corresponds a unique ${\bar{r}}$ such that $\bar{\Xi}({\bar{r}},\bar{b})=0$. The maximum value ${\bar{r}}=1$ ($\lambda=\sqrt{m}$) is achieved only for $\bar{b}=2$ (i.e. for Polyak’s (1.4)) and values ${\bar{r}}\in(0,1)$ correspond to two different real values of $\bar{b}$.

We now have the following result that is proved as in the discrete case.

Consider the differential equation (1.5) (or the equivalent system (4.1)) with parameter $\bar{b}>0$ and assume that $f$ is $m$-strongly convex. Let $\lambda=\sqrt{m}{\bar{r}}$, where ${\bar{r}}>0$ is the value determined by the relation $\bar{\Xi}({\bar{r}},\bar{b})=0$ (see (4.6)) and define the positive semi-definite matrix ${\bar{P}}$ by (4.2) and (4.5). Then the matrix $\bar{T}$ in (4.3) is negative semi-definite.

As a result, if $x(t)$ is a solution of (1.5), the function

decreases monotonically as $t$ increases, which implies

For $\bar{b}=0$, the construction leading to the theorem yields $r=0$, i.e. $\lambda=0$, and,

In addition, $\bar{T}=0$ and therefore the factor in round brackets in (4.7) is an invariant of motion. In this case the system (4.1) is Hamiltonian and the invariant we have found equals $\sqrt{m}$ times the corresponding Hamiltonian function.

The value $\bar{b}=2$, in addition to maximizing the decay rate in $f(x(t))$ in Theorem 4.3 for arbitrary $m$-strongly convex $f$, has another optimality property in the simple one-dimensional case with $f(x)=mx^{2}/2$, when (1.5) or (4.1) describe a damped harmonic oscillator. An elementary computation (see e.g. ) shows that $\bar{b}=2$ is the value of the friction coefficient that ensures the fastest dissipation of the energy $(\dot{x})^{2}/2+mx^{2}/2$.

It will be proved in the Appendix that if $f$, in addition to being strongly convex has Lipschitz continuous gradient, then better decay rates in $f(x(t))$ may be obtained by choosing $\bar{b}$ to be larger than $2$. Therefore $(\dot{x})^{2}/2+mx^{2}/2$ is not the best Lyapunov function to study the rate of decay of $f(x)$ in the damped harmonic oscillator. This is in agreement with Theorem 4.6 below.

Reference gives a Lyapunov function for (1.5) or (4.1) that includes a cross-term $v^{T}\nabla f(x)$ and does not require the strong convexity of $f$. However, the presence of the gradient in the Lyapunov function makes it necessary that $f$ be demanded to be twice-differentiable (the Hessian of $f$ appears when differentiating the Lyapunov function with respect to $t$).

2 Optimality

Steps 2 and 4 in the construction above imply a degree of arbitrariness and it is of interest to ask whether there are alternative choices of $\lambda$ and $\widehat{\bar{P}}\succeq 0$ that, while ensuring $\widehat{\bar{T}}\preceq 0$, furnish better decay rates. We conclude this section by proving that this is not the case.

In the theorem below we use the notation $\bar{r}^{\star}$ and $\widehat{\bar{P}}^{\star}$ for the values obtained, for given $\bar{b}>0$, in the construction leading to Theorem 4.3. (These are functions $\bar{r}^{\star}=\bar{r}^{\star}(b)$ and $\widehat{\bar{P}}^{\star}=\widehat{\bar{P}}^{\star}(b)$, but the dependence on $\bar{b}$ will be dropped from the notation.) In particular, $\bar{p}_{22}^{\star}={m\bar{r}^{\star}}^{2}/2$ and $\bar{\Xi}(\bar{r}^{\star},\bar{b})=0$. The symbols $\lambda$ and $\widehat{\bar{P}}$ are used in the theorem to refer to an arbitrary real number and an arbitrary $2\times 2$ symmetric matrix. Finally, we set $\lambda^{\star}=\sqrt{m}\>\bar{r}^{\star}$ and $\lambda=\sqrt{m}\>\bar{r}$.

With the notation as described, for each fixed $\bar{b}>0$, $\lambda^{\star}={\rm max}\>\lambda$, subject to the constraints $\widehat{\bar{T}}(\lambda,\widehat{\bar{P}})\preceq 0$, $\widehat{\bar{P}}\succeq 0$.

Since we are solving a convex optimization problem, it is sufficient to show that $(\lambda^{\star},\widehat{\bar{P}}^{\star})$ provides a local maximum.

We observed in step 1 above that $\widehat{\bar{T}}\preceq 0$ determines the values of $\bar{p}_{11}$, $\bar{p}_{12}$ as in (4.4). This leaves us with $\lambda$ (or equivalently $\bar{r}$) and $\bar{p}_{22}$ as decision variables. For simplicity we hereafter omit the subindices in $\bar{p}_{22}$.

The constraint $\widehat{\bar{P}}\succeq 0$, implies ${\rm det}(\widehat{\bar{P}})\geq 0$ or (after using the values of $\bar{p}_{11}$, $\bar{p}_{12}$) $\bar{p}\geq(m/2){\bar{r}}^{2}$. The constraint $\widehat{\bar{T}}\preceq 0$ implies $\bar{t}_{11}\bar{t}_{22}-{\bar{t}_{12}}^{2}\geq 0$. We use (4.4), to write $\bar{t}_{11}\bar{t}_{22}-\bar{t}_{12}^{2}\geq 0$ as a function $\Delta(\bar{r},\bar{p})$; tedious algebra leads to the expression:

We will be done if we prove that the pair $(\bar{r}^{\star},\bar{p}^{\star})$ is a local maximum for the problem

At the point $({\bar{r}}^{\star},{\bar{p}}^{\star})$ both constraints are active (in fact they were chosen to be so at steps 2 and 4). If we define the Lagrangian

where $\zeta_{1}$, $\zeta_{2}$ are the multipliers, the proof concludes by showing that the gradient of $\mathcal{L}$ at $({\bar{r}}^{\star},{\bar{p}}^{\star})$ may be annihilated for a suitable choice of positive multipliers.

($|^{\star}$ means evaluation at at $({\bar{r}}^{\star},{\bar{p}}^{\star})$) and

(which implies that $\zeta_{1}$ and $\zeta_{2}$ have the same sign) and eliminate $\zeta_{1}$ to get

In this way we are left with the task of proving that

or, after using the expression for $\Delta$ and some simplification,

Let us denote by $\Lambda=\Lambda({\bar{r}}^{\star},\bar{b})$ the left hand-side of this inequality. When $\bar{b}=2$ and ${\bar{r}}^{\star}=1$, we have $\Lambda=-1$. On the other hand, we know that

and this relation makes it impossible for $\Lambda$ to change sign as $\bar{b}>0$ and the corresponding ${\bar{r}}^{\star}(b)\in(0,1]$ vary. In fact, if $\Lambda$ were to vanish, we would have

something that cannot happen because ${\bar{r}}^{\star}<1$ for $\bar{b}\neq 2$. ∎

Connecting the differential equations with optimization algorithms

The second-order differential equation (1.5) provides a limit for the algorithm (3.1) when $\beta$ changes smoothly with $h=\sqrt{\alpha}$ in such a way that $\beta_{h}=1-\bar{b}\sqrt{m}h+o(h)$ as $h\downarrow 0$. In this section we study this limit when $\bar{b}>0$. As in (3.8) write $\beta_{h}=1-b_{h}\delta=1-b_{h}\sqrt{m}h$. Clearly, $b_{h}\rightarrow\bar{b}$ and, in addition, for $h$ sufficiently small $b_{h}\in(b_{\rm min}^{h},b_{\rm max}^{h})$ (see (3.14)). The application of Theorem 3.3 then gives a rate $\rho^{2}_{h}=1-r_{h}\delta=1-r_{h}\sqrt{m}h$. As noted before, the polynomial $\bar{\Xi}$ in (4.6) is the limit of $\Xi_{\delta}$ in (3.12) as $h$ (or $\delta$) approaches zero, and, accordingly, $r_{h}\rightarrow{\bar{r}}$, where ${\bar{r}}$ solves $\bar{\Xi}({\bar{r}},\bar{b})=0$. Then Theorem 3.3 guarantees that, over one step $k\mapsto k+1$ of the algorithm, $f(x_{k})-f(x^{\star})$ decays by a factor $\rho^{2}_{h}=1-\sqrt{m}{\bar{r}}h+{o}(h)$. Over $k$ steps the decay factor will be $(1-\sqrt{m}{\bar{r}}h+{o}(h))^{k}$, a quantity that in the limit $kh\rightarrow t$ converges to $\exp(-\sqrt{m}{\bar{r}}t)=\exp(-\lambda t)$. This is exactly the decay guaranteed by Theorem 4.3 for $f(x(t))-f(x^{\star})$ over an interval of length $t$.

In addition, the matrices $P_{h}$ in the discrete Lyapunov function converge to the matrix $\widehat{P}$ in the differential equation, because from the expression for the entries in (3.11) and (4.5)

The above discussion and standard results on the convergence of discretizations of ordinary differential equations imply the following result.

Fix the parameter $\bar{b}>0$ and the initial conditions $x(0)$, $\dot{x}(0)$ for the differential equation (1.5). For small $h>0$, consider the optimization algorithm (3.1) with parameters $\alpha=h^{2}$ and $\beta=\beta_{h}=1-\bar{b}\sqrt{m}h+o(h)$. Assume that the initial points $x_{-1}$, $x_{0}$ are such that, as $h\downarrow 0$, $x_{0}\rightarrow x(0)$ and $(1/h)(x_{0}-x_{-1})\rightarrow\dot{x}(0)$. Then, in the limit $kh\rightarrow t$,

$x_{k}\rightarrow x(t)$ and $(1/h)(x_{k+1}-x_{k})\rightarrow\dot{x}(t)$.

The discrete Lyapunov function in (3.16) converges to the Lyapunov function in (4.7).

As a consequende of this theorem, the Lyapunov function of the differential equation could have been derived alternatively by first finding the Lyapunov function for the discrete optimization algorithm and then taking limits. In our research we first investigated the discrete case and then studied the differential equations; in hindsight we saw it would have been easier to first deal with the differential equation and then carry out the analysis of the algorithm by mimicking the treatment of the continuous case. References find Lyapunov functions for different optimization algorithms by first constructing Lyapunov functions for suitable so-called high-resolution differential equations. In our context, this would mean perturbing (4.1) with suitable $h$-dependent terms so as to obtain an ($h$-dependent) differential equation for which the algorithm has a high order of consistency. The idea behind those high-resolution equations is very old in the numerical analysis of ordinary and partial differential equations, where they are known as modified equations, see e.g. or [24, Chapter 10] and, for the stochastic case, .

Heavy Ball and other methods

The paper has given rise to a number of contributions that aim to understand the behaviour of optimization methods by seeing them as discretizations of differential equations. However it is well known that the long-time properties of a differential equation are not automatically inherited by their discretizations, regardless of the value of the step-size chosen. A very simple example is provided by the application of Euler’s rule to the harmonic oscillator: for all step-sizes the discrete trajectories grow while the continuous solutions stay bounded. A more relevant example in an optimization context may be seen in . On the other hand properties of the discretizations may often be extrapolated to the continuous limit; a general discussion of these points in different settings may be seen in .

In the setting of the preceding section, it is not true that discretizing a dissipative differential equation with a known a Lyapunov function will always yield an optimization algorithm with a “suitable” Lyapunov function. We now illustrate this fact by means of the Heavy Ball algorithm obtained by choosing $\gamma=0$ and $\beta\neq 0$ in (2.2).

We proceed as in Section 3, rewrite the algorithm in terms of $d_{k}$ and $x_{k}$ and then cast it in the general format (2.1). We will presently prove that a discrete Lyapunov with properties similar to the Lyapunov function for Nesterov’s method in Theorem 3.3 does not exist. We argue by contradiction. With the notation as in Section 3, we consider

$p_{ij}=m\,\phi_{ij}(\beta,\delta)$, $(i,j)=(1,1),(1,2),(2,2)$, such that $\widehat{P}\succeq 0$,

and suppose that the corresponding $T(\lambda,P)$ is $\preceq 0$ for each $\delta<c/\sqrt{\kappa}$. As in Remark 3.4 to ensure equivariance with respect to changes of scale, the number $c$ and functions $\phi_{ij}$ and $\psi$ are assumed to be independent of the constants $m$ and $L$ associated with $f$ and the values of the parameters $\alpha$ and $\beta$ in the Heavy Ball algorithm.

For future reference, the element $t_{11}$ is found to have the expression:

This has to be $\leq 0$ for $\delta<c/\sqrt{\kappa}$.

Next, as in the preceding section, we assume that $\beta$ changes smoothly with $h$ in such a way that, for some $\bar{b}>0$, $\beta=\beta_{h}=1-\bar{b}\delta+o(h)=1-\bar{b}\sqrt{m}h+o(h)$. Clearly the algorithm is then a consistent discretization of the differential equation (1.5), and we assume that $r_{h}$, $p_{ij}^{h}$ converge to their differential equation counterparts $\bar{r}$ and ${\bar{p}}_{ij}$.This hypothesis is not necessarily in the argument that follows. It is enough to suppose that $r_{h}$, $p_{ij}^{h}$ have finite limits.

This cannot happen because $L$ may be arbitrarily large.

The Heavy Ball algorithm is a “more natural” discretization of (1.5) than Nesterov’s, in that, as conventional linear multistep methods, it does not evaluate $\nabla f$ at a linear combination of $x_{k}$, $x_{k-1}$ (cf. Remark 4.1).

The contradiction in (6.1) arises because we insisted in $T$ being $\preceq 0$ for “large” non-dimensional stepsizes $\delta=\sqrt{m}h<c/\sqrt{\kappa}$. For optimization algorithms that, in the limit $h\downarrow 0$, approximate a differential equation with decay $\exp(-\lambda h)=\exp(-\bar{r}\delta)$ in a time-interval of length $h$, such large stepsizes seem to be necessary to achieve accelerated rates $1-\mathcal{O}(\sqrt{\kappa})$ rather than rates $1-\mathcal{O}(\kappa)$.

The reference constructs a Lyapunov function for the Heavy Ball method, but it only operates for $\delta=\mathcal{O}(1/\kappa)$ and, while useful in showing convergence, does not provide acceleration. For an additional convergence proof of the Heavy Ball algorithm see ; again this reference does not prove acceleration.

The three-parameter family of methods (2.2) contains algorithms, like Nesterov’s, that “inherit” the ODE Lyapunov function for stepsizes $\delta<c/\sqrt{\kappa}$ and algorithms, like the Heavy Ball, that do not. In fact the situation for the Heavy Ball is arguably the rule rather than the exception. For (2.2),

where we observe the unwelcome presence of the factor $L-m$ that created the difficulties in the analysis of the Heavy Ball algorithm. If we look at a situation where $\beta$ changes with $h$ as above and in addition $\gamma$ is also allowed to change with $h$ and approaches a limit, a Lyapunov function that has the form envisaged and works for $\delta<c/\sqrt{\kappa}$ may only exist if $\beta_{h}-\gamma_{h}$ vanishes (at least in the limit $h\downarrow 0$) to offset the factor, i.e. if the algorithm is not far away from Nesterov’s.

Acknowledgement. We are thankful to an anonymous referee for helping us to improve the discussion of our results.

References

Appendix

In Theorem 4.6 we proved that, for each $\bar{b}>0$, the rate of decay $\lambda$ provided by Theorem 4.3 is the best one may obtain by using Theorem 2.3 if one chooses $\sigma=0$. In this Appendix we investigate whether $\lambda$ may be improved by a suitable choice of $\sigma>0$. Since for $\sigma\neq 0$, the matrix $\bar{M}^{(3)}$ that contains the constant $L$ contributes to $T$, the following results require that $f$, in addition to being $m$-strongly convex (as in Theorem 4.3) is $L$-smooth, i.e. they hold for $f\in\mathcal{F}_{m,L}$.

When $\sigma\neq 0$ the expressions for the $t_{ij}$ in Section 4 have to be replaced by:

As in Section 4, we set $\lambda=\sqrt{m}\>{\bar{r}}$ and, in addition, $\sigma=m\bar{s}$ (the variable $\bar{s}$ is, as $\bar{r}$, non-dimensional). We shall show that it is possible, for given $m$ and $L$, to find values of the six parameters ${\bar{p}}_{11}$, ${\bar{p}}_{12}$, ${\bar{p}}_{22}$, $\bar{b}$, $\bar{s}$, $\bar{r}$, in such a way that the constraints $\widehat{\bar{T}}\preceq 0$, $\widehat{\bar{P}}\succeq 0$, $\bar{s}\geq 0$ are satisfied and, at the same time, $\bar{r}>1$, so that by using the matrix ${\bar{M}}^{(3)}$ it is possible to improve on the best value $\bar{r}=1$ (associated with $\bar{b}=2$ and leading to $\lambda=\sqrt{m}$) that may be achieved in Theorem 4.3.

For given $m$ and $L$, we determine the values of the six parameters as follows:

First step. We impose ${\bar{t}}_{22}=0$, a requirement that leads to the relation

Second step. We impose ${\bar{t}}_{23}=0$ and get

Third step. We require ${\rm det}(\widehat{\bar{P}})=0$. Therefore

Note that for $\bar{r},\bar{s}\geq 0$ we have ${\bar{p}}_{22}>0$ and thus the third step guarantees that $\widehat{\bar{P}}\succeq 0$.

Fourth step. We next demand that ${\bar{t}}_{12}=0$ and obtain

The four preceding displayed formulas allow us to express the parameters ${\bar{p}}_{12}$, ${\bar{p}}_{22}$, and $\bar{b}$ as known functions of $\bar{s}$ and $\bar{r}$.

Fifth step. At this stage, we have ensublue that ${\bar{t}}_{12}$, ${\bar{t}}_{22}$, ${\bar{t}}_{23}$ vanish. As a result, the condition $\widehat{\bar{T}}\preceq 0$ is equivalent to $\widehat{\bar{T}}^{13}\preceq 0$ where $\widehat{\bar{T}}^{13}$ is the $2\times 2$ matrix obtained by suppressing from $\widehat{\bar{T}}$ its second row and column. Furthermore ${\bar{t}}_{33}<0$ for $\bar{s}>0$ and then we shall have $\widehat{\bar{T}}^{13}\preceq 0$ if we impose that ${\rm det}(\widehat{\bar{T}}^{13})=0$, or

By using the displayed formulas above, the last equation becomes a relation $F(\bar{r},\bar{s})=0$, between $\bar{r}$ and $\bar{s}$, with

We next show that the rational curve $F(\bar{r},\bar{s})=0$ in the $(\bar{r},\bar{s})$ real plane has points with $\bar{s}>0$ and $\bar{r}>1$.

It is easily checked that the point $\bar{r}=1$, $\bar{s}=0$ lies on the curve $F=0$ and has $\bar{b}=0$. This could have been anticipated because, if $\bar{s}=0$ and $\bar{b}=2$, the construction in this appendix just reproduces the construction in Section 4, which yields $\bar{r}=1$.

By removing the denominator in the rational function $F$ so as to have a polynomial equation for the curve and looking at the Newton diagram at $\bar{r}=1$, $\bar{s}=0$, one sees that in the neighbourhood of this point the curve consists of a single branch that may be parameterized by $\bar{r}$. A Taylor expansion reveals that

In this way, choosing a sufficiently small value of the parameter $\bar{s}>0$, there are two possible values of the rate $\bar{r}$

one of which is $>1$. In conclusion we have proved analytically that the introduction of $\sigma$ and $\bar{M}^{(3)}$ in $T$ makes it possible to achieve rates $\bar{r}>1$ (or $\lambda>\sqrt{m}$).

We next determined the value of $\bar{s}$ that leads to the largest possible $\bar{r}$ on the curve $F=0$. In view of the involved expression of $F$, we proceeded numerically and found this largest value by continuation along the curve, starting from $\bar{r}=1$, $\bar{s}=0$. The results, for different values of $\kappa$, are given in Table 1. For the small condition number $\kappa=10$, the table shows that it is possible to achieve a decay $\approx\exp(-1.086\sqrt{m}t)$ by fixing the dissipation coefficient at the value $\bar{b}\approx 2.35$ rather than at $\bar{b}=2$ as in Polyak’s (1.4)—this is a marginal improvement on the best decay $\exp(-\sqrt{m}t)$ that one may insure without using $\bar{M}^{(3)}$. In addition the improvement quickly decreases as the condition number grows: for $\kappa=10^{3}$ the decay is $\exp(-1.0039\sqrt{m}t)$. In fact, we observe in the table that, as $\kappa\uparrow\infty$, $\bar{r}\approx 1+0.38\kappa^{-2/3}$. Of course as $\kappa$ increases, $\bar{r}$ and $\bar{b}$ approach the values $1$ and $2$ that correspond to the situation studied in Section 4, where $f$ is not assumed to possess Lipschitz gradients. A similar convergence obtains for the matrix $\widehat{\bar{P}}\succeq 0$. Also note that $\bar{s}\approx 0.50\kappa^{-1/3}$: as the condition number increases the parameter $\sigma=\sqrt{m}\bar{s}$ that multiplies $\bar{M}^{(3)}$ decreases, as it may have been expected.