A systematic approach to Lyapunov analyses of continuous-time models in convex optimization

Introduction

Convex optimization is an important tool in the numerical analyst toolbox. It serves, among others, for framing modeling problems in data science and signal processing. A number of convex optimization problem take the form:

where $f$ is convex, differentiable, and $x\in\mathbf{R}^{d}$ contains the variables of the model. First-order methods are very popular to solve these problems, due to their attractive low cost per iteration, and to the fact that data science applications typically do not require very accurate solutions . Gradient descent is a common first-order method, which starts from a point $x_{0}\in\mathbf{R}^{d}$, its iterates are given by the simple recursion

where $\gamma>0$ is a step size. Denoting $\dot{X}_{t}=\frac{d}{dt}X_{t}$, gradient descent with small step sizes $\gamma$ is directly related to the so-called gradient flow:

where the solution $X_{t}$ of the ODE verifies $X_{t_{k}}\approx x_{k}$ with the identification $t_{k}=\gamma k$. In numerical integration, gradient descent (2) is also known as the Euler explicit scheme for integrating the gradient flow. Recently, Su et al. have interpreted Nesterov’s accelerated gradient in a similar fashion through its continuous version, paving the way to several continuous analyses of accelerated methods .

where $\sigma(X_{t})$ is a noise parameter connected to parameters of the method, that was further developed by Shi et al. . This connection has raised many questions regarding approximate equivalences between optimization methods and continuous-time models, and tools for analyzing convergence speeds of discrete methods. Usually, gradient flows and first-order methods are studied via worst-case convergence properties, that have to be verified for any function of a given class, and any trajectory generated by the ODE or the optimization method. In many cases, continuous approaches seem to allows for shorter and simpler proofs, together with intuitions on what can be expected from optimization methods.

The analysis of continuous-time models often relies on Lyapunov stability arguments, as in system theory and physics, where energy dissipation plays a crucial role. The existence of such Lyapunov functions provides direct convergence proof for ODEs under consideration. The main challenge in the Lyapunov approach is to find a suitable function that is decreasing along all trajectories generated by an ODE.

From an outsider point of view, these analyses are often seen as complicated and technical to reach. In this work, we remedy this problem by extending the systematic approach based on semidefinite programming (SDP) developed by Drori and Teboulle for certifying convergence of optimization methods. This technique is referred as “performance estimation problems” (PEPs). The main contribution of this work is to provide a tool for simultaneously analyzing convergence of continuous-time models, and constructing Lyapunov functions suited to gradient flows in a systematic way, using small-sized SDPs reformulations. Furthermore, this procedure benefits from tightness properties, meaning that the feasibility of the SDP allows to conclude about the existence of some Lyapunov functions.

When dealing with convergence rates of gradient flows, many proofs typically rely on a Lyapunov function. In control theory, such functions are common for studying stability of dynamical systems .

Given a trajectory $X_{t}$ , we call $\mathcal{V}:\mathbf{R}^{d}\rightarrow\mathbf{R}$ a Lyapunov function if it is differentiable and satisfies the conditions:

$\frac{d}{dt}\mathcal{V}(X_{t})\leqslant 0$.

Given an ODE starting from $x_{0}$ and a class of functions $\mathcal{F}$, a function $\mathcal{V}(\cdot)$ is a Lyapunov function if the inequality $\frac{d}{dt}\mathcal{V}(X_{t})\leqslant 0$ is verified for all trajectories $X_{t}$ generated by ODEs originating from functions $f\in\mathcal{F}$. There exist similar definitions of Lyapunov functions for discrete optimization methods .

Lyapunov functions are suited to deriving both linear and sublinear convergence rates. When looking for linear convergence rates (as we may expect for strongly convex functions), the third condition is typically replaced by $\frac{d}{dt}\mathcal{V}(X_{t})\leqslant-\tau\mathcal{V}(X_{t})$, where $\tau$ depends on the class of functions and on the ODE.

With this approach, convergence guarantees highly depend on the choice of Lyapunov functions. For a specific ODE and class of functions $\mathcal{F}$, there are often multiple choices of valid Lyapunov functions. In this work, we look for Lyapunov functions in a class of quadratic functions that is popular both in the discrete and continuous time literature.

2 Prior works

Lyapunov functions are common for analyzing continuous and discrete time models in convex optimization. Convergence proofs for Nesterov’s accelerated gradient method typically rely on such Lyapunov functions , [10, Theorem 4.8]. In the recent work , the authors proposed Lyapunov-based analyses for many first-order methods, for linear and sublinear convergence rates. Continuous-time versions of optimization methods also often involve Lyapunov arguments, such as Nesterov’s accelerated gradient flow introduced in , and its high-resolution ODEs for strongly convex functions proposed in , or accelerated mirror descent whose continuous dynamics was analyzed in .

Different techniques were developed to compute suitable Lyapunov functions. The authors of proposed an approach based on Bregman Lagrangian for accelerated methods in potentially non-Euclidean settings, further developed in . directly derived Lyapunov functions from Hamiltonian equations describing dynamics of ODEs. Using similar conservation laws in a dilated coordinate system, also generated Lyapunov functions in a principled way.

Given a class of functions and an optimization method, proving a convergence rate mostly consists in combining inequalities characterizing the class of functions at hand. Recently, the automated search for combination of inequalities formulated as semidefinite program was formalized by Drori and Teboulle , and led to the notion of performance estimation problems. Their work was followed up in to provide worst-case bounds in a principled way, and extended to the Lyapunov framework . A competing strategy inspired by control theory was developed by , where Lyapunov functions for discrete-time models are constructed using integral quadratic constraints (IQCs) and semidefinite programming; a similar approach was applied to continuous-time models in . Connections between Lyapunov functions obtained via the IQC framework in continuous and discrete-time, were later highlighted by .

For stochastic differential equations (SDEs), convergence proofs can also be obtained through the Lyapunov approach, together with Ito’s calculus. analyzed both SGD, SAGA , and SVRG , for some well-chosen Lyapunov functions. extended the framework of to the stochastic setting. To the best of our knowledge, a systematic way of verifying a Lyapunov functions in the stochastic setting has not been developed yet.

3 Contributions and organization

In this work, we are concerned with worst-case convergence analyses of ordinary and stochastic differential equations, for modeling optimization methods. We propose a principled approach to worst-case analyses based on Lyapunov functions, SDPs and Ito’s calculus.

In Section 2, we extend the performance estimation approach developed for optimization methods to gradient flows, that originates from a (possibly strongly) convex function. In short, we find Lyapunov functions as feasible points of certain linear matrix inequalities (LMIs). Building on the first part of this work for ODEs, we analyze continuous versions of stochastic optimization algorithms. All codes for numerical results are provided at https://github.com/CMoucer/PEP_ODEs.

Section 3 studies properties of trajectories generated by SDEs, as approximations to stochastic gradient methods. We obtain a simple version of the trade-off between forgetting the initial conditions and diminishing the noise, with and without averaging. It appears that decreasing step sizes, together with a non-uniform version of averaging, allows to reach an optimal trade-off. Our results match those obtained for the stochastic gradient method, in a compact way compared with discrete analyses.

In Section 4, we prove that accelerated gradient flows require diminishing step sizes to converge in the stochastic setting. In contrast to first-order stochastic gradient flow, averaging does not preserve convergence for fixed step size.

4 Assumptions

Throughout this work, functions to be minimized are convex (see Problem 1). Under this assumption, stationary points are global minimizers. We restrict ourselves to continuous-time versions of gradient descent, accelerated gradient descent and stochastic gradient descent. Such methods gather information about functions to be minimized by evaluating its (sub)gradient at past iterates.

