C Éline Moucer
email: ce-line.moucer@inria.fr

Adrien Taylor
email: adrien.taylor@inria.fr

Francis Bach
email: francis.bach@inria.fr

A SYSTEMATIC APPROACH TO LYAPUNOV ANALYSES OF CONTINUOUS-TIME MODELS IN CONVEX OPTIMIZATION

Keywords: Convex optimization, continuous-time models, first-order methods, worst-case analysis, performance estimation, stochastic differential equations, ordinary differential equations

First-order methods are often analyzed via their continuous-time models, where their worst-case convergence properties are usually approached via Lyapunov functions. In this work, we provide a systematic and principled approach to find and verify Lyapunov functions for classes of ordinary and stochastic differential equations. More precisely, we extend the performance estimation framework, originally proposed by Drori and Teboulle [9], to continuous-time models. We retrieve convergence results comparable to those of discrete methods using fewer assumptions and convexity inequalities, and provide new results for stochastic accelerated gradient flows.

1. 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:

(1.1) min

x∈R d f (x),
where f is convex, differentiable, and x ∈ 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 [START_REF] Bottou | The tradeoffs of large scale learning[END_REF]. Gradient descent is a common first-order method, which starts from a point x 0 ∈ R d , its iterates are given by the simple recursion (1.2)

x k+1 = x k -γ∇f (x k ),
where γ > 0 is a step size. Denoting Ẋt = d dt X t , gradient descent with small step sizes γ is directly related to the so-called gradient flow:

(1.3) Ẋt = -∇f (X t ), X 0 = x 0 ∈ R d ,
where the solution X t of the ODE verifies X t k ≈ x k with the identification t k = γk. In numerical integration, gradient descent (1.2) is also known as the Euler explicit scheme for integrating the gradient flow. Recently, Su et al. [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] have interpreted Nesterov's accelerated gradient [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF] in a similar fashion through its continuous version, paving the way to several continuous analyses of accelerated methods [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Wilson | A Lyapunov analysis of accelerated methods in optimization[END_REF][START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF].

Many applications entail some randomness and require a stochastic modeling of the function f , which is often defined in terms of an expectation f (x) = E ξ [ f (x, ξ)]. The function f is the expectation over some random variable ξ, and accounts some random modeling. When ξ is drawn uniformly from a finite set of possible samples (ξ 1 , ..., ξ n ), we have a finite sum f (x) = 1 n n k=1 f (x, ξ k ). As soon as the number of data points n is large, computing the gradient of a finite sum, as it is done in gradientbased methods, is possibly expensive (computing the gradient of each element of the sum, which is possibly very large). Stochastic gradient descent (SGD) provides an alternative with lower computational burden per iteration, by evaluating only the gradient of a single f (•, ξ i k ) per iteration,

x k+1 = x k -γ∇ f (x k , ξ i k ),
where γ > 0 is the step size, ξ i k is drawn uniformly at random in (ξ 1 , ..., ξ n ), and thereby E ξ [∇ f (x k , ξ i k )] = f (x k ) is an unbiased estimate of the full gradient. Li et al. [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF][START_REF] Li | Stochastic modified equations and dynamics of stochastic gradient algorithms i: Mathematical foundations[END_REF] have proven a connection with stochastic differential equations (SDE):

dX t = -∇f (X t )dt + σ(X t )dB t ,
where σ(X t ) is a noise parameter connected to parameters of the method, that was further developed by Shi et al. [START_REF] Shi | On learning rates and schr\" odinger operators[END_REF][START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF]. This connection has raised many questions regarding approximate equivalences between optimization methods and continuoustime 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 [START_REF] Drori | Performance of first-order methods for smooth convex minimization: a novel approach[END_REF] 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.

1.1. 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 [START_REF] Kalman | Control system analysis and design via the "second method" of lyapunov: I-continuous-time systems[END_REF].

Definition 1.1. Given a trajectory X t , we call V : R d → R a Lyapunov function if it is differentiable and satisfies the conditions:

1. V(x) = 0 ⇐⇒ x = x ⋆ , 2. V(X t ) 0, 3. d dt V(X t ) 0. Given an ODE starting from x 0 and a class of functions F , a function V(•) is a Lyapunov function if the inequality d dt V(X t ) 0 is verified for all trajectories X t generated by ODEs originating from functions f ∈ F . There exist similar definitions of Lyapunov functions for discrete optimization methods [START_REF] Taylor | Lyapunov functions for first-order methods: Tight automated convergence guarantees[END_REF][START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF][START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF].

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 d dt V(X t ) -τ V(X t ), where τ 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 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 [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF][START_REF] Aspremont | Acceleration methods[END_REF] and continuous time [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] literature.

1.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 [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], [10, Theorem 4.8]. In the recent work [START_REF] Bansal | Potential-function proofs for gradient methods[END_REF], the authors proposed Lyapunov-based analyses for many first-order methods, for linear and sublinear convergence rates. Continuoustime versions of optimization methods also often involve Lyapunov arguments, such as Nesterov's accelerated gradient flow introduced in [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF], and its high-resolution ODEs for strongly convex functions proposed in [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF], or accelerated mirror descent whose continuous dynamics was analyzed in [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF].

Different techniques were developed to compute suitable Lyapunov functions. The authors of [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] proposed an approach based on Bregman Lagrangian for accelerated methods in potentially non-Euclidean settings, further developed in [START_REF] Wilson | A Lyapunov analysis of accelerated methods in optimization[END_REF]. [START_REF] Diakonikolas | Generalized momentumbased methods: A hamiltonian perspective[END_REF] directly derived Lyapunov functions from Hamiltonian equations describing dynamics of ODEs. Using similar conservation laws in a dilated coordinate system, [START_REF] Suh | Continuous-time analysis of agm via conservation laws in dilated coordinate systems[END_REF] 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 [START_REF] Drori | Performance of first-order methods for smooth convex minimization: a novel approach[END_REF], and led to the notion of performance estimation problems. Their work was followed up in [START_REF] Taylor | Smooth strongly convex interpolation and exact worst-case performance of first-order methods[END_REF][START_REF] Taylor | Exact worst-case performance of first-order methods for composite convex optimization[END_REF] to provide worst-case bounds in a principled way, and extended to the Lyapunov framework [START_REF] Taylor | Lyapunov functions for first-order methods: Tight automated convergence guarantees[END_REF]. A competing strategy inspired by control theory was developed by [START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF][START_REF] Hu | Dissipativity theory for Nesterov's accelerated method[END_REF], 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 [START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF]. Connections between Lyapunov functions obtained via the IQC framework in continuous and discrete-time, were later highlighted by [START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF].

For stochastic differential equations (SDEs), convergence proofs can also be obtained through the Lyapunov approach, together with Ito's calculus. [START_REF] Orvieto | Continuous-time models for stochastic optimization algorithms[END_REF] analyzed both SGD, SAGA [START_REF] Defazio | SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], and SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF], for some well-chosen Lyapunov functions. [START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF][START_REF] Xu | Continuous and discrete-time accelerated stochastic mirror descent for strongly convex functions[END_REF] extended the framework of [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] 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.

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 ma-trix 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.

1.4. Assumptions. Throughout this work, functions to be minimized are convex (see Problem 1.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 : R d → R is convex if for all x, y ∈ R d , and for all λ ∈ [0, 1], f (λx+(1-λ)y) λf (x)+(1-λ)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 [START_REF] Bolte | Characterizations of Lojasiewicz inequalities: Subgradient flows, talweg, convexity[END_REF]Section 3.2]. Then, f is convex if and only if for all x, y ∈ R d , f (y) f (x) + ∇f (x), y -x . 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,

∇f (x) -∇f (y) L x -y . A differentiable function f is µ-strongly convex if for any x, y ∈ R d it satisfies ∇f (x) -∇f (y) µ x -y .
Strong convexity ensures the function is not too flat, and the unicity of the minimizer x ⋆ . We denote by F µ,L the family of L-smooth µ-strongly closed convex proper functions, with 0 µ L +∞. Weaker assumptions than strong convexity are also encountered in the literature for analyzing gradient algorithms, and leads to similar convergence guarantees. Among them, Lojasiewicz [START_REF] Lojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF] introduced the Lojasiewicz inequality, under which Polyak [START_REF] Polyak | Gradient methods for minimizing functionals[END_REF]Theorem 4] showed linear convergence of gradient descent. Other relaxed versions of strong convexity followed [START_REF] Bolte | From error bounds to the complexity of first-order descent methods for convex functions[END_REF][START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF].

2.

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.

Gradient flow. We consider the gradient flow

Ẋt = -∇f (X t ), X 0 = x 0 ∈ R d .
Without further assumptions, the function f is decreasing along the trajectory X t solution to the gradient flow. The Lyapunov function

V(X t ) = f (X t ) has indeed a negative time-derivative d dt V(X t ) = Ẋt ⊤ ∇f (X t ) = -∇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.

2.1.1. Minimizing strongly convex functions. Let f be µ-strongly convex (µ > 0), admitting thus a unique minimizer x ⋆ such that f (x ⋆ ) = f ⋆ , and consider the Problem (1.1). Thanks to strong convexity, it is possible to prove linear convergence of the gradient flow to its stationary point.

f ∈ F µ,∞ , (2.1) f (X t ) -f ⋆ e -2µt (f (x 0 ) -f ⋆ ).
This convergence guarantee follows directly from the time-derivative of the Lya-

punov function V(X t ) = f (X t ) -f ⋆ , together with strong convexity (or Lojasiewicz inequality): d dt V(X t ) = Ẋt ⊤ ∇f (X t ) = -∇f (X t ) 2 -2µ(f (X t )-f ⋆ ) -2µ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:

(2.2) V a,c (X t ) = a • (f (X t ) -f ⋆ ) + c • X t -x ⋆ 2 ,
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 (2.1). Given a Lyapunov function V a,c , the idea is to find the smallest value of τ such that the condition

(2.3) d dt V a,c (X t ) -τ V a,c (X t ),
holds for any functions f ∈ F µ,∞ and any trajectory X t generated by the gradient flow. After integrating, a convergence guarantee in function values is given by V a,c (X t ) e -τ t V a,c (X 0 ). Given a certain Lyapunov function V a,c and a time t, we get that the largest acceptable τ is a solution to:

(2.4)

