{'fraction_non_alphanumeric': 0.06504065040650407, 'fraction_numerical': 0.01883602994445786, 'mean_word_length': 4.004027926960258, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast convex optimization via time scale and averaging of the steepest descent

In a Hilbert setting, we develop a gradient-based dynamic approach for fast solving convex optimization problems. By applying time scaling, averaging, and perturbation techniques to the continuous steepest descent (SD), we obtain high-resolution ODEs of the Nesterov and Ravine methods. These dynamics involve asymptotically vanishing viscous damping and Hessian driven damping (either in explicit or implicit form). Mathematical analysis does not require developing a Lyapunov analysis for inertial systems. We simply exploit classical convergence results for (SD) and its external perturbation version, then use tools of differential and integral calculus, including Jensen's inequality. By way of illustration of the flexibility of the method we consider the case where the initial dynamics is the regularized Newton method. Then, we show that the technique can be naturally extended to the case of a monotone cocoercive operator. Our approach leads to parallel algorithmic results, which we study in the case of proximal algorithms.


Introduction.
In a real Hilbert space H, we develop a general dynamic approach whose objective is the rapid resolution of convex optimization problems via first-order methods. We consider the minimization problem where, throughout the paper, we make the following assumptions on the function f to be minimized ∇f is Lipschitz continuous on the bounded sets of H.
Our study is part of the close links between dissipative dynamical systems and optimization algorithms, the latter being obtained by temporal discretization of the continuous dynamics. The dynamic system approach makes it possible to harmoniously combine the ingredients of our approach, namely time scaling, averaging, and perturbation techniques to finally obtain inertial dynamics with proven fast optimization properties. Let us recall the mainstream for the study of accelerated gradient methods. Then, in the next section, we will describe our approach, which is novel in many respects. The classical approach can be traced back to Gelfand and Tsetlin [36], then Polyak [41,42], and Attouch, Goudou and Redont [16], where fast first-order optimization algorithms are naturally linked to damped inertial dynamics via temporal discretization. In this approach, which is very inspired by mechanics, PDE's, and control theory the focus is on the design of the damping term of the inertial dynamic, and its asymptotic stabilization effect.
In recent years, for the minimization of general convex differentiable functions, an in-depth study has been carried out linking the accelerated gradient method of Nesterov ( [40,39]) to inertial dynamics with vanishing viscosity (the viscous damping coefficient tends to zero as t → +∞). That is the Su-Boyd-Candès dynamical system [45] (1.3)ẍ(t) + α tẋ (t) + ∇f (x(t)) = 0, where fast optimization is obtained by taking α ≥ 3. In addition to viscous damping, taking into account geometric damping (in our situation it is driven by the Hessian of the function to be minimized) makes it possible to improve the performance of these methods by attenuating the oscillations. That is the dynamic (1.4)ẍ(t) + α tẋ (t) + β∇ 2 f (x(t))ẋ(t) + ∇f (x(t)) = 0, considered by Attouch, Peypouquet, and Redont in [20]. The central role played by these dynamics has been confirmed by their recent interpretation as the low and high resolution of the Nesterov and Ravine accelerated gradient methods, see Shi, Du, Jordan and Su [44], Attouch, Chbani, Fadili, and Riahi [10,9]. The study of the fast convergence properties of these dynamics and algorithms is based on Lyapunov analysis. But this approach does not give a simple explanation of the choice of dynamics and algorithms. Moreover the Lyapunov analysis is not trivial, and there is no general rule to find the Lyapunov functions. Still, these techniques have allowed a deep understanding of dynamics and accelerated algorithms for optimization, see [3], [4], [6,7], [10], [11], [12], [14], [19], [26], [28], [32], [33], [34], [45], [44] for some recent contributions.
2. Time scaling and averaging approach. We propose a new approach to these questions which is based on a combination of time scaling, averaging, and perturbation techniques. The basis property which underlies our approach is the fact that by time scaling one can always accelerated a dynamic. The kinematic interpretation is clear. By changing the time clock the trajectories are travelled more or less quickly. The averaging technique shifts from a dynamic of order n to a dynamic of the immediately higher order n + 1. While preserving the asymptotic convergence properties, this makes it possible to achieve a dynamic where the coefficient in front of the gradient is fixed, and therefore suitable for developing an accelerated gradient algorithm. The perturbation technique gives flexibility to these methods, exploiting the fact that the convergence properties are preserved when adding an external perturbation which vanishes asymptotically sufficiently fast.
Under the standing assumption (A) on f we know that, for any z 0 ∈ H, there exists a unique classical global solution z ∈ C 1 ([t 0 , +∞[: H) of (SD) satisfying z(t 0 ) = z 0 , see [5,Theorem 17.1.1]. We fix t 0 > 0 as the origin of time (because of the singularity at the origin of coefficient α t which enter into our analysis).
2.1. Classical facts concerning the continuous steepest descent. We will have to consider the gradient system with an external perturbation. The following theorem investigates its asymptotic properties under appropriate conditions on the perturbation term.
Then the following statements are true: i) (convergence of gradients towards zero) ∇f (z (t)) = o 1 √ t as t → +∞; ii) (integral estimate of the gradients) as t → +∞; iv) the solution trajectory z(t) converges weakly as t → +∞, and its limit belongs to S = argmin f .
Proof. For sake of completeness, let us recall some of the arguments used in the asymptotic analysis of the perturbed steepest descent system when the perturbation g satisfies (2.3). Given z * ∈ S, let T > t 0 be fixed and for every T ≥ t ≥ t 0 consider and thus lim inf t→+∞ t (f (z (t)) − inf H f ) = 0. It remains to show that the limit exists. To this end we compute the time and apply once again [1, Lemma 5.1]. Statement iii) follows from assumption (2.3) and (2.8).
Let L be the Lipschitz constant of ∇f on a ball containing the trajectory z(·). It follows from [10, Lemma 1] that for every t ≥ t 0 which implies that lim t→+∞ t ∇f (z(t)) 2 = 0 and proves i).
The convergence rate becomes (take t = τ (s) in iii) of Theorem 2.1) By introducing a function τ (s) that grows faster than the identity (namely τ (s) ≥ s), we have accelerated the dynamic, passing from the asymptotic convergence rate 1/t for (SD) to 1/τ (s) for (2.12). The price to pay is that we no longer have an autonomous dynamic in (2.12), with as major drawback the fact that the coefficient in front of the gradient term tends to infinity as s → +∞. This prevents from using gradient methods to discretize it. Recall that for gradient methods the step size has to be less than or equal to twice the inverse of the Lipschitz constant of the gradient. To overcome this we come with the second step of our method which is averaging.
Several recent papers have been devoted to the role of the averaging techniques in the acceleration of optimization algorithms. In [43] Scieur, D'Aspremont, and Bach developed a universal acceleration method based on averaging via polynomial techniques. In [37] Poveda and Li developed another universal acceleration method via averaging technique for singularly perturbed hybrid dynamical systems with restarting mechanisms. Let us explain some simple and, as far as we know, new ideas concerning the averaging methods based on differential equations. The use of differential equations allows us to use differential and integral calculus, which turns out to be a flexible approach.
Let us attach to y(·) the new function In doing so, we passed from the first-order differential equation (2.12) to the second-order differential equation (2.18), with the advantage that now the coefficient in front of the gradient is fixed. Let us now particularize the time scale τ (·). Set θ(s) =τ (s). We are led to solve the fist-order linear differential equation After multiplication by s −α , we get equivalently By integrating we get According to θ(s) =τ (s), we finally get This leads to the choice α > 1. A particular and relatively simple situation is to take C 4 = C 5 = 0, which gives .
In doing so, we passed from the first-order differential equation to the second-order differential equation (2.28), with the advantage that now the coefficient in front of the gradient is fixed. Of course, we have to prove that the fast convergence rates are preserved. The above dynamic (2.28) is a particular case of the Inertial System with Implicit Hessian Damping considered by Alecsa, László, and Pinţa in [2], see also Attouch, Fadili, and Kungurtsev [15] in the perturbed case. The rationale justifying the use of the term "implicit" comes from the observation that by Taylor expansion (as s → +∞ we haveẋ(s) → 0) one has ∇f x(s) + β(s)ẋ(s) ≈ ∇f (x(s)) + β(s)∇ 2 f (x(s))ẋ(s), hence making the Hessian damping appear indirectly. We can therefore expect the dynamic (2.28) to have an asymptotic behavior similar to that of (2.31)ẍ(s) + α sẋ (s) + s α − 1 ∇ 2 f (x(s))ẋ(s) + ∇f (x(s)) = 0.

