# Continuous Newton-like Methods featuring Inertia and Variable Mass

Camille Castera, Hedy Attouch, Jalal Fadili, Peter Ochs

<figure><img src="image_1.png" /></figure>

## To cite this version:

Camille Castera, Hedy Attouch, Jalal Fadili, Peter Ochs. Continuous Newton-like Methods featuring Inertia and Variable Mass. 2023. hal-03951620

HAL Id: hal-03951620

https://hal.science/hal-03951620

Preprint submitted on 23 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Continuous Newton-like Methods featuring Inertia and Variable Mass

Camille Castera, Hedy Attouch, Jalal Fadili, Peter Ochs

► To cite this version:

Camille Castera, Hedy Attouch, Jalal Fadili, Peter Ochs. Continuous Newton-like Methods featuring Inertia and Variable Mass. SIAM Journal on Optimization, 2024, 34 (1), pp.251-277. 10.1137/23M1549675. hal-03951620v2

HAL Id: hal-03951620

https://hal.science/hal-03951620v2

Submitted on 19 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Continuous Newton-like Methods featuring Inertia and Variable Mass

**Camille Castera***

Department of Mathematics
University of Tübingen
Germany

**Hedy Attouch**

IMAG

Université Montpellier
CNRS, France

**Peter Ochs**

Department of Mathematics
University of Tübingen
Germany

**Jalal Fadili**

ENSICAEN

Normandie Université
CNRS, GREYC, France

## ABSTRACT

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.

# 1 Introduction

## 1.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 \in \mathbb{R}^n} f(x). $$

<footnote>*Corresponding author: camille.castera@protonmail.com</footnote>
# Continuous Newton-like Methods featuring Inertia and Variable Mass

**Camille Castera***

Department of Mathematics
University of Tübingen
Germany

**Hedy Attouch**

IMAG
Université Montpellier
CNRS, France

**Jalal Fadili**

ENSICAEN
Normandie Université
CNRS, GREYC, France

**Peter Ochs**

Departments of Computer Science and Mathematics
Saarland University
Germany

*Dedicated to the memory of Hedy Attouch, outstanding mathematician and beloved collaborator.*

## ABSTRACT

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.

# 1 Introduction

## 1.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 \in \mathbb{R}^n} f(x). $$

<footnote>*Corresponding author: camille.castera@protonmail.com
Published in SIAM Journal on Optimization 34(1):251–277</footnote>
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 [37, 33, 34]. 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 trade-off. 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 [18, 23, 25, 39, 31], 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 [32, 13]. 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:

$$\varepsilon(t)\ddot{x}(t) + \alpha(t)\dot{x}(t) + \beta\nabla^2 f(x(t))\dot{x}(t) + \nabla f(x(t)) = 0, \quad t \ge 0, \qquad (\text{VM-DIN-AVD})$$

where $f: \mathbb{R}^n \to \mathbb{R}$ is a smooth convex twice continuously differentiable function, with gradient $\nabla f$ and Hessian $\nabla^2 f$ defined on $\mathbb{R}^n$ equipped with scalar product $\langle \cdot, \cdot \rangle$, and induced norm $\|\cdot\|$. Additionally, $f$ is assumed to be coercive, and strongly convex on bounded subsets of $\mathbb{R}^p$. The functions $\varepsilon, \alpha: \mathbb{R}_+ \to \mathbb{R}_+$ (where $\mathbb{R}_+ = [0, +\infty]$) are differentiable, non-increasing, and $\varepsilon(t) > 0$ for all $t \ge 0$. Together with $\beta > 0$, they are control parameters that define the type of dynamics that drives the trajectory (or solution) $x: \mathbb{R}_+ \to \mathbb{R}^n$, whose first- and second-order derivatives are denoted $\dot{x}$ and $\ddot{x}$ 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 [2], where $\varepsilon$ and $\alpha$ were constant. DIN was then extended to the case of non-constant asymptotically vanishing dampings (AVD) $\alpha$ [9]. In this work we introduce the non-constant parameter $\varepsilon$ called variable mass (VM) in front of the acceleration $\ddot{x}$, in the same way that $\alpha$ 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 $\nabla f$ but not $\nabla^2 f$, significantly reducing computational costs, hence enabling the design of practical algorithms, see e.g., [20, 10, 21]. The key idea behind this is the relation $\nabla^2 f(x(t))\dot{x}(t) = \frac{d}{dt}\nabla 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 $\varepsilon = \alpha = 0$, one obtains<fnref n="1" /> the Continuous Newton (CN) method [24]

$$\beta\nabla^2 f(x_N(t))\dot{x}_N(t) + \nabla f(x_N(t)) = 0, \quad t \ge 0, \qquad (\text{CN})$$

<footnote n="1">¹CN is usually considered with $\beta = 1$, we put $\beta$ in the system to ease the discussions below.</footnote>
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 [43, 38, 39]. 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 trade-off. It is therefore essential to develop algorithms that take the best of both worlds. Several quasi-Newton algorithms partly address this issue, for example BFGS methods [22, 28, 30, 45, 36], 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 [37, 16]. 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:

$$ \varepsilon(t)\ddot{x}(t) + \alpha(t)\dot{x}(t) + \beta\nabla^2 f(x(t))\dot{x}(t) + \nabla f(x(t)) = 0, \quad t \ge 0, \qquad (\text{VM-DIN-AVD}) $$

