Hedy Attouch
email: hedy.attouch@umontpellier.fr.

Jalal Fadili
email: jalal.fadili@greyc.ensicaen.fr.

FROM THE RAVINE METHOD TO THE NESTEROV METHOD AND VICE VERSA: A DYNAMICAL SYSTEM PERSPECTIVE

Keywords: Ravine method, Nesterov accelerated gradient method, Hessian driven damping, high resolution ODE, convergence rates, Lyapunov analysis, proximal algorithms AMS subject classifications. 37N40, 46N10, 49M30, 65B99, 65K05, 65K10, 90B50, 90C25

We revisit the Ravine method of Gelfand and Tsetlin from a dynamical system perspective, study its convergence properties, and highlight its similarities and differences with the Nesterov accelerated gradient method. The two methods are closely related. They can be deduced from each other by reversing the order of the extrapolation and gradient operations in their definitions. They benefit from similar fast convergence of values and convergence of iterates for general convex objective functions. We will also establish the high resolution ODE of the Ravine and Nesterov methods, and reveal an additional geometric damping term driven by the Hessian for both methods. This will allow us to prove fast convergence towards zero of the gradients not only for the Ravine method but also for the Nesterov method for the first time. We also highlight connections to other algorithms stemming from more subtle discretization schemes, and finally describe a Ravine version of the proximal-gradient algorithms for general structured smooth + non-smooth convex optimization problems.

1. Introduction. In a real Hilbert space H, we revisit the Ravine method of Gelfand and Tsetlin [START_REF] Gelfand | Printszip nelokalnogo poiska v sistemah avtomatich, Optimizatsii[END_REF] from a dynamical system perspective, study its fast convergence properties, and compare it with the Nesterov accelerated gradient method [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF], which we coin here NAG for short. We first consider the case of smooth convex optimization (1.1) min {f (x) : x ∈ H} , where f : H → R is a convex function of class C 1 , whose gradient ∇f is Lipschitz continuous, and which satisfies argmin H (f ) = ∅. We will unveil the close links between the Ravine method and NAG which are sometimes confused in the literature. Indeed, the two methods stem from different discretizations of similar continuous dynamics and can be deduced from each other by reversing the order of the extrapolation and gradient update operations in their definitions. Thus, they benefit from similar fast convergence properties.

On the other hand, the high resolution ODE of the Ravine and Nesterov methods reveal an additional geometric Hessian-driven damping term. The Hessian damping, which is a special case of strong damping in PDE's, plays an important role in attenuating the oscillations. This paves the way to proving new results on fast convergence towards zero of the gradients for both methods. To achieve even better attenuation of the oscillations, we also highlight connections to other algorithms stemming from more subtle discretization schemes of the high resolution ODE. We finally examine the case of "smooth + nonsmooth" structured convex optimization problems, and introduce a first-order inertial proximal gradient algorithm which is based on the Ravine method.

2. Damped inertial dynamical systems for fast optimization. Damped inertial dynamics have a natural mechanical and physical interpretation. Asymptotically, they tend to stabilize the system at a minimizer of the global energy function. As such, they offer an intuitive way to develop fast optimization methods. Let us briefly describe the main damped inertial dynamics used in optimization, their mechanical interpretation, and how the Ravine method stands among them. The Ravine method was a precursor of the accelerated gradient methods. It has long been ignored and, surprisingly enough, is at the forefront of current research. According to the notes in [START_REF] Sra | Optimization for Machine Learning: Subgradient method; Accelerated gradient[END_REF] :

"Ravine method worked well and sparked numerous heuristics for selecting its parameters and improving its behavior. However, its convergence was never proved. It inspired Polyak's heavy-ball method, which seems to have inspired Nesterov's accelerated gradient method ".

2.1. Heavy Ball with Friction (HBF). The heavy ball with friction method was introduced by Polyak in 1964 [START_REF] Polyak | Introduction to optimization[END_REF][START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF]. It describes the movement of a material point of positive mass subjected to a driving force governed by the gradient of the function to be minimized and a viscous friction force. According to the fundamental equation of mechanics, and having normalized the mass equal to one, it is written as follows (HBF) ẍ(t) + γ ẋ(t) + ∇f (x(t)) = 0, where γ > 0 is a fixed viscous damping parameter. The (HBF) method proves to be a useful tool for exploring the local minima of a smooth non-convex function [START_REF] Attouch | The heavy ball with friction method. The continuous dynamical system, global exploration of the local minima of a real-valued function by asymptotical analysis of a dissipative dynamical system[END_REF]. For convex optimization, it is especially interesting in the strongly convex case, where an appropriate choice of the damping coefficient provides linear convergence with an optimal rate. Unfortunately, in the case of a general convex function, it only provides a sublinear rate of convergence of values of order O 1 t . A decisive step to improve it, and to pass from the rate O 1 t to the faster rate O 1 t 2 , has been accomplished by considering algorithms associated with inertial dynamics with asymptotically vanishing damping coefficients. This is the motivation behind the dynamic (AVD α ) and associated algorithm (NAG α ) described hereafter.

2.2. Nesterov Accelerated Gradient method (NAG α ). In recent years, an in-depth study was carried out linking the NAG method to inertial dynamics with vanishing viscosity, see [START_REF] Attouch | Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity[END_REF][START_REF] Attouch | Convergence rates of inertial forward-backward algorithms[END_REF][START_REF] Attouch | Rate of convergence of inertial gradient dynamics with time-dependent viscous damping coefficient[END_REF][START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Cabot | On the long time behavior of second order differential equations with asymptotically small dissipation[END_REF][START_REF] Cabot | Second order differential equations with asymptotically small dissipation and piecewise flat potentials[END_REF][START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF]. Given α a positive parameter, the following second-order ODE (AVD α ) ẍ(t) + α t ẋ(t) + ∇f (x(t)) = 0, was introduced in [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF]. An appropriate temporal discretized version of this ODE with step size s > 0 gives the scheme (NAG α ) which reads

(NAG α ) y k = x k + 1 -α k (x k -x k-1 ) x k+1 = y k -s∇f (y k ).
The scheme (NAG α ) performs a gradient step at y k , which is an extrapolated point obtained from x k and the previous iterate x k-1 .

The method depends in an essential way on the tuning of the extrapolation parameter α k which takes the form α k = 1 -α k in the (NAG α ) scheme. So α k tends to one from below in a subtle controlled way. The historical version of the accelerated gradient method of Nesterov corresponds to α = 3, with the asymptotic convergence rate f (x(t)) -min H f = O 1 t 2 for the continuous dynamic (AVD α ), and

f (x k ) -min H f = O 1
k 2 for the corresponding scheme (NAG α ). Taking α > 3 provides convergence of the trajectories and the improved convergence rate, with small o instead of capital O in the above convergence rates. These results are obtained by Lyapunov analysis [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Attouch | The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than 1/k 2[END_REF][START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF], as summarized below.

Theorem 2.1. Suppose that f : H → R is a convex differentiable function such that ∇f is L-Lipschitz continuous, and S = argmin H (f ) = ∅. Take α ≥ 3, and s ∈]0, 1/L]. Let (x k ) k∈N be a sequence generated by the (NAG α ) algorithm. Set t k = k-1 α-1 , and define, for each integer k ≥ 1

E k := t 2 k (f (x k ) -f (z)) + 1 2s x k-1 -z + t k (x k -x k-1 ) 2 .
Then, the sequence (E k ) k∈N is non-increasing, and as k → +∞

f (x k ) -min H f = O 1 k 2 , x k -x k-1 = O 1 k .
In addition, when α > 3,

f (x k ) -min H f = o 1 k 2 , x k -x k-1 = o 1 k
, and w-lim

x k = x ∈ S,
where w-lim stands for the weak limit.