Closed loop control.
The continuous steepest descent system (SD) provides several quantities which are strictly increasing and converge to +∞, so they are eligible for time scaling. Since ż(t) is monotonically decreasing to zero, taking t = τ (s) = 1 ż(s) p is eligible for any p > 0. According to (SD) we have that the same holds true for Since τ comes into the scaling technique via its derivative, we can as well consider for p > 0 t =τ (s) = 1 ż(s) p . Alternatively, we may also consider building τ with the help of the function values since f (z(·)) is decreasing.
All these choices will lead to closed-loop dynamical systems which have attracted of many researchers recently, see for instance [4,38].
2.5. Classical convergence results for the system with explicit Hessian-driven damping. The above dynamic (2.31) is a particular case of the general dynamic on [s 0 , +∞[ studied by Attouch, Chbani, Fadili, Riahi in [10]. Let us recall the convergence results obtained in [10] by performing a Lyapunov analysis for (DIN − AVD) α,β,b . They will serve as a comparison with our results. Then, for every solution trajectory x : [s 0 , +∞[→ H of (DIN − AVD) α,β,b , the following statements are true: Let us specialize the above result to (2.31) by taking b(s) = 1 and β(s) = s α−1 for every s ≥ s 0 . We get w(·) ≡ α−3 α−1 . Thus (G 2 ) is satisfied for α > 3, where as (G 3 ) reduces to 0 ≤ (α−3) 2 α−1 , which is also satisfied. We thus get the convergence rates for the trajectories of (2.31) : 2.6. Convergence rates. The above results suggest that similar results are valid for (2.18). We will examine the convergence rates satisfied by the trajectories of (2.18) simply by using time scale and averaging arguments. We do not make any other Lyapunov analysis, we just use the Lyapunov analysis of the steepest descent, which was recalled in Subsection 2.1. Let us state our first result which concerns the implicit Hessian driven damping dynamics. Then there exists a positive constant C such that the following statements are true: iii) (integral estimate of the gradients) +∞ s0 iv) just assume that α > 1: then x(s) converges weakly as s → +∞, and its limit belongs to S = argmin f .
Proof. First observe that the time scaling and averaging operations allow us to obtain any solution trajectory of (2.32). To this end we need to choose adequately the initial data in the first order evolution equations (2.29) and (2.14) which are respectively attached to the scaling and averaging operations. So let us give a solution trajectory x(·) of (2.32) which satisfies the Cauchy data x(s 0 ) = x 0 andẋ(s 0 ) = x 1 . Let us verify that this solution is reached by taking successively , which, by the second equation of (2.33), leads toẋ(s 0 ) = x 1 . i) Let us first interpret the passage from y to x as an averaging process. To this end we rewrite (2.34) as After multiplication of (2.35) by s α−2 , we get equivalently By integrating (2.37) from s 0 to s, and according to x(s 0 ) = x 0 , we obtain where the last equality follows from the choice of y(s 0 ) as given by (2.33). Then, observe that x(s) can be simply written as follows where µ s is the positive Radon measure on [s 0 , s] defined by is considered as a small perturbation for large values of s. First, let us get rid of the perturbation term. Since y(s) converges weakly as s → +∞ (a property of (SD)), we have that the trajectory y(·) remains bounded.
From (2.38) we deduce that x(·) is also bounded. Since ∇f is Lipschitz continuous on the bounded sets, it follows that ∇f is bounded on the bounded sets, which immediately implies that f is also Lipschitz continuous on the bounded sets. Therefore there exists L f > 0 such that For α > 3 we have α − 1 > 2, and the above estimate gives We are therefore reduced to examining the convergence rate of f s s0 y(u)dµ s (u) − inf f towards zero as s → +∞. We have that µ s is a positive Radon measure on [s 0 , s] whose total mass is equal to 1. It is therefore a probability measure, and s s0 y(u) dµ s (u) is obtained by averaging the trajectory y(·) on [s 0 , s] with respect to µ s . From there, we can deduce fast convergence properties for the solution trajectories of (2.18). According to the convexity of f , and Jensen's inequality, we obtain that where the last inequality above comes from (2.27). According to the definition of µ s , it yields For α > 3 we have We therefore obtain that We have already observed in i) that the trajectory x(·) is bounded. Let L > 0 be the Lipschitz constant of ∇f on a ball that contains the trajectory x(·). It follows from [39, Theorem 2.1.5] and the fact that ∇f (x * ) = 0 for x * ∈ S that for every s ≥ s 0 A recent version of this result which is valid in a general Hilbert space can be found in [10, Lemma 1]. From here we obtain ii).
iii) According to the energy estimates for the continuous steepest descent system, see Theorem 2.1, we have that for any solution trajectory of After making the change of time variable , and setting z(u(s)) = y(s), the above integral estimate becomes Equivalently y(s) = x(s) + s α−1ẋ (s). Replacing y by this equivalent formulation in (2.47) gives the claim iv) For trajectory convergence, following the approach of the paper, we do not perform a Lyapunov analysis, but take advantage of the fact that the trajectory z (s) of the steepest descent dynamics converges weakly towards a solution x * ∈ S as s → +∞. This immediately implies that y(s) = z(u(s)) converges weakly to x * as s → +∞. In other words, for each v ∈ H To pass from the convergence of y to that of x, we use the interpretation of x as an average of y plus a negligeable term. The convergence then results from the general property which says that convergence entails ergodic convergence. Let us make this precise. Since the perturbation term ξ is negligible as s → +∞, we have by definition of x(s) After elementary calculus, we just need to prove that if a(·) is a positive real valued function which verifies lim u→+∞ a(u) = 0, then lim s→+∞ A(s) = 0, where Given > 0, let us take T such that s 0 < T and a(u) ≤ for t ≥ T . For s > T , let us write Letting s go to +∞ we get lim sup This being true for any > 0, we infer that lim s→+∞ A(s) = 0, which gives the claim.
Remark 2.4. Let us emphasize the fact that the convergence of the solution trajectories of the damped inertial system (2.32) is valid under the condition α > 1. This contrasts with the stronger condition α > 3 which is required to obtain the convergence of the trajectories for the system without Hessian driven damping and which is the low-resolution ODE of Nesterov's accelerated gradient method. This property is related to the presence of the Hessian driven damping and the particular form of the corresponding coefficient s α−1 .
2.7. From capital "O" to small "o" convergence rates. It has been established in [19] that taking α > 3 provides the convergence rate of the values o 1/s 2 , thus improving the convergence rate of Nesterov O 1/s 2 . We will see that the time scale and averaging approach is flexible enough to also provide the convergence rate o 1/s 2 in the case α > 3. This is made precise in the following theorem.
Then, as s → +∞ Proof. i) From the convergence properties of the steepest descent, as recalled in Theorem 2.3, we have t for a positive function ε 0 (·) which satisfies lim t→+∞ ε 0 (t) = 0 as t → +∞. From the definition of y(·) we get with ε(s) = ε 0 (τ (s)). Since τ (s) → +∞ as s → +∞, we have that ε(s) tends to zero as s → +∞. Let us show that the property (2.51) is inherited by x(·). For simplicity, we suppose that x(s 0 ) = y(s 0 ), in which case is obtained by averaging the trajectory y(·) on [s 0 , s] with respect to µ s . According to the convexity of f , and Jensen's inequality, we deduce that, for all s ≥ s 0 Let us explicit the last above integral. Recall that τ (s) = s 2 2(α − 1) and Equivalently Therefore, for α > 3 we get So, to prove our claim, we just need to show that lim sup We follow an argument similar to that developed in the proof of Theorem 2.3 iv). Given η > 0, let us take T η such that s 0 < T η and ε(u) ≤ η for u ≥ T η . For s > T η , let us write Letting s go to +∞ we get lim sup This being true for any η > 0, we obtain the claim i).
ii) We first follow an argument similar to the one above. From , we get with ε(s) = 2ε 0 (τ (s)) which tends to zero as s → +∞. Let us now establish a similar estimate for ẋ(s) . According to (2.38) Let us derivate this expression. We geṫ Let us reformulate this expression in terms ofẏ(s), which is the quantity whose speed of convergence is known by (2.54). Indeed, by integration by parts we have After multiplication of the above expression by α−1 s α , and according to y(s 0 ) = x 0 , we get Since ε(s) tends to zero as s → +∞, an argument similar to that used in i) gives lim s→+∞ sẋ(s) = 0. Let us reformulate (2.61) as a differential equation with y(s) as the state variable. According to (2.60) we have which after derivation gives . Combining (2.61) and (2.63) we obtain By replacingẋ(s) by its formulation given in (2.62) we obtain Putting together (2.62) and (2.66) we obtain the system Let us now consider time scaling of the above system. Set A similar calculation as in Subsection 2.2 gives This is a perturbation of the steepest descent dynamical system in the product space H × H. Precisely, setting Z(t) = (Y (t), X(t)) ∈ H × H, we have where Ψ(Z) = f (Y ) and A : H × H → H × H is the linear operator whose matrix is given by The case β 0 = 2 is of particular interest since then A is a symmetric operator which is positive semidefinite.
The study of this system then follows from classical results concerning gradient dynamical systems with multiscale aspects, see for example the work of Attouch and Czarnecki [13].
The dynamical system (2.69) for β 0 = 2 fits in the dynamics (2.70)Ż(t) + β 0 ∇Ψ (Z(t)) + ε (t) ∇Φ (Z(t)) = 0, with ε (t) := (α−1) 2β0t and Φ (·) = 1 2 · 2 A , where · A is the seminorm induced by the symmetric and positive semidefinite operator A. The differential equation (2.70) can be seen as a continuous time approach of the bilevel optimization problem Let Z * be a solution of (2.71), meaning that Φ (Z * ) = inf D Φ and Ψ (Z * ) = inf H×H Ψ. Define for t ≥ t 0 Indeed the above quantity is nonnegative because we are in the particular simple case where Φ(Y, This shows that lim t→+∞ Z (t) − Z * ∈ R exists and that This yields lim inf t→+∞ tΨ Φ (t) = 0. On the other hand, we have d dt This means that Ψ Φ is decreasing and furthermore d dt Since Ψ and Φ are convex and weakly lower semicontinuous, the second condition of Opial's lemma is fulfilled. Thus Z(t) converges weakly to a solution of (2.71) as t → +∞. Moreover, according to the definition of Φ and A, it holds According to the definitions of Ψ and Φ we conclude that there exists x * ∈ argmin H f such that (Y (t), X(t)) converges weakly to (x * , x * ) as t → +∞, This means that x(s) converges weakly to x * as s → +∞, 3. Nonsmooth case. As a basic property, the couple of variables (x, y) defined in Subsection 2.6 satisfies the following differential system, which involves only first-order derivatives in time and space (3.1)   ẏ (s) + s α−1 ∇f (y(s)) = 0 x(s) + α−1 s (x(s) − y(s)) = 0. This naturally suggests the extension of the above results to the nonsmooth case (replace the gradient of f by its subdifferential). Given f : H → R ∪ {+∞} a convex, lower semicontinuous, and proper function, we thus get the differential inclusion system x(s) + α−1 s (x(s) − y(s)) = 0. Solving this system gives generalized solutions to the second-order differential inclusion (s) 0 whose direct study raises several difficulties. To extend the results of the previous section to this nonsmooth case, we need to avoid the arguments using the Lipschitz continuity of ∇f . So, we are led to consider the Cauchy problem for (3.3) with initial data x 0 ∈ dom f and x 1 = 0, that is with initial velocity equal to zero. So doing we have y(t 0 ) = x(t 0 ) = x 0 , which allows to interpret x as an average of y without correcting term (ξ = 0 with the notations of Subsection 2.6). Then, the existence and uniqueness of a strong solution to the associated Cauchy problem relies on the equivalent formulation of (3. , we obtain the existence and uniqueness of a global strong solution of the Cauchy problem associated with (3.2). As a major advantage of the time scaling and averaging techniques, the arguments used in the previous section still work in this more general nonsmooth situation. The rules of differential calculus are still valid for strong solutions, see [31,chapter VIII.2], and Jensen's inequality is still valid for a nonsmooth function f . Indeed Jensen's inequality is classical for a smooth convex function f . Its extension to the nonsmooth convex case can be obtained by first writing it for the Moreau-Yosida regularization f λ of f , then passing to the limit when λ ↓ 0. According to the monotone convergence of f λ towards f , we can pass to the limit in the integral term thanks to the Beppo-Levy monotone convergence theorem. Therefore we obtain the following theorem.  ∀s ≥ s 0 ; ii) just assume that α > 1: then x(s) converges weakly as s → +∞, and its limit belongs to S = argmin f . Remark 3.2. As a consequence of i), we have that x(s) remains in the domain of f for all s ≥ t 0 . This viability property strongly depends on the fact that the initial position belongs to the domain of f , and that the initial velocity has been taken equal to zero. Remark 3.3. The above theorem is valid in an infinite dimensional setting, which makes it applicable to nonlinear PDEs. According to the classical theory for the steepest descent dynamic in the nonsmooth convex case, we have that at each s > s 0 the right derivative of y(s) exists and satisfies − dy ds where (∂f (y(s))) 0 is the element of minimum norm of ∂f (y(s)). Thus, y(·) may exhibit shocks. By contrast, x, which is an average of y, has a continuous derivative, hence does not exhibit shocks. Still x is not twice differentiable.

A proximal type algorithm for nonsmooth minimization.
Recall the second-order differential inclusion (3.3) associated with (1.1) in Section 3 or equivalently formulated as a first order system (see (3.2)) (4.2)   ẏ (s) + s α−1 ∂f (y(s)) 0 x(s) + α−1 s (x(s) − y(s)) = 0. We consider the following implicit discretization in time of (4.2), written for every k ≥ 0 as This gives rise to the following numerical algorithm , where y 0 , x 0 ∈ H are arbitrary initial points. This algorithm has some analogy with the relaxed inertial proximal algorithm for maximally monotone operators in [8].
For the step size sequence (s k ) k≥0 we impose the following recurrence condition which is nothing but the classic step size rule of Nesterov's accelerated gradient method [40,39,24]. We will nevertheless continue to work with the sequence (s k ) k≥0 in order to emphasize that all the convergence statements provided in this section are valid for α > 1, contrary to the stronger condition α > 3 for the low resolution case (see Remark 2.4). We have that s k ∼ k, precisely, s k ≥ (k + 1) (α − 1) /2 for every k ≥ 0 (cf. [24,Lemma 4.3]). Furthermore, by a telescoping sum argument, we have from (4.5) Since s 0 = 0, we can conclude that for every k ≥ 0 (4.6) Let us also notice that (4.4) can be written only in terms of the sequence (x k ) k≥0 which can be seen as a direct discretization of (4.1).
In the following theorem we analyze the convergence properties of the sequence (y k ) k≥0 generated by (4.4). One can easily observe the similarity between the properties of this sequence and those of the trajectory y(·) generated by the continuous steepest descent system after time scaling.
In particular, if f is differentiable, then iv) the sequence of iterates (y k ) k≥0 converges weakly as k → +∞, and its limit belongs to S = argmin H f .
Proof. Let k ≥ 0 be fixed. Take z * ∈ S = argmin f . According to (4.4) there exists η k+1 ∈ ∂f (y k+1 ) such that y k+1 − y k + s k+1 α−1 η k+1 = 0. By convexity of f we deduce that The statements i) and ii) follow from [23,Lemma 5.31]. In addition, the limit lim k→+∞ y k − z * ∈ R exists, which means that the first condition of the discrete Opial's lemma is fulfilled. The According to [7,Lemma 22] we get which proves iii), thanks to (4.6). Finally, since lim k→+∞ f (y k ) = inf H f and f is convex and lower semicontinuous, the second condition of Opial's lemma is also fulfilled. This gives the weak convergence of the sequence (y k ) k≥0 to an element in S = argmin f .
Let us now consider the convergence properties of (x k ) k≥0 which are associated with the averaging process.
Theorem 4.2. Let (x k ) k≥0 be the sequence generated by (4.4). The following statements are true: ii) the sequence of iterates (x k ) k≥0 converges weakly as k → +∞, and its limit belongs to S = argmin H f . Proof. We begin by interpreting (x k ) k≥0 as an average of (y k ) k≥0 . Multiplying the second equation in (4.4) by s 2 k+1 and using (4.5), we obtain for every k ≥ 0 s 2 k+1 x k+1 = s 2 k+1 − (α − 1) s k+1 x k + (α − 1) s k+1 y k+1 = s 2 k x k + (α − 1) s k+1 y k+1 . By a telescopic sum argument and using (4.6), we obtain from here the following expression for each k ≥ 0 Then, Jensen's inequality gives for every k ≥ 0 where in the last equation we use (4.6). In Theorem 4.1 we have shown i≥0 s i (f (y i ) − inf H f ) < +∞, which gives the announced convergence rates. The fact that (x k ) k≥0 converges weakly to an element in S = argmin f as k → +∞ follows from the fact that convergence entails ergodic convergence.

