# 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
Department of Mathematics
University of Tübingen
Germany

###### 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.

11footnotetext: Corresponding author<EMAIL_ADDRESS>
##  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).$

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\geq 0,$ (VM-DIN-AVD)

where $f\colon\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\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ (where
$\mathbb{R}_{+}=[0,+\infty[$) are differentiable, non-increasing, and
$\varepsilon(t)>0$ for all $t\geq 0$. Together with $\beta>0$, they are
control parameters that define the type of dynamics that drives the trajectory
(or solution) $x\colon\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{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\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 obtains111CN is
usually considered with $\beta=1$, we put $\beta$ in the system to ease the
discussions below. 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\geq 0,$
(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), 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\geq 0,$ (LM)

also known as regularized Newton method since it stabilizes (CN). In the rest
of the paper, the solutions of (CN) and (LM) will always be denoted by $x_{N}$
and $x_{LM}$ respectively. Alvarez et al. [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 improvement222DIN-like systems were thought to yield faster
vanishing of the gradient compared to inertial first-order dynamics, until
recently [4]. 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.

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\geq\alpha$ (above) and $\varepsilon\leq\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).

###  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 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 $\varepsilon$, and “not too
large” positive dampings $\alpha$, and derive bounds that (formally) take the
form $\|x(t)-x_{N}(t)\|=O(\sqrt{\varepsilon(t)})$. We then extend these
results to larger dampings $\alpha$ 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 $\varepsilon$ and $\alpha$ 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.

Table 1: Informal summary of Section 4. Comparison of (VM-DIN-AVD) with other dynamics Parameters of (VM-DIN-AVD) |  | Speed of convergence
---|---|---
Dominant parameter |  Integrability in $+\infty$ |  |  w.r.t. (CN) |  w.r.t. (LM)
variable mass $\varepsilon$ | yes |  | as fast | as fast
no |  | faster | faster
viscous damping $\alpha$ | yes |  | as fast |  only depends on $\varepsilon$
no |  | slower

###  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
$\varepsilon=1$ and constant $\alpha$ 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
$\alpha$, with in particular the DIN-AVD system with $\alpha(t)=\alpha_{0}/t$
as introduced in [9]. The reason for this popular choice for $\alpha$ is its
link with Nesterov’s method [41]. Non-constant choices for $\beta$ have been
considered [10, 1, 29, 11]. We keep it constant here, and rather focus on non-
constant $\varepsilon$, unlike prior work that used constant $\varepsilon=1$.
The mass $\varepsilon$ was only considered in the original work [2], but only
for fixed $\varepsilon$, $\beta=1$ and constant $\alpha=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
$\varepsilon(t)$. Our approach—which consists in studying the connections with
other second-order dynamics as $\varepsilon$ 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 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^{\prime}(0)=\varepsilon^{\prime}_{0}\leq 0$ and
$\alpha(0)=\alpha_{0}\geq 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^{\prime}(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^{\prime}(t)x(t)+\frac{\nu(t)}{\beta}x(t)+\frac{1}{\beta}y(t)&=0\end{cases}$
(gVM-DIN-AVD)

with initial conditions
$(x(0),y(0))=\left(x_{0},-\varepsilon_{0}\dot{x}_{0}-\beta\nabla
f(x_{0})-(\alpha_{0}-\varepsilon^{\prime}_{0}-\frac{1}{\beta}\varepsilon_{0})x_{0}\right)$.
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\geq 0$ and
$(u,v)\in\mathbb{R}^{n}\times\mathbb{R}^{n}$, define the mapping

$G\left(t,(u,v)\right)=\begin{pmatrix}\frac{1}{\varepsilon(t)}(-\beta\nabla
f(u)-\nu(t)u-v)\\\
-\nu^{\prime}(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\left(t,(x(t),y(t))\right)$ for all $t\geq 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^{\star}\in\mathbb{R}^{n}$. Consequently, any two trajectories that converge
to $x^{\star}$ 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\geq 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}\geq 0$ such that for all $t\geq 0$,
$|\varepsilon^{\prime}(t)|\leq c_{1}\varepsilon(t)$ and $\alpha(t)\leq
c_{2}\varepsilon(t)$.

The assumption states that $\alpha$ must decrease at least as fast as
$\varepsilon$ (up to a constant).333Assumption 1 can actually hold only after
some large-enough $t_{0}\geq 0$, we take $t_{0}=0$ for the sake of simplicity.
The reason for assuming $|\varepsilon^{\prime}(t)|\leq 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.444This is a
consequence of Gronwall’s lemma, see e.g., [22]. 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}\geq 0$. There exist
$C_{0},C_{1},C_{2}\geq 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\geq 0$,

$\|x(t)-x_{N}(t)\|\leq
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)}\mathop{}\\!\mathrm{d}s.$
(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\geq 0$.

######  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\geq 0$, define the function,

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

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

$\frac{\mathop{}\\!\mathrm{d}U}{\mathop{}\\!\mathrm{d}t}(t)=\frac{\varepsilon^{\prime}(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\leq
0.$

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

###### Proof.

Let $t\geq 0$, since $x$ is twice differentiable, $U$ is differentiable and,

$\frac{\mathop{}\\!\mathrm{d}U}{\mathop{}\\!\mathrm{d}t}(t)=\frac{\varepsilon^{\prime}(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{\mathop{}\\!\mathrm{d}U}{\mathop{}\\!\mathrm{d}t}(t)=\frac{\varepsilon^{\prime}(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^{\prime}(t)\leq 0$. Furthermore $f$ is convex so
$\langle\nabla^{2}f(x(t))\dot{x}(t),\dot{x}(t)\rangle\geq 0$. Hence $U$ is
non-increasing. ∎

We then state the following bound.

######  Lemma 3.5.

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

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

###### Proof.

Let $t\geq 0$, according to Lemma 3.4, $U$ is non-increasing so $U(t)\leq
U(0)$, or equivalently,

$\frac{\varepsilon(t)}{2}\|\dot{x}(t)\|^{2}+f(x(t))-f(x^{\star})\leq\frac{\varepsilon_{0}}{2}\|\dot{x}_{0}\|^{2}+f(x_{0})-f(x^{\star}).$

This implies in particular that,

$\varepsilon(t)\|\dot{x}(t)\|^{2}\leq\varepsilon_{0}\|\dot{x}_{0}\|^{2}+2(f(x_{0})-f(x^{\star})),$

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

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

where $C=\sqrt{\varepsilon_{0}\|\dot{x}_{0}\|^{2}+2(f(x_{0})-f(x^{\star}))}$.
∎

###### 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\geq 0$, $U(t)\leq U(0)$, so in particular

$f(x(t))\leq 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 $\mathsf{K}_{0}=\left\\{y\in\mathbb{R}^{n}\mid f(y)\leq M_{0}\right\\}$ is
bounded, since $f$ is coercive ($\lim_{\|y\|\to+\infty}f(y)=+\infty$). So for
all $t\geq 0$, $x(t)\in\mathsf{K}_{0}$. Since $M_{0}$ (and hence
$\mathsf{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 $\mathsf{K}_{0}$ at all
time. Let $x_{N}$ be the solution of (CN). One can see similarly that for all
$t\geq 0$, $f(x_{N}(t))\leq f(x_{N}(0))=f(x_{0})\leq M_{0}$. So we also have
$x_{N}(t)\in\mathsf{K}_{0}$ for all $t\geq 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 $\mathsf{K}_{0}$. We denote $\mu>0$ the strong-convexity
parameter of $f$ on $\mathsf{K}_{0}$. Equivalently, we have that $\nabla f$ is
strongly monotone on $\mathsf{K}_{0}$, that is, $\forall
y_{1},y_{2}\in\mathsf{K}_{0}$,

$\langle\nabla f(y_{1})-\nabla
f(y_{2}),y_{1}-y_{2}\rangle\geq\mu\|y_{1}-y_{2}\|^{2}.$ (2)

Let $t\geq 0$, since $x(t)\in\mathsf{K}_{0}$ and $x_{N}(t)\in\mathsf{K}_{0}$,
by combining (2) with the Cauchy–Schwarz inequality, we deduce that

$\|x(t)-x_{N}(t)\|\leq\frac{1}{\mu}\|\nabla f(x(t))-\nabla f(x_{N}(t))\|.$ (3)

Therefore, it is sufficient to bound the difference of gradients in order to
bound $\|x(t)-x_{N}(t)\|$. First, remark that (CN) can be rewritten as
follows: $\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\nabla
f(x_{N}(t))+\frac{1}{\beta}\nabla f(x_{N}(t))=0$. So we can integrate, for all
$t\geq 0$,

$\nabla f(x_{N}(t))=e^{-\frac{t}{\beta}}\nabla f(x_{0}).$ (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{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\left[\varepsilon(t)\dot{x}(t)+\beta\nabla
f(x(t))\right]+\frac{1}{\beta}\varepsilon(t)\dot{x}(t)+\nabla
f(x(t))=\left(\frac{1}{\beta}\varepsilon(t)+\varepsilon^{\prime}(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^{\prime}(t)-\alpha(t)\right)\dot{x}(t),\quad
t\geq 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^{\prime}(s)-\alpha(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s.$

We thus have the following expression for $\nabla f(x)$, for all $t\geq 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^{\prime}(s)-\alpha(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s.$
(5)

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

$\|x(t)-x_{N}(t)\|\leq\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^{\prime}(s)-\alpha(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s\right\|.$

Using the triangle inequality, we obtain,

$\|x(t)-x_{N}(t)\|\leq\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^{\prime}(s)-\alpha(s)\right|\|\dot{x}(s)\|\mathop{}\\!\mathrm{d}s.$
(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\geq 0$,
$\frac{\varepsilon(t)\|\dot{x}(t)\|}{\beta\mu}\leq 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\geq 0$,

$\left|\frac{1}{\beta}\varepsilon(s)+\varepsilon^{\prime}(s)-\alpha(s)\right|\|\dot{x}(s)\|\leq\left(\frac{1}{\beta}+c_{1}+c_{2}\right)\varepsilon(s)\|\dot{x}(s)\|\leq\left(\frac{1}{\beta}+c_{1}+c_{2}\right)C\sqrt{\varepsilon(s)}.$

This proves the theorem with
$C_{2}=\frac{C}{\beta\mu}\left(\frac{1}{\beta}+c_{1}+c_{2}\right)$. ∎

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\geq
0$,

$\|x(t)-x_{N}(t)\|\leq
C_{0}e^{-\frac{t}{\beta}}\varepsilon_{0}\|\dot{x}_{0}\|+C_{3}\sqrt{\varepsilon(t)}.$

###### Proof of Corollary 3.6.

For all $t\geq 0$, define
$J(t)=\int_{0}^{t}e^{\frac{s}{\beta}}\sqrt{\varepsilon(s)}\mathop{}\\!\mathrm{d}s$.
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^{\prime}(s)}{2\sqrt{\varepsilon(s)}}\mathop{}\\!\mathrm{d}s=\beta
e^{\frac{t}{\beta}}\sqrt{\varepsilon(t)}-\beta\varepsilon_{0}+\int_{s=0}^{t}\beta
e^{\frac{s}{\beta}}\frac{-\varepsilon^{\prime}(s)}{2\varepsilon(s)}\sqrt{\varepsilon(s)}\mathop{}\\!\mathrm{d}s.$
(7)

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

$J(t)\leq\beta
e^{\frac{t}{\beta}}\sqrt{\varepsilon(t)}+c_{1}\frac{\beta}{2}\int_{s=0}^{t}e^{\frac{s}{\beta}}\sqrt{\varepsilon(s)}\mathop{}\\!\mathrm{d}s=\beta
e^{\frac{t}{\beta}}\sqrt{\varepsilon(t)}+c_{1}\frac{\beta}{2}J(t).$

So, $\left(1-c_{1}\frac{\beta}{2}\right)J(t)\leq\beta
e^{t}\sqrt{\varepsilon(t)}$. By assumption $1-c_{1}\frac{\beta}{2}>0$,
therefore,
$J(t)\leq\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. ∎

######  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.

###  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}\geq 0$ such that for all $t\geq 0$,
$|\varepsilon^{\prime}(t)|\leq c_{1}\varepsilon(t)$ and
$|\alpha^{\prime}(t)|\leq 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}\geq 0$. There exist constants $C,\tilde{C}\geq 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\geq 0$,

$\|x(t)-x_{N}(t)\|\leq
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))\mathop{}\\!\mathrm{d}s\right],$
(8)

and,

$\|x(t)-x_{LM}(t)\|\leq\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))\mathop{}\\!\mathrm{d}s\right].$
(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)\leq\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^{\star}$ 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}$,
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^{\star}=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}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
eigenspace555This 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$. 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\geq 0.$ (Q1-VM-DIN-AVD)

Here (and throughout what follows) $\lambda>0$ denotes any eigenvalue of
$A^{T}A$ and $x\colon\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\geq 0,$
(Q1-LM)

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

$\beta\dot{x}_{N}(t)+x_{N}(t)=0,\quad t\geq 0,$ (Q1-CN)

where again, $x_{N}\colon\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\geq 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}\mathop{}\\!\mathrm{d}s\right).$
(10)

Since the minimizer is $x^{\star}=0$, we directly see that $x_{N}$ converges
exponentially fast to $x^{\star}$, 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\geq 0,$ (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\geq 0$, $r(t)\neq 0$, we can use the changes of
variables $\tau=\int_{0}^{t}\sqrt{r(s)}\mathop{}\\!\mathrm{d}s$ and
$w=r^{1/4}z$ and show that $w$ is solution to

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

where666We express $\psi(\tau)$ using $t$ for the sake of readability, using
the one-to-one correspondence between $\tau$ and $t$.
$\psi(\tau)=\frac{4r(t)r^{\prime\prime}(t)-5r^{\prime}(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)}\mathop{}\\!\mathrm{d}s\right)\quad\text{and}\quad\hat{z}_{2}(t)=r(t)^{-1/4}\exp\left(\int_{0}^{t}-\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\right).$
(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\colon\mathbb{R}_{+}\to\mathbb{R}$ be a real, positive, twice
continuously differentiable function, and define
$\varphi(t)=\frac{4r(t)r^{\prime\prime}(t)-5r^{\prime}(t)^{2}}{16r(t)^{5/2}}$
for all $t\geq 0$. Then for any $T>0$, the differential equation,

$\ddot{z}(t)-r(t)z(t)=0,\quad t\in[0,T],$ (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)}\mathop{}\\!\mathrm{d}s\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)}\mathop{}\\!\mathrm{d}s\right),$

where
$\displaystyle|\delta_{1}(t)|\leq\exp\left(\frac{1}{2}\int_{0}^{t}|\varphi(s)|\mathop{}\\!\mathrm{d}s\right)-1$
and
$\displaystyle|\delta_{2}(t)|\leq\exp\left(-\frac{1}{2}\int_{t}^{T}|\varphi(s)|\mathop{}\\!\mathrm{d}s\right)-1$.
If in addition
$\displaystyle\int_{0}^{+\infty}|\varphi(s)|\mathop{}\\!\mathrm{d}s<+\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\geq 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.

###  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\geq 0$,
$\varepsilon_{0}<\frac{(\beta\lambda)^{2}}{2|\alpha^{\prime}(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^{\prime}(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^{\prime}(s)|\mathop{}\\!\mathrm{d}s\leq\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\geq 0$, define

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

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

$\ddot{y}(t)-r(t)y(t)=0,\quad t\geq 0,$ (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\geq 0$, $r(t)$ is positive. Let $t>0$,

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

Since $\varepsilon^{\prime}(t)\leq 0$ and $\alpha^{\prime}(t)\leq 0$, one can
check that a sufficient condition for (17) to hold is,

$r(t)>0\impliedby\frac{(\alpha(t)+\beta\lambda)^{2}}{4}>\left(\frac{|\alpha^{\prime}(t)|}{2}+\lambda\right)\varepsilon(t)\impliedby\frac{(\beta\lambda)^{2}}{2|\alpha^{\prime}(t)|+4\lambda}>\varepsilon_{0}.$

So under Assumption 3, for all $t\geq 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}\mathop{}\\!\mathrm{d}s\right)+\dot{x}(t)\exp\left(\int_{0}^{t}\frac{p(s)}{2}\mathop{}\\!\mathrm{d}s\right),$

and

$\ddot{y}(t)=\exp\left(\int_{0}^{t}\frac{p(s)}{2}\mathop{}\\!\mathrm{d}s\right)\left[\left(\frac{p(t)^{2}}{4}+\frac{p^{\prime}(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}\mathop{}\\!\mathrm{d}s\right)\left(\frac{p(t)^{2}}{4}+\frac{p^{\prime}(t)}{2}-\frac{\lambda}{\varepsilon(t)}\right)x(t)\\\
=\left(\frac{p(t)^{2}}{4}+\frac{p^{\prime}(t)}{2}-\frac{\lambda}{\varepsilon(t)}\right)y(t)=r(t)y(t).$

∎

Lemma 4.5 gives a reformulation of (Q1-VM-DIN-AVD) suited to apply Theorem
4.2. To use the theorem for all $t\geq 0$, we need to ensure that
$\varphi(t)=\frac{4r(t)r^{\prime\prime}(t)-5r^{\prime}(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\limits_{t\to\infty}{\varepsilon(t)}=0$ and
$\varepsilon^{\prime}(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)|\mathop{}\\!\mathrm{d}s<+\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\geq 0$, the
solution of (Q1-VM-DIN-AVD) is

$\displaystyle\begin{split}&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))\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\right),\end{split}$
(18)

where for all $t\geq 0$,

$|\delta_{1}(t)|\leq\exp\left(\frac{1}{2}\int_{0}^{t}|\varphi(s)|\mathop{}\\!\mathrm{d}s\right)-1\quad\text{and}\quad|\delta_{2}(t)|\leq\exp\left(-\frac{1}{2}\int_{t}^{+\infty}|\varphi(s)|\mathop{}\\!\mathrm{d}s\right)-1.$
(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{{\scriptstyle\textrm{def}}}{{=}}x(t)\exp\left(\int_{0}^{t}\frac{p(s)}{2}\mathop{}\\!\mathrm{d}s\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)|\mathop{}\\!\mathrm{d}s<+\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\geq 0$,

$y(t)=A\frac{1+\delta_{1}(t)}{r(t)^{1/4}}\exp\left(\int_{0}^{t}\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\right)+B\frac{1+\delta_{2}(t)}{r(t)^{1/4}}\exp\left(\int_{0}^{t}-\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\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}\mathop{}\\!\mathrm{d}s\right)$,
we obtain that for all $t\geq 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)}\mathop{}\\!\mathrm{d}s\right)+B\frac{1+\delta_{2}(t)}{r(t)^{1/4}}\exp\left(\int_{0}^{t}-\frac{p(s)}{2}-\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\right).$
(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\geq 0$,

$\displaystyle\begin{split}\sqrt{r(s)}&=\frac{p(s)}{2}\sqrt{1+\frac{2p^{\prime}(s)}{p(s)^{2}}-\frac{4\lambda}{\varepsilon(s)p(s)^{2}}}\\\
&=\frac{p(s)}{2}\left(1+\frac{p^{\prime}(s)}{p(s)^{2}}-\frac{2\lambda}{\varepsilon(s)p(s)^{2}}-\frac{1}{8}\left(\frac{2p^{\prime}(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^{\prime}(s)}{2p(s)}-\frac{\lambda}{\varepsilon(s)p(s)}-\frac{1}{16}\left(\frac{2p^{\prime}(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^{\prime}(s)\varepsilon(s)}{2(\alpha(s)+\beta\lambda)}-\frac{\lambda}{\alpha(s)+\beta\lambda}-\frac{1}{16}\left(\frac{2p^{\prime}(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^{\prime}(s)}{2(\alpha(s)+\beta\lambda)}-\frac{\varepsilon^{\prime}(s)}{2\varepsilon(s)}-\frac{\lambda}{\alpha(s)+\beta\lambda}-\frac{1}{16}\left(\frac{2p^{\prime}(s)}{p(s)^{3/2}}-\frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s)+\beta\lambda)^{3/2}}\right)^{2}+o(\varepsilon(s))\end{split}$
(21)

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

$\displaystyle\begin{split}&\exp\left(\int_{0}^{t}-\frac{p(s)}{2}+\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\right)\\\
&=\exp\left(\int_{0}^{t}\frac{\alpha^{\prime}(s)/2}{\alpha(s)+\beta\lambda}-\frac{\varepsilon^{\prime}(s)}{2\varepsilon(s)}-\frac{\lambda}{\alpha(s)+\beta\lambda}-\frac{1}{16}\left(\frac{2p^{\prime}(s)}{p(s)^{3/2}}-\frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s)+\beta\lambda)^{3/2}}\right)^{2}+o(\varepsilon(s))\mathop{}\\!\mathrm{d}s\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^{\prime}(s)}{p(s)^{3/2}}-\frac{4\lambda\sqrt{\varepsilon(s)}}{(\alpha(s)+\beta\lambda)^{3/2}}\right)^{2}+o(\varepsilon(s))\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\right),\end{split}$

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\left(\int_{0}^{t}-\frac{p(s)}{2}-\sqrt{r(s)}\mathop{}\\!\mathrm{d}s\right)$,
and up to redefining $A$ and $B$ so as to encompass all the constants, we
obtain (18) and the result is proved. ∎

###  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 equivalence777Two 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$. 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 3 and 4,
there exists $C>0$ such that the following asymptotic equivalences hold:

$\displaystyle\begin{split}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))\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\right).\end{split}$
(22)

