# continuous_time_analysis_of_momentum_methods__538416cd.pdf

Journal of Machine Learning Research 22 (2021) 1-40 Submitted 6/19; Revised 9/20; Published 1/21

Continuous Time Analysis of Momentum Methods

Nikola B. Kovachki nkovachki@caltech.edu Computing and Mathematical Sciences California Institute of Technology Pasadena, CA 91125, USA

Andrew M. Stuart astuart@caltech.edu Computing and Mathematical Sciences California Institute of Technology Pasadena, CA 91125, USA

Editor: Suvrit Sra

Gradient descent-based optimization methods underpin the parameter training of neural networks, and hence comprise a signiﬁcant component in the impressive test results found in a number of applications. Introducing stochasticity is key to their success in practical problems, and there is some understanding of the role of stochastic gradient descent in this context. Momentum modiﬁcations of gradient descent such as Polyak s Heavy Ball method (HB) and Nesterov s method of accelerated gradients (NAG), are also widely adopted. In this work our focus is on understanding the role of momentum in the training of neural networks, concentrating on the common situation in which the momentum contribution is ﬁxed at each step of the algorithm. To expose the ideas simply we work in the deterministic setting.

Our approach is to derive continuous time approximations of the discrete algorithms; these continuous time approximations provide insights into the mechanisms at play within the discrete algorithms. We prove three such approximations. Firstly we show that standard implementations of ﬁxed momentum methods approximate a time-rescaled gradient descent ﬂow, asymptotically as the learning rate shrinks to zero; this result does not distinguish momentum methods from pure gradient descent, in the limit of vanishing learning rate. We then proceed to prove two results aimed at understanding the observed practical advantages of ﬁxed momentum methods over gradient descent, when implemented in the non-asymptotic regime with ﬁxed small, but non-zero, learning rate. We achieve this by proving approximations to continuous time limits in which the small but ﬁxed learning rate appears as a parameter; this is known as the method of modiﬁed equations in the numerical analysis literature, recently rediscovered as the high resolution ODE approximation in the machine learning context. In our second result we show that the momentum method is approximated by a continuous time gradient ﬂow, with an additional momentum-dependent second order time-derivative correction, proportional to the learning rate; this may be used to explain the stabilizing eﬀect of momentum algorithms in their transient phase. Furthermore in a third result we show that the momentum methods admit an exponentially attractive invariant manifold on which the dynamics reduces, approximately, to a gradient ﬂow with respect to a modiﬁed loss function, equal to the original loss function plus a small perturbation proportional to the learning rate; this small correction provides convexiﬁcation of the loss function and encodes additional robustness present in momentum methods, beyond the transient phase.

c 2021 Nikola B. Kovachki and Andrew M. Stuart.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided at http://jmlr.org/papers/v22/19-466.html.

Kovachki and Stuart

Keywords: Optimization, Machine Learning, Deep Learning, Gradient Flows, Momentum Methods, Modiﬁed Equation, Invariant Manifold

1. Introduction

1.1 Background and Literature Review

At the core of many machine learning tasks is solution of the optimization problem

arg min u Rd Φ(u) (1)

where Φ : Rd R is an objective (or loss) function that is, in general, non-convex and diﬀerentiable. Finding global minima of such objective functions is an important and challenging task with a long history, one in which the use of stochasticity has played a prominent role for many decades, with papers in the early development of machine learning Geman and Geman (1987); Styblinski and Tang (1990), together with concomitant theoretical analyses for both discrete Bertsimas et al. (1993) and continuous problems Kushner (1987); Kushner and Clark (2012). Recent successes in the training of deep neural networks have built on this older work, leveraging the enormous computer power now available, together with empirical experience about good design choices for the architecture of the networks; reviews may be found in Goodfellow et al. (2016); Le Cun et al. (2015). Gradient descent plays a prominent conceptual role in many algorithms, following from the observation that the equation du

dt = Φ(u) (2)

will decrease Φ along trajectories. The most widely adopted methods use stochastic gradient decent (SGD), a concept introduced in Robbins and Monro (1951); the basic idea is to use gradient decent steps based on a noisy approximation to the gradient of Φ. Building on deep work in the convex optimization literature, momentum-based modiﬁcations to stochastic gradient decent have also become widely used in optimization. Most notable amongst these momentum-based methods are the Heavy Ball Method (HB), due to Polyak (1964), and Nesterov s method of accelerated gradients (NAG) Nesterov (1983). To the best of our knowledge, the ﬁrst application of HB to neural network training appears in Rumelhart et al. (1986). More recent work, such as Sutskever et al. (2013), has even argued for the indispensability of such momentum based methods for the ﬁeld of deep learning. From these two basic variants on gradient decent, there have come a plethora of adaptive methods, incorporating momentum-like ideas, such as Adam Kingma and Ba (2014), Adagrad Duchi et al. (2011), and RMSProp Tieleman and Hinton (2012). There is no consensus on which method performs best and results vary based on application. The recent work of Wilson et al. (2017) argues that the rudimentary, non-adaptive schemes SGD, HB, and NAG result in solutions with the greatest generalization performance for supervised learning applications with deep neural network models. There is a natural physical analogy for momentum methods, namely that they relate to a damped second order Hamiltonian dynamic with potential Φ:

dt2 + γ(t)du

dt + Φ(u) = 0. (3)

Continuous Time Analysis of Momentum Methods

This perspective goes back to Polyak s original work Polyak (1964, 1987) and was further expanded on in Qian (1999), although no proof was given. For NAG, the work of Su et al. (2014) proves that the method approximates a damped Hamiltonian system of precisely this form, with a time-dependent damping coeﬃcient. The analysis in Su et al. (2014) holds when the momentum factor is chosen according to the rule

λ = λn = n n + 3, (4)

where n is the iteration count; this choice was proposed in the original work of Nesterov (1983) and results in a choice of λ which is asymptotic to 1. In the setting where Φ is µ-strongly convex, it is proposed in Nesterov (2014) that the momentum factor is ﬁxed and chosen close to 1; speciﬁcally it is proposed that

where h > 0 is the time-step (learning rate). In Wilson et al. (2016), a limiting equation for both HB and NAG of the form

u + 2 µ u + Φ(u) = 0

is derived under the assumption that λ is ﬁxed with respect to iteration number n, and dependent on the time-step h as speciﬁed in (5); convergence is obtained to order O(h1/2). Using insight from this limiting equation it is possible to choose the optimal value of µ to maximize the convergence rate in the neighborhood of a locally strongly convex objective function. Further related work is developed in Shi et al. (2018) where separate limiting equations for HB and NAG are derived both in the cases of λ given by (4) and (5), obtaining convergence to order O(h3/2). Much work has also gone into analyzing these methods in the discrete setting, without appeal to the continuous time limits, see Hu and Lessard (2017); Lessard et al. (2016), as well as in the stochastic setting, establishing how the eﬀect on the generalization error, for example, Gadat et al. (2018); Loizou and Richt arik (2017); Yang et al. (2016). In this paper, however, our focus is on the use of continuous time limits as a methodology to explain optimization algorithms. In many machine learning applications, especially for deep learning, NAG and HB are often used with a constant momentum factor λ that is chosen independently of the iteration count n (contrary to (4)) and independently of the learning rate h (contrary to (5)). In fact, popular books on the subject such as Goodfellow et al. (2016) introduce the methods in this way, and popular articles, such as He et al. (2016) to name one of many, simply state the value of the constant momentum factor used in their experiments. Widely used deep learning libraries such as Tensorﬂow Abadi et al. (2015) and Py Torch Paszke et al. (2017) implement the methods with a ﬁxed choice of momentum factor. Momentum based methods used in this way, with ﬁxed momentum, have not been carefully analyzed. We will undertake such an analysis, using ideas from numerical analysis, and in particular the concept of modiﬁed equations Griﬃths and Sanz-Serna (1986); Chartier et al. (2007) and from the theory of attractive invariant manifolds Hirsch et al. (2006); Wiggins (2013); both ideas are explained in the text Stuart and Humphries (1998). It is noteworthy that

Kovachki and Stuart

the high resolution ODE approximation described in Shi et al. (2018) may be viewed as a rediscovery of the method of modiﬁed equations. We emphasize the fact that our work is not at odds with any previous analyses of these methods, rather, we consider a setting which is widely adopted in deep learning applications and has not been subjected to continuous time analysis to date.

1.2 Our Contribution

We study momentum-based optimization algorithms for the minimization task (1), with learning rate independent momentum, ﬁxed at every iteration step, focusing on deterministic methods for clarity of exposition. Our approach is to derive continuous time approximations of the discrete algorithms; these continuous time approximations provide insights into the mechanisms at play within the discrete algorithms. We prove three such approximations. The ﬁrst shows that the asymptotic limit of the momentum methods, as learning rate approaches zero, is simply a rescaled gradient ﬂow (2). The second two approximations include small perturbations to the rescaled gradient ﬂow, on the order of the learning rate, and give insight into the behavior of momentum methods when implemented with momentum and ﬁxed learning rate. Through these approximation theorems, and accompanying numerical experiments, we make the following contributions to the understanding of momentum methods as often implemented within machine learning:

We show that momentum-based methods with a ﬁxed momentum factor, satisfy, in the continuous-time limit obtained by sending the learning rate to zero, a rescaled version of the gradient ﬂow equation (2).

We show that such methods also approximate a damped Hamiltonian system of the form (3), with small mass m (on the order of the learning rate) and constant damping γ(t) = γ; this approximation has the same order of accuracy as the approximation of the rescaled equation (2) but provides a better qualitative understanding of the ﬁxed learning rate momentum algorithm in its transient phase.