-τ = max Xt∈R d ,d∈N,f ∈Fµ,∞ d dt V a,c (X t ), subject to V a,c (X t ) = 1, Ẋt = -∇f (X t ).
This minimization problem is invariant in t. It is established in [START_REF] Taylor | Smooth strongly convex interpolation and exact worst-case performance of first-order methods[END_REF][START_REF] Drori | Performance of first-order methods for smooth convex minimization: a novel approach[END_REF] 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 ∈ F µ,∞ , the maximization problem (2.4) is infinite-dimensional. Recall that a differentiable function f ∈ F µ,∞ verifies for all points x, y

∈ R d , f (x) -f (y) -∇f (y), x -y µ 2 x -y 2 .
Introducing alternate variables f t , f ⋆ , g t and g ⋆ (informally:

f t = f (X t ), f ⋆ = f (x ⋆ ), g t = ∇f (X t ) and g ⋆ = ∇f (x ⋆ ) = 0), it holds that (2.5) -τ = max Xt∈R d ,d∈N d dt V a,c (X t ) subject to V a,c (X t ) = 1, Ẋt = -g t , f i -f j -g j , X i -X j µ 2 X i -X j 2 , ∀i, j = t, ⋆.
The fact (2.5) produces an upper bound on -τ directly follows from the fact that any sampled strongly convex function satisfy these inequalities at the sampled points X t and x ⋆ . Thereby, any feasible point to (2.4) corresponds to a feasible for (2.5) with the same objective value. In the other direction, [45, Corollary 2] (which provides a constructive way to obtain some f ∈ F µ,∞ that interpolates the triplets (X i , g i , f i ) t,⋆ ) ensures that any feasible point to (2.5) can be translated to a feasible point to (2.4) with also the same objective value, thereby reaching the equivalence between formulations (2.4) and (2.5).

In a second stage, we introduce

G = X t -x ⋆ 2 X t -x ⋆ , g t X t -x ⋆ , g t g t 2
0 a Gram matrix and a vector F = [f t , f ⋆ ], thereby obtaining a semidefinite reformulation:

(2.6)

-τ = max G 0,F ∈R 2 ,d∈N b ⊤ 0 F + Tr(A 0 G) subject to b ⊤ 1 F + Tr(A 1 G) 0, b ⊤ 2 F + Tr(A 2 G) 0.
where

A 0 = c • τ -c -c -a , A 1 = -µ/2 1/2 1/2 0 , A 2 = -µ/2 0 0 0 , b 0 = a • τ [1, -1] ⊤ b 1 = [-1, 1] ⊤ and b 2 = [1, -1] ⊤ .
Those developments allow arriving to the following worst-case results.

Theorem 2.1. Let V a,c be a quadratic Lyapunov function (2.2) with a, c 0, τ 0, and functions f ∈ F µ,∞ be strongly convex with parameter µ > 0. We consider gradient flows (1.3) originating from functions f ∈ F µ,∞ and starting from an initial point x 0 ∈ R d , where d ∈ N is the dimension. The following assertions are equivalent:

• The inequality d dt V a,c (X t ) -τ V a,c (X t ) is satisfied for all f ∈ F µ,∞
, all trajectories X t solutions to the gradient flow and all dimensions d ∈ N.

• There exist λ 1 , λ 2 0 such that

(2.7) S = τ c -µ 2 (λ 1 + λ 2 ) -c + λ1 2 -c + λ1 2 -a 0, τ a = λ 1 -λ 2 .
Proof. The procedure to obtain the LMI is fully detailed in Appendix A.1, and consists in taking the standard Lagrangian dual of the SDP (2.6).

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 V a,c to decrease at a specific rate τ 0 for all functions f ∈ F µ,∞ . The infeasibility of the LMI allows to conclude about the existence of a Lyapunov function such that the rate τ 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 τ , the LMI is indeed jointly convex in λ 1 , λ 2 , a, c. Finally, thanks to linearity of the feasibility problem in τ , 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 1a, 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µ from Scieur et al. [START_REF] Scieur | Integration methods and accelerated optimization algorithms[END_REF]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 1b and the method in [START_REF] Taylor | Interpolation and performance estimation of first-order methods for convex optimization[END_REF]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.

2.1.2. Minimizing convex functions. In the case where f ∈ F 0,∞ , 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

V(X t , t) = t(f (X t )-f ⋆ )+ 1 2 X t -x ⋆ 2 from [37, p. 7] ver- ifies d
dt V(X t , t) 0 for any functions f ∈ F 0,∞ , and any trajectories X t generated by the gradient flow (1.3) (proof :

d dt V(X t ) = t ∇f (X t ), Ẋt +f (X t )-f ⋆ + Ẋt , X t -x ⋆ = -t ∇f (X t ) 2 + f (X t ) -f ⋆ -∇f (X t ), X t -x ⋆
-t ∇f (X t ) 2 using convexity). After integrating, we recover a convergence bound in function values from the literature [37, p.7], [START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF]Section 6.3

.1], f (X t ) -f ⋆ x 0 -x ⋆ 2 2t .
We consider the family of quadratic Lyapunov functions:

(2.8) V at,ct (X t , t) = a t (f (X t ) -f ⋆ ) + c t X t -x ⋆ 2 ,
where a t , c t are differentiable functions from R + to R + . When the Lyapunov function is decreasing along the trajectory X t , that is d dt V(X t ) 0, a convergence guarantee in function values is given by

f (X t ) -f ⋆ V(x 0 , 0) a t = a 0 (f (x 0 ) -f ⋆ ) + c 0 x 0 -x ⋆ 2 a t .
Looking for a worst-case guarantee with a Lyapunov approach can be cast as:

0 max Xt∈R d ,d∈N,f ∈F0,∞ d dt V at,ct subject to Ẋt = -∇f (X t ).
Remark 2.2. The strongly convex case as defined above is a particular case of the convex one, using a specific Lyapunov function Φ(•), such that V(X t , t) = e τ t Φ(X t ) where Φ(

X t ) = a • (f (X t ) -f ⋆ ) + c • X t -x ⋆ 2 . Then, d dt V(X t , t) 0 is equivalent to d dt Φ(X t ) -τ Φ(X t ).
Theorem 2.3. Let V at,ct defined in (2.8) be a quadratic Lyapunov function for a t , c t nonnegative differentiable functions and functions f ∈ F 0,∞ be convex. Given the gradient flow (1.3) originating from convex functions f ∈ F 0,∞ and starting from x 0 ∈ R d , where d ∈ N is the dimension. The following assertions are equivalent:

• The inequality d dt V at,ct (X t , t) 0 is satisfied for all functions f ∈ F 0,∞ ,
for all trajectories X t generated by gradient flows, all functions f ∈ F µ,∞ and for all dimensions d ∈ N.

• There exist λ

(1) t , λ (2) t 0 such that S = ċt -c t + λ (1) t 2 -c t + λ (1) t 2 -a t 0, ȧt = λ (1) t -λ (2)
t .

If these assertions are satisfied, c t is decreasing, and a t tc t .

Proof. The LMI is obtained following the previous methodology (Appendix A.1), and the inequality a t tc t directly from the LMI.

Choosing c t = 1 2 and a t = t for the Lyapunov function parameters, and λ

(1) t = 1, λ (2) t = 0, we retrieve the Lyapunov function V(x, t) = t(f (x) -f ⋆ ) + 1 2 x -x ⋆ 2 from [37, p. 7], that verifies f (X t ) -f ⋆ 1 2t x 0 -x ⋆
2 for all convex functions f ∈ F 0,∞ , for all trajectories X t generated by the gradient flow, and all dimensions d ∈ N.

The differential LMI equivalence from Theorem 2.3 is jointly convex in λ (1) t , λ (2)
t , c t , a t , ȧt , ċ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 c t t together with ċt 0 implies that the convergence cannot be faster than 1 t . 2.2. Accelerated gradient flows. A major improvement to gradient descent dates back to Nesterov in 1983 [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], with an accelerated gradient method (AGM), (2.9)

x k+1 = y k -γ∇f (y k ), y k+1 = x k+1 + α k (x k+1 -x k ),
where γ, α 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 [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] with the heavy-ball method, starting from x 0 , x 1 ∈ R d , (2.10)

x k+2 = x k+1 + α k (x k+1 -x k ) -γ∇f (x k+1 ),
where α 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 γ, these schemes happen to be closely related to second-order differential equations, for β t 0 a function that depends on α k ,

Ẍt + β t Ẋ + ∇f (X t ) = 0.
Recently, accelerated gradient methods have been analyzed using second-order differential equations [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF][START_REF] Wilson | A Lyapunov analysis of accelerated methods in optimization[END_REF][START_REF] Hu | Dissipativity theory for Nesterov's accelerated method[END_REF]. 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 [START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF][START_REF] Shi | On learning rates and schr\" odinger operators[END_REF][START_REF] Wilson | A Lyapunov analysis of accelerated methods in optimization[END_REF]. Taking integration theory's point of view, Scieur et al. [START_REF] Scieur | Integration methods and accelerated optimization algorithms[END_REF] 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 V at,Pt be the family of quadratic Lyapunov functions for second-order gradient flows, 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

(2.11) V at,Pt (X t , t) = a t (f (X t ) -f ⋆ ) + X t -X ⋆ Ẋt ⊤ P t X t -
(X t ) -f ⋆ V(x0)
at .

2.2.1. Minimizing strongly convex functions. Let f ∈ F µ,L , with strong convexity parameter µ > 0, and scheme parameters be defined by γ = 1

L and α =

1- √ µγ 1+
√ µγ in Nesterov's accelerated gradient methods (2.9). When reducing the step size γ, the continuous-time limit of y k in (2.9) is exactly the Polyak damped oscillator [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF],

(2.12) Ẍt + 2 √ µ Ẋt + ∇f (X t ) = 0,
where t = γk, as it has already been highlighted in previous works [START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF][START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF][START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF]. This ODE is also the limit of the heavy-ball method (2.10). Shi et al. [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF] proved a convergence guarantee in f

(X t ) -f ⋆ = O(e - √ µt 4
) using a Lyapunov-based approach. This bound was improved to f (X t )-f ⋆ = O(e -√ µt ), by Wilson et al. [47,Appendix B] using the Bregman-Lagrangian approach, and by Sanz-Serna and Zygalakis using the IQC framework [START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF][START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF]. We compute linear convergence guarantees using quadratic Lyapunov functions with constant parameters (2.11).

Theorem 2.4. Let V a,P be a quadratic Lyapunov function (2.11), where a 0, and P a symmetric matrix. Let τ 0, µ 0, d ∈ N be the dimension, and the Polyak damped oscillator (2.12) be starting from x 0 ∈ R d , and originating from strongly convex functions f ∈ F µ,∞ . The following assertions are equivalent:

• The inequality d dt V a,P (X t ) -τ V a,P (X t )
is satisfied for all functions f ∈ F µ,∞ , all trajectories X t generated by Polyak damped oscillator, and all di-

mensions d ∈ R d . • There exist λ 1 , λ 2 , ν 1 , ν 2 0, such that (2.13)   -µ 2 (λ 1 + λ 2 ) + τ p 11 p 11 -2 √ µp 12 + τ p 12 -p 12 + λ1 2 * 2(p 12 -2 √ µp 22 ) + τ p 22 -p 22 + a 2 * * 0   0, τ a = λ 1 -λ 2 , P 0 0 0 +   µ 2 (ν 1 + ν 2 ) 0 -ν1 2 0 0 0 -ν1 2 0 0   0, a = ν 2 -ν 1 .
Proof. The LMI equivalence is obtained introducing the Gram matrix G = P ⊤ P , where P = ( Ẋt , X t -x ⋆ , 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 λ 1 , λ 2 , ν 1 , ν 2 0 and in Lyapunov parameters a, P . Given parameters a, P and a rate τ , it is a verification tool for V a,P to be a Lyapunov function decreasing at rate τ . As for gradient flows, we can perform a bisection search over τ to find the fastest linear convergence rate that can be verified using quadratic Lyapunov functions (see Figure 2a). It provides a numerical tool for choosing Lyapunov parameters (Figure 2b) for which the worst-case linear convergence rate is achieved. Corollary 2.5. Let us consider the Polyak damped oscillator (2.12) starting from x 0 ∈ R d where d ∈ N is the dimension, and originating from strongly convex functions f ∈ F µ,∞ where µ > 0. The Lyapunov function This class of quadratic Lyapunov functions is inspired by [START_REF] Taylor | Lyapunov functions for first-order methods: Tight automated convergence guarantees[END_REF] in discrete time, where a stricter positivity condition on P 0 hindered proving tight convergence of Nesterov's accelerated gradient. Similarly, the Lyapunov function from Corollary 2.12 is defined by

V(X t ) = f (X t ) -f ⋆ + X t -X ⋆ Ẋt ⊤ 4/9µ 2/3 √ µ 2/3 √ µ 1/2 X t -X ⋆ Ẋt , verifies d dt V(X t ) -4/3 √ µV(X t ) for all functions f ∈ F µ,∞ ,
P = 4/9µ 2/3 √ µ 2/3 √ µ 1/2
, 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 [START_REF] Serna | The connections between lyapunov functions for some optimization algorithms and differential equations[END_REF][START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF]Theorem 4.3], Fig. 2: Comparison between numerical worst-case rate guarantee obtained numerically with PEP, and with references. Trajectories X t are generated by the damped gradient flow (2.12) for µ-strongly convex functions.

V(X t ) = f (X t ) -f ⋆ + 1 2 X t -X ⋆ Ẋt ⊤ µ √ µ √ µ 1 X t -X ⋆ Ẋt that verifies d dt V(X t ) - √ µV(X t ) for all functions f ∈ F µ,∞ ,
= √ µ, λ 1 = √ µ, λ 2 = 0, ν 1 = 0, ν 2 = a =
Figure 2b provides a numerical help for computing the Lyapunov function from Corollary 2.5. Figure 2a together with Corollary 2.5 allows to conclude that this bound cannot be improved when changing the Lyapunov function among the class of quadratic functions (2.11).

Minimizing convex functions.

As for gradient flow, rates are sublinear when the accelerated gradient flow originates from convex functions. Let f ∈ F 0,L , the step size be γ 1 L , and α = k-1 k+2 be the scheme parameter in Nesterov's accelerated method. Su et al. [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]Section 2] proved the connection between the first-order scheme and a second order ODE known as accelerated gradient flow (AGF):

(2.14) Ẍt + 3 t Ẋt + ∇f (X t ) = 0.
Su et al. [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]Theorem 3] proved the following inequality is verified for all functions f and all trajectories X t generated by AGF,

f (X t ) -f ⋆ ≤ 2 x 0 -x ⋆ 2 t 2 .
Their proof exhibits a Lyapunov function

V(X t , t) = t 2 (f (X t ) -f ⋆ ) + (X t -x ⋆ ) + t 2t(f (X t ) -f ⋆ ) + t 2 Ẋt , ∇f (X t ) + 2 X t -x ⋆ + t 2 Ẋt , 3 Ẋt + t Ẍt = 2t(f (X t ) -f ⋆ ) - ∇f (X t ), X t -x ⋆ 0 by convexity of f ).
The following theorem provides a systematic condition for a quadratic function V (2.11) to be a Lyapunov function for AGF (2.14).

Theorem 2.6. Let V at,Pt be a quadratic Lyapunov function (2.11), where a t 0 is a differentiable function, and P t 0 with differentiable parameters. Given the accelerated gradient flow (2.14) starting from x 0 ∈ R d where d ∈ N is the dimension, and originating from convex functions f ∈ F 0,∞ , the following assertions are equivalent:

• The inequality d dt V at,Pt (X t , t) 0 is satisfied for all functions f ∈ F 0,∞ ,
+ at 2 * * 0    0, ȧt = λ (1) t -λ (2)
t .

Proof. The proof follows those from Theorem 2.4 and 2.1.

The LMI from Theorem 2.6 is differentiable and convex in λ

(1) t , λ (2)
t . Even when studying second-order gradient flows, the number of interpolation inequalities (or of dual values λ

(i) t )
= 0, we prove that

that the Lyapunov function V(X t , t) = t 2 (f (X t ) -f ⋆ ) + 2 (X t -x ⋆ ) + t 2 Ẋ 2 verifies d dt
V(X t , t) 0, for all functions f ∈ F 0,∞ , for all trajectories X t generated accelerated gradient flows, and all dimensions d ∈ 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.

2.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.

2.3.1.

Non-autonomous first-order gradient flow. Let f be a convex function, and X t be the solution to the non-autonomous first-order gradient flow,

(2.15) Ẋt = -α t ∇f (X t ), X 0 = x 0 ∈ R d ,
where α t 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 α t . A change of variable connects this ODE to the gradient flow (1.3), for which α t = 1. Let Y t be the solution to the gradient flow, and τ t = t 0 α s ds be a time change variable. Then, the variable

X t = Y τt verifies Ẋt = d dt Y τt = α t Ẏτt = -α t ∇f (Y τt ) = -α t ∇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.3.

Corollary 2.7. Let us consider non-autonomous gradient flows (2.15) starting from x 0 ∈ R d where d ∈ N is the dimension, and originating from possibly strongly convex functions f ∈ F µ,∞ , where µ 0.

• If µ = 0, the Lyapunov function V(X t , t) = t 0 α s ds (f (X t ) -f ⋆ ) + 1 2 X t - x ⋆ 2 verifies d dt V(X t , t)
0 for all functions f ∈ F 0,∞ , for all trajectories X t generated by non-autonomous gradient flows and all dimensions d ∈ N.

A convergence guarantee is given by f

(X t ) -f ⋆ 1 2 t 0 α s ds x 0 -x ⋆ 2 . • If µ > 0, the Lyapunov function V(X t , t) = e 2µ t 0 αsds (f (X t ) -f ⋆ ) verifies d dt V(X t , t) -2µα t V(X t , t) for all strongly convex functions f ∈ F µ,∞ (µ > 0)
, all trajectories X t generated by the non-autonomous gradient flow and all dimensions d ∈ N.

A convergence guarantee is given by f (X t ) -f ⋆ e -2µ t 0 αsds (f (x 0 ) -f ⋆ ). Remark 2.8. When considering α t = 1 above, that is τ t = t, we recover exactly the results from Theorem 2.1 and Theorem 2.3.

As mentioned by Orvieto and Lucchi [START_REF] Orvieto | Continuous-time models for stochastic optimization algorithms[END_REF] in the stochastic setting, and for accelerated methods by Wibisono et al. [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF], one can thus either work with Y t generated by the non-autonomous gradient flow (2.15), 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 (2.15) with α s = α > 0 and originating from f ∈ F 0,L , a condition on step sizes h > 0 arises 0 h 2 Lα . 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 (1.3). However, non-autonomous gradient flows (2.15) 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 α t = 1).

2.3.2.

A non-autonomous second-order gradient flows. Nesterov's accelerated gradient flow reaches an (2) t 0, there is no time change formula that connects their associated ODE. In other words, they are in the same time-scale. Theorem 2.10 provides an LMI equivalence for analyzing convergence of trajectories X t generated by second-order gradient flows (2.16) using quadratic Lyapunov functions (2.11).

Theorem 2.10. Let us consider V at,Pt quadratic Lyapunov functions (2.11), and non-autonomous second-order gradient flows (2.16) starting from x 0 ∈ R d where d ∈ N is the dimension, and originating from possibly non-strongly convex functions f ∈ F µ,∞ where µ 0. The following assertions are equivalent:

• The inequality d dt V at,Pt (X t , t) 0 is satisfied for all functions f ∈ F µ,∞ , for X t generated by non-autonomous second-order gradient flows, and all dimensions d ∈ N.

• There exist λ

(1) t , λ (2) t 0 such that,    -µ 2 (λ (1)
+ at 2 * * 0    0, ȧt = λ (1) t -λ (2)
t .

For a given function β t 0, the LMI remains convex in λ

(1)

t , λ (2)
t 0. Compared with the LMI derived from Nesterov's accelerated gradient method in Theorem 2.6, the LMI is parametrized by β t .

Corollary 2.11. Let us consider non-autonomous second-order gradient flows starting from x 0 ∈ R d (2.16) where d ∈ N is the dimension, and originating from possibly non-strongly convex functions f ∈ F µ,∞ where µ 0. The Lyapunov function

V(X t , t) = a t (f (X t ) -f ⋆ ) + 1 2a t a t Ẋt + ȧt (X t -x ⋆ ) 2 ,
with a t defined by:

• If µ > 0, a t = min( √ µ, 2 3 β t ), • If µ = 0, a t = min(( √ a 0 + ( p (11) 0
/2)t) 2 , lim ǫ→0,ǫ>0 a ǫ e t ǫ 2 3 βsds ), verifies d dt V(X t , t) 0 for all functions f ∈ F µ,∞ , all X t generated the second order gradient flows and all dimensions d ∈ N.

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

Non-autonomous second-order gradient flows (2.16) cannot converge faster than Nesterov's accelerated gradient flow in function values, that is to say not faster than 1 t 2 , when using quadratic Lyapunov functions parametrized by β t 0 from Corollary 2.11. To analyze Nesterov's accelerated gradient methods using ODEs, Su et al. [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] introduced parametrized second-order gradient flows, that fit the model (2.16), (2.17)

Ẍt + r t Ẋt + ∇f (X t ) = 0. When r 3, the guarantee f (X t ) -f ⋆ (r-1) 2 x0-x⋆ 2 2t 2
holds for all convex functions f and all trajectories X t generated by accelerated gradient flows [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]Theorem 5]