As a consequence, the convergence of $x$ to $x^{\star}$ is:

  1. (i)

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

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

  3. (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\geq 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))\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\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}\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\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))\mathop{}\\!\mathrm{d}s\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}\mathop{}\\!\mathrm{d}s=\int_{0}^{t}-\frac{1}{\beta}+\frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)}\mathop{}\\!\mathrm{d}s=-\frac{t}{\beta}+\int_{0}^{t}\frac{\alpha(s)}{\beta(\alpha(s)+\beta\lambda)}\mathop{}\\!\mathrm{d}s$
(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))\mathop{}\\!\mathrm{d}s$
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)\geq 0$ for
all $s\geq 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))\mathop{}\\!\mathrm{d}s$
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\geq 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)}\mathop{}\\!\mathrm{d}s<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). ∎

##  5 Numerical Experiments

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^{\star}$
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.

Figure 3:  Similar experiment and figures as those described in Figure 2, but
for the function
$f(x)=\log\left(\sum_{i=1}^{n}e^{x_{i}}+e^{-x_{i}}\right)+\frac{1}{2}\|Ax\|^{2}$.

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.

Figure 5:  Numerical validation of Theorem 4.9: distance to the optimum
$x^{\star}$ 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.

We present two set 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.

###  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)\|\leq\|x(t)-x^{\star}\|+\|x_{N}(t)-x^{\star}\|$, 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^{\star}$ 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
$\varepsilon$, 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
$\varepsilon>\alpha$, $x$ may be closer to $x_{N}$ and when
$\varepsilon\leq\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
$\varepsilon$. 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
$\varepsilon$ 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 $\varepsilon(t)=\varepsilon_{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 $\varepsilon$ 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 corresponding888That is, the solution of (LM) for the same
$\alpha$ as that considered for (VM-DIN-AVD). $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.

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].