Let us recall a few basic definitions and properties characterizing the classes of functions under consideration within the next sections. A function $f:\mathbf{R}^{d}\rightarrow\mathbf{R}$ is convex if for all $x,y\in\mathbf{R}^{d}$, and for all $\lambda\in$, $f(\lambda x+(1-\lambda)y)\leqslant\lambda f(x)+(1-\lambda)f(y)$. We consider in particular the class of convex closed proper (CCP) functions (i.e., functions whose epigraphs are non-empty closed convex sets). For simplicity, we assume in addition differentiability of $f$, even if results do not require it for convex gradient differential inclusions [4, Section 3.2]. Then, $f$ is convex if and only if for all $x,y\in\mathbf{R}^{d}$, $f(y)\geqslant f(x)+\langle\nabla f(x),y-x\rangle$. Smoothness is a common assumption for analyzing optimization methods, that limits the growth rate of the function. A function $f$ is $L$-smooth if the gradient is $L$-Lipschitz, that is if for any $x,y$,

A differentiable function $f$ is $\mu$-strongly convex if for any $x,y\in\mathbf{R}^{d}$ it satisfies

Strong convexity ensures the function is not too flat, and the unicity of the minimizer $x_{\star}$. We denote by $\mathcal{F}_{\mu,L}$ the family of $L$-smooth $\mu$-strongly closed convex proper functions, with $0\leqslant\mu\leqslant L\leqslant+\infty$. Weaker assumptions than strong convexity are also encountered in the literature for analyzing gradient algorithms, and leads to similar convergence guarantees. Among them, Łojasiewicz introduced the Łojasiewicz inequality, under which Polyak [28, Theorem 4] showed linear convergence of gradient descent. Other relaxed versions of strong convexity followed .

A principled approach to Lyapunov functions for gradient flows

In this section, we study convergence properties of the gradient flow and its accelerated versions, via quadratic Lyapunov functions. We prove that verifying such a Lyapunov function can be formulated as a LMI. This framework allows to search for Lyapunov functions, and to derive convergence bounds for non-autonomous gradient flows.

Without further assumptions, the function $f$ is decreasing along the trajectory $X_{t}$ solution to the gradient flow. The Lyapunov function $\mathcal{V}(X_{t})=f(X_{t})$ has indeed a negative time-derivative $\frac{d}{dt}\mathcal{V}(X_{t})=\dot{X_{t}}^{\top}\nabla f(X_{t})=-\|\nabla f(X_{t}))\|^{2}$.

Under additional convexity assumptions, it is possible to derive Lyapunov functions in a principled way, and deduce worst-case convergence speeds of the gradient flow. As a first stage, let us consider gradient flows originating from strongly convex functions and establish linear (or exponential) convergence of the flows.

Let $f$ be $\mu$-strongly convex ($\mu>0$), admitting thus a unique minimizer $x_{\star}$ such that $f(x_{\star})=f_{\star}$, and consider the Problem (1). Thanks to strong convexity, it is possible to prove linear convergence of the gradient flow to its stationary point. Scieur et al. proved in [33, Proposition 1.1] a convergence bound in function values for gradient flows originating from $f\in\mathcal{F}_{\mu,\infty}$,

This convergence guarantee follows directly from the time-derivative of the Lyapunov function $\mathcal{V}(X_{t})=f(X_{t})-f_{\star}$, together with strong convexity (or Łojasiewicz inequality): $\frac{d}{dt}\mathcal{V}(X_{t})=\dot{X_{t}}^{\top}\nabla f(X_{t})=-\|\nabla f(X_{t})\|^{2}\leqslant-2\mu(f(X_{t})-f_{\star})\leqslant-2\mu\mathcal{V}(X_{t})$.

Given the specific gradient flows studied in this work, it is reasonable to search for Lyapunov functions made of linear combinations of function values, and a quadratic form in the trajectory $X_{t}$. We simply refer to them as quadratic Lyapunov functions:

where $a,c$ are fixed nonnegative constants that do not depend on $t$. Such Lyapunov functions are common for proving convergence of gradient flows (and of optimization methods), and cover for instance the Lyapunov used to prove convergence of the gradient flow under strong convexity (4). Given a Lyapunov function $\mathcal{V}_{a,c}$, the idea is to find the smallest value of $\tau$ such that the condition

holds for any functions $f\in\mathcal{F}_{\mu,\infty}$ and any trajectory $X_{t}$ generated by the gradient flow. After integrating, a convergence guarantee in function values is given by $\mathcal{V}_{a,c}(X_{t})\leqslant e^{-\tau t}\mathcal{V}_{a,c}(X_{0})$. Given a certain Lyapunov function $\mathcal{V}_{a,c}$ and a time $t$, we get that the largest acceptable $\tau$ is a solution to:

This minimization problem is invariant in $t$. It is established in that these so-called performance estimation problems (PEP) can be formulated as SDPs (details are provided in Appendix A.1 for completeness). Because of the condition $f\in\mathcal{F}_{\mu,\infty}$, the maximization problem (7) is infinite-dimensional. Recall that a differentiable function $f\in\mathcal{F}_{\mu,\infty}$ verifies for all points $x,y\in\mathbf{R}^{d}$, $f(x)-f(y)-\langle\nabla f(y),x-y\rangle\geqslant\frac{\mu}{2}\|x-y\|^{2}$. Introducing alternate variables $f_{t}$, $f_{\star}$, $g_{t}$ and $g_{\star}$ (informally: $f_{t}=f(X_{t})$, $f_{\star}=f(x_{\star})$, $g_{t}=\nabla f(X_{t})$ and $g_{\star}=\nabla f(x_{\star})=0$), it holds that

The fact (8) produces an upper bound on $-\tau$ directly follows from the fact that any sampled strongly convex function satisfy these inequalities at the sampled points $X_{t}$ and $x_{\star}$. Thereby, any feasible point to (7) corresponds to a feasible for (8) with the same objective value. In the other direction, [45, Corollary 2] (which provides a constructive way to obtain some $f\in\mathcal{F}_{\mu,\infty}$ that interpolates the triplets $(X_{i},g_{i},f_{i})_{t,\star}$) ensures that any feasible point to (8) can be translated to a feasible point to (7) with also the same objective value, thereby reaching the equivalence between formulations (7) and (8).

In a second stage, we introduce $G=\begin{pmatrix}\|X_{t}-x_{\star}\|^{2}&\langle X_{t}-x_{\star},g_{t}\rangle\\ \langle X_{t}-x_{\star},g_{t}\rangle&\|g_{t}\|^{2}\end{pmatrix}\succeq 0$ a Gram matrix and a vector $F=[f_{t},f_{\star}]$, thereby obtaining a semidefinite reformulation:

where $A_{0}=\begin{pmatrix}c\cdot\tau&-c\\ -c&-a\end{pmatrix}$, $A_{1}=\begin{pmatrix}-\mu/2&1/2\\ 1/2&0\end{pmatrix}$, $A_{2}=\begin{pmatrix}-\mu/2&0\\ 0&0\end{pmatrix}$, $b_{0}=a\cdot\tau[1,\ -1]^{\top}$ $b_{1}=[-1,\ 1]^{\top}$ and $b_{2}=[1,\ -1]^{\top}$. Those developments allow arriving to the following worst-case results.

Let $\mathcal{V}_{a,c}$ be a quadratic Lyapunov function (5) with $a,c\geqslant 0$, $\tau\geqslant 0$, and functions $f\in\mathcal{F}_{\mu,\infty}$ be strongly convex with parameter $\mu>0$. We consider gradient flows (3) originating from functions $f\in\mathcal{F}_{\mu,\infty}$ and starting from an initial point $x_{0}\in\mathbf{R}^{d}$, where $d\in\mathbf{N}$ is the dimension. The following assertions are equivalent:

The inequality $\frac{d}{dt}\mathcal{V}_{a,c}(X_{t})\leqslant-\tau\mathcal{V}_{a,c}(X_{t})$ is satisfied for all $f\in\mathcal{F}_{\mu,\infty}$, all trajectories $X_{t}$ solutions to the gradient flow and all dimensions $d\in\mathbf{N}$.

There exist $\lambda_{1},\lambda_{2}\geqslant 0$ such that

The procedure to obtain the LMI is fully detailed in Appendix A.1, and consists in taking the standard Lagrangian dual of the SDP (9).

A few conclusions can be drawn from the LMI equivalence from Theorem 2.1. First, it provides a necessary and sufficient condition for a quadratic Lyapunov function $\mathcal{V}_{a,c}$ to decrease at a specific rate $\tau\geqslant 0$ for all functions $f\in\mathcal{F}_{\mu,\infty}$. The infeasibility of the LMI allows to conclude about the existence of a Lyapunov function such that the rate $\tau$ is achieved for all functions in the class. Second, it is possible to search over the class of quadratic Lyapunov functions that certify linear convergence. Given a rate $\tau$, the LMI is indeed jointly convex in $\lambda_{1},\lambda_{2},a,c$. Finally, thanks to linearity of the feasibility problem in $\tau$, a bisection search allows to optimize over the convergence rate. In other words, we can optimize jointly over the Lyapunov function and the worst-case guarantee.

In Figure 1(a), we obtain the fastest linear convergence rate that can be achieved using quadratic Lyapunov functions. Together with Theorem 2.1 we retrieve the known linear worst-case convergence speed in $e^{-2\mu}$ from Scieur et al. [33, Proposition 1.1], without improvement. The numerical approach though allows to ensure tightness with a numerical function $f$ that matches this convergence guarantee (see Figure 1(b) and the method in [43, Chapter 3]). Let us build on these results to analyze the gradient flow originating from a (possibly non strongly) convex function, where the difficulty comes from the time-dependence of Lyapunov functions. Following theorems are obtained using the same methodology.

1.2 Minimizing convex functions

In the case where $f\in\mathcal{F}_{0,\infty}$, worst-case convergence rates often are sublinear. Again, as in discrete time, it is possible to obtain convergence guarantees using time-dependent quadratic Lyapunov functions.

The Lyapunov function $\mathcal{V}(X_{t},t)=t(f(X_{t})-f_{\star})+\frac{1}{2}\|X_{t}-x_{\star}\|^{2}$ from [37, p. 7] verifies $\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant 0$ for any functions $f\in\mathcal{F}_{0,\infty}$, and any trajectories $X_{t}$ generated by the gradient flow (3) (proof: $\frac{d}{dt}\mathcal{V}(X_{t})=t\langle\nabla f(X_{t}),\dot{X}_{t}\rangle+f(X_{t})-f_{\star}+\langle\dot{X}_{t},X_{t}-x_{\star}\rangle=-t\|\nabla f(X_{t})\|^{2}+f(X_{t})-f_{\star}-\langle\nabla f(X_{t}),X_{t}-x_{\star}\rangle\leqslant-t\|\nabla f(X_{t})\|^{2}$ using convexity). After integrating, we recover a convergence bound in function values from the literature [37, p.7], [11, Section 6.3.1],

We consider the family of quadratic Lyapunov functions:

where $a_{t},c_{t}$ are differentiable functions from $\mathbf{R}^{+}$ to $\mathbf{R}^{+}$. When the Lyapunov function is decreasing along the trajectory $X_{t}$, that is $\frac{d}{dt}\mathcal{V}(X_{t})\leqslant 0$, a convergence guarantee in function values is given by

Looking for a worst-case guarantee with a Lyapunov approach can be cast as:

The strongly convex case as defined above is a particular case of the convex one, using a specific Lyapunov function $\Phi(\cdot)$, such that $\mathcal{V}(X_{t},t)=e^{\tau t}\Phi(X_{t})$ where $\Phi(X_{t})=a\cdot(f(X_{t})-f_{\star})+c\cdot\|X_{t}-x_{\star}\|^{2}$. Then, $\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant 0$ is equivalent to $\frac{d}{dt}\Phi(X_{t})\leqslant-\tau\Phi(X_{t})$.

Let $\mathcal{V}_{a_{t},c_{t}}$ defined in (11) be a quadratic Lyapunov function for $a_{t},c_{t}$ nonnegative differentiable functions and functions $f\in\mathcal{F}_{0,\infty}$ be convex. Given the gradient flow (3) originating from convex functions $f\in\mathcal{F}_{0,\infty}$ and starting from $x_{0}\in\mathbf{R}^{d}$, where $d\in\mathbf{N}$ is the dimension. The following assertions are equivalent:

The inequality $\displaystyle\frac{d}{dt}\mathcal{V}_{a_{t},c_{t}}(X_{t},t)\leqslant 0$ is satisfied for all functions $f\in\mathcal{F}_{0,\infty}$, for all trajectories $X_{t}$ generated by gradient flows, all functions $f\in\mathcal{F}_{\mu,\infty}$ and for all dimensions $d\in\mathbf{N}$.

There exist $\lambda^{(1)}_{t},\lambda^{(2)}_{t}\geqslant 0$ such that

If these assertions are satisfied, $c_{t}$ is decreasing, and $a_{t}\geqslant tc_{t}$.

The LMI is obtained following the previous methodology (Appendix A.1), and the inequality $a_{t}\geqslant tc_{t}$ directly from the LMI.

Choosing $c_{t}=\frac{1}{2}$ and $a_{t}=t$ for the Lyapunov function parameters, and $\lambda^{(1)}_{t}=1$, $\lambda^{(2)}_{t}=0$, we retrieve the Lyapunov function $\mathcal{V}(x,t)=t(f(x)-f_{\star})+\frac{1}{2}\|x-x_{\star}\|^{2}$ from [37, p. 7], that verifies $f(X_{t})-f_{\star}\leqslant\frac{1}{2t}\|x_{0}-x_{\star}\|^{2}$ for all convex functions $f\in\mathcal{F}_{0,\infty}$, for all trajectories $X_{t}$ generated by the gradient flow, and all dimensions $d\in\mathbf{N}$.

The differential LMI equivalence from Theorem 2.4 is jointly convex in $\lambda^{(1)}_{t}$, $\lambda^{(2)}_{t}$, $c_{t}$, $a_{t}$, $\dot{a}_{t}$, $\dot{c}_{t}$, introducing derivatives (to be compared with iterates for discrete optimization methods). In contrast with the minimization of strongly convex functions, numerically solving this LMI is not available because of the dependence in $t$. Other quadratic Lyapunov functions could be derived. However, the condition $a_{t}\geqslant c_{t}t$ together with $\dot{c}_{t}\leqslant 0$ implies that the convergence cannot be faster than $\frac{1}{t}$.

2 Accelerated gradient flows

A major improvement to gradient descent dates back to Nesterov in 1983 , with an accelerated gradient method (AGM),

where $\gamma,\alpha_{k}$ are nonnegative parameters depending on the class of functions to minimize. The current iterate is computed thanks to a so-called momentum. This combination of past iterates allows more control over the accumulated error. This idea was first introduced by Polyak with the heavy-ball method, starting from $x_{0},x_{1}\in\mathbf{R}^{d}$,

where $\alpha_{k}>0$ is the momentum method. Yet, compared with Nesterov’s accelerated gradient method, the heavy-ball method lacks global acceleration beyond quadratics.