where $f: \mathbb{R}^n \to \mathbb{R}$ is a smooth function, with gradient $\nabla f$ and Hessian $\nabla^2 f$ defined on $\mathbb{R}^n$ equipped with scalar product $\langle \cdot, \cdot \rangle$, and induced norm $\| \cdot \|$. The system is controlled by two functions $\varepsilon, \alpha: \mathbb{R}_+ \to \mathbb{R}_+$ (where $\mathbb{R}_+ = [0, +\infty[$) and a parameter $\beta > 0$ that define the type of dynamics that drives the trajectory (or solution) $x: \mathbb{R}_+ \to \mathbb{R}^n$, whose first- and second-order derivatives are denoted $\dot{x}$ and $\ddot{x}$ 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 [3], where $\varepsilon$ and $\alpha$ were constant. DIN was then extended to the case of non-constant asymptotically vanishing dampings (AVD) $\alpha$ [11]. In this work we introduce the non-constant parameter $\varepsilon$ called variable mass (VM) in front of the acceleration $\ddot{x}$, in the same way that $\alpha$ 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 $\nabla f$ but not $\nabla^2 f$, significantly reducing computational costs, hence enabling the design of practical algorithms, see e.g., [24, 12, 25]. The key idea behind this is the relation $\nabla^2 f(x(t))\dot{x}(t) = \frac{d}{dt}\nabla 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 special cases. Indeed, taking $\varepsilon = \alpha = 0$, one obtains<fnref n="1" /> the Continuous Newton (CN) method [29]

$$ \beta\nabla^2 f(x_N(t))\dot{x}_N(t) + \nabla f(x_N(t)) = 0, \quad t \ge 0, \qquad (\text{CN}) $$

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),

<footnote n="1">¹CN is usually considered with $\beta = 1$, we put $\beta$ in the system to ease the discussions below.</footnote>
<figure><img src="image_3.png" /><figcaption>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 $\varepsilon \ge \alpha$ (above) and $\varepsilon \le \alpha$ (below). Right: 2D illustration of the trajectories of (VM-DIN-AVD) for several choices of $\varepsilon$ on a quadratic function. Using fast-vanishing $\varepsilon(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).</figcaption></figure>

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 $\alpha \ne 0$ but $\varepsilon = 0$, we recover the Levenberg-Marquardt (LM) method,

$$ \alpha(t)\dot{x}_{LM}(t) + \beta\nabla^2 f(x_{LM}(t))\dot{x}_{LM}(t) + \nabla f(x_{LM}(t)) = 0, \quad t \ge 0, \qquad (\text{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. [2] showed that for $\alpha = 0, \beta = 1$, and $\varepsilon$ 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 $\sqrt{\varepsilon}$ at all time. Yet, despite the benefits of this class of ODEs, such as stabilization properties, see e.g., [9, 10], no improvement<fnref n="2" /> of the rate of convergence (in values) has been shown compared to inertial first-order dynamics [37, 41]. This raises the question:

“are these ODEs really of Newton type?”,

which is crucial in view of designing faster algorithms from them.

<footnote n="2">²DIN-like systems were thought to yield faster vanishing of the gradient compared to inertial first-order dynamics, until recently [4].</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Left: phase diagram on distances from (VM-DIN-AVD) to (CN) and (LM) (see Section 3). The color of each patch indicates which distance is considered, and the scaling of a corresponding upper-bound on this distance is written (in white for prior work, in black for our contributions). The green line separates the cases $\varepsilon \ge \alpha$ (above) and $\varepsilon \le \alpha$ (below). Right: 2D illustration of the trajectories of (VM-DIN-AVD) for several choices of $\varepsilon$ on a quadratic function. Fast-vanishing $\varepsilon(t)$ (dark-blue solid curves) bring solutions of (VM-DIN-AVD) close to that of (CN), making them, more robust to bad conditioning compared to first-order dynamics (e.g., gradient descent).</figcaption></figure>

which also justifies the terminology “Newton-like” in DIN. When $\alpha \neq 0$ but $\varepsilon = 0$, we recover the Levenberg-Marquardt (LM) method,

$$ \alpha(t)\dot{x}_{LM}(t) + \beta\nabla^2 f(x_{LM}(t))\dot{x}_{LM}(t) + \nabla f(x_{LM}(t)) = 0, \quad t \ge 0, \qquad (\text{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. Since the introduction of DIN, it is known (see [3]) that for $\alpha = 0, \beta = 1$, and $\varepsilon$ 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 $\sqrt{\varepsilon}$ at all times. Yet, despite the benefits of this class of ODEs, such as stabilization properties [11, 12], no improvement<fnref n="2" /> of the rate of convergence (in values) has been shown compared to inertial first-order dynamics [43, 47]. This raises the question:

“are these ODEs really of Newton type?”,

which is crucial in view of designing faster algorithms from them.

## 1.2 Main Contributions

We show that the answer to this question is partially positive, and closely related to the choices of $\varepsilon$ and $\alpha$. We provide general results on the role played by these two control parameters and how they

<footnote n="2">²DIN-like systems were thought to yield faster vanishing of the gradient compared to inertial first-order dynamics, until recently [5].</footnote>
<table><caption>Table 1: Informal summary of Section 4. Comparison of (VM-DIN-AVD) with other dynamics</caption><thead><tr><th colspan="2">Parameters of (VM-DIN-AVD)</th><th colspan="2">Speed of convergence</th></tr><tr><th>Dominant parameter</th><th>Integrability in +∞</th><th>w.r.t. (CN)</th><th>w.r.t. (LM)</th></tr></thead><tbody><tr><td rowspan="2">variable mass ε</td><td>yes</td><td>as fast</td><td>as fast</td></tr><tr><td>no</td><td><span style="color:green;">faster</span></td><td><span style="color:green;">faster</span></td></tr><tr><td rowspan="2">viscous damping α</td><td>yes</td><td>as fast</td><td>only depends</td></tr><tr><td>no</td><td><span style="color:red;">slower</span></td><td>on ε</td></tr></tbody></table>

## 1.2 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(\sqrt{ε(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.

## 1.3 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 [2] and called DIN (for Dynamical Inertial Newton-like system). Except in a few cases [20, 19], 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 [9]. The reason for this popular choice for $α$ is its link with Nesterov’s method [41]. Non-constant choices for $β$ have been considered [10, 1, 29, 11]. 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 [2], but only for fixed $ε$, $β = 1$ and constant $α = 0$. VM-DIN-AVD is however closely related to the IGS system considered in [11] 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 [11], which is of independent interest. The literature on DIN is rich, let us mention further connections with Nesterov’s method [40, 1], extensions with Tikhonov
<table><caption>Table 1: Informal summary of Section 4. Comparison of (VM-DIN-AVD) with other dynamics</caption><thead><tr><th colspan="2">Parameters of (VM-DIN-AVD)</th><th colspan="2">Speed of convergence</th></tr><tr><th>Dominant parameter</th><th>Integrability in +∞</th><th>w.r.t. (CN)</th><th>w.r.t. (LM)</th></tr></thead><tbody><tr><td rowspan="2">variable mass ε</td><td>yes</td><td>as fast</td><td>as fast</td></tr><tr><td>no</td><td><span style="color:green;">faster</span></td><td><span style="color:green;">faster</span></td></tr><tr><td rowspan="2">viscous damping α</td><td>yes</td><td>as fast</td><td>only depends</td></tr><tr><td>no</td><td><span style="color:red;">slower</span></td><td>on ε</td></tr></tbody></table>

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 show the 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(\sqrt{ε(t)})$. We then extend these results to larger $α$ 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.

## 1.3 Related work

The importance and challenges of finding systems that approximate and preserve the benefits of (CN) were highlighted by [2]. 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 [3] and called DIN (for Dynamical Inertial Newton-like system). Except in a few cases [24, 23], 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 [11]. The reason for this popular choice for $α$ is its link with Nesterov’s method [47]. Non-constant choices for $β$ have been considered [12, 1, 34, 13]. We keep it constant, and rather focus on non-constant $ε$, unlike prior work that used constant $ε = 1$. The mass $ε$ was only considered in the original work [3], but only for fixed $ε, β = 1$ and constant $α = 0$. VM-DIN-AVD is however related to the IGS system considered in [13] 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 [13], which is of independent interest. The literature on DIN is rich, let us mention further connections with Nesterov’s method [46, 1], extensions with Tikhonov regularization [19] and closed-loop damping [14, 34]. The non-smooth and possibly non-convex cases have been considered in [6, 7, 24]. Finally, avoidance of strict saddle points in smooth non-convex optimization has been shown in [23].
regularization [16] and closed-loop damping [12, 29]. The non-smooth and possibly non-convex cases have been considered in [5, 6, 20]. Finally, avoidance of strict saddle points in smooth non-convex optimization has been shown in [19].

The influence of the damping $\alpha$ on the (LM) dynamics has been studied in [7, 8]. Interestingly, the conditions enforced on $\alpha$ in these papers (formally a sub-exponential decay) are very similar to those we make on $\varepsilon$ and $\alpha$ 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. [10] provided closed-form solutions to (VM-DIN-AVD) for $\varepsilon \equiv 1$ and special choices of $\alpha$. Our work rather deals with approximate solutions which allows considering a wide class of functions $\varepsilon$ and $\alpha$. We rely on the Liouville-Green (LG) method [30, 26] presented in Section 4. Generalizations of LG are also often referred to as WKB methods [44, 28, 17] 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.

## 1.4 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 $\varepsilon$ and $\alpha$ is then carried out on quadratic functions in Section 4. Finally, numerical experiments are presented in Section 5, and some conclusions are then drawn.

# 2 Existence and Uniqueness of Solutions

In the sequel, we fix $x_0 \in \mathbb{R}^n$ and $\dot{x}_0 \in \mathbb{R}^n$, such that, unless stated otherwise, (VM-DIN-AVD) is always considered with initial condition $(x(0), \dot{x}(0)) = (x_0, \dot{x}_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 $\varepsilon(0) = \varepsilon_0 > 0$, $\varepsilon'(0) = \varepsilon'_0 \le 0$ and $\alpha(0) = \alpha_0 \ge 0$. In addition to the definitions of $\varepsilon$ and $\alpha$ in Section 1.1, we assume that $\varepsilon$ 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: \mathbb{R}_+ \to \mathbb{R}^n$. Notably, our reformulation does not involve $\nabla^2 f$, in the same fashion as [2, 9]. For all $t$, defining $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$, we show in Appendix A that (VM-DIN-AVD) is equivalent to

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + \nu(t)x(t) + y(t) = 0 \\ \dot{y}(t) + \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t) = 0 \end{cases} \quad (gVM-DIN-AVD) $$

with initial conditions $(x(0), y(0)) = (x_0, -\varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) - (\alpha_0 - \varepsilon'_0 - \frac{1}{\beta}\varepsilon_0)x_0)$. One can notice that in the special case where $\varepsilon$ is taken constant and equal to 1 (that is when (VM-DIN-AVD) is simply the DIN-AVD system [9]), we recover the same first-order formulation as that in [9].

We can then apply the Cauchy-Lipschitz Theorem. For all $t \ge 0$ and $(u, v) \in \mathbb{R}^n \times \mathbb{R}^n$, define the mapping

$$ G(t, (u, v)) = \begin{pmatrix} \frac{1}{\varepsilon(t)}(-\beta\nabla f(u) - \nu(t)u - v) \\ -\nu'(t)u - \frac{\nu(t)}{\beta}u - \frac{1}{\beta}v \end{pmatrix}, $$
The influence of the damping $\alpha$ on the (LM) dynamics has been studied in [8, 9]. Interestingly, the conditions enforced on $\alpha$ in these papers (formally a sub-exponential decay) are very similar to those we make on $\varepsilon$ and $\alpha$ for (VM-DIN-AVD) (see Assumptions 2 and 3). The (LM) dynamics is also related to the system in [10].

Regarding the second part of our analysis, which deals with the case where $f$ is quadratic, a recent work [12] provided closed-form solutions to (VM-DIN-AVD) for $\varepsilon \equiv 1$ and special choices of $\alpha$. Our work rather deals with approximate solutions which allows considering a wide class of functions $\varepsilon$ and $\alpha$. We rely on the Liouville-Green (LG) method [35, 31] presented in Section 4. Generalizations of LG are also often referred to as WKB methods [50, 33, 21] 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.

## 1.4 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 $\varepsilon$ and $\alpha$ is then carried out on quadratic functions in Section 4. Finally, numerical experiments are presented in Section 5, and some conclusions are then drawn.

# 2 Existence and Uniqueness of Solutions

We always assume the following on the system (VM-DIN-AVD) and the functions $f, \varepsilon, \alpha$ defined in the introduction.

**Assumption 1.** $\bullet$ $f: \mathbb{R}^n \to \mathbb{R}$ is convex, twice continuously differentiable, coercive, and strongly convex on bounded subsets of $\mathbb{R}^n$. We fix $x_0, \dot{x}_0 \in \mathbb{R}^n$, such that, unless stated otherwise, (VM-DIN-AVD) has initial condition $(x(0), \dot{x}(0)) = (x_0, \dot{x}_0)$, and (CN) and (LM) have initial conditions $x_N(0) = x_{LM}(0) = x_0$.

$\bullet$ $\alpha: \mathbb{R}_+ \to \mathbb{R}$ is non-increasing, non-negative and differentiable. $\varepsilon: \mathbb{R}_+ \to \mathbb{R}$ is non-increasing, positive, and twice differentiable with bounded second derivative. We fix initial values: $\varepsilon(0) = \varepsilon_0 > 0$, $\varepsilon'(0) = \varepsilon'_0 \le 0$ and $\alpha(0) = \alpha_0 \ge 0$.

**Theorem 2.1.** Under Assumption 1 there exists a unique global solution $x: \mathbb{R}_+ \to \mathbb{R}^n$ to (VM-DIN-AVD).

The proof relies on the Cauchy-Lipschitz Theorem, we sketch the main elements. We reformulate (VM-DIN-AVD) into a first-order (in time) system by introducing an auxiliary variable $y: \mathbb{R}_+ \to \mathbb{R}^n$. Notably, our reformulation does not involve $\nabla^2 f$, in the same fashion as [3, 11]. For all $t$, defining $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$, we show in Appendix A that (VM-DIN-AVD) is equivalent to

$$ \left\{ \begin{array}{l} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + \nu(t)x(t) + y(t) = 0 \\ \dot{y}(t) + \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t) = 0 \end{array} \right. \qquad (g\text{VM-DIN-AVD}) $$

with initial conditions $(x(0), y(0)) = (x_0, -\varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) - (\alpha_0 - \varepsilon'_0 - \frac{1}{\beta}\varepsilon_0)x_0)$. One can notice that in the special case where $\varepsilon$ is taken constant and equal to 1 (that is when (VM-DIN-AVD) is simply the DIN-AVD system [11]), we recover the same first-order formulation as that in [11]. For all $t \ge 0$
so that (gVM-DIN-AVD) rewrites $(\dot{x}(t), \dot{y}(t)) = G(t, (x(t), y(t)))$ for all $t \ge 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 $\alpha$ and $\varepsilon$. Since $f$ is coercive and strongly convex on bounded sets, it has a unique minimizer $x^* \in \mathbb{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)\| \underset{t \to +\infty}{\longrightarrow} 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 \ge 0$, and which typically depend on $\varepsilon$ and/or $\alpha$. We first present the case where $\alpha$ is small relative to $\varepsilon$ and then generalize.

## 3.1 Comparison with (CN) under Moderate Viscous Dampings

When the viscous damping $\alpha$ remains moderate w.r.t. the variable mass $\varepsilon$, 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 \ge 0$ such that for all $t \ge 0$, $|\varepsilon'(t)| \le c_1\varepsilon(t)$ and $\alpha(t) \le c_2\varepsilon(t)$.

The assumption states that $\alpha$ must decrease at least as fast as $\varepsilon$ (up to a constant).<fnref n="3" /> The reason for assuming $|\varepsilon'(t)| \le c_1\varepsilon(t)$ is technical and will appear more clearly in the proofs below. It formally means that $\varepsilon$ can decrease at most exponentially fast.<fnref n="4" /> This is a relatively mild assumption that holds, for example, for any polynomial decay $\varepsilon_0/(t+1)^a$, $a \in \mathbb{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 \ge 0$. There exist $C_0, C_1, C_2 \ge 0$, depending on $c_1, c_2$, such that for all $(\varepsilon, \alpha)$ 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 \ge 0$,

$$ \|x(t) - x_N(t)\| \le C_0 e^{-\frac{t}{\beta}} \varepsilon_0 \|\dot{x}_0\| + C_1 \sqrt{\varepsilon(t)} + C_2 \int_{s=0}^{t} e^{\frac{1}{\beta}(s-t)} \sqrt{\varepsilon(s)} \, ds. \quad (1) $$

This extends a previous result from [2, Proposition 3.1] which states a similar bound for constant $\varepsilon$, $\alpha \equiv 0$ and $\beta = 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 $\varepsilon$ and $\alpha$, and that the result is non-asymptotic. This allows in particular choosing $(\varepsilon, \alpha)$ to control the distance from $x$ to $x_N$, for all time $t \ge 0$.

<footnote n="3">³Assumption 1 can actually hold only after some large-enough $t_0 \ge 0$, we take $t_0 = 0$ for the sake of simplicity.</footnote>

<footnote n="4">⁴This is a consequence of Gronwall's lemma, see e.g., [22].</footnote>
and $(u, v) \in \mathbb{R}^n \times \mathbb{R}^n$, define the mapping $G(t, (u, v)) = \begin{pmatrix} \frac{1}{\varepsilon(t)}(-\beta\nabla f(u) - \nu(t)u - v) \\ -\nu'(t)u - \frac{\nu(t)}{\beta}u - \frac{1}{\beta}v \end{pmatrix}$, so that (gVM-DIN-AVD) rewrites $(\dot{x}(t), \dot{y}(t)) = G(t, (x(t), y(t)))$ for all $t \ge 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.2) 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. We conclude with the following important remark.

**Remark 2.2.** *Thanks to the first-order reformulation (gVM-DIN-AVD), Theorem 2.1 does not require the Lipschitz continuity of $\nabla^2 f$ and not even that of $\nabla f$ to obtain global existence and uniqueness. In contrast, these smoothness assumptions would have been required if the standard first-order reformulation of (VM-DIN-AVD) in phase-space (position-velocity) was used. In fact, (gVM-DIN-AVD) allows for a natural extension of (VM-DIN-AVD) to non-smooth convex functions as we detail in Appendix A.3, and which is not the case for (CN) and (LM). This allows defining an inertial Newton-like dynamics in the non-smooth setting which is important for many applications; e.g. [24].*

# 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 $\alpha$ and $\varepsilon$. Since $f$ is coercive and strongly convex on bounded sets, it has a unique minimizer $x^* \in \mathbb{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)\| \xrightarrow[t \to +\infty]{} 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 \ge 0$, and which typically depend on $\varepsilon$ and/or $\alpha$. We first present the case where $\alpha$ is small relative to $\varepsilon$ and then generalize.

## 3.1 Comparison with (CN) under Moderate Viscous Dampings

When the damping $\alpha$ remains moderate w.r.t. the variable mass $\varepsilon$, one expects the solutions of (VM-DIN-AVD) to be close to that of (CN). We make the following assumptions.

**Assumption 2.** *Assumption 1 holds and there exists $c_1, c_2 \ge 0$ such that for all $t \ge 0$, $|\varepsilon'(t)| \le c_1\varepsilon(t)$ and $\alpha(t) \le c_2\varepsilon(t)$.*

The assumption states that $\alpha$ must decrease at least as fast as $\varepsilon$ (up to a constant).<fnref n="3" /> The reason for assuming $|\varepsilon'(t)| \le c_1\varepsilon(t)$ is technical and will appear more clearly in the proofs below. It formally means that $\varepsilon$ can decrease at most exponentially fast.<fnref n="4" /> This is a relatively mild assumption that holds, for example, for any polynomial decay $\varepsilon_0/(t+1)^a$, $a \in \mathbb{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 \ge 0$. There exist $C_0, C_1, C_2 \ge 0$, depending on $c_1, c_2$, such that for all $(\varepsilon, \alpha)$ for which Assumption 2 holds with constants $c_1$ and $c_2$, the

<footnote n="3">³Assumption 2 can hold only after some $t_0 \ge 0$, we take $t_0 = 0$ for the sake of simplicity.</footnote>

<footnote n="4">⁴This is a consequence of Gronwall's lemma, see e.g., [27].</footnote>
**Remark 3.3.** Under Assumption 1, the dynamics (VM-DIN-AVD) is dominated by the variable mass. The damping $\alpha$ 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 $(\varepsilon, \alpha)$, and let $x$ be the corresponding solution of (VM-DIN-AVD). For all $t \ge 0$, define the function,

$$U(t) = \frac{\varepsilon(t)}{2} \|\dot{x}(t)\|^2 + f(x(t)) - f(x^*).$$

Then $U$ is differentiable and for all $t > 0$,

$$\frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 - \alpha(t) \|\dot{x}(t)\|^2 - \beta \langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle \le 0.$$

Therefore, in particular, $U$ is non-increasing.

*Proof.* Let $t \ge 0$, since $x$ is twice differentiable, $U$ is differentiable and,

$$\frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 + \varepsilon(t) \langle \dot{x}(t), \ddot{x}(t) \rangle + \langle \dot{x}(t), \nabla f(x(t)) \rangle.$$

We use the fact that $x$ is solution of (VM-DIN-AVD), to substitute $\varepsilon(t)\ddot{x}(t)$ by its expression,

$$\frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 - \alpha(t) \|\dot{x}(t)\|^2 - \beta \langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle.$$

By assumption $\varepsilon$ is non-increasing so for all $t > 0$, $\varepsilon'(t) \le 0$. Furthermore $f$ is convex so $\langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle \ge 0$. Hence $U$ is non-increasing. $\square$

We then state the following bound.

**Lemma 3.5.** There exists $C \ge 0$ such that for any $(\varepsilon, \alpha)$ and the corresponding solution $x$ of (VM-DIN-AVD), for all $t \ge 0$ it holds that,

$$\varepsilon(t) \|\dot{x}(t)\| \le C \sqrt{\varepsilon(t)}.$$

*Proof.* Let $t \ge 0$, according to Lemma 3.4, $U$ is non-increasing so $U(t) \le U(0)$, or equivalently,

$$\frac{\varepsilon(t)}{2} \|\dot{x}(t)\|^2 + f(x(t)) - f(x^*) \le \frac{\varepsilon_0}{2} \|\dot{x}_0\|^2 + f(x_0) - f(x^*).$$

This implies in particular that,

$$\varepsilon(t) \|\dot{x}(t)\|^2 \le \varepsilon_0 \|\dot{x}_0\|^2 + 2(f(x_0) - f(x^*)),$$

and hence by multiplying both sides by $\varepsilon(t)$ and composing with the square-root we obtain that,

$$\varepsilon(t) \|\dot{x}(t)\| \le C \sqrt{\varepsilon(t)},$$

$$\text{where } C = \sqrt{\varepsilon_0 \|\dot{x}_0\|^2 + 2(f(x_0) - f(x^*))}.$$

$\square$
corresponding solution $x$ of (VM-DIN-AVD) is such that for all $t \ge 0$,

$$ \|x(t) - x_N(t)\| \le C_0 e^{-\frac{t}{\beta}} \varepsilon_0 \|\dot{x}_0\| + C_1 \sqrt{\varepsilon(t)} + C_2 \int_{s=0}^{t} e^{\frac{1}{\beta}(s-t)} \sqrt{\varepsilon(s)} \, ds. \quad (1) $$

This extends a previous result from [3, Proposition 3.1] which states a similar bound for constant $\varepsilon$, $\alpha \equiv 0$ and $\beta = 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 result comes from the fact that $C_0, C_1, C_2$ do not depend on $\varepsilon$ and $\alpha$, and that the result is non-asymptotic. This allows in particular for choosing $(\varepsilon, \alpha)$ to control the distance from $x$ to $x_N$, for all time $t \ge 0$.

**Remark 3.3.** Under Assumption 2, the dynamics (VM-DIN-AVD) is dominated by the variable mass. The damping $\alpha$ 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 $(\varepsilon, \alpha)$, and let $x$ be the corresponding solution of VM-DIN-AVD. For all $t \ge 0$, define the function, $U(t) = \frac{\varepsilon(t)}{2} \|\dot{x}(t)\|^2 + f(x(t)) - f(x^*)$. Then $U$ is differentiable and for all $t > 0$,

$$ \frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 - \alpha(t) \|\dot{x}(t)\|^2 - \beta \langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle \le 0. $$

Therefore, in particular, $U$ is non-increasing.

*Proof.* Let $t \ge 0$, since $x$ is twice differentiable, $U$ is differentiable and,

$$ \frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 + \varepsilon(t) \langle \dot{x}(t), \ddot{x}(t) \rangle + \langle \dot{x}(t), \nabla f(x(t)) \rangle. $$

We use the fact that $x$ is solution of (VM-DIN-AVD), to substitute $\varepsilon(t)\ddot{x}(t)$ by its expression, $\frac{dU}{dt}(t) = \frac{\varepsilon'(t)}{2} \|\dot{x}(t)\|^2 - \alpha(t) \|\dot{x}(t)\|^2 - \beta \langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle$. By assumption $\varepsilon$ is non-increasing so for all $t > 0$, $\varepsilon'(t) \le 0$. Furthermore $f$ is convex so $\langle \nabla^2 f(x(t)) \dot{x}(t), \dot{x}(t) \rangle \ge 0$. Hence $U$ is non-increasing. $\square$

We then state the following bound.

**Lemma 3.5.** For any $(\varepsilon, \alpha)$, the corresponding solution $x$ of (VM-DIN-AVD) is such that for all $t \ge 0$,

$$ \varepsilon(t) \|\dot{x}(t)\| \le \sqrt{2U(0)} \sqrt{\varepsilon(t)}. $$

From Lemma 3.4, $U$ is non-increasing so $\forall t \ge 0$, $U(t) \le U(0)$. Then in particular $\frac{\varepsilon(t)}{2} \|\dot{x}(t)\|^2 \le U(0)$ and the proof follows by multiplying both sides by $\varepsilon(t)$.

*Proof of Theorem 3.1.* Let $(\varepsilon, \alpha)$ 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 \ge 0$, $U(t) \le U(0)$, so in particular $f(x(t)) \le f(x_0) + \frac{\varepsilon_0}{2} \|\dot{x}_0\|^2$. Denoting $M_0 = f(x_0) + \frac{\varepsilon_0}{2} \|\dot{x}_0\|^2$, the set $K_0 = \{y \in \mathbb{R}^n | f(y) \le M_0\}$ is bounded, since $f$ is coercive ($\lim_{\|y\| \to +\infty} f(y) = +\infty$). So for all $t \ge 0$, $x(t) \in K_0$. Since $M_0$ (and hence $K_0$) depends only on $\varepsilon_0$, $x_0$ and $\dot{x}_0$, we have proved that for any choice $(\varepsilon, \alpha)$, the corresponding
*Proof of Theorem 3.1.* Let $(\varepsilon, \alpha)$ 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 \ge 0$, $U(t) \le U(0)$, so in particular

$$f(x(t)) \le f(x_0) + \frac{\varepsilon_0}{2} \|\dot{x}_0\|^2.$$

Denoting $M_0 = f(x_0) + \frac{\varepsilon_0}{2} \|\dot{x}_0\|^2$, the set $K_0 = \{y \in \mathbb{R}^n \mid f(y) \le M_0\}$ is bounded, since $f$ is coercive ($\lim_{\|y\| \to +\infty} f(y) = +\infty$). So for all $t \ge 0$, $x(t) \in K_0$. Since $M_0$ (and hence $K_0$) depends only on $\varepsilon_0$, $x_0$ and $\dot{x}_0$, we have proved that for any choice $(\varepsilon, \alpha)$, 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 \ge 0$, $f(x_N(t)) \le f(x_N(0)) = f(x_0) \le M_0$. So we also have $x_N(t) \in K_0$ for all $t \ge 0$.

Now, fix $c_1, c_2 > 0$, and let $(\varepsilon, \alpha)$ 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 $\mu > 0$ the strong-convexity parameter of $f$ on $K_0$. Equivalently, we have that $\nabla f$ is strongly monotone on $K_0$, that is, $\forall y_1, y_2 \in K_0$,

$$\langle \nabla f(y_1) - \nabla f(y_2), y_1 - y_2 \rangle \ge \mu \|y_1 - y_2\|^2. \quad (2)$$

Let $t \ge 0$, since $x(t) \in K_0$ and $x_N(t) \in K_0$, by combining (2) with the Cauchy–Schwarz inequality, we deduce that

$$\|x(t) - x_N(t)\| \le \frac{1}{\mu} \|\nabla f(x(t)) - \nabla f(x_N(t))\|. \quad (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: $\frac{d}{dt}\nabla f(x_N(t)) + \frac{1}{\beta}\nabla f(x_N(t)) = 0$. So we can integrate, for all $t \ge 0$,

$$\nabla f(x_N(t)) = e^{-\frac{t}{\beta}} \nabla f(x_0). \quad (4)$$

We now turn our attention to $\nabla f(x(t))$, for which we cannot find a closed-form solution in general. We rewrite (VM-DIN-AVD) in the following equivalent form

$$\frac{d}{dt} [\varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t))] + \frac{1}{\beta}\varepsilon(t)\dot{x}(t) + \nabla f(x(t)) = \left(\frac{1}{\beta}\varepsilon(t) + \varepsilon'(t) - \alpha(t)\right) \dot{x}(t).$$

Introducing the variable $\omega(t) = \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t))$, the latter is thus solution to

$$\begin{cases} \dot{\omega}(t) + \frac{1}{\beta}\omega(t) = \left(\frac{1}{\beta}\varepsilon(t) + \varepsilon'(t) - \alpha(t)\right) \dot{x}(t), & t \ge 0, \\ \omega(0) = \varepsilon_0\dot{x}_0 + \beta\nabla f(x_0). \end{cases}$$

This is a non-homogeneous first-order ODE in $\omega$, whose solution can be expressed using the integrating factor

$$\omega(t) = e^{-\frac{t}{\beta}}(\varepsilon_0\dot{x}_0 + \beta\nabla f(x_0)) + e^{-\frac{t}{\beta}}\int_0^t e^{\frac{s}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s) - \alpha(s)\right) \dot{x}(s) \,ds.$$

We thus have the following expression for $\nabla f(x)$, for all $t \ge 0$,

$$\beta\nabla f(x(t)) = \beta e^{-\frac{t}{\beta}}\nabla f(x_0) + e^{-\frac{t}{\beta}}\varepsilon_0\dot{x}_0 - \varepsilon(t)\dot{x}(t) + e^{-\frac{t}{\beta}}\int_0^t e^{\frac{s}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s) - \alpha(s)\right) \dot{x}(s) \,ds. \quad (5)$$
solution $x$ of (VM-DIN-AVD) is inside $K_0$ at all times. Let $x_N$ be the solution of (CN). One can similarly see that for all $t \ge 0$, $f(x_N(t)) \le f(x_N(0)) = f(x_0) \le M_0$. So we also have $x_N(t) \in K_0$ for all $t \ge 0$.

Now, fix $c_1, c_2 > 0$, and let $(\varepsilon, \alpha)$ such that Assumption 2 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 $\mu_{K_0} > 0$ the strong-convexity parameter of $f$ on $K_0$. Equivalently, we have that $\nabla f$ is strongly monotone on $K_0$, that is, $\forall y_1, y_2 \in K_0$,

$$ \langle \nabla f(y_1) - \nabla f(y_2), y_1 - y_2 \rangle \ge \mu_{K_0} \|y_1 - y_2\|^2. \quad (2) $$

Let $t \ge 0$, since $x(t) \in K_0$ and $x_N(t) \in K_0$, by combining (2) with the Cauchy–Schwarz inequality, we deduce that

$$ \|x(t) - x_N(t)\| \le \frac{1}{\mu_{K_0}} \|\nabla f(x(t)) - \nabla f(x_N(t))\|. \quad (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: $\frac{d}{dt}\nabla f(x_N(t)) + \frac{1}{\beta}\nabla f(x_N(t)) = 0$. So we can integrate, for all $t \ge 0$,

$$ \nabla f(x_N(t)) = e^{-\frac{t}{\beta}} \nabla f(x_0). \quad (4) $$

We now turn our attention to $\nabla f(x(t))$, for which we cannot find a closed-form solution in general. We rewrite (VM-DIN-AVD) in the following equivalent form

$$ \frac{d}{dt} [\varepsilon(t) \dot{x}(t) + \beta \nabla f(x(t))] + \frac{1}{\beta} \varepsilon(t) \dot{x}(t) + \nabla f(x(t)) = \left( \frac{1}{\beta} \varepsilon(t) + \varepsilon'(t) - \alpha(t) \right) \dot{x}(t). $$

Introducing the variable $\omega(t) = \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t))$, the latter is thus solution to

$$ \dot{\omega}(t) + \frac{1}{\beta}\omega(t) = \left( \frac{1}{\beta}\varepsilon(t) + \varepsilon'(t) - \alpha(t) \right) \dot{x}(t), \quad t \ge 0, $$

with $\omega(0) = \varepsilon_0\dot{x}_0 + \beta\nabla f(x_0)$. This is a non-homogeneous first-order ODE in $\omega$, whose solution can be expressed using the integrating factor

$$ \omega(t) = e^{-\frac{t}{\beta}} (\varepsilon_0 \dot{x}_0 + \beta \nabla f(x_0)) + e^{-\frac{t}{\beta}} \int_0^t e^{\frac{s}{\beta}} \left( \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) - \alpha(s) \right) \dot{x}(s) \, ds. $$

We thus have the following expression for $\nabla f(x)$, for all $t \ge 0$,

$$ \beta \nabla f(x(t)) = \beta e^{-\frac{t}{\beta}} \nabla f(x_0) + e^{-\frac{t}{\beta}} \varepsilon_0 \dot{x}_0 - \varepsilon(t) \dot{x}(t) \\ + e^{-\frac{t}{\beta}} \int_0^t e^{\frac{s}{\beta}} \left( \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) - \alpha(s) \right) \dot{x}(s) \, ds. \quad (5) $$

We can now use (4) and (5) in (3) to get

$$ \|x(t) - x_N(t)\| \le \frac{1}{\beta \mu_{K_0}} \left\| e^{-\frac{t}{\beta}} \varepsilon_0 \dot{x}_0 - \varepsilon(t) \dot{x}(t) + e^{-\frac{t}{\beta}} \int_0^t e^{\frac{s}{\beta}} \left( \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) - \alpha(s) \right) \dot{x}(s) \, ds \right\|. $$
We can now use (4) and (5) in (3) to get

$$ \|x(t) - x_N(t)\| \le \frac{1}{\beta\mu} \left\| e^{-\frac{t}{\beta}}\varepsilon_0 \dot{x}_0 - \varepsilon(t)\dot{x}(t) + e^{-\frac{t}{\beta}} \int_0^t e^{\frac{s}{\beta}} \left( \frac{1}{\beta}\varepsilon(s) + \varepsilon'(s) - \alpha(s) \right) \dot{x}(s) \, ds \right\|. $$

Using the triangle inequality, we obtain,

$$ \|x(t) - x_N(t)\| \le \frac{\varepsilon_0 \|\dot{x}_0\|}{\beta\mu} e^{-\frac{t}{\beta}} + \frac{\varepsilon(t) \|\dot{x}(t)\|}{\beta\mu} + \frac{1}{\beta\mu} \int_0^t e^{\frac{1}{\beta}(s-t)} \left| \frac{1}{\beta}\varepsilon(s) + \varepsilon'(s) - \alpha(s) \right| \|\dot{x}(s)\| \, ds. \quad (6) $$

The first term in (6) corresponds to the first one in (1) with $C_0 = 1/(\beta\mu)$. As for the second-one, by direct application of Lemma 3.5, there exists $C > 0$ such that for all $t \ge 0$, $\frac{\varepsilon(t)\|\dot{x}(t)\|}{\beta\mu} \le C\sqrt{\varepsilon(t)}$, so we set $C_1 = C/(\beta\mu)$. Regarding the last term in (6), using Assumption 1 and again Lemma 3.5, it holds that, for all $s \ge 0$,

$$ \left| \frac{1}{\beta}\varepsilon(s) + \varepsilon'(s) - \alpha(s) \right| \|\dot{x}(s)\| \le \left( \frac{1}{\beta} + c_1 + c_2 \right) \varepsilon(s) \|\dot{x}(s)\| \le \left( \frac{1}{\beta} + c_1 + c_2 \right) C\sqrt{\varepsilon(s)}. $$

This proves the theorem with $C_2 = \frac{C}{\beta\mu} (\frac{1}{\beta} + c_1 + c_2)$.

<figure><img src="image_4.png" /></figure>

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 $\dot{x}_0 = 0$, the second-one decays like $\sqrt{\varepsilon(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 < \frac{2}{\beta}$, then there exists $C_3 > 0$ such that for all $t \ge 0$,

$$ \|x(t) - x_N(t)\| \le C_0 e^{-\frac{t}{\beta}} \varepsilon_0 \|\dot{x}_0\| + C_3 \sqrt{\varepsilon(t)}. $$

*Proof of Corollary 3.6.* For all $t \ge 0$, define $J(t) = \int_0^t e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \, ds$. An integration by parts yields

$$ J(t) = \left[ \beta e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \right]_{s=0}^t - \int_{s=0}^t \beta e^{\frac{s}{\beta}} \frac{\varepsilon'(s)}{2\sqrt{\varepsilon(s)}} \, ds = \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} - \beta \varepsilon_0 + \int_{s=0}^t \beta e^{\frac{s}{\beta}} \frac{-\varepsilon'(s)}{2\varepsilon(s)} \sqrt{\varepsilon(s)} \, ds. \quad (7) $$

By assumption, $0 \le c_1 < \frac{2}{\beta}$ such that for all $s > 0$, $|\varepsilon'(s)| \le c_1\varepsilon(s)$, which in our setting is equivalent to $\frac{-\varepsilon'(s)}{\varepsilon(s)} \le c_1$. So we deduce from (7) that

$$ J(t) \le \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} + c_1 \frac{\beta}{2} \int_{s=0}^{t} e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \, ds = \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} + c_1 \frac{\beta}{2} J(t). $$

So, $(1 - c_1 \frac{\beta}{2}) J(t) \le \beta e^{t} \sqrt{\varepsilon(t)}$. By assumption $1 - c_1 \frac{\beta}{2} > 0$, therefore, $J(t) \le \frac{2}{2 - c_1 \beta} e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)}$. Finally, using this in (1) and setting $C_3 = C_1 + C_2 \frac{2}{2 - c_1 \beta}$, we obtain the result. □

<figure><img src="image_5.png" /></figure>

**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 $\alpha$ is “not too large” w.r.t. $\varepsilon$, and do not study (LM). We now state a more general result that covers these cases.
Using the triangle inequality, we obtain,

$$ \|x(t) - x_N(t)\| \le \frac{\varepsilon_0 \|\dot{x}_0\|}{\beta \mu_{\kappa_0}} e^{-\frac{t}{\beta}} + \frac{\varepsilon(t) \|\dot{x}(t)\|}{\beta \mu_{\kappa_0}} \\ + \frac{1}{\beta \mu_{\kappa_0}} \int_0^t e^{\frac{1}{\beta}(s-t)} \left| \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) - \alpha(s) \right| \|\dot{x}(s)\| \, ds. \quad (6) $$