Appendices

##  Appendix 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,

$\displaystyle\begin{split}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\left[\varepsilon(t)\dot{x}(t)\right]+\beta\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\nabla
f(x(t))+\alpha(t)\dot{x}(t)-\varepsilon^{\prime}(t)\dot{x}(t)+\nabla
f(x(t))=0,\quad t\geq 0.\end{split}$ (24)

We then integrate (24) for all $t\geq 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^{\prime}(s))\dot{x}(s)+\nabla
f(x(s))\mathop{}\\!\mathrm{d}s=0.$ (25)

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

$z(t)=\int_{0}^{t}(\alpha(s)-\varepsilon^{\prime}(s))\dot{x}(s)+\nabla
f(x(s))\mathop{}\\!\mathrm{d}s-\varepsilon_{0}\dot{x}_{0}-\beta\nabla
f(x_{0}).$

We differentiate $z$, for all $t>0$,
$\dot{z}(t)=(\alpha(t)-\varepsilon^{\prime}(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^{\prime}(t))\dot{x}(t)-\nabla
f(x(t))=0\end{cases},\quad t\geq 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^{\prime}(t)-\frac{1}{\beta}\varepsilon(t))\dot{x}(t)+z(t)=0\end{cases},\quad
t\geq 0.$ (26)

To ease the readability, we recall the notation
$\nu(t)=\alpha(t)-\varepsilon^{\prime}(t)-\frac{1}{\beta}\varepsilon(t)$ from
Section 2. Then define for all $t\geq 0$, $y(t)=z(t)-\nu(t)x(t)$, and
differentiate, $\dot{y}(t)=\dot{z}(t)-\nu(t)\dot{x}(t)-\nu^{\prime}(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^{\prime}(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^{\prime}_{0}-\frac{1}{\beta}\varepsilon_{0})x_{0}.$

######  Remark A.1.

Notice that the quantity
$\nu(t)=\alpha(t)-\varepsilon^{\prime}(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.

###  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{\varepsilon(t)}{2}\|\dot{x}(t)\|^{2}+f(x(t))-f(x^{\star})$ 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\geq 0$, $U(t)\geq
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^{\prime}(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^{\prime}(s))x(s)\mathop{}\\!\mathrm{d}s.$

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

$\|y(t)\|\leq e^{-\frac{t}{\beta}}\|y_{0}\|+\sup_{s\geq
0}\|(\nu(s)+\beta\nu^{\prime}(s))x(s)\|e^{-\frac{t}{\beta}}\int_{0}^{t}\frac{1}{\beta}e^{\frac{s}{\beta}}\mathop{}\\!\mathrm{d}s\leq\|y_{0}\|+\sup_{s\geq
0}\|(\nu(s)+\beta\nu^{\prime}(s))x(s)\|.$ (27)

Using the definition of $\varepsilon$ and $\alpha$ from Sections 1.1 and 2,
observe that $\varepsilon$, $\alpha$, $\varepsilon^{\prime}$ and
$\alpha^{\prime}$ are all bounded on $\mathbb{R}_{+}$, and
$\varepsilon^{\prime\prime}$ is assumed to be bounded. So $\nu$ and
$\nu^{\prime}$ 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.

##  Appendix 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 $\varepsilon$ 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\geq 0$, $x(t)$, $x_{N}(t)$ and $x_{LM}(t)$ belong to
the bounded set $\mathsf{K}_{0}$ defined in that proof. Since $f$ is
$\mu$-strongly convex on $\mathsf{K}_{0}$, the proof relies again on bounding
difference of gradients, indeed, for all $t\geq 0$,

$\|x(t)-x_{N}(t)\|\leq\frac{1}{\mu}\|\nabla f(x(t))-\nabla f(x_{N}(t))\|\
\text{and}\ \|x(t)-x_{LM}(t)\|\leq\frac{1}{\mu}\|\nabla f(x(t))-\nabla
f(x_{LM}(t))\|.$ (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\geq 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^{\prime}(s)-\alpha(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s.$

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\geq 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^{\prime}(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s-\int_{0}^{t}e^{\frac{s-t}{\beta}}\alpha(s)\dot{x}(s)\mathop{}\\!\mathrm{d}s.$
(29)

Then, for all $s\geq 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{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}s}(e^{\frac{s}{\beta}}x(s))-\frac{1}{\beta}e^{\frac{s}{\beta}}x(s),$
(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)\mathop{}\\!\mathrm{d}s=\left[\alpha(s)e^{\frac{s}{\beta}}x(s)\right]_{0}^{t}-\int_{0}^{t}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)e^{\frac{s}{\beta}}x(s)\mathop{}\\!\mathrm{d}s.$

Therefore,

$e^{-\frac{t}{\beta}}\int_{0}^{t}e^{\frac{s}{\beta}}\alpha(s)\dot{x}(s)\mathop{}\\!\mathrm{d}s=\alpha(t)x(t)-e^{-\frac{t}{\beta}}\alpha_{0}x_{0}-\int_{0}^{t}e^{\frac{s-t}{\beta}}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)x(s)\mathop{}\\!\mathrm{d}s,$
(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^{\prime}(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s\\\
-\alpha(t)x(t)+e^{-\frac{t}{\beta}}\alpha_{0}x_{0}+\int_{0}^{t}e^{\frac{s-t}{\beta}}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)x(s)\mathop{}\\!\mathrm{d}s.$
(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\geq 0}\|x(t)\|\leq R.$ (33)

We are now in position to prove (8).

#### Distance from $x$ to $x_{N}$.

We first gather (4) and (32). For all $t\geq 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^{\prime}(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s\\\
+e^{-\frac{t}{\beta}}\alpha_{0}x_{0}-\alpha(t)x(t)+\int_{0}^{t}e^{\frac{s-t}{\beta}}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)x(s)\mathop{}\\!\mathrm{d}s.$

We then use (28) and triangle inequalities,

$\beta\mu\|x(t)-x_{N}(t)\|\leq
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^{\prime}(s)\right|\|\dot{x}(s)\|\mathop{}\\!\mathrm{d}s\\\
+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^{\prime}(s)\right|\|x(s)\|\mathop{}\\!\mathrm{d}s.$
(34)

By Assumption 2, for all $s\geq 0$,
$|\frac{\varepsilon(s)}{\beta}+\varepsilon^{\prime}(s)|\leq(\frac{1}{\beta}+c_{1})\varepsilon(s)$
and
$|\frac{\alpha(s)}{\beta}+\alpha^{\prime}(s)|\leq(\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)\|\leq
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)}\mathop{}\\!\mathrm{d}s\\\
+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)\mathop{}\\!\mathrm{d}s.$