We also show that, for the approximate Hamiltonian system, the dynamics admit an exponentially attractive invariant manifold, locally representable as a graph mapping co-ordinates to their velocities. The map generating this graph describes a gradient ﬂow in a potential which is a small (on the order of the learning rate) perturbation of Φ see (21); the correction to the potential is convexifying, does not change the global minimum, and provides insight into the ﬁxed learning rate momentum algorithm beyond its initial transient phase.

We provide numerical experiments which illustrate the foregoing considerations, for simple linear test problems, and for the MNIST digit classiﬁcation problem; in the latter case we consider SGD and thereby demonstrate that the conclusions of our theory have relevance for understanding the stochastic setting as well.

Taken together our results are interesting because they demonstrate that the popular belief that (ﬁxed) momentum methods resemble the dynamics induced by (3) is misleading. Whilst it is true, the mass in the approximating equation is small and as a consequence

Continuous Time Analysis of Momentum Methods

understanding the dynamics as gradient ﬂows (2), with modiﬁed potential, is more instructive. In fact, in the ﬁrst application of HB to neural networks described in Rumelhart et al. (1986), the authors state that [their] experience has been that [one] get[s] the same solutions by setting [the momentum factor to zero] and reducing the size of [the learning rate]. However our theorems should not be understood to imply that there is no practical diﬀerence between momentum methods (with ﬁxed learning rate) and SGD. There is indeed a practical diﬀerence as has been demonstrated in numerous papers throughout the machine learning literature, and our experiments in Section 5 further conﬁrm this. We show that while these methods have the same transient dynamics, they are approximated diﬀerently. Our results demonstrate that, although momentum methods behave like a gradient descent algorithm, asymptotically, this algorithm has a modiﬁed potential. Furthermore, although this modiﬁed potential (20) is on the order of the learning rate, the fact that the learning rate is often chosen as large as possible, constrained by numerical stability, means that the correction to the potential may be signiﬁcant. Our results may be interpreted as indicating that the practical success of momentum methods stems from the fact that they provide a more stable discretization to (2) than the forward Euler method employed in SGD. The damped Hamiltonian dynamic (11), as well the modiﬁed potential, give insight into how this manifests. Our work gives further theoretical justiﬁcation for the exploration of the use of diﬀerent numerical integrators for the purposes of optimization such as those performed in Scieur et al. (2017); Betancourt et al. (2018); Zhang et al. (2018).

While our analysis is conﬁned to the non-stochastic case to simplify the exposition, the results will, with some care, extend to the stochastic setting using ideas from averaging and homogenization Pavliotis and Stuart (2008) as well as continuum analyses of SGD as in Li et al. (2017); Feng et al. (2018); indeed, in the stochastic setting, sharp uniform in time error estimates are to be expected for empirical averages Mattingly et al. (2010); Dieuleveut et al. (2017). To demonstrate that our analysis is indeed relevant in the stochastic setting, we train a deep autoencoder with mini-batching (stochastic) and verify that our convergence results still hold. The details of this experiment are given in section 5. Furthermore we also conﬁne our analysis to ﬁxed learning rate, and impose global bounds on the relevant derivatives of Φ; this further simpliﬁes the exposition of the key ideas, but is not essential to them; with considerably more analysis the ideas exposed in this paper will transfer to adaptive time-stepping methods and much less restrictive classes of Φ.

The paper is organized as follows. Section 2 introduces the optimization procedures and states the convergence result to a rescaled gradient ﬂow. In section 3 we derive the modiﬁed, second-order equation and state convergence of the schemes to this equation. Section 4 asserts the existence of an attractive invariant manifold, demonstrating that it results in a gradient ﬂow with respect to a small perturbation of Φ. In section 5, we train a deep autoencoder, showing that our results hold in a stochastic setting with Assumption 1 violated. We conclude in section 6. All proofs of theorems are given in the appendices so that the ideas of the theorems can be presented clearly within the main body of the text.

1.3 Notation

We use | | to denote the Euclidean norm on Rd. We deﬁne f : Rd Rd by f(u) := Φ(u) for any u Rd. Given parameter λ [0, 1) we deﬁne λ := (1 λ) 1.

Kovachki and Stuart

For two Banach spaces A, B, and A0 a subset in A, we denote by Ck(A0; B) the set of k-times continuously diﬀerentiable functions with domain A0 and range B. For a function u Ck(A0; B), we let Dju denote its j-th (total) Fr echet derivative for j = 1, . . . , k. For a function u Ck([0, ), Rd), we denote its derivatives by du

dt2 , etc. or equivalently by u, u, etc. To simplify our proofs, we make the following assumption about the objective function.

Assumption 1 Suppose Φ C3(Rd; R) with uniformly bounded derivatives. Namely, there exist constants B0, B1, B2 > 0 such that

Dj 1f = DjΦ Bj 1

for j = 1, 2, 3 where denotes any appropriate operator norm.

We again stress that this assumption is not key to developing the ideas in this work, but is rather a simpliﬁcation used to make our results global. Without Assumption 1, and no further assumption on Φ such as convexity, one could only hope to give local results i.e. in the neighborhood of a critical point of Φ. Such analysis could indeed be carried out (see for example Carr (2012)), but we choose not to do so here for the sake of clarity of exposition. In section 5, we give a practical example where this assumption is violated and yet the behavior is as predicted by our theory. Finally we observe that the nomenclature learning rate is now prevalent in machine learning, and so we use it in this paper; it refers to the object commonly referred to as time-step in the ﬁeld of numerical analysis.

2. Momentum Methods and Convergence to Gradient Flow

In subsection 2.1 we state Theorem 2 concerning the convergence of a class of momentum methods to a rescaled gradient ﬂow. Subsection 2.2 demonstrates that the HB and NAG methods are special cases of our general class of momentum methods, and gives intuition for proof of Theorem 2; the proof itself is given in Appendix A. Subsection 2.3 contains a numerical illustration of Theorem 2.

2.1 Main Result

The standard Euler discretization of (2) gives the discrete time optimization scheme

un+1 = un + hf(un), n = 0, 1, 2, . . . . (6)

Implementation of this scheme requires an initial guess u0 Rd. For simplicity we consider a ﬁxed learning rate h > 0. Equation (2) has a unique solution u C3([0, ); Rd) under Assumption 1 and for un = u(nh)

sup 0 nh T |un un| C(T)h;

see Stuart and Humphries (1998), for example.

Continuous Time Analysis of Momentum Methods

In this section we consider a general class of momentum methods for the minimization task (1) which can be written in the form, for some a 0 and λ (0, 1),

un+1 = un + λ(un un 1) + hf(un + a(un un 1)), n = 0, 1, 2, . . . ,

u1 = u0 + hf(u0) . (7)

Again, implementation of this scheme requires an an initial guess u0 Rd. The parameter choice a = 0 gives HB and a = λ gives NAG. In Appendix A we prove the following:

Theorem 2 Suppose Assumption 1 holds and let u C3([0, ); Rd) be the solution to

dt = (1 λ) 1 Φ(u)

u(0) = u0 (8)

with λ (0, 1). For n = 0, 1, 2, . . . let un be the sequence given by (7) and deﬁne un := u(nh). Then for any T 0, there is a constant C = C(T) > 0 such that

sup 0 nh T |un un| Ch.

Note that (8) is simply a sped-up version of (2): if v solves (2) and w solves (8) then v(t) = w((1 λ)t) for any t [0, ). This demonstrates that introduction of momentum in the form used within both HB and NAG results in numerical methods that do not diﬀer substantially from gradient descent.

2.2 Link to HB and NAG

The HB method is usually written as a two-step scheme taking the form (Sutskever et al. (2013))

vn+1 = λvn + hf(un)

un+1 = un + vn+1

with v0 = 0, λ (0, 1) the momentum factor, and h > 0 the learning rate. We can re-write this update as

un+1 = un + λvn + hf(un)

= un + λ(un un 1) + hf(un)

hence the method reads

un+1 = un + λ(un un 1) + hf(un)

u1 = u0 + hf(u0). (9)

Similarly NAG is usually written as (Sutskever et al. (2013))

vn+1 = λvn + hf(un + λvn)

un+1 = un + vn+1

Kovachki and Stuart

with v0 = 0. Deﬁne wn := un + λvn then

wn+1 = un+1 + λvn+1 = un+1 + λ(un+1 un)

un+1 = un + λvn + hf(un + λvn)

= un + (wn un) + hf(wn)

= wn + hf(wn).

Hence the method may be written as

un+1 = un + λ(un un 1) + hf(un + λ(un un 1))

u1 = u0 + hf(u0). (10)

It is clear that (9) and (10) are special cases of (7) with a = 0 giving HB and a = λ giving NAG. To intuitively understand Theorem 2, re-write (8) as

If we discretize the du/dt term using forward diﬀerences and the λdu/dt term using backward diﬀerences, we obtain

u(t + h) u(t)

h λu(t) u(t h)

h f(u(t)) f u(t) + hau(t) u(t h)

with the second approximate equality coming from the Taylor expansion of f. This can be rearranged as

u(t + h) u(t) + λ(u(t) u(t h)) + hf(u(t) + a(u(t) u(t h)))

which has the form of (7) with the identiﬁcation un u(nh).

2.3 Numerical Illustration

Figure 1 compares trajectories of the momentum numerical method (7) with the rescaled gradient ﬂow (8), for the two-dimensional problem Φ(u) = 1

2 u, Qu . We pick Q to be positive-deﬁnite so that the minimum is achieved at the point (0, 0)T and make it diagonal so that we can easily control its condition number. In particular, the condition number of Q is given as

κ = max{Q11, Q22}

min{Q11, Q22} .

We see that, as the condition number is increased, both HB and NAG exhibit more pronounced transient oscillations and are thus further away from the trajectory of (8), however, as the learning rate h is decreased, the oscillations dampen and the trajectories match more