When reducing the step size $\gamma$, these schemes happen to be closely related to second-order differential equations, for $\beta_{t}\geqslant 0$ a function that depends on $\alpha_{k}$,

Recently, accelerated gradient methods have been analyzed using second-order differential equations . Reversely, the accelerated gradient method and the heavy-ball method may be seen as discretization schemes of these second order ODEs, as many other schemes. Discretization techniques are, among others, discussed by . Taking integration theory’s point of view, Scieur et al. proved that these multi-step methods may even be seen as discretization schemes of the gradient flow (for quadratics).

Again, ODEs and multi-step first-order methods as defined above are often handled using quadratic Lyapunov functions. We extend the systematic Lyapunov approach developed previously to accelerated gradient flows. Let $\mathcal{V}_{a_{t},P_{t}}$ be the family of quadratic Lyapunov functions for second-order gradient flows,

where $P=\begin{pmatrix}p^{(11)}_{t}&p^{(12)}_{t}\\ p^{(12)}_{t}&p^{(22)}_{t}\end{pmatrix}$ is a symmetric matrix with differentiable parameters, and $a_{t}$ is a differentiable function, such that the Lyapunov is nonnegative when evaluated on the gradient flow. After integrating, this approach leads to convergence bounds for instance in function values, such that $f(X_{t})-f_{\star}\leqslant\frac{\mathcal{V}(x_{0})}{a_{t}}$.

Let $f\in\mathcal{F}_{\mu,L}$, with strong convexity parameter $\mu>0$, and scheme parameters be defined by $\gamma=\frac{1}{L}$ and $\alpha=\frac{1-\sqrt{\mu\gamma}}{1+\sqrt{\mu\gamma}}$ in Nesterov’s accelerated gradient methods (12). When reducing the step size $\gamma$, the continuous-time limit of $y_{k}$ in (12) is exactly the Polyak damped oscillator ,

where $t=\gamma k$, as it has already been highlighted in previous works . This ODE is also the limit of the heavy-ball method (13). Shi et al. proved a convergence guarantee in $f(X_{t})-f_{\star}={O}(e^{-\frac{\sqrt{\mu}t}{4}})$ using a Lyapunov-based approach. This bound was improved to $f(X_{t})-f_{\star}={O}(e^{-\sqrt{\mu}t})$, by Wilson et al. [47, Appendix B] using the Bregman-Lagrangian approach, and by Sanz-Serna and Zygalakis using the IQC framework . We compute linear convergence guarantees using quadratic Lyapunov functions with constant parameters (14).

Let $\mathcal{V}_{a,P}$ be a quadratic Lyapunov function (14), where $a\geqslant 0$, and $P$ a symmetric matrix. Let $\tau\geqslant 0$, $\mu\geqslant 0$, $d\in\mathbf{N}$ be the dimension, and the Polyak damped oscillator (15) be starting from $x_{0}\in\mathbf{R}^{d}$, and originating from strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$. The following assertions are equivalent:

The inequality $\displaystyle\frac{d}{dt}\mathcal{V}_{a,P}(X_{t})\leqslant-\tau\mathcal{V}_{a,P}(X_{t})$ is satisfied for all functions $f\in\mathcal{F}_{\mu,\infty}$, all trajectories $X_{t}$ generated by Polyak damped oscillator, and all dimensions $d\in\mathbf{R}^{d}$.

There exist $\lambda_{1},\lambda_{2},\nu_{1},\nu_{2}\geqslant 0$, such that

The LMI equivalence is obtained introducing the Gram matrix $G=P^{\top}P$, where $P=(\dot{X}_{t},X_{t}-x_{\star},g_{t})$ and using the PEP methodology. The first LMI refers to the non-increasing condition for the Lyapunov function, and the second LMI to the positivity of the Lyapunov for all trajectories and convex functions.

The LMI is a feasibility problem, jointly convex in $\lambda_{1},\lambda_{2},\nu_{1},\nu_{2}\geqslant 0$ and in Lyapunov parameters $a,P$. Given parameters $a,P$ and a rate $\tau$, it is a verification tool for $\mathcal{V}_{a,P}$ to be a Lyapunov function decreasing at rate $\tau$. As for gradient flows, we can perform a bisection search over $\tau$ to find the fastest linear convergence rate that can be verified using quadratic Lyapunov functions (see Figure 2(a)). It provides a numerical tool for choosing Lyapunov parameters (Figure 2(b)) for which the worst-case linear convergence rate is achieved.

Let us consider the Polyak damped oscillator (15) starting from $x_{0}\in\mathbf{R}^{d}$ where $d\in\mathbf{N}$ is the dimension, and originating from strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$ where $\mu>0$. The Lyapunov function

verifies $\frac{d}{dt}\mathcal{V}(X_{t})\leqslant-{4}/{3}\sqrt{\mu}\mathcal{V}(X_{t})$ for all functions $f\in\mathcal{F}_{\mu,\infty}$, all trajectories $X_{t}$ generated by Polyak damped oscillator, and all dimensions $d\in\mathbf{N}$.

Taking $\lambda_{1}={4}/{3}\sqrt{\mu}$, $\lambda_{2}=0$, $\nu_{1}=0$ and $\nu_{2}=1$, we verify the LMI for this Lyapunov function $\mathcal{V}$, with $\tau={4}/{3}\sqrt{\mu}$.

This class of quadratic Lyapunov functions is inspired by in discrete time, where a stricter positivity condition on $P\succeq 0$ hindered proving tight convergence of Nesterov’s accelerated gradient. Similarly, the Lyapunov function from Corollary 15 is defined by $P=\begin{pmatrix}{4}/{9}\mu&{2}/{3}\sqrt{\mu}\\ {2}/{3}\sqrt{\mu}&{1}/{2}\end{pmatrix}$, which is not positive semidefinite. In the continuous-time models’ literature, we usually choose matrices $P$ positive semidefinite, such as in the Lyapunov function from [32, 34, Theorem 4.3],

that verifies $\displaystyle\frac{d}{dt}\mathcal{V}(X_{t})\leqslant-\sqrt{\mu}\mathcal{V}(X_{t})$ for all functions $f\in\mathcal{F}_{\mu,\infty}$, all trajectories $X_{t}$ generated by Polyak damped oscillator, and all dimensions $d\in\mathbf{N}$. This Lyapunov is a feasible point of the LMI (LABEL:LMI_acc) from Theorem 2.6, with $\tau=\sqrt{\mu}$, $\lambda_{1}=\sqrt{\mu}$, $\lambda_{2}=0$, $\nu_{1}=0$, $\nu_{2}=a=1$. Relaxing the condition $P\succeq 0$ improves results from Sanz-Serna and Zygalakis and Wilson et al. by a factor $4/3$ in Corollary 2.8.

Figure 2(b) provides a numerical help for computing the Lyapunov function from Corollary 2.8. Figure 2(a) together with Corollary 2.8 allows to conclude that this bound cannot be improved when changing the Lyapunov function among the class of quadratic functions (14).

2.2 Minimizing convex functions

As for gradient flow, rates are sublinear when the accelerated gradient flow originates from convex functions. Let $f\in\mathcal{F}_{0,L}$, the step size be $\gamma\leqslant\frac{1}{L}$, and $\alpha=\frac{k-1}{k+2}$ be the scheme parameter in Nesterov’s accelerated method. Su et al. [37, Section 2] proved the connection between the first-order scheme and a second order ODE known as accelerated gradient flow (AGF):

Su et al. [37, Theorem 3] proved the following inequality is verified for all functions $f$ and all trajectories $X_{t}$ generated by AGF,