This proves (8).

#### Expressing $\nabla f(x_{LM})$.

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

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\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\geq 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)\mathop{}\\!\mathrm{d}s.$

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)\mathop{}\\!\mathrm{d}s=\alpha(t)x_{LM}(t)-e^{-\frac{t}{\beta}}\alpha_{0}x_{0}-e^{-\frac{t}{\beta}}\int_{0}^{t}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)e^{\frac{s}{\beta}}x_{LM}(s)\mathop{}\\!\mathrm{d}s.$

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

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}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\leq
0.$

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

$\sup_{t\geq 0}\|x_{LM}(t)\|\leq R.$ (35)

#### Distance from $x$ to $x_{LM}$.

We substract 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^{\prime}(s)\right)\dot{x}(s)\mathop{}\\!\mathrm{d}s\\\
-\alpha(t)(x(t)-x_{LM}(t))-\int_{0}^{t}e^{\frac{s-t}{\beta}}\left(\alpha^{\prime}(s)+\frac{\alpha(s)}{\beta}\right)(x(s)-x_{LM}(s))\mathop{}\\!\mathrm{d}s,$

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

$\beta\mu\|x(t)-x_{LM}(t)\|\leq
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)}\mathop{}\\!\mathrm{d}s\\\
+\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)\|\mathop{}\\!\mathrm{d}s.$