Continuous Time Analysis of Momentum Methods

(a) HB: κ = 5

(b) HB: κ = 10

(c) HB: κ = 20

(d) NAG: κ = 5

(e) NAG: κ = 10

(f) NAG: κ = 20

Figure 1: Comparison of trajectories for HB and NAG with the gradient ﬂow (8) on the twodimensional problem Φ(u) = 1

2 u, Qu with λ = 0.9 ﬁxed. We vary the condition number of Q as well as the learning rate h.

Figure 2: The numerical rate of convergence, as a function of the learning rate h, of HB and NAG to the gradient ﬂow (8) for the problem described in Figure 1.

and more closely. This observation from Figure 1 is quantiﬁed in Figure 2 where we estimate the rate of convergence, as a function of h, which is deﬁned as

= log2 u(h) u u(h/2) u

where u(α) is the numerical solution using time-step α. The ﬁgure shows that the rate of convergence is indeed close to 1, as predicted by our theory. In summary the behavior of the momentum methods is precisely that of a rescaled gradient ﬂow, but with initial transient oscillations which capture momentum eﬀects, but disappear as the learning rate is decreased. We model these oscillations in the next section via use of a modiﬁed equation.

Kovachki and Stuart

3. Modiﬁed Equations

The previous section demonstrates how the momentum methods approximate a time rescaled version of the gradient ﬂow (2). In this section we show how the same methods may also be viewed as approximations of the damped Hamiltonian system (3), with mass m on the order of the learning rate, using the method of modiﬁed equations. In subsection 3.1 we state and discuss the main result of the section, Theorem 3. Subsection 3.2 gives intuition for proof of Theorem 3; the proof itself is given in Appendix B. And the section also contains comments on generalizing the idea of modiﬁed equations. In subsection 3.3 we describe a numerical illustration of Theorem 3.

3.1 Main Result

The main result of this section quantiﬁes the sense in which momentum methods do, in fact, approximate a damped Hamiltonian system; it is proved in Appendix B.

Theorem 3 Fix λ (0, 1) and assume that a 0 is chosen so that α := 1

2(1+λ 2a(1 λ)) is strictly positive. Suppose Assumption 1 holds and let u C4([0, ); Rd) be the solution to

dt2 + (1 λ)du

u(0) = u0, du

dt (0) = u 0. (11)

Suppose further that h (1 λ)2/2αB1. For n = 0, 1, 2, . . . let un be the sequence given by (7) and deﬁne un := u(nh). Then for any T 0, there is a constant C = C(T) > 0 such that sup 0 nh T |un un| Ch.

Theorem 2 demonstrates the same order of convergence, namely O(h), to the rescaled gradient ﬂow equation (8), obtained from (11) simply by setting h = 0. In the standard method of modiﬁed equations the limit system (here (8)) is perturbed by small terms (in terms of the assumed small learning rate) and an increased rate of convergence is obtained to the modiﬁed equation (here (11)). In our setting however, because the small modiﬁcation is to a higher derivative (here second) than appears in the limit equation (here ﬁrst order), an increased rate of convergence is not obtained. This is due to the nature of the modiﬁed equation, whose solution has derivatives that are inversely proportional to powers of h; this fact is quantiﬁed in Lemma 8 from Appendix B. It is precisely because the modiﬁed equation does not lead to a higher rate of convergence that the initial parameter u 0 is arbitrary; the same rate of convergence is obtained no matter what value it takes. It is natural to ask, therefore, what is learned from the convergence result in Theorem 3. The answer is that, although the modiﬁed equation (11) is approximated at the same order as the limit equation (8), it actually contains considerably more qualitative information about the dynamics of the system, particularly in the early transient phase of the algorithm; this will be illustrated in subsection 3.3. Indeed we will make a speciﬁc choice of u 0 in our

Continuous Time Analysis of Momentum Methods

numerical experiments, namely

dt (0) = 1 2α 2α λ + 1f(u0), (12)

to better match the transient dynamics.

3.2 Intuition and Wider Context

3.2.1 Idea Behind The Modified Equations

In this subsection, we show that the scheme (7) exhibits momentum, in the sense of approximating a momentum equation, but the size of the momentum term is on the order of the step size h. To see this intuitively, we add and subtract un un 1 to the right hand size of (7) then we can rearrange it to obtain

hun+1 2un + un 1

h2 + (1 λ)un un 1

h = f(un + a(un un 1)).

This can be seen as a second order central diﬀerence and ﬁrst order backward diﬀerence discretization of the momentum equation

dt2 + (1 λ)du

noting that the second derivative term has size of order h.

3.2.2 Higher Order Modified Equations For HB

We will now show that, for HB, we may derive higher order modiﬁed equations that are consistent with (9). Taking the limit of these equations yields an operator that agrees with with our intuition for discretizing (8). To this end, suppose Φ C b (Rd, R) and consider the ODE(s),

hk 1(1 + ( 1)kλ)

dtk = f(u) (13)

noting that p = 1 gives (8) and p = 2 gives (11). Let u C ([0, ), Rd) be the solution to (13) and deﬁne un := u(nh), u(k) n := dku

dtk (nh) for n = 0, 1, 2, . . . and k = 1, 2, . . . , p. Taylor expanding yields

un 1 = un +

k! u(k) n + hp+1I n

I n = ( 1)p+1

0 (1 s)p dp+1u

dtp+1 ((n s)h)ds.

Kovachki and Stuart

un+1 un λ(un un 1) =

k! u(k) n + λ

k! u(k) n + hp+1(I+ n λI n )

hk 1(1 + ( 1)kλ)

k! u(k) n + hp+1(I+ n λI n )

= hf(un) + hp+1(I+ n λI n )

showing consistency to order p + 1. As is the case with (11) however, the I n terms will be inversely proportional to powers of h hence global accuracy will not improve. We now study the diﬀerential operator on the l.h.s. of (13) as p . Deﬁne the sequence of diﬀerential operators Tp : C ([0, ), Rd) C ([0, ), Rd) by

hk 1(1 + ( 1)kλ)

dtk , u C ([0, ), Rd).

Taking the Fourier transform yields

F(Tpu)(ω) =

hk 1(1 + ( 1)kλ)(iω)k

where i = 1 denotes the imaginary unit. Suppose there is a limiting operator Tp T as p then taking the limit yields

F(Tu)(ω) = 1

h(eihω + λe ihω λ 1)F(u)(ω).

Taking the inverse transform and using the convolution theorem, we obtain

(Tu)(t) = 1

h F 1(eihω + λe ihω λ 1)(t) u(t)

h ( (1 + λ)δ(t) + λδ(t + h) + δ(t h)) u(t)

( (1 + λ)δ(t τ) + λδ(t τ + h) + δ(t τ h)) u(τ) dτ

h ( (1 + λ)u(t) + λu(t h) + u(t + h))

= u(t + h) u(t)

h λ u(t) u(t h)

where δ( ) denotes the Dirac-delta distribution and we abuse notation by writing its action as an integral. The above calculation does not prove convergence of Tp to T, but simply conﬁrms our intuition that (9) is a forward and backward discretization of (8).

Continuous Time Analysis of Momentum Methods

(a) HB: κ = 5

(b) HB: κ = 10

(c) HB: κ = 20

(d) NAG: κ = 5

(e) NAG: κ = 10

(f) NAG: κ = 20

Figure 3: Comparison of trajectories for HB and NAG with the Hamiltonian dynamic (11) on the two-dimensional problem Φ(u) = 1

2 u, Qu with λ = 0.9 ﬁxed. We vary the condition number of Q as well as the learning rate h.

Figure 4: The numerical rate of convergence, as a function of the learning rate h, of HB and NAG to the momentum equation (11) for the problem described in Figure 3.

3.3 Numerical Illustration

Figure 3 shows trajectories of (7) and (11) for diﬀerent values of a and h on the twodimensional problem Φ(u) = 1

2 u, Qu , varying the condition number of Q. We make the speciﬁc choice of u 0 implied by the initial condition (12). Figure 4 shows the numerical order of convergence as a function of h, as deﬁned in Section 2.3, which is near 1, matching our theory. We note that the oscillations in HB are captured well by (11), except for a slight shift when h and κ are large. This is due to our choice of initial condition which cancels the maximum number of terms in the Taylor expansion initially, but the overall rate of convergence remains O(h) due to Lemma 8. Other choices of u 0 also result in O(h)

Kovachki and Stuart

convergence and can be picked on a case-by-case basis to obtain consistency with diﬀerent qualitative phenomena of interest in the dynamics. Note also that α|a=λ < α|a=0. As a result the transient oscillations in (11) are more quickly damped in the NAG case than in the HB case; this is consistent with the numerical results. However panels (d)-(f) in Figure 1 show that (11) is not able to adequately capture the oscillations of NAG when h is relatively large. We leave for future work, the task of ﬁnding equations that are able to appropriately capture the oscillations of NAG in the large h regime.

4. Invariant Manifold

The key lessons of the previous two sections are that the momentum methods approximate a rescaled gradient ﬂow of the form (2) and a damped Hamiltonian system of the form (3), with small mass m which scales with the learning rate, and constant damping γ. Both approximations hold with the same order of accuracy, in terms of the learning rate, and numerics demonstrate that the Hamiltonian system is particularly useful in providing intuition for the transient regime of the algorithm. In this section we link the two theorems from the two preceding sections by showing that the Hamiltonian dynamics with small mass from section 3 has an exponentially attractive invariant manifold on which the dynamics is, to leading order, a gradient ﬂow. That gradient ﬂow is a small, in terms of the learning rate, perturbation of the time-rescaled gradient ﬂow from section 2.

4.1 Main Result

Deﬁne vn := (un un 1)/h (14)

noting that then (7) becomes