. When r < 3, Attouch et al. [1, Theorem 2.1] proved a convergence bound in f (X t ) - f ⋆ = O( 1 t 2r/3
), improving the results from [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]Theorem 7] that required additional assumptions on f . Using Corollary 2.11, we retrieve a similar bound in function values

f (X t ) -f ⋆ 2 x0-x⋆ 2 r 2 9t 2r 3
. Remark 2.12. Polynomial convergence can be achieved up to a change of variable, as it was shown by Wisobono et al. [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] with τ t = t p/2 (α t = p 2 t p/2-1 ). Nesterov's accelerated gradient flow (r = 3) has an ODE Ẍt + p+1 t Ẋt + p 2 4 t p-2 ∇f (X t ) = 0 for p 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 ⋆ x0-x⋆ 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 ⋆ ) and (x ⋆ , X t ), see Theorem 2.6) 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.

3. 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) [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF][START_REF] Orvieto | Continuous-time models for stochastic optimization algorithms[END_REF][START_REF] Xu | Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms[END_REF]. 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) (3.1)

x k+1 = x k -γ∇ f (x k , ξ i k ),
where γ > 0 is the step size, ξ i k are uniformly drawn in (ξ 1 , ..., ξ n ), and where ∇ f (x k , ξ i k ) is an unbiased estimate of full gradient ∇f (x k ). Li et al. [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF] introduced stochastic modified equations (SMEs) to model SGD, rewriting it as

x k+1 = x k -γ f (x k , ξ i k ) + √ hV k (x k ),
where

V k (x) = √ γ(∇f (x k )-∇ f (x k , ξ i k )
has zero mean and a covariance matrix equal

to γΣ(x k ) = γ( n i=1 (∇f (x k ) -∇ f (x k , ξ i k ))(∇f (x k ) -∇ f (x k , ξ i k )) ⊤ ). The corresponding SDE is given by (3.2) dX t = -∇f (X t )dt + (γΣ(X t )) 1/2 dB t ,
where B t is a standard Brownian motion, is an (order-1 weak) approximation of SGD ([21, Theorem 1], [START_REF] Li | Stochastic modified equations and dynamics of stochastic gradient algorithms i: Mathematical foundations[END_REF]). 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 E[x k ] -E[X(kγ)] Cγ for k ∈ [0, T γ ].
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 [START_REF] Shi | Acceleration via symplectic discretization of high-resolution differential equations[END_REF][START_REF] Orvieto | Continuous-time models for stochastic optimization algorithms[END_REF][START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF].

Remark 3.1. The SDE (3.2) is an approximation of SGD for small step sizes γ. When taking the step size to zero, the noise term actually disappears, and the limiting ODE of SGD is exactly the gradient flow (1.3). Similarly, the stochastic Langevin dynamics x k+1 = x k -γ∇f (x k ) -√ γξ k has the limiting ODE dX t = -∇f (X t )dt + √ 2dB 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 [START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF][START_REF] Taylor | Stochastic first-order methods: non-asymptotic and computer-aided analyses via potential functions[END_REF]. Convergence to a stationary point requires diminishing step sizes such as γ k = 1

x k+1 = x k -γh k ∇f (x k ),
where γ is the maximum allowed learning rate and h k ∈ [0, 1] is the varying part. For h t 0 a continuous function corresponding to discrete h k , the SDE is given by

(3.3) dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t .
We treat the covariance matrix Σ(X t ) as symmetric, already implied by the notation Σ(X t ) 1/2 , but unstructured with bounded variance Σ(X t ) Σ along any trajectory X t generated by the approximating SDE (3.3). Compared to ODEs, functions f ∈ F 0,L to be optimized using SDEs are in addition assumed to be possibly L-smooth with L ∈ (0, ∞], 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.

Lyapunov functions do not always extend to the stochastic setting.

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 γ > 0, an approximating SDE of SGD is

dX t = -∇f (X t )dt + (γΣ(X t )) 1/2 dB t .
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]. Lemma 3.2 (Ito's Lemma). For g a twice differentiable function, and X t a stochastic process solution to the SDE (3.2),

dg(X t , t) = ∂ ∂t g(X t , t)dt + ∂ ∂x g(X t , t)dX t + 1 2 γTr( ∂ 2 ∂x 2 g(X t , t)Σ(X t ))dt.
3.1.1. Minimizing strongly convex functions. When the SDE originates from (possibly non-smooth) strongly convex functions f ∈ F µ,∞ , 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 V(x, t) = e 2µt (f (x)f ⋆ ) is a Lyapunov function in the worst-case. Applying Ito's formula to this Lyapunov function with X t a solution to the SDE (3.2), we obtain d dt EV(X t , t)

1 2 e 2µt γETr(∇ 2 xx f (X t )Σ(X t )). After integrating between 0 and t, we have

E(f (X t ) -f ⋆ ) e -2µt f (x 0 ) -f ⋆ + 1 2 γ t 0 e 2µs ETr(∇ 2 xx f (X s )Σ(X s ))ds .
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 < ∞ and using the bounded covariance assumption Σ(X t ) Σ, the variance term is bounded by 1 2 LγTr(Σ). Under constant step sizes, the SDE approximating SGD converges to a diffusion, and cannot get to a stationary point x ⋆ in the worst-case, but the extra term is linear in the step size γ. As in the deterministic case, the forgetting of initial conditions remains of order O(e -2µt ).

3.1.2. Minimizing convex functions. When f ∈ F 0,∞ , Lyapunov functions induce convergence bounds with a possibly diverging variance term. When considering deterministic gradient flows (1.3) originating from convex functions, a Lyapunov function followed from Theorem 2.3:

V(x, t) = t(f (x) -f ⋆ ) + 1 2 x -x ⋆ 2
. Assuming X t are solutions X t to SDEs (3.2) originating from convex functions f ∈ F 0,∞ , we compute the derivative to this Lyapunov function thanks to Ito's formula

d dt EV(X t , t) -tE ∇f (X t ) 2 +E 1 2 Tr((t∇ 2 x f (X t )+I)Σ(X t ))
. Thanks to twice differentiability, bounded covariance Σ(X t ) Σ, and assuming in addition L-smoothness of f , a convergence bound is given by

E(f (X t ) -f ⋆ ) x 0 -x ⋆ 2 t + 1 2 (L t 2 + 1)Tr(Σ).
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 [START_REF] Ruppert | Efficient estimations from a slowly convergent Robbins Monroe process[END_REF] and diminishing step sizes.

3.2.

Diminishing the step size is a key to success. We study convergence of SDEs with varying step sizes (3.3). 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 [START_REF] Orvieto | Continuous-time models for stochastic optimization algorithms[END_REF]Theorem 5]). Our problem can be formulated as follows, where V are Lyapunov functions, max

Xt∈R d ,d∈N, f ∈F0,∞ d dt EV(X t , t), subject to dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , where d dt EV(X t , t) = E[ ∂ ∂t V(X t , t) + ∂ ∂x V(X t , t) dXt dt ] + γ 2 ETr( ∂ 2 ∂x 2 V(X t , t)Σ(X t )
) is computed with Ito's formula. The first two terms corresponds exactly to taking the derivative in trajectories generated by ODEs (2.15) (or the SDEs (3.3) with γ = 0), and the last term to a variance term. Because of the trace in a second-order derivative and in the covariance matrix Σ(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.

Corollary 3.3. Let f ∈ F 0,∞ be twice differentiable functions, and let us consider SDEs with varying step sizes (3.3) starting from X 0 ∈ R d , where d ∈ N is the dimension. The quadratic Lyapunov function from the family (2.8)