Finally, using (33) and (35), for all $s\geq 0$, $\|x(s)-x_{LM}(s)\|\leq 2R$,
which concludes the proof.

##  Appendix 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.

$\displaystyle p^{\prime}$
$\displaystyle=\frac{\alpha^{\prime}\varepsilon-(\alpha+\beta\lambda)\varepsilon^{\prime}}{\varepsilon^{2}},$
$\displaystyle p^{\prime\prime}$
$\displaystyle=\frac{\alpha^{\prime\prime}\varepsilon^{2}-2\alpha^{\prime}\varepsilon^{\prime}\varepsilon-(\alpha+\beta\lambda)\varepsilon^{\prime\prime}\varepsilon+2(\alpha+\beta\lambda)(\varepsilon^{\prime})^{2}}{\varepsilon^{3}}.$

Then,

$\displaystyle\begin{split}r&=\frac{p^{2}}{4}\left(1+\frac{2p^{\prime}}{p^{2}}-\frac{4\lambda}{\varepsilon
p^{2}}\right)=\frac{(\alpha+\beta\lambda)^{2}}{4\varepsilon^{2}}\left(1+\frac{2p^{\prime}\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^{\prime}\varepsilon}{(\alpha+\beta\lambda)^{2}}-\frac{2\varepsilon^{\prime}}{(\alpha+\beta\lambda)}-\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}\right).\end{split}$
(36)