un+1 = un + hλvn + hf(un + havn)

and vn+1 = un+1 un

h = λvn + f(un + havn).

Hence we can re-write (7) as

un+1 = un + hλvn + hf(un + havn)

vn+1 = λvn + f(un + havn). (15)

Note that if h = 0 then (15) shows that un = u0 is constant in n, and that vn converges to (1 λ) 1f(u0). This suggests that, for h small, there is an invariant manifold which is a small perturbation of the relation vn = λf(un) and is representable as a graph. Motivated by this, we look for a function g : Rd Rd such that the manifold

v = λf(u) + hg(u) (16)

is invariant for the dynamics of the numerical method:

vn = λf(un) + hg(un) vn+1 = λf(un+1) + hg(un+1). (17)

Continuous Time Analysis of Momentum Methods

We will prove the existence of such a function g by use of the contraction mapping theorem to ﬁnd ﬁxed point of mapping T deﬁned in subsection 4.2 below. We seek this ﬁxed point in set Γ which we now deﬁne:

Deﬁnition 4 Let γ, δ > 0 be as in Lemmas 9, 10. Deﬁne Γ := Γ(γ, δ) to be the closed subset of C(Rd; Rd) consisting of γ-bounded functions:

g Γ := sup ξ Rd |g(ξ)| γ, g Γ

that are δ-Lipshitz:

|g(ξ) g(η)| δ|ξ η|, g Γ, ξ, η Rd.

Theorem 5 Fix λ (0, 1). Suppose that h is chosen small enough so that Assumption 11 holds. For n = 0, 1, 2, . . ., let un, vn be the sequences given by (15). Then there is a τ > 0 such that, for all h (0, τ), there is a unique g Γ such that (17) holds. Furthermore,

|vn λf(un) hg(un)| (λ + h2λδ)n|v0 λf(u0) hg(u0)|

where λ + h2λδ < 1.

The statement of Assumption 11, and the proof of the preceding theorem, are given in Appendix C. The assumption appears somewhat involved at ﬁrst glance but inspection reveals that it simply places an upper bound on the learning rate h, as detailed in Lemmas 9, 10. The proof of the theorem rests on the Lemmas 13, 14 and 15 which establish that the operator T is well-deﬁned, maps Γ to Γ, and is a contraction on Γ. The operator T is deﬁned, and expressed in a helpful form for the purposes of analysis, in the next subsection. In the next subsection we obtain the leading order approximation for g, given in equation (31). Theorem 5 implies that the large-time dynamics are governed by the dynamics on the invariant manifold. Substituting the leading order approximation for g into the invariant manifold (16) and using this expression in the deﬁnition (14) shows that

vn = (1 λ) 1 Φ(un) + 1

2h λ( λ a)| Φ(un)|2 , (18a)

un = un 1 h(1 λ) 1 Φ(un) + 1

2h λ( λ a)| Φ(un)|2 . (18b)

c = λ λ a + 1

we see that for large time the dynamics of momentum methods, including HB and NAG, are approximately those of the modiﬁed gradient ﬂow

dt = (1 λ) 1 Φh(u) (20)

Kovachki and Stuart

with Φh(u) = Φ(u) + 1

2hc| Φ(u)|2. (21)

To see this we proceed as follows. Note that from (20)

2(1 λ) 2 | Φ(u)|2 + O(h)

then Taylor expansion shows that, for un = u(nh),

un = un 1 + h un h2

2 un + O(h3)

= un 1 h λ Φ(un) + 1

2hc | Φ(un)|2 + 1

4h2 λ2 | Φ(un)|2 + O(h3)

where we have used that Df(u)f(u) = 1

2 | Φ(u)|2 .

Choosing c = λ( λ a + 1/2) we see that

un = un 1 h(1 λ) 1 Φ(un) + 1

2h λ( λ a)| Φ(un)|2 + O(h3). (22)

Notice that comparison of (18b) and (22) shows that, on the invariant manifold, the dynamics are to O(h2) the same as the equation (20); this is because the truncation error between (18b) and (22) is O(h3). Thus we have proved:

Theorem 6 Suppose that the conditions of Theorem 5 hold. Then for initial data started on the invariant manifold and any T 0, there is a constant C = C(T) > 0 such that

sup 0 nh T |un un| Ch2,

where un = u(nh) solves the modiﬁed equation (20) with c = λ( λ a + 1/2).

4.2 Intuition

We will deﬁne mapping T : C(Rd; Rd) C(Rd; Rd) via the equations

p = ξ + hλ λf(ξ) + hg(ξ) + hf ξ + ha λf(ξ) + hg(ξ)

λf(p) + h(Tg)(p) = λ λf(ξ) + hg(ξ) + f ξ + ha λf(ξ) + hg(ξ) . (23)

A ﬁxed point of the mapping g 7 Tg will give function g so that, under (23), identity (17) holds. Later we will show that, for g in Γ and all h suﬃciently small, ξ can be found from (23a) for every p, and that thus (23b) deﬁnes a mapping from g Γ into Tg C(Rd; Rd). We will then show that, for h suﬃciently small, T : Γ 7 Γ is a contraction.

Continuous Time Analysis of Momentum Methods

For any g C(Rd; Rd) and ξ Rd deﬁne

wg(ξ) := λf(ξ) + hg(ξ) (24)

zg(ξ) := λwg(ξ) + f ξ + hawg(ξ) . (25)

With this notation the ﬁxed point mapping (23) for g may be written

p = ξ + hzg(ξ), λf(p) + h(Tg)(p) = zg(ξ). (26)

Then, by Taylor expansion,

f ξ + ha λf(ξ) + hg(ξ) = f ξ + hawg(ξ)

= f(ξ) + ha Z 1

0 Df ξ + shawg(ξ) wg(ξ)ds

= f(ξ) + ha I(1) g (ξ)

where the last line deﬁnes I(1) g . Similarly

f(p) = f(ξ + hzg(ξ))

= f(ξ) + h Z 1

0 Df ξ + shzg(ξ) zg(ξ)ds

= f(ξ) + h I(2) g (ξ),

where the last line now deﬁnes I(2) g . Then (23b) becomes

λ f(ξ) + h I(2) g (ξ) + h(Tg)(p) = λ λf(ξ) + hλg(ξ) + f(ξ) + ha I(1) g (ξ)

and we see that (Tg)(p) = λg(ξ) + a I(1) g (ξ) λI(2) g (ξ).

In this light, we can rewrite the deﬁning equations (23) for T as

p = ξ + hzg(ξ), (29)

(Tg)(p) = λg(ξ) + a I(1) g (ξ) λI(2) g (ξ). (30)

for any ξ Rd. Perusal of the above deﬁnitions reveals that, to leading order in h,

wg(ξ) = zg(ξ) = λf(ξ), I(1) g (ξ) = I(2) g (ξ) = λDf(ξ)f(ξ).

Thus setting h = 0 in (29), (30) shows that, to leading order in h,

g(p) = λ2(a λ)Df(p)f(p). (31)

Note that since f(p) = Φ(p), Df is the negative Hessian of Φ and is thus symmetric. Hence we can write g in gradient form, leading to

2 λ2(a λ) | Φ(p)|2 . (32)

Remark 7 This modiﬁed potential (21) also arises in the construction of Lyapunov functions for the one-stage theta method see Corollary 5.6.2 in Stuart and Humphries (1998).

Kovachki and Stuart

(a) HB: un given by (15)

(b) HB: vn given by (15)

(c) HB: en given by (34)

(d) NAG: un given by (15)

(e) NAG: vn given by (15)

(f) NAG: en given by (34)

Figure 5: Invariant manifold for HB and NAG with h = 2 6 and λ = 0.9 on the twodimensional problem Φ(u) = 1

2 u, Qu , varying the condition number of Q. Panels (c), (f) show the distance from the invariant manifold for the largest condition number κ = 20.

4.3 Numerical Illustration

In Figure 5 panels (a),(b),(d),(e), we plot the components un and vn found by solving (15) with initial conditions u0 = (1, 1)T and vn = (0, 0)T in the case where Φ(u) = 1

2 u, Qu . These initial conditions correspond to initializing the map oﬀthe invariant manifold. To leading order in h the invariant manifold is given by (see equation (18))

v = (1 λ) 1 Φ(u) + 1

2h λ( λ a)| Φ(u)|2 . (33)

To measure the distance of the trajectory shown in panels (a),(b),(d),(e) from the invariant manifold we deﬁne

en = vn + (1 λ) 1 Φ(un) + 1

2h λ( λ a)| Φ(un)|2 . (34)

Panels (c),(f) show the evolution of en as well as the (approximate) bound on it found from substituting the leading order approximation of g into the following upper bound from Theorem 5: (λ + h2λδ)n|v0 λf(u0) hg(u0)|.

5. Deep Learning Example

Our theory is developed under quite restrictive assumptions, in order to keep the proofs relatively simple and to allow a clearer conceptual development. The purpose of the numerical experiments in this section is twofold: ﬁrstly to demonstrate that our theory sheds

Continuous Time Analysis of Momentum Methods

h = 20 h = 2 1 h = 2 2 h = 2 3 h = 2 4 h = 2 5 h = 2 6

GF n/a 4.3948 4.5954 5.6769 7.0049 8.6468 10.6548 HB 3.6775 4.0157 4.5429 5.6447 7.0720 8.7070 10.6848 NAG 3.2808 3.7166 4.4579 5.6087 7.0557 8.6987 10.6814 Wilson 6.7395 7.5177 8.3491 9.2543 10.2761 11.3776 12.4123 HB-µ 5.7099 6.6146 7.6202 8.6629 9.7838 11.0039 12.1743 NAG-µ 5.6867 6.6033 7.6131 8.6556 9.7783 11.0015 12.1738

