{'fraction_non_alphanumeric': 0.05681629948649992, 'fraction_numerical': 0.016730163988736127, 'mean_word_length': 4.097952290479206, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic behavior for a gradient algorithm with energy and momentum

This paper investigates a novel gradient algorithm, using both energy and momentum (called AGEM), for solving general non-convex optimization problems. The solution properties of the AGEM algorithm, including uniformly boundedness and convergence to critical points, are examined. The dynamic behavior is studied through analysis of a high-resolution ODE system. Such ODE system is nonlinear and obtained by taking the limit of the discrete scheme while keeping the momentum effect through a rescale of the momentum parameter. In particular, we show global well-posedness of the ODE system, time-asymptotic convergence of solution trajectories, and further establish a linear convergence rate for objective functions satisfying the Polyak-Lojasiewicz condition.


Introduction
In this paper we will be considering the unconstrained optimization problem, (1.1) min θ∈R n f (θ), where the objective function f : R n → R is assumed to be differentiable and bounded from below, i.e., f + c > 0 for some c ∈ R. The most prominent applications such as training machine learning models have involved large-scale problems in which the objective function is typically non-convex and in the form of a finite sum over terms associated with individual data. In this setting, first-order gradient methods such as gradient descent (GD) and its variants are often applied since they are computationally efficient and has satisfying performance [5]; consequently, there has been much effort on the theory and practices of speeding up the convergence of first-order methods, among which adaptive step size [26,12] and momentum [22] are two widely-used techniques.
The essential idea in many adaptive gradient methods is using the past gradient values to adjust the effective step size, with renowned works including Adam [12] and RMSProp [26]. AEGD (adaptive gradient descent with energy) is another type of adaptive gradient methods that was first introduced by us in [16]. In this work, we study a variant of AEGD with momentum (termed as AGEM), which takes the following form: starting form θ 0 and r 0 = F (θ 0 ) with F = √ f + c, inductively define v k = βv k−1 + (1 − β)∇F (θ k ), v −1 = 0, (1.2a) where β ∈ (0, 1) is the momentum parameter. This when taking β = 0 reduces to the original AEGD in [16]. Distinguished from other adaptive gradient methods, AEGD and its variants use an energy variable r k to adjust the effective learning rate. The AGEM (1.2) is comparable to the SGEM (Stochastic Gradient with energy and momentum) algorithms introduced in [13]. In sharp contrast to GD, all these energy-adaptive gradient methods feature the unconditional energy stability property, irrespective of the base step size η. The excellent performance of AEGD type schemes has been illustrated with various optimization tasks [16,13]. On the theoretical side, the convergence rates obtained so far depend on (kr k ) −1 , instead of k −1 ; this could degrade the convergence behavior when the energy variable decays too fast. This issue consists of one of the main challenges in better understanding this class of energy-adaptive gradient methods.
On the other hand, the connection between ODEs and numerical optimization is often established via taking step sizes to be very small so that the trajectory or solution path converges to a curve modeled by an ODE. The well-established theory and dynamic properties of ODEs can provide deeper insights into optimization; see e.g., [1,3,2,25].
In this work, we therefore derive an ODE system for (1.2), which is the exact limit of (1.2) by taking small step sizes while keeping = βη/(1 − β) fixed. This ODE system with unknown U := (v, r, θ) reads v = −v + ∇F (θ), for t > 0, with initial conditions U (0) = (0, F (θ 0 ), θ 0 ); here, θ 0 is the starting point in the AGEM scheme,U denotes the time derivative. This work is the first to model AEGDtype schemes by ODE systems. With this connection, the convergence properties of (1.2) will be studied from the perspective of continuous dynamics for general smooth objective functions.
For a gradient method, the geometric property of f often affects the convergence rates. For non-convex f , we consider an old condition originally introduced by Polyak [20], who showed that it is a sufficient condition for GD to achieve a linear convergence rate; see [11] for a recent convergence study under this condition. Such condition is often called the Polyak-Lojasiewicz (PL) inequality in the literature since it is a special case of the inequality introduced by Lojasiewicz in 1963 [17]. A number of important nonconvex problems in machine learning, including some classes of neural networks, have been recently shown to satisfy the PL condition [4,7,24,15]. It is believed that the PL condition provides a relevant and attractive setting for many machine learning problems, particularly in the over-parametrized regime.

Contributions
We study the convergence behavior of both the discrete scheme (1.2) and its continuous counterpart (1.3). Our main contributions are summarized as follows.
1. For (1.2), we prove that the iterates are uniformly bounded and converge to the set of critical points of the objective function when the step size is suitably small, also lim k→∞ r k = r * > 0.
2. We derive (1.3) as a high resolution ODE system of the discrete scheme (1.2). For this ODE system, we first show the global well-posedness by constructing a suitable Lyapunov function. We then establish the time-asymptotic convergence of the solution trajectory toward a local minimum of the objective function using the Lasalle invariance principle. Moreover, a positive lower bound of the energy variable is provided.
3. For objective functions that satisfy the PL condition (see Section 2.2), we are able to show the linear convergence rate for some α > 0, where θ * is the global minimizer of f (not necessarily convex).
For non-convex optimization problems, obtaining convergence rates is often a challenging task. The technique used in this work to obtain the linear convergence rate for system (1.3) is inspired by that used in [2], where the authors derived a linear convergence rate for the Heavy Ball dynamical system, that is (1.4) with a(t) being a constant. For the more complicated nonlinear system (1.3), we are able to construct a class of control functions which play a similar role to that in [2]. The proof of the convergence rate results for (1.3) is significantly different and more involved. In particular, we explore a subtle optimization argument to identify an optimal control function to achieve the desired decay rate.
Further related works A recent line of research on the theoretical analysis of accelerated gradient algorithms has been on the study of their continuous-time limits (ODEs). With well-established theory and techniques such as Lyapunov functions, ODEs prove to be very useful when studying the convergence behavior of a discrete algorithm. The convergence behavior of the heavy ball method has been studied via the following second-order ODE: (1.4)θ(t) + a(t)θ + ∇f (θ) = 0, θ(0) = θ 0 ,θ = 0 with a(t) being a smooth, positive and non-increasing function. The Lyapunov function E = 1 2 |θ| 2 + f (θ) for (1.4) plays an essential role in analysis for both convex setting [1] and non-convex setting [3]. The Nesterov algorithm is shown to be connected to (1.4) with a(t) = 3 t in [25], where the authors recover the optimal convergence rate of O(1/t 2 ) of the discrete scheme in the convex setting [19]. A non-autonomous system of differential equations was derived and analyzed in [8] as the continuous time limit of adaptive optimization methods including Adam and RMSprop. In addition to convergence behavior, such dynamical approach also allows for better understanding of qualitative features of discrete algorithms. For example, [23] investigated the difference between the acceleration phenomena of the heavy ball method and the Nesterov algorithm by analyzing some high-resolution ODEs. By analyzing the sign GD of formθ = − ∇f (θ) |∇f (θ)| and its variants, [18] identified three behavior patterns of Adam and RMSprop: fast initial progress, small oscillations, and large spikes. From a continuous dynamical system satisfying nice converging properties, one may discretize it to obtain novel algorithms; for example, a new second-order inertial optimization method (called INNA algorithm) is obtained from a continuous inclusion system when the objective function is non-smooth [6]. To our best knowledge, (1.3) is different from the existing ODEs for optimization methods, and also subtle to analyze.
The rest of the paper is organized as follows. In Section 2, we begin with the problem setup, recall the existing energy-adaptive gradient methods and propose a new method. The PL condition and related properties are discussed. In Section 3, we analyze the solution properties of the proposed method. In Section 4, we derive a continuous dynamical system and present a convergence result to show that the dynamical system is faithful model for the discrete scheme. In Section 5, we analyze the global existence and asymptotic convergence of the solution to the dynamical system. Finally, we provide a convergence rate for objective functions satisfy PL condition in Section 6.

Energy-adaptive gradient algorithms
For the objective function in (1.1), we make the following assumptions: (2) f is bounded from below so that F = √ f + c is well-defined for some c ∈ R.
Throughout this work, we denote We use ∇f and ∇ 2 f to denote the gradient and Hessian of f ; and ∂ i f to denote the i-th element of ∇f . For x ∈ R n , we use |x| to denote its length, i.e. |x| = x 2 1 + ... + x 2 n . For x, y ∈ R n , we use x ⊗ y to denote their outer product. For an integer m, we use [m] to denote the set {1, 2, ..., m}.
Under Assumption 2.1, note that f + c = F 2 , we have ∇f = 2F ∇F and The fact that F ≥ F * > 0 implies that we also have F ∈ C 2 (R n ), and F is coercive.
2.1. Energy-adaptive gradient algorithms. We recall the AEGD (Adaptive gradient descent with energy) method This is the global and non-stochastic version of those introduced in [16]. The most striking feature is that it is unconditionally energy stable in the sense that r k , which serves as an approximation to F (θ k ) = f (θ k ) + c, is strictly decreasing in k unless θ k+1 = θ k . A natural way to speed up AEGD is to add momentum. The following scheme (called SGEM) introduced in [13] is an attempt in this direction: where β ∈ (0, 1) is the momentum parameter. It is shown to achieve a faster convergence speed than AEGD while generalizing better than stochastic gradient descent with momentum in training benchmark deep neural networks [13].
In this work, we propose a new variant of (2.2), called AGEM: All the above three methods feature the same unconditional energy stability property, stated in the following.  .3) are all unconditionally energy stable in the sense that for any step size η > 0, we have This implies that r k is strictly decreasing and converges to r * ≥ 0 as k → ∞, and lim k→∞ |θ k+1 − θ k | 2 = 0. Moreover, The excellent performance of (2.2) has been illustrated in [13] with various optimization tasks. Here we present a simple experiment on testing a non-convex function as shown in Figure 1, from which we see that AGEM exhibits faster convergence than SGEM, AEGD and GDM (GD with momentum).
The main objective of this work is to present a detailed analysis of AGEM using both its discrete and continuous formulations. These results are expected to also hold true for (2.2).
We are particularly interested in the convergence properties of AGEM. This will be studied from the perspective of continuous dynamics for general smooth objective functions. Furthermore, based on a structural condition on f -the so-called Polyak-Lojasiewicz (PL) condition, we derive a convergence rate in Section 6.
2.2. PL condition for non-convex objectives. Here we present a short account on the PL condition and its relevance to the loss function in the context of deep neural networks.
For a differentiable function f : R n → R with argminf (θ) = ∅, we say f satisfies the PL condition if there exists some constant µ > 0, such that the following inequality holds for any θ * ∈ argminf (θ): This means that a local minimizer θ * is also a global minimizer. We note that strongly convex functions satisfy the PL condition, while a function satisfies the PL condition is not necessary to be convex. For example, the function is not convex but satisfies the PL condition with µ = 1 32 , f * = 0. As we mentioned in the introduction, PL condition has attracted increasing attention in connection with the convergence of training deep neural network tasks. Here we conduct an informal investigation on how a loss function commonly used in training deep learning tasks can satisfy the PL condition. Consider the following loss function are training data points, g is a deep neural network prediction parameterized by θ. The gradient of such loss function is where H is a m × m matrix with (i, j) entry defined by From (2.7) it can be seen that as long as the smallest eigenvalue of H is lower bounded by µ, we have This is indeed the PL condition. In fact, using over-parameterization and random initialization, [9] has proved that for a two-layer neural network with rectified linear unit (ReLU) activation, the smallest eigenvalue of H is indeed strictly positive and it plays an essential role in deriving the linear convergence rate of gradient descent in [9]. A similar convergence result for implicit networks is also obtained in [10].

