Camille Castera
email: camille.castera@protonmail.com

Hedy Attouch

Jalal Fadili

Peter Ochs

Continuous Newton-like Methods featuring Inertia and Variable Mass

We introduce a new dynamical system, at the interface between second-order dynamics with inertia and Newton's method. This system extends the class of inertial Newton-like dynamics by featuring a time-dependent parameter in front of the acceleration, called variable mass. For strongly convex optimization, we provide guarantees on how the Newtonian and inertial behaviors of the system can be non-asymptotically controlled by means of this variable mass. A connection with the Levenberg-Marquardt (or regularized Newton's) method is also made. We then show the effect of the variable mass on the asymptotic rate of convergence of the dynamics, and in particular, how it can turn the latter into an accelerated Newton method. We provide numerical experiments supporting our findings. This work represents a significant step towards designing new algorithms that benefit from the best of both first-and second-order optimization methods.

Introduction 1.Problem Statement

A major challenge in modern unconstrained convex optimization consists in building fast algorithms while maintaining low computational cost and memory footprint. This plays a central role in many key applications such as large-scale machine learning problems or data processing. The problems we are aiming to study are of the form min x∈R n f (x).

Large values of n demand for algorithms at the interface of first-and second-order optimization. Limited computational capabilities explain why gradient-based (first-order) algorithms remain prominent in practice. Unfortunately, they often require many iterations which is true even for the provably best algorithms for certain classes of optimization problems; for example that of convex and strongly convex functions with Lipschitz continuous gradient [START_REF] Boris | Some methods of speeding up the convergence of iteration methods[END_REF][START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence O(1/k 2 )[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. On the other hand, algorithms using second-order information (the Hessian of f )-with Newton's method as prototype-adapt locally to the geometry of the objective, allowing them to progress much faster towards a solution. However, each iteration comes with high computational and memory costs, which highlights a challenging tradeoff. It is therefore essential to develop algorithms that take the best of both worlds. They are commonly referred to as (limited-memory) quasi-Newton methods. Several quasi-Newton algorithms partly address this issue, for example BFGS methods [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF][START_REF] Fletcher | A new approach to variable metric algorithms[END_REF][START_REF] Goldfarb | A family of variable-metric methods derived by variational means[END_REF][START_REF] David | Conditioning of quasi-Newton methods for function minimization[END_REF][START_REF] Dong | On the limited memory BFGS method for large scale optimization[END_REF]], yet, in very large-scale applications, first-order algorithms often remain the preferred choice.

In order to reach a new level of efficiency, deep insights into the mechanism and relations between algorithms are required. To that aim, an insightful approach is to see optimization algorithms as discretization of ordinary differential equations (ODEs): for small-enough step-sizes, iterates can be modeled by a continuous-time trajectory [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF][START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF]. Obtaining a fast algorithm following this strategy depends on two ingredients: choosing an ODE for which rapid convergence to a solution can be proved, and discretizing it with an appropriate scheme that preserves the favorable properties of the ODE.

Both steps are highly challenging, our work focuses on the ODE matter. We study the following second-order dynamical system in a general setting:

ε(t)ẍ(t) + α(t) ẋ(t) + β∇ 2 f (x(t)) ẋ(t) + ∇f (x(t)) = 0, t ≥ 0, (VM-DIN-AVD)
where f : R n → R is a smooth convex twice continuously differentiable function, with gradient ∇f and Hessian ∇ 2 f defined on R n equipped with scalar product •, • , and induced norm • . Additionally, f is assumed to be coercive, and strongly convex on bounded subsets of R P . The functions ε, α : R + → R + (where R + = [0, +∞[) are differentiable, non-increasing, and ε(t) > 0 for all t ≥ 0. Together with β > 0, they are control parameters that define the type of dynamics that drives the trajectory (or solution) x : R + → R n , whose first-and second-order derivatives are denoted ẋ and ẍ respectively. We call the above dynamics (VM-DIN-AVD), which stands for "Variable Mass Dynamical Inertial Newton-like system with Asymptotically Vanishing Damping" since it generalizes a broad class of ODEs whose original member is DIN [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF], where ε and α were constant. DIN was then extended to the case of non-constant asymptotically vanishing dampings (AVD) α [START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF]. In this work we introduce the non-constant parameter ε called variable mass (VM) in front of the acceleration ẍ, in the same way that α is called (viscous) damping by analogy with classical mechanics. A key feature of these ODEs, that positions them at the interface of first-and second-order optimization, is that they possess equivalent forms involving only ∇f but not ∇ 2 f , significantly reducing computational costs, hence enabling the design of practical algorithms, see e.g., [START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF][START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF][START_REF] Chen | First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow[END_REF]. The key idea behind this is the relation ∇ 2 f (x(t)) ẋ(t) = d dt ∇f (x(t)), see Section 2 for an equivalent formulation of (VM-DIN-AVD) exploiting this. This paper emphasizes the relation between (VM-DIN-AVD) and well-studied special cases. Indeed, taking ε = α = 0, one obtains 1 the Continuous Newton (CN) method [START_REF] Konstantinovich | Nonlinear functional equations and continuous analogues of iteration methods[END_REF] β∇ 2 f (x N (t)) ẋN (t) + ∇f (x N (t)) = 0, t ≥ 0, (CN) For each patch, the color indicates which of the distances x N (t) -x(t) and x LM (t) -x(t) is considered, the scaling of a corresponding upper-bound on this distance is written; in white for prior work and in black for our contributions. The green line separates the cases ε ≥ α (above) and ε ≤ α (below). Right: 2D illustration of the trajectories of (VM-DIN-AVD) for several choices of ε on a quadratic function. Using fast-vanishing ε(t) (dark-blue solid curves), one can bring the solution of (VM-DIN-AVD) close to that of (CN), making it, for example, more robust to bad conditioning compared to first-order dynamics (such as gradient descent).

known notably for being invariant to affine transformations and yielding fast vanishing of the gradient (see Section 3). In fact, this observation shows that (VM-DIN-AVD) is a singular perturbation of (CN), which also justifies the terminology "Newton-like" in DIN. When α = 0 but ε = 0, we recover the Levenberg-Marquardt (LM) method,

α(t) ẋLM (t) + β∇ 2 f (x LM (t)) ẋLM (t) + ∇f (x LM (t)) = 0, t ≥ 0, (LM)
also known as regularized Newton method since it stabilizes (CN). In the rest of the paper, the solutions of (CN) and (LM) will always be denoted by x N and x LM respectively. Alvarez et al. [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] showed that for α = 0, β = 1, and ε constant and small, (VM-DIN-AVD) is a "perturbed" Newton method since the distance between the solutions of (VM-DIN-AVD) and (CN) is at most proportional to √ ε at all time. Yet, despite the benefits of this class of ODEs, such as stabilization properties, see e.g., [START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF][START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF], no improvement2 of the rate of convergence (in values) has been shown compared to inertial first-order dynamics [START_REF] Boris | Some methods of speeding up the convergence of iteration methods[END_REF][START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]. This raises the question: "are these ODEs really of Newton type?", which is crucial in view of designing faster algorithms from them.

Main Contributions

We show that the answer to this question is partially positive, and closely related to the choices of ε and α. We provide general results on the role played by these two control parameters and how they can be chosen to control (VM-DIN-AVD), and make it close to (CN) for all time, as illustrated on the right-hand side of Figure 1, but also to obtain fast convergence. This represents a first step towards building new fast practical algorithms. Our main contributions are the following:

-We provide a first-order equivalent formulation of (VM-DIN-AVD), and address the questions of existence and uniqueness of the solutions of (VM-DIN-AVD) under mild assumptions.

-We generalize the perturbed Newtonian property discussed above to non-constant and possibly vanishing variable masses ε, and "not too large" positive dampings α, and derive bounds that (formally) take the form x(t) -x N (t) = O( ε(t)). We then extend these results to larger dampings α and make the connection between (VM-DIN-AVD) and (LM). This contribution is summarized in the phase diagram of Figure 1.

-Using quadratic functions as a model for strongly convex functions, we shed light on techniques to efficiently approximate solutions of (VM-DIN-AVD). We then show how ε and α affect the speed of convergence. Depending on their setting, the solutions of (VM-DIN-AVD) may either converge as fast as that of (CN), faster, or rather have a (LM) nature, as summarized in Table 1.

-We provide numerical experiments supporting our theoretical findings.

Related work

The system (VM-DIN-AVD) belongs to the class of inertial systems with viscous and geometric ("Hessian-driven") dampings, initially introduced with constant ε = 1 and constant α in [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] and called DIN (for Dynamical Inertial Newton-like system). Except in a few cases [START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF][START_REF] Castera | Inertial Newton algorithms avoiding strict saddle points[END_REF], most of the follow-up work then considered extensions of DIN with non-constant AVD α, with in particular the DIN-AVD system with α(t) = α 0 /t as introduced in [START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF]. The reason for this popular choice for α is its link with Nesterov's method [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights[END_REF]. Non-constant choices for β have been considered [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF][START_REF] Daniel Alecsa | An extension of the second order dynamical system that models Nesterov's convex gradient method[END_REF][START_REF] Lin | A control-theoretic perspective on optimal high-order optimization[END_REF][START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF]. We keep it constant here, and rather focus on non-constant ε, unlike prior work that used constant ε = 1. The mass ε was only considered in the original work [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF], but only for fixed ε, β = 1 and constant α = 0. VM-DIN-AVD is however closely related to the IGS system considered in [START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF] as it is actually equivalent to the latter after dividing both sides of (VM-DIN-AVD) by ε(t). Our approach-which consists in studying the connections with other second-order dynamics as ε vanishes asymptotically-is however different from the one followed in [START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF], which is of independent interest. The literature on DIN is rich, let us mention further connections with Nesterov's method [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF][START_REF] Daniel Alecsa | An extension of the second order dynamical system that models Nesterov's convex gradient method[END_REF], extensions with Tikhonov regularization [START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian driven damping[END_REF] and closed-loop damping [START_REF] Attouch | Fast optimization via inertial dynamics with closed-loop damping[END_REF][START_REF] Lin | A control-theoretic perspective on optimal high-order optimization[END_REF]. The non-smooth and possibly non-convex cases have been considered in [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[END_REF][START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF]. Finally, avoidance of strict saddle points in smooth non-convex optimization has been shown in [START_REF] Castera | Inertial Newton algorithms avoiding strict saddle points[END_REF].

The influence of the damping α on the (LM) dynamics has been studied in [START_REF] Attouch | A continuous dynamical Newton-like approach to solving monotone inclusions[END_REF][START_REF] Attouch | Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in hilbert spaces[END_REF]. Interestingly, the conditions enforced on α in these papers (formally a sub-exponential decay) are very similar to those we make on ε and α for (VM-DIN-AVD) (see Assumptions 1 and 2).

Regarding the second part of our analysis, which deals with the case where f is quadratic, Attouch et al. [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF] provided closed-form solutions to (VM-DIN-AVD) for ε ≡ 1 and special choices of α. Our work rather deals with approximate solutions which allows considering a wide class of functions ε and α. We rely on the Liouville-Green (LG) method [START_REF] Liouville | Mémoire sur le développement des fonctions ou parties de fonctions en séries dont les divers termes sont assujétis à satisfaire à une même équation différentielle du second ordre, contenant un paramètre variable[END_REF][START_REF] Green | On the motion of waves in a variable canal of small depth and width[END_REF] presented in Section 4. Generalizations of

LG are also often referred to as WKB methods [START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik[END_REF][START_REF] Anthony | Wellenmechanik und halbzahlige Quantisierung[END_REF][START_REF] Brillouin | Remarques sur la mécanique ondulatoire[END_REF] and seem to be mostly used in physics so far. To the best of the authors' knowledge, the current work seem to be one of the first to use the LG method in optimization, and the first for DIN-like systems.

Organization

The paper is organized as follows. We discuss the existence of solutions in Section 2. Our main results, from a non-asymptotic control perspective, are then presented in Section 3. An analysis of the role played asymptotically by ε and α is then carried out on quadratic functions in Section 4. Finally, numerical experiments are presented in Section 5, and some conclusions are then drawn.

Existence and Uniqueness of Solutions

In the sequel, we fix x 0 ∈ R n and ẋ0 ∈ R n , such that, unless stated otherwise, (VM-DIN-AVD) is always considered with initial condition (x(0), ẋ(0)) = (x 0 , ẋ0 ), and (CN) and (LM) with initial condition x N (0) = x LM (0) = x 0 . We also fix initial values for the control parameters ε(0) = ε 0 > 0, ε (0) = ε 0 ≤ 0 and α(0) = α 0 ≥ 0. In addition to the definitions of ε and α in Section 1.1, we assume that ε is twice differentiable, with bounded second derivative. In order to use the Cauchy-Lipschitz Theorem, we reformulate (VM-DIN-AVD) into a first-order (in time) system by introducing an auxiliary variable y : R + → R n . Notably, our reformulation does not involve ∇ 2 f , in the same fashion as [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF][START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF]. For all t, defining ν(t) = α(t) -ε (t) -1 β ε(t), we show in Appendix A that (VM-DIN-AVD) is equivalent to

ε(t) ẋ(t) + β∇f (x(t)) + ν(t)x(t) + y(t) = 0 ẏ(t) + ν (t)x(t) + ν(t) β x(t) + 1 β y(t) = 0 (gVM-DIN-AVD) with initial conditions (x(0), y(0)) = x 0 , -ε 0 ẋ0 -β∇f (x 0 ) -(α 0 -ε 0 -1 β ε 0 )x 0 .
One can notice that in the special case where ε is taken constant and equal to 1 (that is when (VM-DIN-AVD) is simply the DIN-AVD system [START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF]), we recover the same first-order formulation as that in [START_REF] Attouch | Fast convex optimization via inertial dynamics with Hessian driven damping[END_REF].

We can then apply the Cauchy-Lipschitz Theorem. For all t ≥ 0 and (u, v) ∈ R n × R n , define the mapping

G (t, (u, v)) = 1 ε(t) (-β∇f (u) -ν(t)u -v) -ν (t)u -ν(t) β u -1 β v
, so that (gVM-DIN-AVD) rewrites ( ẋ(t), ẏ(t)) = G (t, (x(t), y(t))) for all t ≥ 0. Since f is twice continuously differentiable, one can see that G is continuously differentiable w.r.t. its second argument (u, v). Consequently G is locally Lipschitz continuous w.r.t. (u, v) and by the Cauchy-Lipschitz Theorem, for each initial condition, there exists a unique local solution to (gVM-DIN-AVD) and thus to (VM-DIN-AVD). We then show that this solution is actually global (in Appendix A) by proving the boundedness of (x, y). We omit the existence and uniqueness of the solutions of (CN) and (LM) since these are standard results, obtained with similar arguments.

3 Non-asymptotic Control of (VM-DIN-AVD)

The purpose of this section is to understand how close x might be to x N and x LM , as a function of α and ε. Since f is coercive and strongly convex on bounded sets, it has a unique minimizer x ∈ R n . Consequently, any two trajectories that converge to x will eventually be arbitrarily close to each other. Thus, asymptotic results of the form x(t) -x N (t) ----→ t→+∞ 0 are not precise enough to claim, for example, that x has a "Newtonian behavior". Instead, we will derive upper bounds on the distance between trajectories that hold for all time t ≥ 0, and which typically depend on ε and/or α. We first present the case where α is small relative to ε and then generalize.

Comparison with (CN) under Moderate Viscous Dampings

When the viscous damping α remains moderate w.r.t. the variable mass ε, one can expect the solutions of (VM-DIN-AVD) to be close to that of (CN). We make the following assumptions. Assumption 1. There exists c 1 , c 2 ≥ 0 such that for all t ≥ 0, |ε (t)| ≤ c 1 ε(t) and α(t) ≤ c 2 ε(t).

The assumption states that α must decrease at least as fast as ε (up to a constant). 3 The reason for assuming |ε (t)| ≤ c 1 ε(t) is technical and will appear more clearly in the proofs below. It formally means that ε can decrease at most exponentially fast. 4 This is a relatively mild assumption that holds, for example, for any polynomial decay ε 0 /(t + 1) a , a ∈ N.

We start with the main result of this section. Theorem 3.1. Let x N be the solution of (CN), and let c 1 , c 2 ≥ 0. There exist C 0 , C 1 , C 2 ≥ 0, depending on c 1 , c 2 , such that for all (ε, α) for which Assumption 1 holds with constants c 1 and c 2 , the corresponding solution x of (VM-DIN-AVD) is such that for all t ≥ 0,

x(t) -x N (t) ≤ C 0 e -t β ε 0 ẋ0 + C 1 ε(t) + C 2 t s=0 e 1 β (s-t) ε(s) ds. (1)
This extends a previous result from [2, Proposition 3.1] which states a similar bound for constant ε, α ≡ 0 and β = 1. Theorem 3.1 corresponds to the blue parts in the phase diagram of Figure 1 (see also Corollary 3.6 below). Remark 3.2. The strength of the above result comes from the fact that the constants C 0 , C 1 , C 2 do not depend on ε and α, and that the result is non-asymptotic. This allows in particular choosing (ε, α) to control the distance from x to x N , for all time t ≥ 0.

Remark 3.3. Under Assumption 1, the dynamics (VM-DIN-AVD) is dominated by the variable mass. The damping α does not appear in Theorem 3.1.

The above theorem and remarks emphasize the "Newtonian nature" of (VM-DIN-AVD). We present two lemmas before proving Theorem 3.1, and then state a simpler bound than (1), see Corollary 3.6. Lemma 3.4. Let (ε, α), and let x be the corresponding solution of (VM-DIN-AVD). For all t ≥ 0, define the function,

U (t) = ε(t) 2 ẋ(t) 2 + f (x(t)) -f (x ).
Then U is differentiable and for all t > 0,

dU dt (t) = ε (t) 2 ẋ(t) 2 -α(t) ẋ(t) 2 -β ∇ 2 f (x(t)) ẋ(t), ẋ(t) ≤ 0.
Therefore, in particular, U is non-increasing.

Proof. Let t ≥ 0, since x is twice differentiable, U is differentiable and,

dU dt (t) = ε (t) 2 ẋ(t) 2 + ε(t) ẋ(t), ẍ(t) + ẋ(t), ∇f (x(t)) .
We use the fact that x is solution of (VM-DIN-AVD), to substitute ε(t)ẍ(t) by its expression,

dU dt (t) = ε (t) 2 ẋ(t) 2 -α(t) ẋ(t) 2 -β ∇ 2 f (x(t)) ẋ(t), ẋ(t) .
By assumption ε is non-increasing so for all t > 0, ε (t) ≤ 0. Furthermore f is convex so

∇ 2 f (x(t)) ẋ(t), ẋ(t) ≥ 0. Hence U is non-increasing.
We then state the following bound. Lemma 3.5. There exists C ≥ 0 such that for any (ε, α) and the corresponding solution x of (VM-DIN-AVD), for all t ≥ 0 it holds that,

ε(t) ẋ(t) ≤ C ε(t).
Proof. Let t ≥ 0, according to Lemma 3.4, U is non-increasing so U (t) ≤ U (0), or equivalently,

ε(t) 2 ẋ(t) 2 + f (x(t)) -f (x ) ≤ ε 0 2 ẋ0 2 + f (x 0 ) -f (x ).
This implies in particular that,

ε(t) ẋ(t) 2 ≤ ε 0 ẋ0 2 + 2(f (x 0 ) -f (x )),
and hence by multiplying both sides by ε(t) and composing with the square-root we obtain that,

ε(t) ẋ(t) ≤ C ε(t),
where

C = ε 0 ẋ0 2 + 2(f (x 0 ) -f (x )).
Proof of Theorem 3.1. Let (ε, α) as defined in Sections 1.1 and 2, and let x be the corresponding solution of (VM-DIN-AVD). Then, according to Lemma 3.4, for all t ≥ 0, U (t) ≤ U (0), so in particular

f (x(t)) ≤ f (x 0 ) + ε 0 2 ẋ0 2 . Denoting M 0 = f (x 0 ) + ε 0 2 ẋ0 2 , the set K 0 = {y ∈ R n | f (y) ≤ M 0 } is bounded, since f is coercive (lim y →+∞ f (y) = +∞).
So for all t ≥ 0, x(t) ∈ K 0 . Since M 0 (and hence K 0 ) depends only on ε 0 , x 0 and ẋ0 , we have proved that for any choice (ε, α), the corresponding solution x of (VM-DIN-AVD) is inside K 0 at all time. Let x N be the solution of (CN). One can see similarly that for all t ≥ 0,

f (x N (t)) ≤ f (x N (0)) = f (x 0 ) ≤ M 0 . So we also have x N (t) ∈ K 0 for all t ≥ 0.
Now, fix c 1 , c 2 > 0, and let (ε, α) such that Assumption 1 is satisfied with constants c 1 , c 2 . Let x be the corresponding solution of (VM-DIN-AVD). Since f is strongly convex on bounded sets, it is strongly convex on K 0 . We denote µ > 0 the strong-convexity parameter of f on K 0 . Equivalently, we have that ∇f is strongly monotone on K 0 , that is,

∀y 1 , y 2 ∈ K 0 , ∇f (y 1 ) -∇f (y 2 ), y 1 -y 2 ≥ µ y 1 -y 2 2 . ( 2
)
Let t ≥ 0, since x(t) ∈ K 0 and x N (t) ∈ K 0 , by combining ( 2) with the Cauchy-Schwarz inequality, we deduce that

x(t) -x N (t) ≤ 1 µ ∇f (x(t)) -∇f (x N (t)) . (3)
Therefore, it is sufficient to bound the difference of gradients in order to bound x(t) -x N (t) . First, remark that (CN) can be rewritten as follows:

d dt ∇f (x N (t)) + 1 β ∇f (x N (t)) = 0. So we can integrate, for all t ≥ 0, ∇f (x N (t)) = e -t β ∇f (x 0 ). (4)
We now turn our attention to ∇f (x(t)), for which we cannot find a closed-form solution in general.

We rewrite (VM-DIN-AVD) in the following equivalent form

d dt [ε(t) ẋ(t) + β∇f (x(t))] + 1 β ε(t) ẋ(t) + ∇f (x(t)) = 1 β ε(t) + ε (t) -α(t) ẋ(t).
Introducing the variable ω(t) = ε(t) ẋ(t) + β∇f (x(t)), the latter is thus solution to

ω(t) + 1 β ω(t) = 1 β ε(t) + ε (t) -α(t) ẋ(t), t ≥ 0, ω(0) = ε 0 ẋ0 + β∇f (x 0 ).
This is a non-homogeneous first-order ODE in ω, whose solution can be expressed using the integrating factor

ω(t) = e -t β (ε 0 ẋ0 + β∇f (x 0 )) + e -t β t 0 e s β 1 β ε(s) + ε (s) -α(s) ẋ(s) ds.
We thus have the following expression for ∇f (x), for all t ≥ 0,

β∇f (x(t)) = βe -t β ∇f (x 0 ) + e -t β ε 0 ẋ0 -ε(t) ẋ(t) + e -t β t 0 e s β 1 β ε(s) + ε (s) -α(s) ẋ(s) ds. (5)
We can now use ( 4) and ( 5) in (3) to get

x(t) -x N (t) ≤ 1 βµ e -t β ε 0 ẋ0 -ε(t) ẋ(t) + e -t β t 0 e s β 1 β ε(s) + ε (s) -α(s) ẋ(s) ds .
Using the triangle inequality, we obtain,

x(t) -x N (t) ≤ ε 0 ẋ0 βµ e -t β + ε(t) ẋ(t) βµ + 1 βµ t 0 e 1 β (s-t) 1 β ε(s) + ε (s) -α(s) ẋ(s) ds. ( 6
)
The first term in [START_REF] Attouch | Continuous Newton-like inertial dynamics for monotone inclusions[END_REF] corresponds to the first one in (1) with C 0 = 1/(βµ). As for the second-one, by direct application of Lemma 3.5, there exists C > 0 such that for all t ≥ 0, ε(t) ẋ(t) βµ ≤ C ε(t), so we set C 1 = C/(βµ). Regarding the last term in (6), using Assumption 1 and again Lemma 3.5, it holds that, for all s ≥ 0,

1 β ε(s) + ε (s) -α(s) ẋ(s) ≤ 1 β + c 1 + c 2 ε(s) ẋ(s) ≤ 1 β + c 1 + c 2 C ε(s).
This proves the theorem with

C 2 = C βµ 1 β + c 1 + c 2 .
Let us analyze the bound in Theorem 3.1. The first term in (1) decays exponentially fast and can even be zero if the initial speed is ẋ0 = 0, the second-one decays like ε(t), however, the rate at which the last term decreases is less obvious. The following corollary gives a less-tight but easier-to-understand estimate. Corollary 3.6. Consider the same assumptions and variables as in Theorem 3.1. If furthermore c 1 < 2 β , then there exists C 3 > 0 such that for all t ≥ 0,

x(t) -x N (t) ≤ C 0 e -t β ε 0 ẋ0 + C 3 ε(t).
Proof of Corollary 3.6. For all t ≥ 0, define J(t) = t 0 e s β ε(s) ds. An integration by parts yields

J(t) = βe s β ε(s) t s=0 - t s=0 βe s β ε (s) 2 ε(s) ds = βe t β ε(t) -βε 0 + t s=0 βe s β -ε (s) 2ε(s) ε(s) ds. (7)
By assumption, 0

≤ c 1 < 2 β such that for all s > 0, |ε (s)| ≤ c 1 ε(s), which in our setting is equivalent to -ε (s) ε(s) ≤ c 1 .
So we deduce from (7) that

J(t) ≤ βe t β ε(t) + c 1 β 2 t s=0 e s β ε(s) ds = βe t β ε(t) + c 1 β 2 J(t). So, 1 -c 1 β 2 J(t) ≤ βe t ε(t). By assumption 1 -c 1 β 2 > 0, therefore, J(t) ≤ 2 2-c 1 β e t β
ε(t). Finally, using this in (1) and setting

C 3 = C 1 + C 2 2 2-c 1 β
, we obtain the result. Remark 3.7. The above proofs suggest that an extension to the case where f is non-smooth but strongly convex is possible using regularization techniques. This is left for future work.

So far our results only cover the case where α is "not too large" w.r.t. ε, and do not study (LM). We now state a more general result that covers these cases.

Generalization to Arbitrary Viscous Dampings with Sub-exponential Decay

This time we do not assume a link between ε and α but only sub-exponential decays. Assumption 2. There exists c 1 , c 2 ≥ 0 such that for all t ≥ 0, |ε (t)| ≤ c 1 ε(t) and |α (t)| ≤ c 2 α(t).

We are now in position to state the main result of this section. Theorem 3.8. Let x N and x LM be the solution of (CN) and (LM) respectively, and let c 1 , c 2 ≥ 0. There exist constants C, C ≥ 0, depending on c 1 , c 2 , such that for all ε and α for which Assumption 2 holds with c 1 and c 2 , the corresponding solution x of (VM-DIN-AVD) is such that for all t ≥ 0,

x(t) -x N (t) ≤ C e -t β + ε(t) + α(t) + t s=0 e 1 β (s-t) ( ε(s) + α(s)) ds , (8)
and,

x(t) -x LM (t) ≤ C e -t β + ε(t) + α(t) + t s=0 e 1 β (s-t) ( ε(s) + α(s)) ds . ( 9
)
The proof is postponed to Appendix B. Although it follows a similar reasoning as that of Theorem 3.1, more involved estimates are needed.

Let us comment on these results. The bound (8) generalizes Theorem 3.1, although the constant involved will, in general, be larger than those in (1) (see the proof of Theorem 3.8 in appendix). Theorem 3.8 allows for far more flexibility in the choice of ε and α in order to control x and make it possibly close to x N . The bound in ( 9) is the same as that in (8) (up to a constant), but this time w.r.t. x LM , thus connecting (VM-DIN-AVD) to (LM). We make the following two important remarks. First (9) involves α, suggesting that making x close to x LM requires not only ε but also α to vanish asymptotically. Additionally, Theorem 3.8 does not state to which of x N and x LM the solution of (VM-DIN-AVD) is the closest.

It remains an open question to know whether one can make (9) independent of α, and to state to which trajectory x is the closest. Yet, the numerical experiments in Section 5 suggest that neither are possible. Indeed, we observe that for some functions f , x is sometimes closer to x N than to x LM , even when ε(t) ≤ α(t).

Nevertheless, Theorem 3.8 answers the question asked in the introduction: yes, (VM-DIN-AVD) is really of second-order nature since it can be brought close to the second-order dynamics (CN) and (LM). Doing so, it benefits from the good properties of these methods, such as the robustness to bad conditioning, as previously illustrated on the right of Figure 1. This concludes the analysis from a control perspective. We will now derive an approximation of the solution x in order to study the impact that ε and α have on the speed of convergence of x to x compared to the speeds of convergence of x N and x LM .

Approximate Solutions and Asymptotic Analysis on Quadratic Functions

We consider the particular case where f is a strongly-convex quadratic function in order to study the asymptotic behavior of (VM-DIN-AVD) w.r.t. (CN) and (LM). Quadratic functions are the prototypical example of strongly-convex functions. In particular, any strongly-convex function can be locally approximated by a quadratic one around its minimizer, making the latter a good model for understanding the local behavior of dynamics. In this section, f is quadratic:

f (y) = 1 2 Ay -b 2 2 for all y ∈ R n ,
where A ∈ R n×n is symmetric positive definite and b ∈ R n . Without loss of generality, we take b = 0, so that the unique minimum is x = 0.

Setting: the Special Case of Quadratic Functions

Quadratic functions are particularly interesting in our setting since DIN-like ODEs take a simpler form (as observed in [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF][START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF]). Indeed, ∀y ∈ R n , ∇f (y) = A T Ay and ∇ 2 f (y) = A T A. Since ∇ 2 f (y) is independent of y we can rewrite (VM-DIN-AVD) in an eigenspace5 of A T A. That is, we can study the ODE coordinate-wise by looking at one-dimensional problems of the form

ε(t)ẍ(t) + (α(t) + βλ) ẋ(t) + λx(t) = 0, t ≥ 0. (Q1-VM-DIN-AVD)
Here (and throughout what follows) λ > 0 denotes any eigenvalue of A T A and x : R + → R now denotes the corresponding coordinate (function) of the solution of (VM-DIN-AVD) in an eigenspace of A T A. The dynamics (Q1-VM-DIN-AVD) is a linear second-order ODE in x with non-constant coefficients. Similarly, (LM) can be rewritten coordinate-wise as

(α(t) + βλ) ẋLM (t) + λx LM (t) = 0, t ≥ 0, (Q1-LM)
where x LM : R + → R, and (CN) becomes

β ẋN (t) + x N (t) = 0, t ≥ 0, (Q1-CN)
where again, x N : R + → R is one-dimensional. Observe in particular that (CN) and (LM) are now first-order linear ODEs, whose solutions have the closed forms, for all t ≥ 0,

x N (t) = x 0 e -t β and x LM (t) = x 0 exp -

t 0 λ α(s) + βλ ds . ( 10
)
Since the minimizer is x = 0, we directly see that x N converges exponentially fast to x , with a rate independent of λ. The rate of x LM depends however on λ and how fast α vanishes.

Unfortunately, except for some special choices of ε and α (see [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF]), one cannot solve the second-order linear ODE (Q1-VM-DIN-AVD) in closed form in general. Additionally, it is hopeless to circumvent the difficulty by finding a closed form for ∇f (x), accordingly to what we did in Section 3, since here ∇f (x) = λx. In order to study the speed of convergence of x despite not having access to a closed form, we will approximate it with a controlled error, via a method that we now present.

The Liouville-Green Method

In what follows, we rely on the Liouville-Green method [START_REF] Liouville | Mémoire sur le développement des fonctions ou parties de fonctions en séries dont les divers termes sont assujétis à satisfaire à une même équation différentielle du second ordre, contenant un paramètre variable[END_REF][START_REF] Green | On the motion of waves in a variable canal of small depth and width[END_REF], a technique for obtaining nonasymptotic approximations to solutions of linear second-order ODEs with non-constant coefficients. First, we give the intuition behind the method, following the presentation of [START_REF] Olver | Asymptotics and special functions[END_REF]. Consider the differential equation

z(t) -r(t)z(t) = 0, t ≥ 0, ( 11
)
where r is real-valued, positive, and twice continuously differentiable. Any linear second-order ODE can be reformulated in the form [START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF], see Lemma 4.5 below. Since for all t ≥ 0, r(t) = 0, we can use the changes of variables τ = t 0 r(s) ds and w = r 1/4 z and show that w is solution to

ẅ(τ ) -(1 + ψ(τ ))w(τ ) = 0, t ≥ 0, (12)
where6 ψ(τ ) = 4r(t)r (t)-5r (t) 2 16r(t) 3 . The LG method consists in neglecting the term ψ(τ ) in [START_REF] Attouch | Fast optimization via inertial dynamics with closed-loop damping[END_REF], which simply yields two approximate solutions ŵ1 (τ ) = e τ and ŵ2 (τ ) = e -τ . Expressing this in terms of z and t, we obtain

ẑ1 (t) = r(t) -1/4 exp t 0 r(s) ds and ẑ2 (t) = r(t) -1/4 exp t 0 -r(s) ds . ( 13
)
Those are the LG approximations of the solutions of [START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF]. They are formally valid on any interval [0, T ], T > 0 when ψ is "not too large", provided that √ r is integrable on [0, T ].

Remark 4.1. There exists other (but less intuitive) ways to derive the approximations above, which allow for generalization to higher-order linear ODEs, see e.g., [START_REF] Carl | Asymptotic methods and perturbation theory[END_REF]Chapter 10].

The advantage of this approach is the possibility to estimate the error made using (13) w.r.t. the true solutions of [START_REF] Attouch | Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling[END_REF]. This is expressed in the following theorem which gathers results from [START_REF] Blumenthal | Über asymptotische Integration linearer Differentialgleichungen mit Anwendung auf eine asymptotische Theorie der Kugelfunktionen[END_REF][START_REF] William | Error bounds for the Liouville-Green (or WKB) approximation[END_REF][START_REF] James | Improved error bounds for the Liouville-Green (or WKB) approximation[END_REF].

Theorem 4.2 (Olver [START_REF] Olver | Asymptotics and special functions[END_REF]). Let r : R + → R be a real, positive, twice continuously differentiable function, and define ϕ(t) = 4r(t)r (t)-5r (t) 2 16r(t) 5/2 for all t ≥ 0. Then for any T > 0, the differential equation,

z(t) -r(t)z(t) = 0, t ∈ [0, T ], (14)
has two real and twice continuously differentiable solutions defined for all t ∈ [0, T ] by, -Since we assumed r to be twice continuously differentiable and positive, ϕ is continuous, so it is integrable except maybe for t → +∞.

z 1 (t) = 1 + δ 1 (t) r(t)
-For the sake of simplicity, the formulation of Theorem 4.2 slightly differs from that in [START_REF] Olver | Asymptotics and special functions[END_REF], the original formulation can be recovered by a change of variable.

Liouville-Green Approximation of (VM-DIN-AVD)

We now proceed to make use of the LG method for approximating the solutions of (Q1-VM-DIN-AVD).

The reader only interested in the result can jump directly to the Section 4.4. We first make the following assumption. Assumption 3. The functions α and ε are three times continuously differentiable, and ε 0 is such that ∀t ≥ 0, ε 0 < (βλ) 2 2|α (t)|+4λ . Remark 4.4. The condition on ε 0 in Assumption 3 is only technical, so that r defined below is positive. It can be easily satisfied since |α (t)| is uniformly bounded. Indeed, α is non-increasing and nonnegative (see Section 1.1), from which one can deduce that

+∞ 0 |α (s)| ds ≤ α 0 .
We now rewrite (Q1-VM-DIN-AVD) in the form [START_REF] Carl | Asymptotic methods and perturbation theory[END_REF]. Lemma 4.5. Suppose that Assumption 3 holds, and let x be the solution of (Q1-VM-DIN-AVD). For all t ≥ 0, define

p(t) = α(t) + βλ ε(t) and r(t) = p(t) 2 4 + p (t) 2 - λ ε(t) . ( 15
)
Then, p and r are twice continuously differentiable, r is positive and the function y defined for all t ≥ 0 by y(t) = x(t) exp

t 0 p(s) 2 ds is a solution to ÿ(t) -r(t)y(t) = 0, t ≥ 0, (16)
with initial condition (y(0), ẏ(0)) = (x 0 , ẋ0 + p(0) 2 x 0 ).

Proof. We first check that for all t ≥ 0, r(t) is positive. Let t > 0,

r(t) > 0 ⇐⇒ (α(t) + βλ) 2 4ε(t) 2 + α (t) 2ε(t) - (α(t) + βλ)ε (t) ε(t) 2 - λ ε(t) > 0. (17)
Since ε (t) ≤ 0 and α (t) ≤ 0, one can check that a sufficient condition for [START_REF] Brillouin | Remarques sur la mécanique ondulatoire[END_REF] to hold is,

r(t) > 0 ⇐= (α(t) + βλ) 2 4 > |α (t)| 2 + λ ε(t) ⇐= (βλ) 2 2|α (t)| + 4λ > ε 0 .
So under Assumption 3, for all t ≥ 0, r(t) > 0. We then check that y is indeed solution to [START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian driven damping[END_REF] Lemma 4.5 gives a reformulation of (Q1-VM-DIN-AVD) suited to apply Theorem 4.2. To use the theorem for all t ≥ 0, we need to ensure that ϕ(t) = 4r(t)r (t)-5r (t) 2 16r(t) 5/2 is integrable. To this aim we make the following assumption. Assumption 4. The functions ε and α have first, second and third-order derivatives that are integrable on [0, +∞[. In addition, lim t→∞ ε(t) = 0 and ε (t) 2 /ε(t) is integrable on [0, +∞[. Remark 4.6. Assumption 4 holds for most decays used in practice, with in particular any polynomial decay of the form ε 0 (t+1) a and α 0 (t+1) b , a ∈ N \ {0} and b ∈ N. Note that ε and α need not be integrable and α can even be constant. The proof of this lemma, relies on relatively simple arguments but involves long computations and is thus postponed to Appendix C. We can now use Theorem 4.2 to obtain an exact form for the solution of (Q1-VM-DIN-AVD) based on the LG approximations. Theorem 4.8. Suppose that Assumptions 3 and 4 hold. There exists A, B ∈ R such that x(0) = x 0 , ẋ(0) = ẋ0 and for all t ≥ 0, the solution of (Q1-VM-DIN-AVD) is

x(t) = A 1 + δ 1 (t) r(t) 1/4 α(t) + βλ ε(t) exp t 0 - λ α(s) + βλ - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds +B 1 + δ 2 (t) r(t) 1/4 ε(t) α(t) + βλ exp t 0 - α(s) + βλ ε(s) + λ α(s) + βλ + λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds , (18)
where for all t ≥ 0,

|δ 1 (t)| ≤ exp 1 2 t 0 |ϕ(s)| ds -1 and |δ 2 (t)| ≤ exp - 1 2 +∞ t |ϕ(s)| ds -1. (19)
Thanks to the bounds [START_REF] Castera | Inertial Newton algorithms avoiding strict saddle points[END_REF], we now have an approximation of x. We will use it in particular to compare x asymptotically to the solutions of (Q1-LM) and (Q1-CN). Before this, we prove Theorem 4.8.

Proof of Theorem 4.8. Let x be the solution of (Q1-VM-DIN-AVD) define p, r as in [START_REF] Blumenthal | Über asymptotische Integration linearer Differentialgleichungen mit Anwendung auf eine asymptotische Theorie der Kugelfunktionen[END_REF]. Let us also define y(t

) def = x(t) exp t 0 p(s)
2 ds . According to Lemma 4.5, r is positive and y is solution to [START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian driven damping[END_REF]. Then, from Lemma 4.7, T t |ϕ(s)| ds < +∞, so we can apply Theorem 4.2 to y on [0, +∞[. Therefore, there exists A, B ∈ R, such that ∀t ≥ 0,

y(t) = A 1 + δ 1 (t) r(t) 1/4 exp t 0 r(s) ds + B 1 + δ 2 (t) r(t) 1/4 exp t 0 -r(s) ds ,
where A and B are determined by the initial conditions, and δ 1 , δ 2 are such that (19) holds.

Going back to x(t) = y(t) exp t 0 -p(s) 2 ds , we obtain that for all t ≥ 0,

x(t) = A 1 + δ 1 (t) r(t) 1/4 exp t 0 - p(s) 2 + r(s) ds + B 1 + δ 2 (t) r(t) 1/4 exp t 0 - p(s) 2 -r(s) ds . (20)
It now remains to expand the terms in the two exponentials in [START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF] in order to obtain [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF]. To this aim, we approximate r(s), let s ≥ 0,

r(s) = p(s) 2 1 + 2p (s) p(s) 2 - 4λ ε(s)p(s) 2 = p(s) 2 1 + p (s) p(s) 2 - 2λ ε(s)p(s) 2 - 1 8 2p (s) p(s) 2 - 4λ ε(s)p(s) 2 2 + o(ε(s) 2 ) = p(s) 2 + p (s) 2p(s) - λ ε(s)p(s) - 1 16
2p (s) p(s) 3/2 - 4λ ε(s)p(s) 3/2 2 + o(ε(s)) = p(s) 2 + p (s)ε(s) 2(α(s) + βλ) - λ α(s) + βλ - 1 16
2p (s) p(s) 3/2 - 4λ ε(s) (α(s) + βλ) 3/2 2 + o(ε(s)) = p(s) 2 + α (s) 2(α(s) + βλ) - ε (s) 2ε(s) - λ α(s) + βλ - 1 16 2p (s) p(s) 3/2 - 4λ ε(s) (α(s) + βλ) 3/2 2 + o(ε(s)) (21)
Focusing on the first exponential term in [START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF], we deduce from ( 21) that for all t ≥ 0,

exp t 0 - p(s) 2 + r(s) ds = exp   t 0 α (s)/2 α(s) + βλ - ε (s) 2ε(s) - λ α(s) + βλ - 1 16 2p (s) p(s) 3/2 - 4λ ε(s) (α(s) + βλ) 3/2 2 + o(ε(s)) ds   = α(t) + βλ) √ α 0 + βλ √ ε 0 ε(t) exp   t 0 -λ α(s) + βλ - 1 16 2p (s) p(s) 3/2 - 4λ ε(s) (α(s) + βλ) 3/2 2 + o(ε(s)) ds   = α(t) + βλ) √ α 0 + βλ √ ε 0 ε(t) exp t 0 -λ α(s) + βλ - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds ,
where the last line relies on further computations postponed to Lemma C.1 in Appendix C. Performing the exact same type of computations on exp t 0 -p(s) 2 -r(s) ds , and up to redefining A and B so as to encompass all the constants, we obtain [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF] and the result is proved.

Comparison of x with x LM and x N

We now have an expression for x which is almost explicit: we do not know δ 1 and δ 2 in closed form, but they are uniformly bounded (by Lemma 4.7). We will now compare the asymptotic behavior of [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF] with those of the solutions of (Q1-LM) and (Q1-CN) that we denoted x LM and x N respectively. Our main result of Section 4 is the following, where ∼ +∞ denotes the asymptotic equivalence 7 between two functions as t → ∞. Theorem 4.9. Let x be the solution of (Q1-VM-DIN-AVD), given in [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF], and x LM and x N whose closed forms are stated in [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF]. Under Assumptions 3 and 4, there exists C > 0 such that the following asymptotic equivalences hold:

x(t) ∼ +∞ x LM (t)C exp t 0 - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds , and x(t) ∼ +∞ x N (t)C exp t 0 α(s) β(α(s) + βλ) - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds . ( 22
)
As a consequence, the convergence of x to x is: (i) Faster than that of x LM if ε is non-integrable and as fast otherwise.

(ii) Slower than that of x N if α is non-integrable and as fast if α is integrable, in the case where ∀t ≥ 0, α(t) > ε(t).

(iii) Faster than that of x N if ε is non-integrable and as fast if ε is integrable, in the case where ∀t ≥ 0, α(t) < ε(t).

While the results of Section 3 were related to the closeness of (VM-DIN-AVD) w.r.t. (CN) and (LM) from a control perspective, Theorem 4.9 provides a different type of insight. First, the results are asymptotic, so they only allow to control (VM-DIN-AVD) for large t. They provide however a clear understanding of the nature of the solutions of (VM-DIN-AVD) and their convergence. The conclusions (summarized in Table 1) are in accordance with what we would expect: when the viscous damping is larger than the variable mass, (VM-DIN-AVD) behaves more like the Levenberg-Marquardt method than the Newton one, but it actually becomes an accelerated Levenberg-Marquardt dynamics when ε is non-integrable but vanishing. However, when the variable mass ε is larger than α, the dynamics is closer to the one of the Newton method, and can actually be an accelerated Newton dynamics, again for non-integrable ε. This is analogous to the necessary condition that α must be non-integrable in order to accelerate first-order methods in convex optimization (see [START_REF] Attouch | Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity[END_REF]). We conclude this section by proving Theorem 4.9.

Proof of Theorem 4.9. Thanks to Assumptions 3 and 4, Theorem 4.8 tells us that x has the form [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF].

We now analyze the two terms in [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF].

First, we know from Theorem 4.8 that δ 1 (0) = 0 and lim t→+∞ δ 2 (t) = 0. In addition, by Lemma 4.7, δ 1 and δ 2 are uniformly bounded by some positive constant. Then r(t) -1/4 decays asymptotically like

ε(t) and α is bounded. So A 1+δ 1 (t) r(t) 1/4 √ α(t)+βλ √ ε(t)
is asymptotically equivalent to some constant c 1 ∈ R as

t → +∞. Similarly, B 1+δ 2 (t) r(t) 1/4 √ ε(t) √ α(t)+βλ is equivalent to c 2 ε(t), with c 2 ∈ R.
We now analyze the "exponential factors" in [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF]. On the one hand,

λ α(s)+βλ + λ 2 ε(s) (α(s)+βλ) 3 + o(ε(s)) converges to 1 β as s → ∞, while α(s)+βλ ε(s)
diverges to +∞. Therefore, we deduce that,

exp t 0 - α(s) + βλ ε(s) + λ α(s) + βλ + λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds = o exp t 0 - λ α(s) + βλ - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds .
As a consequence, the second term in (18) will decrease to 0 faster than the first-one (let alone the additional ε(t) decay that we have just discussed). The asymptotic behavior of x will thus be governed by the first term in [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF].

Let us now focus on the first term in [START_REF] George | The convergence of a class of double-rank minimization algorithms[END_REF]. Observe that exp t 0 -λ α(s)+βλ ds is exactly the decay of x LM in [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF]. Thus, we have proved that there exists C > 0, such that the following asymptotic equivalence holds,

A 1 + δ 1 (t) r(t) 1/4 α(t) + βλ ε(t) exp t 0 - λ α(s) + βλ - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds ∼ +∞ x LM (t)C exp t 0 - λ 2 ε(s) (α(s) + βλ) 3 + o(ε(s)) ds ,
which proves the first part of [START_REF] Emmrich | Discrete versions of Gronwall's lemma and their application to the numerical analysis of parabolic problems[END_REF]. The second equivalence in ( 22) is obtained using the following identity,

t 0 - λ α(s) + βλ ds = t 0 - 1 β + α(s) β(α(s) + βλ) ds = - t β + t 0 α(s) β(α(s) + βλ) ds ( 23
)
and e -t/β is precisely the rate at which x N decreases. So ( 22) holds.

It finally remains to deduce the conclusions of the theorem from [START_REF] Emmrich | Discrete versions of Gronwall's lemma and their application to the numerical analysis of parabolic problems[END_REF].

-Regarding the comparison with x LM , the integral

t 0 -λ 2 ε(s) (α(s)+βλ) 3 + o(ε(s))
ds converges if and only if ε is integrable on [0, +∞[, and diverges to -∞ when ε is not. So x converges to 0 at least as fast as x LM and faster when ε is not integrable.

-As for the comparison with x N , if α(s) > ε(s) ≥ 0 for all s ≥ 0, then the integral

t 0 α(s) β(α(s)+βλ) - λ 2 ε(s) (α(s)+βλ) 3 + o(ε(s)) ds is convergent in +∞ if
and only if α is integrable and diverges to +∞ when α is non-integrable. So when α is integrable, the speed of convergence of x is the same as that of x N . When α is not integrable, the convergence to 0 is slower but still holds. Indeed, for all s ≥ 0

α(s) β(α(s)+ 1 β ) < α(s) βα(s) = 1 β . Thus for all t > 0, -t β + t 0 α(s)
β(α(s)+βλ) ds < 0. -Finally, the comparison with x N in the case ε(s) > α(s) is exactly the same as the comparison with x LM using [START_REF] Fletcher | A new approach to variable metric algorithms[END_REF].

Numerical Experiments

We present two set of experiments that illustrate our main results from Sections 3 and 4. We first detail the general methodology.

Methodology

We compare the solutions of (CN), (LM) and (VM-DIN-AVD) obtained for strongly-convex functions in dimension n = 100. Since closed-form solutions are not available, they are estimated via discretization schemes with small step-sizes γ = 10 -1 . We used Euler semi-explicit schemes, where a linear system is solved at each iteration, for the sake of stability. The resulting algorithms are detailed in Appendix D. Distance to CN (t) = 0/log t, (t) = 0/t (t) = 0/ t, (t) = 0/t (t) = 0/t, (t) = 0/t (t) = 0/t 2 , (t) = 0/t (t) = 0e t/ , (t) = 0/t

(t) = 0/log t, (t) = 0/t 2 (t) = 0/ t, (t) = 0/t 2 (t) = 0/t, (t) = 0/t 2 (t) = 0/t 2 , (t) = 0/t 2 (t) = 0e t/ , (t) = 0/t 2
(x) = e -x 2 + 1 2 Ax 2 .
Left figures: distance x(t) -x N (t) versus time t, each curve corresponds to a different choice of ε; middle figures: distance x(t) -x LM (t) , again for several ε. Right figures: distance to the optimum x for reference, x N and x LM are in dotted and dashed lines, other curves correspond to (VM-DIN-AVD) for several choices of ε. The brown curve is often hidden behind the purple (and sometimes the pink) curve. Top and bottom rows show results respectively for non-integrable and integrable viscous dampings α. The theoretical bounds from Theorem 3.8 are only displayed on Figure 4 below, for the sake of readability.

First Experiment: Distance between Trajectories

We begin with an empirical validation of the results of Section 3 on the distance between x, x LM and x N . Each of Figures 2, 3 and4 (as well as Figures 6 in Appendix D) corresponds to a different strongly-convex function, specified below its corresponding figure. In order to ensure strong convexity, each function contains a quadratic term of the form Ax 2 , where A is symmetric positive definite.

Several observations can be made from the numerical results, but we first note on the right plots of Figures 2 to 4 that x N always converges asymptotically linearly (i.e., exponentially fast). This is also the case for x and x LM in some (but not all) cases. This is important because x(t) -x N (t) ≤

x(t) -x + x N (t) -x , so if both x and x N converge linearly, then the bounds of Theorems 3.1 and 3.8 need not be tight asymptotically. That being said, the strength of these results is to be nonasymptotic and this is highlighted by the experiments as we now explain.

Looking at the left and middle plots of Figures 2, 3 and4, we observe that Theorems 3.1 and 3.8 seem empirically validated, since the distances x(t) -x N (t) and x(t) -x LM (t) decrease relatively fast to zero. Again, when x converges rapidly to x this is not very insightful, however, the main interest Distance to CN (t) = 0/log t, (t) = 0/t (t) = 0/ t, (t) = 0/t (t) = 0/t, (t) = 0/t (t) = 0/t 2 , (t) = 0/t (t) = 0e t/ , (t) = 0/t Distance to CN

(t) = 0/log t, (t) = 0/t 2 (t) = 0/ t, (t) = 0/t 2 (t) = 0/t, (t) = 0/t 2 (t) = 0/t 2 , (t) = 0/t 2 (t) = 0e t/ , (t) = 0/t 2
(x) = log ( n i=1 e x i + e -x i ) + 1 2 Ax 2 .
of our theorems appears on the left of Figures 3 and4: the blue and green curves, corresponding to slowly decaying choices of ε, converge more slowly than other trajectories. However, when taking faster decays, we recover fast convergence and closeness to x N (this is particularly true for the purple curve). Very similar observations are made w.r.t. x LM on the middle plots. Despite not being stated in the theorems of Section 3, the experiments match the intuition that when ε > α, x may be closer to x N and when ε ≤ α, x would rather be closer to x LM . This is more noticeable on the top rows of the figures, where α is not integrable.

Figure 4 suggests that the bounds in Theorem 3.8 are rather tight for small t, since, for example, the blue and green curves on the left show a relatively slow vanishing of x(t)-x N (t) for slowly decaying ε. The experiments show that the bounds seem however often too pessimistic for large t, for which the second part of our study provides better insights (see Section 4 and below). Interestingly, slow decays of ε may sometimes result in faster convergence for x than fast decays (and also faster convergence than x LM ), notably on Figure 2. We also note that ε(t) = ε 0 /t combined either with α(t) = α 0 /t or α(t) = α 0 /t 2 seems to always yield fast convergence on these experiments (and sometimes the fastest of all dynamics).

Second Experiment: Empirical Validation of Theorem 4.9

We now turn our attention to the solutions x, x N and x LM for a quadratic function of the form f (y) = 2 Ax 2 . The thin "dash dotted" curves represent the theoretical bounds from Theorem 3.8 for each choice of (ε, α) considered. as the corresponding 8 x LM when ε is not integrable and regardless of α, and x is faster when ε is nonintegrable. Then on the left-hand side, when comparing to x N , x is slower in settings where α is larger than ε and non-integrable (red curves), or almost as fast when α is integrable (pink curve). However, acceleration w.r.t. to x N is indeed achieved for non-integrable ε regardless of α (first-two blue curves), and the rate is the same as that of x N when ε is integrable (third blue curve).

Conclusions and Perspectives

We introduced a general ODE (VM-DIN-AVD) featuring variable mass, and provided a deep understanding on the behavior of its solutions w.r.t. time dependent control parameters ε and α, both, asymptotically and non-asymptotically. We can conclude that (VM-DIN-AVD) is indeed of (regularized) Newton type, since it can be controlled to be close to both (CN) and (LM). Yet we also showed that (VM-DIN-AVD) fundamentally differs from the other two dynamics in its nature. In particular, Theorem 4.9 and the numerical experiments emphasized that ε and α can accelerate (or slow down) (VM-DIN-AVD) w.r.t. (CN) and (LM). We also note that our bounds in Theorems 3.1 and 3.8 seem relatively tight, in particular for functions with large gradients (see Figure 4). Our contribution yields a complete and satisfying picture on the relation between the three systems, which was only partially understood. We believe that our results build a strong foundation for the development of algorithms that combine the best properties of first-and second-order optimization methods. As for future work, we showed that (VM-DIN-AVD) is promising from an optimization perspective. So far we approximated solutions of (VM-DIN-AVD) via schemes that required solving a linear system at each iteration (this is also true for (CN) and (LM)). Our new understanding on (ε, α) paves the way towards designing new Newton-like algorithms with a significantly reduced computational cost, which is crucial for large-scale optimization. Another open question is whether it is possible to preserve the properties evidenced in this work when ε is defined in a closed-loop manner (formally depending on x rather than on t). Finally, it would be worth investigating how the current work can be extended to general convex and/or non-smooth functions.

A.2 Local Solutions are Global

Using the formulation (gVM-DIN-AVD), we proved local existence and uniqueness of solutions of (VM-DIN-AVD) in Section 2. Using the same notations, we justify that the local solution (x, y) actually exists globally. According to Lemma 3.4, the Lyapunov function

U (t) = ε(t) 2 ẋ(t) 2 + f (x(t)) -f (x )
is non-negative and decreasing. Thus, it is uniformly bounded on R + and the same holds for t → f (x(t)) since for all t ≥ 0, U (t) ≥ f (x(t)). Then, f is coercive by assumption, so x is uniformly bounded on R + (otherwise f (x) could not remain bounded). We now prove that y is also uniformly bounded. From (gVM-DIN-AVD), for all t > 0, ẏ(t) = -

Using the definition of ε and α from Sections 1.1 and 2, observe that ε, α, ε and α are all bounded on R + , and ε is assumed to be bounded. So ν and ν are bounded, and since we also proved that x is uniformly bounded on R + , we deduce from ( 27) that y is uniformly bounded as well. Hence, the unique local solution (x, y) is global.

B Proof of Theorem 3.8

This section is devoted to proving the general result of Section 3. Fix some constants c 1 , c 2 > 0 and let ε and α such that Assumption 2 is satisfied with these constants. Let x be the corresponding solution of (VM-DIN-AVD), x N , and x LM that of (CN) and (LM), respectively. Following the same arguments as in the beginning of the proof of Theorem 3.1, for all t ≥ 0, x(t), x N (t) and x LM (t) belong to the bounded set K 0 defined in that proof. Since f is µ-strongly convex on K 0 , the proof relies again on bounding difference of gradients, indeed, for all t ≥ 0,

x(t) -x N (t) ≤ 1 µ ∇f (x(t)) -∇f (x N (t)) and x(t) -x LM (t) ≤ 1 µ ∇f (x(t)) -∇f (x LM (t)) . (28
) Recall also that the closed form of ∇f (x N ) is given in (4).

Expressing ∇f (x). We follow the exact same steps as in the proof of Theorem 3.1 to obtain the expression of ∇f (x) given in ( 5), which we recall, for all t ≥ 0,

β∇f (x(t)) = βe -t β ∇f (x 0 ) + e -t β ε 0 ẋ0 -ε(t) ẋ(t) + t 0 e s-t β 1 β ε(s) + ε (s) -α(s) ẋ(s) ds.
Here we do not assume any relation between ε and α, and we thus need to find a more suitable expression for ∇f (x(t)). We first expand the terms in the integral, for all t ≥ 0,

β∇f (x(t)) = βe -t β ∇f (x 0 ) + e -t β ε 0 ẋ0 -ε(t) ẋ(t) + t 0 e s-t β 1 β ε(s) + ε (s) ẋ(s) ds - t 0 e s-t β α(s) ẋ(s) ds. (29)
Then, for all s ≥ 0, we have the identity,

which we use to perform an integration by part on the last integral in (29), Therefore,

e -t β t 0 e s β α(s) ẋ(s) ds = α(t)x(t) -e -t β α 0 x 0 - t 0 e s-t β α (s) + α(s) β x(s) ds, (31)
and we can substitute in [START_REF] Lin | A control-theoretic perspective on optimal high-order optimization[END_REF],

β∇f (x(t)) = βe -t β ∇f (x 0 ) + e -t β ε 0 ẋ0 -ε(t) ẋ(t) + t 0 e s-t β 1 β ε(s) + ε (s) ẋ(s) ds -α(t)x(t) + e -t β α 0 x 0 + t 0 e s-t β α (s) + α(s) β x(s) ds. ( 32
)
Uniform boundedness. In view of exploiting [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF], we recall that for all (ε, α), x is uniformly bounded. So there exists R > 0 such that for all (ε, α), the corresponding solution x of (VM-DIN-AVD) is such that sup

t≥0 x(t) ≤ R. (33)
We are now in position to prove [START_REF] Attouch | Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in hilbert spaces[END_REF].

Distance from x to x N . We first gather ( 4) and [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF]. For all t ≥ 0,

β∇f (x(t)) -β∇f (x N (t)) = e -t β ε 0 ẋ0 -ε(t) ẋ(t) + t 0 e s-t β 1 β ε(s) + ε (s) ẋ(s) ds + e -t β α 0 x 0 -α(t)x(t) + t 0 e s-t β α (s) + α(s) β x(s) ds.
We then use [START_REF] Anthony | Wellenmechanik und halbzahlige Quantisierung[END_REF] and triangle inequalities,

βµ x(t) -x N (t) ≤ e -t β ε 0 ẋ0 + ε(t) ẋ(t) + t 0 e s-t β ε(s) β + ε (s) ẋ(s) ds + e -t β α 0 x 0 + α(t) x(t) + t 0 e s-t β α(s) β + α (s) x(s) ds. ( 34
)
By Assumption 2, for all s ≥ 0,

| ε(s) β + ε (s)| ≤ ( 1 β + c 1 )ε(s) and | α(s) β + α (s)| ≤ ( 1 β + c 2 )α(s)
. We then use Lemma 3.5 (denoting by C > 0 the constant stated in the lemma) on the first line of (34), and we use the boundedness (33) on the second line to obtain,

βµ x(t) -x N (t) ≤ e -t β ε 0 ẋ0 + C ε(t) + C 1 β + c 1 t 0 e s-t β ε(s) ds + e -t β α 0 x 0 + Rα(t) + R 1 β + c 2 t 0 e s-t β α(s) ds.
This proves [START_REF] Attouch | Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in hilbert spaces[END_REF].

Expressing ∇f (x LM ). We now repeat previous arguments but for (LM). First, (LM) is equivalent to

d dt ∇f (x LM (t)) + 1 β ∇f (x LM (t)) = -α(t) ẋLM (t).
So using an integrating factor one can check that for all t ≥ 0,

∇f (x LM (t)) = e -t β ∇f (x 0 ) -e -t β t 0 1 β e s β α(s) ẋLM (s) ds.
We can then follow exactly steps [START_REF] Liouville | Mémoire sur le développement des fonctions ou parties de fonctions en séries dont les divers termes sont assujétis à satisfaire à une même équation différentielle du second ordre, contenant un paramètre variable[END_REF] to (31) so as to obtain,

e -t β t 0 e s β α(s) ẋLM (s) ds = α(t)x LM (t) -e -t β α 0 x 0 -e -t β t 0 α (s) + α(s) β e s β x LM (s) ds.
Finally, remark that for all t ≥ 0,

d dt f (x LM (t)) = -α(t) ẋLM (t) 2 -β ẋLM (t), ∇ 2 f (x LM (t)) ẋLM (t) ≤ 0.
So f (x LM (t)) ≤ f (x 0 ) and using the coercivity of f as before we deduce that for all choices α, sup t≥0

x LM (t) ≤ R.

Distance from x to x LM . We substract gradients,

β∇f (x(t)) -β∇f (x LM (t))) = e -t β ε 0 ẋ0 -ε(t) ẋ(t) + t 0 e s-t β 1 β ε(s) + ε (s) ẋ(s) ds -α(t)(x(t) -x LM (t)) - t 0 e s-t β α (s) + α(s) β (x(s) -x LM (s)) ds,
and we proceed as before using (28), Assumption 2 and Lemma 3.5. It holds that,

βµ x(t) -x LM (t) ≤ e -t β ε 0 ẋ0 + C ε(t) + C 1 β + c 1 t 0 e 1 β (s-t) ε(s) ds + α(t) x(t) -x LM (t) + 1 β + c 2 t 0 e 1 β (s-t) x(s) -x LM (s) ds.
Finally, using ( 33) and ( 35), for all s ≥ 0, x(s) -x LM (s) ≤ 2R, which concludes the proof.

On the one hand, it holds that, On the other hand, we have,

r 2 ε 5 =   - (α + βλ) 2 -4λε (α+βλ) 2 + 1 + (-2(α+βλ)ε +2α ε) (α+βλ) 2 ε 2 √ ε + (α + βλ) 2 √ ε 4 - 4λε (α + βλ) 2 + 8λα ε (α + βλ) 3
r ε(t) 3 = -λε ε + 4λα ε ε α + βλ + 2λα ε 2 α + βλ - 6λα 2 ε 2 (α + βλ) 2 + (α + βλ) 2 2 3ε 2 ε + ε 4λε (α + βλ) 2 -1 + 2 ((α + βλ)ε -α ε) (α + βλ
(α + βλ)ε ε -6α ε 2 + 3α ε ε + 3α ε ε -α ε 2 + 1 α + βλ 4 (α + βλ)ε -α ε α ε + (α + βλ)ε -α ε α ε + 2 -2(α + βλ)ε 2 + (α + βλ)ε ε + 2α ε ε -α ε 2 α - 2 α + βλ 2λε - 4λα ε α + βλ + 2 (α + βλ)ε ε -α ε + - 2(α + βλ)ε 2 ε + (α + βλ)ε + 2α ε -α ε - 2 ((α + βλ)ε -α ε) α α + βλ α ε - 3 (α + βλ)ε ε -α ε 2 α 2 (α + βλ) 2 .
We then analyze the integrability of each of the terms above. By Assumption 4, ε , ε and ε are integrable and the same goes for α , α and α , which is enough to justify the integrability of almost all the terms above. We finally see that we also need (ε ) 2 ε and (ε ) 3 ε to be integrable, which holds by Assumption 4. Overall, ϕ is integrable on R + .

We now state and prove the following result which was used at the end of the proof of Theorem 4.8. Proof. We omit the time dependence on s ≥ 0 for the sake of readability. Using Assumption 3 we can define and expand the following quantity, Assumption 4, implies in particular that |ε (t)| = o(ε(t)) and that α (t) → 0, which we use in the equality above to obtain the desired conclusion.

D Additional Experiments and Details

We first detail the discretization that we used for approximating the solutions of the three ODEs considered in Section 5. We use Euler discretization schemes with fixed step-size γ > 0 and approximate the solutions at times t k = γk, for all k ∈ N. For a trajectory x, we use the notation x(t k ) def = x (k) . The approximation of (CN) is obtained by explicit discretization, so that for all k ∈ N, we have,

x (k+1) N = x (k) N -γ β∇ 2 f (x (k) N ) -1 ∇f (x (k) N ). (38)
Then, defining ε k = ε(t k ) and α k = α(t k ), (LM) and (VM-DIN-AVD) are obtained via Euler semiimplicit discretization. The solution of (LM) is approximated by,

x (k+1) LM = x (k) LM -γ α k I n + β∇ 2 f (x (k) LM ) -1 ∇f (x (k) LM ), (39)
where I n is the identity matrix on R n . The solution of (VM-DIN-AVD) is obtained similarly,

x (k+1) = x (k) + (ε k + γα k )I n + γβ∇ 2 f (x (k) ) -1 ε k (x (k) -x (k-1) ) -γ 2 ∇f (x (k) ) . (40)
As safety check, one can see that for ε k = 0, (40) is equivalent to [START_REF] David | Conditioning of quasi-Newton methods for function minimization[END_REF], which is itself equivalent to [START_REF] Rossum | Python reference manual[END_REF] when α k = 0. Distance to CN (t) = 0/log t, (t) = 0/t (t) = 0/ t, (t) = 0/t (t) = 0/t, (t) = 0/t (t) = 0/t 2 , (t) = 0/t (t) = 0e t/ , (t) = 0/t Distance to CN (t) = 0/log t, (t) = 0/t 2 (t) = 0/ t, (t) = 0/t 2 (t) = 0/t, (t) = 0/t 2 (t) = 0/t 2 , (t) = 0/t 2 (t) = 0e t/ , (t) = 0/t 2 Distance to CN (t) = 0/log t, (t) = 0/t (t) = 0/ t, (t) = 0/t (t) = 0/t, (t) = 0/t (t) = 0/t 2 , (t) = 0/t (t) = 0e t/ , (t) = 0/t Distance to CN (t) = 0/log t, (t) = 0/t 2 (t) = 0/ t, (t) = 0/t 2 (t) = 0/t, (t) = 0/t 2 (t) = 0/t 2 , (t) = 0/t 2 (t) = 0e t/ , (t) = 0/t 2

Figure 1 :

 1 Figure 1: Left: phase diagram on distances from (VM-DIN-AVD) to (CN) and (LM) (see Section 3).For each patch, the color indicates which of the distances x N (t) -x(t) and x LM (t) -x(t) is considered, the scaling of a corresponding upper-bound on this distance is written; in white for prior work and in black for our contributions. The green line separates the cases ε ≥ α (above) and ε ≤ α (below). Right: 2D illustration of the trajectories of (VM-DIN-AVD) for several choices of ε on a quadratic function. Using fast-vanishing ε(t) (dark-blue solid curves), one can bring the solution of (VM-DIN-AVD) close to that of (CN), making it, for example, more robust to bad conditioning compared to first-order dynamics (such as gradient descent).

  The next lemma states the integrability of ϕ on [0, +∞[. Lemma 4.7. Under Assumption 3 and 4, +∞ 0 |ϕ(s)| ds < +∞.

Figure 2 :

 2 Figure2: Comparison of the solutions x N , x LM and x of (CN), (LM) and (VM-DIN-AVD) respectively, for a strongly-convex function of the form f (x) = e -x 2 + 1 2 Ax 2 . Left figures: distance x(t) -x N (t) versus time t, each curve corresponds to a different choice of ε; middle figures: distance x(t) -x LM (t) , again for several ε. Right figures: distance to the optimum x for reference, x N and x LM are in dotted and dashed lines, other curves correspond to (VM-DIN-AVD) for several choices of ε. The brown curve is often hidden behind the purple (and sometimes the pink) curve. Top and bottom rows show results respectively for non-integrable and integrable viscous dampings α. The theoretical bounds from Theorem 3.8 are only displayed on Figure4below, for the sake of readability.

Figure 3 :

 3 Figure 3: Similar experiment and figures as those described in Figure 2, but for the function f(x) = log ( n i=1 e x i + e -x i ) + 1 2 Ax 2 .

Figure 4 :

 4 Figure4: Similar experiment and figures as those described in Figure2, but for the function f (x) = x 50 + 1 2 Ax 2 . The thin "dash dotted" curves represent the theoretical bounds from Theorem 3.8 for each choice of (ε, α) considered.

Figure 5 : 1 2

 51 Figure 5: Numerical validation of Theorem 4.9: distance to the optimum x as a function of time. on a quadratic function f (x) = 1 2 Ax 2 . Left: speed comparison w.r.t. (CN) for several choices of ε and α. Right: Comparison with LM in for α integrable or not and several choices of ε. Shades of blue represent cases where ε(t) > α(t) while shades of red represent the opposite setting.

- 2 (√ ε 2 .

 22 α + βλ)ε -4α ε + 2α ε (α + βλ) 2 -2 (-2(α + βλ)ε + 2α ε) α (α + βλ) βλ) 2 + 1 + (-2(α + βλ)ε + 2α ε) (α + βλ) 2 α

) 2 + 2 ( 2 6

 222 α + βλ) 4λε (α + βλ) 2 -1 + 2 ((α + βλ)ε -α ε) (α + βλ) 2 α ε -(α + βλ)α + α 2 2 4λε (α + βλ) 2 -1 + 2 ((α + βλ)ε -α ε) (α + βλ) 2 ε -(α + βλ)ε ε -α ε 2 -(α + βλ)ε -α ε ε -2 -2(α + βλ)ε 2 ε + (α + βλ)ε + 2α ε -α ε ε + 2 2λε -4λα ε α + βλ + 2 (α + βλ)ε ε -α ε + -2(α + βλ)ε 2 ε + (α + βλ)ε + 2α ε -α ε -2 ((α + βλ)ε -α ε) α α + βλ ε -1 (α + βλ)ε 3 ε -6(α + βλ)ε ε +

Lemma C. 1 . 2 = λ 2

 122 Under Assumptions 3 and 4, for all s ≥ 0ε(s) (α(s) + βλ)3 + o(ε(s)).

3 = λ 2 ε

 32 (α + βλ) 3 -p λ √ ε p 3/2 (α + βλ) 3/2 + (p ) 2 4p (α + βλ) 3 -λ(α ε -(α + βλ)ε ) + 1 4(α + βλ) 3 (α ) 2 ε + (ε ) 2 ε (α + βλ) 2 -2α ε (α + βλ) .

Figure 6 :

 6 Figure6: Similar experiment and figures as those described in Figure2, but for a poorly conditioned quadratic f (x) =1 2 Ax 2 (first-two rows) and the function f (x) = log ( n i=1 e x i ) + 1 2 Ax 2 (last-two rows).

Table 1 :

 1 Informal summary of Section 4. Comparison of (VM-DIN-AVD) with other dynamics

	Parameters of (VM-DIN-AVD)	Speed of convergence
	Dominant parameter Integrability in +∞	w.r.t. (CN) w.r.t. (LM)
	variable mass ε	yes no	as fast faster	as fast faster
	viscous damping α	yes no	as fast slower	only depends on ε

  1/4 exp Remark 4.3. We make the following remarks regarding the above result. -Note that z 1 and z 2 in Theorem 4.2 are exact solutions to[START_REF] Carl | Asymptotic methods and perturbation theory[END_REF]. The LG approximations ẑ1 and ẑ2 are obtained by neglecting the unknown functions δ 1 and δ 2 in z 1 and z 2 . The theorem gives a nonasymptotic bound for the errors |z 1 (t) -ẑ1 (t)| and |z 2 (t) -ẑ2 (t)|, t ≥ 0.

		0	t	r(s) ds	and z 2 (t) =	1 + δ 2 (t) r(t) 1/4 exp -	0	t	r(s) ds ,
	where |δ 1 (t)| ≤ exp	1 2

t 0 |ϕ(s)| ds -1 and |δ 2 (t)| ≤ exp -1 2 T t |ϕ(s)| ds -1. If in addition +∞ 0 |ϕ(s)| ds < +∞, then the results above also hold for T = +∞.

1 2

 1 Ay 2 , y ∈ R n , and for several choices of ε and α. The results in Figure5exactly match the expected behavior summarized in Table1. Indeed, looking first at the right-hand side of Figure5, x is as fast

	10 4 10 3			Distance to CN				Distance to LM			Speed of convergence Levenberg-Marquardt Continuous Newton
	10 11
	10 18
	10 25
	10 46 10 39 10 32	(t) = 0/log t, (t) = 0/t (t) = 0/ t, (t) = 0/t (t) = 0/t, (t) = 0/t (t) = 0/t 2 , (t) = 0/t (t) = 0e t/ , (t) = 0/t
	0	20	40	60	t	80 100 120 140	0	20	40	60	t	80 100 120 140	0	20	40	60	t	80 100 120 140
	10 4			Distance to CN				Distance to LM			Speed of convergence
	10 5
	10 14
	10 23
	10 32
	10 59 10 50 10 41	(t) = 0/log t, (t) = 0/t 2 (t) = 0/ t, (t) = 0/t 2 (t) = 0/t, (t) = 0/t 2 (t) = 0/t 2 , (t) = 0/t 2 (t) = 0e t/ , (t) = 0/t 2								Levenberg-Marquardt Continuous Newton
	0	20	40	60	t	80 100 120 140	0	20	40	60	t	80 100 120 140	0	20	40	60	t	80 100 120 140

CN is usually considered with β = 1, we put β in the system to ease the discussions below.

DIN-like systems were thought to yield faster vanishing of the gradient compared to inertial first-order dynamics, until recently[START_REF] Attouch | From the ravine method to the Nesterov method and vice versa: A dynamical system perspective[END_REF].

Assumption 1 can actually hold only after some large-enough t 0 ≥ 0, we take t 0 = 0 for the sake of simplicity.

This is a consequence of Gronwall's lemma, see e.g.,[START_REF] Emmrich | Discrete versions of Gronwall's lemma and their application to the numerical analysis of parabolic problems[END_REF].

This can be generalized to the case where A T A is only semi-definite by considering orthogonal projections on an eigenspace spanned by the positive eigenvalues of A T A.

We express ψ(τ ) using t for the sake of readability, using the one-to-one correspondence between τ and t.

Two real-valued functions g 1 and g 2 are asymptotically equivalent in +∞ if and only if lim t→∞ g1(t) g2(t) = 1.

That is, the solution of (LM) for the same α as that considered for (VM-DIN-AVD).

Acknowledgment

C. Castera, J. Fadili and P. Ochs are supported by the ANR-DFG joint project TRINOM-DS under number ANR-20-CE92-0037-01. The numerical experiments were made thanks to the development teams of the following libraries: Python [38], Numpy [43] and Matplotlib [27].

Appendices

A Equivalent First-order System and Global Existence of Solutions A.1 First-order Equivalent Formulation

We reformulate (VM-DIN-AVD) as a system of ODE involving only first-order time derivatives and the gradient of f . For this purpose, notice that for all t > 0 (VM-DIN-AVD) can be rewritten as,

We then integrate [START_REF] Konstantinovich | Nonlinear functional equations and continuous analogues of iteration methods[END_REF] for all t ≥ 0,

For all t ≥ 0, we define the variable,

We differentiate z, for all t > 0, ż(t) = (α(t) -ε (t)) ẋ(t) + ∇f (x(t)), so that we can rewrite (25) as,

We substitute the first line in the second-one,

To ease the readability, we recall the notation

which is (gVM-DIN-AVD). Finally, the initial condition on y is

Remark A.1. Notice that the quantity ν(t) = α(t) -ε (t) -1 β ε(t) involved in (gVM-DIN-AVD) also plays a key role in our analysis of Section 3, see e.g., [START_REF] Attouch | Continuous Newton-like inertial dynamics for monotone inclusions[END_REF]. In particular the sign of ν(t) changes the nature of (VM-DIN-AVD) and is related to Assumption 1.

C Integrability of ϕ and Additional Asymptotic Computations

Below we prove Lemma 4.7.

Proof. We suppose that Assumptions 3 and 4 hold. As stated in Remark 4.3, since ϕ is continuous, we only need to check its integrability when t tends to +∞. Let t > 0, we first establish some useful identities, we omit the dependence on t for the sake of readability.

Then,

Then, to justify that ϕ is integrable, we prove that r r 3/2 and (r ) 2 r 5/2 are integrable. Since we know that 1/r decays at the same speed as ε 2 , we can equivalently show that ε 3 r and ε 5 (r ) 2 are integrable. To this aim we fully expand all the terms in [START_REF] Boris | Some methods of speeding up the convergence of iteration methods[END_REF], and compute (r ) 2 , which is extremely long and involved.