The first term in (6) corresponds to the first one in (1) with $C_0 = 1/(\beta\mu_{\kappa_0})$. As for the second-one, by direct application of Lemma 3.5, for all $t \ge 0$, $\varepsilon(t)\|\dot{x}(t)\| \le \sqrt{2U(0)}\sqrt{\varepsilon(t)}$, so we set $C_1 = \sqrt{2U(0)}/(\beta\mu_{\kappa_0})$. Regarding the last term in (6), using Assumption 2 and again Lemma 3.5, it holds that, for all $s \ge 0$,

$$ \left| \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) - \alpha(s) \right| \|\dot{x}(s)\| \le \left( \frac{1}{\beta} + c_1 + c_2 \right) \varepsilon(s) \|\dot{x}(s)\| \le \left( \frac{1}{\beta} + c_1 + c_2 \right) \sqrt{2U(0)\varepsilon(s)}. $$

This proves the theorem with $C_2 = \frac{\sqrt{2U(0)}}{\beta\mu_{\kappa_0}} (\frac{1}{\beta} + c_1 + c_2)$. $\square$

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 $\dot{x}_0 = 0$, the second-one decays like $\sqrt{\varepsilon(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.** Let the same assumptions and variables as in Theorem 3.1. If furthermore $c_1 < \frac{2}{\beta}$, then there exists $C_3 > 0$ such that for all $t \ge 0$,

$$ \|x(t) - x_N(t)\| \le C_0 e^{-\frac{t}{\beta}} \varepsilon_0 \|\dot{x}_0\| + C_3 \sqrt{\varepsilon(t)}. $$

*Proof of Corollary 3.6.* For all $t \ge 0$, define $J(t) = \int_0^t e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \, ds$. An integration by parts yields

$$ \begin{aligned} J(t) &= \left[ \beta e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \right]_{s=0}^t - \int_{s=0}^t \beta e^{\frac{s}{\beta}} \frac{\varepsilon'(s)}{2\sqrt{\varepsilon(s)}} \, ds \\ &= \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} - \beta \varepsilon_0 + \int_{s=0}^t \beta e^{\frac{s}{\beta}} \frac{-\varepsilon'(s)}{2\varepsilon(s)} \sqrt{\varepsilon(s)} \, ds. \quad (7) \end{aligned} $$

By assumption, $0 \le c_1 < \frac{2}{\beta}$ such that for all $s > 0$, $|\varepsilon'(s)| \le c_1\varepsilon(s)$, which in our setting is equivalent to $\frac{-\varepsilon'(s)}{\varepsilon(s)} \le c_1$. So we deduce from (7) that

$$ J(t) \le \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} + c_1 \frac{\beta}{2} \int_{s=0}^{t} e^{\frac{s}{\beta}} \sqrt{\varepsilon(s)} \, ds = \beta e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)} + c_1 \frac{\beta}{2} J(t). $$

So, $(1 - c_1 \frac{\beta}{2}) J(t) \le \beta e^{t} \sqrt{\varepsilon(t)}$. By assumption $1 - c_1 \frac{\beta}{2} > 0$, therefore, $J(t) \le \frac{2}{2-c_1\beta} e^{\frac{t}{\beta}} \sqrt{\varepsilon(t)}$. We use this in (1) and set $C_3 = C_1 + C_2 \frac{2}{2-c_1\beta}$ to get the result. $\square$

**Remark 3.7.** Observe that local strong convexity is only required to get (2) and (3). In fact, local strong convexity can be greatly weakened and our claims above can be generalized to a large sub-class of strictly convex functions as we explain in Appendix B. We did not directly consider the strictly convex
## 3.2 Generalization to Arbitrary Viscous Dampings with Sub-exponential Decay

This time we do not assume a link between $\varepsilon$ and $\alpha$ but only sub-exponential decays.