Solution properties of the AGEM algorithm
In this section, we prove that the iterates generated by (2.3) are uniformly bounded and convergent when η is suitably small.
we have the following results: (1) (Uniformly boundedness) θ k and v k are bounded for all k ≥ 1. Moreover, (2) (Lower bound on energy) The lower bound of r k is strictly positive: (3) (Convergence) As k → ∞, we have where = ηβ 1−β , then we can claim: in which, where r k+1 < r k is used in the last inequality. Connecting (1) To show the uniformly boundedness, we take a summation on (3.2) from 0 to k and ignore the negative terms on the RHS to have Using Q 0 = F (θ 0 ) and the bound in (2.5), we further get ).
Now we show that when η is suitably small, {θ k } is uniformly bounded. This further ensures the boundedness of {v k }, since To proceed, we use the notation: andΣ k for the convex hull * of Σ k . By the assumption, Σ 2 ⊂ Σ 3 are bounded domains. We want to identify conditions on η such that {θ k } ∈ Σ 2 for all k ≥ 1. This can be carried out by induction. Obviously θ 0 ∈ Σ 2 . It suffices to show that under suitable conditions on η, we have The remaining task is to find the condition on η, for which we argue in two steps: On the other hand, |∇F (θ)|. * The convex hull of a given set S is the (unique) minimal convex set containing S. This ensures (3.9) if we take This gives .
Remark 3.2. If f is assumed to be L-smooth, then L F is uniformly bounded. In such case, the boundedness of the solution is achieved without any restriction on η, as can be seen from (3.7). However, the L-smoothness assumption would exclude a large class of objective functions such as |θ| 4 .