V(X t , t) = a (1) t (f (X t ) -f ⋆ ) + 1 2 X t -x ⋆ 2 , with ȧ (1)
E[f (X t ) -f ⋆ ] X 0 -X ⋆ 2 a (1) t + γ 2a (1) t t 0 h 2 s Tr((∇ 2 xx f (X s )a (1) s + 1 2 I d )Σ(X s ))ds.
Proof. The Lyapunov function V directly comes from Corollary 2.7, for nonautonomous first-order gradient flows. The bound on E[f (X t ) -f ⋆ ] is obtained using Ito's formula on V along trajectories X t generated by approximating SDEs (3.3).

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 ( ȧ(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 = 1 (t+1) α , where α 0. The Lyapunov function from Corollary 3.3 is given by V(X t , t) = a

(1)

t (f (X t ) -f ⋆ ) + 1 2 X t -x ⋆ 2 with a (1) t = (t + 1)
   γTr(Σ) (t+1) 1-α L (t+1) 2-3α -1 2-3α + 1 2 (t+1) 1-2α -1 1-2α if 0 α < 1, γTr(Σ) log(t) L t 0 log(s+1) (s+1) 2 ds + 1 2 (1 -1 t+1 ) if α = 1.
Because of the variance term, the Lyapunov converges if and only if α 1 2 . In other words, convergence is not guaranteed for constant step sizes (α = 0). If α ∈ (1/2, 2/3), the convergence in function value is bounded by O( 1 t 2α-1 ). If α ∈ (2/3, 1), the convergence in function value is bounded by O( 1 t 1-α ). As for SGD [2, Theorem 5], the convergence regime changes at α = 2 3 with a global convergence rate in 1 t 1/3 , for which the variance term and the term that forgets the initial conditions converges at the same rate (up to log(t)). It is therefore possible to reach convergence with diminishing step sizes. Other techniques have been developed to improve the tradeoff between faster convergence and larger step sizes.

3.3.

Averaging for larger step sizes. Polyak-Ruppert averaging [START_REF] Ruppert | Efficient estimations from a slowly convergent Robbins Monroe process[END_REF][START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF] is a standard way to improve convergence of SGD. In the discrete setting, convergence guarantees are considered at an averaged sequence defined by: (3.4)

x k+1 = x k -γ∇f i k (x k ), xk = 1 k k i=1 x i .
Other averaging techniques were later developed, such as primal averaging [START_REF] Tao | Primal averaging: A new gradient evaluation step to attain the optimal individual convergence[END_REF] and averaging with respect to some nonnegative function [START_REF] Roux | Anytime tail averaging[END_REF]. Primal averaging is detailed for comparison in Appendix C.

Polyak-Ruppert averaging.

In this section, we analyze convergence properties of SDEs (3.3) with varying step sizes under Polyak-Ruppert averaging. Taylor and Bach [START_REF] Taylor | Stochastic first-order methods: non-asymptotic and computer-aided analyses via potential functions[END_REF]Theorem 6] provided a decreasing Lyapunov function along (x k , x k ) (3.4), and a condition on the step size for convergence to the optimum. An approximating SDE for Polyak-Ruppert averaging is:

(3.5) dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , d Xt = X t -Xt t dt,
with step size γ > 0 that is taken close to zero, and a variable term h t ∈ [0, 1]. We introduce the family of quadratic Lyapunov functions:

(3.6) V a (1) t ,a (2)
t ,Pt (X t , Xt , t) = a (1) t a (2) t ⊤ f (X t ) -f ⋆ f ( Xt ) -f ⋆ + X t -X ⋆ Xt -X ⋆ ⊤ P t X t -X ⋆ Xt -X ⋆ ,
where a

(1)

t , a (2)
t ,Pt , and SDEs (3.5), the performance estimation problem is formulated by:

max f ∈F0,∞,d∈N,Xt∈R d d dt EV a (1) t ,a (2)
t ,Pt (X t , Xt , t),

subject to dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , d Xt = X t -Xt t dt.
A possible control of this quantity is presented in Theorem 3.4.

Theorem 3.4. Let f ∈ F 0,∞ be twice differentiable functions, and SDEs be taken under Polyak-Ruppert averaging (3.5) starting from x 0 ∈ R d , where d ∈ N is the dimension, and V be quadratic Lyapunov functions defined by (3.6). The following assertions are equivalent,

• The inequality d dt EV(X t , Xt , t) 1 2
Tr((a

(1) t ∇ xx f (X t ) + 2p (11) t I d )Σ(X t ))h 2 t γ is satisfied for all functions f ∈ F 0,∞ ,
t +λ (3) t 2 * * -a (1) t 0 * * * 0      0, ȧ(1) t + λ (1) t + λ (5) t = λ (4) t + λ (6) t , ȧ(2) t + λ (2) t + λ (6) t = λ (3) t + λ (5) t + ȧ(1) t .
Proof. The proof follows the methodology from Section 2.1.1, and using the Gram matrix G = P ⊤ P , where P = (X t -x ⋆ , Xt -x ⋆ , g t , ḡ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 (V 0,a [START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF] t ,Pt ). tht is a non-increasing function, the Lyapunov function

V(X t , Xt , t) = a (2) t (f ( Xt ) -f ⋆ ) + a (2) t 2th t X t -X ⋆ 2 verifies d dt EV(X t , Xt , t) a (2) t t h t Tr(Σ(X t )
) for all twice differentiable functions f ∈ F 0,∞ , all trajectories X t generated by SDEs (3.5) and all dimensions d ∈ N. Then,

it holds that E[f ( Xt ) -f ⋆ ] x 0 -x ⋆ 2 2a (2) t + γ 2a (2) t t 0 a (2) s s h s Tr(Σ(X s ))ds.
Proof. The proof follows from Theorem 3.4 (see Appendix B.1).

When a

(2) t = t maximal possible value), the step size verifies 0. The variance term does not diverge if and only if h t is constant. Then, a convergence bound is given by E

[f ( Xt ) -f ⋆ ] x0-x⋆ 2 2t
+ 1 2 Tr(σ)γh. The decreasing condition on t → a (2) t tht suggests a trade-off between converging and diminishing step size, as obtained without averaging.

Under the assumptions of Corollary 3.5, let us consider a bounded covariance matrix Σ(X t ) Σ, a step size h t = 1 (t+1) α , and a

(2) t = t β for α 0 and 0 β 1 some parameters. The decreasing condition imposes α + β 1. A different behavior is expected from α and β: on the one hand, an ideal step size should be large (α small), and on the other hand, we aim at converging as fast as possible (β large). The term that forgets the initial conditions is bounded by O( Under weighted averaging, and with

C u t = ut t 0 usds , the SDE is (3.7) dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , d Xu t = (X t -Xu t )C u t dt.
We study convergence of this generalized version of Polyak-Ruppert averaging, and compare it to traditional averaging techniques. Theorem 3.6. Let f ∈ F µ,∞ be twice differentiable functions, possibly strongly convex µ 0, and SDEs be given by (3.7), starting from

x 0 ∈ R d where d ∈ N is the dimension. Assuming ˙ ut 2ht
2µu t , the Lyapunov function

V(X t , Xu t , t) = u t C u t (f ( Xu t ) -f ⋆ ) + u t 2h t X t -x ⋆ 2 , verifies d dt EV(X t , Xu t , t) 1 2 u t h t γTr(Σ(X t )
) for all twice differentiable functions f ∈ F 0,∞ , all trajectories (X t , Xt ) generated by SDEs with averaging (3.7), and all dimensions d ∈ N.

Then, it holds that E[f ( Xu t ) -f ⋆ ] x 0 -x ⋆ 2 u 0 2h 0 t 0 u s ds + γ 4 t 0 u s ds t 0 u s h s Tr(Σ(X s ))ds.
Proof. These results are obtained by replacing 1 t → C u t (see Appendix B.2). Let us consider convex functions f ∈ F 0,∞ , a bounded covariance matrix Σ(X t ) Σ, a step size h t = 1 (t+1) α , and an averaging function u t = 1 (t+1) β , with α, β 0. From Theorem 3.6, the terms that forgets the initial conditions is in For strongly convex functions f ∈ F µ,∞ , polynomial convergence can be reached for the term that contains the initial conditions. However, the variance term cannot converge faster than the step size. Given a bounded covariance matrix Σ(X t ) Σ, the variance term behaves after integrating by part t 0 ushsds t 0 usds = h t -t 0 ( u) ḣsds t 0 usds . For diminishing step sizes, t 0 ushsds t 0 usds h t . Hence, averaging allows a better convergence for the terms containing initial conditions, but does not play a role in the variance term. To conclude, weighted averaging does not improve convergence results obtained under Polyak-Ruppert or primal averaging. The trade-off between the forgetting of the initial conditions and the noise term mostly relies on step sizes.

We have analyzed convergence of SGD together with averaging techniques in continuous-time, using approximating SDEs (3.3). Compared with discrete time, the continuous-time analysis leads to similar convergence results, while benefiting from simpler formulations, fewer assumptions especially on step sizes. Using this approach, we analyzed the trade-off between non-uniform averaging and step sizes, paving the way to a better understanding of averaging techniques. In the next session, we explore new convergence analyses for stochastic accelerated methods.

4. Accelerating the gradient flow. For both stochastic and deterministic models, we have retrieved known convergence results for continuous-time models approximating optimization methods. In this section, we provide convergence guarantees for second order gradient flows, and in particular for AGF (2.14).

In the deterministic setting, convergence of gradient descent was improved using a momentum. In this section, let f ∈ F 0,∞ be twice differentiable function and γ > 0 be constant step sizes. Li [START_REF] Li | Stochastic modified equations and dynamics of stochastic gradient algorithms i: Mathematical foundations[END_REF]Theorem 16,Section 4.4] proved that Nesterov accelerated gradient has the approximating SDE (for order-1 weak approximations), (4.1)

d 2 X t + 3 t dX t + ∇f (X t )dt + γΣ(X t )dB t = 0.
As for SGD, the Lyapunov function V(x, t) = t 2 (f (x)-f ⋆ )+2 (x-x ⋆ )+ t 2 ẋ 2 obtained from Theorem 2.6 does not allow to conclude about convergence to a stationary point of trajectories X t generated by stochastic accelerated gradient flows (4.1). Assuming bounded covariance Σ(X t ) Σ and applying Ito's formula to V

(•) along X t , E[f (X t ) -f ⋆ ] x 0 -x ⋆ 2 t 2 + γTr(Σ)3t.
In the following, we explore Polyak-Ruppert averaging together with diminishing step sizes, to analyze convergence of second-order SDEs.

4.1. Averaging does not preserve convergence rates. Averaging was a key to success for improving convergence of SGD (see Section 3). It is natural to wonder if averaging preserves the acceleration of Nesterov's gradient flow [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]. Let us define stochastic sedond-order gradient flows with Polyak-Ruppert averaging:

(4.2) d 2 X t + β t dX t + ∇f (X t )dt + γΣdB t = 0, d Xt = X t -Xt t dt,
where β t 0 is a function, and a family of quadratic Lyapunov functions,

(4.3) V a (1) t ,a (2)
t ,Pt (X t , Xt , t) = a (1) t a (2) t ⊤ f (X t ) -f ⋆ f ( Xt ) -f ⋆ +   Ẋt X t -X ⋆ Xt -X ⋆   ⊤ P t   Ẋt X t -X ⋆ Xt -X ⋆   ,
where P t 0 and a

(1) t , a (2)
t 0 are differentiable functions.

Theorem 4.1. Let f ∈ F 0,∞ be twice differentiable functions, V = V a (1) t ,a (2)
t ,Pt be quadratic Lyapunov functions, and stochastic accelerated gradient flows under Polyak-Ruppert averaging (4.2) with constant step sizes h t = 1, starting from x 0 ∈ R d , where d ∈ N is the dimension. Then it holds that:

• When a

t = 0, if the function V verifies d dt EV(X t , Xt , t) Tr(2p (1)
γΣ(X t )) for all twice-differentiable functions f ∈ F 0,∞ , all trajectories (X t , Xt ) generated by SDEs (4.2), and all dimensions d ∈ N, then V = 0.

• When a (2) t = 0, the Lyapunov function V(X t , t) = a (1) t (f (X t ) -f ⋆ ) + 1 2a (1) t a (1) t Ẋt + ȧ(1) t (X t -x ⋆ ) 2 , verifies d dt EV(X t , t) Tr(2a (1)
t γΣ(X t )), with a

(1) t t 2 for all functions f ∈ F 0,∞ , all trajectories X t generated by stochastic second-order gradient flows (4.2) and all dimensions d ∈ N. Then, it holds that

E[f (X t ) -f ⋆ ] 2 x 0 -x ⋆ 2 a (1) t + 1 2a (1) t γ t 0 a (1)
s Tr(Σ(X s ))ds.

Proof. The proof follow from Theorem 2.10 (see Appendix B.3).

To conclude, with constant step sizes, there is no Lyapunov functions that allows to forget the initial conditions while reducing the variance term. Primal averaging leads to similar results (see Appendix C). Therefore, averaging plays a different role in second-order SDEs than in first-order SDEs under Polyak-Ruppert averaging.

4.2. Diminishing step sizes. Averaging was not conclusive for finding a convergence guarantee of Nesterov's accelerated gradient flow with a diffusion term. Let us consider second-order stochastic gradient flows with varying step sizes h t 0, (4.4)

d 2 X t + β t dX t + h t ∇f (X t )dt + h t γΣ(X t )dB t = 0.
We propose a Lyapunov function among the class of quadratic functions (2.11).

Theorem 4.2. Let f ∈ F 0,∞ be twice differentiable functions, X t be a trajectory generated by stochastic second-order gradient flows (4.4) with varying step sizes h t 0 and γ > 0, starting from x 0 ∈ R d , where d ∈ N is the dimension. Assuming ȧt a t 2 3 (β t + 1 2 ḣt ht ) and t → ( ȧt) 2 2htat a decreasing function, the Lyapunov function

V(X t ) = a t (f (X t ) -f ⋆ ) + 1 2h t a t a t Ẋt + ȧt (X t -x ⋆ ) 2 ,
verifies d dt EV(X t , t) 1 4 Tr(a t h t Σ(X t ))γ for all functions f ∈ F 0,∞ , all trajectories X t generated by stochastic second-order gradient flows (4.4) and all dimensions d ∈ N.

Proof. This result is obtained extending the LMI for ODEs from Theorem 2.10 and Corollary 2.11 to varying step sizes.

Let us now assume that h t = and α + β 2, and

E[f (X t ) -f ⋆ ] β 2 t β x 0 -x ⋆ 2 + γ 4t β t 0 s β (s + 1) α Tr(Σ(X s ))ds.
We assume bounded covariance of Σ(X t ) Σ, then β min( 2b-α 3 , 2 -α). On the one hand, the smaller the step sizes, the better the convergence for the term that contains the initial conditions. On the other hand, the variance term behaves as 1

t β if β
α -1, and as 1 t α-1 otherwise (convergence requiring then α 1 and β 1). For Nesterov's accelerated gradient flow with b = 3, we have β = 2 -α α -1, and therefore α 3 2 . Taking α = 3 2 , a convergence bound is given by:

E[f (X t ) -f ⋆ ] 9 4 √ t x 0 -x ⋆ 2 + log t √ t γTr(Σ).
We retrieve result from Corollary 3.5 for the SDE approximating SGD with Polyak-Ruppert averaging, with smaller step sizes. It does not seem possible to accelerate SGD when diminishing the step size. Ghadimi and Lan [31, Corollary 3] proved a convergence bound for a randomized stochastic accelerated gradient method with β = 2 and α = 1 2 , that we do not retrieved. However, in their approach, the function is minimized over a compact, convex domain, whereas our approach focuses on an unbounded domain.

Conclusion and future work.

We have developed a systematic approach for generating Lyapunov functions for families of ODEs and SDEs. Verifying such a Lyapunov function can be formulated as an LMI. From this formulation, it is possible to derive Lyapunov functions and the associated convergence bounds. We retrieve results from discrete optimization methods with shorter proofs, and fewer assumptions.

While obtaining guarantees for stochastic optimization methods might be tedious, the SDE approach allows for simpler analyses of the trade-off between the variance term and term that forgets the initial conditions. A shortcoming of this approach is that this analysis does not include approximation guarantees between optimization methods and their continuous-time counterparts. In the deterministic setting, stability techniques are often developed to quantify this approximation efficiency [START_REF] Gautschi | Numerical Analysis[END_REF]. In stochastic analysis, stochastic modified equation have been introduced by Li et al. [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF]Theorem 1] to better approximate SGD, and stochastic methods with momentum. For non-convex functions, some analysis have also been done by Shi et al [START_REF] Shi | On learning rates and schr\" odinger operators[END_REF]. However, these approximation theorems often require many assumptions on the class of functions, which we believe could be further simplified using computer-assisted proofs.