Figure 6: Final training errors for the autoencoder on MNIST for six training methods over diﬀerent learning rates. GF refers to equation (35) while HB and NAG to (7) all with ﬁxed λ = 0.9.

(a) HB, NAG to (35)

(b) HB-µ, NAG-µ to (36)

Figure 7: The numerical rate of convergence for the parameters of the autoencoder, as a function of the learning rate h, of HB and NAG to (35) (a), as well as of HB-µ and NAG-µ to (36) (b).

light on a stochastic version of gradient descent applied, furthermore, to a setting in which the objective function does not satisfy the global assumptions which facilitate our analysis; and second to show that methods implemented as we use them here (with learning-rate independent momentum, ﬁxed at every step of the iteration) can out-perform other choices on speciﬁc problems.

Our numerical experiments in this section are undertaken with in the context of the example given in Sutskever et al. (2013). We train a deep autoencoder, using the architecture of Hinton and Salakhutdinov (2006) on the MNIST dataset Le Cun and Cortes (2010). Since our work is concerned only with optimization and not generalization, we present our results only on the training set of 60,000 images and ignore the testing set. We ﬁx an initialization of the autoencoder following Glorot and Bengio (2010) and use it to test every optimization method. Furthermore, we ﬁx a batch size of 200 and train for 500 epochs, not shuﬄing the data set during training so that each method sees the same realization of the noise. We use the mean-squared error as our loss function.

Kovachki and Stuart

We compare HB and NAG given by (7) to the re-scaled gradient ﬂow (8) which we discretize in the standard way to yield the numerical method

un+1 = un h (1 λ) Φ(un), (35)

hence the momentum term λ only acts to re-scale the learning rate. We do not test against equation (11) because, to discretize it faithfully, we would need to use a time-step much lower than h (because (11) contains a term of order h), but doing so would mean that we need to train for many more epochs compared to HB and NAG so that the same ﬁnal time is reached. This, in turn, implies that the methods would see diﬀerent realization of the noise. Thus, to compare them well, we would need to perform a Monte Carlo simulation, however, since we do not state any of our results in a stochastic setting, we leave this for future work. We also compare our results to those of Wilson et al. (2016) which analyze HB and NAG in the setting where Φ is µ-strongly convex and λ is given by (5) that is

They obtain the limiting equation

u + 2 µ u + Φ(u) = 0

which we discretize via a split-step method to yield

un+1 = un + 1 2 µ

1 e 2 µh vn

vn+1 = e 2 µhvn

h Φ(un+1) (36)

where we have mapped the the time-step h in HB and NAG to

h as in done in Wilson et al. (2016). We choose this discretization because it allows us to directly solve for the linear parts of the ODE (in the enlarged state-space), yielding a more accurate approximation than the forward-Euler method used to obtain (35). A detailed derivation is given in Appendix D. We will refer to the method in equation (36) as Wilson. Further we refer to equation (7) with λ given by (5) and a = 0 as HB-µ and equation (7) with λ given by (5) and a = λ as NAG-µ. Since deep neural networks are not strongly convex, there is no single optimal choice of µ; we simply set µ = 1 in our experiments. Figure 6 gives the ﬁnal training errors for each method for several learning rates. We were unable to train the autoencoder using (35) with h = 1 since λ = 0.9 implies an eﬀective learning rate of 10 for which the system blows up. In general, NAG is the best performing method for relatively large h which is an observation that is consistently made in the deep learning literature. Further, we note that as the learning rate decreases, the ﬁnal errors become closer indicating convergence to the appropriate limiting equations. Figure 6 showcases the practical eﬀectiveness of momentum methods as they provide a way of discretizing the gradient ﬂow (2) with a large eﬀective learning rate that forward Euler cannot accommodate. From this perspective, we can view momentum methods as

Continuous Time Analysis of Momentum Methods

providing a more stable discretization to gradient ﬂows in a manner illustrated by (20). Such a viewpoint informs the works Scieur et al. (2017); Betancourt et al. (2018); Zhang et al. (2018). To further illustrate the point of convergence to the limiting equation, we compute the numerical rate of convergence, deﬁned in Section 2.3, as a function of h for the neural network parameters between (35) and HB and NAG as well as between (36) and HB-µ and NAG-µ. Figure 7 gives the results. We note that this rate is around 1 as predicted by our theory while the rate for (36) is around 0.5 which is also consistent with the theory in Wilson et al. (2016).

6. Conclusion

Together, equations (8), (11) and (20) describe the dynamical systems which are approximated by momentum methods, when implemented with ﬁxed momentum, in a manner made precise by the four theorems in this paper. The insight obtained from these theorems sheds light on how momentum methods perform optimization tasks.

Acknowledgments

Both authors are supported, in part, by the US National Science Foundation (NSF) grant DMS 1818977, the US Oﬃce of Naval Research (ONR) grant N00014-17-1-2079, and the US Army Research Oﬃce (ARO) grant W911NF-12-2-0022. Both authors are also grateful to the anonymous reviewers for their invaluable suggestions which have helped to signiﬁcantly strengthen this work.

Kovachki and Stuart

Proof [of Theorem 2] Taylor expanding yields

un+1 = un + h λf(un) + O(h2)

and un = un 1 + h λf(un) + O(h2).

Hence (1 + λ)un λun 1 = un + hλ λf(un) + O(h2).

Subtracting the third identity from the ﬁrst, we ﬁnd that

un+1 ((1 + λ)un λun 1) = hf(un) + O(h2)

by noting λ λλ = 1. Similarly,

a(un un 1) = ha λf(un) + O(h2)

hence Taylor expanding yields

f(un + a(un un 1)) = f(un) + a Df(un)(un un 1)

0 (1 s)D2f(un + sa(un un 1))[un un 1]2ds

= f(un) + ha λDf(un)f(un) + O(h2).

From this, we conclude that

hf(un + a(un un 1)) = hf(un) + O(h2)

hence un+1 = (1 + λ)un λun 1 + hf(un + a(un un 1)) + O(h2).

Deﬁne the error en := un un then

en+1 = (1 + λ)en λen 1 + h (f(un + a(un un 1)) f(un + a(un un 1))) + O(h2)

= (1 + λ)en λen 1 + h Mn((1 + a)en aen 1) + O(h2)

where, from the mean value theorem, we have

0 Df s un + a(un un 1) + 1 s un + a(un un 1) ds.

Now deﬁne the concatenation En+1 := [en+1, en] R2d then

En+1 = A(λ)En + h A(a) n En + O(h2)

where A(λ), A(a) n R2d 2d are the block matrices

A(λ) := (1 + λ)I λI I 0I

, A(a) n := (1 + a)Mn a Mn 0I 0I

Continuous Time Analysis of Momentum Methods

with I Rd d the identity. We note that A(λ) has minimal polynomial

µA(λ)(z) = (z 1)(z λ)

and is hence diagonalizable. Thus there is a norm on on R2d such that its induced matrix norm m satiﬁes A(λ) m = ρ(A(λ)) where ρ : R2d 2d R+ maps a matrix to its spectral radius. Hence, since λ (0, 1), we have A(λ) m = 1. Thus

En+1 (1 + h A(a) n m) En + O(h2).

Then, by ﬁnite dimensional norm equivalence, there is a constant α > 0, independent of h, such that

1 + a a 0 0

2a2 + 2a + 1 Mn 2

where 2 denotes the spectral 2-norm. Using Assumption 1, we have

thus, letting c := α

2a2 + 2a + 1B1, we ﬁnd

En+1 (1 + hc) En + O(h2).

Then, by Gr onwall lemma,

En+1 (1 + hc)n E1 n + (1 + hc)n+1 1

= (1 + hc)n E1 n + O(h)

noting that the constant in the O(h) term is bounded above in terms of T, but independently of h. Finally, we check the initial condition

E1 = u1 u1 u0 u0

= h( λ 1)f(u0) + O(h2) 0

as desired.

Proof [of Theorem 3] Taylor expanding yields

un 1 = un h un + h2

0 (1 s)2...u((n s)h)ds.

Kovachki and Stuart

Then using equation (11)

un+1 un λ(un un 1) = h(1 λ) un + h2

2 (1 + λ) un + h3

2 (I+ n λI n )

= hf(un) + h2a(1 λ) un + h3

2 (I+ n λI n ). (37)

a(un un 1) = ha un h2

2 a un + h3

f(un + a(un un 1)) = f(un) + ha Df(un) un Df(un) h2

If n = a2 Z 1

0 (1 s)D2f(un + sa(un un 1))[un un 1]2ds.

Diﬀerentiating (11) yields

dt3 + (1 λ)d2u

dt2 = Df(u)du

hf(un + a(un un 1)) = hf(un) + h2a (hα...u n + (1 λ) un) Df(un) h3

= hf(un) + h2a(1 λ) un + h3aα...u n Df(un) h3

Rearranging this we obtain an expression for hf(un) which we plug into equation (37) to yield un+1 un λ(un un 1) = hf(un + a(un un 1)) + LTn

2 (I+ n λI n ) | {z }

O hexp (1 λ)

h3aα...u n | {z }

O hexp (1 λ)

2α n + Df(un) h3

| {z } O(h2)

h If n |{z} O(h3)

The bounds (in braces) on the four terms above follow from employing Assumption 1 and Lemma 8. From them we deduce the existence of constants K1, K2 > 0 independent of h such that

|LTn| h K1exp (1 λ)

2α n + h2K2.

We proceed similarly to the proof of Theorem 2, but with a diﬀerent truncation error structure, and ﬁnd the error satsiﬁes

En+1 (1 + hc) En + h K1exp (1 λ)

2α n + h2K2

Continuous Time Analysis of Momentum Methods