Discrete to continuous-time-limits
Scheme (2.3) can be studied by considering the limiting ODEs obtained by taking the step size η → 0. However, different ODEs can be obtained when the hyper-parameters are scaled differently. Here we shall derive an ODE system which is simpler for analysis.
For any k ≥ 0, let t k = kη, and assume v k = v(t k ), r k = r(t k ), θ k = θ(t k ) for some sufficiently smooth curve (v(t), r(t), θ(t)). Performing the Taylor expansion on v k−1 in powers of η, we get Similarly, from (2.3c) we have We first discard O(η 2 ) and keep O(η) terms in (4.1) (4.2), (4.3), which leads to the following ODE system: Let η → 0, then (4.4) reduces to the following ODE system: This says that the ODE system (4.5) is equivalent to the gradient flow (4.6).
To study the momentum effect, we keep = βη 1−β unchanged while letting η → 0. This way, we obtain the following system of ODEs: Compared with (4.7), the high-resolution ODEs (4.4) are more accurate continuous-time counterparts for the corresponding discrete algorithm (2.3). While the analysis of (4.4) would be more involved. In the rest of this work, we focus on (4.7) and provide a detailed analysis of it. This result follows from a standard ODE analysis. In the present setting, it is guaranteed by the boundedness of U k as shown in Theorem 3.1 and the boundedness of U (t) that will be shown in Theorem 5.1.