There has been an active literature devoted to study these questions, from various perspectives, which have given an in-depth understanding of the (NAG α ) method; see for example [START_REF] Apidopoulos | The differential inclusion modeling the FISTA algorithm and optimality of convergence rate in the case b ≤ 3[END_REF][START_REF] Attouch | Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity[END_REF][START_REF] Attouch | Convergence rates of inertial forward-backward algorithms[END_REF][START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Attouch | Rate of convergence of the Nesterov accelerated gradient method in the subcritical case α ≤ 3[END_REF][START_REF] Aujol | Optimal convergence rates for Nesterov acceleration[END_REF][START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF][START_REF] Kim | Optimized first-order methods for smooth convex minimization[END_REF][START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF][START_REF] May | Asymptotic for a second-order evolution equation with convex potential and vanishing damping term[END_REF][START_REF] Muehlebach | A Dynamical Systems Perspective on Nesterov Acceleration[END_REF][START_REF] Siegel | Accelerated first-order methods: Differential equations and Lyapunov functions[END_REF][START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF][START_REF] Villa | Accelerated and inexact forward-backward[END_REF].

Ravine method.

The Ravine method was introduced by Gelfand and Tsetlin [START_REF] Gelfand | Printszip nelokalnogo poiska v sistemah avtomatich, Optimizatsii[END_REF] in 1961. It mimics the flow of water in the mountains which first flows rapidly downhill through small, steep ravines and then flows along the main river in the valley. It also models the transmission of nerve impulses. It has been recently brought to the fore by Polyak [START_REF] Polyak | Accelerated gradient methods revisited[END_REF]. According to the above mechanical interpretation, the Ravine Accelerated Gradient method (coined RAG for short) generates sequences (y k ) k∈N which satisfy

(RAG α ) w k = y k -s∇f (y k ) y k+1 = w k + 1 -α k+1 (w k -w k-1 ) .
Historically, the Ravine method was introduced with a fixed extrapolation coefficient. Taking the extrapolation coefficient equal to 1 -α k+1 makes the Ravine method in accordance with (NAG α ) and is crucial to obtain an accelerated method. A geometric view of (RAG α ) is given in Figure 2.1.

• • • • • y k-1 y k w k-1 = y k-1 -s∇f (y k-1 ) w k = y k -s∇f (y k ) y k+1 = w k + 1 -α k+1 (w k -w k-1 ) S = argmin H (f )
Fig. 2.1: A geometrical illustration of (RAG α ).

2.4. From Ravine to Nesterov and vice versa. The Ravine method has been ignored for a long time, and sometimes confused with (NAG α ) in the literature. Indeed, the two methods can be deduced from each other by reversing the order of the extrapolation and gradient operations. Even more confusing, they come within the same equations. Specifically, the variable y k which enters the definition of (NAG α ) follows the (RAG α ) algorithm. Despite the elementary proof of this result, we state it as a theorem, because of its importance.

Theorem 2.2. (i) Let (x k ) k∈N be a sequence generated by the algorithm (NAG α ). Let (y k ) k∈N be the associated sequence given by

y k = x k + 1 -α k (x k -x k-1
). Then, (y k ) k∈N follows the algorithm (RAG α ). (ii) Conversely, if (y k ) k∈N is a sequence generated by (RAG α ), then the sequence (x k ) k∈N defined by

x k+1 = y k -s∇f (y k ) obeys the algorithm (NAG α ).
Proof.

(i) Suppose that the iterates (x k ) k∈N follow (NAG α ). According to the definition of y k

y k+1 = x k+1 + 1 - α k + 1 (x k+1 -x k ) = y k -s∇f (y k ) + 1 - α k + 1 y k -s∇f (y k ) -(y k-1 -s∇f (y k-1 )) .
Set w k = y k -s∇f (y k ) (which is nothing but x k+1 ). We obtain that the sequence (y k ) k∈N obeys (RAG α ).

(ii) Conversely, from the definition of y k and w k in (RAG α ), we have

y k+1 = y k -s∇f (y k ) + 1 - α k + 1 y k -s∇f (y k ) -(y k-1 -s∇f (y k-1 )) .
Setting x k+1 := y k -s∇f (y k ) as devised, we deduce that

y k+1 = x k+1 + 1 - α k + 1 (x k+1 -x k ) .
Equivalently

y k = x k + 1 - α k (x k -x k-1 )
.

Putting together the above relation with the definition of x k+1 , we obtain that (x k ) k∈N follows (NAG α ). This completes the proof.

Remark 2.1. The order of the two operations, gradient and extrapolation is reversed in the two algorithms. In (NAG α ) first the extrapolation operation is performed, followed by a gradient step. In the Ravine method (RAG α ), it is the reverse order. Although different, this is reminiscent of the approach followed in [START_REF] Attouch | Backward-forward algorithms for structured monotone inclusions in Hilbert spaces[END_REF] which makes it possible to switch from forward-backward algorithms to backward-forward algorithms.

Equipped with Theorem 2.2, we now transfer convergence properties known for (NAG α ) to (RAG α ). In particular, we will show that (NAG α ) and (RAG α ) share the same asymptotic convergence rates.

Theorem 2.3. Let f : H → R be a C 1 convex function whose gradient is L-Lipschitz continuous, and S = argmin H (f ) = ∅. Let (y k ) k∈N be the sequence generated by (RAG α ) with α ≥ 3 and sL ≤ 1. Then, the following properties hold:

(i) f (y k ) -min H f = O 1 k 2 as k → +∞. If, in addition, α > 3, then (ii) f (y k ) -min H f = o 1 k 2 as k → +∞. Let (x k
) k∈N be the associated trajectory generated by (NAG α ), i.e. x k+1 = y k -s∇f (y k ). Then (iii) w-lim y k = w-lim x k ∈ S.

Proof.

(i) According to Theorem 2.2, the sequence (x k ) k∈N defined by (2.1)

x k+1 = y k -s∇f (y k )
follows (NAG α ). Let us take advantage of the convergence properties of (NAG α ), as described in Theorem 2.1. We thus have

(2.2) f (x k ) -min H f = O 1 k 2 , x k -x k-1 = O 1 k .
According to (2.1), we have -1 s (x k+1 -y k ) = ∇f (y k ). Using successively the convex subdifferential inequality, the Cauchy-Schwarz inequality, and the triangle inequality, we obtain

f (y k ) -min H f ≤ f (x k ) -min H f + 1 s x k+1 -y k , x k -y k ≤ f (x k ) -min H f + 1 s x k+1 -y k y k -x k , ≤ f (x k ) -min H f + 1 s ( x k+1 -x k + y k -x k ) y k -x k . (2.3)
Using Theorem 2.2 on the equivalence between (RAG α ) and (NAG α ), we have

y k = x k + 1 - α k (x k -x k-1 ) . Therefore y k -x k ≤ x k -x k-1 . (2.4) Combining (2.3) with (2.4) we obtain (2.5) f (y k ) -min H f ≤ f (x k ) -min H f + 1 s ( x k+1 -x k + x k -x k-1 ) x k -x k-1 .
According to (2.2) we conclude. (ii) Similar arguments when α > 3 yield the o 1

k 2 rate. (iii) Since y k -x k ≤ x k -x k-1 , and x k -x k-1 = O 1
k → 0, we have y k -x k → 0, i.e. y k -x k converges strongly to zero. Combining this with the fact that w-lim x k = x ∈ S, see Theorem 2.1, it follows that the sequence (y k ) k∈N converges weakly to the same limit x .

2.5. Inertial dynamics with Hessian driven damping. It the light of recent work in the study of the acceleration of first order algorithms through the lens of dynamical systems, we will show that the high resolution ODE's of (RAG α ) and (NAG α ) contain an additional Hessian-driven geometric damping term. The underlying dynamic, called (DIN-AVD α,β,b ) is the subject of this section.

The Hessian driven damping plays an important role in various domains. As such, it can be introduced from various perspectives: geometric (also called strong) damping of oscillating systems, which is a central theme in PDE's and mechanics, high resolution ODE of the Ravine method (this will be analyzed in section 4), and regularization of the Newton method. In addition to the viscous damping already present in (AVD α ), taking into account the geometric damping which is driven by the Hessian of the function to be minimized makes it possible to improve the performance of these methods by attenuating their oscillations (see Figure 2.2). In the rest of this section, we follow the lines of [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF], see also [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF] for a special case and motivation. When f is twice continuously differentiable, the dynamic is written as follows

(DIN-AVD α,β,b ) ẍ(t) + α t ẋ(t) + β∇ 2 f (x(t)) ẋ(t) + b(t)∇f (x(t)) = 0,
where b(t) takes into account the temporal scaling effect. The prefix DIN, which stands shortly for Dynamical Inertial Newton, refers to the interpretation of this dynamic as a regularized continuous Newton method, see [START_REF] Attouch | A continuous dynamical Newton-Like approach to solving monotone inclusions[END_REF][START_REF] Attouch | Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators[END_REF][START_REF] Attouch | Continuous Newton-like Inertial Dynamics for Monotone Inclusions, Set Valued and Variational Analysis[END_REF]. At first glance, the presence of the Hessian may seem to entail numerical difficulties. Fortunately, this is not the case as the Hessian intervenes in the form ∇ 2 f (x(t)) ẋ(t), which is nothing but the derivative with respect to time of the function t → ∇f (x(t)).

The temporal discretization of the dynamic (DIN-AVD α,β,b ) with b(t) = 1 + β/t, provides the first-order algorithm proposed in [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF] (IGAHD)

y k = x k + 1 -α k (x k -x k-1 ) -β √ s (∇f (x k ) -∇f (x k-1 )) -β √ s k ∇f (x k-1 ) x k+1 = y k -s∇f (y k ),
where (IGAHD) stands for Inertial Gradient Algorithm with Hessian driven Damping. By comparison with (NAG α ), (IGAHD) has a correction term which contains the difference of the gradients at two consecutive steps. While preserving the convergence properties of (NAG α ), (IGAHD) provides fast convergence to zero of the gradients, and reduces the oscillatory aspects. This is made precise in the following theorem, see [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF][START_REF] Attouch | Convergence of iterates for first-order optimization algorithms with inertia and Hessian driven damping[END_REF].

Theorem 2.4. Let f : H → R be a C 1 convex function whose gradient is L-Lipschitz continuous, and S = argmin H (f ) = ∅. Suppose that α ≥ 3, 0 < β < 2 √ s, sL ≤ 1. Let (x k
) k∈N be a sequence generated by (IGAHD). Then, the following holds:

(i) f (x k ) -min H f = O 1 k 2 , and x k -x k-1 = O 1 k as k → +∞; (ii) k k 2 ∇f (y k ) 2 < +∞ and k k 2 ∇f (x k ) 2 < +∞.
In addition, when α > 3, (iii

) f (x k ) -min H f = o 1 k 2 and x k -x k-1 = o 1 k ; (iv) w-lim x k ∈ S.
A number of other recent papers have contributed to this subject or closely related ones, see [START_REF] Adly | Finite convergence of proximal-gradient inertial algorithms combining dry friction with Hessiandriven damping[END_REF][START_REF] Adly | Finite time stabilization of continuous inertial dynamics combining dry friction with Hessiandriven damping[END_REF][START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian damping[END_REF][START_REF] Kim | Accelerated Proximal Point Method for Maximally Monotone Operators[END_REF][START_REF] Lin | A Control-Theoretic Perspective on Optimal High-Order Optimization[END_REF].

Let us mention another important advantage of (DIN-AVD α,β,b ), which confirms its natural connection with first-order methods. When β > 0, the presence of the Hessian driven damping in the dynamics (DIN-AVD α,β,b ) allows to formulate (DIN-AVD α,β,b ) as an equivalent first-order system both in time and in space, without explicit evaluation of the Hessian. This makes it possible to extend the existence of trajectories and the convergence results to the non-smooth case f ∈ Γ 0 (H) (the class of proper, lower semicontinuous and convex functions on H), by simply replacing the gradient of f by the subdifferential ∂f . This approach was initiated in [START_REF] Álvarez | A second-order gradient-like dissipative dynamical system with Hessiandriven damping. Application to optimization and mechanics[END_REF] and [START_REF] Attouch | A dynamical approach to an inertial forward-backward algorithm for convex minimization[END_REF][START_REF] Attouch | Fast convex minimization via inertial dynamics with Hessian driven damping[END_REF], and used in the perturbed case in [START_REF] Attouch | On the effect of perturbations, errors in first-order optimization methods with inertia and Hessian driven damping[END_REF]. From a mechanical perspective, non-smooth f permits to model non-elastic shocks in unilateral mechanics, see [START_REF] Attouch | A second-order differential system with Hessian-driven damping; Application to non-elastic shock laws[END_REF].

In (DIN-AVD α,β,b ), the Hessian appears explicitly. A closely related ODE is obtained by considering an approach where the Hessian driven damping appears in an implicit form. This was initiated in [START_REF] Alecsa | An extension of the second order dynamical system that models Nesterov's convex gradient method[END_REF], see also [START_REF] Muehlebach | A Dynamical Systems Perspective on Nesterov Acceleration[END_REF] for a related autonomous system in the case of a strongly convex function f . This ODE, coined (ISIHD) for Inertial System with Implicit Hessian Damping, takes the form

(ISIHD) ẍ(t) + α t ẋ(t) + ∇f x(t) + β(t) ẋ(t) = 0,
where α ≥ 3 and

β(t) = γ + β t , γ, β ≥ 0.
The rationale justifying the use of the term "implicit" comes from the observation that by a Taylor expansion (as t → +∞ we have ẋ(t) → 0 which justifies using Taylor expansion), one has

∇f x(t) + β(t) ẋ(t) ≈ ∇f (x(t)) + β(t)∇ 2 f (x(t)) ẋ(t),
hence making the Hessian damping appear indirectly in (ISIHD). As for (DIN-AVD α,β,b ), this ODE was found to have a smoothing effect on the oscillations.

3. The dynamical system perspective of (NAG α ). This sections reviews the close ties between the (NAG α ) algorithm and the associated (AVD α ) system. They will serve as a basis for exploring similar questions for (RAG α ). In doing so, we highlight general methods and tools that allow moving from continuous dynamics to algorithms via temporal discretization, and vice versa.

From continuous dynamics to algorithms and vice versa.

A general and successful recipe to pass from continuous gradient dynamics to gradient algorithms is to follow the following two-step procedure: i) First consider the implicit discretization of the continuous dynamic, and so obtain a proximal algorithm. It is a well known fact that the implicit discretization usually preserves the asymptotic convergence properties of the continuous dynamic. This fact has been well documented for first-order evolution dynamics associated with convex optimization problems [START_REF] Peypouquet | Evolution Equations for Maximal Monotone Operators: Asymptotic Analysis in Continuous and Discrete Time[END_REF], and explains the importance of the proximal algorithm. This type of property is also directly linked with the exponential formula and the Trotter-Lie-Kato formula for the generation of contraction semigroups generated by maximally monotone operators. However, the proximal operator prox sf := (I + s∇f ) -1 , s > 0, may not be easy to compute, which justifies the next step. ii) In the so obtained proximal algorithm, replace the proximal step associated with the operator prox sf by a gradient step associated with the operator I -s∇f . By taking s sufficiently small (typically less than or equal to the inverse of the Lipschitz constant of the gradient of the function which is to be minimized), one can expect to preserve the convergence properties. A major advantage of this procedure is that the proximal and gradient steps have a similar structure and are therefore likely to be combined in proximal gradient algorithms for structured optimization.