An important consequence of Assumption 4 is that
$|\varepsilon^{\prime}(t)|=o(\varepsilon(t))$,
$|\varepsilon^{\prime\prime}(t)|=o(\varepsilon^{\prime}(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$,

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

and

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

Then, to justify that $\varphi$ is integrable, we prove that
$\frac{r^{\prime\prime}}{r^{3/2}}$ and $\frac{(r^{\prime})^{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^{\prime\prime}$ and $\varepsilon^{5}(r^{\prime})^{2}$ are
integrable. To this aim we fully expand all the terms in (37), and compute
$(r^{\prime})^{2}$, which is extremely long and involved.  On the one hand, it
holds that, $\displaystyle\small r^{\prime 2}\varepsilon^{5}=$
$\displaystyle\left[-\frac{(\alpha+\beta\lambda)^{2}\left(-\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}+1+\frac{\left(-{2(\alpha+\beta\lambda)\varepsilon^{\prime}}+{2\alpha^{\prime}\varepsilon}\right)}{(\alpha+\beta\lambda)^{2}}\right)\varepsilon^{\prime}}{2\sqrt{\varepsilon}}+\frac{(\alpha+\beta\lambda)^{2}\sqrt{\varepsilon}}{4}\left(-\frac{4\lambda\varepsilon^{\prime}}{(\alpha+\beta\lambda)^{2}}+\frac{8\lambda\alpha^{\prime}\varepsilon}{(\alpha+\beta\lambda)^{3}}\right.\right.$
$\displaystyle\left.\left.+\frac{2\left(-\frac{2(\alpha+\beta\lambda)\varepsilon^{\prime}}{\varepsilon}+{2\alpha^{\prime}\varepsilon}\right)\varepsilon^{\prime}}{(\alpha+\beta\lambda)^{2}}+\frac{\left(\frac{4(\alpha+\beta\lambda)\varepsilon^{\prime
2}}{\varepsilon}-{2(\alpha+\beta\lambda)\varepsilon^{\prime\prime}}-{4\alpha^{\prime}\varepsilon^{\prime}}+{2\alpha^{\prime\prime}\varepsilon}\right)}{(\alpha+\beta\lambda)^{2}}-\frac{2\left(-{2(\alpha+\beta\lambda)\varepsilon^{\prime}}+{2\alpha^{\prime}\varepsilon}\right)\alpha^{\prime}}{(\alpha+\beta\lambda)^{3}}\right)\right.$
$\displaystyle\left.+\frac{(\alpha+\beta\lambda)}{2}{\left(-\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}+1+\frac{\left(-{2(\alpha+\beta\lambda)\varepsilon^{\prime}}+{2\alpha^{\prime}\varepsilon}\right)}{(\alpha+\beta\lambda)^{2}}\right)}\alpha^{\prime}\sqrt{\varepsilon}\phantom{\frac{\frac{\beta}{\beta}}{\beta}}\right]^{2}.$
On the other hand, we have, $\displaystyle
r^{\prime\prime}\varepsilon(t)^{3}=$
$\displaystyle-\lambda\varepsilon^{\prime\prime}\varepsilon+\frac{4\lambda\alpha^{\prime}\varepsilon^{\prime}\varepsilon}{\alpha+\beta\lambda}+\frac{2\lambda\alpha^{\prime\prime}\varepsilon^{2}}{\alpha+\beta\lambda}-\frac{6\lambda\alpha^{\prime
2}\varepsilon^{2}}{(\alpha+\beta\lambda)^{2}}+\frac{(\alpha+\beta\lambda)^{2}}{2}\left(-\frac{3\varepsilon^{\prime
2}}{\varepsilon}+\varepsilon^{\prime\prime}\right)\left(\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}-1+\frac{2\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)}{(\alpha+\beta\lambda)^{2}}\right)$
$\displaystyle+{2(\alpha+\beta\lambda)\left(\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}-1+\frac{2\left((\alpha+\beta\lambda)\varepsilon^{\prime}-\alpha^{\prime}\varepsilon\right)}{(\alpha+\beta\lambda)^{2}}\right)\alpha^{\prime}\varepsilon^{\prime}}-\frac{\left((\alpha+\beta\lambda)\alpha^{\prime\prime}+\alpha^{\prime
2}\right)}{2}\left(\frac{4\lambda\varepsilon}{(\alpha+\beta\lambda)^{2}}-1+\frac{2\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)}{(\alpha+\beta\lambda)^{2}}\right)\varepsilon$
$\displaystyle-{\left(\frac{(\alpha+\beta\lambda)\varepsilon^{\prime}}{\varepsilon}-\alpha^{\prime}\right)\varepsilon^{\prime
2}}-\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)\varepsilon^{\prime\prime}-2\left(-\frac{2(\alpha+\beta\lambda)\varepsilon^{\prime
2}}{\varepsilon}+{(\alpha+\beta\lambda)\varepsilon^{\prime\prime}}+{2\alpha^{\prime}\varepsilon^{\prime}}-\alpha^{\prime\prime}\varepsilon\right)\varepsilon^{\prime}$
$\displaystyle+{2\left(2\lambda\varepsilon^{\prime}-\frac{4\lambda\alpha^{\prime}\varepsilon}{\alpha+\beta\lambda}+2\left(\frac{(\alpha+\beta\lambda)\varepsilon^{\prime}}{\varepsilon}-\alpha^{\prime}\right)\varepsilon^{\prime}+\left(-\frac{2(\alpha+\beta\lambda)\varepsilon^{\prime
2}}{\varepsilon}+{(\alpha+\beta\lambda)\varepsilon^{\prime\prime}}+{2\alpha^{\prime}\varepsilon^{\prime}}-\alpha^{\prime\prime}\varepsilon\right)-\frac{2\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)\alpha^{\prime}}{\alpha+\beta\lambda}\right)\varepsilon^{\prime}}$
$\displaystyle-\frac{1}{2}\left(\frac{6(\alpha+\beta\lambda)\varepsilon^{\prime
3}}{\varepsilon}-{6(\alpha+\beta\lambda)\varepsilon^{\prime}\varepsilon^{\prime\prime}}+{(\alpha+\beta\lambda)\varepsilon^{\prime\prime\prime}\varepsilon}-{6\alpha^{\prime}\varepsilon^{\prime
2}}+{3\alpha^{\prime}\varepsilon^{\prime\prime}\varepsilon}+{3\alpha^{\prime\prime}\varepsilon^{\prime}\varepsilon}-\alpha^{\prime\prime\prime}\varepsilon^{2}\right)$
$\displaystyle+\frac{1}{\alpha+\beta\lambda}\left({4\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)\alpha^{\prime}\varepsilon^{\prime}}+{\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)\alpha^{\prime\prime}\varepsilon}+{2\left(-{2(\alpha+\beta\lambda)\varepsilon^{\prime
2}}+{(\alpha+\beta\lambda)\varepsilon^{\prime\prime}\varepsilon}+{2\alpha^{\prime}\varepsilon^{\prime}\varepsilon}-\alpha^{\prime\prime}\varepsilon^{2}\right)\alpha^{\prime}}\right)$
$\displaystyle-\frac{2}{\alpha+\beta\lambda}\left(2\lambda\varepsilon^{\prime}-\frac{4\lambda\alpha^{\prime}\varepsilon}{\alpha+\beta\lambda}+2\left(\frac{(\alpha+\beta\lambda)\varepsilon^{\prime}}{\varepsilon}-\alpha^{\prime}\right)\varepsilon^{\prime}+\left(-\frac{2(\alpha+\beta\lambda)\varepsilon^{\prime
2}}{\varepsilon}+{(\alpha+\beta\lambda)\varepsilon^{\prime\prime}}+{2\alpha^{\prime}\varepsilon^{\prime}}-\alpha^{\prime\prime}\varepsilon\right)-\frac{2\left({(\alpha+\beta\lambda)\varepsilon^{\prime}}-\alpha^{\prime}\varepsilon\right)\alpha^{\prime}}{\alpha+\beta\lambda}\right)\alpha^{\prime}\varepsilon$
$\displaystyle-\frac{3\left({(\alpha+\beta\lambda)\varepsilon^{\prime}\varepsilon}-\alpha^{\prime}\varepsilon^{2}\right)\alpha^{\prime
2}}{(\alpha+\beta\lambda)^{2}}.$ We then analyze the integrability of each of
the terms above. By Assumption 4, $\varepsilon^{\prime}$,
$\varepsilon^{\prime\prime}$ and $\varepsilon^{\prime\prime\prime}$ are
integrable and the same goes for $\alpha^{\prime}$, $\alpha^{\prime\prime}$
and $\alpha^{\prime\prime\prime}$, which is enough to justify the
integrability of almost all the terms above. We finally see that we also need
$\frac{(\varepsilon^{\prime})^{2}}{\varepsilon}$ and
$\frac{(\varepsilon^{\prime})^{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\geq 0$,

$\frac{1}{16}\left(\frac{2p^{\prime}(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\geq 0$ for the sake of readability. Using
Assumption 3 we can define and expand the following quantity,

$\displaystyle\frac{1}{16}\left(\frac{2p^{\prime}}{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^{\prime}\lambda\sqrt{\varepsilon}}{p^{3/2}(\alpha+\beta\lambda)^{3/2}}+\frac{(p^{\prime})^{2}}{4p^{3}}$
$\displaystyle=\frac{\lambda^{2}\varepsilon}{(\alpha+\beta\lambda)^{3}}-\lambda(\alpha^{\prime}\varepsilon-(\alpha+\beta\lambda)\varepsilon^{\prime})+\frac{1}{4(\alpha+\beta\lambda)^{3}}\left((\alpha^{\prime})^{2}\varepsilon+\frac{(\varepsilon^{\prime})^{2}}{\varepsilon}(\alpha+\beta\lambda)^{2}-2\alpha^{\prime}\varepsilon^{\prime}(\alpha+\beta\lambda)\right).$

Assumption 4, implies in particular that
$|\varepsilon^{\prime}(t)|=o(\varepsilon(t))$ and that $\alpha^{\prime}(t)\to
0$, which we use in the equality above to obtain the desired conclusion. ∎

##  Appendix 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{{\scriptstyle\textrm{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)}).$ (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)}),$ (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}\left(\varepsilon_{k}(x^{(k)}-x^{(k-1)})-\gamma^{2}\nabla
f(x^{(k)})\right).$ (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$.

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\left(\sum_{i=1}^{n}e^{x_{i}}\right)+\frac{1}{2}\|Ax\|^{2}$ (last-
two rows).

## References

  * Alecsa et al. [2021] 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.
  * Alvarez et al. [2002] 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.
  * Attouch and Cabot [2017] Hedy Attouch and Alexandre Cabot.  Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity.  _Journal of Differential Equations_ , 263(9):5412–5458, 2017.
  * Attouch and Fadili [2022] 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.
  * Attouch and László [2020] 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.
  * Attouch and László [2021] 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.
  * Attouch and Svaiter [2011] 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.
  * Attouch et al. [2013] 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.
  * Attouch et al. [2016] 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.
  * Attouch et al. [2020] 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.
  * Attouch et al. [2022a] 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, 2022a.
  * Attouch et al. [2022b] 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)_ , 2022b.
  * Benaïm et al. [2005] Michel Benaïm, Josef Hofbauer, and Sylvain Sorin.  Stochastic approximations and differential inclusions.  _SIAM Journal on Control and Optimization_ , 44(1):328–348, 2005.
  * Bender and Orszag [1999] Carl M Bender and Steven A Orszag.  _Asymptotic methods and perturbation theory_.  Springer, 1999.
  * Blumenthal [1912] Otto Blumenthal.  Über asymptotische Integration linearer Differentialgleichungen mit Anwendung auf eine asymptotische Theorie der Kugelfunktionen.  _Archiv der Mathematik und Physik_ , 19:136–174, 1912.
  * Boţ et al. [2021] 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.
  * Brillouin [1926] Léon Brillouin.  Remarques sur la mécanique ondulatoire.  _Journal de Physique et Le Radium_ , 7(12):353–368, 1926.
  * Broyden [1970] Charles George Broyden.  The convergence of a class of double-rank minimization algorithms.  _IMA Journal of Applied Mathematics_ , 6(1):76–90, 1970.
  * Castera [2021] Camille Castera.  Inertial Newton algorithms avoiding strict saddle points.  _arXiv:2111.04596_ , 2021.
  * Castera et al. [2021] 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.
  * Chen and Luo [2019] Long Chen and Hao Luo.  First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow.  _arXiv:1912.09276_ , 2019.
  * Emmrich [1999] Etienne Emmrich.  _Discrete versions of Gronwall’s lemma and their application to the numerical analysis of parabolic problems_.  TU Fachbereich, 1999.
  * Fletcher [1970] Roger Fletcher.  A new approach to variable metric algorithms.  _The computer journal_ , 13(3):317–322, 1970\.
  * Gavurin [1958] Mark Konstantinovich Gavurin.  Nonlinear functional equations and continuous analogues of iteration methods.  _Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika_ , 1(5):18–31, 1958.
  * Goldfarb [1970] Donald Goldfarb.  A family of variable-metric methods derived by variational means.  _Mathematics of computation_ , 24(109):23–26, 1970.
  * Green [1838] 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.
  * Hunter [2007] John D Hunter.  Matplotlib: A 2D graphics environment.  _Computing in science & engineering_, 9(3):90–95, 2007.
  * Kramers [1926] Hendrik Anthony Kramers.  Wellenmechanik und halbzahlige Quantisierung.  _Zeitschrift für Physik_ , 39(10):828–840, 1926.
  * Lin and Jordan [2022] Tianyi Lin and Michael I Jordan.  A control-theoretic perspective on optimal high-order optimization.  _Math. Program._ , 195(1):929–975, 2022.
  * Liouville [1836] 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.
  * Liu and Nocedal [1989] Dong C Liu and Jorge Nocedal.  On the limited memory BFGS method for large scale optimization.  _Math. Program._ , 45(1):503–528, 1989.
  * Ljung [1977] Lennart Ljung.  Analysis of recursive stochastic algorithms.  _IEEE transactions on automatic control_ , 22(4):551–575, 1977.
  * Nesterov [1983] 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.
  * Nesterov [2003] Yurii Nesterov.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media, 2003.
  * Olver [1997] Frank Olver.  _Asymptotics and special functions_.  Academic Press, 1997.
  * Olver [1961] 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.
  * Polyak [1964] Boris T Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, 1964.
  * Rossum [1995] Guido Rossum.  _Python reference manual_.  CWI (Centre for Mathematics and Computer Science), 1995.
  * Shanno [1970] David F Shanno.  Conditioning of quasi-Newton methods for function minimization.  _Mathematics of computation_ , 24(111):647–656, 1970.
  * Shi et al. [2022] 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.
  * Su et al. [2014] 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.
  * Taylor [1982] James G Taylor.  Improved error bounds for the Liouville-Green (or WKB) approximation.  _Journal of Mathematical Analysis and Applications_ , 85(1):79–89, 1982.
  * Walt et al. [2011] 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.
  * Wentzel [1926] Gregor Wentzel.  Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik.  _Zeitschrift für Physik_ , 38(6):518–529, 1926.