Dynamic solution behavior
This section is devoted to the analysis of (4.7), including global existence, asymptotic convergence to equilibria using the celebrated LaSalle invariance principle [14], as well as convergence rates for non-convex objective functions satisfying the PL property.

Global Existence and Uniqueness.
Theorem 5.1. Under Assumption 2.1, the ODE (4.7) with initial condition (θ 0 , r 0 > 0, v 0 ) has a unique global solution (θ(t), v(t), r(t)) ∈ C 1 ([0, ∞]; R 2n+1 ). In particular, for any t > 0, we have Proof. For f ∈ C 1 (R n ), due to Picard and Lindelöf, there exists a unique local solution to (4.7) with initial condition (θ 0 , r 0 , v 0 ). According to the extension theorem, the only way in which global existence can fail is if the solution 'escapes' to infinity. Hence, to establish the global existence, it suffices to show that (θ(t), r(t), v(t)) remains bounded for all t > 0. To do this, we introduce a suitable control function. A candidate based on the discrete counterpart is Using the chain rule, we find that if (θ, r, v) satisfies (4.7), then Hence Q is guaranteed to decrease. In particular, we have which implies the boundedness of F (θ). This together with the coerciveness of f , hence of F , guarantees that θ(t) is bounded for all t > 0. This also ensures the boundedness of v(t). In fact, from (4.7a), using a comparison argument we have On the other hand which is bounded for bounded θ. From (4.7), we see that r decreases monotonically, thus 0 ≤ r(t) ≤ r 0 . This concludes the proof.

Asymptotic behavior of solutions.
Theorem 5.3. Under Assumption 2.1, let (θ(t), r(t), v(t)) be the global solution to the ODE (4.7) as stated in Theorem 5.1 with initial condition (θ 0 , r 0 , v 0 ). If Proof. The upper bound on F follows directly from (5.4). We prove the convergence result in three steps.
Integrating on both sides from 0 to t and using (5.5), we have From this lower bound, we obtain Step 2. We proceed to analyze the asymptotic behavior of θ and v by using LaSalle's invariance principle [14]. Recall the function Q defined in Theorem 5.1, This form and the fact r(t) ≥ r * > 0 implies that is a compact positively invariant set with respect to (4.7), and for any (θ, r, v) ∈ Ω, since r(t) ≥ r * > 0 for all t > 0, thus Next we identify the largest invariant set in Σ. Suppose that at some t 1 , v(t 1 ) = 0, then we haveθ = 0,ṙ = 0, but v(t 1 ) = ∇F (θ(t 1 )). As a result, the trajectory will not stay in Σ unless ∇F (θ(t 1 )) = 0. That is, the largest invariant set in Σ must be By LaSalle's invariance principle, every trajectory of (4.7) starting in Ω tends to Σ 0 as t → ∞. From Step 1, r * is the only limit of r(t). Assume that ∇f (θ) = 0 contains only finite number of critical points of f , then all trajectory admits a limit in since ∇f (θ) = 2F (θ)∇F (θ) and ∇F (θ) = 0.
Step 3. Finally we show that f (θ * ) is a local minimum of f , rather than a local maximum. It suffices to show that F (θ * ) is a local minimum of F . We prove this by contradiction. Suppose F (θ * ) is a local maximum. Recall that F (θ) → F (θ * ), r → r * and |v| → 0 as t → ∞, thus for any δ > 0, there exists t 1 = t 1 (δ) such that for any t ≥ t 1 , Since F is continuous, there exists t 2 ≥ t 1 such that where we have used the assumption that F (θ * ) is a local maximum. From Step 2, we know that Q is non-increasing in t, hence for any t ≥ t 2 , we have provided δ is small so that (r * + δ)δ < 1 2 . Now letting t → ∞,we have F (θ * ) = lim t→∞ Q(t) < F (θ * ). This is a contradiction, hence F (θ * ) has to be a local minimum of F .