program,

0 max G 0,F ∈R 2 b ⊤ 0 F + Tr(A 0 G), subject to b ⊤ 1 F + Tr(A 1 G) 0, b ⊤ 2 F + Tr(A 2 G) 0,
where

A 0 = ċt -c t -c t -a t , A 1 = -µ/2 1/2 1/2 0 , A 2 = -µ/2 0 0 0 , b 0 = ȧt [1, -1] ⊤ b 1 = [-1, 1] ⊤ and b 2 = [1, -1] ⊤ .
In this convex case, strong duality holds via Slater's conditions [START_REF] Taylor | Smooth strongly convex interpolation and exact worst-case performance of first-order methods[END_REF]Theorem 6]. Consider the Lagrangian dual of the SDP, min

λ (1) t ,λ (2) t 0 0 subject to A 0 + λ (1) t A 1 + λ (2) t A 2 0 λ (2) t b 2 + λ (1) t b 1 + b 0 = 0.
This problem is a feasibility problem, min

λ (1) t ,λ (2) t 0 0 subject to S = ċt -µ 2 (λ (1)
t + λ (2) t ) -c t + λ (1) t 2 -c t + λ (1) t 2 -a t 0, ȧt = λ (1) t -λ (2)
t .

When µ = 0, this formulation is exactly the LMI formulation obtained in Theorem 2.3. When µ > 0, the LMI corresponds to Theorem 2.1 for ȧt = e 2µt a and ċt = e 2µt c.

A.

+ at 2 * * 0    0, ȧt = λ (1) t -λ (2)
t .

Convex functions. Let us first assume the functions f to minimize are convex (µ = 0). While determining a suitable Lyapunov function V at,Pt , we try to find a large parameter a t at time t. To this end, we saturate ȧt = λ

(1) t (λ

(2) t = 0). The LMI can then be simplified into: Assuming for instance a 0 = 0 and p (11) 0

S = ṗ (
= 2, we obtain:

a t min(t 2 , lim ǫ→0 a ǫ e ⊤ ǫ 2 3 βsds ).
Strongly convex functions. Similar result can be obtain assuming an exponential form for the parameters a t = ae tτ , and P t = P e tτ , where τ > 0 is a convergence parameter to determine.

Taking a generic form for β t = r t ), the condition lim ǫ→0 a ǫ e ⊤ ǫ 2 3 βsds takes more sense. Indeed, a t = t 2r/3 is a solution to a t = lim ǫ→0 a ǫ e Proof for Theorem 3.4. We consider the stochastic gradient flow under Polyak-Ruppert averaging (3.5), that is recalled here:

dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , d Xt = X t -Xt t dt,
where h t is the varying part of the step size, and γ > 0 the maximum allowed step size.

Let us consider the family of quadratic Lyapunov functions defined by

(B.1) V a (1) t ,a (2)
t ,Pt (X t , Xt , t) = a (1) t a (2) t ⊤ f (X t ) -f ⋆ f ( Xt ) -f ⋆ + X t -X ⋆ Xt -X ⋆ ⊤ P t X t -X ⋆ Xt -X ⋆ ,
where a

t , a

t ,Pt (X t , Xt , t), subject to

dX t = -h t ∇f (X t )dt + h t (γΣ(X t )) 1/2 dB t , d Xt = X t -Xt t dt.
Let us rewrite the SDE,

dX t d Xt = 0 0 1 t -1 t X t Xt dt + -h t 0 0 0 ∇f (X t ) ∇f ( Xt ) dt + h t (γΣ(X t )) 1/2 0 dB t .
Let us denote Y t = X t Xt , and V a (1)

t ,a (2)
t ,Pt (X t , Xt , t) = V(Y t , t), then d dt V(Y t , t) = ∂ ∂t V(Y t , t) + ∂ ∂y V dYt dt + 1 2 γh 2 t Tr[ ∂ 2 ∂Y 2 t V(Y t , t) ⊤ Σ(X t ) 0 ]
. Then, we have,

d dt EV a (1) t ,a (2)
t ,Pt (X t , Xt , t) = d dt EV(Y t , t), = E[ ∂ ∂t V(Y t , t) + ∂ ∂y V dY t dt ] + 1 2 γh 2 t Tr( ∂ 2 ∂Y 2 t V(Y t , t) ⊤ Σ(X t ) 0 ).
The noise term only depends on the second derivative of f in X t and contains a noise term Σ(X t ). We propose to not take the noise term into account in the SDP formulation. The inequality

(B.2) d dt EV a (1) t ,a (2)
t ,Pt (X t , Xt , t)

0 max Xt∈R d ,d∈N, f ∈F0,∞ d dt V a (1) t ,a (2)
t ,Pt (X t , Xt , t)

subject to dX t = -h t ∇f ( Xt )dt, d Xt = X t -Xt t dt.
Using convex interpolation on triplets {(X t , g t , f t ), ( Xt , ḡt , ft ), (X ⋆ , 0, f (X ⋆ ))} such that there exists f ∈ F 0,∞ verifying f t = f (X t ), g t = ∇f (X t ), f ( Xt ) = ft and ft = ∇f ( Xt ), the infinite dimensional problem has a quadratic program formulation.

Introducing the Gram matrix