where we abuse notation and continue to write K1, K2 when, in fact, the constants have changed by use of ﬁnite-dimensional norm equivalence. Deﬁne K3 := K2/c then summing this error, we ﬁnd

En+1 (1 + hc)n E1 + h K3((1 + hc)n+1 1) + h K1

j=0 (1 + hc)jexp (1 λ)

= (1 + hc)n E1 + h K3((1 + hc)n+1 1) + h K1Sn.

Sn = exp (1 λ)

(1 + hc)n+1exp (1 λ)

2α (n + 1) 1

(1 + hc)exp 1 λ

Let T = nh then

Sn (1 + hc)n+1exp 1 λ

(1 + hc)exp 1 λ

2exp c T + 1 λ

From this we deduce that

En+1 (1 + hc)n E1 + O(h)

noting that the constant in the O(h) term is bounded above in terms of T, but independently of h. For the initial condition, we check

u1 u1 = h(u 0 f(u0)) + h2

which is O(h) by Lemma 8. Putting the bounds together we obtain

sup 0 nh T En C(T)h.

Lemma 8 Suppose Assumption 1 holds and let u C3([0, ); Rd) be the solution to

dt2 + (1 λ)du

u(0) = u0, du

dt (0) = v0

for some u0, v0 Rd and α > 0 independent of h. Suppose h (1 λ)2/2αB1 then there are constants C(1), C(2) 1 , C(2) 2 , C(3) 1 , C(3) 2 > 0 independent of h such that for any t [0, ),

| u(t)| C(1),

| u(t)| C(2) 1 h exp (1 λ)

2hα t + C(2) 2 ,

|...u(t)| C(3) 1 h2 exp (1 λ)

2hα t + C(3) 2 .

Kovachki and Stuart

One readily veriﬁes that the result of Lemma 8 is tight by considering the one-dimensional case with f(u) = u. This implies that the result of Theorem 3 cannot be improved without further assumptions. Proof [of Lemma 8] Deﬁne v := u then

hα ((1 λ)v f(u)) .

Deﬁne w := (1 λ)v f(u) hence v = (1/hα)w and u = v = λ(w + f(u)). Thus

w = (1 λ) v Df(u) u

hα w Df(u)( λ(w + f(u))).

Hence we ﬁnd

1 2 d dt|w|2 = (1 λ)

hα |w|2 λ w, Df(u)w λ w, Df(u)f(u)

hα |w|2 + λ| w, Df(u)w | + λ| w, Df(u)f(u) |

hα |w|2 + λB1|w|2 + λB0B1|w|

hα |w|2 + (1 λ)

2hα |w|2 + λB0B1|w|

2hα |w|2 + λB0B1|w|

by noting that our assumption h (1 λ)2/2αB1 implies λB1 (1 λ)/2hα. Hence

d dt|w| (1 λ)

2hα |w| + λB0B1

so, by Gr onwall lemma,

|w(t)| exp (1 λ)

2hα t |w(0)| + 2h λ2αB0B1

1 exp (1 λ)

2hα t |w(0)| + hβ1

where we deﬁne β1 := 2 λ2αB0B1. Hence

| u(t)| = | v(t)|

hαexp (1 λ)

2hα t |w(0)| + β1

= |(1 λ)v0 f(u0)|

hα exp (1 λ)

Continuous Time Analysis of Momentum Methods

thus setting C(2) 1 = |(1 λ)v0 f(u0)|/α and C(2) 1 = β1/α gives the desired result. Further,

| u(t)| = |v(t)|

λ(|w(t)| + |f(u(t))|)

λ(|w(0)| + hβ1 + B0)

hence we deduce the existence of C(1). Now deﬁne z := w then

hα z λDf(u)z + G(u, v, w)

where we deﬁne G(u, v, w) := λ(Df(u)(Df(u)v) + D2f(u)[v, w] + D2f(u)[Df(u)v, f(u)]). Using Assumption 1 and our bounds on w and v, we deduce that there is a constant C > 0 independent of h such that |G(u, v, w)| C

1 2 d dt|z|2 = (1 λ)

hα |z|2 λ z, Df(u)z + z, G(u, v, w)

hα |z|2 + λB1|z|2 + C|z|

2hα |z|2 + C|z|

as before. Thus we ﬁnd d dt|z| (1 λ)

2hα |z| + C

so, by Gr onwall lemma,

|z(t)| exp (1 λ)

2hα t |z(0)| + hβ2

where we deﬁne β2 := 2 λαC. Recall that

...u = v = 1

|z(0)| (1 λ)|(1 λ)v0 f(u0)|

hα + B1|v0|

hence we ﬁnd

|...u(t)| (1 λ)|(1 λ)v0 f(u0)

h2α2 + B1|v0|

Thus we deduce that there is a constant C(3) 1 > 0 independent of h such that

|...u(t)| C(3) 1 h2 exp (1 λ)

2hα t + C(3) 2

as desired where C(3) 2 = β2/α.

Kovachki and Stuart

Appendix C.

For the results of Section 4 we make the following assumption on the size of h. Recall ﬁrst that by Assumption 1 there are constants B0, B1, B2 > 0 such that

Dj 1f = DjΦ Bj 1

for j = 1, 2, 3.

Lemma 9 Suppose h > 0 is small enough such that

λ + h B1(a + λ λ) < 1

then there is a τ1 > 0 such that for any γ [τ1, )

(λ + h B1(a + λ λ))γ + λB0B1(a + λ) γ. (38)

Using Lemma 9 ﬁx γ [τ1, ) and deﬁne the constants

K1 := λB0 + hγ

K3 := B0 + λK1 α2 := h2(λ + ha B1),