Their proof exhibits a Lyapunov function $\mathcal{V}(X_{t},t)=t^{2}(f(X_{t})-f_{\star})+\|(X_{t}-x_{\star})+\frac{t}{2}\dot{X}_{t}\|^{2}$, whose derivative is decreasing along trajectories $X_{t}$ (proof: $\frac{d}{dt}\mathcal{V}(X_{t},t)=2t(f(X_{t})-f_{\star})+t^{2}\langle\dot{X}_{t},\nabla f(X_{t})\rangle+2\langle X_{t}-x_{\star}+\frac{t}{2}\dot{X}_{t},3\dot{X}_{t}+t\ddot{X}_{t}\rangle=2t(f(X_{t})-f_{\star})-\langle\nabla f(X_{t}),X_{t}-x_{\star}\rangle\leqslant 0$ by convexity of $f$). The following theorem provides a systematic condition for a quadratic function $\mathcal{V}$ (14) to be a Lyapunov function for AGF (17).

Let $\mathcal{V}_{a_{t},P_{t}}$ be a quadratic Lyapunov function (14), where $a_{t}\geqslant 0$ is a differentiable function, and $P_{t}\succeq 0$ with differentiable parameters. Given the accelerated gradient flow (17) starting from $x_{0}\in\mathbf{R}^{d}$ where $d\in\mathbf{N}$ is the dimension, and originating from convex functions $f\in\mathcal{F}_{0,\infty}$, the following assertions are equivalent:

The inequality $\displaystyle\frac{d}{dt}\mathcal{V}_{a_{t},P_{t}}(X_{t},t)\leqslant 0$ is satisfied for all functions $f\in\mathcal{F}_{0,\infty}$, all trajectories $X_{t}$ generated by the accelerated gradient flow and all dimensions $d\in\mathbf{N}$.

There exist $\lambda^{(1)}_{t},\lambda^{(2)}_{t}\geqslant 0$ such that

The proof follows those from Theorem 2.6 and 2.1.

The LMI from Theorem 2.10 is differentiable and convex in $\lambda^{(1)}_{t}$, $\lambda^{(2)}_{t}$. Even when studying second-order gradient flows, the number of interpolation inequalities (or of dual values $\lambda_{t}^{(i)}$) is bounded by two, enforcing short proofs.

Theorem 2.10 allows to retrieve the Lyapunov function exhibited by Su et al. [37, Theorem 3], and its associated convergence guarantee. Choosing $a_{t}=t^{2}$ and $P_{t}=2\begin{pmatrix}1&t/2\\ t/2&t^{2}/4\end{pmatrix}$ for the Lyapunov parameters, and $\lambda^{(1)}_{t}=t$, $\lambda^{(2)}_{t}=0$, we prove that that the Lyapunov function \mathcal{V}(X_{t},t)=t^{2}(f(X_{t})-f_{\star})+2\big{\|}(X_{t}-x_{\star})+\frac{t}{2}\dot{X}\big{\|}^{2} verifies $\displaystyle\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant 0$, for all functions $f\in\mathcal{F}_{0,\infty}$, for all trajectories $X_{t}$ generated accelerated gradient flows, and all dimensions $d\in\mathbf{N}$. As for gradient flows originating from convex functions, we did not compute numerically worst-case scenarios because of the time-dependence of the Lyapunov function.

3 Higher-order convergence and time dilation

In this section, we analyze convergence of non-autonomous gradient flows, and provide convergence rates depending on their parameters. It appears that higher order convergence of gradient flows is highly connected to time dilation.

Let $f$ be a convex function, and $X_{t}$ be the solution to the non-autonomous first-order gradient flow,

where $\alpha_{t}\geqslant 0$ is a continuous function (such that the flow is converging). It is natural to wonder if it is possible to accelerate such gradient flows when changing $\alpha_{t}$. A change of variable connects this ODE to the gradient flow (3), for which $\alpha_{t}=1$. Let $Y_{t}$ be the solution to the gradient flow, and $\tau_{t}=\int_{0}^{t}\alpha_{s}ds$ be a time change variable. Then, the variable $X_{t}=Y_{\tau_{t}}$ verifies $\dot{X}_{t}=\frac{d}{dt}Y_{\tau_{t}}=\alpha_{t}\dot{Y}_{\tau_{t}}=-\alpha_{t}\nabla f(Y_{\tau_{t}})=-\alpha_{t}\nabla f(X_{t})$, which is exactly the non-autonomous gradient flow. The following corollaries can be obtained by performing the appropriate change of variable in LMI from Theorem 2.4.

Let us consider non-autonomous gradient flows (18) starting from $x_{0}\in\mathbf{R}^{d}$ where $d\in\mathbf{N}$ is the dimension, and originating from possibly strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$, where $\mu\geqslant 0$.

If $\mu=0$, the Lyapunov function $\mathcal{V}(X_{t},t)=\left(\int_{0}^{t}\alpha_{s}ds\right)(f(X_{t})-f_{\star})+\frac{1}{2}\|X_{t}-x_{\star}\|^{2}$ verifies $\displaystyle\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant 0$ for all functions $f\in\mathcal{F}_{0,\infty}$, for all trajectories $X_{t}$ generated by non-autonomous gradient flows and all dimensions $d\in\mathbf{N}$.

A convergence guarantee is given by $\displaystyle f(X_{t})-f_{\star}\leqslant\frac{1}{2\int_{0}^{t}\alpha_{s}ds}\|x_{0}-x_{\star}\|^{2}$.

If $\mu>0$, the Lyapunov function $\mathcal{V}(X_{t},t)=e^{2\mu\int_{0}^{t}\alpha_{s}ds}(f(X_{t})-f_{\star})$ verifies $\displaystyle\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant-2\mu\alpha_{t}\mathcal{V}(X_{t},t)$ for all strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$ ($\mu>0$), all trajectories $X_{t}$ generated by the non-autonomous gradient flow and all dimensions $d\in\mathbf{N}$.

A convergence guarantee is given by $\displaystyle f(X_{t})-f_{\star}\leqslant e^{-2\mu\int_{0}^{t}\alpha_{s}ds}(f(x_{0})-f_{\star})$.

When considering $\alpha_{t}=1$ above, that is $\tau_{t}=t$, we recover exactly the results from Theorem 2.1 and Theorem 2.4.

As mentioned by Orvieto and Lucchi in the stochastic setting, and for accelerated methods by Wibisono et al. , one can thus either work with $Y_{t}$ generated by the non-autonomous gradient flow (18), or with $X_{t}$ generated by the gradient flow. However, the acceleration on $Y_{t}$ is not preserved after discretizing the flow. Applying explicit Euler scheme to a non-autonomous gradient flow (18) with $\alpha_{s}=\alpha>0$ and originating from $f\in\mathcal{F}_{0,L}$, a condition on step sizes $h>0$ arises $0\leqslant h\leqslant\frac{2}{L\alpha}$.

When focusing on continuous-time models for analyzing explicit optimization methods, that only calls for past gradient iterates, we prefer working with the gradient flow (3). However, non-autonomous gradient flows (18) may be useful for analyzing other methods such as proximal methods. More generally, and in the next section, we analyze gradient flows without adjusting the time scale (taking $\alpha_{t}=1$).

3.2 A non-autonomous second-order gradient flows

Nesterov’s accelerated gradient flow reaches an $\frac{1}{t^{2}}$ convergence (see Theorem 2.10) in function values. Considering the family of quadratic Lyapunov functions (14), we study convergence of non-autonomous second-order gradient flows and draw comparison with the accelerated gradient flow. Let $\beta_{t}\geqslant 0$ be a continuous function, and a second-order non-autonomous ODE,