Convergence rate
How fast can the solution converge to the equilibrium point depends on the geometry of the objective function near θ * . Indeed, under the PL condition, we can show the flow can reach the equilibrium exponentially fast.
The proof presented below is in three steps: (1) We introduce a candidate control function, with parameters a, λ > 0 to be determined.
(2) We identify admissible pairs (a, λ) such that E(t) decays exponentially fast. and derive the convergence rate of W based on the convergence rate of E. Proof.
Step 1: Decay of the control function. Define where parameters a, λ > 0 will be determined later so that along the trajectory of (4.7) E = E(t) has an exponential decay rate.
For each term in E we proceed to evaluate their time derivative along the trajectory of (4.7). For the first term, we get For the second term, we have Here we used ∇F (θ) = ∇f (θ)/(2F (θ)) and the PL property for f . Hence (6.5) reduces to For the third term, we have d dt Combining the above three estimates together we obtain Note that (6.4) can be rewritten as This upon substitution into (6.6) gives where b(t) := −2ar(t) + 1 + λr(t) F (θ(t)) . This leads to as long as we can identify a, λ so that for all t > 0.
Step 2: Admissible choice for λ, a. Using the solution bounds r * ≤ r(t) ≤ r 0 and In order to ensure (6.10c) we must have b * < 2. Putting these together, condition (6.10) can be ensured by the following stronger constraints: A direct calculation shows that for small , such admissible pair exists, with a = O( ), and λ ≤ F * r 0 which is induced from b * < 2.
The constraint b * < 2 would impose a lower bound for a ∼ (which we should avoid to stay consistency with the discrete algorithm) unless λ is chosen to satisfy b * +2ar * ≤ 2−δ. This is equivalent to requiring With the above two constraints, it is safe to replace b * by 2 − δ in (6.11b) and (6.11c), respectively, so that The convergence rate estimate in the next step requires a to be large as possible, we thus simply take The second relation in (6.12) now reduces to Note that the constraint b * ≥ δ is met if 2ar 0 ≤ 1 − δ + λr * F (θ 0 ) , imposing another upper bound on . To be concrete, we take (6.14) This is (6.2) with r 0 = F (θ 0 ). Hence, for suitably small , we have obtained a set of admissible pairs (a, λ) with as δ varies in (0, 1).

Discussion
This paper studies a new energy-adaptive gradient algorithm with momentum (AGEM). Empirically, we see that the use of momentum can speed up the convergence significantly in a non-convex example. Theoretically, we show that AGEM still enjoys the unconditional energy stability property as AEGD does, and prove convergence to critical points for general non-convex objective functions when the step size is suitably small.
To model the dynamical behavior of AGEM, we derive a high resolution ODE system in which the momentum effect is kept. This system with infinitely many steady states is subtle to analyze. Nonetheless, with several analytical tools, we are able to establish a series of results: (i) the global well-posedness by constructing a suitable Lyapunov function; (ii) the time-asymptotic convergence of the solution trajectory toward a local minimum by using the LaSalle invariance principle; and (iii) a linear convergence rate for non-convex objective functions that satisfy the PL condition.
We believe that the analysis conducted in this article can be extended to study other variants of the AEGD. For instance, the corresponding ODE system for scheme (2.2) of form can be shown to share similar solution properties. For large-scale optimization problems, a stochastic version of AGEM should be the right choice in practice, thus it would be interesting to study the convergence behavior of AGEM in the stochastic setting.