α1 := λ 1 + h B1( λ + a(1 + h λB1)) + λ λ(B1 + h B2K3) + ha(a B2K1 + B1 λ(B1 + h B2K3) ,

α0 := a B2K1(1 + ha λB1) + λ(a B2 1 + B2K3) + λ2B1(1 + ha B1)(B1 + h B2K3). (39)

Lemma 10 Suppose h > 0 is small enough such that

α2 1 > 4α2α0, α1 < 0

then there are τ 2 > 0 such that for any δ (τ 2 , τ + 2 ]

α2δ2 + α1δ + α0 0. (40)

Using Lemma 10 ﬁx δ (τ 2 , τ + 2 ]. We make the following assumption on the size of the learning rate h which is achievable since λ (0, 1).

Assumption 11 Let Assumption 1 hold and suppose h > 0 is small enough such that the assumptions of Lemmas 9, 10 hold. Deﬁne K2 := λB1 + hδ and suppose h > 0 is small enough such that c := h(λK2 + B1(1 + ha K2)) < 1. (41)

Deﬁne constants

Q1 := λδ + a(B1K2 + B2K1(1 + ha K2)) + λ((B1 + h B2K3)(λK2 + B1(1 + ha K2)) + B2K3),

Q2 := h(a(B1 + ha B2K1) + λ(λ + ha B1)(B1 + h B2K3)),

Q3 := h(λK2 + B1(1 + ha K2)),

µ := λ + Q2 + h2(λ + ha B1)Q1

Continuous Time Analysis of Momentum Methods

Suppose h > 0 is small enough such that

Q3 < 1, µ < 1. (43)

Lastly assume h > 0 is small enough such that

λ + h2λδ < 1. (44)

Proof [of Lemma 9.] Since λ + h B1(a + λ λ) < 1 and λB0B1(a + λ) > 0 the line deﬁned by

(λ + h B1(a + λ λ))γ + λB0B1(a + λ)

will intersect the identity line at a positive γ and lie below it thereafter. Hence setting

τ1 = λB0B1(a + λ) 1 λ + h B1(a + λ λ)

completes the proof.

Proof [of Lemma 10.] Note that since α2 > 0, the parabola deﬁned by

α2δ2 + α1δ + α0

is upward-pointing and has roots

α2 1 4α2α0 2α2 .

Since α2 1 > 4α2α0, ζ R with ζ+ = ζ . Since α1 < 0, ζ+ > 0 hence setting τ + 2 = ζ+ and τ 2 = max{0, ζ } completes the proof.

The following proof refers to four lemmas whose statement and proof follow it. Proof [of Theorem 5.] Deﬁne τ > 0 as the maximum h such that Assumption 11 holds. The contraction mapping principle together with Lemmas 13, 14, and 15 show that the operator T deﬁned by (29) and (30) has a unique ﬁxed point in Γ. Hence, from its deﬁnition and equation (23b), we immediately obtain the existence result. We now show exponential attractivity. Recall the deﬁnition of the operator T namely equations (29), (30):

p = ξ + hzg(ξ)

(Tg)(p) = λg(ξ) + a I(1) g (ξ) λI(2) g (ξ).

Let g Γ be the ﬁxed point of T and set

p = un + hzg(un)

g(p) = λg(un) + a I(1) g (un) λI(2) g (un).

Kovachki and Stuart

|vn+1 λf(un+1) hg(un+1)| |vn+1 λf(un+1) hg(p)| + h|g(p) g(un+1)|

|vn+1 λf(un+1) hg(p)| + hδ|p un+1|

since g Γ. Since, by deﬁnition,

vn+1 = λvn + f(un + havn)

|vn+1 λf(un+1) hg(p)| = |λvn + f(un + havn) λf(un+1) h(λg(un) + a I(1) g (un) λI(2) g (un))|

= λ|vn λf(un) hg(un)|

by noting that

f(un + havn) = f(un) + ha I(1) g (un)

f(un+1) = f(un) + h I(2) g (un).

From deﬁnition,

un+1 = un + hλvn + hf(un + havn)

|p un+1| = |un + hzg(un) un hλvn hf(un + havn)|

= h|λ( λf(un) + hg(un)) + f(un + havn) λvn f(un + havn)|

= hλ|vn λf(un) hg(un)|.

|vn+1 λf(un+1) hg(un+1)| (λ + h2λδ)|vn λf(un) hg(un)|

as desired. By Assumption 11, λ + h2λδ < 1.

The following lemma gives basic bounds which are used in the proof of Lemmas 13, 14, 15.

Continuous Time Analysis of Momentum Methods

Lemma 12 Let g, q Γ and ξ, η Rd then the quantities deﬁned by (24), (25), (27), (28) satisfy the following:

|wg(ξ)| K1,

|wg(ξ) wg(η)| K2|ξ η|,

|wg(ξ) wq(ξ)| h|g(ξ) q(ξ)|,

|zg(ξ)| K3,

|zg(ξ) zg(η)| (λK2 + B1 (1 + ha K2)) |ξ η|,

|zg(ξ) zq(ξ)| h (λ + ha B1) |g(ξ) q(ξ)|,

|I(1) g (ξ)| B1K1,

|I(1) g (ξ) I(1) g (η)| (B1K2 + B2K1(1 + ha K2))|ξ η|,

|I(1) g (ξ) I(1) q (ξ)| h(B1 + ha B2K1)|g(ξ) q(ξ)|,

|I(2) g (ξ)| B1K3

|I(2) g (ξ) I(2) g (η)| ((B1 + h B2K3)(λK2 + B1(1 + ha K2)) + B2K3)|ξ η|,

|I(2) g (ξ) I(2) q (ξ)| h(λ + h B1a)(B1 + h B2K3)|g(ξ) q(ξ)|.

Proof These bounds relay on applications of the triangle inequality together with boundedness of f and its derivatives as well as the fact that functions in Γ are bounded and Lipschitz. To illustrate the idea, we will prove the bounds for wg, wq, I(1) g , and I(1) q . To that end,

|wg(ξ)| = | λf(ξ) + hg(ξ)|

λ|f(ξ)| + h|g(ξ)|

establishing the ﬁrst bound. For the second,

|wg(ξ) wg(η)| λ|f(ξ) f(η)| + h|g(ξ) g(η)|

λB1|ξ η| + hδ|ξ η|

as desired. Finally,

|wg(ξ) wq(ξ)| = | λf(ξ) + hg(ξ) λf(ξ) hq(ξ)|

= h|g(ξ) q(ξ)|

as desired. We now turn to the bounds for I(1) g , I(1) q ,

|I(1) g (ξ)| Z 1

0 |Df(ξ + shawg(ξ))||wg(ξ)|ds

Kovachki and Stuart

establishing the ﬁrst bound. For the second bound,

|I(1) g (ξ) I(1) g (η)| Z 1

0 |Df(ξ + shawg(ξ))wg(ξ) Df(η + shawg(η))wg(ξ)|ds

0 |Df(η + shawg(η))wg(ξ) Df(η + shawg(η))wg(η)|ds

0 (|ξ η| + sha|wg(ξ) wg(η)|)ds + B1|wg(ξ) wg(η)|

K1B2(|ξ η| + ha K2|ξ η|) + B1K2|ξ η|

= (B1K2 + B2K1(1 + ha K2))|ξ η|

as desired. Finally

|I(1) g (ξ) I(1) q (ξ)| Z 1

0 |Df(ξ + shawg(ξ))wg(ξ) Df(ξ + shawg(ξ))wq(ξ)|ds

0 |Df(ξ + shawg(ξ))wq(ξ) Df(ξ + shawq(ξ))wq(ξ)|ds

0 |wg(ξ) wq(ξ)|ds + K1B2

0 |ξ + shawg(ξ) ξ shawq(ξ)|ds

h B1|g(ξ) q(ξ)| + h2a B2K1|g(ξ) q(ξ)|

= h(B1 + ha B2K1)|g(ξ) q(ξ)|

as desired. The bounds for zg, zq, I(2) g , and I(2) q follow similarly.

We also need the following three lemmas:

Lemma 13 Suppose Assumption 11 holds. For any g Γ and p Rd there exists a unique ξ Rd satisfying (29).

Lemma 14 Suppose Assumption 11 holds. The operator T deﬁned by (30) satisﬁes T : Γ Γ.

Lemma 15 Suppose Assumption 11 holds. For any g1, g2 Γ, we have

Tg1 Tg2 Γ µ g1 g2 Γ

where µ < 1.

Now we prove these three lemmas. Proof [of Lemma 13.] Consider the iteration of the form

ξk+1 = p hzg(ξk).

For any two sequences {ξk}, {ηk} generated by this iteration we have, by Lemma 12,

|ξk+1 ηk+1| h|zg(ηk) zg(ξk)|

h(λK2 + B1(1 + ha K2))|ξk ηk|

Continuous Time Analysis of Momentum Methods

which is a contraction by (41).

Proof [of Lemma 14.] Let g Γ and p Rd then by Lemma 13 there is a unique ξ Rd

such that (29) is satisﬁed. Then

|(Tg)(p)| λ|g(ξ)| + a|I(1) g (ξ)| + λ|I(2) g (ξ)|

λγ + a B1( λB0 + hγ) + λB1(λ( λB0 + hγ) + B0)

= (λ + h B1(a + λ λ))γ + λB0B1(a + λ)

with the last inequality following from (38). Let p1, p2 Rd then, by Lemma 13, there exist ξ1, ξ2 Rd such that (29) is satisﬁed with p = {p1, p2}. Hence, by Lemma 12,

|(Tg)(p1) (Tg)(p2)| λ|g(ξ1) g(ξ2)| + a|I(1) g (ξ1) I(1) g (ξ2)| + λ|I(2) g (ξ1) I(2) g (ξ2)|

where we deﬁne

K := λδ + a(B1K2 + B2K1(1 + ha K2)) + λ((B1 + h B2K3)(λK2 + B1(1 + ha K2)) + B2K3).

Now, using (29) and the proof of Lemma 13,

|ξ1 ξ2| |p1 p2| + h|zg(ξ1) zg(ξ2)|

|p1 p2| + c|ξ1 ξ2|.

Since c < 1 by (41), we obtain

|ξ1 ξ2| 1 1 c|p1 p2|

thus |(Tg)(p1) (Tg)(p2)| K 1 c|p1 p2| δ|p1 p2|.

To see the last inequality, we note that

K 1 c δ K δ(1 c) 0

and K δ(1 c) = α2δ2 + α1δ + α0 by (39) hence (40) gives the desired result.

Proof [of Lemma 15.] By Lemma 13, for any p Rd and g1, g2 Γ, there are ξ1, ξ2 Rd

p = ξj + hzgj(ξj)

(Tgj)(p) = λgj(ξj) + a I(1) gj (ξj) λI(2) gj (ξj)

Kovachki and Stuart

for j = 1, 2. Then

|(Tg1)(p) (Tg2)(p)| λ|g1(ξ1) g2(ξ2)| + a|I(1) g1 (ξ1) I(1) g2 (ξ2)| + λ|I(2) g1 (ξ1) I(2) g2 (ξ2)|.

|g1(ξ1) g2(ξ2)| = |g1(ξ1) g2(ξ2) g2(ξ1) + g2(ξ1)|

|g1(ξ1) g2(ξ1)| + δ|ξ1 ξ2|.

Similarly, by Lemma 12,

|I(1) g1 (ξ1) I(1) g2 (ξ2)| = |I(1) g1 (ξ1) I(1) g2 (ξ2) I(1) g2 (ξ1) + I(1) g2 (ξ1)|

|I(1) g1 (ξ1) I(1) g2 (ξ1)| + |I(1) g2 (ξ1) I(1) g2 (ξ2)|

h(B1 + ha B2K1)|g1(ξ1) g2(ξ1)| + (B1K2 + B2K1(1 + ha K2))|ξ1 ξ2|

|I(2) g1 (ξ1) I(2) g2 (ξ2)| = |I(2) g1 (ξ1) I(2) g2 (ξ2) I(2) g2 (ξ1) + I(2) g2 (ξ1)|

|I(2) g1 (ξ1) I(2) g2 (ξ1)| + |I(2) g2 (ξ1) I(2) g2 (ξ2)|

h(λ + h B1a)(B1 + h B2K3)|g1(ξ1) g2(ξ1)|+

+ ((B1 + h B2K3)(λK2 + B1(1 + ha K2)) + B2K3)|ξ1 ξ2|

Putting these together and using (42), we obtain

|(Tg1)(p) (Tg2)(p)| (λ + Q2)|g1(ξ1) g2(ξ1)| + Q1|ξ1 ξ2|.

Now, by Lemma 12,

|ξ1 ξ2| h|zg1(ξ1) zg2(ξ2) zg2(ξ1) + zg2(ξ1)|

h(|zg1(ξ1) zg2(ξ1)| + |zg2(ξ1) zg2(ξ2)|)

h2(λ + ha B1)|g1(ξ) g2(ξ1)| + h(λK2 + B1(1 + ha K2))|ξ1 ξ2|

= h2(λ + ha B1)|g1(ξ) g2(ξ1)| + Q3|ξ1 ξ2|

using (42). Since, by (43), Q3 < 1, we obtain

|ξ1 ξ2| h2(λ + ha B1)

1 Q3 |g1(ξ1) g2(ξ1)|

|(Tg1)(p) (Tg2)(p)| λ + Q2 + h2(λ + ha B1)Q1

|g1(ξ1) g2(ξ1)|

= µ|g1(ξ1) g2(ξ1)|

by (42). Taking the supremum over ξ1 then over p gives the desired result. Since µ < 1 by (43), we obtain that T is a contraction on Γ.

Continuous Time Analysis of Momentum Methods

We consider the equation

u + 2 µ u + Φ(u) = 0

u(0) = u0, u(0) = v0.

Set v = u then we have u v

= v 2 µv Φ(u)

Deﬁne the maps

f1(u, v) := v 2 µv

, f2(u, v) := 0 Φ(u)

= f1(u, v) + f2(u, v).

We ﬁrst solve the system u v

= f1(u, v).

Clearly v(t) = e 2 µtv0

u(t) = u0 + Z t

0 e 2 µsv0 ds

= u0 + 1 2 µ

1 e 2 µt v0.

This gives us the ﬂow map

ψ1(u, v; t) =

" u + 1 2 µ 1 e 2 µt v e 2 µtv

We now solve the system u v

= f2(u, v).

Clearly u(t) = u0

hence v(t) = v0 t Φ(u0).

This gives us the ﬂow map

ψ2(u, v; t) = u v t Φ(u)

Kovachki and Stuart

The composition of the ﬂow maps is then

(ψ2 ψ1)(u, v; t) =

" u + 1 2 µ 1 e 2 µt v

e 2 µtv t Φ u + 1 2 µ 1 e 2 µt v #

Mapping t to the time-step

h gives the numerical method (36).

Mart ın Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeﬀrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoﬀrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Man e, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Vi egas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. Tensor Flow: Large-scale machine learning on heterogeneous systems, 2015. URL https://www.tensorflow.org/. Software available from tensorﬂow.org.

Dimitris Bertsimas, John Tsitsiklis, et al. Simulated annealing. Statistical science, 8(1): 10 15, 1993.

Michael Betancourt, Michael I. Jordan, and Ashia C. Wilson. On symplectic optimization, 2018.

Jack Carr. Applications of centre manifold theory, volume 35. Springer Science & Business Media, 2012.

Philippe Chartier, Ernst Hairer, and Gilles Vilmart. Numerical integrators based on modiﬁed diﬀerential equations. Mathematics of computation, 76(260):1941 1953, 2007.

Aymeric Dieuleveut, Alain Durmus, and Francis Bach. Bridging the gap between constant step size stochastic gradient descent and markov chains. ar Xiv preprint ar Xiv:1707.06386, 2017.

John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. J. Mach. Learn. Res., 12:2121 2159, July 2011. ISSN 1532-4435. URL http://dl.acm.org/citation.cfm?id=1953048.2021068.

Yuanyuan Feng, Lei Li, and Jian-Guo Liu. Semigroups of stochastic gradient descent and online principal component analysis: properties and diﬀusion approximations. Communications in Mathematical Sciences, 16(3):777 789, 2018.

S ebastien Gadat, Fabien Panloup, Soﬁane Saadane, et al. Stochastic heavy ball. Electronic Journal of Statistics, 12(1):461 529, 2018.

Continuous Time Analysis of Momentum Methods

Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. In Readings in computer vision, pages 564 584. Elsevier, 1987.

Xavier Glorot and Yoshua Bengio. Understanding the diﬃculty of training deep feedforward neural networks. In In Proceedings of the International Conference on Artiﬁcial Intelligence and Statistics (AISTATS 10). Society for Artiﬁcial Intelligence and Statistics, 2010.

Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016. http://www.deeplearningbook.org.

DF Griﬃths and JM Sanz-Serna. On the scope of the method of modiﬁed equations. SIAM Journal on Scientiﬁc and Statistical Computing, 7(3):994 1008, 1986.

Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770 778, 2016.

Geoﬀrey Hinton and Ruslan Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504 507, 2006.

Morris W Hirsch, Charles Chapman Pugh, and Michael Shub. Invariant manifolds, volume 583. Springer, 2006.

Bin Hu and Laurent Lessard. Dissipativity theory for nesterov s accelerated method. In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1549 1557. JMLR. org, 2017.

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. Co RR, abs/1412.6980, 2014. URL http://arxiv.org/abs/1412.6980.

Harold J Kushner. Asymptotic global behavior for stochastic approximation and diﬀusions with slowly decreasing noise eﬀects: global minimization via monte carlo. SIAM Journal on Applied Mathematics, 47(1):169 185, 1987.

Harold Joseph Kushner and Dean S Clark. Stochastic approximation methods for constrained and unconstrained systems, volume 26. Springer Science & Business Media, 2012.

Yann Le Cun and Corinna Cortes. MNIST handwritten digit database. 2010. URL http: //yann.lecun.com/exdb/mnist/.

Yann Le Cun, Yoshua Bengio, and Geoﬀrey E. Hinton. Deep learning. Nature, 521 (7553):436 444, 2015. doi: 10.1038/nature14539. URL https://doi.org/10.1038/ nature14539.

Laurent Lessard, Benjamin Recht, and Andrew Packard. Analysis and design of optimization algorithms via integral quadratic constraints. SIAM Journal on Optimization, 26(1): 57 95, 2016.

Kovachki and Stuart

Qianxiao Li, Cheng Tai, and Weinan E. Stochastic modiﬁed equations and adaptive stochastic gradient algorithms. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 2101 2110, International Convention Centre, Sydney, Australia, 06 11 Aug 2017. PMLR.

Nicolas Loizou and Peter Richt arik. Linearly convergent stochastic heavy ball method for minimizing generalization error. ar Xiv preprint ar Xiv:1710.10737, 2017.

Jonathan C Mattingly, Andrew M Stuart, and Michael V Tretyakov. Convergence of numerical time-averaging and stationary measures via poisson equations. SIAM Journal on Numerical Analysis, 48(2):552 577, 2010.

Yurii Nesterov. A method of solving a convex programming problem with convergence rate o(1/k2). Soviet Mathematics Doklady, 27(2):372 376, 1983.

Yurii Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Springer Publishing Company, Incorporated, 1 edition, 2014. ISBN 1461346916, 9781461346913.

Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary De Vito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In NIPS-W, 2017.

Grigorios Pavliotis and Andrew Stuart. Multiscale Methods: Averaging and Homogenization, volume 53. 01 2008. doi: 10.1007/978-0-387-73829-1.

Boris Polyak. Some methods of speeding up the convergence of iteration methods. Ussr Computational Mathematics and Mathematical Physics, 4:1 17, 12 1964. doi: 10.1016/ 0041-5553(64)90137-5.

Boris T. Polyak. Introduction to optimization. New York: Optimization Software, Inc., 1987.

Ning Qian. On the momentum term in gradient descent learning algorithms. Neural Netw., 12(1):145 151, January 1999. ISSN 0893-6080. doi: 10.1016/S0893-6080(98)00116-6. URL http://dx.doi.org/10.1016/S0893-6080(98)00116-6.

Herbert Robbins and Sutton Monro. A stochastic approximation method. Ann. Math. Statist., 22(3):400 407, 09 1951. doi: 10.1214/aoms/1177729586. URL https://doi. org/10.1214/aoms/1177729586.

D. E. Rumelhart, G. E. Hinton, and R. J. Williams. Parallel distributed processing: Explorations in the microstructure of cognition, vol. 1. chapter Learning Internal Representations by Error Propagation, pages 318 362. MIT Press, Cambridge, MA, USA, 1986. ISBN 0-262-68053-X. URL http://dl.acm.org/citation.cfm?id=104279.104293.

Damien Scieur, Vincent Roulet, Francis Bach, and Alexandre d Aspremont. Integration methods and optimization algorithms. In Advances in Neural Information Processing Systems, pages 1109 1118, 2017.

Continuous Time Analysis of Momentum Methods

Bin Shi, Simon S Du, Michael I Jordan, and Weijie J Su. Understanding the acceleration phenomenon via high-resolution diﬀerential equations. ar Xiv preprint ar Xiv:1810.08907, 2018.

Andrew Stuart and Anthony R Humphries. Dynamical systems and numerical analysis, volume 2. Cambridge University Press, 1998.

MA Styblinski and T-S Tang. Experiments in nonconvex optimization: stochastic approximation with function smoothing and simulated annealing. Neural Networks, 3(4):467 483, 1990.

Weijie Su, Stephen Boyd, and Emmanuel Candes. A diﬀerential equation for modeling nesterov s accelerated gradient method: Theory and insights. In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 27, pages 2510 2518. Curran Associates, Inc., 2014. URL http://papers.nips.cc/paper/ 5322-a-differential-equation-for-modeling-nesterovs-accelerated-gradient-method-theory-an pdf.

Ilya Sutskever, James Martens, George Dahl, and Geoﬀrey Hinton. On the importance of initialization and momentum in deep learning. In Proceedings of the 30th International Conference on International Conference on Machine Learning - Volume 28, ICML 13, pages III 1139 III 1147. JMLR.org, 2013. URL http://dl.acm.org/citation.cfm? id=3042817.3043064.

T. Tieleman and G. Hinton. Lecture 6.5 Rms Prop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural Networks for Machine Learning, 2012.

Stephen Wiggins. Normally hyperbolic invariant manifolds in dynamical systems, volume 105. Springer Science & Business Media, 2013.

Ashia C. Wilson, Benjamin Recht, and Michael I. Jordan. A lyapunov analysis of momentum methods in optimization. Co RR, abs/1611.02635, 2016. URL http://arxiv.org/abs/ 1611.02635.

Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nati Srebro, and Benjamin Recht. The marginal value of adaptive gradient methods in machine learning. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4148 4158. Curran Associates, Inc., 2017. URL http://papers.nips.cc/paper/ 7003-the-marginal-value-of-adaptive-gradient-methods-in-machine-learning. pdf.

Tianbao Yang, Qihang Lin, and Zhe Li. Uniﬁed convergence analysis of stochastic momentum methods for convex and non-convex optimization. ar Xiv preprint ar Xiv:1604.03257, 2016.

Kovachki and Stuart

Jingzhao Zhang, Aryan Mokhtari, Suvrit Sra, and Ali Jadbabaie. Direct runge-kutta discretization achieves acceleration. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS 18, page 3904 3913, Red Hook, NY, USA, 2018. Curran Associates Inc.