For any two functions $\beta_{t}^{(1)},\beta_{t}^{(2)}\geqslant 0$, there is no time change formula that connects their associated ODE. In other words, they are in the same time-scale. Theorem 2.15 provides an LMI equivalence for analyzing convergence of trajectories $X_{t}$ generated by second-order gradient flows (19) using quadratic Lyapunov functions (14).

Let us consider $\mathcal{V}_{a_{t},P_{t}}$ quadratic Lyapunov functions (14), and non-autonomous second-order gradient flows (19) starting from $x_{0}\in\mathbf{R}^{d}$ where $d\in\mathbf{N}$ is the dimension, and originating from possibly non-strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$ where $\mu\geqslant 0$. The following assertions are equivalent:

The inequality $\displaystyle\frac{d}{dt}\mathcal{V}_{a_{t},P_{t}}(X_{t},t)\leqslant 0$ is satisfied for all functions $f\in\mathcal{F}_{\mu,\infty}$, for $X_{t}$ generated by non-autonomous second-order gradient flows, and all dimensions $d\in\mathbf{N}$.

There exist $\lambda^{(1)}_{t},\lambda^{(2)}_{t}\geqslant 0$ such that,

For a given function $\beta_{t}\geqslant 0$, the LMI remains convex in $\lambda^{(1)}_{t},\lambda^{(2)}_{t}\geqslant 0$. Compared with the LMI derived from Nesterov’s accelerated gradient method in Theorem 2.10, the LMI is parametrized by $\beta_{t}$.

Let us consider non-autonomous second-order gradient flows starting from $x_{0}\in\mathbf{R}^{d}$ (19) where $d\in\mathbf{N}$ is the dimension, and originating from possibly non-strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$ where $\mu\geqslant 0$. The Lyapunov function

If $\mu>0$, $a_{t}=\min(\sqrt{\mu},\frac{2}{3}\beta_{t})$,

If $\mu=0$, $a_{t}=\min((\sqrt{a_{0}}+(\sqrt{p_{0}^{(11)}}/2)t)^{2},\lim_{\epsilon\rightarrow 0,\\ \epsilon>0}a_{\epsilon}e^{\int_{\epsilon}^{t}\frac{2}{3}\beta_{s}ds})$,

verifies $\displaystyle\frac{d}{dt}\mathcal{V}(X_{t},t)\leqslant 0$ for all functions $f\in\mathcal{F}_{\mu,\infty}$, all $X_{t}$ generated the second order gradient flows and all dimensions $d\in\mathbf{N}$.

The proof follows from Theorem 2.15, and is detailed in Appendix A.2.

Non-autonomous second-order gradient flows (19) cannot converge faster than Nesterov’s accelerated gradient flow in function values, that is to say not faster than $\frac{1}{t^{2}}$, when using quadratic Lyapunov functions parametrized by $\beta_{t}\geqslant 0$ from Corollary 2.16. To analyze Nesterov’s accelerated gradient methods using ODEs, Su et al. introduced parametrized second-order gradient flows, that fit the model (19),

When $r\geqslant 3$, the guarantee $f(X_{t})-f_{\star}\leqslant\frac{(r-1)^{2}\|x_{0}-x_{\star}\|^{2}}{2t^{2}}$ holds for all convex functions $f$ and all trajectories $X_{t}$ generated by accelerated gradient flows [37, Theorem 5]. When $r<3$, Attouch et al. [1, Theorem 2.1] proved a convergence bound in $f(X_{t})-f_{\star}=\mathcal{O}(\frac{1}{t^{2r/3}})$, improving the results from [37, Theorem 7] that required additional assumptions on $f$. Using Corollary 2.16, we retrieve a similar bound in function values $f(X_{t})-f_{\star}\leqslant\frac{2\|x_{0}-x_{\star}\|^{2}r^{2}}{9t^{\frac{2r}{3}}}$.

Polynomial convergence can be achieved up to a change of variable, as it was shown by Wisobono et al. with $\tau_{t}=t^{p/2}$ ($\alpha_{t}=\frac{p}{2}t^{p/2-1}$). Nesterov’s accelerated gradient flow ($r=3$) has an ODE $\ddot{X}_{t}+\frac{p+1}{t}\dot{X}_{t}+\frac{p^{2}}{4}t^{p-2}\nabla f(X_{t})=0$ for $p\geqslant 2$. For all convex functions and all trajectories $X_{t}$ generated by the ODE starting from $x_{0}$, a convergence bound in function value is given by $f(X_{t})-f_{\star}\leqslant\frac{\|x_{0}-x_{\star}\|^{2}}{2t^{p}}$.

We have extended the performance estimation approach to continuous-time models, using Lyapunov functions. Given an ODE and a class of functions, we presented a semidefinite formulation equivalent with the existence of a quadratic Lyapunov function. This LMI provides a principled way to generate Lyapunov functions, even when dealing with non-autonomous parametrized ODEs. It turns out only two convex inequalities (in $(X_{t},x_{\star})$ and $(x_{\star},X_{t})$, see Theorem 2.10) are involved in convergence proofs for continuous-time models. For strongly convex functions, we proved numerically worst-case guarantees from Corollary 2.1 and cannot be improved using a specific family of quadratic Lyapunov functions. Even if the connection between discrete and continuous-time does not allow to directly transfer results to the analysis of first-order methods, their convergence analyses are be closely related.

SDEs for modeling SGD

Convergence results for stochastic convex optimization often require additional assumptions on function classes, refined choices of step sizes and averaged iterates. Their analyses raise challenges and more complex proofs in contrast with deterministic methods. Prior works have been concerned with the connection between stochastic methods and stochastic differential equations (SDEs) . This section is devoted to convergence analyses of SDEs using the Lyapunov framework. As for ODEs, verifying a Lyapunov function will be cast as verifying the feasibility of a small-sized LMI.

Let us recall the stochastic gradient descent method (SGD)

where $B_{t}$ is a standard Brownian motion, is an (order-$1$ weak) approximation of SGD ([21, Theorem 1], ). The SDE approximation of SGD allows to take into account the role of fixed step size in the dynamics of SGD (while keeping them small). Under mild assumptions on $f$, Li et al. proved the weak approximation of SGD by this SDE on a finite interval $[0,T]$: there exists $C>0$ such that $\|\mathbf{E}[x_{k}]-\mathbf{E}[X(k\gamma)]\|\leqslant C\gamma$ for $k\in[0,\frac{T}{\gamma}]$. However, this approach is limited since $C$ depends exponentially on $T$. In the literature, matching rates are often obtained extending Lyapunov functions from continuous to discrete time .

The SDE (22) is an approximation of SGD for small step sizes $\gamma$. When taking the step size to zero, the noise term actually disappears, and the limiting ODE of SGD is exactly the gradient flow (3). Similarly, the stochastic Langevin dynamics $x_{k+1}=x_{k}-\gamma\nabla f(x_{k})-\sqrt{\gamma}\xi_{k}$ has the limiting ODE $dX_{t}=-\nabla f(X_{t})dt+\sqrt{2}dB_{t}$, where the step size is not taken into account.

Compared with the gradient flow, SGD does not converge to a stationary point under fixed step sizes . Convergence to a stationary point requires diminishing step sizes such as $\gamma_{k}=\frac{1}{\sqrt{k}}$. Li et al. [21, Section 4.1] (and later Orvieto and Lucchi [26, Section 2.1]) proposed to include this varying learning rate in the dynamics:

where $\gamma$ is the maximum allowed learning rate and $h_{k}\in$ is the varying part. For $h_{t}\geqslant 0$ a continuous function corresponding to discrete $h_{k}$, the SDE is given by