3.2. Passing from (AVD α ) to (NAG α ) by temporal discretization. Let us illustrate the above procedure in the case of the (AVD α ) dynamic. Implicit time discretization of (AVD α ), with step size h > 0, gives

x k+1 -2x k + x k-1 h 2 + α kh x k -x k-1 h + ∇f (x k+1 ) = 0.
After multiplication by s = h 2 , we obtain

(3.1) (x k+1 -x k ) -(x k -x k-1 ) + α k (x k -x k-1 ) + s∇f (x k+1 ) = 0. Equivalently (3.2) x k+1 + s∇f (x k+1 ) = x k + 1 - α k (x k -x k-1 ),
which gives

(3.3) x k+1 = prox sf x k + 1 - α k (x k -x k-1 ) .
So, we obtain the inertial proximal algorithm

y k = x k + 1 - α k (x k -x k-1 ) x k+1 = prox sf (y k ) .
This algorithm was initiated by Güler in [START_REF] Güler | New proximal point algorithms for convex minimization[END_REF][START_REF] Güler | On the convergence of the proximal point algorithm for convex optimization[END_REF]. It is a key ingredient of the FISTA method [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Replacing the proximal step by a gradient step, we obtain the (NAG α ) method.

Remark 3.1. One may wonder if a full implicit discretization, which also involves the damping term, leads to the same algorithm. So consider

x k+1 -2x k + x k-1 h 2 + α kh x k+1 -x k h + ∇f (x k+1 ) = 0.
Similar calculation as above gives the inertial proximal algorithm

   y k = x k + 1 1+ α k (x k -x k-1 ) x k+1 = prox s 1+ α k f (y k ) .
The gradient version of the above algorithm takes the form

   y k = x k + k k + α (x k -x k-1 ) x k+1 = y k -sk k+α ∇f (y k ) .
The above form of the extrapolation coefficient, namely k k+α is often used by practitioners, though they maintain the step size equal to s rather than sk k+α , which is obviously asymptotically equal to s. This leads to results similar to those with the extrapolation coefficient 1-α k which is clearly asymptotically equivalent. This is well documented in [START_REF] Attouch | Convergence rates of inertial forward-backward algorithms[END_REF], where the case of general damping and extrapolation coefficients is considered.

3.3.

Passing from (NAG α ) to (AVD α ). We use here a standard limiting argument to pass from (NAG α ) to (AVD α ) (see also [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF]). First write (NAG α ) equivalently as

x k+1 = x k + 1 - α k (x k -x k-1 ) -s∇f x k + 1 - α k (x k -x k-1 ) .
Thus, with s = h 2 , this is also equivalent to

(3.4) x k+1 -2x k + x k-1 h 2 + α kh x k -x k-1 h + ∇f (y k ) = 0.
For each k ∈ N, set t k = kh, and we use the ansatz x k = X(t k ) for some smooth curve t → X(t) defined for t ≥ t 0 > 0. Performing a Taylor expansion in powers of h, when h is close to zero, of the different quantities involved in (NAG α ), we obtain

x k+1 = X(t k+1 ) = X(t k ) + h Ẋ(t k ) + 1 2 h 2 Ẍ(t k ) + O(h 3 ) (3.5) x k-1 = X(t k-1 ) = X(t k ) -h Ẋ(t k ) + 1 2 h 2 Ẍ(t k ) + O(h 3 ). (3.6)
By adding (3.5) and (3.6), we obtain

x k+1 -2x k + x k-1 h 2 = Ẍ(t k ) + O(h).
Moreover, (3.6) gives

x k -x k-1 h = Ẋ(t k ) + O(h).
According to the L-Lipschitz continuity property of ∇f , the definition of y k and 1 -α/k ≤ 1, we have

∇f (y k ) -∇f (x k ) ≤ L y k -x k ≤ L x k -x k-1 . Therefore (3.7) ∇f (y k ) = ∇f (X(t k )) + O(h).
Plugging (3.5), (3.5) and (3.7) into (3.4), we obtain

(3.8) Ẍ(t k ) + α t k Ẋ(t k ) + ∇f (X(t k )) + O(h) = 0.
When h is small, we can neglect the O(h) term which, at the limit, gives that X(•) follows the ODE (AVD α ).

3.4. High resolution ODE of (NAG α ). The high resolution method is extensively used in fluid mechanics, where physical phenomena occur at multiple scales, see for example [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF] for a comprehensive presentation of geophysical fluid dynamics. The idea in our context is not to let h → 0, but to take into account the terms of order h = √ s in the asymptotic expansions, and to discard the terms of order h 2 = s and higher. Moreover, to make the Hessian appear explicitly (see also the discussion on the system (ISIHD) above), we will have to refine the Taylor expansion (3.7). By doing so for (NAG α ), we now show that a Hessian-driven damping term appears in the associated continuous inertial ODE. This is a distinctly new feature and we are not aware of any such a result for (NAG α ).

Theorem 3.1. Assume that f is C 2 . The high resolution ODE with temporal step size √ s of (NAG α ) gives the inertial dynamic with Hessian driven damping

(3.9) Ẍ(t) + √ s∇ 2 f (X(t)) Ẋ(t) + α t Ẋ(t) + 1 + α √ s 2t ∇f (X(t)) = 0.
Proof. Recall the equivalent formulations of (NAG α ) in (3.4) with s = h 2 . For each k ∈ N, set t k := h(k + c) for a real parameter c to be adjusted later, and use the ansatz that x k = X(t k ) for some smooth curve t → X(t) defined for t ≥ t 0 > 0. Performing a Taylor expansion in powers of h, when h is close to zero, of the different quantities involved in (NAG α ), we obtain

x k+1 = X(t k+1 ) = X(t k ) + h Ẋ(t k ) + 1 2 h 2 Ẍ(t k ) + + 1 6 h 3 ... X(tk) + O(h 4 ) (3.10) x k-1 = X(t k-1 ) = X(t k ) -h Ẋ(t k ) + 1 2 h 2 Ẍ(t k ) - 1 6 h 3 ... X(tk) + O(h 4 ). (3.11)
By adding (3.10) and (3.11), we obtain

x k+1 -2x k + x k-1 h 2 = Ẍ(t k ) + O(h 2 ).
Moreover, (3.11) gives

x k -x k-1 h = Ẋ(t k ) - 1 2 h Ẍ(t k ) + O(h 2 ).
We also have

∇f (y k ) = ∇f x k + 1 - α k (x k -x k-1 ) = ∇f x k + h 1 - α k x k -x k-1 h = ∇f X(t k ) + h 1 - α k Ẋ(t k ) - 1 2 h Ẍ(t k ) + O(h 2 ) = ∇f X(t k ) + h 1 - α k Ẋ(t k ) + O(h 2 ) = ∇f (X(t k )) + h 1 - α k ∇ 2 f (X(t k )) Ẋ(t k ) + O(h 2 ).
Putting this with (3.10) and (3.11) into (3.4), we obtain

(3.12) Ẍ(t k ) + α kh Ẋ(t k ) - 1 2 h Ẍ(t k ) + ∇f (X(t k )) + h 1 - α k ∇ 2 f (X(t k )) Ẋ(t k ) + O(h 2 ) = 0.
Equivalently,

(3.13) 1 - α 2k Ẍ(t k ) + α kh Ẋ(t k ) + ∇f (X(t k )) + h 1 - α k ∇ 2 f (X(t k )) Ẋ(t k ) + O(h 2 ) = 0.
Dividing by 1 -α 2k gives

Ẍ(t k ) + α h(k -α 2 ) Ẋ(t k ) + 1 + αh 2h(k -α 2 ) ∇f (X(t k )) + h 1 - α 2 k -α 2 ∇ 2 f (X(t k )) Ẋ(t k ) + O(h 2 ) = 0.
Set c = -α 2 and thus t k := h(k -α 2 ). We obtain

Ẍ(t k ) + α t k Ẋ(t k ) + 1 + αh 2t k ∇f (X(t k )) + h 1 - αh 2 t k ∇ 2 f (X(t k )) Ẋ(t k ) + O(h 2 ) = 0.
By neglecting the term of order s = h 2 , and keeping the terms of order h = √ s, we obtain the claimed inertial dynamic with Hessian driven damping. This completes the proof.

4. The dynamical system perspective of the Ravine method. Let us examine successively the low then the high resolution ODE of the Ravine method.

4.1. Low resolution ODE of (RAG α ). According to the definition of the Ravine method, we have

y k+1 = y k -s∇f (y k ) + 1 - α k + 1 y k -s∇f (y k ) -(y k-1 -s∇f (y k-1 )) = y k + 1 - α k + 1 (y k -y k-1 ) -s∇f (y k ) -s 1 - α k + 1 ∇f (y k ) -∇f (y k-1 ) .
After division by s = h 2 , we obtain, equivalently (4.1)

y k+1 -2y k + y k-1 h 2 + α kh + h y k -y k-1 h + ∇f (y k ) + 1 - α k + 1 (∇f (y k ) -∇f (y k-1 )) = 0.
We now follow an argument similar to the one in section 3.3. For each k ∈ N, set t k = kh, and assume that y k = Y (t k ) for some smooth curve t → Y (t) defined for t ≥ t 0 > 0. Performing a Taylor expansion in powers of h, when h is close to zero, of the different quantities involved in (4.1), we obtain

Ÿ (t k ) + α t k Ẏ (t k ) + ∇f (Y (t k )) + O(h) = 0.
Letting h → 0 gives that Y (•) is a solution trajectory of (AVD α ). We therefore obtain the same inertial dynamics as that associated with (NAG α ).

4.2.

High resolution ODE of (RAG α ). By letting h → 0 in (4.1), the term 1 -α k+1 (∇f (y k ) -∇f (y k-1 )) disappears at the limit. Indeed, as we will see, this term is numerically important. To take it into account, we will perform a high resolution of (4.1). The approach will be similar to that developed in [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF], which will make appear the Hessian-driven damping in the associated continuous inertial equation. This is made precise in the following theorem.

Theorem 4.1. Assume that f is C 2 . The high resolution ODE with temporal step size √ s of (RAG α ) gives the inertial dynamic with Hessian driven damping

(4.2) Ÿ (t) + α t Ẏ (t) + √ s∇ 2 f (Y (t)) Ẏ (t) + 1 + α √ s 2t ∇f (Y (t)) = 0.
Proof. Let us start from the equivalent formulation (4.1) of (RAG α ), which we rewrite as follows

(4.3) y k+1 -2y k + y k-1 h 2 + α k + 1 y k -y k-1 h 2 + ∇f (y k ) + k + 1 -α k + 1 (∇f (y k ) -∇f (y k-1 )) = 0.
Let us arrange the above formula, so as to prepare it for its analysis by Taylor expansion. After multiplying (4.3) by k+1 k+1-α , we get

(4.4) k + 1 k + 1 -α y k+1 -2y k + y k-1 h 2 + α k + 1 -α y k -y k-1 h 2 + k + 1 k + 1 -α ∇f (y k ) + ∇f (y k ) -∇f (y k-1 ) = 0.
Notice then that

y k -y k-1 h 2 = y k+1 -y k h 2 - y k+1 -2y k + y k-1 h 2 .
Thus, (4.4) can be formulated equivalently as follows

k + 1 k + 1 -α - α k + 1 -α y k+1 -2y k + y k-1 h 2 + α k + 1 -α y k+1 -y k h 2 + k + 1 k + 1 -α ∇f (y k ) + ∇f (y k ) -∇f (y k-1 ) = 0.
After reduction we obtain, equivalently

(4.5) y k+1 -2y k + y k-1 h 2 + α (k + 1 -α)h y k+1 -y k h + 1 + α k + 1 -α ∇f (y k ) + ∇f (y k ) -∇f (y k-1 ) = 0.
Building on (4.5), we are now following a device similar to the one developed in section 3.3, and which uses Taylor expansions, but now taken at a higher order. For each k ∈ N, set t k := (k + c)h, where c is a real parameter that will be adjusted later. Assume that y k = Y (t k ) for some smooth curve t → Y (t) defined for t ≥ t 0 > 0. Performing a Taylor expansion in powers of h, when h is close to zero, of the different quantities involved in (4.5), we obtain

y k+1 = Y (t k+1 ) = Y (t k ) + h Ẏ (t k ) + 1 2 h 2 Ÿ (t k ) + 1 6 h 3 ... Y (t k ) + O(h 4 ) (4.6) y k-1 = Y (t k-1 ) = Y (t k ) -h Ẏ (t k ) + 1 2 h 2 Ÿ (t k ) - 1 6 h 3 ... Y (t k ) + O(h 4 ). (4.7)
By adding (4.6) and (4.7) we obtain

y k+1 -2y k + y k-1 h 2 = Ÿ (t k ) + O(h 2 ).
Moreover, (4.6) gives

y k+1 -y k h = Ẏ (t k ) + 1 2 h Ÿ (t k ) + O(h 2 ).
By Taylor expansion of ∇f we have

∇f (y k ) -∇f (y k-1 ) = ∇ 2 f (Y (t k )) Ẏ (t k )h + O h 2 .
Plugging all of the above results into (4.5), we obtain

Ÿ (t k ) + O(h 2 ) + α (k + 1 -α)h Ẏ (t k ) + 1 2 h Ÿ (t k ) + O(h 2 ) + k + 1 k + 1 -α ∇f (Y (t k )) + h∇ 2 f (Y (t k )) Ẏ (t k ) + O h 2 = 0.
After multiplication by (k+1-α)h α , and reduction of the terms involving Ÿ (t k ), we obtain

h α k + 1 - α 2 Ÿ (t k ) + Ẏ (t k ) + (k + 1)h α ∇f (Y (t k )) + h (k + 1 -α)h α ∇ 2 f (Y (t k )) Ẏ (t k ) + O(h 3 ) = 0. Dividing by h α (k + 1 -α 2 ) yields Ÿ (t k ) + α (k + 1 -α 2 )h Ẏ (t k ) + 1 + α 2 k + 1 -α 2 ∇f (Y (t k )) + h 1 - α 2 k + 1 -α 2 ∇ 2 f (Y (t k )) Ẏ (t k ) + O(h 2 ) = 0.
Take c = 1 -α 2 and thus t k := (k + 1 -α 2 )h. We obtain

Ÿ (t k ) + α t k Ẏ (t k ) + 1 + αh 2t k ∇f (Y (t k )) + h 1 - αh 2t k ∇ 2 f (Y (t k )) Ẏ (t k ) + O(h 2 ) = 0.
By neglecting the term of order s = h 2 , and keeping the terms of order h = √ s, we obtain the claimed inertial dynamic with Hessian driven damping. This completes the proof.

A few remarks are in order.

Remark 4.1. The high resolution ODE's of (RAG α ) and (NAG α ) have the same structure but have also differences. First, they are given in terms of two different variables: x for (NAG α ) and y for (RAG α ). Observe also that to get the high resolution ODE (3.9), the Hessian appears after applying an extra Taylor expansion on the gradient, which is reminiscent of our discussion on the implicit Hessian damping ODE (ISIHD). On the other hand, the Hessian appears from an explicit discretization in the ODE (4.2) associated to (RAG α ).

Remark 4.2. Recall that the dynamic with Hessian driven damping (DIN-AVD α,β,b ) which supports the inertial gradient algorithm (IGAHD) developed in [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF] is given by

(4.8) ẍ(t) + α t ẋ(t) + β∇ 2 f (x(t)) ẋ(t) + 1 + β t ∇f (x(t)) = 0.
It is in accordance with the high resolution ODE's (3.9) and (4.2) of respectively (NAG α ) and (RAG α ), and allows to interpret the Hessian driven damping coefficient β of (4.8) as a temporal step size. This also paves the way to proving fast convergence to zero of the gradients which is the subject of the next section.

5.

Fast convergence to zero of the gradients for (RAG α ) and (NAG α ). Let us come to another central point of our study, which concerns the fast convergence towards zero of the gradients. Closely related results were obtained in [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF], and [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Shi | Acceleration via symplectic discretization of high-resolution differential equations[END_REF] in different contexts and discretizations. In [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF], the algorithm considered is (IGAHD), whose underlying dynamic is (DIN-AVD α,β,b ), but the discretization is different. It is inspired by the Nesterov scheme, which contrasts with the Ravine method which is based on an explicit discretization scheme. In [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF], the structure of the algorithm is the same as that of (RAG α ), but the extrapolation parameter is different, which requires an independent proof. Other discretizations are also discussed in [START_REF] Shi | Acceleration via symplectic discretization of high-resolution differential equations[END_REF] after a first-order equivalent reformulation of (4.8). The fast convergence rates on the gradients shown in [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Shi | Acceleration via symplectic discretization of high-resolution differential equations[END_REF] turn out to be weaker than ours. Observe also that developing the Ravine method with a general extrapolation parameter, as was done by Attouch and Cabot in [START_REF] Attouch | Convergence rates of inertial forward-backward algorithms[END_REF] for the accelerated gradient method, is an interesting research venue that we leave to a future work.

5.1. The case of (RAG α ). For the convenience of the reader, we give a self-contained proof, which is based on Lyapunov analysis. We will rely on the following equivalent formulation of the Ravine method which was obtained in (4.1), and which gives rise to the dynamic interpretation with the damping driven by the Hessian:

(5.1)

y k+1 -2y k + y k-1 h 2 + α (k + 1 -α)h y k+1 -y k h + 1 + α k + 1 -α ∇f (y k ) + ∇f (y k ) -∇f (y k-1 ) = 0.
To make the notations shorter, it is convenient to introduce the discrete velocity v k which is defined for each k ∈ N by

v k = 1 h (y k+1 -y k ).
So, the constitutive equation (5.1) can be equivalently written as

(5.2) v k -v k-1 = - α (k + 1 -α) v k -h ∇f (y k ) -∇f (y k-1 -h k + 1 k + 1 -α ∇f (y k ).
Given x ∈ argmin H (f ), our Lyapunov analysis is based on the energy sequence (E k ) k∈N defined by

E k := h 2 (k + 2 -α)(k + 1)(f (y k ) -f (x )) + 1 2 z k 2 (5.
3)

z k := (α -1)(y k+1 -x ) + h(k + 2 -α) v k + h∇f (y k ) . (5.4) Theorem 5.1. Let f : H → R be a C 1 convex
function whose gradient is L-Lipschitz continuous. Let (y k ) k∈N be the sequence generated by (RAG α ), where α ≥ 3 and sL < 1. Then the sequence (E k ) k∈N defined by (5.3)-(5.4) is non-increasing, and the following convergence rate is satisfied:

k k 2 ∇f (y k ) 2 < +∞.
In addition, when α > 3,

k k(f (y k ) -f (x )) < +∞.
Proof. By definition of E k , we have

E k+1 -E k = h 2 (k + 2 -α)(k + 1)(f (y k+1 ) -f (y k )) + h 2 (2k + 4 -α)(f (y k+1 ) -f (x )) + 1 2 z k+1 2 - 1 2 z k 2 . (5.5)
Let us compute this last expression with the help of the elementary identity

(5.6) 1 2 z k+1 2 - 1 2 z k 2 = z k+1 -z k , z k+1 - 1 2 z k+1 -z k 2 .
First observe that the constitutive equation ( 5.2) gives

(5.7) (v k + h∇f (y k )) -(v k-1 + h∇f (y k-1 )) = - α (k + 1 -α) v k -h k + 1 k + 1 -α ∇f (y k ).
Therefore,

(k + 1 -α)(v k + h∇f (y k )) -(k + 1 -α)(v k-1 + h∇f (y k-1 ) = -αv k -h(k + 1)∇f (y k ).
Equivalently,

(5.8) (k + 1)(v k + h∇f (y k )) -(k + 1 -α)(v k-1 + h∇f (y k-1 )) = -h(k + 1 -α)∇f (y k ).
Using successively the definition of z k and (5.8), we obtain

z k+1 -z k = (α -1)(y k+2 -y k+1 ) +h(k + 3 -α) v k+1 + h∇f (y k+1 ) -h(k + 2 -α) v k + h∇f (y k ) = h(α -1)v k+1 + h(k + 3 -α) v k+1 + h∇f (y k+1 ) -h(k + 2 -α) v k + h∇f (y k ) = h(k + 2) v k+1 + h∇f (y k+1 ) -h(k + 2 -α) v k + h∇f (y k ) -h 2 (α -1)∇f (y k+1 ) = -h 2 (k + 2 -α)∇f (y k+1 ) -h 2 (α -1)∇f (y k+1 ) = -h 2 (k + 1)∇f (y k+1 ).
Plugging this into (5.6), we obtain

1 2 z k+1 2 - 1 2 z k 2 = - 1 2 h 4 (k + 1) 2 ∇f (y k+1 ) 2 -h 2 (k + 1) ∇f (y k+1 ), (α -1)(y k+1 -x ) + h(k + 2 -α) v k + h∇f (y k ) -h 2 (k + 1)∇f (y k+1 ) = 1 2 h 4 (k + 1) 2 ∇f (y k+1 ) 2 -h 2 (k + 1) ∇f (y k+1 ), (α -1)(y k+1 -x ) + h(k + 2 -α) v k + h∇f (y k ) .
Let us rearrange this last expression so that we have terms involving only ∇f (y k+1 ). For this, we use (5.7) that we write as follows

v k + h∇f (y k )) = v k+1 + h∇f (y k+1 ) + α k + 2 -α v k+1 + h k + 2 k + 2 -α ∇f (y k+1 ).
Therefore,

(α -1)(y k+1 -x ) + h(k + 2 -α) v k + h∇f (y k ) = (α -1)(y k+1 -x ) + h(k + 2 -α) v k+1 + h∇f (y k+1 ) + α k + 2 -α v k+1 + h k + 2 k + 2 -α ∇f (y k+1 ) = (α -1)(y k+1 -x ) + h(k + 2)v k+1 + h 2 (2k + 4 -α)∇f (y k+1 ).
Collecting the above results we obtain

1 2 z k+1 2 - 1 2 z k 2 = 1 2 h 4 (k + 1) 2 ∇f (y k+1 ) 2 -h 2 (k + 1) ∇f (y k+1 ), (α -1)(y k+1 -x ) + h(k + 2)v k+1 + h 2 (2k + 4 -α)∇f (y k+1 ) .
Combining this inequality with (5.5) we get (5.9)

E k+1 -E k = h 2 (k + 2 -α)(k + 1)(f (y k+1 ) -f (y k )) + h 2 (2k + 4 -α)(f (y k+1 ) -f (x )) + 1 2 h 4 (k + 1) 2 ∇f (y k+1 ) 2 -h 2 (k + 1) ∇f (y k+1 ), (α -1)(y k+1 -x ) + h(k + 2)v k+1 + h 2 (2k + 4 -α)∇f (y k+1 ) .
According to the basic gradient inequality for convex differentiable functions whose gradient is L-Lipschitz continuous, we have

f (y k ) ≥ f (y k+1 ) + ∇f (y k+1 ), y k -y k+1 + 1 2L ∇f (y k+1 ) -∇f (y k ) 2 , f (x ) ≥ f (y k+1 ) + ∇f (y k+1 ), x -y k+1 .
Combining the above inequalities with (5.9), we obtain

E k+1 -E k ≤ -h 2 (k + 2 -α)(k + 1) ∇f (y k+1 ), y k -y k+1 + 1 2L ∇f (y k+1 ) -∇f (y k ) 2 +h 2 (2k + 4 -α)(f (y k+1 ) -f (x )) + h 2 (k + 1)(α -1)(f (x ) -f (y k+1 )) + 1 2 h 4 (k + 1) 2 ∇f (y k+1 ) 2 -h 2 (k + 1) ∇f (y k+1 ), h(k + 2)v k+1 + h 2 (2k + 4 -α)∇f (y k+1 ) .
Equivalently,

E k+1 -E k ≤ h 2 (k + 2 -α)(k + 1) ∇f (y k+1 ), hv k - 1 2L ∇f (y k+1 ) -∇f (y k ) 2 -h 2 k(α -3) + 2α -5 (f (y k+1 ) -f (x )) -h 2 (k + 1) ∇f (y k+1 ), h(k + 2)v k+1 - 1 2 h 4 (k + 1)(3k + 7 -2α) ∇f (y k+1 ) 2 .
Let us put together the terms involving the scalar product with ∇f (y k+1 ). We get

E k+1 -E k ≤ -h 2 (k + 2 -α)(k + 1) 1 2L ∇f (y k+1 ) -∇f (y k ) 2 +h 3 (k + 1) ∇f (y k+1 ), (k + 2 -α)v k -(k + 2)v k+1 -h 2 k(α -3) + 2α -5 (f (y k+1 ) -f (x )) - 1 2 h 4 (k + 1)(3k + 7 -2α) ∇f (y k+1 ) 2 .
According to (5.8) we have

(5.10) (k + 2 -α)(v k + h∇f (y k )) -(k + 2)(v k+1 + h∇f (y k+1 )) = h(k + 2 -α)∇f (y k+1 ). Therefore (5.11) (k + 2 -α)v k -(k + 2)v k+1 = -h(k + 2 -α)∇f (y k ) + h(2k + 4 -α)∇f (y k+1 ).
Combining the above results we get

E k+1 -E k + h 2 k(α -3) + 2α -5 (f (y k+1 ) -f (x )) ≤ - h 2 2L (k + 2 -α)(k + 1) ∇f (y k+1 ) -∇f (y k ) 2 +h 3 (k + 1) ∇f (y k+1 ), -h(k + 2 -α)∇f (y k ) + h(2k + 4 -α)∇f (y k+1 ) - 1 2 h 4 (k + 1)(3k + 7 -2α) ∇f (y k+1 ) 2 .

Equivalently

(5.12)

E k+1 -E k + h 2 k(α -3) + 2α -5 (f (y k+1 ) -f (x )) + R(∇f (y k ), ∇f (y k+1 )) ≤ 0, where R(X, Y ) = h 2 2L (k + 2 -α)(k + 1) Y -X 2 + 1 2 h 4 (k + 1)(3k + 7 -2α) Y 2 -h 3 (k + 1) Y, -h(k + 2 -α)X + h(2k + 4 -α)Y .
To conclude, we just need to prove that the quadratic form R is positive definite. A simple procedure consists in computing inf X R(X, Y ). For fixed Y , the minimum of

R(•, Y ) is achieved at X with X -Y = -h 2 LY . Therefore inf X R(X, Y ) = h 4 L 2 (k + 2 -α)(k + 1)h 2 Y 2 + 1 2 h 4 (k + 1)(3k + 7 -2α) Y 2 -h 3 (k + 1) Y, -h(k + 2 -α)(1 -h 2 L)Y + h(2k + 4 -α)Y .
After reduction, we get

inf X R(X, Y ) = h 4 (k + 1) 2 Y 2 k(1 -Lh 2 ) + Lh 2 (α -2) + 3 -2α .
Returning to (5.12) we obtain

E k+1 -E k + h 2 k(α -3) + 2α -5 (f (y k+1 ) -f (x )) + h 4 (k + 1) 2 ∇f (y k+1 ) 2 k(1 -Lh 2 ) + Lh 2 (α -2) + 3 -2α ≤ 0.
So, when α ≥ 3 and Lh 2 ∈]0, 1[, we obtain that (E k ) k∈N is a non-negative non-increasing sequence, hence convergent. By summing the above inequalities over k, we finally obtain

k k 2 ∇f (y k+1 ) 2 < +∞.
Also note that when α > 3 we obtain

k k(f (y k+1 ) -f (x )) < +∞.
Remark 5.1. It is easy to see that our result entails

min 1≤i≤k ∇f (y i ) 2 = O 1 k 3 ,
which recovers the rate found in [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Shi | Acceleration via symplectic discretization of high-resolution differential equations[END_REF].

5.2. The case of (NAG α ). We have the following fast convergence of gradients to zero for (NAG α ). Surprisingly, this result has never been established before in the literature.

Theorem 5.2. Let f : H → R be a C 1 convex function whose gradient is L-Lipschitz continuous. Let (x k ) k∈N be the sequence generated by (NAG α ), where α ≥ 3 and sL < 1. Then the following convergence rate is satisfied:

k k 2 ∇f (x k ) 2 < +∞.
In addition, when α > 3,

k k(f (x k ) -f (x )) < +∞.
Proof. The analysis relies on the energy function

E k (t) = t 2 k (f (x k ) -min H f ) + 1 2s x k-1 -x + t k (x k -x k-1 ) 2 ,
where t k := k-1 α-1 . We will also need the following refined version of the descent lemma

(5.13) f (y -s∇f (y)) ≤ f (x) + ∇f (y), y -x - s 2 ∇f (y) 2 - s 2 ∇f (x) -∇f (y) 2
valid for any (x, y) ∈ H 2 and s ∈]0, 1/L]. Let us write successively (5.13) at y = y k and x = x k , then at y = y k and x = x . Recalling that x k+1 = y k -s∇f (y k ), we obtain the two inequalities

f (x k+1 ) ≤ f (x k ) + ∇f (y k ), y k -x k - s 2 ∇f (y k ) 2 - s 2 ∇f (x k ) -∇f (y k ) 2 , (5.14) f (x k+1 ) ≤ min H f + ∇f (y k ), y k -x -s ∇f (y k ) 2 .
(5.15) Multiplying (5.14) by t k+1 -1 ≥ 0, then adding (5.15), we derive that

t k+1 f (x k+1 ) ≤ (t k+1 -1)f (x k ) + min H f + ∇f (y k ), (t k+1 -1)(y k -x k ) + y k -x - s 2 (t k+1 + 1) ∇f (y k ) 2 - s 2 (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 . (5.16) Observing that t k verifies the identity α k := 1 - α k = t k -1 t k+1
, we obtain

(t k+1 -1)(y k -x k ) + y k = t k+1 y k -(t k+1 -1)x k = x k + t k+1 α k (x k -x k-1 ) = x k-1 + (1 + t k+1 α k )(x k -x k-1 ) = x k-1 + t k (x k -x k-1 ). Setting z k = x k-1 + t k (x k -x k-1
), we then deduce from (5.16) that

t k+1 (f (x k+1 ) -min H f ) ≤ (t k+1 -1)(f (x k ) -min H f ) + ∇f (y k ), z k -x - s 2 (t k+1 + 1) ∇f (y k ) 2 - s 2 (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 . (5.17)
On the other hand, observe that

z k+1 -z k = x k + t k+1 (x k+1 -x k ) -x k-1 -t k (x k -x k-1 ) = t k+1 (x k+1 -x k ) -(t k -1)(x k -x k-1 ) = t k+1 (x k+1 -x k -α k (x k -x k-1 )) = t k+1 (x k+1 -y k ) = -s t k+1 ∇f (y k ).
It ensues that (5.18)

z k+1 -x 2 = z k -x 2 -2st k+1 ∇f (y k ), z k -x + s 2 t 2 k+1 ∇f (y k ) 2 .
By using this equality in (5.17), we find

t k+1 (f (x k+1 ) -min H f ) ≤ (t k+1 -1)(f (x k ) -min H f ) + 1 2st k+1 ( z k -x 2 -z k+1 -x 2 ) + st k+1 2 ∇f (y k ) 2 - s 2 (t k+1 + 1) ∇f (y k ) 2 - s 2 (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 . (5.19)
which is equivalent to

t 2 k+1 (f (x k+1 ) -min H f ) + 1 2s z k+1 -x 2 ≤ (t 2 k+1 -t k+1 )(f (x k ) -min H f ) + 1 2s z k -x 2 - st 2 k+1 2 ∇f (y k ) 2 + (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 .
Using the expression of the sequence (E k ) k∈N , we obtain

E k+1 ≤ E k + (t 2 k+1 -t 2 k -t k+1 )(f (x k ) -min H f ) - st 2 k+1 2 ∇f (y k ) 2 + (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 .
Elementary algebraic computation gives

∇f (y k ) 2 + (t k+1 -1) ∇f (x k ) -∇f (y k ) 2 ≥ t k+1 -1 t k+1 ∇f (x k ) 2 .
We finally obtain (5.20)

E k+1 + s 2 t k+1 (t k+1 -1) ∇f (x k ) 2 ≤ E k + (t 2 k+1 -t 2 k -t k+1 )(f (x k ) -min H f ).
By definition of t k , we have t 2 k+1 -t 2 k -t k+1 = -k(α-3)+1 (α-1) 2 ≤ 0 for α ≥ 3. Thus, (E k ) k∈N is a non-increasing non-negative function, and summing (5.20), we obtain leading to a rate on the gradients

k k(k + 1 -α) ∇f (x k ) 2 < +∞.
∇f (x k ) 2 = O log k k 4 .
This is in agreement with our result in Theorem 5.2. On the other hand, one infers from our result that ∇f (x k ) 2 must decrease at least as fast as O 1 k 3 (log k) ν , for ν > 1. 6. The Ravine accelerated proximal gradient method. Let us now extend the Ravine method (RAG α ) to the case of additively structured "smooth + non-smooth" convex minimization problems (6.1) min

x∈H {θ(x) := f (x) + g(x)} ,
where we make the following assumptions (H)

     f : H → R is a C 1 convex
(θ) = ∅.
Note that by the above assumptions, θ : H → R ∪ {+∞} is a proper, convex and lower semicontinuous function. A natural extension of (NAG α ) to this setting is (6.2)

y k = x k + 1 -α k (x k -x k-1 ) x k+1 = prox sg (y k -s∇f (y k ))
which generalizes FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Let us recall the convergence properties of this algorithm, which are valid under the assumption sL ≤ 1:

• α = 3: According to [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], one has

θ(x k ) -min H θ = O 1 k 2 as k → +∞.
• α > 3: According to [START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF] and [START_REF] Attouch | The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than 1/k 2[END_REF] θ(x k ) -min

H θ = o 1 k 2 , x k -x k-1 = o 1 k
as k → +∞, and w-lim x k ∈ S.

• α ≤ 3: According to [START_REF] Apidopoulos | Convergence rate of inertial Forward-Backward algorithm beyond Nesterov's rule[END_REF] and [START_REF] Attouch | Rate of convergence of the Nesterov accelerated gradient method in the subcritical case α ≤ 3[END_REF] θ(x k ) -min

H θ = O k -2α 3 .
Before presenting our algorithm, let us introduce the operator T s : H → H defined by (6.3) T s (y) = 1 s y -prox sg (y -s∇f (y)) .

Note that solving (6.1) is equivalent to find a zero of T s . In addition, the operator T s reduces to the gradient operator ∇f when g = 0. Thus the algorithm (6.2) can be formulated in an equivalent way in the following form

y k = x k + 1 -α k (x k -x k-1 ) x k+1 = y k -sT s (y k ).
As a consequence, all the algebraic developments concerning the Ravine method (RAG α ) can be extended to the structured additive setting, by just replacing the gradient operator ∇f by T s . Indeed, one can easily show that for sL ≤ 1, the two operators share the following properties: monotonicity, co-coercivity and Lipschitz continuity which play a central role in the Lyapunov analysis.

With the help of this analogy, we are now in position to introduce the Ravine Accelerated Proximal Gradient algorithm ((RAPG α ) for short):

(RAPG α ) w k = y k -sT s (y k ) y k+1 = w k + 1 -α k (w k -w k-1 ).
Figure 6.1 gives some geometric insight into the scheme (RAPG α ).

• • • • • y k-1 y k w k-1 = prox sg (y k-1 -s∇f (y k-1 )) w k = prox sg (y k -s∇f (y k )) y k+1 = w k + 1 -α k+1 (w k -w k-1 )
S = argmin H (f + g) Fig. 6.1: Geometrical illustration of the (RAPG α ) algorithm.

By following an argument similar to that of Theorem 4.1, the high resolution ODE of the algorithm (RAPG α ) gives (6.4)

Ÿ (t) + α t Ẏ (t) + √ s d dt T s (Y (t)) + 1 + α √ s 2t T s (Y (t)) = 0,
where the term d dt T s (Y (t)) is interpreted as the distributional derivative of the absolutely continuous function t → T s (Y (t)). The ODE (6.4) is a Regularized Inertial Newton dynamic which has been recently studied in [START_REF] Adly | Asymptotic behavior of Newton-like inertial dynamics involving the sum of potential and nonpotential terms[END_REF][START_REF] Adly | Newton-type inertial algorithms for solving monotone equations governed by sums of potential and nonpotential operators[END_REF] and [START_REF] Attouch | Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators[END_REF][START_REF] Attouch | Continuous Newton-like Inertial Dynamics for Monotone Inclusions, Set Valued and Variational Analysis[END_REF]. The existence and uniqueness of a strong solution to the Cauchy problem associated with (6.4) has been proved in [START_REF] Adly | Asymptotic behavior of Newton-like inertial dynamics involving the sum of potential and nonpotential terms[END_REF]Theorem 2.1]. It is based on the equivalent reformulation of (6.4) as a first-order system.

Let us establish some fast convergence properties of (RAPG α ) which can be deduced from the well established results concerning (6.2). Theorem 6.1. Suppose that (H) holds, α ≥ 3, and sL ∈]0, 1]. Let (y k ) k∈N and (w k ) k∈N be the sequences generated by (RAPG α ). Then the following properties are satisfied:

θ(w k ) -min H θ = O 1 k 2 as k → +∞; When α > 3, θ(w k ) -min H θ = o 1 k 2 and w-lim y k = w-lim w k ∈ S.
Proof. As a key ingredient, we use the following analog of the gradient descent lemma for composite optimization, see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], [START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF]: for any (x, y) ∈ H 2 (6.5)

θ(y -sT s (y)) ≤ θ(x) + T s (y), y -x - s 2 T s (y) 2 .
Taking x = w k-1 y = y k , we deduce that

(6.6) θ(w k ) ≤ θ(w k-1 ) + T s (y k ), y k -w k-1 .
By Cauchy-Schwarz inequality, we get

(6.7) θ(w k ) ≤ θ(w k-1 ) + T s (y k ) y k -w k-1 .
By definition of (RAPG α ), and by using the link with (6.2) (recall that w k = x k+1 where (x k ) k∈N are the iterates generated by (6.2)), we have

T s (y k ) = 1 s w k -y k ≤ 1 s ( w k -w k-1 + y k -w k-1 ) ≤ 1 s ( w k -w k-1 + w k-1 -w k-2 ) = 1 s ( x k+1 -x k + x k -x k-1 ) = O 1 k ,
where the latter rate is known from [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Attouch | The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than 1/k 2[END_REF]. In a similar way, (6.8)

y k -w k-1 ≤ w k-1 -w k-2 = x k -x k-1 = O 1 k .
By plugging the above estimates into (6.7), we obtain

θ(w k ) -min H θ ≤ θ(x k ) -min H θ + C k 2 = O 1 k 2 .
where the rate O(1/k 2 ) on θ(x k ) -min H is known from [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

For α > 3, we know from [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Attouch | The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than 1/k 2[END_REF]

that θ(x k ) -min H = o 1 k 2 and x k -x k-1 = o 1 k .
We thus argue as above to obtain the claimed o 1 k 2 rate. In addition, we have from (6.8) that y k -x k → 0, i.e. y k -x k converges strongly to zero. Since the sequence (x k ) k∈N converges weakly when α > 3, see [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF][START_REF] Attouch | The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than 1/k 2[END_REF], it follows that the sequence (y k ) k∈N converges weakly to the same limit as (x k ) k∈N .

7. The strongly convex case. In this section, we briefly discuss the strongly convex case. Recall that f : H → R is said to be µ-strongly convex for some µ > 0 if f -µ 2 • 2 is convex. In this case, a proper tuning of the viscous damping coefficient in the dynamic (HBF) of Polyak provides exponential convergence rate with optimal rate. Theorem 7.1. Suppose that f : H → R is a C 1 and µ-strongly convex function for some µ > 0. Let x(•) : [t 0 , +∞[→ H be a solution trajectory of (HBF) with γ = 2 √ µ, i.e.

(7.1) ẍ(t) + 2 √ µ ẋ(t) + ∇f (x(t)) = 0, with initial condition (x(t 0 ), ẋ(t 0 )), t 0 ≥ 0. Then, for all t ≥ t 0 f (x(t)) -min

H f ≤ Ce -√ µ(t-t0)
where

C := f (x(t 0 )) -min H f + µ x(t 0 ) -x 2 + ẋ(t 0 ) 2 .
According to the procedure described in section 3.1, we consider three different discretizations of (7.1) inspired by the inertial proximal algorithm, then the Nesterov method, and finally the Ravine method. Let h > 0 be the temporal step size. Inertial proximal algorithm. Implicit time discretization of (7.1) gives

x k+1 -2x k + x k-1 h 2 + √ µ x k+1 -x k-1 h + ∇f (x k+1 ) = 0.
After multiplication by s = h 2 , we obtain

(7.2) (1 + h √ µ)(x k+1 -x k ) + s∇f (x k+1 = (1 -h √ µ)(x k -x k-1 ), which gives (7.3) x k+1 = prox s 1+ √ µs f x k + 1 - √ µs 1 + √ µs (x k -x k-1 ) .
So, we obtain the inertial proximal algorithm (7.4)

y k = x k + 1- √ µs 1+
√ µs (x k -x k-1 ) x k+1 = prox s 1+ √ µs f (y k ) Nesterov method. Replacing the proximal step by a gradient step in (7.4), we obtain

y k = x k + 1- √ µs 1+ √ µs (x k -x k-1 ) x k+1 = y k -s 1+ √ µs ∇f (y k ).
Ravine method. Interverting the role of the variables x k and y k we obtain the Ravine method

     w k = y k - s 1 +
√ µs ∇f (y k )

y k+1 = w k + 1 - √ µs 1 + √ µs (w k -w k-1 ).
This scheme is closely related to the classical form of the Nesterov accelerated gradient method for strongly convex minimization. Note however that our approach makes appear a gradient step size s 1+ √ µs which is slightly different from the most usually used step size s. A detailed Lyapunov analysis for this scheme is an interesting subject that we leave for a future work.

8. Comparison to related algorithms. The following table gives a bird's eye view of the relationships between (NAG α ), (RAG α ). Comparison with (IGAHD). Recall that the algorithm (IGAHD), introduced by the authors in [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF], is based on the dynamic (4.8), a special case of (DIN-AVD α,β,b ), with damping parameters α ≥ 3 and β ≥ 0. This dynamic is essentially the same as the high resolution ODE (4.2) associated to (RAG α ) (the same holds for the ODE resp. (3.9) of (RAG α )). The two dynamics can be deduced from each other by a linear temporal reparameterization, which preserves their convergence properties. However, an important message to keep in mind here is that the algorithms (RAG α ) and (IGAHD) markedly differ in the type of temporal discretization used to obtain them. The algorithm (RAG α ) is obtained by explicit discretization of (4.2), namely 2). In a nutshell, from a theoretical point of view, these methods behave similarly. Nevertheless, it was shown in [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF] that, numerically, (IGAHD) enjoys much less oscillations that (NAG α ) (and thus (RAG α )). We conjecture that this is a consequence of the more subtle discretization underlying (IGAHD). So far, this lacks clear theoretical justification and we believe that it is a nice research program to undertake in the future.

9. Conclusion, Perspective. This work was intended to unveil the relationship between the Nesterov accelerated method and the Ravine method, which as been ignored for a long time and sometimes confused with the Nesterov. We have shed light on these connections through the perspective of dynamical systems. We believe that this work paves the way to many important questions that remain to be answered. Among them, we mention the following ones:

• Design better structure-preserving discretization schemes/algorithms for inertial systems, and understand their fundamental limits/performance.

• For additively structured "smooth + nonsmooth" convex minimization problems, develop a Lyapunov analysis showing the fast convergence to zero of the operators (which correspond to the gradients for the Ravine accelerated gradient method). • Develop a Ravine accelerated method for linearly constrained optimization problems, which is based on the augmented Lagrangian approach, and the ADMM algorithm. • Study the introduction of perturbation, errors into the Ravine method, so as to prepare the stochastic versions of this algorithm. • Generalization and tuning of the extrapolation parameter.

• The case of monotone inclusions.

Fig. 2 . 2 :

 22 Fig. 2.2: Evolution of the objective (left) and trajectories (right) for (AVD α ) (α = 3.1) and (DIN-AVD α,β,b ) (α = 3.1, β = 1) on an ill-conditioned quadratic problem in R 2 .

When α > 3 ,

 3 summing again we get the second claimed estimatek (k(α -3) + 1)(f (x k ) -min H f ) < +∞.Remark 5.2. Another way to show Theorem 5.2 is to use Theorem 5.1 and the equivalence result in Theorem 2.2 between (NAG α ) and (RAG α ).

Remark 5 . 3 . 3 ,

 533 Theorem 5.2 yields the rate min 1≤i≤k ∇f (x i ) 2 = O 1 k which matches the complexity bound in [44, Item 2.]. In [44, Item 3.], a better complexity bound is obtained by applying (NAG α ) with α = 3 to a Tikhonov regularization of f with an asymptotically vanishing parameter. From those complexity bounds, one can straightforwardly show that this parameter has to scale as O log k k 2

(8. 1 )

 1 y k+1 -2y k + y k-1 h 2 + α (k + 1 -α)h y k+1 -y k h + 1 + α k + 1 -α ∇f (y k ) + ∇f (y k ) -∇f (y k-1 ) = 0.On the other hand (IGAHD) is obtained by the time discretization of (4.8)1 s (x k+1 -2x k + x k-1 ) + α ks (x k -x k-1 ) + β √ s (∇f (x k ) -∇f (x k-1 )) + β k √ s ∇f (x k-1 ) + ∇f (y k ) = 0,with y k is an extrapolated point inspired by Nesterov's scheme. The convergence properties of (IGAHD), recalled in Theorem 2.4 are similar to those (RAG α ) and (NAG α ) (see Theorem 2.1, Theorem 2.3, Theorem 5.1 and Theorem 5.

  function and ∇f is L-Lipschitz continuous; g : H → R ∪ {+∞} is proper, convex and lower semicontinuous; S := argmin H

Acknwoledgements. We would like to thank N. Boumal who brought the reference [44] to our attention.