Explicit Hessian damping.
In this section we consider an explicit Hessian driven damping version of the Su-Boyd-Cándes dynamical system (1.3). As explained in (2.31) the following dynamic (∇f (y (s))) + ∇f (y (s)) = 0 comes naturally by performing a Taylor expansion in (2.32). Just note that to match the general writing of inertial dynamics for optimization, the above formula is obtained by taking α + 1 instead of α − 1. This dynamic falls into the setting of (DIN − AVD) α,β,b described in Theorem 2.2. Let us show that this dynamic can be reached by time scaling of a perturbed version of the continuous steepest descent dynamical system. As for the implicit case, we make in the latter the change of time variable τ (s) = s 2 2(α + 1) .
Then there exists a positive constant C such that the following statemnets are true: ∇f (x(s)) ≤ C s ∀s ≥ s 0 ; iii) (integral estimate of the gradients) +∞ s0 s 3 ∇f (y(s)) 2 ds < +∞; iv) the solution trajectory y(s) converges weakly as s → +∞, and its limit belongs to S = argmin f .
6. Combining implicit and explicit Hessian driven damping. As a starting dynamic to develop our time scaling and averaging technique, let us consider in this section the regularized Newton dynamical system It is a special case of the regularized Newton dynamic This system has been studied by Attouch and Svaiter [22], Attouch, Redont, and Svaiter [21], Attouch, Marques Alves, and Svaiter [17] to solve monotone inclusions (that is, with a general maximally monotone operator A instead of ∇f ). In this approach, the central question is the adjustment of the Levenberg-Marquardt regularization parameter λ(·) in front of the velocity term. Indeed, taking lim t→+∞ λ(t) = 0 allows to be asymptotically close to the Newton method. Our situation concerns the simpler case λ(t) ≡ λ > 0 constant, which fits with the convergence properties proved in these papers, and the fact that this dynamic is well-posed. For precise statements concerning the asymptotic behavior of this system see [22,Theorem 4.1,Remark 4.2] and [22,Theorem 3.9]. For the sake of completeness, let us state and prove its convergence properties useful to our study.
Theorem 6.1. Let f : H → R be a convex C 1 function whose gradient is Lipschitz continuous on the bounded sets of H. Let z : [t 0 , +∞[→ H be a solution trajectory of (6.1). Then the following statements are true: iii) (integral estimates of velocities and gradients) +∞ t0 ż (t) 2 dt < +∞ and +∞ t0 v (t) 2 dt < +∞, which lead to v) the solution trajectory z(t) converges weakly as t → +∞, and its limit belongs to S = argmin f .
Proof. Given z * ∈ S, we define the following function for every t ≥ t 0 , which is nonnegative due to the convexity of f and the fact that v (t) = ∇f (z (t)) According to (6.1) and the convexity of f we have for every We integrate from t 0 to t > t 0 and let t tend to +∞, and obtain i).
The inequality above also implies that the function φ is nonincreasing on [t 0 , +∞[, from which we deduce that φ(t) converges as t → +∞, and that the trajectory z(·) is bounded. Let L > 0 be the Lipschitz constant of ∇f on a ball that contains the trajectory z(·). It follows from (6.1) that for every Since t → f (z (t)) − inf H f is nonincreasing, by applying Lemma 9.3 in the Appendix we obtain ii).
In addition, integrating again from t 0 to t > t 0 and letting t tend to +∞, we obtain iii).
Consequently, for item iv) we can use again Lemma 9.3. Precisely, observe that for every t ≥ t 0 d dt which means t → 1 2 v (t) 2 is nonincreasing. The conclusion follows by using (6.3). Finally, let us recall that lim t→+∞ φ (t) ∈ R exists. Moreover, according to ii), iv), and the fact that the trajectory z(·) is bounded, we conlcude that This yields the first condition of Opial's lemma. Now let z be a weak sequential cluster point of the trajectory.
Since v (t) = ∇f (z (t)) for every t ≥ t 0 , it follows from the sequential closedness property of ∇f in the weak × strong topology on H × H that z ∈ S, which proves that the second condition of Opial's lemma is also satisfied. This completes the proof.
The convergence rate becomes s 2 as s → +∞. So, we have accelerated the dynamic for the values, passing from the convergence rate 1/t to 1/s 2 .
Then there exists a positive constant C such that the following statements are true: ∀s ≥ s 0 ; ii) (integral estimate of the gradients) +∞ t0 s w (s) 2 ds < +∞; iii) (convergence of gradients towards zero) w (t) = o 1 s as s → +∞; iv) the solution trajectory x(s) converges weakly as s → +∞, and its limit belongs to S = argmin f .
Proof. Let x(·) be a solution trajectory of (6.15) which satisfies the Cauchy data x(s 0 ) = x 0 ,ẋ(s 0 ) = x 1 , and w (s 0 ) := ∇f x 0 + s0 α−1 x 1 . In the same way as in the proof of Theorem 2.3, we first show that such a solution is reached by considering first the solution y(·) of the system (6.16) Indeed, (6.17) givesẋ(s 0 ) = α−1 s0 (y(s 0 ) − x(s 0 )), which, by the second equation of (6.16), equals to x 1 . We have already seen that x can be seen as an averaging process (6.17) of y as follows (6.18) x where µ s is the measure on [s 0 , s] defined by Since the averaging process is the same as in Theorem 2.3, the convergence rate of 1/s 2 for the values follows by the same arguments.
The statements ii) and iii) follow from Theorem 6.1 after time rescaling and integration by substitution. In order to show iv), we need to pass from the convergence of y to that of x by using the interpretation of x as an average of y plus a negligible term. Moreover, we know that y converges weakly to an element in S thanks to Theorem 6.1. The convergence of x follows from similar arguments as in Theorem 2.3.
7. Operator case. What underlies our study is the first-order in time evolution system taken as a starting point and its convergence properties. In this regard, it is natural to consider also the evolution equation governed by a cocoercive operator. Indeed, gradients of convex functions and cocoercive operators are the two basic examples of monotone operators ensuring the convergence of the associated semi-groups of contractions. We will successively examine the case of a general cocoercive operator, then the additively structured case. In the case of monotone operators, the convergence rates are expressed in terms of velocities and residuals.