We treat the covariance matrix $\Sigma(X_{t})$ as symmetric, already implied by the notation $\Sigma(X_{t})^{1/2}$, but unstructured with bounded variance $\Sigma(X_{t})\preceq\Sigma$ along any trajectory $X_{t}$ generated by the approximating SDE (23). Compared to ODEs, functions $f\in\mathcal{F}_{0,L}$ to be optimized using SDEs are in addition assumed to be possibly $L$-smooth with $L\in(0,\infty]$, and to be twice differentiable.

We propose to analyze approximating SDEs with averaging techniques, and to include later varying step sizes. Verifying Lyapunov functions thanks to small-sizes LMI, we retrieve convergence results from discrete optimization methods, using appropriate choices of step sizes.

The analysis of the gradient flow in the deterministic case provides Lyapunov functions that are decreasing along trajectories generated by ODEs. The direct transfer of these Lyapunov functions to the stochastic setting is not always suited to the variance term, as detailed below. Under constant step sizes $\gamma>0$, an approximating SDE of SGD is

In SDE theory, a time-differential of a function of a solution to a stochastic process is given by Ito’s Lemma [39, Theorem 4.2].

For $g$ a twice differentiable function, and $X_{t}$ a stochastic process solution to the SDE (22),

When the SDE originates from (possibly non-smooth) strongly convex functions $f\in\mathcal{F}_{\mu,\infty}$, the Lyapunov function from the deterministic setting extends well to SDEs. In the deterministic setting, we have shown in Theorem 2.1 that the function $\mathcal{V}(x,t)=e^{2\mu t}(f(x)-f_{\star})$ is a Lyapunov function in the worst-case. Applying Ito’s formula to this Lyapunov function with $X_{t}$ a solution to the SDE (22), we obtain $\frac{d}{dt}\mathbf{E}\mathcal{V}(X_{t},t)\leqslant\frac{1}{2}e^{2\mu t}\gamma\mathbf{E}{\rm Tr}(\nabla^{2}_{xx}f(X_{t})\Sigma(X_{t}))$. After integrating between and $t$, we have

Additional requirements on $f$, such as smoothness and twice differentiability, are needed for convergence. For example, if $f$ is in addition $L$-smooth with $L<\infty$ and using the bounded covariance assumption $\Sigma(X_{t})\preceq\Sigma$, the variance term is bounded by $\frac{1}{2}L\gamma\text{Tr}(\Sigma)$. Under constant step sizes, the SDE approximating SGD converges to a diffusion, and cannot get to a stationary point $x_{\star}$ in the worst-case, but the extra term is linear in the step size $\gamma$. As in the deterministic case, the forgetting of initial conditions remains of order $O(e^{-2\mu t})$.

1.2 Minimizing convex functions

When $f\in\mathcal{F}_{0,\infty}$, Lyapunov functions induce convergence bounds with a possibly diverging variance term. When considering deterministic gradient flows (3) originating from convex functions, a Lyapunov function followed from Theorem 2.4: $\mathcal{V}(x,t)=t(f(x)-f_{\star})+\frac{1}{2}\|x-x_{\star}\|^{2}$. Assuming $X_{t}$ are solutions $X_{t}$ to SDEs (22) originating from convex functions $f\in\mathcal{F}_{0,\infty}$, we compute the derivative to this Lyapunov function thanks to Ito’s formula $\frac{d}{dt}\mathbf{E}\mathcal{V}(X_{t},t)\leqslant-t\mathbf{E}\|\nabla f(X_{t})\|^{2}+\mathbf{E}\frac{1}{2}\text{Tr}((t\nabla_{x}^{2}f(X_{t})+I)\Sigma(X_{t}))$. Thanks to twice differentiability, bounded covariance $\Sigma(X_{t})\preceq\Sigma$, and assuming in addition $L$-smoothness of $f$, a convergence bound is given by

Taylor and Bach proved a comparable convergence bound for SGD [41, Theorem 5], applying the Lyapunov performance estimation approach under similar assumptions (bounded variance, smoothness of $f$). Tough, we cannot conclude about convergence of SGD in the worst-case without further assumptions. Optimization methods have been developed to ensure the global convergence of SGD to the optimum, among them averaging and diminishing step sizes.

2 Diminishing the step size is a key to success

We study convergence of SDEs with varying step sizes (23). In contrast to the deterministic setting, in which varying step sizes correspond to a time rescaling, the time change plays a direct role on the variance term (explicit formula by Orvieto and Lucchi in [26, Theorem 5]). Our problem can be formulated as follows, where $\mathcal{V}$ are Lyapunov functions,

where $\frac{d}{dt}\mathbf{E}\mathcal{V}(X_{t},t)=\mathbf{E}[\frac{\partial}{\partial t}\mathcal{V}(X_{t},t)+\frac{\partial}{\partial x}\mathcal{V}(X_{t},t)\frac{dX_{t}}{dt}]+\frac{\gamma}{2}\mathbf{E}{\rm Tr}(\frac{\partial^{2}}{\partial x^{2}}\mathcal{V}(X_{t},t)\Sigma(X_{t}))$ is computed with Ito’s formula. The first two terms corresponds exactly to taking the derivative in trajectories generated by ODEs (18) (or the SDEs (23) with $\gamma=0$), and the last term to a variance term. Because of the trace in a second-order derivative and in the covariance matrix $\Sigma(X_{t})$, we do not take this term into account in an LMI reformulation. Instead, we propose to first derive a family of Lyapunov functions using LMIs from the deterministic setting, and then to optimize their parameters so that the variance term converges conveniently.

Let $f\in\mathcal{F}_{0,\infty}$ be twice differentiable functions, and let us consider SDEs with varying step sizes (23) starting from $X_{0}\in\mathbf{R}^{d}$, where $d\in\mathbf{N}$ is the dimension. The quadratic Lyapunov function from the family (11)