**Assumption 2.** There exists $c_1, c_2 \ge 0$ such that for all $t \ge 0$, $|\varepsilon'(t)| \le c_1\varepsilon(t)$ and $|\alpha'(t)| \le c_2\alpha(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 \ge 0$. There exist constants $C, \tilde{C} \ge 0$, depending on $c_1, c_2$, such that for all $\varepsilon$ and $\alpha$ 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 \ge 0$,

$$ \|x(t) - x_N(t)\| \le C \left[ e^{-\frac{t}{\beta}} + \sqrt{\varepsilon(t)} + \alpha(t) + \int_{s=0}^{t} e^{\frac{1}{\beta}(s-t)} (\sqrt{\varepsilon(s)} + \alpha(s)) \, ds \right], \quad (8) $$

and,

$$ \|x(t) - x_{LM}(t)\| \le \tilde{C} \left[ e^{-\frac{t}{\beta}} + \sqrt{\varepsilon(t)} + \alpha(t) + \int_{s=0}^{t} e^{\frac{1}{\beta}(s-t)} (\sqrt{\varepsilon(s)} + \alpha(s)) \, ds \right]. \quad (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 $\varepsilon$ and $\alpha$ 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 $\alpha$, suggesting that making $x$ close to $x_{LM}$ requires not only $\varepsilon$ but also $\alpha$ 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 $\alpha$, 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 $\varepsilon(t) \le \alpha(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 $\varepsilon$ and $\alpha$ have on the speed of convergence of $x$ to $x^*$ compared to the speeds of convergence of $x_N$ and $x_{LM}$.

## 4 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) = \frac{1}{2}\|Ay - b\|_2^2$ for all $y \in \mathbb{R}^n$,
case to emphasize the main ideas and ease the reading. Following Remark 2.2, it seems also that a non-smooth extension is possible using regularization techniques. Unlike the strictly convex setting, the non-smooth one would however require further investigations since the vanilla Newton method is not applicable to non-smooth functions.

So far our results only cover the case where $\alpha$ is “not too large” w.r.t. $\varepsilon$, and do not study (LM). We now state a more general result that covers these cases.

## 3.2 Generalization to Sub-exponentially Decaying Viscous Dampings

This time we do not assume a link between $\varepsilon$ and $\alpha$ but only sub-exponential decays.

**Assumption 3.** Assumption 1 holds and there exists $c_1, c_2 \ge 0$ such that for all $t \ge 0$, $|\varepsilon'(t)| \le c_1\varepsilon(t)$ and $|\alpha'(t)| \le c_2\alpha(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 solutions of (CN) and (LM) respectively, and let $c_1, c_2 \ge 0$. There exist constants $C, \tilde{C} \ge 0$, depending on $c_1, c_2, \varepsilon_0, \alpha_0$ and the initial conditions, such that for all $\varepsilon$ and $\alpha$ for which Assumption 3 holds with $c_1, c_2$, the corresponding solution $x$ of (VM-DIN-AVD) is such that for all $t \ge 0$,

$$ \|x(t) - x_N(t)\| \le C \left[ e^{-\frac{t}{\tilde{\beta}}} + \sqrt{\varepsilon(t)} + \alpha(t) + \int_{s=0}^{t} e^{\frac{1}{\tilde{\beta}}(s-t)} (\sqrt{\varepsilon(s)} + \alpha(s)) \, ds \right], \quad (8) $$

and,

$$ \|x(t) - x_{LM}(t)\| \le \tilde{C} \left[ e^{-\frac{t}{\tilde{\beta}}} + \sqrt{\varepsilon(t)} + \alpha(t) + \int_{s=0}^{t} e^{\frac{1}{\tilde{\beta}}(s-t)} (\sqrt{\varepsilon(s)} + \alpha(s)) \, ds \right]. \quad (9) $$

The proof is omitted but key elements are presented in Appendix C. 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 C). Theorem 3.8 allows for far more flexibility in the choice of $\varepsilon$ and $\alpha$ 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 $\alpha$, suggesting that making $x$ close to $x_{LM}$ requires not only $\varepsilon$ but also $\alpha$ 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 $\alpha$, 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 $\varepsilon(t) \le \alpha(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 $\varepsilon$ and $\alpha$ have on the speed of convergence of $x$ to $x^*$ compared to the speeds of convergence of $x_N$ and $x_{LM}$.
where $A \in \mathbb{R}^{n \times n}$ is symmetric positive definite and $b \in \mathbb{R}^n$. Without loss of generality, we take $b = 0$, so that the unique minimum is $x^* = 0$.

## 4.1 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 [10, 40]). Indeed, $\forall y \in \mathbb{R}^n$, $\nabla f(y) = A^T A y$ and $\nabla^2 f(y) = A^T A$. Since $\nabla^2 f(y)$ is independent of $y$ we can rewrite (VM-DIN-AVD) in an eigenspace$^5$ of $A^T A$. That is, we can study the ODE coordinate-wise by looking at one-dimensional problems of the form

$$ \varepsilon(t)\ddot{x}(t) + (\alpha(t) + \beta\lambda)\dot{x}(t) + \lambda x(t) = 0, \quad t \ge 0. \qquad (\text{Q1-VM-DIN-AVD}) $$

Here (and throughout what follows) $\lambda > 0$ denotes any eigenvalue of $A^T A$ and $x: \mathbb{R}_+ \to \mathbb{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

$$ (\alpha(t) + \beta\lambda)\dot{x}_{LM}(t) + \lambda x_{LM}(t) = 0, \quad t \ge 0, \qquad (\text{Q1-LM}) $$

where $x_{LM}: \mathbb{R}_+ \to \mathbb{R}$, and (CN) becomes

$$ \beta\dot{x}_N(t) + x_N(t) = 0, \quad t \ge 0, \qquad (\text{Q1-CN}) $$

where again, $x_N: \mathbb{R}_+ \to \mathbb{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 \ge 0$,

$$ x_N(t) = x_0 e^{-\frac{t}{\beta}} \quad \text{and} \quad x_{LM}(t) = x_0 \exp \left( - \int_0^t \frac{\lambda}{\alpha(s) + \beta\lambda} ds \right). \qquad (10) $$

Since the minimizer is $x^* = 0$, we directly see that $x_N$ converges exponentially fast to $x^*$, with a rate independent of $\lambda$. The rate of $x_{LM}$ depends however on $\lambda$ and how fast $\alpha$ vanishes.

Unfortunately, except for some special choices of $\varepsilon$ and $\alpha$ (see [10]), 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 $\nabla f(x)$, accordingly to what we did in Section 3, since here $\nabla f(x) = \lambda 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.

## 4.2 The Liouville–Green Method

In what follows, we rely on the Liouville–Green method [30, 26], a technique for obtaining non-asymptotic approximations to solutions of linear second-order ODEs with non-constant coefficients. First, we give the intuition behind the method, following the presentation of [35]. Consider the differential equation

$$ \ddot{z}(t) - r(t)z(t) = 0, \quad t \ge 0, \qquad (11) $$

<footnote>$^5$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$.</footnote>
# 4 Approximate Solutions and Asymptotic Analysis on Quadratics

In addition to Assumption 1, we consider the 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) = \frac{1}{2}\|Ay - b\|_2^2$ for all $y \in \mathbb{R}^n$, where $A \in \mathbb{R}^{n \times n}$ is symmetric positive definite and $b \in \mathbb{R}^n$. Without loss of generality, we take $b = 0$, so that the unique minimum is $x^* = 0$.

## 4.1 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 [12, 46]). Indeed, $\forall y \in \mathbb{R}^n$, $\nabla f(y) = A^T Ay$ and $\nabla^2 f(y) = A^T A$. Since $\nabla^2 f(y)$ is independent of $y$ we can rewrite (VM-DIN-AVD) in an eigenspace<fnref n="5" /> of $A^T A$. That is, we can study the ODE coordinate-wise by looking at one-dimensional problems of the form

$$ \varepsilon(t)\ddot{x}(t) + (\alpha(t) + \beta\lambda)\dot{x}(t) + \lambda x(t) = 0, \quad t \ge 0. \qquad (\text{Q1-VM-DIN-AVD}) $$

Here (and throughout what follows) $\lambda > 0$ denotes any eigenvalue of $A^T A$ and $x: \mathbb{R}_+ \to \mathbb{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

$$ (\alpha(t) + \beta\lambda)\dot{x}_{LM}(t) + \lambda x_{LM}(t) = 0, \quad t \ge 0, \qquad (\text{Q1-LM}) $$

where $x_{LM}: \mathbb{R}_+ \to \mathbb{R}$, and (CN) becomes

$$ \beta\dot{x}_N(t) + x_N(t) = 0, \quad t \ge 0, \qquad (\text{Q1-CN}) $$

where again, $x_N: \mathbb{R}_+ \to \mathbb{R}$ is one-dimensional. Observe in particular that (CN) and (LM) are now first-order linear ODEs, whose solutions have closed forms: $\forall t \ge 0$,

$$ x_N(t) = x_0 e^{-\frac{t}{\beta}} \quad \text{and} \quad x_{LM}(t) = x_0 \exp \left( - \int_0^t \frac{\lambda}{\alpha(s) + \beta\lambda} ds \right). \qquad (10) $$

Since the minimizer is $x^* = 0$, we see that $x_N$ converges exponentially fast to $x^*$, with a rate independent of $\lambda$ while the rate of $x_{LM}$ depends on $\lambda$ and how fast $\alpha$ vanishes.

Unfortunately, except for some special choices of $\varepsilon$ and $\alpha$ (see [12]), 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 $\nabla f(x)$, accordingly to what we did in Section 3, since here $\nabla f(x) = \lambda 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.

<footnote n="5">⁵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$.</footnote>
where $r$ is real-valued, positive, and twice continuously differentiable. Any linear second-order ODE can be reformulated in the form (11), see Lemma 4.5 below. Since for all $t \ge 0$, $r(t) \ne 0$, we can use the changes of variables $\tau = \int_0^t \sqrt{r(s)} ds$ and $w = r^{1/4}z$ and show that $w$ is solution to

$$ \ddot{w}(\tau) - (1 + \psi(\tau))w(\tau) = 0, \quad t \ge 0, \qquad (12) $$

where<fnref n="6" /> $\psi(\tau) = \frac{4r(t)r''(t)-5r'(t)^2}{16r(t)^3}$. The LG method consists in neglecting the term $\psi(\tau)$ in (12), which simply yields two approximate solutions $\hat{w}_1(\tau) = e^\tau$ and $\hat{w}_2(\tau) = e^{-\tau}$. Expressing this in terms of $z$ and $t$, we obtain

$$ \hat{z}_1(t) = r(t)^{-1/4} \exp \left( \int_0^t \sqrt{r(s)} \, ds \right) \quad \text{and} \quad \hat{z}_2(t) = r(t)^{-1/4} \exp \left( \int_0^t -\sqrt{r(s)} \, ds \right). \qquad (13) $$

Those are the LG approximations of the solutions of (11). They are formally valid on any interval $[0, T]$, $T > 0$ when $\psi$ is “not too large”, provided that $\sqrt{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., [14, Chapter 10].

The advantage of this approach is the possibility to estimate the error made using (13) w.r.t. the true solutions of (11). This is expressed in the following theorem which gathers results from [15, 36, 42].

**Theorem 4.2** (Olver [35]). Let $r: \mathbb{R}_+ \to \mathbb{R}$ be a real, positive, twice continuously differentiable function, and define $\varphi(t) = \frac{4r(t)r''(t)-5r'(t)^2}{16r(t)^{5/2}}$ for all $t \ge 0$. Then for any $T > 0$, the differential equation,

$$ \ddot{z}(t) - r(t)z(t) = 0, \quad t \in [0, T], \qquad (14) $$

has two real and twice continuously differentiable solutions defined for all $t \in [0, T]$ by,

$$ z_1(t) = \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t \sqrt{r(s)} \, ds \right) \quad \text{and} \quad z_2(t) = \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( -\int_0^t \sqrt{r(s)} \, ds \right), $$

where $|\delta_1(t)| \le \exp\left(\frac{1}{2}\int_0^t |\varphi(s)| ds\right) - 1$ and $|\delta_2(t)| \le \exp\left(-\frac{1}{2}\int_t^T |\varphi(s)| ds\right) - 1$. If in addition $\int_0^{+\infty} |\varphi(s)| ds < +\infty$, then the results above also hold for $T = +\infty$.

**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 (14). The LG approximations $\hat{z}_1$ and $\hat{z}_2$ are obtained by neglecting the unknown functions $\delta_1$ and $\delta_2$ in $z_1$ and $z_2$. The theorem gives a non-asymptotic bound for the errors $|z_1(t) - \hat{z}_1(t)|$ and $|z_2(t) - \hat{z}_2(t)|$, $t \ge 0$.

- Since we assumed $r$ to be twice continuously differentiable and positive, $\varphi$ is continuous, so it is integrable except maybe for $t \to +\infty$.

- For the sake of simplicity, the formulation of Theorem 4.2 slightly differs from that in [35], the original formulation can be recovered by a change of variable.

<footnote n="6">⁶We express $\psi(\tau)$ using $t$ for the sake of readability, using the one-to-one correspondence between $\tau$ and $t$.</footnote>
## 4.2 The Liouville–Green Method

In what follows, we rely on the Liouville–Green method [35, 31], a technique for obtaining non-asymptotic approximations to solutions of linear second-order ODEs with non-constant coefficients. First, we give the intuition behind the method, following the presentation of [42]. Consider the differential equation

$$ \ddot{z}(t) - r(t)z(t) = 0, \quad t \ge 0, \qquad (11) $$

where $r$ is real-valued, positive, and twice continuously differentiable. Any linear second-order ODE can be reformulated in the form (11), see Lemma 4.5 below. Since for all $t \ge 0$, $r(t) \ne 0$, we can use the changes of variables $\tau = \int_0^t \sqrt{r(s)} ds$ and $w = r^{1/4}z$ and show that $w$ is solution to

$$ \ddot{w}(\tau) - (1 + \psi(\tau))w(\tau) = 0, \quad t \ge 0, \qquad (12) $$

where$^6$ $\psi(\tau) = \frac{4r(t)r''(t)-5r'(t)^2}{16r(t)^3}$. The LG method consists in neglecting the term $\psi(\tau)$ in (12), which simply yields two approximate solutions $\hat{w}_1(\tau) = e^\tau$ and $\hat{w}_2(\tau) = e^{-\tau}$. Expressing this in terms of $z$ and $t$, we obtain

$$ \hat{z}_1(t) = r(t)^{-1/4} \exp \left( \int_0^t \sqrt{r(s)} \, ds \right) \quad \text{and} \quad \hat{z}_2(t) = r(t)^{-1/4} \exp \left( \int_0^t -\sqrt{r(s)} \, ds \right). \qquad (13) $$

Those are the LG approximations of the solutions of (11). They are formally valid on any $[0, T]$, $T > 0$ when $\psi$ is “not too large” and if $\sqrt{r}$ is integrable on $[0, T]$.

**Remark 4.1.** There exists other (but less intuitive) ways to derive the LG approximations which allow for generalization to higher-order linear ODEs [17, Chapter 10].

The advantage of this approach is the possibility to estimate the error made using (13) w.r.t. the true solutions of (11). This is expressed in the following theorem which gathers results from [18, 41, 48].

**Theorem 4.2** ([42]). Let $r: \mathbb{R}_+ \to \mathbb{R}$ be a real, positive, twice continuously differentiable function, and define $\varphi(t) = \frac{4r(t)r''(t)-5r'(t)^2}{16r(t)^{5/2}}$ for all $t \ge 0$. Then for any $T > 0$, the differential equation,

$$ \ddot{z}(t) - r(t)z(t) = 0, \quad t \in [0, T], \qquad (14) $$

has two real and twice continuously differentiable solutions defined $\forall t \in [0, T]$ by,

$$ z_1(t) = \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t \sqrt{r(s)} \, ds \right) \quad \text{and} \quad z_2(t) = \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( - \int_0^t \sqrt{r(s)} \, ds \right), $$

where $|\delta_1(t)| \le \exp\left(\frac{1}{2} \int_0^t |\varphi(s)| ds\right) - 1$ and $|\delta_2(t)| \le \exp\left(-\frac{1}{2} \int_t^T |\varphi(s)| ds\right) - 1$. If in addition $\int_0^{+\infty} |\varphi(s)| ds < +\infty$, then the results above also hold for $T = +\infty$.

**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 (14). The LG approximations $\hat{z}_1$ and $\hat{z}_2$ are obtained by neglecting the unknown functions $\delta_1$ and $\delta_2$ in $z_1$ and $z_2$. The theorem gives a non-asymptotic bound for the errors $|z_1(t) - \hat{z}_1(t)|$ and $|z_2(t) - \hat{z}_2(t)|$, $t \ge 0$.

<footnote>$^6$We express $\psi(\tau)$ via $t$ using the one-to-one correspondence between $\tau$ and $t$ to ease readability.</footnote>
## 4.3 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 $\alpha$ and $\varepsilon$ are three times continuously differentiable, and $\varepsilon_0$ is such that $\forall t \ge 0$, $\varepsilon_0 < \frac{(\beta\lambda)^2}{2|\alpha'(t)|+4\lambda}$.

**Remark 4.4.** The condition on $\varepsilon_0$ in Assumption 3 is only technical, so that $r$ defined below is positive. It can be easily satisfied since $|\alpha'(t)|$ is uniformly bounded. Indeed, $\alpha$ is non-increasing and non-negative (see Section 1.1), from which one can deduce that $\int_0^{+\infty} |\alpha'(s)| ds \le \alpha_0$.

We now rewrite (Q1-VM-DIN-AVD) in the form (14).

**Lemma 4.5.** Suppose that Assumption 3 holds, and let $x$ be the solution of (Q1-VM-DIN-AVD). For all $t \ge 0$, define

$$p(t) = \frac{\alpha(t) + \beta\lambda}{\varepsilon(t)} \quad \text{and} \quad r(t) = \frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)}. \qquad (15)$$

Then, $p$ and $r$ are twice continuously differentiable, $r$ is positive and the function $y$ defined for all $t \ge 0$ by $y(t) = x(t) \exp \left( \int_0^t \frac{p(s)}{2} ds \right)$ is a solution to

$$\ddot{y}(t) - r(t)y(t) = 0, \quad t \ge 0, \qquad (16)$$

with initial condition $(y(0), \dot{y}(0)) = (x_0, \dot{x}_0 + \frac{p(0)}{2}x_0)$.

*Proof.* We first check that for all $t \ge 0$, $r(t)$ is positive. Let $t > 0$,

$$r(t) > 0 \iff \frac{(\alpha(t) + \beta\lambda)^2}{4\varepsilon(t)^2} + \frac{\alpha'(t)}{2\varepsilon(t)} - \frac{(\alpha(t) + \beta\lambda)\varepsilon'(t)}{\varepsilon(t)^2} - \frac{\lambda}{\varepsilon(t)} > 0. \qquad (17)$$

Since $\varepsilon'(t) \le 0$ and $\alpha'(t) \le 0$, one can check that a sufficient condition for (17) to hold is,

$$r(t) > 0 \iff \frac{(\alpha(t) + \beta\lambda)^2}{4} > \left(\frac{|\alpha'(t)|}{2} + \lambda\right) \varepsilon(t) \iff \frac{(\beta\lambda)^2}{2|\alpha'(t)| + 4\lambda} > \varepsilon_0.$$

So under Assumption 3, for all $t \ge 0$, $r(t) > 0$. We then check that $y$ is indeed solution to (16). Let $t > 0$,

$$\dot{y}(t) = \frac{p(t)}{2}x(t) \exp \left( \int_{0}^{t} \frac{p(s)}{2} ds \right) + \dot{x}(t) \exp \left( \int_{0}^{t} \frac{p(s)}{2} ds \right),$$

and

$$\ddot{y}(t) = \exp \left( \int_{0}^{t} \frac{p(s)}{2} ds \right) \left[ \left( \frac{p(t)^2}{4} + \frac{p'(t)}{2} \right) x(t) + p(t) \dot{x}(t) + \ddot{x}(t) \right].$$

Since $x$ solves (Q1-VM-DIN-AVD), it holds that $\ddot{x}(t) = -p(t)\dot{x}(t) - \frac{\lambda}{\varepsilon(t)}x(t)$, so,

$$\ddot{y}(t) = \exp \left( \int_{0}^{t} \frac{p(s)}{2} ds \right) \left( \frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)} \right) x(t) \\ = \left( \frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)} \right) y(t) = r(t)y(t). \quad \square$$
– Since we assumed $r$ to be twice continuously differentiable and positive, $\varphi$ is continuous, so it is integrable except maybe for $t \to +\infty$.

– For the sake of simplicity, the formulation of Theorem 4.2 slightly differs from that in [42], the original formulation can be recovered by a change of variable.

## 4.3 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 4.** The functions $\alpha$ and $\varepsilon$ are three times continuously differentiable, and $\varepsilon_0$ is such that $\forall t \ge 0, \varepsilon_0 < \frac{(\beta\lambda)^2}{2|\alpha'(t)|+4\lambda}$.

**Remark 4.4.** The condition on $\varepsilon_0$ in Assumption 4 is only technical, so that $r$ defined below is positive. It can be easily satisfied since $|\alpha'(t)|$ is uniformly bounded. Indeed, $\alpha$ is non-increasing and non-negative (see Section 1.1), from which one can deduce that $\int_0^{+\infty} |\alpha'(s)| ds \le \alpha_0$.

We now rewrite (Q1-VM-DIN-AVD) in the form (14).

**Lemma 4.5.** Suppose that Assumption 4 holds, and let $x$ be the solution of Q1-VM-DIN-AVD. For all $t \ge 0$, define

$$p(t) = \frac{\alpha(t) + \beta\lambda}{\varepsilon(t)} \quad \text{and} \quad r(t) = \frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)}. \qquad (15)$$

Then, $p$ and $r$ are twice continuously differentiable, $r$ is positive and the function $y$ defined for all $t \ge 0$ by $y(t) = x(t) \exp \left( \int_0^t \frac{p(s)}{2} ds \right)$ is a solution to

$$\ddot{y}(t) - r(t)y(t) = 0, \quad t \ge 0, \qquad (16)$$

with initial condition $(y(0), \dot{y}(0)) = (x_0, \dot{x}_0 + \frac{p(0)}{2}x_0)$.

*Proof.* We first check that for all $t \ge 0$, $r(t)$ is positive. Let $t > 0$,

$$r(t) > 0 \iff \frac{(\alpha(t) + \beta\lambda)^2}{4\varepsilon(t)^2} + \frac{\alpha'(t)}{2\varepsilon(t)} - \frac{(\alpha(t) + \beta\lambda)\varepsilon'(t)}{\varepsilon(t)^2} - \frac{\lambda}{\varepsilon(t)} > 0. \qquad (17)$$

Since $\varepsilon'(t) \le 0$ and $\alpha'(t) \le 0$, one can check that a sufficient condition for (17) to hold is,

$$r(t) > 0 \iff \frac{(\alpha(t) + \beta\lambda)^2}{4} > \left( \frac{|\alpha'(t)|}{2} + \lambda \right) \varepsilon(t) \iff \frac{(\beta\lambda)^2}{2|\alpha'(t)| + 4\lambda} > \varepsilon_0.$$

So under Assumption 4, for all $t \ge 0$, $r(t) > 0$. We then check that $y$ is indeed solution to (16). Let $t > 0$,

$$\dot{y}(t) = \frac{p(t)}{2}x(t) \exp\left(\int_{0}^{t} \frac{p(s)}{2} ds\right) + \dot{x}(t) \exp\left(\int_{0}^{t} \frac{p(s)}{2} ds\right), \quad \text{and,}$$

$$\ddot{y}(t) = \exp\left(\int_{0}^{t} \frac{p(s)}{2} ds\right) \left[ \left( \frac{p(t)^2}{4} + \frac{p'(t)}{2} \right) x(t) + p(t) \dot{x}(t) + \ddot{x}(t) \right].$$
Lemma 4.5 gives a reformulation of (Q1-VM-DIN-AVD) suited to apply Theorem 4.2. To use the theorem for all $t \ge 0$, we need to ensure that $\varphi(t) = \frac{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 $\varepsilon$ and $\alpha$ have first, second and third-order derivatives that are integrable on $[0, +\infty]$. In addition, $\lim_{t \to \infty} \varepsilon(t) = 0$ and $\varepsilon'(t)^2/\varepsilon(t)$ is integrable on $[0, +\infty[$.

**Remark 4.6.** Assumption 4 holds for most decays used in practice, with in particular any polynomial decay of the form $\frac{\varepsilon_0}{(t+1)^a}$ and $\frac{\alpha_0}{(t+1)^b}$, $a \in \mathbb{N} \setminus \{0\}$ and $b \in \mathbb{N}$. Note that $\varepsilon$ and $\alpha$ need not be integrable and $\alpha$ can even be constant.

The next lemma states the integrability of $\varphi$ on $[0, +\infty[$.

**Lemma 4.7.** Under Assumption 3 and 4, $\int_0^{+\infty} |\varphi(s)| ds < +\infty$.

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 \in \mathbb{R}$ such that $x(0) = x_0$, $\dot{x}(0) = \dot{x}_0$ and for all $t \ge 0$, the solution of (Q1-VM-DIN-AVD) is

$$x(t) = A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) ds \right) \\ + B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \frac{\sqrt{\varepsilon(t)}}{\sqrt{\alpha(t) + \beta\lambda}} \exp \left( \int_0^t -\frac{\alpha(s) + \beta\lambda}{\varepsilon(s)} + \frac{\lambda}{\alpha(s) + \beta\lambda} + \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) ds \right), \quad (18)$$

where for all $t \ge 0$,

$$|\delta_1(t)| \le \exp\left(\frac{1}{2} \int_0^t |\varphi(s)| ds\right) - 1 \quad \text{and} \quad |\delta_2(t)| \le \exp\left(-\frac{1}{2} \int_t^{+\infty} |\varphi(s)| ds\right) - 1. \quad (19)$$

Thanks to the bounds (19), 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 (15). Let us also define $y(t) \stackrel{\text{def}}{=} x(t) \exp \left( \int_0^t \frac{p(s)}{2} ds \right)$. According to Lemma 4.5, $r$ is positive and $y$ is solution to (16).

Then, from Lemma 4.7, $\int_t^T |\varphi(s)| ds < +\infty$, so we can apply Theorem 4.2 to $y$ on $[0, +\infty[$. Therefore, there exists $A, B \in \mathbb{R}$, such that $\forall t \ge 0$,

$$y(t) = A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t \sqrt{r(s)} ds \right) + B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\sqrt{r(s)} ds \right),$$

where $A$ and $B$ are determined by the initial conditions, and $\delta_1, \delta_2$ are such that (19) holds.

Going back to $x(t) = y(t) \exp \left( \int_0^t -\frac{p(s)}{2} ds \right)$, we obtain that for all $t \ge 0$,

$$x(t) = A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\frac{p(s)}{2} + \sqrt{r(s)} ds \right) + B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\frac{p(s)}{2} - \sqrt{r(s)} ds \right). \quad (20)$$
Since $x$ solves (Q1-VM-DIN-AVD), it holds that $\ddot{x}(t) = -p(t)\dot{x}(t) - \frac{\lambda}{\varepsilon(t)}x(t)$, so,

$$ 
\begin{aligned}
\ddot{y}(t) &= \exp\left(\int_{0}^{t} \frac{p(s)}{2} ds\right) \left(\frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)}\right) x(t) \\
&= \left(\frac{p(t)^2}{4} + \frac{p'(t)}{2} - \frac{\lambda}{\varepsilon(t)}\right) y(t) = r(t)y(t).
\end{aligned}
$$

<figure><img src="image_3.png" /></figure>

Lemma 4.5 gives a reformulation of (Q1-VM-DIN-AVD) suited to apply Theorem 4.2. To use the theorem for all $t \ge 0$, we need to ensure that $\varphi(t) = \frac{4r(t)r''(t)-5r'(t)^2}{16r(t)^{5/2}}$ is integrable. To this aim we make the following assumption.

**Assumption 5.** The functions $\varepsilon$ and $\alpha$ have first, second and third-order derivatives that are integrable on $[0, +\infty[$. In addition, $\lim_{t \to \infty} \varepsilon(t) = 0$ and $\varepsilon'(t)^2/\varepsilon(t)$ is integrable on $[0, +\infty[$.

**Remark 4.6.** Assumption 5 holds for most decays used in practice, with in particular any polynomial decay of the form $\frac{\varepsilon_0}{(t+1)^a}$ and $\frac{\alpha_0}{(t+1)^b}$, $a \in \mathbb{N} \setminus \{0\}$ and $b \in \mathbb{N}$. Note that $\varepsilon$ and $\alpha$ need not be integrable and $\alpha$ can even be constant.

The next lemma states the integrability of $\varphi$ on $[0, +\infty[$.

**Lemma 4.7.** Under Assumption 4 and 5, $\int_{0}^{+\infty} |\varphi(s)| ds < +\infty$.

The proof of this lemma involves long computations and is thus postponed to Appendix D. 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 4 and 5 hold. There exists $A, B \in \mathbb{R}$ such that $x(0) = x_0$, $\dot{x}(0) = \dot{x}_0$ and for all $t \ge 0$, the solution of (Q1-VM-DIN-AVD) is

$$ 
\begin{aligned}
x(t) = &A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) ds \right) \\
&+ B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \frac{\sqrt{\varepsilon(t)}}{\sqrt{\alpha(t) + \beta\lambda}} \\
&\exp \left( \int_0^t -\frac{\alpha(s) + \beta\lambda}{\varepsilon(s)} + \frac{\lambda}{\alpha(s) + \beta\lambda} + \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) ds \right),
\end{aligned}
$$

where for all $t \ge 0$,

$$ |\delta_1(t)| \le \exp\left(\frac{1}{2} \int_0^t |\varphi(s)| ds\right) - 1 \quad \text{and} \quad |\delta_2(t)| \le \exp\left(-\frac{1}{2} \int_t^{+\infty} |\varphi(s)| ds\right) - 1. \quad (19) $$

Thanks to the bounds (19), 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 (15). Let us also define $y(t) \stackrel{\text{def}}{=} x(t) \exp\left(\int_0^t \frac{p(s)}{2} ds\right)$. According to Lemma 4.5, $r$ is positive and $y$ is solution to (16).
It now remains to expand the terms in the two exponentials in (20) in order to obtain (18). To this aim, we approximate $\sqrt{r(s)}$, let $s \ge 0$,

$$
\begin{align*}
\sqrt{r(s)} &= \frac{p(s)}{2} \sqrt{1 + \frac{2p'(s)}{p(s)^2} - \frac{4\lambda}{\varepsilon(s)p(s)^2}} \\
&= \frac{p(s)}{2} \left( 1 + \frac{p'(s)}{p(s)^2} - \frac{2\lambda}{\varepsilon(s)p(s)^2} - \frac{1}{8} \left( \frac{2p'(s)}{p(s)^2} - \frac{4\lambda}{\varepsilon(s)p(s)^2} \right)^2 + o(\varepsilon(s)^2) \right) \\
&= \frac{p(s)}{2} + \frac{p'(s)}{2p(s)} - \frac{\lambda}{\varepsilon(s)p(s)} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda}{\varepsilon(s)p(s)^{3/2}} \right)^2 + o(\varepsilon(s)) \\
&= \frac{p(s)}{2} + \frac{p'(s)\varepsilon(s)}{2(\alpha(s) + \beta\lambda)} - \frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s)) \\
&= \frac{p(s)}{2} + \frac{\alpha'(s)}{2(\alpha(s) + \beta\lambda)} - \frac{\varepsilon'(s)}{2\varepsilon(s)} - \frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s)) \tag{21}
\end{align*}
$$

Focusing on the first exponential term in (20), we deduce from (21) that for all $t \ge 0$,

$$
\begin{align*}
& \exp \left( \int_0^t -\frac{p(s)}{2} + \sqrt{r(s)} \, ds \right) \\
&= \exp \left( \int_0^t \frac{\alpha'(s)}{\alpha(s) + \beta\lambda} \frac{\varepsilon'(s)}{2\varepsilon(s)} - \frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s)) \, ds \right) \\
&= \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\alpha_0 + \beta\lambda}} \frac{\sqrt{\varepsilon_0}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t \frac{-\lambda}{\alpha(s) + \beta\lambda} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s)) \, ds \right) \\
&= \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\alpha_0 + \beta\lambda}} \frac{\sqrt{\varepsilon_0}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t \frac{-\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right),
\end{align*}
$$

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(\int_0^t -\frac{p(s)}{2} - \sqrt{r(s)} \, ds)$, and up to redefining A and B so as to encompass all the constants, we obtain (18) and the result is proved. $\square$

## 4.4 Comparison of $x$ with $x_{LM}$ and $x_N$

We now have an expression for $x$ which is almost explicit: we do not know $\delta_1$ and $\delta_2$ in closed form, but they are uniformly bounded (by Lemma 4.7). We will now compare the asymptotic behavior of (18) 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 $\sim_{+\infty}$ denotes the asymptotic equivalence$^7$ between two functions as $t \to \infty$.

<footnote>$^7$Two real-valued functions $g_1$ and $g_2$ are asymptotically equivalent in $+\infty$ if and only if $\lim_{t \to \infty} \frac{g_1(t)}{g_2(t)} = 1$.</footnote>
Then, from Lemma 4.7, $\int_t^T |\varphi(s)| ds < +\infty$, so we can apply Theorem 4.2 to $y$ on $[0, +\infty[$. Therefore, there exists $A, B \in \mathbb{R}$, such that $\forall t \ge 0$,

$$y(t) = A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t \sqrt{r(s)} \, ds \right) + B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\sqrt{r(s)} \, ds \right),$$

where $A$ and $B$ are determined by the initial conditions, and $\delta_1, \delta_2$ are such that (19) holds. Going back to $x(t) = y(t) \exp \left( \int_0^t -\frac{p(s)}{2} ds \right)$, we obtain that for all $t \ge 0$,

$$x(t) = A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\frac{p(s)}{2} ds + \sqrt{r(s)} ds \right) + B \frac{1 + \delta_2(t)}{r(t)^{1/4}} \exp \left( \int_0^t -\frac{p(s)}{2} ds - \sqrt{r(s)} ds \right). \quad (20)$$

It now remains to expand the terms in the two exponentials in (20) in order to obtain (18). To this aim, we approximate $\sqrt{r(s)}$, let $s \ge 0$,

$$\begin{aligned}
\sqrt{r(s)} &= \frac{p(s)}{2} \sqrt{1 + \frac{2p'(s)}{p(s)^2} - \frac{4\lambda}{\varepsilon(s)p(s)^2}} \\
&= \frac{p(s)}{2} \left( 1 + \frac{p'(s)}{p(s)^2} - \frac{2\lambda}{\varepsilon(s)p(s)^2} - \frac{1}{8} \left( \frac{2p'(s)}{p(s)^2} - \frac{4\lambda}{\varepsilon(s)p(s)^2} \right)^2 + o(\varepsilon(s)^2) \right) \\
&= \frac{p(s)}{2} + \frac{p'(s)}{2p(s)} - \frac{\lambda}{\varepsilon(s)p(s)} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda}{\varepsilon(s)p(s)^{3/2}} \right)^2 + o(\varepsilon(s)) \quad (21) \\
&= \frac{p(s)}{2} + \frac{p'(s)\varepsilon(s)}{2(\alpha(s) + \beta\lambda)} - \frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s)) \\
&= \frac{p(s)}{2} + \frac{\alpha'(s)}{\alpha(s) + \beta\lambda} - \frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\varepsilon'(s)}{2\varepsilon(s)} - \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 + o(\varepsilon(s))
\end{aligned}$$

To ease the readability, we denote $h(t) = \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s)+\beta\lambda)^{3/2}} \right)^2$. Focusing on the first exponential term in (20), we deduce from (21) that for all $t \ge 0$,

$$\begin{aligned}
& \exp \left( \int_0^t -\frac{p(s)}{2} ds + \sqrt{r(s)} ds \right) \\
&= \exp \left( \int_0^t \frac{\alpha'(s)}{\alpha(s) + \beta\lambda} ds - \frac{\varepsilon'(s)}{2\varepsilon(s)} ds - \frac{\lambda}{\alpha(s) + \beta\lambda} ds - \frac{1}{16} h(t) + o(\varepsilon(s)) ds \right) \\
&= \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\alpha_0 + \beta\lambda}} \frac{\sqrt{\varepsilon_0}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t \frac{-\lambda}{\alpha(s) + \beta\lambda} ds - \frac{1}{16} h(t) + o(\varepsilon(s)) ds \right) \\
&= \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\alpha_0 + \beta\lambda}} \frac{\sqrt{\varepsilon_0}}{\sqrt{\varepsilon(t)}} \exp \left( \int_0^t \frac{-\lambda}{\alpha(s) + \beta\lambda} ds - \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) ds \right),
\end{aligned}$$

where the last line relies on further computations postponed to Lemma D.1 in Appendix D. Performing the exact same type of computations on $\exp \left( \int_0^t -\frac{p(s)}{2} ds - \sqrt{r(s)} ds \right)$, and up to redefining $A$ and $B$ so as to encompass all the constants, we obtain (18) and the result is proved. $\square$
**Theorem 4.9.** Let $x$ be the solution of (Q1-VM-DIN-AVD), given in (18), and $x_{LM}$ and $x_N$ whose closed forms are stated in (10). Under Assumptions 3 and 4, there exists $C > 0$ such that the following asymptotic equivalences hold:

$$ x(t) \sim_{+\infty} x_{LM}(t)C \exp \left( \int_{0}^{t} \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right), \quad \text{and} \\ x(t) \sim_{+\infty} x_N(t)C \exp \left( \int_{0}^{t} \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} - \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right). \quad (22) $$