G = Y ⊤ Y 0, Y = [X t -x ⋆ , Xt -x ⋆ ,
t + λ

(4) t = λ (5) t + λ (1) t + ȧ(1) t , λ (2) t + λ (6) t + ȧ(2) t = λ (3) t + λ (5)
t , λ

t + λ (2) t + ȧ(2) t = λ (4) t + λ (5) t + ȧ(1) t , λ (i) t 0, i ∈ {1, ..., 6}, (1)
where λ (i) t , i ∈ {1, ..., 6} are dual values associated with interpolation inequalities. These three inequalities can be reduced to two independent equalities.

Proof for Corollary 3.5. Using the LMI formulation of Theorem 3.4, the Lya-

punov function V(X t , Xt , t) = a (2) t (f ( Xt ) -f ⋆ ) + p (11) t X t -X ⋆ 2 satisfies the LMI:      ṗ(11) t 0 λ (6) t +λ (4) t 2 -h t p (11) t a
(2) t 2t - λ (5) t 2 * 0 - λ (6) t 2 λ (3)
t +λ

t 2t * * 0 0 * * * 0      0, λ (5) t 2 - a (2)
t + λ

(4) t = λ (5) t + λ (1) t + ȧ(1) t , λ (2) t + λ (6) t + ȧ(2) t = λ (3) t + λ (5)
t , λ (i) t 0, i ∈ {1, ..., 6}.

From the LMI, we conclude that λ Proof. Following the same reasoning as in the proof for Theorem 2.3, the worstcase guarantee can be formulated as the maximization problem:

(3) t = λ (6) t = 0, λ (5) t = a (2) t t , λ (4)
0 max Xt∈R d ,d∈N, f ∈F0,∞ d dt V a (1) t ,a (2)
t ,Pt (X t , Xt , t) subject to dX t = -h t ∇f (X t )dt + h t (γΣ(X t ))

- λ (6) t 2 λ (3) t +λ (5) t 2 -a (2) t C u t * * -a (11) t 0 * * * 0         0, λ (6)
t + λ

(4) t = λ (5) t + λ (1) t + ȧ(1) t , λ (2) t + λ (6) t + ȧ(2) t = λ (3) t + λ (5)
t , λ Proof. We follow the same reasoning as in Theorem 2.3 and Theorem 3.4, the worst-case guarantee can be formulated as the maximization problem:

t + λ (2) t + ȧ(2) t = λ (4) t + λ (5) t + ȧ(1) t , λ (i) t 0, i ∈ {1, ..., 6}. Considering the Lyapunov function V(X t , Xu t , t) = a (2) t (f ( Xu t )-f ⋆ )+p (11) t X t -x ⋆ 2 , it follows that ȧ(2) t C u t a (1)
0 max f ∈F0,∞,d∈N,Xt∈R d d dt V a (1) t ,a (2)
t ,Pt (X t , Xt , t), subject to d 2 X t + β t dX t + ∇f (X t )dt + γΣdB t = 0, d Xt = X t -Xt t dt, where V a (1) t ,a (2)
t ,Pt (X t , Xt , t) = a (1) t a (2) t ⊤ f (X t ) -f ⋆ f ( Xt ) -f ⋆ +   Ẋt X t -X ⋆ Xt -X ⋆   ⊤ P t   Ẋt X t -X ⋆ Xt -X ⋆   with P t
0, with differentiable coefficients, and a

t , a

(2) t are nonnegative differentiable functions. Introducing the Gram matrix G = P ⊤ P with P = [ Ẋt , X t -X ⋆ , Xt -x ⋆ , ∇f (X t ), ∇f ( Xt )], the LMI equivalent reformulation is given by:

             ṗ(
t -a

(2)

t /t 2 * * * 0 0 * * * * 0              0, λ (6)
t + λ

(4) t = λ (5) t + λ (1) t + ȧ(1) t , λ (2) t + λ (6) t + ȧ(2) t = λ (3) t + λ (5)
t , λ

t ,Pt = 0. Otherwise, the LMI can be simplified to the LMI of Theorem 2.10.

Appendix C. Primal averaging. In this appendix, we analyze primal averaging, that is to compare with Polyak-Ruppert averaging. Both continuous-time models happens to have very similar behavior.

C.1. First-order gradient flow. Let us consider primal averaging, introduced by Tao et al. [START_REF] Tao | Primal averaging: A new gradient evaluation step to attain the optimal individual convergence[END_REF], that led to a simpler analysis of averaging in SGD [START_REF] Taylor | Stochastic first-order methods: non-asymptotic and computer-aided analyses via potential functions[END_REF]Theorem 7] (C.1) dX t = -h t ∇f ( Xt )dt + h t (γΣ(X t )) Primal averaging involves at most two interpolation inequalities whereas Polyak-Ruppert averaging may involve up to four inequalities to obtain the same convergence guarantee. Compared with Polyak-Ruppert averaging, primal averaging has similar convergence guarantees but requires less interpolation inequalities.

Remark C.2. In discrete time, convergence of primal averaging is often proven using its connection to the heavy ball method [START_REF] Flammarion | From averaging to acceleration, there is only a step-size[END_REF][START_REF] Ghadimi | Global convergence of the heavy-ball method for convex optimization[END_REF]. In continuous time, primal averaging also allows to describe the dynamics of the averaged sequence by a secondorder SDE Ẍt + 2 t Ẋt + 1 t ∇f ( Xt )+ √ γΣ(Xt) t dBt dt = 0. The convergence bound above can be obtained directly from this SDE. This analysis allows to understand the relationship between averaging and accelerated gradient flows.

Under primal averaging, the convergence in function values is exactly the same as under Polyak-Ruppert averaging, but involves only one interpolation inequality. When considering optimization methods, [START_REF] Taylor | Stochastic first-order methods: non-asymptotic and computer-aided analyses via potential functions[END_REF]Theorem 7] also pointed out that primal averaging leads to shorter proofs. C.2. Second order-gradient flow. We have seen primal averaging for stochastic first-order gradient flows (3.3) leads to a second-order SDE. For the accelerated gradient flow under primal averaging, it is possible to obtain a third order SDE, A natural quadratic Lyapunov is given by V at,Pt (X t ) = a t (f (X t ) -f ⋆ ) + U ⊤ t P t U t , where U t = X t -X ⋆ , Ẋt , Ẍt ⊤ , and a t is a nonnegative differentiable function, and P t 0 has differentiable parameters.

d
Theorem C.3. Let f ∈ F 0,∞ be twice differentiable functions, β t 0 be differentiable functions, and consider the following ODE starting from x 0 ∈ R d where d ∈ N is the dimension,

d 3 X t + α t d 2 X t + β t dX t + γ t ∇f (X t )dt = 0.
If a quadratic Lyapunov function V at,Pt as defined above verifies d dt V at,Pt (X t , t) 0 for all twice-differentiable functions f ∈ F 0,∞ , all solutions X t to ODEs defined above, and all dimensions d ∈ N, then V at,Pt = 0.

Proof. Following the reasoning from Appendix A. and a t = 0, and the Lyapunov V at,Pt = 0 is the unique solution to this LMI.

Theorem C.3 concludes in an other way that it is not possible to obtain a convergence guarantee on an averaged sequence of an accelerated method. However, we can wonder if there exists other transformations that ensures convergence of both the averaged sequence and the noise.

  Strong convexity parameter µWorst-case convergence rate τ 2µ PEP (a) Worst-case rate τ on the class of quadratic Lyapunov functions for the gradient flow.

1 Fig. 1 :

 11 Fig. 1: Comparison between numerical values for τ obtained by solving the LMI (2.7) and the reference established in the literature [33, Proposition 1.1], for trajectories X t generated by gradient flow (1.3) on µ-strongly convex functions.

  all trajectories X t generated by Polyak damped oscillator, and all dimensions d ∈ N.Proof. Taking λ 1 = 4/3 √ µ, λ 2 = 0, ν 1 = 0 and ν 2 = 1, we verify the LMI for this Lyapunov function V, with τ = 4/3 √ µ.

1 .

 1 Relaxing the condition P 0 improves results from Sanz-Serna and Zygalakis and Wilson et al. by a factor 4/3 in Corollary 2.5.

  Worst-case guarantee optimized over the class of quadratic Lyapunov functions.

  Lyapunov parameters P for τ = 4/3 √ µ and a = 1, as a function of the condition number µ.

√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:

Corollary 3 . 5 (

 35 Averaging and diminishing step sizes). Let f ∈ F 0,∞ be twice differentiable functions, SDEs be taken under Polyak-Ruppert averaging (3.5) starting from x 0 ∈ R d where d ∈ N is the dimension. Assuming ȧ

1 ( 2 (

 12 t+1) 1-β , and the variance term in 1 (t+1) α . Both terms converge at same rate for α = β = 1 up to log (t + 1)). We retrieve convergence results for SGD under Polyak-Ruppert averaging [2,Theorem 6].

⊤ ǫ 2 3

 2 βsds (which requires a 0 = 0). Appendix B. Proofs for SDEs. B.1. Proof for Theorem 3.4 and Corollary 3.5.

  Scieur et al. proved in [33, Proposition 1.1] a convergence bound in function values for gradient flows originating from

  all trajectories X t generated by Polyak damped oscillator, and all dimensions d ∈ N. This Lyapunov is a feasible point of the LMI (2.13) from Theorem 2.4, with τ

  all trajectories X t generated by the accelerated gradient flow and all dimensions d ∈ N.• There exist λ

				(1) t , λ (2) t	0 such that
	S =	  	ṗ(11) t *	p 2(p (11) t (12) t	-3 t p -3 t p (12) t (22) t	+ ) + ṗ(12) t ṗ(22) t	-p -p (12) t (22) + t	(1) λ t 2

  is bounded by two, enforcing short proofs. Theorem 2.6 allows to retrieve the Lyapunov function exhibited by Su et al.[START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF] Theorem 3], and its associated convergence guarantee. Choosing a t = t 2 and P t =

	2	1 t/2 t 2 /4 t/2	for the Lyapunov parameters, and λ (1) t	= t, λ

  1 t 2 convergence (see Theorem 2.6) in function values. Considering the family of quadratic Lyapunov functions (2.11), we study convergence of non-autonomous second-order gradient flows and draw comparison with the accelerated gradient flow. Let β t Remark 2.9. Wibisono et al. [46, Theorem 2.2] proved this ODE is related to the family of ODEs defined by Ẏt + βt Y t + αt ∇(Y t ) = 0. For α t a strictly nonnegative differentiable function, a time change formula τ t =

		0 be a continuous function, and a second-order
	non-autonomous ODE,
	(2.16) (2.16), the trajectory Y t = X τt is solution to Ÿt + (β τt Ẍt + β t 0 √ √ α s ds, and X t a solution to α t -αt 2αt ) Ẏt + α t ∇f (Y t ) = 0.
	For any two functions β	(1) t , β