with $\dot{a}_{t}^{(1)}=2h_{t}$ verifies $\displaystyle\frac{d}{dt}\mathbf{E}(\mathcal{V}(X_{t},t)\leqslant h_{t}^{2}\text{Tr}((\nabla_{xx}^{2}f(X_{t})a^{(1)}_{t}+\frac{1}{2}I_{d})\Sigma(X_{t}))$ for all twice differentiable functions $f\in\mathcal{F}_{0,\infty}$, all trajectories $X_{t}$ generated by SDEs (23), and all dimensions $d\in\mathbf{N}$. Then, it holds that $\displaystyle\mathbf{E}[f(X_{t})-f_{\star}]\leqslant\frac{\|X_{0}-X_{\star}\|^{2}}{a^{(1)}_{t}}+\frac{\gamma}{2a^{(1)}_{t}}\int_{0}^{t}h_{s}^{2}\text{Tr}((\nabla_{xx}^{2}f(X_{s})a^{(1)}_{s}+\frac{1}{2}I_{d})\Sigma(X_{s}))ds$.

The Lyapunov function $\mathcal{V}$ directly comes from Corollary 2.12, for non-autonomous first-order gradient flows. The bound on $\mathbf{E}[f(X_{t})-f_{\star}]$ is obtained using Ito’s formula on $\mathcal{V}$ along trajectories $X_{t}$ generated by approximating SDEs (23).

The convergence bound from Corollary 3.3 is divided into two terms: a term that forgets the initial conditions and a variance term due to noise. Convergence is mostly controlled by the step size ($\dot{a}^{(1)}_{t}=2h_{t}$). Bach and Moulines [2, Theorem 5] provided a comparable but much more complex analysis for stochastic gradient descent, for a specific family of step sizes. Let us compare our results, assuming $h_{t}=\frac{1}{(t+1)^{\alpha}}$, where $\alpha\geqslant 0$. The Lyapunov function from Corollary 3.3 is given by $\mathcal{V}(X_{t},t)=a^{(1)}_{t}(f(X_{t})-f_{\star})+\frac{1}{2}\|X_{t}-x_{\star}\|^{2}$ with $a^{(1)}_{t}=(t+1)^{1-\alpha}$. The forgetting of the initial condition is thus bounded by $\frac{\|x_{0}-x_{\star}\|^{2}}{(t+1)^{1-\alpha}}$. Provided $\Sigma(X_{t})\preceq\Sigma$, the variance term is bounded by:

Because of the variance term, the Lyapunov converges if and only if $\alpha\geqslant\frac{1}{2}$. In other words, convergence is not guaranteed for constant step sizes ($\alpha=0$). If $\alpha\in(1/2,2/3)$, the convergence in function value is bounded by $O(\frac{1}{t^{2\alpha-1}})$. If $\alpha\in(2/3,1)$, the convergence in function value is bounded by $O(\frac{1}{t^{1-\alpha}})$. As for SGD [2, Theorem 5], the convergence regime changes at $\alpha=\frac{2}{3}$ with a global convergence rate in $\frac{1}{t^{1/3}}$, for which the variance term and the term that forgets the initial conditions converges at the same rate (up to $\text{log}(t)$). It is therefore possible to reach convergence with diminishing step sizes. Other techniques have been developed to improve the trade-off between faster convergence and larger step sizes.

3 Averaging for larger step sizes

Polyak-Ruppert averaging is a standard way to improve convergence of SGD. In the discrete setting, convergence guarantees are considered at an averaged sequence defined by:

Other averaging techniques were later developed, such as primal averaging and averaging with respect to some nonnegative function . Primal averaging is detailed for comparison in Appendix C.

In this section, we analyze convergence properties of SDEs (23) with varying step sizes under Polyak-Ruppert averaging. Taylor and Bach [41, Theorem 6] provided a decreasing Lyapunov function along $(\bar{x}_{k},x_{k})$ (24), and a condition on the step size for convergence to the optimum. An approximating SDE for Polyak-Ruppert averaging is:

with step size $\gamma>0$ that is taken close to zero, and a variable term $h_{t}\in$. We introduce the family of quadratic Lyapunov functions:

where $a^{(1)}_{t},a^{(2)}_{t}\geqslant 0$ are differentiable real functions, and $P_{t}=\begin{pmatrix}p^{(11)}_{t}&p^{(12)}_{t}\\ p^{(12)}_{t}&p^{(22)}_{t}\end{pmatrix}\succeq 0$ has differentiable parameters. Given a Lyapunov function $\mathcal{V}_{a^{(1)}_{t},a^{(2)}_{t},P_{t}}$, and SDEs (25), the performance estimation problem is formulated by:

A possible control of this quantity is presented in Theorem 3.5.

Let $f\in\mathcal{F}_{0,\infty}$ be twice differentiable functions, and SDEs be taken under Polyak-Ruppert averaging (25) starting from $x_{0}\in\mathbf{R}^{d}$, where $d\in\mathbf{N}$ is the dimension, and $\mathcal{V}$ be quadratic Lyapunov functions defined by (26). The following assertions are equivalent,

The inequality $\displaystyle\frac{d}{dt}\mathbf{E}\mathcal{V}(X_{t},\bar{X}_{t},t)\leqslant\frac{1}{2}\text{Tr}((a^{(1)}_{t}\nabla_{xx}f(X_{t})+2p^{(11)}_{t}I_{d})\Sigma(X_{t}))h_{t}^{2}\gamma$ is satisfied for all functions $f\in\mathcal{F}_{0,\infty}$, all trajectories $(X_{t},\bar{X}_{t})$ generated by SDEs under Polyak-Ruppert averaging (25), and all dimensions $d\in\mathbf{N}$.

There exist $\lambda^{(1)}_{t},...,\lambda^{(6)}_{t}\geqslant 0$ such that,

The proof follows the methodology from Section 2.1.1, and using the Gram matrix $G=P^{\top}P$, where $P=(X_{t}-x_{\star},\bar{X}_{t}-x_{\star},g_{t},\bar{g}_{t})$ (see Appendix B.1).

The variance term increases with $a^{(1)}_{t}$, and its convergence requires an additional smoothness assumption on $f$. For this reason, we propose to analyze the convergence based on Lyapunov functions on the averaged sequence only ($\mathcal{V}_{0,a^{(2)}_{t},P_{t}}$).

Let $f\in\mathcal{F}_{0,\infty}$ be twice differentiable functions, SDEs be taken under Polyak-Ruppert averaging (25) starting from $x_{0}\in\mathbf{R}^{d}$ where $d\in\mathbf{N}$ is the dimension. Assuming $\dot{a}^{(2)}_{t}\leqslant\frac{a^{(2)}_{t}}{t}$ and $t\rightarrow\frac{a^{(2)}_{t}}{th_{t}}$ is a non-increasing function, the Lyapunov function

verifies $\frac{d}{dt}\mathbf{E}\mathcal{V}(X_{t},\bar{X}_{t},t)\leqslant\frac{a^{(2)}_{t}}{t}h_{t}\text{Tr}(\Sigma(X_{t}))$ for all twice differentiable functions $f\in\mathcal{F}_{0,\infty}$, all trajectories $X_{t}$ generated by SDEs (25) and all dimensions $d\in\mathbf{N}$. Then, it holds that $\displaystyle\mathbf{E}[f(\bar{X}_{t})-f_{\star}]\leqslant\frac{\|x_{0}-x_{\star}\|^{2}}{2a^{(2)}_{t}}+\frac{\gamma}{2a^{(2)}_{t}}\int_{0}^{t}\frac{a^{(2)}_{s}}{s}h_{s}\text{Tr}(\Sigma(X_{s}))ds$.

The proof follows from Theorem 3.5 (see Appendix B.1).

When $a^{(2)}_{t}=t$ (its maximal possible value), the step size verifies $\dot{h}_{t}\geqslant 0$. The variance term does not diverge if and only if $h_{t}$ is constant. Then, a convergence bound is given by $\mathbf{E}[f(\bar{X}_{t})-f_{\star}]\leqslant\frac{\|x_{0}-x_{\star}\|^{2}}{2t}+\frac{1}{2}\text{Tr}(\sigma)\gamma h$. The decreasing condition on $t\rightarrow\frac{a^{(2)}_{t}}{th_{t}}$ suggests a trade-off between converging and diminishing step size, as obtained without averaging.

Under the assumptions of Corollary 3.7, let us consider a bounded covariance matrix $\Sigma(X_{t})\preceq\Sigma$, a step size $h_{t}=\frac{1}{(t+1)^{\alpha}}$, and $a_{t}^{(2)}=t^{\beta}$ for $\alpha\geqslant 0$ and $0\leqslant\beta\leqslant 1$ some parameters. The decreasing condition imposes $\alpha+\beta\leqslant 1$. A different behavior is expected from $\alpha$ and $\beta$: on the one hand, an ideal step size should be large ($\alpha$ small), and on the other hand, we aim at converging as fast as possible ($\beta$ large). The term that forgets the initial conditions is bounded by $O(\frac{1}{t^{\beta}})$, and the variance term by $O(\frac{1}{t^{\alpha}})$ if $\beta\neq\alpha$. When $\alpha=\beta$, the variance term is in $\frac{\text{log}(t)}{2t^{\beta}}$. Hence, a natural choice is $\alpha=\beta=\frac{1}{2}$, retrieving results from [2, Theorem 4] [41, Table 2] in discrete optimization.