General cocoercive case. Recall that a single-valued operator
We mention mention the following classical properties of cocoercive operators: a) If M is ρ-cocoercive, then it is monotone and 1 ρ -Lipschitz continuous, and hence maximally monotone. b) If f is a convex differentiable function whose gradient is L-Lipschitz continuous, then the operator ∇f is 1/L-cocoercive (Baillon-Haddad theorem).
Let us come to our study. Given M : H → H a ρ-cocoercive operator, consider the monotone equation Therefore, we can develop a strategy parallel to that developed in [18], which consists in replacing a maximally monotone operator by its Yosida approximation. This leads to the inertial dynamics (7.19)ÿ (s) + α sẏ (s) + s α + 1 d ds (M A,B,µ (y (s))) + M A,B,µ (y (s)) = 0.
It is immediate to apply Theorem 7.2 to (7.19). This open new perspectives for forward-backward accelerated regularized Newton methods. Ay where A is a linear operator from R n to R m , g : R n → R ∪ {+∞} is a convex, lower semicontinuous and proper function which acts as a regularizer. This class of problems covers many interesting situations arising in the signal and image processing such as the LASSO problem. We are in the setting of (7.18) with M A,B,µ (y) := 1 µ y − prox µg (y + µA * (b − Ay)) , where 0 < µ < 1 A 2 . Following [10], we can also interpret our problem as minimizing the Moreau envelope of f computed for an ad hoc metric. Precisely, if W := 1 µ I−A * A and f W (y) := argmin u∈R n f (u) + 1 2 y − u 2 W , then ∇f W (y) = M A,B,µ (y) and inf H f = inf H f W . So, we can use this example to illustrate also the convergence properties of the dynamic with explicit Hessian damping considered in Section 5. Thus, we can expect a convergence rate of 1/s 2 for f (y W (s)) − inf H f , where y W (s) := prox µg (y (s) + γA * (b − Ay (s))).
We conduct the experiment for a randomly generated matrix A and by taking g to be the 1 norm with various choices for α ∈ {1.001, 2, 3, 5}. Notice that we do not need need to solve (7.19) directly but rather exploit the fact that it can be rewritten as the first order equatioṅ y (s) + s α + 1 M A,B,µ (y (s)) = s α 0 s α y 1 + s 0 α + 1 y 0 , as shown previously. We then use the Runge-Kutta method to solve this equation numerically. Below we plot for the different choices of α the values of f at y W (s) and of the norm of the operator at y (s), respectively.

Conclusion and perspectives.
Our study has shed new light on the acceleration of gradient methods for general convex differentiable optimization. Based on a dynamical approach, and using relatively elementary tools, namely time scaling and time averaging, we were able to introduce a large class of damped inertial dynamics with proven fast convergence rates. The method is very flexible and can be adapted to various initial start-up dynamics. We paid particular attention to the steepest continuous descent and the regularized Newton method. Our approach confirms the important role played by the Hessian driven damping role in these questions. Our results open the door to a systematic study of the corresponding fast proximal and gradient methods. Importantly, our approach also works when we take a cocoercive operator instead of a gradient operator. Applications to inertial forward-backward methods, ADMM, can therefore be expected, to cite some of the most popular algorithms for convex structured optimization. Application to stochastic differential equations can be considered as well.
9. Appendix. Opial's lemma in continuous form is a basic ingredient of the asymptotic analysis of dynamical systems. (ii) every weak sequential limit point of x(t), as t → ∞, belongs to S. Then x(t) converges weakly as t → +∞, and its limit belongs to S.