t Ẋt + ∇f (X t ) = 0.

  all trajectories (X t , Xt ) generated by SDEs under Polyak-Ruppert averaging (3.5), and all dimensions d ∈ N.• There exist λ

			(1) t , ..., λ (6) t	0 such that,
	    	ṗ11 t + * 2p (12) t t	ṗ12 t -ṗ22 t -p (12) t t + 2p (22) p (22) t t t t	λ (6) t +λ (4) t 2 -λ (6) t 2 -h t p -h t p (12) (11) t t	-	a 2t -(2) t (2) t 2t + a λ (5) (5) λ t 2

  Weighted averaging. Polyak-Ruppert performs uniform averaging of trajectories X t over the time step. We introduce weighted averaging to analyze SGD, that is defined with respect to a function u t 0[START_REF] Roux | Anytime tail averaging[END_REF],

	1 t β ), and the variance term t t) by O( 1 2t β . Hence, a natural choice is α = β = 1 2 , retrieving results from [2, Theorem 4] [41, Table 2] in discrete optimization.
	3.3.2. xu t =	1 0 u s ds t

α ) if β = α. When α = β, the variance term is in log(t 0 u s x s ds.

  2. Proof for Corollary 2.11.

	Proof. From Theorem 2.3, we get a LMI reformulation to finding a suitable Lya-
	punov function: There exist λ (1) t , λ (2) t	0 such that,
	   -µ 2 (λ (1) t + λ (2) t ) + ṗ 11 *	p 2(p (11) t (12) t	-β t p -β t p (12) t (22) t	+ ) + ṗ(12) t ṗ(22) t	-p -p (12) t (22) t	+ λ1 2

  [START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF] From linear algebra, the matrix S is semi-definite negative if and only if the diagonal terms are negative. We obtain two conditions on a t : t a t . For all ǫ > 0, after integrating between ǫ, and t, that a t a ǫ e βsds . Therefore, a t lim ǫ→0 a ǫ e

						ṗ(11) t	0,
						ȧt -β t a t +	1 2	ȧt 0.
	In addition, we assumed P t 0 to enforce positivity of the Lyapunov function V at,Pt , which is equivalent to p (11) t (p (12) t ) 2 p (22) t = ( ȧt) 2 2at .
	On the one hand, we conclude that the function ȧt 2 √ at	p	(11) 0	, and after integrating
	between 0 and t, that conclude that ȧt 2 3 β ⊤ √ a t ǫ 2 3 ⊤ √ a 0 + (11) p 0 2 t (a t = O(t 2 )). On the other hand, we ǫ 2 3 βsds , and a bound on a t is given by:
	a t min(( √ a 0 + ( p	(11) 0	/2)t) 2 , lim ǫ→0	a ǫ e	⊤ ǫ	2 3 βsds ).
	p	(22) t	=	t * a t 2	,	p 2(p (11) t (12) t	-β t p -β t p (12) t (22) t	+ ) + ṗ(12) t ṗ(22) t	0,
	p	(12) t	=	ȧt 2	.

  t are differentiable real functions, and P t = p

	that can be cast as a minimization problem at time t,
	max Xt∈R d ,d∈N, f ∈F0,∞	d dt	EV a (1) t ,a
			p	(11) t (12) t	p p	(12) t t (22)	0 has dif-
	ferentiable parameters. We look for a worst-case guarantee in the Lyapunov approach,

  g t , ḡt ], the maximization problem is reformulated as an SDP program, whose dual is the following LMI,

	    	ṗ(11) t	+ * * *	2p (12) t t	ṗ(12) t ṗ(22) -t	p (12) t t + -2p (22) p (22) t t t t * *	λ (6) t +λ (4) t 2 -h t p (12) -h t p t -λ (6) (11) t t 2 * -a t (11)	a 2t -(2) t (3) λ t +λ (5) t 2 0 0	λ (5) t 2 -a 2t (2) t	    	0,
	λ

  1/2 dB t ,

	where xu t = given by:	1 0 usds t	⊤ 0 u s x s ds and C u t = ut t 0 usds . Its LMI equivalent reformulation is
	      	ṗ(11) t +2p (12) t *	C u t	ṗ(12) t (12) t -p -2p (22) t ṗ(22) +(p t (22) t C u )C u t t	λ (6) t +λ (4) t 2 (12) -h t p -h t p t	(11) t	a	(2) t C u t -	(5) λ t 2
	

d Xu t = (X t -Xu t )C u t dt,

  [START_REF] Fazlyab | Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems[END_REF]

	t -2β t p * *	(11) t	ṗ(12) + p (13) t t + p ṗ(22) t + 2 p (22) (12) t (22) t t t *	ṗ(13) t p (12) t β t + p -p (13) t t -(23) t p (23) t -2 2p (22) t t + p (33) t t ṗ(33) t -2p (33) t t	-p (1) t -p (11) t (12) t λ (4) t +λ (6) t + a 2 -λ (6) t 2 -p (13) t	0 2t -a (2) t (5) λ t +λ	(5) λ t 2

t -β t p

  = 0, and because P t is positive semidefinite, p = 0, the unique feasible Lyapunov function is V a[START_REF] Attouch | Rate of convergence of the nesterov accelerated gradient method in the subcritical case α 3[END_REF]

	t + λ (2) t +	ȧ(2) t = λ (4) t + λ (5) t +	ȧ(1) t ,
	λ (i) t	0, i ∈ {1, ..., 6}.
	Thus, p a (1) t = 0. If in addition a (11) t	(2)	(12) t	= p	(13) t	= 0. Then,

t t ,a

  1/2 dB t , Unlike Polyak-Ruppert averaging, the gradient is evaluated on the averaged sequence Xt , that drives the dynamics on X t . The family of Lyapunov functions V a (1) As for Polyak-Ruppert averaging, we consider the family of Lyapunov functions V 0,a (2) t ,Pt . It is possible to obtain an LMI condition on the Lyapunov, as done in Theorem 3.4.Theorem C.1. Let f ∈ F 0,∞ be twice differentiable functions, and SDEs (C.1) starting from x 0 ∈ R d , where d ∈ N is the dimension and originating from convex ∈ F 0,∞ , for all trajectories (X t , Xt ) generated by the SDEs (C.1), and all dimensions d ∈ N.Then, it holds that E[f ( Xt ) -f ⋆ ] x 0 -x ⋆Proof. This results is obtained following the same reasoning as for Polyak-Ruppert averaging in Theorem 2.3. The equivalent LMI to inequalityd dt EV(X t , t) Tr(Σ)γh t ,for any twice-differentiable function f ∈ F 0,∞ and any trajectory (X t , Xt ) generated by the SDE (C.1), is given by ∈ {1, ..., 6}.The Lyapunov function defined by V(X t , Xt , t) = a

								2a	(2) t h 0	2	+	γ 2a (2) t	0	t	Tr(Σ(X s ))h s	a	(2) s s	ds.
	     λ (6) ṗ(11) t t + λ + * * * (4) 2p (12) t t t = λ (5) t + λ ṗ(12) t ṗ(22) -t (1) t + p (12) t t + -2p (22) p (22) t t t t * * ȧ(1) t ,			λ (6) t +λ (4) t 2 -λ (6) t 2 (11) t -a *			a 2t -(2) t (3) λ t +λ (5) t 2	λ (5) t 2 -h t p -a (2) t 2t -h t p (11) t (12) t 0 0	    	0,
	λ (2) t + λ (6) t +	ȧ(2) t = λ (3) t + λ (5) t ,
	λ (1) t + λ (2) t +	ȧ(2) t = λ (4) t + λ (5) t +	ȧ(1) t ,
	λ (i) t verifies the following LMI 0, i (2)	(11) t	X t -x ⋆	2
	d Xt = ṗ(12) t -p (12) X t -Xt t t t + p (22) dt. t t ṗ(22) t -2p (22) t t * * t + λ = λ 2p (12) t t (4) + * * * t + λ      ṗ(11) t λ (6) t (5) (1) t ,	λ (6) t +λ (4) t 2 -λ (6) t 2 0 *			a 2t -(2) t (3) λ t +λ (5) t 2	λ (5) t 2 -h t p -a (2) t 2t -h t p (11) t (12) t 0 0	    	0,
	λ (2) t + λ (6) t +	ȧ(2) t = λ (3) t + λ (5) t ,
	t ,a tht is a non-increasing function, the (2) t ,Pt I d )Σ)h 2 ȧ(2) (11) t = 0, and t ∇ xx f ( Xt )+2p (1) 0, i ∈ {1, ..., 6}. causes a noise 1 2 Tr((a λ (i) t Thus, λ (6) t = λ (4) t = λ (5) t = λ (1) t ȧ(2) t λ (2) t + t = λ (3) t = a (2) t 2t , t γ. functions f . Assuming ȧ(2) t a (2) t t and t → a (2) t following Lyapunov function ṗ(11) t = ˙ (2) a t 2th t 0.
		V(X t , t) = a	(2) t (f ( Xt ) -f ⋆ ) +	(2) t 2th t a	X t -X ⋆	2 ,
	verifies d dt EV(X t , t)	γ 2 Tr(Σ(X t ))h t	a	(2) t

t for all twice-differentiable functions f

t (f ( Xt) -f ⋆ ) + p

  1, the problem of finding a suitable Lyapunov function can be formulated as a maximization problem subject to d 3 X t + α t d 2 X t + β t dX t + γ t ∇f (X t )dt = 0,

	0	max f ∈F0,∞,d∈N,Xt∈R d	d dt	V at,Pt ,
	has the following LMI equivalence,
	          ȧt = λ ṗ(11) +2p * (12) t * * (1) t -λ (2) (11) t t .	ṗ(12) ṗ(22)	+p	(11) t (22) t + p -α t p (13) t (12) t + 2p (12) t (23) t * *	ṗ(13) +p ṗ(23) +p ṗ(33) (23) t (33) t t *	(13) t (13) t	(11) t (12) 2 -h t p -h t p at t λ (1) t (13) t 2 -h t p 0	         	0,
	It follows directly from the LMI that p Thus λ (1) t (2) = 0 = λ t	(11) t	= 0. Since P t 0, p	(11) t	= p	(12) t	= p	(13) t	= 0.

t -2α t p t -β t p t -α t p t -2β t p t -β t p

, whose derivative is decreasing along trajectories X t (proof: d dt V(X t , t) =

Acknowledgments. This work was funded by MTE and the Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grant SEQUOIA 724063).

Appendix A. Proof for ODEs.

A.1. Proofs for Theorems 2.3 and 2.1.

Proof. Let us consider X t a solution to the gradient flow (1.3) starting from x 0 ∈ R d , and a quadratic Lyapunov function of the form V at,ct

2 such that a t , c t are differentiable and V at,ct is nonnegative over the trajectory X t . The family of function f is supposed to be µ-strongly convex, with µ 0. We are going to prove the LMI equivalence of Theorem 2.3 and Theorem 2.1.

Let us first rephrase our objective, that is computing a Lyapunov function that is decreasing along X t . 0 max

This problem remains infinite-dimensional, but we can reformulate it thanks to an interpolation theorem of the class of functions F µ,L on (X i , g i , f i ), with 0 µ L ∞. A set {(X i , g i , f i ) i∈I } is said to be F µ,L -interpolable if there exist f ∈ F µ,L such that f i = f (X i ) and g i = ∇f (X i ) for all i ∈ I. It turns out that the formulation for this class of functions is quite simple.

interpolable if and only if the following set of conditions holds for every pair of indices i ∈ I and j ∈ J

.

This theorem allows to replace the problem above by a quadratic program, linear in f t , f ⋆ ) and quadratic in (x t , g t , x ⋆ ). 0 max

Let G be a Gram matrix defined by G = X t -x ⋆ 2 X t -x ⋆ , g t X t -x ⋆ , g t g t 2 0 and the vector F = [f t , f ⋆ ]. The quadratic program can thus be formulated into a semidefinite