As a consequence, the convergence of $x$ to $x^*$ is:

(i) Faster than that of $x_{LM}$ if $\varepsilon$ is non-integrable and as fast otherwise.

(ii) Slower than that of $x_N$ if $\alpha$ is non-integrable and as fast if $\alpha$ is integrable, in the case where $\forall t \ge 0, \alpha(t) > \varepsilon(t)$.

(iii) Faster than that of $x_N$ if $\varepsilon$ is non-integrable and as fast if $\varepsilon$ is integrable, in the case where $\forall t \ge 0, \alpha(t) < \varepsilon(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 $\varepsilon$ is non-integrable but vanishing. However, when the variable mass $\varepsilon$ is larger than $\alpha$, the dynamics is closer to the one of the Newton method, and can actually be an accelerated Newton dynamics, again for non-integrable $\varepsilon$. This is analogous to the necessary condition that $\alpha$ must be non-integrable in order to accelerate first-order methods in convex optimization (see [3]). 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 (18). We now analyze the two terms in (18).

First, we know from Theorem 4.8 that $\delta_1(0) = 0$ and $\lim_{t \to +\infty} \delta_2(t) = 0$. In addition, by Lemma 4.7, $\delta_1$ and $\delta_2$ are uniformly bounded by some positive constant. Then $r(t)^{-1/4}$ decays asymptotically like $\sqrt{\varepsilon(t)}$ and $\alpha$ is bounded. So $A \frac{1+\delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t)+\beta\lambda}}{\sqrt{\varepsilon(t)}}$ is asymptotically equivalent to some constant $c_1 \in \mathbb{R}$ as $t \to +\infty$. Similarly, $B \frac{1+\delta_2(t)}{r(t)^{1/4}} \frac{\sqrt{\varepsilon(t)}}{\sqrt{\alpha(t)+\beta\lambda}}$ is equivalent to $c_2\varepsilon(t)$, with $c_2 \in \mathbb{R}$.

We now analyze the “exponential factors” in (18). On the one hand, $\frac{\lambda}{\alpha(s)+\beta\lambda} + \frac{\lambda^2\varepsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\varepsilon(s))$ converges to $\frac{1}{\beta}$ as $s \to \infty$, while $\frac{\alpha(s)+\beta\lambda}{\varepsilon(s)}$ diverges to $+\infty$. Therefore, we deduce that,

$$ \exp \left( \int_{0}^{t} -\frac{\alpha(s) + \beta\lambda}{\varepsilon(s)} + \frac{\lambda}{\alpha(s) + \beta\lambda} + \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right) = o \left( \exp \left( \int_{0}^{t} -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right) \right). $$
## 4.4 Comparison of $x$ with $x_{LM}$ and $x_N$

We now have an expression for $x$ which is almost explicit: we do not know $\delta_1$ and $\delta_2$ in closed form, but they are uniformly bounded (by Lemma 4.7). We will now compare the asymptotic behavior of (18) 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 $\sim_{+\infty}$ denotes the asymptotic equivalence<fnref n="7" /> between two functions as $t \to \infty$.

**Theorem 4.9.** Let $x$ be the solution of (Q1-VM-DIN-AVD), given in (18), and $x_{LM}$ and $x_N$ whose closed forms are stated in (10). Under Assumptions 4 and 5, there exists $C > 0$ such that the following asymptotic equivalences hold:

$$ 
\begin{aligned}
x(t) &\sim_{+\infty} x_{LM}(t)C \exp \left( \int_0^t \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right), \quad \text{and} \\
x(t) &\sim_{+\infty} x_N(t)C \exp \left( \int_0^t \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} - \frac{\lambda^2 \varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds \right).
\end{aligned}
\quad (22) $$

As a consequence, the convergence of $x$ to $x^*$ is:

(i) Faster than that of $x_{LM}$ if $\varepsilon$ is non-integrable and as fast otherwise.

(ii) Slower than that of $x_N$ if $\alpha$ is non-integrable and as fast if $\alpha$ is integrable, in the case where $\forall t \ge 0, \alpha(t) > \varepsilon(t)$.

(iii) Faster than that of $x_N$ if $\varepsilon$ is non-integrable and as fast if $\varepsilon$ is integrable, in the case where $\forall t \ge 0, \alpha(t) < \varepsilon(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 controlling (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 $\varepsilon$ is non-integrable but vanishing. However, when the variable mass $\varepsilon$ is larger than $\alpha$, the dynamics is closer to the one of the Newton method, and can actually be an accelerated Newton dynamics, again for non-integrable $\varepsilon$. This is analogous to the necessary condition that $\alpha$ must be non-integrable in order to accelerate first-order methods in convex optimization (see [4]). We conclude this section by proving Theorem 4.9.

**Proof of Theorem 4.9.** Thanks to Assumptions 4 and 5, Theorem 4.8 tells us that $x$ has the form (18). We now analyze the two terms in (18).

First, we know from Theorem 4.8 that $\delta_1(0) = 0$ and $\lim_{t \to +\infty} \delta_2(t) = 0$. In addition, by Lemma 4.7, $\delta_1$ and $\delta_2$ are uniformly bounded by some positive constant. Then $r(t)^{-1/4}$ decays asymptotically like

<footnote n="7">⁷Two real-valued functions $g_1$ and $g_2$ are asymptotically equivalent in $+\infty$ if and only if $\lim_{t \to +\infty} \frac{g_1(t)}{g_2(t)} = 1$.</footnote>
As a consequence, the second term in (18) will decrease to 0 faster than the first-one (let alone the additional $\epsilon(t)$ decay that we have just discussed). The asymptotic behavior of $x$ will thus be governed by the first term in (18).

Let us now focus on the first term in (18). Observe that $\exp(\int_0^t -\frac{\lambda}{\alpha(s)+\beta\lambda} ds)$ is exactly the decay of $x_{LM}$ in (10). Thus, we have proved that there exists $C > 0$, such that the following asymptotic equivalence holds,

$$A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\epsilon(t)}} \exp \left( \int_0^t -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2 \epsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\epsilon(s)) ds \right) \\ \sim_{+\infty} x_{LM}(t) C \exp \left( \int_0^t -\frac{\lambda^2 \epsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\epsilon(s)) ds \right),$$

which proves the first part of (22). The second equivalence in (22) is obtained using the following identity,

$$\int_0^t -\frac{\lambda}{\alpha(s) + \beta\lambda} ds = \int_0^t -\frac{1}{\beta} + \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} ds = -\frac{t}{\beta} + \int_0^t \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} ds \quad (23)$$

and $e^{-t/\beta}$ is precisely the rate at which $x_N$ decreases. So (22) holds.

It finally remains to deduce the conclusions of the theorem from (22).

– Regarding the comparison with $x_{LM}$, the integral $\int_0^t -\frac{\lambda^2\epsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\epsilon(s)) ds$ converges if and only if $\epsilon$ is integrable on $[0, +\infty[$, and diverges to $-\infty$ when $\epsilon$ is not. So $x$ converges to 0 at least as fast as $x_{LM}$ and faster when $\epsilon$ is not integrable.

– As for the comparison with $x_N$, if $\alpha(s) > \epsilon(s) \ge 0$ for all $s \ge 0$, then the integral $\int_0^t \frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)} - \frac{\lambda^2\epsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\epsilon(s)) ds$ is convergent in $+\infty$ if and only if $\alpha$ is integrable and diverges to $+\infty$ when $\alpha$ is non-integrable. So when $\alpha$ is integrable, the speed of convergence of $x$ is the same as that of $x_N$. When $\alpha$ is not integrable, the convergence to 0 is slower but still holds. Indeed, for all $s \ge 0$ $\frac{\alpha(s)}{\beta(\alpha(s)+\frac{1}{\beta})} < \frac{\alpha(s)}{\beta\alpha(s)} = \frac{1}{\beta}$. Thus for all $t > 0$, $-\frac{t}{\beta} + \int_0^t \frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)} ds < 0$.

– Finally, the comparison with $x_N$ in the case $\epsilon(s) > \alpha(s)$ is exactly the same as the comparison with $x_{LM}$ using (23). $\square$

# 5 Numerical Experiments

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

## 5.1 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 $\gamma = 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.
$\sqrt{\varepsilon(t)}$ and $\alpha$ is bounded. So $A \frac{1+\delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t)+\beta\lambda}}{\sqrt{\varepsilon(t)}}$ is asymptotically equivalent to some constant $c_1 \in \mathbb{R}$ as $t \to +\infty$. Similarly, the factor $B \frac{1+\delta_2(t)}{r(t)^{1/4}} \frac{\sqrt{\varepsilon(t)}}{\sqrt{\alpha(t)+\beta\lambda}}$ is equivalent to $c_2\varepsilon(t)$, with $c_2 \in \mathbb{R}$.

We now analyze the “exponential factors” in (18). On the one hand, $\frac{\lambda}{\alpha(s)+\beta\lambda} + \frac{\lambda^2\varepsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\varepsilon(s))$ converges to $\frac{1}{\beta}$ as $s \to \infty$, while $\frac{\alpha(s)+\beta\lambda}{\varepsilon(s)}$ diverges to $+\infty$. Therefore, we deduce that,

$$ \exp\left(\int_{0}^{t} -\frac{\alpha(s) + \beta\lambda}{\varepsilon(s)} + \frac{\lambda}{\alpha(s) + \beta\lambda} + \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds\right) = o\left(\exp\left(\int_{0}^{t} -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds\right)\right). $$

As a consequence, the second term in (18) will decrease to 0 faster than the first-one (let alone the additional $\varepsilon(t)$ decay that we have just discussed). The asymptotic behavior of $x$ will thus be governed by the first term in (18).

Let us now focus on the first term in (18). Observe that $\exp\left(\int_{0}^{t} -\frac{\lambda}{\alpha(s)+\beta\lambda} \, ds\right)$ is exactly the decay of $x_{LM}$ in (10). Thus, we have proved that there exists $C > 0$, such that the following asymptotic equivalence holds,

$$ A \frac{1 + \delta_1(t)}{r(t)^{1/4}} \frac{\sqrt{\alpha(t) + \beta\lambda}}{\sqrt{\varepsilon(t)}} \exp\left(\int_0^t -\frac{\lambda}{\alpha(s) + \beta\lambda} - \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds\right) \sim_{+\infty} x_{LM}(t)C \exp\left(\int_0^t -\frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)) \, ds\right), $$

which proves the first part of (22). The second equivalence in (22) is obtained using the following identity,

$$ \int_{0}^{t} -\frac{\lambda}{\alpha(s) + \beta\lambda} \, ds = \int_{0}^{t} -\frac{1}{\beta} + \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} \, ds = -\frac{t}{\beta} + \int_{0}^{t} \frac{\alpha(s)}{\beta(\alpha(s) + \beta\lambda)} \, ds \quad (23) $$

and $e^{-t/\beta}$ is precisely the rate at which $x_N$ decreases. So (22) holds.

It finally remains to deduce the conclusions of the theorem from (22).

– Regarding the comparison with $x_{LM}$, the integral $\int_{0}^{t} -\frac{\lambda^2\varepsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\varepsilon(s)) \, ds$ converges if and only if $\varepsilon$ is integrable on $[0, +\infty[$, and diverges to $-\infty$ when $\varepsilon$ is not. So $x$ converges to 0 at least as fast as $x_{LM}$ and faster when $\varepsilon$ is not integrable.

– As for the comparison with $x_N$, if $\alpha(s) > \varepsilon(s) \ge 0$ for all $s \ge 0$, then the integral $\int_{0}^{t} \frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)} - \frac{\lambda^2\varepsilon(s)}{(\alpha(s)+\beta\lambda)^3} + o(\varepsilon(s)) \, ds$ is convergent in $+\infty$ if and only if $\alpha$ is integrable and diverges to $+\infty$ when $\alpha$ is non-integrable. So when $\alpha$ is integrable, the speed of convergence of $x$ is the same as that of $x_N$. When $\alpha$ is not integrable, the convergence to 0 is slower but still holds. Indeed, for all $s \ge 0$ $\frac{\alpha(s)}{\beta(\alpha(s)+\frac{1}{\beta})} < \frac{\alpha(s)}{\beta\alpha(s)} = \frac{1}{\beta}$. Thus for all $t > 0$, $-\frac{t}{\beta} + \int_{0}^{t} \frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)} \, ds < 0$.

– Finally, the comparison with $x_N$ in the case $\varepsilon(s) > \alpha(s)$ is exactly the same as the comparison with $x_{LM}$ using (23). $\square$
<figure><img src="image_6.png" /><figcaption>Figure 2: 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} + \frac{1}{2}\|Ax\|^2$. Left figures: distance $\|x(t) - x_N(t)\|$ versus time $t$, each curve corresponds to a different choice of $\varepsilon$; middle figures: distance $\|x(t) - x_{LM}(t)\|$, again for several $\varepsilon$. 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 $\varepsilon$. 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 $\alpha$. The theoretical bounds from Theorem 3.8 are only displayed on Figure 4 below, for the sake of readability.</figcaption></figure>

## 5.2 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 and 4 (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)\| \le \|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 non-asymptotic and this is highlighted by the experiments as we now explain.

Looking at the left and middle plots of Figures 2, 3 and 4, 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
**Remark 4.10.** As discussed in the beginning of this section, the main motivation for studying quadratic functions is that strongly convex functions with Lipschitz continuous gradient can be tightly both upper- and lower-bounded by quadratic functions around their minimizers. The intuitive attempt to extend our analysis to this class of functions consists therefore in applying Theorem 4.8 to the upper and lower quadratic bounds. Such analysis is hence only possible for initial conditions close-enough to the minimizer, and even in that case we did not manage to derive precise-enough asymptotic estimates for the solution of (VM-DIN-AVD) using this strategy. Whether this is possible or not remains an open question.

# 5 Numerical Experiments

<figure><img src="image_4.png" /><figcaption>Figure 2: 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} + \frac{1}{2}\|Ax\|^2$. Left figures: distance $\|x(t) - x_N(t)\|$ versus time $t$, each curve corresponds to a different choice of $\varepsilon$; middle figures: distance $\|x(t) - x_{LM}(t)\|$, again for several $\varepsilon$. 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 $\varepsilon$. 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 $\alpha$. The theoretical bounds from Theorem 3.8 are only displayed on Figure 4 below, for the sake of readability.</figcaption></figure>

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

## 5.1 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
<figure><img src="image_7.png" /><figcaption>Figure 3: Similar experiment and figures as those described in Figure 2, but for the function $f(x) = \log (\sum_{i=1}^{n} e^{x_i} + e^{-x_i}) + \frac{1}{2} \|Ax\|^2$.</figcaption></figure>

of our theorems appears on the left of Figures 3 and 4: the blue and green curves, corresponding to slowly decaying choices of $\epsilon$, 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 $\epsilon > \alpha$, $x$ may be closer to $x_N$ and when $\epsilon \le \alpha$, $x$ would rather be closer to $x_{LM}$. This is more noticeable on the top rows of the figures, where $\alpha$ 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 $\epsilon$. 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 $\epsilon$ 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 $\epsilon(t) = \epsilon_0/t$ combined either with $\alpha(t) = \alpha_0/t$ or $\alpha(t) = \alpha_0/t^2$ seems to always yield fast convergence on these experiments (and sometimes the fastest of all dynamics).

## 5.3 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) = \frac{1}{2}\|Ay\|^2$, $y \in \mathbb{R}^n$, and for several choices of $\epsilon$ and $\alpha$. The results in Figure 5 exactly match the expected behavior summarized in Table 1. Indeed, looking first at the right-hand side of Figure 5, $x$ is as fast
<figure><img src="image_5.png" /><figcaption>Figure 3: Similar experiment and figures as those described in Figure 2, but for the function $f(x) = \log (\sum_{i=1}^{n} e^{x_i} + e^{-x_i}) + \frac{1}{2} \|Ax\|^2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Similar experiment and figures as those described in Figure 2, but for the function $f(x) = \sum_{i=1}^{n} x_i^{50} + \frac{1}{2} \|Ax\|^2$. The thin “dash dotted” curves represent approximations of the theoretical bounds from Theorem 3.8 for each choice of $(\epsilon, \alpha)$ considered.</figcaption></figure>

schemes with small$^8$ step-sizes $\gamma = 10^{-1}$. All ODEs are initialized at $x_0 \in \mathbb{R}^n$, where each coordinate of $x_0$ belongs to $\{-1, 1\}$. For all functions this means that $x_0$ is approximately at distance $\sqrt{n}$ of $x^*$, hence not very close to $x^*$. We then approximated the ODEs on time intervals that are long-enough to observe both non-asymptotic and asymptotic convergence behaviors. We used Euler semi-explicit

<footnote>$^8$In practice Newton’s method would converge faster when used with $\gamma = \beta = 1$, here we take a smaller step-size $\gamma$ to more accurately approximate the solutions of the ODEs.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 4: Similar experiment and figures as those described in Figure 2, but for the function $f(x) = \|x\|^{50} + \frac{1}{2}\|Ax\|^{2}$. The thin “dash dotted” curves represent the theoretical bounds from Theorem 3.8 for each choice of $(\varepsilon, \alpha)$ considered.</figcaption></figure>

as the corresponding<fnref n="8" /> $x_{LM}$ when $\varepsilon$ is not integrable and regardless of $\alpha$, and $x$ is faster when $\varepsilon$ is non-integrable. Then on the left-hand side, when comparing to $x_N$, $x$ is slower in settings where $\alpha$ is larger than $\varepsilon$ and non-integrable (red curves), or almost as fast when $\alpha$ is integrable (pink curve). However, acceleration w.r.t. to $x_N$ is indeed achieved for non-integrable $\varepsilon$ regardless of $\alpha$ (first-two blue curves), and the rate is the same as that of $x_N$ when $\varepsilon$ is integrable (third blue curve).

# 6 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 $\varepsilon$ and $\alpha$, 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 $\varepsilon$ and $\alpha$ 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.

<footnote n="8">⁸That is, the solution of (LM) for the same $\alpha$ as that considered for (VM-DIN-AVD).</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5: Numerical validation of Theorem 4.9: distance to the optimum $x^*$ as a function of time on a quadratic function $f(x) = \frac{1}{2}\|Ax\|^2$. Left: speed comparison w.r.t. (CN) for several choices of $\varepsilon$ and $\alpha$. Right: Comparison with LM for $\alpha$ integrable or not and several choices of $\varepsilon$. Shades of blue represent cases where $\varepsilon(t) > \alpha(t)$ while shades of red represent the opposite setting.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Same experiments as in Figure 5 but in the setting $\gamma = \beta = 1$ (so that Newton's method converges in one iteration), for large values of $\varepsilon_0$ (left figures) and small $\varepsilon_0$ (right). The use of a large step-size $\gamma$ makes this setting out of the scope of our theoretical results. Due to fast convergence, different behaviors can be observed for large values of $\varepsilon_0$ and smaller ones.</figcaption></figure>

schemes obtained by solving linear systems, for the sake of stability. The resulting algorithms are detailed in Appendix E.

## 5.2 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 and 4 corresponds to a different strongly convex function, specified below its corresponding figure. 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)\| \le \|x(t) - x^*\| + \|x_N(t) - x^*\|$, so if both $x$ and $x_N$ converge linearly, then the bounds of Theorems 3.1
<figure><img src="image_9.png" /><figcaption>Figure 5: Numerical validation of Theorem 4.9: distance to the optimum $x^*$ as a function of time. On a quadratic function $f(x) = \frac{1}{2}\|Ax\|^2$. Left: speed comparison w.r.t. (CN) for several choices of $\varepsilon$ and $\alpha$. Right: Comparison with LM in for $\alpha$ integrable or not and several choices of $\varepsilon$. Shades of blue represent cases where $\varepsilon(t) > \alpha(t)$ while shades of red represent the opposite setting.</figcaption></figure>

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 $(\varepsilon, \alpha)$ 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 $\varepsilon$ 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.

## 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].
and 3.8 need not be asymptotically tight. That being said, the strength of these results is to be non-asymptotic and this is highlighted by the experiments as we now explain.

The left and middle plots of Figures 2, 3 and 4 are consistent with Theorems 3.1 and 3.8, 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 of our theorems appears on the left of Figures 3 and 4: the blue and green curves, corresponding to slowly decaying choices of $\epsilon$, 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 $\epsilon > \alpha$, $x$ may be closer to $x_N$ and when $\epsilon \le \alpha$, $x$ would rather be closer to $x_{LM}$. This is more noticeable on the top rows of the figures, where $\alpha$ 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 $\epsilon$. 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 $\epsilon$ might result in faster convergence for $x$ than fast decays (and also faster convergence than $x_{LM}$), notably on Figure 2. We also note that $\epsilon(t) = \epsilon_0/t$ combined either with $\alpha(t) = \alpha_0/t$ or $\alpha(t) = \alpha_0/t^2$ seems to very often yield fast convergence in these experiments.

## 5.3 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) = \frac{1}{2}||Ay||^2$, $y \in \mathbb{R}^n$, and for several choices of $\epsilon$ and $\alpha$. The results in Figure 5 exactly match the expected behavior summarized in Table 1. Indeed, looking first at the right-hand side of Figure 5, $x$ is as fast as the corresponding<fnref n="9" /> $x_{LM}$ when $\epsilon$ is not integrable and regardless of $\alpha$, and $x$ is faster when $\epsilon$ is non-integrable. Then on the left-hand side, when comparing to $x_N$, $x$ is slower in settings where $\alpha$ is larger than $\epsilon$ and non-integrable (red curves), or almost as fast when $\alpha$ is integrable (pink curve). However, acceleration w.r.t. to $x_N$ is indeed achieved for non-integrable $\epsilon$ regardless of $\alpha$ (first-two blue curves), and the rate is the same as that of $x_N$ when $\epsilon$ is integrable (third blue curve). Finally, on Figure 6 we empirically investigate the case of large step-sizes $\gamma$ and take $\gamma = \beta = 1$ (optimal for Newton's method on quadratics). Due to large step-sizes our theoretical results do not hold and the choice of $\epsilon_0$ matters more. This evidences a trade-off between $\epsilon_0$ and $\gamma$ in the case of discrete algorithms which is not present in the analysis of ODEs and is worth investigating in a future work.

# 6 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 $\epsilon$ and $\alpha$, 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 $\epsilon$ and $\alpha$ 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

<footnote n="9">⁹That is, the solution of (LM) for the same $\alpha$ as that considered for (VM-DIN-AVD).</footnote>
# 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,

$$ \frac{d}{dt}[\varepsilon(t)\dot{x}(t)] + \beta \frac{d}{dt}\nabla f(x(t)) + \alpha(t)\dot{x}(t) - \varepsilon'(t)\dot{x}(t) + \nabla f(x(t)) = 0, \quad t \ge 0. \quad (24) $$

We then integrate (24) for all $t \ge 0$,

$$ \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) - \varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) + \int_0^t (\alpha(s) - \varepsilon'(s))\dot{x}(s) + \nabla f(x(s)) ds = 0. \quad (25) $$

For all $t \ge 0$, we define the variable,

$$ z(t) = \int_{0}^{t} (\alpha(s) - \varepsilon'(s))\dot{x}(s) + \nabla f(x(s)) ds - \varepsilon_{0}\dot{x}_{0} - \beta\nabla f(x_{0}). $$

We differentiate $z$, for all $t > 0$, $\dot{z}(t) = (\alpha(t) - \varepsilon'(t))\dot{x}(t) + \nabla f(x(t))$, so that we can rewrite (25) as,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + z(t) = 0 \\ \dot{z}(t) - (\alpha(t) - \varepsilon'(t))\dot{x}(t) - \nabla f(x(t)) = 0 \end{cases}, \quad t \ge 0. $$

We substitute the first line in the second-one,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + z(t) = 0 \\ \beta\dot{z}(t) - \beta(\alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t))\dot{x}(t) + z(t) = 0 \end{cases}, \quad t \ge 0. \quad (26) $$

To ease the readability, we recall the notation $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$ from Section 2. Then define for all $t \ge 0$, $y(t) = z(t) - \nu(t)x(t)$, and differentiate, $\dot{y}(t) = \dot{z}(t) - \nu(t)\dot{x}(t) - \nu'(t)x(t)$. We finally rewrite (26) as,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + \nu(t)x(t) + y(t) = 0 \\ \dot{y}(t) + \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t) = 0 \end{cases}. $$

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

$$ y(0) = z(0) - \nu(0)x(0) = -\varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) - (\alpha_0 - \varepsilon'_0 - \frac{1}{\beta}\varepsilon_0)x_0. $$

**Remark A.1.** Notice that the quantity $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$ involved in (gVM-DIN-AVD) also plays a key role in our analysis of Section 3, see e.g., (6). In particular the sign of $\nu(t)$ changes the nature of (VM-DIN-AVD) and is related to Assumption 1.
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 linear systems (this is also true for (CN) and (LM)). Our new understanding on $(\varepsilon, \alpha)$ 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 $\varepsilon$ 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.

## Acknowledgment

C. Castera, J. Fadili and P. Ochs are supported by the ANR-DFG joint project TRINOM-DS under numbers ANR-20-CE92-0037-01 and OC150/5-1. The numerical experiments were made thanks to the development teams of the following libraries: Python [44], Numpy [49] and Matplotlib [32].
## 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) = \frac{\epsilon(t)}{2} \|\dot{x}(t)\|^2 + f(x(t)) - f(x^*)$ is non-negative and decreasing. Thus, it is uniformly bounded on $\mathbb{R}_+$ and the same holds for $t \mapsto f(x(t))$ since for all $t \ge 0$, $U(t) \ge f(x(t))$. Then, $f$ is coercive by assumption, so $x$ is uniformly bounded on $\mathbb{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$, $\dot{y}(t) = -\frac{1}{\beta}y(t) - (\frac{\nu(t)}{\beta} + \nu'(t))x(t)$ so we can use the following integrating factor,

$$y(t) = e^{-\frac{t}{\beta}}y_0 - e^{-\frac{t}{\beta}}\int_0^t \frac{1}{\beta}e^{\frac{s}{\beta}}(\nu(s) + \beta\nu'(s))x(s) \, ds.$$

Using triangle inequalities, for all $t \ge 0$,

$$\|y(t)\| \le e^{-\frac{t}{\beta}} \|y_0\| + \sup_{s \ge 0} \|(\nu(s) + \beta\nu'(s))x(s)\| e^{-\frac{t}{\beta}} \int_0^t \frac{1}{\beta} e^{\frac{s}{\beta}} \, ds \le \|y_0\| + \sup_{s \ge 0} \|(\nu(s) + \beta\nu'(s))x(s)\|. \quad (27)$$

Using the definition of $\epsilon$ and $\alpha$ from Sections 1.1 and 2, observe that $\epsilon, \alpha, \epsilon'$ and $\alpha'$ are all bounded on $\mathbb{R}_+$, and $\epsilon''$ is assumed to be bounded. So $\nu$ and $\nu'$ are bounded, and since we also proved that $x$ is uniformly bounded on $\mathbb{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 $\epsilon$ and $\alpha$ 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 \ge 0$, $x(t)$, $x_N(t)$ and $x_{LM}(t)$ belong to the bounded set $K_0$ defined in that proof. Since $f$ is $\mu$-strongly convex on $K_0$, the proof relies again on bounding difference of gradients, indeed, for all $t \ge 0$,

$$\|x(t) - x_N(t)\| \le \frac{1}{\mu} \|\nabla f(x(t)) - \nabla f(x_N(t))\| \quad \text{and} \quad \|x(t) - x_{LM}(t)\| \le \frac{1}{\mu} \|\nabla f(x(t)) - \nabla f(x_{LM}(t))\|. \quad (28)$$

Recall also that the closed form of $\nabla f(x_N)$ is given in (4).

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

$$\beta \nabla f(x(t)) = \beta e^{-\frac{t}{\beta}} \nabla f(x_0) + e^{-\frac{t}{\beta}} \epsilon_0 \dot{x}_0 - \epsilon(t) \dot{x}(t) + \int_0^t e^{\frac{s-t}{\beta}} \left( \frac{1}{\beta} \epsilon(s) + \epsilon'(s) - \alpha(s) \right) \dot{x}(s) \, ds.$$
# Appendices

## A Equivalent First-order System and 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,

$$ \frac{d}{dt} [\varepsilon(t)\dot{x}(t)] + \beta \frac{d}{dt} \nabla f(x(t)) + \alpha(t)\dot{x}(t) - \varepsilon'(t)\dot{x}(t) + \nabla f(x(t)) = 0, \quad t \ge 0. \quad (24) $$

We then integrate (24) for all $t \ge 0$,

$$ \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) - \varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) + \int_0^t (\alpha(s) - \varepsilon'(s))\dot{x}(s) + \nabla f(x(s)) ds = 0. \quad (25) $$

For all $t \ge 0$, we define the variable,

$$ z(t) = \int_{0}^{t} (\alpha(s) - \varepsilon'(s))\dot{x}(s) + \nabla f(x(s)) ds - \varepsilon_{0}\dot{x}_{0} - \beta\nabla f(x_{0}). $$

We differentiate $z$, for all $t > 0$, $\dot{z}(t) = (\alpha(t) - \varepsilon'(t))\dot{x}(t) + \nabla f(x(t))$, so that we can rewrite (25) as,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + z(t) = 0 \\ \dot{z}(t) - (\alpha(t) - \varepsilon'(t))\dot{x}(t) - \nabla f(x(t)) = 0 \end{cases}, \quad t \ge 0. $$

We substitute the first line in the second-one,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + z(t) = 0 \\ \beta\dot{z}(t) - \beta(\alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t))\dot{x}(t) + z(t) = 0 \end{cases}, \quad t \ge 0. \quad (26) $$

To ease the readability, we recall the notation $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$ from Section 2. Then define for all $t \ge 0$, $y(t) = z(t) - \nu(t)x(t)$, and differentiate, $\dot{y}(t) = \dot{z}(t) - \nu(t)\dot{x}(t) - \nu'(t)x(t)$. We finally rewrite (26) as,

$$ \begin{cases} \varepsilon(t)\dot{x}(t) + \beta\nabla f(x(t)) + \nu(t)x(t) + y(t) = 0 \\ \dot{y}(t) + \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t) = 0 \end{cases}. $$

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

$$ y(0) = z(0) - \nu(0)x(0) = -\varepsilon_0\dot{x}_0 - \beta\nabla f(x_0) - (\alpha_0 - \varepsilon'_0 - \frac{1}{\beta}\varepsilon_0)x_0. $$

**Remark A.1.** Notice that the quantity $\nu(t) = \alpha(t) - \varepsilon'(t) - \frac{1}{\beta}\varepsilon(t)$ involved in (gVM-DIN-AVD) also plays a key role in our analysis of Section 3, see e.g., (6). In particular the sign of $\nu(t)$ changes the nature of (VM-DIN-AVD) and is related to Assumption 2.
Here we do not assume any relation between $\varepsilon$ and $\alpha$, and we thus need to find a more suitable expression for $\nabla f(x(t))$. We first expand the terms in the integral, for all $t \ge 0$,

$$ \beta\nabla f(x(t)) = \beta e^{-\frac{t}{\beta}}\nabla f(x_0) + e^{-\frac{t}{\beta}}\varepsilon_0\dot{x}_0 - \varepsilon(t)\dot{x}(t) + \int_0^t e^{\frac{s-t}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s)\right)\dot{x}(s)\,ds - \int_0^t e^{\frac{s-t}{\beta}}\alpha(s)\dot{x}(s)\,ds. \quad (29) $$

Then, for all $s \ge 0$, we have the identity,

$$ e^{\frac{s}{\beta}}\dot{x}(s) = e^{\frac{s}{\beta}}\dot{x}(s) + \frac{1}{\beta}e^{\frac{s}{\beta}}x(s) - \frac{1}{\beta}e^{\frac{s}{\beta}}x(s) = \frac{d}{ds}(e^{\frac{s}{\beta}}x(s)) - \frac{1}{\beta}e^{\frac{s}{\beta}}x(s), \quad (30) $$

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

$$ \int_0^t e^{\frac{s}{\beta}}\alpha(s)\dot{x}(s)\,ds = \left[\alpha(s)e^{\frac{s}{\beta}}x(s)\right]_0^t - \int_0^t \left(\alpha'(s) + \frac{\alpha(s)}{\beta}\right)e^{\frac{s}{\beta}}x(s)\,ds. $$

Therefore,

$$ e^{-\frac{t}{\beta}} \int_0^t e^{\frac{s}{\beta}} \alpha(s) \dot{x}(s) \, ds = \alpha(t)x(t) - e^{-\frac{t}{\beta}} \alpha_0 x_0 - \int_0^t e^{\frac{s-t}{\beta}} \left( \alpha'(s) + \frac{\alpha(s)}{\beta} \right) x(s) \, ds, \quad (31) $$

and we can substitute in (29),

$$ \beta\nabla f(x(t)) = \beta e^{-\frac{t}{\beta}}\nabla f(x_0) + e^{-\frac{t}{\beta}}\varepsilon_0\dot{x}_0 - \varepsilon(t)\dot{x}(t) + \int_0^t e^{\frac{s-t}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s)\right)\dot{x}(s)\,ds \\ - \alpha(t)x(t) + e^{-\frac{t}{\beta}}\alpha_0 x_0 + \int_0^t e^{\frac{s-t}{\beta}}\left(\alpha'(s) + \frac{\alpha(s)}{\beta}\right)x(s)\,ds. \quad (32) $$

**Uniform boundedness.** In view of exploiting (32), we recall that for all $(\varepsilon, \alpha)$, $x$ is uniformly bounded. So there exists $R > 0$ such that for all $(\varepsilon, \alpha)$, the corresponding solution $x$ of (VM-DIN-AVD) is such that

$$ \sup_{t \ge 0} \|x(t)\| \le R. \quad (33) $$

We are now in position to prove (8).

**Distance from $x$ to $x_N$.** We first gather (4) and (32). For all $t \ge 0$,

$$ \beta\nabla f(x(t)) - \beta\nabla f(x_N(t)) = e^{-\frac{t}{\beta}}\varepsilon_0\dot{x}_0 - \varepsilon(t)\dot{x}(t) + \int_0^t e^{\frac{s-t}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s)\right)\dot{x}(s)\,ds \\ + e^{-\frac{t}{\beta}}\alpha_0 x_0 - \alpha(t)x(t) + \int_0^t e^{\frac{s-t}{\beta}}\left(\alpha'(s) + \frac{\alpha(s)}{\beta}\right)x(s)\,ds. $$

We then use (28) and triangle inequalities,

$$ \beta\mu\|x(t) - x_N(t)\| \le e^{-\frac{t}{\beta}}\varepsilon_0\|\dot{x}_0\| + \varepsilon(t)\|\dot{x}(t)\| + \int_0^t e^{\frac{s-t}{\beta}}\left|\frac{\varepsilon(s)}{\beta} + \varepsilon'(s)\right|\|\dot{x}(s)\|\,ds \\ + e^{-\frac{t}{\beta}}\alpha_0\|x_0\| + \alpha(t)\|x(t)\| + \int_0^t e^{\frac{s-t}{\beta}}\left|\frac{\alpha(s)}{\beta} + \alpha'(s)\right|\|x(s)\|\,ds. \quad (34) $$
## 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) = \frac{\epsilon(t)}{2} \|\dot{x}(t)\|^2 + f(x(t)) - f(x^*)$ is non-negative and decreasing. Thus, it is uniformly bounded on $\mathbb{R}_+$ and the same holds for $t \mapsto f(x(t))$ since for all $t \ge 0$, $U(t) \ge f(x(t))$. Then, $f$ is coercive by assumption, so $x$ is uniformly bounded on $\mathbb{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$, $\dot{y}(t) = -\frac{1}{\beta}y(t) - (\frac{\nu(t)}{\beta} + \nu'(t))x(t)$ so we can use the following integrating factor,

$$y(t) = e^{-\frac{t}{\beta}}y_0 - e^{-\frac{t}{\beta}}\int_0^t \frac{1}{\beta}e^{\frac{s}{\beta}}(\nu(s) + \beta\nu'(s))x(s) \, ds.$$

Using triangle inequalities, for all $t \ge 0$,

$$\begin{aligned} \|y(t)\| &\le e^{-\frac{t}{\beta}}\|y_0\| + \sup_{s \ge 0} \|(\nu(s) + \beta\nu'(s))x(s)\| e^{-\frac{t}{\beta}} \int_0^t \frac{1}{\beta} e^{\frac{s}{\beta}} \, ds \\ &\le \|y_0\| + \sup_{s \ge 0} \|(\nu(s) + \beta\nu'(s))x(s)\|. \quad (27) \end{aligned}$$

Using the definition of $\epsilon$ and $\alpha$ from Sections 1.1 and 2, observe that $\epsilon, \alpha, \epsilon'$ and $\alpha'$ are all bounded on $\mathbb{R}_+$, and $\epsilon''$ is assumed to be bounded. So $\nu$ and $\nu'$ are bounded, and since we also proved that $x$ is uniformly bounded on $\mathbb{R}_+$, we deduce from (27) that $y$ is uniformly bounded as well. Hence, the unique local solution $(x, y)$ is global.

## A.3 Existence and Uniqueness for Non-smooth Functions

Following Remark 2.2, we here provide the main arguments to extend (VM-DIN-AVD) to the non-smooth setting. Assume, in this section only, that $f$ is convex, proper, lower semi-continuous, and (without loss of generality) that $\text{dom}(f) = \mathbb{R}^n$. Consider

$$\begin{cases} \epsilon(t)\dot{x}(t) + \beta\partial f(x(t)) + \nu(t)x(t) + y(t) & \ni 0, \\ \dot{y}(t) + \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t) & = 0, \end{cases} \quad \text{for a.e. } t \ge 0, \quad (28)$$

where $\partial f$ is the sub-differential of $f$. Note that when $f$ is differentiable, (28) boils down to (gVM-DIN-AVD). An absolutely continuous mapping $(x, y) : \mathbb{R}_+ \to \mathbb{R}^n \times \mathbb{R}^n$ that satisfies (28) is called solution. Remark that defining $\tilde{G}(t, x(t), y(t)) = (\frac{\nu(t)x(t)+y(t)}{\epsilon(t)}, \nu'(t)x(t) + \frac{\nu(t)}{\beta}x(t) + \frac{1}{\beta}y(t))$ and $F(t, x(t), y(t)) = (\frac{\beta}{\epsilon(t)}f(x(t)), 0)$, (28) rewrites $(\dot{x}(t), \dot{y}(t)) + \partial F(t, x(t), y(t)) + \tilde{G}(t, x(t), y(t)) \ni 0$. Then under our assumptions on $\epsilon$ and $\alpha$, $\tilde{G}(t, \cdot, \cdot)$ is Lipschitz continuous which allows showing existence of a solution via the theory of non-linear semigroups; see [20, Proposition 3.12]. Note that the chain rule on $(x, y)$ holds for a.e. $t \ge 0$ thanks to absolute continuity, which allows extending our Lyapunov analysis as well.

## B Generalization to Strictly Convex Functions

Following Remark 3.7, we detail how to generalize Section 3 beyond strongly convex functions. Assume, here only, that $f$ is strictly convex. Since $f$ is continuous and the trajectories are contained in
By Assumption 2, for all $s \ge 0$, $|\frac{\varepsilon(s)}{\beta} + \varepsilon'(s)| \le (\frac{1}{\beta} + c_1)\varepsilon(s)$ and $|\frac{\alpha(s)}{\beta} + \alpha'(s)| \le (\frac{1}{\beta} + c_2)\alpha(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,

$$ \beta\mu\|x(t) - x_N(t)\| \le e^{-\frac{t}{\beta}}\varepsilon_0\|\dot{x}_0\| + C\sqrt{\varepsilon(t)} + C\left(\frac{1}{\beta} + c_1\right) \int_0^t e^{\frac{s-t}{\beta}}\sqrt{\varepsilon(s)}\,ds \\ + e^{-\frac{t}{\beta}}\alpha_0\|x_0\| + R\alpha(t) + R\left(\frac{1}{\beta} + c_2\right) \int_0^t e^{\frac{s-t}{\beta}}\alpha(s)\,ds. $$

This proves (8).

**Expressing** $\nabla f(x_{LM})$. We now repeat previous arguments but for (LM). First, (LM) is equivalent to

$$ \frac{d}{dt}\nabla f(x_{LM}(t)) + \frac{1}{\beta}\nabla f(x_{LM}(t)) = -\alpha(t)\dot{x}_{LM}(t). $$

So using an integrating factor one can check that for all $t \ge 0$,

$$ \nabla f(x_{LM}(t)) = e^{-\frac{t}{\beta}}\nabla f(x_0) - e^{-\frac{t}{\beta}}\int_0^t \frac{1}{\beta}e^{\frac{s}{\beta}}\alpha(s)\dot{x}_{LM}(s)\,ds. $$

We can then follow exactly steps (30) to (31) so as to obtain,

$$ e^{-\frac{t}{\beta}} \int_{0}^{t} e^{\frac{s}{\beta}} \alpha(s) \dot{x}_{LM}(s) \, ds = \alpha(t) x_{LM}(t) - e^{-\frac{t}{\beta}} \alpha_0 x_0 - e^{-\frac{t}{\beta}} \int_{0}^{t} \left( \alpha'(s) + \frac{\alpha(s)}{\beta} \right) e^{\frac{s}{\beta}} x_{LM}(s) \, ds. $$

Finally, remark that for all $t \ge 0$,

$$ \frac{d}{dt}f(x_{LM}(t)) = -\alpha(t)\|\dot{x}_{LM}(t)\|^2 - \beta\langle\dot{x}_{LM}(t), \nabla^2 f(x_{LM}(t))\dot{x}_{LM}(t)\rangle \le 0. $$

So $f(x_{LM}(t)) \le f(x_0)$ and using the coercivity of $f$ as before we deduce that for all choices $\alpha$,

$$ \sup_{t \ge 0} \|x_{LM}(t)\| \le R. \qquad (35) $$

**Distance from $x$ to $x_{LM}$**. We subtract gradients,

$$ \beta\nabla f(x(t)) - \beta\nabla f(x_{LM}(t)) = e^{-\frac{t}{\beta}}\varepsilon_0\dot{x}_0 - \varepsilon(t)\dot{x}(t) + \int_0^t e^{\frac{s-t}{\beta}}\left(\frac{1}{\beta}\varepsilon(s) + \varepsilon'(s)\right)\dot{x}(s)\,ds \\ - \alpha(t)(x(t) - x_{LM}(t)) - \int_0^t e^{\frac{s-t}{\beta}}\left(\alpha'(s) + \frac{\alpha(s)}{\beta}\right)(x(s) - x_{LM}(s))\,ds, $$

and we proceed as before using (28), Assumption 2 and Lemma 3.5. It holds that,

$$ \beta\mu\|x(t) - x_{LM}(t)\| \le e^{-\frac{t}{\beta}}\varepsilon_0\|\dot{x}_0\| + C\sqrt{\varepsilon(t)} + C\left(\frac{1}{\beta} + c_1\right) \int_0^t e^{\frac{1}{\beta}(s-t)}\sqrt{\varepsilon(s)}\,ds \\ + \alpha(t)\|x(t) - x_{LM}(t)\| + \left(\frac{1}{\beta} + c_2\right) \int_0^t e^{\frac{1}{\beta}(s-t)}\|x(s) - x_{LM}(s)\|\,ds. $$

Finally, using (33) and (35), for all $s \ge 0$, $\|x(s) - x_{LM}(s)\| \le 2R$, which concludes the proof.
some compact set $K_0$ (see the proof of Theorem 3.1), it follows that $f$ is uniformly convex on $K_0$; see [15, Proposition 10.15]. This implies that $\nabla f$ is uniformly monotone [15, Example 22.3] on $K_0$. That is, there exist an increasing function $\phi: \mathbb{R}_+ \to \mathbb{R}_+$ such that $\phi(0) = 0$, and $\forall y_1, y_2 \in K_0$, $\langle \nabla f(y_1) - \nabla f(y_2), y_1 - y_2 \rangle \ge \phi(\|y_1 - y_2\|)$. Note that when $f$ is strongly convex, $\phi(s) = \mu_{K_0}s^2$, $\forall s \ge 0$. We deduce that $\frac{\phi(\|y_1 - y_2\|)}{\|y_1 - y_2\|} \le \|\nabla f(y_1) - \nabla f(y_2)\|$, which allows extending (3) and our analysis to strictly convex $f$ whose associated $\phi$ is such that $\phi(s)/s \xrightarrow{s \to 0} 0$. In that case, $\phi$ is called a growth function, and this covers many applications; see [26, 40].

### C Proof of Theorem 3.8: Key arguments

This section is devoted to proving the general result of Section 3. Fix some constants $c_1, c_2 > 0$ and let $\varepsilon$ and $\alpha$ such that Assumption 3 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 \ge 0$, $x(t)$, $x_N(t)$ and $x_{LM}(t)$ belong to the bounded set $K_0$ defined in that proof. Since $f$ is $\mu$-strongly convex on $K_0$, the proof relies again on bounding differences of gradients like in (3). Following the exact same steps as in the proof of Theorem 3.1, we know the closed form of $\nabla f(x_N)$ given in (4), and an expression for $\nabla f(x)$ given in (5). For all $s \ge 0$, we have the identity,

$$e^{\frac{s}{\beta}} \dot{x}(s) = e^{\frac{s}{\beta}} \dot{x}(s) + \frac{1}{\beta} e^{\frac{s}{\beta}} x(s) - \frac{1}{\beta} e^{\frac{s}{\beta}} x(s) = \frac{d}{ds} (e^{\frac{s}{\beta}} x(s)) - \frac{1}{\beta} e^{\frac{s}{\beta}} x(s), \quad (29)$$

which we use to perform an integration by parts

$$\int_{0}^{t} e^{\frac{s}{\beta}} \alpha(s) \dot{x}(s) \, ds = \left[ \alpha(s) e^{\frac{s}{\beta}} x(s) \right]_{0}^{t} - \int_{0}^{t} \left( \alpha'(s) + \frac{\alpha(s)}{\beta} \right) e^{\frac{s}{\beta}} x(s) \, ds.$$

Therefore,

$$e^{-\frac{t}{\beta}} \int_{0}^{t} e^{\frac{s}{\beta}} \alpha(s) \dot{x}(s) \, ds = \alpha(t)x(t) - e^{-\frac{t}{\beta}} \alpha_0 x_0 - \int_{0}^{t} e^{\frac{s-t}{\beta}} \left( \alpha'(s) + \frac{\alpha(s)}{\beta} \right) x(s) \, ds, \quad (30)$$

and we can substitute in (5),

$$\beta \nabla f(x(t)) = \beta e^{-\frac{t}{\beta}} \nabla f(x_0) + e^{-\frac{t}{\beta}} \varepsilon_0 \dot{x}_0 - \varepsilon(t) \dot{x}(t) + \int_{0}^{t} e^{\frac{s-t}{\beta}} \left( \frac{1}{\beta} \varepsilon(s) + \varepsilon'(s) \right) \dot{x}(s) \, ds \\ - \alpha(t)x(t) + e^{-\frac{t}{\beta}} \alpha_0 x_0 + \int_{0}^{t} e^{\frac{s-t}{\beta}} \left( \alpha'(s) + \frac{\alpha(s)}{\beta} \right) x(s) \, ds. \quad (31)$$

One then proves (8) by inserting (31) in (3) and by using the uniform boundedness of $x$ and Assumption 3. The proof of (9) follows the same arguments based on the fact that $\forall t \ge 0$,

$$\nabla f(x_{LM}(t)) = e^{-\frac{t}{\beta}} \nabla f(x_0) - e^{-\frac{t}{\beta}} \int_{0}^{t} \frac{1}{\beta} e^{\frac{s}{\beta}} \alpha(s) \dot{x}_{LM}(s) \, ds.$$
C Integrability of $\varphi$ 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 $\varphi$ is continuous, we only need to check its integrability when $t$ tends to $+\infty$. Let $t > 0$, we first establish some useful identities, we omit the dependence on $t$ for the sake of readability.

$$ 
\begin{aligned}
p' &= \frac{\alpha'\varepsilon - (\alpha + \beta\lambda)\varepsilon'}{\varepsilon^2}, \\
p'' &= \frac{\alpha''\varepsilon^2 - 2\alpha'\varepsilon'\varepsilon - (\alpha + \beta\lambda)\varepsilon''\varepsilon + 2(\alpha + \beta\lambda)(\varepsilon')^2}{\varepsilon^3}.
\end{aligned}
$$

Then,

$$ 
\begin{aligned}
r &= \frac{p^2}{4} \left( 1 + \frac{2p'}{p^2} - \frac{4\lambda}{\varepsilon p^2} \right) = \frac{(\alpha + \beta\lambda)^2}{4\varepsilon^2} \left( 1 + \frac{2p'\varepsilon^2}{(\alpha + \beta\lambda)^2} - \frac{4\lambda\varepsilon}{(\alpha + \beta\lambda)^2} \right) \\
&= \frac{(\alpha + \beta\lambda)^2}{4\varepsilon^2} \left( 1 + \frac{2\alpha'\varepsilon}{(\alpha + \beta\lambda)^2} - \frac{2\varepsilon'}{(\alpha + \beta\lambda)} - \frac{4\lambda\varepsilon}{(\alpha + \beta\lambda)^2} \right).
\end{aligned}
\quad (36) $$

An important consequence of Assumption 4 is that $|\varepsilon'(t)| = o(\varepsilon(t))$, $|\varepsilon''(t)| = o(\varepsilon'(t))$ (and the same holds for $\alpha$ w.r.t. to its derivatives). Therefore, we deduce from (36) that

$$ r(t) \sim_{+\infty} \frac{(\alpha(t) + \beta\lambda)^2}{4\varepsilon(t)^2}, $$

and we note that $1/r$ decays at the same speed as $\varepsilon^2$, which will be useful later. In order to study $\varphi$, we now differentiate $r$,

$$ 
\begin{aligned}
r' &= \frac{p'p}{2} \left( 1 + \frac{2p'}{p^2} - \frac{4\lambda}{\varepsilon p^2} \right) + \frac{1}{4} \left( 2p'' - \frac{4(p')^2}{p} + \frac{8\lambda p'}{\varepsilon p} + \frac{4\lambda\varepsilon'}{\varepsilon^2} \right) \\
&= \frac{2p'}{p}r + \frac{1}{4} \left( 2p'' - \frac{4(p')^2}{p} + \frac{8\lambda p'}{\varepsilon p} + \frac{4\lambda\varepsilon'}{\varepsilon^2} \right),
\end{aligned}
$$

and

$$ 
\begin{aligned}
r'' &= 2 \frac{p''p - (p')^2}{p^2} r + \frac{2p'}{p} r' \\
&\quad + \frac{1}{4} \left( 2p''' + 4 \frac{(p')^3 - 2p''p'p}{p^2} + 8\lambda \frac{p''p\varepsilon - (p')^2\varepsilon - p'p\varepsilon'}{\varepsilon^2 p^2} + \frac{4\lambda\varepsilon''}{\varepsilon^2} - \frac{8\lambda(\varepsilon')^2}{\varepsilon^3} \right).
\end{aligned}
\quad (37) $$

Then, to justify that $\varphi$ is integrable, we prove that $\frac{r''}{r^{3/2}}$ and $\frac{(r')^2}{r^{5/2}}$ are integrable. Since we know that $1/r$ decays at the same speed as $\varepsilon^2$, we can equivalently show that $\varepsilon^3 r''$ and $\varepsilon^5 (r')^2$ are integrable. To this aim we fully expand all the terms in (37), and compute $(r')^2$, which is extremely long and involved.
# D Integrability of $\varphi$ and Additional Asymptotic Computations

Below we prove Lemma 4.7.

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

$$p' = \frac{\alpha'\varepsilon - (\alpha + \beta\lambda)\varepsilon'}{\varepsilon^2}, \quad \text{and} \quad p'' = \frac{\alpha''\varepsilon^2 - 2\alpha'\varepsilon'\varepsilon - (\alpha + \beta\lambda)\varepsilon''\varepsilon + 2(\alpha + \beta\lambda)(\varepsilon')^2}{\varepsilon^3}.$$

Then,

$$\begin{aligned} r &= \frac{p^2}{4} \left( 1 + \frac{2p'}{p^2} - \frac{4\lambda}{\varepsilon p^2} \right) = \frac{(\alpha + \beta\lambda)^2}{4\varepsilon^2} \left( 1 + \frac{2p'\varepsilon^2}{(\alpha + \beta\lambda)^2} - \frac{4\lambda\varepsilon}{(\alpha + \beta\lambda)^2} \right) \\ &= \frac{(\alpha + \beta\lambda)^2}{4\varepsilon^2} \left( 1 + \frac{2\alpha'\varepsilon}{(\alpha + \beta\lambda)^2} - \frac{2\varepsilon'}{(\alpha + \beta\lambda)} - \frac{4\lambda\varepsilon}{(\alpha + \beta\lambda)^2} \right). \end{aligned} \quad (32)$$

An important consequence of Assumption 5 is that $|\varepsilon'(t)| = o(\varepsilon(t))$, $|\varepsilon''(t)| = o(\varepsilon'(t))$ (and the same holds for $\alpha$ w.r.t. to its derivatives). Therefore, we deduce from (32) that

$$r(t) \sim_{+\infty} \frac{(\alpha(t) + \beta\lambda)^2}{4\varepsilon(t)^2},$$

and we note that $1/r$ decays at the same speed as $\varepsilon^2$, which will be useful later. In order to study $\varphi$, we now differentiate $r$,

$$\begin{aligned} r' &= \frac{p'p}{2} \left( 1 + \frac{2p'}{p^2} - \frac{4\lambda}{\varepsilon p^2} \right) + \frac{1}{4} \left( 2p'' - \frac{4(p')^2}{p} + \frac{8\lambda p'}{\varepsilon p} + \frac{4\lambda\varepsilon'}{\varepsilon^2} \right) \\ &= \frac{2p'}{p}r + \frac{1}{4} \left( 2p'' - \frac{4(p')^2}{p} + \frac{8\lambda p'}{\varepsilon p} + \frac{4\lambda\varepsilon'}{\varepsilon^2} \right), \quad \text{and,} \end{aligned} \quad (33)$$

$$\begin{aligned} r'' = 2 \frac{p''p - (p')^2}{p^2} r + \frac{2p'}{p} r' \\ + \frac{1}{4} \left( 2p''' + 4 \frac{(p')^3 - 2p''p'}{p^2} + 8\lambda \frac{p''p\varepsilon - (p')^2\varepsilon - p'p\varepsilon'}{\varepsilon^2 p^2} + \frac{4\lambda\varepsilon''}{\varepsilon^2} - \frac{8\lambda(\varepsilon')^2}{\varepsilon^3} \right). \end{aligned} \quad (34)$$

Then, to justify that $\varphi$ is integrable, we prove that $\frac{r''}{r^{3/2}}$ and $\frac{(r')^2}{r^{5/2}}$ are integrable. Since we know that $1/r$ decays at the same speed as $\varepsilon^2$, we can equivalently show that $\varepsilon^3 r''$ and $\varepsilon^5 (r')^2$ are integrable. To this
On the one hand, it holds that,

$$
\begin{aligned}
r'^2 \epsilon^5 = & \left[ -\frac{(\alpha + \beta\lambda)^2 \left(-\frac{4\lambda\epsilon}{(\alpha+\beta\lambda)^2} + 1 + \frac{-2(\alpha+\beta\lambda)\epsilon' + 2\alpha'\epsilon}{(\alpha+\beta\lambda)^2}\right) \epsilon'}{2\sqrt{\epsilon}} + \frac{(\alpha + \beta\lambda)^2 \sqrt{\epsilon}}{4} \left( -\frac{4\lambda\epsilon'}{(\alpha + \beta\lambda)^2} + \frac{8\lambda\alpha'\epsilon}{(\alpha + \beta\lambda)^3} \right) \right. \\
& \left. + \frac{2\left(-\frac{2(\alpha+\beta\lambda)\epsilon'}{\epsilon} + 2\alpha'\epsilon\right) \epsilon'}{(\alpha + \beta\lambda)^2} + \frac{\left(\frac{4(\alpha+\beta\lambda)\epsilon'^2}{\epsilon} - 2(\alpha + \beta\lambda)\epsilon'' - 4\alpha'\epsilon' + 2\alpha''\epsilon\right)}{(\alpha + \beta\lambda)^2} - \frac{2(-2(\alpha + \beta\lambda)\epsilon' + 2\alpha'\epsilon) \alpha'}{(\alpha + \beta\lambda)^3} \right) \\
& + \frac{(\alpha + \beta\lambda)}{2} \left( -\frac{4\lambda\epsilon}{(\alpha + \beta\lambda)^2} + 1 + \frac{(-2(\alpha + \beta\lambda)\epsilon' + 2\alpha'\epsilon)}{(\alpha + \beta\lambda)^2} \right) \alpha' \sqrt{\epsilon} \right]^2
\end{aligned}
$$

On the other hand, we have,

$$
\begin{aligned}
r''\epsilon(t)^3 = & -\lambda\epsilon''\epsilon + \frac{4\lambda\alpha'\epsilon'\epsilon}{\alpha + \beta\lambda} + \frac{2\lambda\alpha''\epsilon^2}{\alpha + \beta\lambda} - \frac{6\lambda\alpha'^2\epsilon^2}{(\alpha + \beta\lambda)^2} + \frac{(\alpha + \beta\lambda)^2}{2} \left(-\frac{3\epsilon'^2}{\epsilon} + \epsilon''\right) \left(\frac{4\lambda\epsilon}{(\alpha + \beta\lambda)^2} - 1 + \frac{2((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)}{(\alpha + \beta\lambda)^2}\right) \\
& + 2(\alpha + \beta\lambda) \left(\frac{4\lambda\epsilon}{(\alpha + \beta\lambda)^2} - 1 + \frac{2((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)}{(\alpha + \beta\lambda)^2}\right) \alpha'\epsilon' - \frac{((\alpha + \beta\lambda)\alpha'' + \alpha'^2)}{2} \left(\frac{4\lambda\epsilon}{(\alpha + \beta\lambda)^2} - 1 + \frac{2((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)}{(\alpha + \beta\lambda)^2}\right) \epsilon \\
& - \left(\frac{(\alpha + \beta\lambda)\epsilon'}{\epsilon} - \alpha'\right) \epsilon'^2 - ((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon) \epsilon'' - 2 \left(-\frac{2(\alpha + \beta\lambda)\epsilon'^2}{\epsilon} + (\alpha + \beta\lambda)\epsilon'' + 2\alpha'\epsilon' - \alpha''\epsilon\right) \epsilon' \\
& + 2 \left(2\lambda\epsilon' - \frac{4\lambda\alpha'\epsilon}{\alpha + \beta\lambda} + 2 \left(\frac{(\alpha + \beta\lambda)\epsilon'}{\epsilon} - \alpha'\right) \epsilon' + \left(-\frac{2(\alpha + \beta\lambda)\epsilon'^2}{\epsilon} + (\alpha + \beta\lambda)\epsilon'' + 2\alpha'\epsilon' - \alpha''\epsilon\right) - \frac{2((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)\alpha'}{\alpha + \beta\lambda}\right) \epsilon' \\
& - \frac{1}{2} \left(\frac{6(\alpha + \beta\lambda)\epsilon'^3}{\epsilon} - 6(\alpha + \beta\lambda)\epsilon'\epsilon'' + (\alpha + \beta\lambda)\epsilon'''\epsilon - 6\alpha'\epsilon'^2 + 3\alpha'\epsilon''\epsilon + 3\alpha''\epsilon'\epsilon - \alpha'''\epsilon^2\right) \\
& + \frac{1}{\alpha + \beta\lambda} \left(4((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)\alpha'\epsilon' + ((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)\alpha''\epsilon + 2(-2(\alpha + \beta\lambda)\epsilon'^2 + (\alpha + \beta\lambda)\epsilon''\epsilon + 2\alpha'\epsilon'\epsilon - \alpha''\epsilon^2)\alpha'\right) \\
& - \frac{2}{\alpha + \beta\lambda} \left(2\lambda\epsilon' - \frac{4\lambda\alpha'\epsilon}{\alpha + \beta\lambda} + 2 \left(\frac{(\alpha + \beta\lambda)\epsilon'}{\epsilon} - \alpha'\right) \epsilon' + \left(-\frac{2(\alpha + \beta\lambda)\epsilon'^2}{\epsilon} + (\alpha + \beta\lambda)\epsilon'' + 2\alpha'\epsilon' - \alpha''\epsilon\right) - \frac{2((\alpha + \beta\lambda)\epsilon' - \alpha'\epsilon)\alpha'}{\alpha + \beta\lambda}\right) \alpha'\epsilon \\
& - \frac{3((\alpha + \beta\lambda)\epsilon'\epsilon - \alpha'\epsilon^2)\alpha'^2}{(\alpha + \beta\lambda)^2}.
\end{aligned}
$$
aim we fully expand all the terms in (33) and (34). We first deal with (33), it holds that:

$$r'^2 \varepsilon^5 = \left[ -\frac{(\alpha + \beta \lambda)^2 \left( -\frac{4\lambda\varepsilon}{(\alpha + \beta \lambda)^2} + 1 + \frac{-2(\alpha + \beta \lambda)\varepsilon' + 2\alpha'\varepsilon}{(\alpha + \beta \lambda)^2} \right) \varepsilon'}{2\sqrt{\varepsilon}} \right. \\ \left. + \frac{(\alpha + \beta \lambda)^2 \sqrt{\varepsilon}}{4} \left( -\frac{4\lambda\varepsilon'}{(\alpha + \beta \lambda)^2} + \frac{8\lambda\alpha'\varepsilon}{(\alpha + \beta \lambda)^3} + \frac{2\left( -\frac{2(\alpha + \beta \lambda)\varepsilon'}{\varepsilon} + 2\alpha'\varepsilon \right) \varepsilon'}{(\alpha + \beta \lambda)^2} \right. \right. \\ \left. + \frac{\left( \frac{4(\alpha + \beta \lambda)\varepsilon'^2}{\varepsilon} - 2(\alpha + \beta \lambda)\varepsilon'' - 4\alpha'\varepsilon' + 2\alpha''\varepsilon \right)}{(\alpha + \beta \lambda)^2} - \frac{2(-2(\alpha + \beta \lambda)\varepsilon' + 2\alpha'\varepsilon)\alpha'}{(\alpha + \beta \lambda)^3} \right) \\ \left. + \frac{(\alpha + \beta \lambda)}{2} \left( -\frac{4\lambda\varepsilon}{(\alpha + \beta \lambda)^2} + 1 + \frac{(-2(\alpha + \beta \lambda)\varepsilon' + 2\alpha'\varepsilon)}{(\alpha + \beta \lambda)^2} \right) \alpha' \sqrt{\varepsilon} \right]^2 .$$

The computations for $\varepsilon^2 r''$ are omitted since they are longer but very similar.

We analyze the integrability of each of the terms above (and those of $\varepsilon^2 r''$). By Assumption 5, $\varepsilon'$, $\varepsilon''$ and $\varepsilon'''$ are integrable, as well as $\alpha'$, $\alpha''$ and $\alpha'''$, this justifies the integrability of most of the terms. We also need $\frac{(\varepsilon')^2}{\varepsilon}$ and $\frac{(\varepsilon')^3}{\varepsilon}$ to be integrable, which holds by Assumption 5. Overall, $\varphi$ is integrable on $\mathbb{R}_+$.

We now state and prove a result used at the end of the proof of Theorem 4.8.

**Lemma D.1.** Under Assumptions 4 and 5, for all $s \ge 0$,

$$\frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 = \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)).$$

*Proof.* We omit the time dependence on $s \ge 0$ for the sake of readability. Using Assumption 4 we can define and expand the following quantity,

$$\frac{1}{16} \left( \frac{2p'}{p^{3/2}} - \frac{4\lambda\sqrt{\varepsilon}}{(\alpha + \beta\lambda)^{3/2}} \right)^2 = \frac{\lambda^2\varepsilon}{(\alpha + \beta\lambda)^3} - \frac{p'\lambda\sqrt{\varepsilon}}{p^{3/2}(\alpha + \beta\lambda)^{3/2}} + \frac{(p')^2}{4p^3} \\ = \frac{\lambda^2\varepsilon}{(\alpha + \beta\lambda)^3} - \frac{\lambda(\alpha'\varepsilon - (\alpha + \beta\lambda)\varepsilon')}{(\alpha + \beta\lambda)^3} + \frac{(\alpha')^2\varepsilon + \frac{(\varepsilon')^2}{\varepsilon}(\alpha + \beta\lambda)^2 - 2\alpha'\varepsilon'(\alpha + \beta\lambda)}{4(\alpha + \beta\lambda)^3}.$$

Assumption 5, implies in particular that $|\varepsilon'(t)| = o(\varepsilon(t))$ and that $\alpha'(t) \to 0$, which we use in the equality above to obtain the desired conclusion.

## E Additional Experimental Details

We used Euler discretization schemes with fixed step-size $\gamma > 0$, for approximating the solutions of the three ODEs considered in Section 5. For a trajectory $x$, at times $t_k = \gamma k, k \in \mathbb{N}$, we denote
We then analyze the integrability of each of the terms above. By Assumption 4, $\varepsilon'$, $\varepsilon''$ and $\varepsilon'''$ are integrable and the same goes for $\alpha'$, $\alpha''$ and $\alpha'''$, which is enough to justify the integrability of almost all the terms above. We finally see that we also need $\frac{(\varepsilon')^2}{\varepsilon}$ and $\frac{(\varepsilon')^3}{\varepsilon}$ to be integrable, which holds by Assumption 4. Overall, $\varphi$ is integrable on $\mathbb{R}_+$.
□

We now state and prove the following result which was used at the end of the proof of Theorem 4.8.

**Lemma C.1.** Under Assumptions 3 and 4, for all $s \ge 0$,

$$ \frac{1}{16} \left( \frac{2p'(s)}{p(s)^{3/2}} - \frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s) + \beta\lambda)^{3/2}} \right)^2 = \frac{\lambda^2\varepsilon(s)}{(\alpha(s) + \beta\lambda)^3} + o(\varepsilon(s)). $$

*Proof.* We omit the time dependence on $s \ge 0$ for the sake of readability. Using Assumption 3 we can define and expand the following quantity,

$$ \frac{1}{16} \left( \frac{2p'}{p^{3/2}} - \frac{4\lambda\sqrt{\varepsilon}}{(\alpha + \beta\lambda)^{3/2}} \right)^2 = \frac{\lambda^2\varepsilon}{(\alpha + \beta\lambda)^3} - \frac{p'\lambda\sqrt{\varepsilon}}{p^{3/2}(\alpha + \beta\lambda)^{3/2}} + \frac{(p')^2}{4p^3} \\ = \frac{\lambda^2\varepsilon}{(\alpha + \beta\lambda)^3} - \lambda(\alpha'\varepsilon - (\alpha + \beta\lambda)\varepsilon') + \frac{1}{4(\alpha + \beta\lambda)^3} \left( (\alpha')^2\varepsilon + \frac{(\varepsilon')^2}{\varepsilon}(\alpha + \beta\lambda)^2 - 2\alpha'\varepsilon'(\alpha + \beta\lambda) \right). $$

Assumption 4, implies in particular that $|\varepsilon'(t)| = o(\varepsilon(t))$ and that $\alpha'(t) \to 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 $\gamma > 0$ and approximate the solutions at times $t_k = \gamma k$, for all $k \in \mathbb{N}$. For a trajectory $x$, we use the notation $x(t_k) \stackrel{\text{def}}{=} x^{(k)}$. The approximation of (CN) is obtained by explicit discretization, so that for all $k \in \mathbb{N}$, we have,

$$ x_N^{(k+1)} = x_N^{(k)} - \gamma \left[ \beta \nabla^2 f(x_N^{(k)}) \right]^{-1} \nabla f(x_N^{(k)}). \quad (38) $$

Then, defining $\varepsilon_k = \varepsilon(t_k)$ and $\alpha_k = \alpha(t_k)$, (LM) and (VM-DIN-AVD) are obtained via Euler semi-implicit discretization. The solution of (LM) is approximated by,

$$ x_{LM}^{(k+1)} = x_{LM}^{(k)} - \gamma \left[ \alpha_k I_n + \beta \nabla^2 f(x_{LM}^{(k)}) \right]^{-1} \nabla f(x_{LM}^{(k)}), \quad (39) $$

where $I_n$ is the identity matrix on $\mathbb{R}^n$. The solution of (VM-DIN-AVD) is obtained similarly,

$$ x^{(k+1)} = x^{(k)} + \left[ (\varepsilon_k + \gamma\alpha_k)I_n + \gamma\beta\nabla^2 f(x^{(k)}) \right]^{-1} (\varepsilon_k(x^{(k)} - x^{(k-1)}) - \gamma^2\nabla f(x^{(k)})). \quad (40) $$

As safety check, one can see that for $\varepsilon_k = 0$, (40) is equivalent to (39), which is itself equivalent to (38) when $\alpha_k = 0$.
$x(t_k) \stackrel{\text{def}}{=} x^{(k)}$. We approximated (CN) by explicit discretization, $\forall k \in \mathbb{N}$,

$$x_N^{(k+1)} = x_N^{(k)} - \gamma \left[ \beta \nabla^2 f(x_N^{(k)}) \right]^{-1} \nabla f(x_N^{(k)}). \quad (35)$$

Then, defining $\varepsilon_k = \varepsilon(t_k)$ and $\alpha_k = \alpha(t_k)$, (LM) and (VM-DIN-AVD) are obtained via Euler semi-implicit discretization. The solution of (LM) is approximated by,

$$x_{LM}^{(k+1)} = x_{LM}^{(k)} - \gamma \left[ \alpha_k I_n + \beta \nabla^2 f(x_{LM}^{(k)}) \right]^{-1} \nabla f(x_{LM}^{(k)}), \quad (36)$$

where $I_n$ is the identity matrix on $\mathbb{R}^n$. The solution of (VM-DIN-AVD) is similarly,

$$x^{(k+1)} = x^{(k)} + \left[ (\varepsilon_k + \gamma \alpha_k) I_n + \gamma \beta \nabla^2 f(x^{(k)}) \right]^{-1} (\varepsilon_k(x^{(k)} - x^{(k-1)}) - \gamma^2 \nabla f(x^{(k)})). \quad (37)$$

As sanity check, one can see that for $\varepsilon_k = 0$, (37) is equivalent to (36), which is itself equivalent to (35) when $\alpha_k = 0$.

## References

[1] Cristian Daniel Alecsa, Szilárd Csaba László, and Titus Pința. An extension of the second order dynamical system that models Nesterov's convex gradient method. *Applied Mathematics & Optimization*, 84(2):1687–1716, 2021.

[2] F Alvarez and JM Pérez. A dynamical system associated with Newton's method for parametric approximations of convex minimization problems. *Applied Mathematics and Optimization*, 38 (2):193–217, 1998.

[3] Felipe Alvarez, Hedy Attouch, Jérôme Bolte, and Patrick Redont. A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics. *Journal de Mathématiques Pures et Appliquées*, 81(8):747–779, 2002.

[4] Hedy Attouch and Alexandre Cabot. Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity. *Journal of Differential Equations*, 263(9):5412–5458, 2017.

[5] Hedy Attouch and Jalal Fadili. From the ravine method to the Nesterov method and vice versa: A dynamical system perspective. *SIAM Journal on Optimization*, 32(3):2074–2101, 2022.

[6] Hedy Attouch and Szilárd Csaba László. Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators. *SIAM Journal on Optimization*, 30(4):3252–3283, 2020.

[7] Hedy Attouch and Szilárd Csaba László. Continuous Newton-like inertial dynamics for monotone inclusions. *Set-Valued and Variational Analysis*, 29(3):555–581, 2021.

[8] Hedy Attouch and Benar Fux Svaiter. A continuous dynamical Newton-like approach to solving monotone inclusions. *SIAM Journal on Control and Optimization*, 49(2):574–598, 2011.

[9] Hedy Attouch, Patrick Redont, and Benar Fux Svaiter. Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in hilbert spaces. *Journal of Optimization Theory and Applications*, 157(3):624–650, 2013.

[10] Hedy Attouch, Maicon Marques Alves, and Benar Fux Svaiter. A dynamic approach to a proximal-Newton method for monotone inclusions in Hilbert spaces, with complexity $O(1/n^2)$. *Journal of Convex Analysis*, 23(1):139–180, 2016.
<figure><img src="image_10.png" /><figcaption>Figure 6: Similar experiment and figures as those described in Figure 2, but for a poorly conditioned quadratic $f(x) = \frac{1}{2}\|Ax\|^2$ (first-two rows) and the function $f(x) = \log(\sum_{i=1}^n e^{x_i}) + \frac{1}{2}\|Ax\|^2$ (last-two rows).</figcaption></figure>
[11] Hedy Attouch, Juan Peypouquet, and Patrick Redont. Fast convex optimization via inertial dynamics with Hessian driven damping. *Journal of Differential Equations*, 261(10):5734–5783, 2016.

[12] Hedy Attouch, Zaki Chbani, Jalal Fadili, and Hassan Riahi. First-order optimization algorithms via inertial systems with Hessian driven damping. *Math. Program.*, 194(4):1–43, 2020.

[13] Hedy Attouch, Aïcha Balhag, Zaki Chbani, and Hassan Riahi. Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling. *Evolution Equations and Control Theory*, 11(2):487–514, 2022.

[14] Hedy Attouch, Radu Ioan Boț, and Ernö Robert Csetnek. Fast optimization via inertial dynamics with closed-loop damping. *Journal of the European Mathematical Society* (published online), 2022.

[15] Heinz Bauschke and Patrick Combettes. *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*. Springer, 2011.

[16] Michel Benaïm, Josef Hofbauer, and Sylvain Sorin. Stochastic approximations and differential inclusions. *SIAM Journal on Control and Optimization*, 44(1):328–348, 2005.

[17] Carl M Bender and Steven A Orszag. *Asymptotic methods and perturbation theory*. Springer, 1999.

[18] Otto Blumenthal. Über asymptotische Integration linearer Differentialgleichungen mit Anwendung auf eine asymptotische Theorie der Kugelfunktionen. *Archiv der Mathematik und Physik*, 19:136–174, 1912.

[19] Radu Ioan Boț, Ernö Robert Csetnek, and Szilárd Csaba László. Tikhonov regularization of a second order dynamical system with Hessian driven damping. *Math. Program.*, 189(1):151–186, 2021.

[20] Haim Brezis. *Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert*. Elsevier, 1973.

[21] Léon Brillouin. Remarques sur la mécanique ondulatoire. *Journal de Physique et Le Radium*, 7 (12):353–368, 1926.

[22] Charles George Broyden. The convergence of a class of double-rank minimization algorithms. *IMA Journal of Applied Mathematics*, 6(1):76–90, 1970.

[23] Camille Castera. Inertial Newton algorithms avoiding strict saddle points. *arXiv:2111.04596*, 2021.

[24] Camille Castera, Jérôme Bolte, Cédric Févotte, and Edouard Pauwels. An inertial Newton algorithm for deep learning. *Journal of Machine Learning Research*, 22(134):1–31, 2021.

[25] Long Chen and Hao Luo. First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow. *arXiv:1912.09276*, 2019.

[26] Dmitriy Drusvyatskiy, Alexander D Ioffe, and Adrian S Lewis. Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria. *Mathematical Programming*, 185:357–383, 2021.

[27] Etienne Emmrich. *Discrete versions of Gronwall's lemma and their application to the numerical analysis of parabolic problems*. TU Fachbereich, 1999.
References

[1] Cristian Daniel Alecsa, Szilárd Csaba László, and Titus Pinţa. An extension of the second order dynamical system that models Nesterov's convex gradient method. *Applied Mathematics & Optimization*, 84(2):1687–1716, 2021.

[2] Felipe Alvarez, Hedy Attouch, Jérôme Bolte, and Patrick Redont. A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics. *Journal de Mathématiques Pures et Appliquées*, 81(8):747–779, 2002.

[3] Hedy Attouch and Alexandre Cabot. Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity. *Journal of Differential Equations*, 263(9):5412–5458, 2017.

[4] Hedy Attouch and Jalal Fadili. From the ravine method to the Nesterov method and vice versa: A dynamical system perspective. *SIAM Journal on Optimization*, 32(3):2074–2101, 2022.

[5] Hedy Attouch and Szilárd Csaba László. Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators. *SIAM Journal on Optimization*, 30(4):3252–3283, 2020.

[6] Hedy Attouch and Szilárd Csaba László. Continuous Newton-like inertial dynamics for monotone inclusions. *Set-Valued and Variational Analysis*, 29(3):555–581, 2021.

[7] Hedy Attouch and Benar Fux Svaiter. A continuous dynamical Newton-like approach to solving monotone inclusions. *SIAM Journal on Control and Optimization*, 49(2):574–598, 2011.

[8] Hedy Attouch, Patrick Redont, and Benar Fux Svaiter. Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in hilbert spaces. *Journal of Optimization Theory and Applications*, 157(3):624–650, 2013.

[9] Hedy Attouch, Juan Peypouquet, and Patrick Redont. Fast convex optimization via inertial dynamics with Hessian driven damping. *Journal of Differential Equations*, 261(10):5734–5783, 2016.

[10] Hedy Attouch, Zaki Chbani, Jalal Fadili, and Hassan Riahi. First-order optimization algorithms via inertial systems with Hessian driven damping. *Math. Program.*, 194(4):1–43, 2020.

[11] Hedy Attouch, Aïcha Balhag, Zaki Chbani, and Hassan Riahi. Fast convex optimization via inertial dynamics combining viscous and Hessian-driven damping with time rescaling. *Evolution Equations and Control Theory*, 11(2):487–514, 2022.

[12] Hedy Attouch, Radu Ioan Boţ, and Ernö Robert Csetnek. Fast optimization via inertial dynamics with closed-loop damping. *Journal of the European Mathematical Society (published online)*, 2022.

[13] Michel Benaïm, Josef Hofbauer, and Sylvain Sorin. Stochastic approximations and differential inclusions. *SIAM Journal on Control and Optimization*, 44(1):328–348, 2005.

[14] Carl M Bender and Steven A Orszag. *Asymptotic methods and perturbation theory*. Springer, 1999.

[15] Otto Blumenthal. Über asymptotische Integration linearer Differentialgleichungen mit Anwendung auf eine asymptotische Theorie der Kugelfunktionen. *Archiv der Mathematik und Physik*, 19:136–174, 1912.
[28] Roger Fletcher. A new approach to variable metric algorithms. *The computer journal*, 13(3): 317-322, 1970.

[29] Mark Konstantinovich Gavurin. Nonlinear functional equations and continuous analogues of iteration methods. *Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika*, 1(5):18-31, 1958.

[30] Donald Goldfarb. A family of variable-metric methods derived by variational means. *Mathematics of computation*, 24(109):23-26, 1970.

[31] George Green. On the motion of waves in a variable canal of small depth and width. *Transactions of the Cambridge Philosophical Society*, 6:457-462, 1838.

[32] John D Hunter. Matplotlib: A 2D graphics environment. *Computing in science & engineering*, 9 (3):90-95, 2007.

[33] Hendrik Anthony Kramers. Wellenmechanik und halbzahlige Quantisierung. *Zeitschrift für Physik*, 39(10):828-840, 1926.

[34] Tianyi Lin and Michael I Jordan. A control-theoretic perspective on optimal high-order optimization. *Math. Program.*, 195(1):929-975, 2022.

[35] Joseph 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. *Journal de Mathématiques Pures et Appliquées*, 2:253-265, 1836.

[36] Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. *Math. Program.*, 45(1):503-528, 1989.

[37] Lennart Ljung. Analysis of recursive stochastic algorithms. *IEEE transactions on automatic control*, 22(4):551-575, 1977.

[38] Yurii Nesterov. A method for unconstrained convex minimization problem with the rate of convergence $O(1/k^2)$. In *Doklady an USSR*, volume 269, pages 543-547, 1983.

[39] Yurii Nesterov. *Introductory lectures on convex optimization: A basic course*, volume 87. Springer Science & Business Media, 2003.

[40] Peter Ochs, Jalal Fadili, and Thomas Brox. Non-smooth non-convex Bregman minimization: Unification and new algorithms. *Journal of Optimization Theory and Applications*, 181(1):244-278, 2019.

[41] Frank Olver. Error bounds for the Liouville-Green (or WKB) approximation. In *Mathematical Proceedings of the Cambridge Philosophical Society*, volume 57, pages 790-810. Cambridge University Press, 1961.

[42] Frank Olver. *Asymptotics and special functions*. Academic Press, 1997.

[43] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. *USSR Computational Mathematics and Mathematical Physics*, 4(5):1-17, 1964.

[44] Guido Rossum. *Python reference manual*. CWI (Centre for Mathematics and Computer Science), 1995.

[45] David F Shanno. Conditioning of quasi-Newton methods for function minimization. *Mathematics of computation*, 24(111):647-656, 1970.
[16] Radu Ioan Boț, Ernö Robert Csetnek, and Szilárd Csaba László. Tikhonov regularization of a second order dynamical system with Hessian driven damping. *Math. Program.*, 189(1):151–186, 2021.

[17] Léon Brillouin. Remarques sur la mécanique ondulatoire. *Journal de Physique et Le Radium*, 7 (12):353–368, 1926.

[18] Charles George Broyden. The convergence of a class of double-rank minimization algorithms. *IMA Journal of Applied Mathematics*, 6(1):76–90, 1970.

[19] Camille Castera. Inertial Newton algorithms avoiding strict saddle points. arXiv:2111.04596, 2021.

[20] Camille Castera, Jérôme Bolte, Cédric Févotte, and Edouard Pauwels. An inertial Newton algorithm for deep learning. *Journal of Machine Learning Research*, 22(134):1–31, 2021.

[21] Long Chen and Hao Luo. First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow. arXiv:1912.09276, 2019.

[22] Etienne Emmrich. *Discrete versions of Gronwall's lemma and their application to the numerical analysis of parabolic problems*. TU Fachbereich, 1999.

[23] Roger Fletcher. A new approach to variable metric algorithms. *The computer journal*, 13(3): 317–322, 1970.

[24] Mark Konstantinovich Gavurin. Nonlinear functional equations and continuous analogues of iteration methods. *Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika*, 1(5):18–31, 1958.

[25] Donald Goldfarb. A family of variable-metric methods derived by variational means. *Mathematics of computation*, 24(109):23–26, 1970.

[26] George Green. On the motion of waves in a variable canal of small depth and width. *Transactions of the Cambridge Philosophical Society*, 6:457–462, 1838.

[27] John D Hunter. Matplotlib: A 2D graphics environment. *Computing in science & engineering*, 9 (3):90–95, 2007.

[28] Hendrik Anthony Kramers. Wellenmechanik und halbzahlige Quantisierung. *Zeitschrift für Physik*, 39(10):828–840, 1926.

[29] Tianyi Lin and Michael I Jordan. A control-theoretic perspective on optimal high-order optimization. *Math. Program.*, 195(1):929–975, 2022.

[30] Joseph 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. *Journal de Mathématiques Pures et Appliquées*, 2:253–265, 1836.

[31] Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. *Math. Program.*, 45(1):503–528, 1989.

[32] Lennart Ljung. Analysis of recursive stochastic algorithms. *IEEE transactions on automatic control*, 22(4):551–575, 1977.

[33] Yurii Nesterov. A method for unconstrained convex minimization problem with the rate of convergence $O(1/k^2)$. In *Doklady an USSR*, volume 269, pages 543–547, 1983.
[46] Bin Shi, Simon S Du, Michael I Jordan, and Weijie J Su. Understanding the acceleration phenomenon via high-resolution differential equations. *Math. Program.*, 195:79–148, 2022.

[47] Weijie Su, Stephen Boyd, and Emmanuel Candes. A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 27, pages 2510–2518, 2014.

[48] James G Taylor. Improved error bounds for the Liouville-Green (or WKB) approximation. *Journal of Mathematical Analysis and Applications*, 85(1):79–89, 1982.

[49] Stéfan van der Walt, Chris Colbert, and Gael Varoquaux. The NumPy array: a structure for efficient numerical computation. *Computing in Science & Engineering*, 13(2):22–30, 2011.

[50] Gregor Wentzel. Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik. *Zeitschrift für Physik*, 38(6):518–529, 1926.
[34] Yurii Nesterov. *Introductory lectures on convex optimization: A basic course*, volume 87. Springer Science & Business Media, 2003.

[35] Frank Olver. *Asymptotics and special functions*. Academic Press, 1997.

[36] Frank William John Olver. Error bounds for the Liouville-Green (or WKB) approximation. In *Mathematical Proceedings of the Cambridge Philosophical Society*, volume 57, pages 790–810. Cambridge University Press, 1961.

[37] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. *USSR Computational Mathematics and Mathematical Physics*, 4(5):1–17, 1964.

[38] Guido Rossum. *Python reference manual*. CWI (Centre for Mathematics and Computer Science), 1995.

[39] David F Shanno. Conditioning of quasi-Newton methods for function minimization. *Mathematics of computation*, 24(111):647–656, 1970.

[40] Bin Shi, Simon S Du, Michael I Jordan, and Weijie J Su. Understanding the acceleration phenomenon via high-resolution differential equations. *Math. Program.*, 195:79–148, 2022.

[41] Weijie Su, Stephen Boyd, and Emmanuel Candes. A differential equation for modeling Nesterov's accelerated gradient method: Theory and insights. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 27, pages 2510–2518, 2014.

[42] James G Taylor. Improved error bounds for the Liouville-Green (or WKB) approximation. *Journal of Mathematical Analysis and Applications*, 85(1):79–89, 1982.

[43] Stéfan van der Walt, Chris Colbert, and Gael Varoquaux. The NumPy array: a structure for efficient numerical computation. *Computing in Science & Engineering*, 13(2):22–30, 2011.

[44] Gregor Wentzel. Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik. *Zeitschrift für Physik*, 38(6):518–529, 1926.