Dynamical Systems

Motivation

Adam (Kingma and Ba, 2015) is the default optimizer for training neural networks. At step $n$ it keeps an exponential moving average $m_n$ of stochastic gradients and an exponential moving average $v_n$ of their squared coordinates, corrects both for their initialization at zero (the bias correction or debiasing step), and moves the iterate by $-\gamma\,\hat m_n/(\varepsilon + \sqrt{\hat v_n})$ coordinatewise. Despite its use, its convergence on nonconvex objectives was poorly understood when Barakat and Bianchi wrote their paper: Reddi, Kale and Kumar (2018) exhibited convex online problems on which constant-parameter Adam fails to converge, and earlier analyses either dropped the bias correction or replaced Adam by a modified algorithm.

Barakat and Bianchi (arXiv:1810.02263v4, SIAM J. Math. Data Sci. 2021) study Adam through the ODE method: as the step size $\gamma \to 0$ with $1-\alpha \sim a\gamma$ and $1-\beta \sim b\gamma$, the interpolated iterates follow a non-autonomous ordinary differential equation, the continuous-time version of Adam. This mission formalizes their convergence theorem for that ODE (Theorem 3.2): trajectories started as Adam is started, with zero moment estimates, approach the critical points of the objective.

Setting

Let $F : \mathbb R^d \to \mathbb R$ be continuously differentiable with locally Lipschitz gradient $\nabla F$, and let $S : \mathbb R^d \to \mathbb R^d$ be locally Lipschitz. In the paper $F(x) = \mathbb E f(x,\xi)$ is the objective and $S(x) = \mathbb E\,\nabla f(x,\xi)^{\odot 2}$ the coordinatewise second moment of the stochastic gradient; §7 only uses the two properties above. Assume $F$ is coercive ($F(x) \to +\infty$ as $\|x\| \to \infty$) and $S(x) > 0$ coordinatewise for every $x$. Fix constants $a, b, \varepsilon > 0$ with $b \le 4a$.

The state is $z = (x, m, v) \in \mathcal Z = \mathbb R^d \times \mathbb R^d \times \mathbb R^d$, and $\mathcal Z_+$ is the part where $v \ge 0$. All vector operations are coordinatewise. The Adam field (3.3) is, for $t > 0$,

$$h(t, z) = \Big( -\frac{(1-e^{-at})^{-1}\, m}{\varepsilon + \sqrt{(1-e^{-bt})^{-1} v}},\ a(\nabla F(x) - m),\ b(S(x) - v) \Big),$$

and (ODE) is $\dot z(t) = h(t, z(t))$. The factors $(1-e^{-at})^{-1}$ and $(1-e^{-bt})^{-1}$ are the continuous-time bias corrections; they blow up as $t \downarrow 0$. A global solution with initial condition $(x_0, 0, 0)$ is a continuous $z : [0, +\infty) \to \mathcal Z_+$, continuously differentiable on $(0,+\infty)$, satisfying (ODE) for every $t > 0$, with $z(0) = (x_0, 0, 0)$. The critical set is $\mathcal S = \{x : \nabla F(x) = 0\}$.

The proof passes through the autonomous field $h_\infty(z) = \lim_{t\to\infty} h(t,z) = (-m/(\varepsilon+\sqrt v),\ a(\nabla F(x)-m),\ b(S(x)-v))$, its semiflow $\Phi$ on $\mathcal Z_+$, and the notions of equilibrium point, strict Lyapunov function and asymptotic pseudotrajectory (APT) of a semiflow.

Formalization targets

Goal: Theorem 3.2

Assume in addition that $F(\mathcal S)$ has empty interior in $\mathbb R$. For every global solution $z = (x, m, v)$ with initial condition $(x_0, 0, 0)$, the set $\mathcal S$ is non-empty and

$$\lim_{t\to\infty} d(x(t), \mathcal S) = 0, \qquad \lim_{t\to\infty} m(t) = 0, \qquad \lim_{t\to\infty} \big(S(x(t)) - v(t)\big) = 0 .$$

The theorem does not claim that $x(t)$ converges; convergence to a single critical point needs a Łojasiewicz condition (Theorem 3.4, a separate mission).

Milestones

The milestones are the numbered results of §7.1–7.3 used in the proof, in the order the proof uses them: positivity of $v$ along solutions (Lemma 7.4); decrease of $V_\infty$ along $h_\infty$ (Lemma 7.5, autonomous part); compactness of trajectories (Propositions 7.6, 7.7) and lower bounds on $v$ (Lemma 7.8); well-posedness of the autonomous system and its semiflow (Proposition 7.13); the abstract limit-set theorem for APTs (Proposition 7.14, quoted from Benaïm); the strict Lyapunov function $W_\delta(x,m,v) = V_\infty(x,m,v) - \delta\langle\nabla F(x), m\rangle + \delta\|S(x)-v\|^2$ (Proposition 7.15); and the APT property of the (ODE) solution (Proposition 7.17).

Significance

Theorem 3.2 is the deterministic core of the paper's analysis of Adam. Combined with the paper's tracking result (constant-step Adam converges in probability to the ODE solution as $\gamma \to 0$, Theorem 4.3) it says that for small steps Adam's iterates spend their time near critical points of $F$; it is also the input of the decreasing-step almost-sure convergence theorem (Theorem 5.2) and of the convergence-rate theorem (Theorem 3.4), which both rely on the same Lyapunov function and limit-set argument. The condition $b \le 4a$ under which the Lyapunov function decreases is satisfied by the default parameters $\alpha = 0.9$, $\beta = 0.999$.

The result is proved in the paper; to the best of current knowledge it has not been machine-checked. The mission produces a formal proof of the theorem, and reusable formal statements about semiflows on metric spaces: strict Lyapunov functions, and the limit set of a relatively compact asymptotic pseudotrajectory.

Difficulty

The obvious route is to show that $V(t, z(t)) = F(x) + \tfrac12\|m\|^2_{U(t,v)^{-1}}$ decreases and apply LaSalle's invariance principle. This fails on two counts. First, (ODE) is non-autonomous, and LaSalle's principle is a statement about autonomous systems: the derivative of $V$ along the flow vanishes only on a set that depends on $t$. Second, the field is singular: $h(\cdot, z)$ blows up at $t = 0$ and $h(t, \cdot)$ is not Lipschitz near $v = 0$, because of $\sqrt v$, so neither existence and uniqueness nor continuous dependence on initial data come from standard theorems. The energy $V_\infty(z) = \lim_{t\to\infty} V(t,z)$ of the autonomous limit $\dot z = h_\infty(z)$ is a Lyapunov function but not a strict one: it is constant along some non-equilibrium orbits, so non-increase of an energy alone does not identify the limit points.

Formalization scope

Each block $x, m, v$ lives in EuclideanSpace ℝ (Fin d) and $\mathcal Z$ is their product; Mathlib's norm on the product is the maximum of the three Euclidean norms, which is equivalent to the Euclidean norm of $\mathbb R^{3d}$ and changes no statement. $F$ and $S$ are abstract, with the hypotheses of the paper's §7 (Assumptions 7.1, 7.2, 2.3, 2.4); this generality implies the printed theorem for the $F$, $S$ of (2.2). Solutions are maps $\mathbb R \to \mathcal Z$ constrained only on $[0, +\infty)$ (or $[0, T)$), and the equation is required at $t > 0$ only. The field is extended off $\mathcal Z_+$ by $v \mapsto |v|$, as on p. 12 of the paper. A semiflow is Mathlib's Flow ℝ≥0 M, on the subtype $\mathcal Z_+$ or on a compact subset of it. The semiflow $\Phi$ in the milestones is always assumed to be the flow of $(\mathrm{ODE}_\infty)$, never an arbitrary flow. The APT and limit-set definitions are the published StochApproxDyn.LimitSet ones.

The goal is stated for every global solution of the specific field (3.3), and keeps "$\mathcal S$ is non-empty" in the conclusion. A statement about some solution of an unspecified field, or with "$F(\mathcal S)$ has empty interior" strengthened to "$\mathcal S$ is finite", is a different theorem and does not count.

A complete development needs: Grönwall-type comparison and chain-rule arguments along $C^1$ trajectories; well-posedness of the autonomous system on $\mathcal Z_+$ and continuity of its flow; and the limit-set theory of asymptotic pseudotrajectories (Benaïm 1999), adapted to strict Lyapunov functions in the paper's sense. The last part is independent of Adam and reusable in stochastic approximation. Contributions to any milestone are welcome; Propositions 7.13 and 7.14 are the largest pieces of infrastructure.

Selected references

• A. Barakat, P. Bianchi, Convergence and Dynamical Behavior of the ADAM Algorithm for Nonconvex Stochastic Optimization, SIAM J. Math. Data Sci. 3(1), 2021; arXiv:1810.02263v4. https://arxiv.org/abs/1810.02263v4
• D. P. Kingma, J. Ba, Adam: A Method for Stochastic Optimization, ICLR 2015. https://arxiv.org/abs/1412.6980
• S. J. Reddi, S. Kale, S. Kumar, On the Convergence of Adam and Beyond, ICLR 2018. https://arxiv.org/abs/1904.09237
• M. Benaïm, Dynamics of stochastic approximation algorithms, Séminaire de Probabilités XXXIII, Lecture Notes in Math. 1709, Springer, 1999. https://doi.org/10.1007/BFb0096509

Motivation

Accelerated first-order methods for convex minimization, starting with Nesterov's 1983 scheme, reach the rate $f(x_k)-\min f=\mathcal O(1/k^2)$, which is optimal for the class of convex functions with Lipschitz gradient. Su, Boyd and Candès (NIPS 2014; JMLR 2016) showed that Nesterov's method is a discretization of the second-order ODE $\ddot x+\frac3t\dot x+\nabla f(x)=0$, and since then continuous-time inertial dynamics have been the main tool for designing and analysing accelerated algorithms: a Lyapunov argument for the ODE is found first, and then transcribed to the discrete scheme.

Nesterov-type methods oscillate. Adding to the dynamic a damping term driven by the Hessian, $\beta(t)\nabla^2f(x(t))\dot x(t)$, attenuates these oscillations, and because $\nabla^2f(x(t))\dot x(t)$ is the time derivative of $\nabla f(x(t))$, time discretizations of such dynamics remain first-order methods. Attouch, Chbani, Fadili and Riahi (arXiv:1907.10536, Math. Program. 2020) develop this programme. This mission formalizes its first result: a convergence-rate theorem for the continuous dynamic with general damping and time-scale functions.

Timeline.

• 2014: Su, Boyd and Candès identify $\ddot x+\frac3t\dot x+\nabla f(x)=0$ as the continuous limit of Nesterov's method and prove $\mathcal O(1/t^2)$.
• 2016: Attouch, Peypouquet and Redont (J. Differential Equations 261; reference [12] of arXiv:1907.10536) add the constant Hessian damping $\beta\nabla^2 f(x)\dot x$, prove $\mathcal O(1/t^2)$ for $\alpha>3$ and $\int t^2\|\nabla f(x(t))\|^2dt<\infty$.
• 2019–2020: Attouch, Chbani, Fadili and Riahi allow time-dependent $\beta(t)$ and a time scale $b(t)$, and isolate the growth conditions $(\mathcal G_2)$, $(\mathcal G_3)$ under which the Lyapunov argument works (Theorem 1 here).

Setting

Let $\mathcal H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let $f:\mathcal H\to\mathbb R$ be a convex function of class $\mathcal C^2$ that attains its minimum at some point $x^\star$. Fix a starting time $t_0>0$, a constant $\alpha\ge1$, and two non-negative continuous functions $\beta,b:[t_0,+\infty[\to\mathbb R_+$. These are the paper's standing assumptions (H).

A solution trajectory of

$$(\mathrm{DIN\text{-}AVD})_{\alpha,\beta,b}\qquad \ddot x(t)+\frac{\alpha}{t}\dot x(t)+\beta(t)\nabla^2 f(x(t))\dot x(t)+b(t)\nabla f(x(t))=0$$

is a twice differentiable curve $x:[t_0,+\infty[\to\mathcal H$ satisfying this equation at every $t\ge t_0$. The coefficient $\alpha/t$ is an asymptotically vanishing viscous damping, $\beta(t)$ weights the Hessian-driven damping, and $b(t)$ is a time-scale factor in front of the gradient.

When $\beta$ is differentiable, set

$$w(t):=b(t)-\dot\beta(t)-\frac{\beta(t)}{t},\qquad \delta(t):=t^2w(t).$$

The rate in Theorem 1 is expressed through $w$.

Formalization targets

Goal: Theorem 1

Assume $w$ is differentiable and the growth conditions

$$(\mathcal G_2)\quad b(t)>\dot\beta(t)+\frac{\beta(t)}{t},\qquad (\mathcal G_3)\quad t\dot w(t)\le(\alpha-3)w(t)$$

hold for all $t\ge t_0$. Then $w(t)>0$ and

$$f(x(t))-\min_{\mathcal H}f=\mathcal O\!\Big(\frac1{t^2w(t)}\Big),\qquad \int_{t_0}^{\infty}t^2\beta(t)w(t)\|\nabla f(x(t))\|^2dt<\infty,$$ $$\int_{t_0}^{\infty}t\big((\alpha-3)w(t)-t\dot w(t)\big)\big(f(x(t))-\min f\big)dt<\infty .$$

The goal leaves $\beta$ and $b$ general. Specific choices are the companions below.

Milestones

The paper's own proof steps, in order:

1. the derivative (3) of the energy $E(t)=\delta(t)(f(x(t))-f(x^\star))+\frac12\|v(t)\|^2$, where $v(t)=(\alpha-1)(x(t)-x^\star)+t(\dot x(t)+\beta(t)\nabla f(x(t)))$;
2. the identity $\dot v(t)=t[\dot\beta(t)+\beta(t)/t-b(t)]\nabla f(x(t))$;
3. the differential inequality (4);
4. the monotonicity of $E$ via (5)–(6);
5. the two integral bounds by $E(t_0)$.

Companions

• Theorem 3: $\beta(t)\equiv\beta$, $b(t)=1+\beta/t$, $\alpha\ge3$ give $\mathcal O(1/t^2)$ and, for $\beta>0$, $\int t^2\|\nabla f(x(t))\|^2<\infty$.
• Theorem 2: the Attouch–Peypouquet–Redont case $b\equiv1$, $\alpha>3$, with the same conclusions.

Significance

The result. Theorem 1 contains the $\mathcal O(1/t^2)$ rate of the continuous Nesterov dynamic as the case $\beta\equiv0$, $b\equiv1$, $\alpha\ge3$. It recovers the Attouch–Peypouquet–Redont theorem and the vanishing-coefficient system of Case 2 as special cases. The time scale $b$ gives faster rates $\mathcal O(1/(t^2w(t)))$ when $w$ grows. For $\beta>0$ it also gives the integral estimate on $\|\nabla f(x(t))\|^2$, which is the continuous counterpart of the fast decay of gradients proved for the algorithms (IPAHD) and (IGAHD) later in the paper.

Formalizing it. The result is proved in the paper, and nothing about it has been machine-checked. The mission produces:

• a reusable predicate for solutions of second-order ODEs with Hessian damping on a half-line;
• the Lyapunov identity and inequality for that predicate;
• a checked version of the integration step that turns a pointwise differential inequality into finite integrals.

Difficulty

The energy estimate is an exact computation that has to match term by term. The cross terms $\langle\nabla f(x(t)),\dot x(t)\rangle$ cancel only because $\delta=t^2w$ is chosen to match $\dot v$. That computation needs the chain rule $\frac{d}{dt}\nabla f(x(t))=\nabla^2 f(x(t))\dot x(t)$ in a Hilbert space.

The naive conclusion "$\dot E\le0$, so integrate (4)" hides two analytic steps:

• monotonicity from a one-sided derivative on a closed half-line;
• integrability on an unbounded interval of non-negative functions that need not be continuous, because $\dot w$ is only assumed to exist.

The (i) rate needs the positivity of $w$ from $(\mathcal G_2)$ to divide by $\delta(t)$.

Formalization scope

Lean conventions:

• $\mathcal H$ is a real inner product space that is complete. $\nabla f$ is Mathlib's gradient f, and $\nabla^2 f(x)u$ is fderiv ℝ (gradient f) x u.
• Trajectories are maps $\mathbb R\to\mathcal H$ with explicit velocity and acceleration maps, tied to them by HasDerivWithinAt on $[t_0,+\infty[$. Values before $t_0$ are irrelevant, and derivatives at $t_0$ are one-sided.
• $\dot\beta$ and $\dot w$ are explicit maps with derivative hypotheses. The paper writes them, so it presupposes they exist.
• $\min f$ is $f(x^\star)$ for a given global minimizer. $\mathcal O$ is =O[atTop], with constants depending on all data.
• "$\int<\infty$" is IntegrableOn on $[t_0,+\infty[$, never a bound on a Bochner integral that could be the junk value $0$.

These choices rule out trivial readings: Lean's zero value for the derivative of a non-differentiable function, an ODE with an unlinked velocity, and an integral that is zero because the integrand is not integrable. The goal assumes nothing about $E$, $v$ or (4)–(6), which appear only in milestones.

Companion items:

• Theorems 2 and 3 add $\beta>0$ to the gradient-integral claim. The paper derives it from Theorem 1 (ii), which carries a factor $\beta$, so it says nothing at $\beta=0$.
• Theorem 2 must cope with growth conditions that hold only for large $t$.

Contributions are welcome: the chain rule for gradient along a curve, a monotonicity lemma for one-sided derivatives on Ici, and the integrability lemma for non-negative functions bounded by a non-increasing energy.

Selected references

• H. Attouch, Z. Chbani, J. Fadili, H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, Math. Program. (2020); arXiv:1907.10536v2. https://arxiv.org/abs/1907.10536
• H. Attouch, J. Peypouquet, P. Redont, Fast convex minimization via inertial dynamics with Hessian driven damping, J. Differential Equations 261, No. 10 (2016) 5734–5783 (as cited in arXiv:1907.10536v2, ref. [12]). https://arxiv.org/abs/1907.10536
• W. Su, S. Boyd, E. J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, NIPS 27 (2014) 2510–2518; journal version J. Mach. Learn. Res. 17 (2016). https://jmlr.org/papers/v17/15-084.html
• Y. Nesterov, A method of solving a convex programming problem with convergence rate O(1/k²), Soviet Math. Dokl. 27 (1983) 372–376.

Motivation

Adam (Kingma and Ba, 2015) is the default optimizer for training neural networks. At each step it keeps an exponential moving average $m$ of stochastic gradients and an exponential moving average $v$ of their squared coordinates, corrects both for their initialization at zero (the bias correction), and moves the iterate by $-\gamma\,\hat m/(\varepsilon+\sqrt{\hat v})$ coordinatewise.

Barakat and Bianchi (arXiv:1810.02263v4, SIAM J. Math. Data Sci. 2021) study Adam through the ODE method: as the step size $\gamma\to0$ with $1-\alpha\sim a\gamma$ and $1-\beta\sim b\gamma$, the interpolated iterates follow a non-autonomous ordinary differential equation. Their Theorem 3.2 shows that the trajectory of this ODE approaches the critical set of the objective. That leaves open whether the trajectory settles at a single critical point, and how fast. For gradient-like systems the classical answer goes through the Łojasiewicz gradient inequality (Łojasiewicz, 1963, for real-analytic functions), and the dissipative-systems framework of Haraux and Jendoubi (2015) turns such an inequality into convergence with a rate. This mission formalizes the paper's Theorem 3.4, which brings that conclusion to the continuous-time Adam dynamics.

Timeline.

• 1963: Łojasiewicz proves the gradient inequality for real-analytic functions and derives convergence of bounded gradient-flow trajectories.
• 2009: Attouch and Bolte show convergence of the proximal algorithm under the Łojasiewicz property, with rates according to the exponent.
• 2014: Bolte, Sabach and Teboulle prove the same kind of result for proximal alternating linearized minimization.
• 2015: Haraux and Jendoubi collect the continuous-time theory for dissipative autonomous systems.
• 2018–2021: Barakat and Bianchi prove the result for the Adam ODE, which is non-autonomous and has a singular field at $t=0$.

Setting

Let $F:\mathbb R^d\to\mathbb R$ be continuously differentiable with locally Lipschitz gradient $\nabla F$, and let $S:\mathbb R^d\to\mathbb R^d$ be locally Lipschitz. In the paper $F(x)=\mathbb E f(x,\xi)$ and $S(x)=\mathbb E\,\nabla f(x,\xi)^{\odot2}$. Its §7 uses only the two properties above, and so does this mission. Assume $F$ is coercive ($F(x)\to+\infty$ as $\|x\|\to\infty$) and $S(x)>0$ coordinatewise for every $x$. Fix constants $a,b,\varepsilon>0$ with $b\le 4a$.

The state is $z=(x,m,v)\in\mathcal Z=\mathbb R^d\times\mathbb R^d\times\mathbb R^d$, and all vector operations are coordinatewise. The Adam field (3.3) is, for $t>0$,

$$h(t,z)=\Big(-\frac{(1-e^{-at})^{-1}m}{\varepsilon+\sqrt{(1-e^{-bt})^{-1}v}},\ a(\nabla F(x)-m),\ b(S(x)-v)\Big),$$

and (ODE) is $\dot z(t)=h(t,z(t))$. A global solution with initial condition $(x_0,0,0)$ is a continuous $z:[0,+\infty)\to\mathbb R^d\times\mathbb R^d\times[0,+\infty)^d$ that is continuously differentiable on $(0,+\infty)$, satisfies (ODE) for every $t>0$, and has $z(0)=(x_0,0,0)$. The critical set is $\mathcal S=\{x:\nabla F(x)=0\}$.

A number $\theta$ is a Łojasiewicz exponent of $F$ at $x^*$ if there are $c>0$ and $\sigma>0$ with

$$\|\nabla F(x)\|\ \ge\ c\,|F(x)-F(x^*)|^{1-\theta}\qquad\text{whenever }\|x-x^*\|\le\sigma. \tag{3.7}$$

Assumption 3.3 (the Łojasiewicz property) requires an exponent $\theta\in(0,\tfrac12]$ at every $x^*\in\mathcal S$. It holds for real-analytic and semialgebraic $F$.

The proof uses the Lyapunov function $V(t,z)=F(x)+\tfrac12\sum_i m_i^2/U_i(t,v)$, where $U_i(t,v)=a(1-e^{-at})(\varepsilon+\sqrt{v_i/(1-e^{-bt})})$ (3.4)–(3.5). It also uses its perturbation $\tilde W_\delta(t,z)=V(t,z)-\delta\langle\nabla F(x),m\rangle+\delta\|S(x)-v\|^2$ (7.11) and $w_\delta(t)=\tilde W_\delta(t,z(t))$.

Formalization targets

Goal: Theorem 3.4

Assume Assumption 3.3 and that $F(\mathcal S)$ has empty interior in $\mathbb R$. Let $z=(x,m,v)$ be a global solution with initial condition $(x_0,0,0)$. Then there is $x^*\in\mathcal S$ with $x(t)\to x^*$ as $t\to\infty$. Moreover, if $\theta\in(0,\tfrac12]$ is a Łojasiewicz exponent of $F$ at $x^*$, then

$$\|x(t)-x^*\|\le C\,t^{-\frac{\theta}{1-2\theta}}\quad(t>0)\ \text{ if }0<\theta<\tfrac12,\qquad \|x(t)-x^*\|\le C\,e^{-\delta t}\quad(t\ge0)\ \text{ if }\theta=\tfrac12,$$

for a constant $C>0$ (and, in the second case, some $\delta>0$) depending on the trajectory, $x^*$ and $\theta$ but not on $t$.

Milestones

The milestones follow the proof in §7.4:

• the decrease of $V$ along the non-autonomous field (Lemma 7.5, $V$ part);
• convergence to the critical set (Theorem 3.2, restated here);
• the upper bound (7.12) on $w_\delta$ by the four residual terms $|F(x)|$, $\|m\|^2$, $\|\nabla F(x)\|^2$, $\|S(x)-v\|^2$;
• the absolute continuity of $w_\delta$ and the derivative bound (7.13);
• the differential inequality $\dot w_\delta\le-c_4\,w_\delta^{2(1-\theta)}$ of step iv), after the limit value of $F(x(t))$ is normalised to $0$.

Significance

Theorem 3.2 only places the limit points of the trajectory in $\mathcal S$. When $\mathcal S$ contains a continuum of critical points, the trajectory could in principle drift along it forever. Theorem 3.4 excludes this under the Łojasiewicz property. It shows that continuous-time Adam has finite length and converges to one critical point, and it makes the speed explicit. The speed is polynomial with exponent $\theta/(1-2\theta)$ in general and exponential when $\theta=\tfrac12$, which holds for instance at a nondegenerate critical point. The theorem places continuous-time Adam in the same class as gradient flow, heavy-ball dynamics and proximal methods.

The result is proved in the paper (§7.4, with the final integration delegated to arguments of Haraux and Jendoubi). It has not, to current knowledge, been machine-checked. The mission produces a formal proof, together with formal statements of the Lyapunov estimates (7.12)–(7.13) for the Adam ODE. The Łojasiewicz-rate step is reusable for other perturbed gradient-like systems.

Difficulty

The standard Łojasiewicz argument for gradient flow uses $F$ itself as the Lyapunov function, whose derivative is $-\|\nabla F\|^2$. For the Adam ODE neither $F$ nor $V$ works. $V$ decreases only at the rate $\|m\|^2$, which says nothing about $\|\nabla F(x)\|$ when the momentum is small. The system is also non-autonomous, with a field that is singular at $t=0$ and only locally Lipschitz in $z$. As a result $w_\delta$ is only absolutely continuous, so all derivative inequalities hold only almost everywhere. A further obstacle is that the Łojasiewicz inequality at a single critical point does not by itself control the trajectory: the limit set is a connected compact subset of $\mathcal S$, and the inequality has to be made uniform on a neighbourhood of it. Finally, the rate has to be stated for an arbitrary exponent at the limit point, not only for the uniform exponent the argument produces.

Formalization scope

$\mathbb R^d$ is EuclideanSpace ℝ (Fin d) and $\mathcal Z$ is the product E × E × E, whose Mathlib norm is the maximum of the three Euclidean norms. This does not change any statement here, which concern $x$, $m$, $v$ separately. $F$ and $S$ are abstract, under the hypotheses of §7 (ContDiff ℝ 1 F, LocallyLipschitz (gradient F), LocallyLipschitz S, $S>0$ coordinatewise, coercivity as Tendsto F (cocompact _) atTop). Under Assumption 2.2 the paper's $F$ and $S$ satisfy these, so the statements imply the printed ones. A solution is a map ℝ → E × E × E constrained only on $[0,\infty)$. The field uses $|v|$ in place of $v$ (the paper's extension, p. 12) and is never evaluated at $t=0$. Statements are made for every global solution; there is exactly one (Theorem 3.1, a separate mission).

The polynomial bound is stated for $t>0$, because Lean's real power gives $0^{-r}=0$ while the paper's right-hand side is $+\infty$ at $t=0$. The constant $C$ is quantified before $t$. A version with $C$ chosen after $t$, or one that assumes $x(t)$ converges, would be trivial and is not the theorem. (7.13) states absolute continuity of $w_\delta$ on every $[1,T]$ together with the almost-everywhere derivative bound, since an almost-everywhere bound alone says nothing about $w_\delta$.

A complete development needs the Lyapunov calculus of §7.1, compactness of the trajectory, the limit-set theory behind Theorem 3.2, the uniformisation of the Łojasiewicz inequality on a compact connected set of critical points, and the Haraux–Jendoubi integration lemmas. Contributions to the last two are reusable beyond this mission and are welcome, as are proofs of the individual milestones.

Selected references

• A. Barakat and P. Bianchi, Convergence and Dynamical Behavior of the ADAM Algorithm for Nonconvex Stochastic Optimization, SIAM J. Math. Data Sci. 3(1), 2021; arXiv:1810.02263v4. https://arxiv.org/abs/1810.02263v4
• D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, ICLR 2015. https://arxiv.org/abs/1412.6980
• S. Łojasiewicz, Une propriété topologique des sous-ensembles analytiques réels, Les Équations aux Dérivées Partielles, CNRS, 1963, pp. 87–89.
• A. Haraux and M. A. Jendoubi, The Convergence Problem for Dissipative Autonomous Systems, SpringerBriefs in Mathematics, 2015. https://doi.org/10.1007/978-3-319-23407-6
• H. Attouch and J. Bolte, On the convergence of the proximal algorithm for nonsmooth functions involving analytic features, Math. Program. 116, 2009, pp. 5–16. https://doi.org/10.1007/s10107-007-0133-5
• J. Bolte, S. Sabach and M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Math. Program. 146, 2014, pp. 459–494. https://doi.org/10.1007/s10107-013-0701-9

Motivation

A differential variational inequality (DVI) couples an ordinary differential equation with a finite-dimensional variational inequality whose solution enters the dynamics as a control. Contact mechanics with friction, electrical circuits with diodes, dynamic traffic equilibria, hybrid engineering systems and the continuous-time limits of Nash games all have this form. Pang and Stewart's survey (Math. Program. 113, 2008; author's version hal-01366027) set up a common framework for these models and studied the most direct numerical method: the Euler-type time-stepping scheme, which replaces the time derivative by a difference quotient and solves one static variational inequality per time step.

For the scheme to converge, its iterates must exist and be bounded uniformly as the step $h$ goes to zero. Section 7 of the paper proves convergence once such bounds are available (Theorem 7.1). Section 8 supplies them for a class of DVIs in which the algebraic part is not affine: the constraint set is a cone and the VI map is a strongly monotone composite. This mission formalizes that section.

Setting

Fix $T>0$, $\theta\in[0,1]$ and dimensions $n,m,\ell$. The data are $f:[0,T]\times\mathbb R^n\to\mathbb R^n$, $B:[0,T]\times\mathbb R^n\to\mathbb R^{n\times m}$, $G:[0,T]\times\mathbb R^n\to\mathbb R^m$, a map $F:\mathbb R^m\to\mathbb R^m$ and a set $K\subseteq\mathbb R^m$. For $\Phi:\mathbb R^m\to\mathbb R^m$,

$$\mathrm{SOL}(K,\Phi)=\{u\in K:\ (u'-u)^\top\Phi(u)\ge0\ \ \forall u'\in K\}.$$

The initial-value DVI (6.2) is $\dot x=f(t,x)+B(t,x)u$, $u(t)\in\mathrm{SOL}(K,G(t,x(t))+F)$, $x(0)=x^0$.

With $h=T/N$ and $t_{h,i}=ih$, the time-stepping scheme (7.2) starts at $x^{h,0}=x^0$ and computes, for $i=0,\dots,N-1$,

$$x^{h,i+1}=x^{h,i}+h\big[f(t_{h,i+1},\theta x^{h,i}+(1-\theta)x^{h,i+1})+B(t_{h,i},x^{h,i})u^{h,i+1}\big],\qquad u^{h,i+1}\in\mathrm{SOL}(K,G(t_{h,i+1},x^{h,i+1})+F).$$

The standing assumptions are:

• (A) $f$, $B$, $G$ are Lipschitz continuous on $[0,T]\times\mathbb R^n$; (B) $\sigma_B=\sup\|B(t,x)\|<\infty$;
• (D) there is $\eta_G>0$ with $(u-u')^\top\big(G(t,r+B(t_{\rm ref},x^{\rm ref})u)-G(t,r+B(t_{\rm ref},x^{\rm ref})u')\big)\ge\eta_G\|u-u'\|^2$ for all $r,u,u',x^{\rm ref}$ and $t,t_{\rm ref}\in[0,T]$;
• (C′) $F=E^\top\circ\Psi\circ E$ with $E\in\mathbb R^{\ell\times m}$ and $\Psi$ Lipschitz continuous and strongly monotone on $E\mathbb R^m$;
• (E) if $K_1$, $K_2$ are the orthogonal projections of $K$ onto $\ker E$ and $(\ker E)^\perp$, then $K_1\oplus K_2\subseteq K$.

Choosing orthonormal bases $Z$ of $\ker E$ and $W$ of $(\ker E)^\perp$, the map $\Upsilon=(EW)^\top\circ\Psi\circ(EW)$ is strongly monotone with some modulus $\eta_\Upsilon$ (8.7).

Formalization targets

Goal: Theorem 8.1

Let $K$ be a closed convex cone, (A), (B), (D), (C′), $\Psi(0)=0$ and (E) hold. There is $\bar h>0$ such that for every $x^0$ with $\mathrm{SOL}(K,G(0,x^0)+F)\neq\emptyset$ and every $h\in(0,\bar h]$ the scheme has a unique run, and for all small $h$

$$\|x^{h,i+1}\|\le c_{0,x}+c_{1,x}\|x^0\|,\quad \|u^{h,i+1}\|\le c_{0,u}+c_{1,u}\|x^0\|,\quad \|Eu^{h,i+1}-Eu^{h,i}\|\le h\,c_{2,u}.$$

These are the hypotheses (7.5)–(7.6) of the convergence theorem of §7.

Milestones

1. Proposition 8.2: under (A)–(D) and the step bound (8.3), one step of the scheme is solvable, uniquely if $F$ is monotone on $K$.
2. Lemma 8.1: (E) is equivalent to an exchange property of the coordinates $(\mu,\lambda)=(Z^\top u,W^\top u)$.
3. Lemma 8.2: for small $h$, $q\mapsto E\,\mathrm{SOL}(K,q+\Phi_h+F)$ is single-valued and Lipschitz with a constant independent of $h$.
4. Proposition 8.3: for a cone $K$, the step is uniquely solvable and $\|u^h\|$ is bounded by the data.
5. (8.14): $\|u^{h,i+1}\|\le\omega(1+\|x^{h,i}\|+h\|x^{h,i+1}-x^{h,i}\|)$ with $\omega$ independent of $h$ and $x^0$.

A companion item states Proposition 8.1, the derivative form of (D).

Significance

Theorem 8.1 is one of the two routes in the paper to convergence of the time-stepping scheme. The other (Theorem 7.4) needs a linear-growth property of the VI solutions that holds for affine or coercive problems; Theorem 8.1 replaces it by structural assumptions on $(K,F,G)$ that cover nonlinear complementarity systems on cones, including the case $F\equiv0$ that yields the initial-value differential complementarity problem (Corollary 8.1). Combined with Theorem 7.1(a), it shows that the Euler trajectories have uniformly convergent subsequences whose limits are weak solutions of the DVI, and so gives existence of weak solutions constructively.

The paper's results are proved on paper but, as far as a search of Mathlib and the platform shows, no part of DVI theory is formalized. The formalization also settles a printed error: the bound (8.10) of Proposition 8.3 carries a factor $h$ that its proof does not justify, and the printed statement is false (a one-dimensional counterexample is recorded on the milestone). The milestone states the bound the proof does give; (8.14) and Theorem 8.1 are unaffected.

Difficulty

The obvious argument fails at the multiplier bound. The VI map $u\mapsto G(t,x(u))+F(u)$ is strongly monotone, but only with modulus of order $h$, so the solution of each step is bounded only by $O(1/h)$ times the data. This loss is harmless for existence and uniqueness but would destroy the uniform bounds (7.5). Recovering an $O(1)$ bound requires separating the components of $u$ in $\ker E$, where $F$ is blind and only the order-$h$ monotonicity of $G$ helps, from the components in $(\ker E)^\perp$, where $F$ itself is strongly monotone; condition (E) and the cone structure of $K$ are what allow the two parts to be tested separately. A second difficulty is that (7.6) needs a Lipschitz bound on $E\,\mathrm{SOL}$ with a constant that does not blow up as $h\to0$.

Formalization scope

• $\mathbb R^k$ is EuclideanSpace ℝ (Fin k); matrices are continuous linear maps; $E^\top$ is the adjoint; $\mathrm{SOL}(K,\Phi)$ is the published definition SolodovSvaiterVI.Alg21.viSol Φ K.
• The data are functions on $\mathbb R\times\mathbb R^n$; only values with $t\in[0,T]$ enter. (A) uses the metric $|t-t'|+\|x-x'\|$.
• Step sizes with $(N_h+1)h=T$ are $h=T/N$; "for $h\in(0,\bar h]$" and "for $h$ small" are "for $N\ge\bar N$" and "for $N\ge N_1$". $G(t_0,x^0)$ is $G(0,x^0)$.
• (E) is defined with orthogonal projections, not with coordinates, so Lemma 8.1 is not a tautology. Step bounds (8.3), (8.8), (8.9) are multiplied out, so a zero denominator imposes no restriction.
• In Theorem 8.1, uniqueness covers the computed iterates $i\ge1$; $u^{h,0}$ is any element of $\mathrm{SOL}(K,G(0,x^0)+F)$. The constants of (7.5) are uniform in $x^0$; $c_{2,u}$ may depend on $x^0$, as in the proof.
• The last sentence of Theorem 8.1 ("Consequently, the conclusion of Theorem 7.1 holds …") is omitted. It follows by combining this theorem with Theorem 7.1(a), which is the goal of the companion mission Differential Variational Inequalities 2.
• A goal asserting only existence of the iterates, bounds stated for one fixed $h$ rather than for all small $h$ with constants independent of $h$, or (D) restricted to $r=0$ would trivialize the result; none of these is the statement here.

The development needs: existence for strongly monotone (coercive) VIs on closed convex sets; the implicit-function argument for the $x$-equation of a step; orthogonal decompositions along $\ker E$; and discrete Gronwall estimates. The VI existence result and the decomposition lemmas are reusable beyond this mission. Proofs of any milestone, and of the existence theorem for strongly monotone VIs as a separate theorem, are welcome.

Selected references

• J.-S. Pang and D. E. Stewart, Differential variational inequalities, Mathematical Programming 113(2), 2008, 345–424. https://doi.org/10.1007/s10107-006-0052-x (author's version: https://hal.science/hal-01366027v1)
• F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, 2003. https://doi.org/10.1007/b97543
• M. V. Solodov and B. F. Svaiter, A new projection method for variational inequality problems, SIAM J. Control Optim. 37(3), 1999, 765–776. https://doi.org/10.1137/S0363012997317475

Motivation

Networks of coupled oscillators model generators in an electric power grid, pacemaker cells, Josephson junction arrays and many other systems in which units with their own rhythm interact through their phase differences. Dörfler and Bullo (arXiv:0910.5673v4; SIAM J. Control Optim. 50(3), 2012, doi:10.1137/110851584) show that, after a singular-perturbation reduction, the transient-stability question for a lossy power network with non-uniform generators becomes a synchronization question for a non-uniform Kuramoto model: oscillators with different damping coefficients, different natural frequencies, an arbitrary coupling graph and phase shifts induced by line losses.

This mission concerns the cleanest regime of that model: lossless coupling and natural frequencies proportional to the damping. There the oscillators do more than agree on a frequency; they agree on a phase. For the classic Kuramoto model (Kuramoto, 1975) with identical frequencies, phase synchronization from inside an open half-circle was obtained by consensus methods for nonlinear coupled systems (Lin, Francis and Maggiore, 2007) and by Jadbabaie, Motee and Barahona (2004). Theorem V.10 of Dörfler–Bullo extends this to non-uniform damping, directed coupling with a globally reachable node, and gives an explicit worst-case rate for symmetric coupling. The proof rests on the contraction theory of time-varying consensus protocols (Moreau, 2004).

Setting

There are $n \ge 2$ oscillators. Oscillator $i$ has a phase $\theta_i$ on the circle, a damping coefficient $D_i > 0$ and a natural frequency $\omega_i \in \mathbb R$. Oscillators interact through coupling weights $P_{ij} \ge 0$ ($i \ne j$, with $P_{ii} = 0$) and phase shifts $\varphi_{ij} \in [0, \pi/2[$. The non-uniform Kuramoto model (8) is

$$D_i\,\dot\theta_i = \omega_i - \sum_{j=1}^n P_{ij}\sin(\theta_i - \theta_j + \varphi_{ij}), \qquad i = 1,\dots,n .$$

The weights define a directed graph with an edge from $i$ to $j$ whenever $P_{ij} > 0$; a globally reachable node is a node $k$ to which every node has a directed path. For $\gamma \in [0,\pi]$ the closed arc set $\bar\Delta(\gamma)$ is the set of configurations whose angles all lie in a closed arc of length $\gamma$.

This mission assumes zero phase shifts, $\varphi_{\max} = \max_{i,j}\varphi_{ij} = 0$, and proportional frequencies, $\omega_i / D_i = \bar\omega$ for all $i$. The phases synchronize exponentially to a trajectory $\theta_\infty(t)$ when there are $C$ and $\lambda > 0$ with $|\theta_i(t) - \theta_\infty(t)| \le C e^{-\lambda t}$ for all $i$ and $t \ge 0$.

For the rate, $L(a_{ij}) = \mathrm{diag}(\sum_j a_{ij}) - A$ is the Laplacian of a weight array, $\lambda_2(L(P_{ij}))$ is its second-smallest eigenvalue (the algebraic connectivity) when $P = P^T$, $\cos\angle(D\mathbf 1,\mathbf 1) = \sum_i D_i / (\sqrt n \sqrt{\sum_i D_i^2})$, $D_{\max} = \max_i D_i$, $D_{\min} = \min_i D_i$, and $\mathrm{sinc}(x) = \sin(x)/x$.

Formalization targets

Goal: Theorem V.10 (Phase synchronization)

For every solution with $\theta(0) \in \bar\Delta(\gamma)$, $\gamma \in [0,\pi[$:

1. there is $c \in [\theta_{\min}(0), \theta_{\max}(0)]$ such that the phases synchronize exponentially to $\theta_\infty(t) = c + \bar\omega t$;
2. if $P = P^T$, they synchronize exponentially to the weighted mean angle $\theta_\infty(t) = \sum_i D_i\theta_i(0)/\sum_i D_i + \bar\omega t$, and $\delta(t) = \theta(t) - \theta_\infty(t)\mathbf 1$ obeys

$$\|\delta(t)\|_2 \le \sqrt{D_{\max}/D_{\min}}\;\|\delta(0)\|_2\; e^{-\lambda_{\mathrm{ps}} t}, \qquad \lambda_{\mathrm{ps}} = \frac{\lambda_2(L(P_{ij}))\,\mathrm{sinc}(\gamma)\,\cos(\angle(D\mathbf 1,\mathbf 1))^2}{D_{\max}} .$$

Milestones, in the order the proof uses them

• Proof of V.10, p. 27: at a configuration in $\bar\Delta(\gamma)$ with extreme indices $m,\ell$, the arc-length derivative $\dot\theta_m - \dot\theta_\ell = -\sum_k\big(\tfrac{P_{mk}}{D_m}\sin(\theta_m-\theta_k) + \tfrac{P_{\ell k}}{D_\ell}\sin(\theta_k-\theta_\ell)\big)$ is $\le 0$.
• Proof of V.10, p. 27: $\bar\Delta(\gamma)$ is positively invariant for every $\gamma \in [0,\pi[$.
• (43): in the rotating frame $\theta \mapsto \theta - \bar\omega t$ the model is the consensus protocol $\dot\theta_i = -\sum_j a_{ij}(t)(\theta_i - \theta_j)$ with $a_{ij}(t) = (P_{ij}/D_i)\,\mathrm{sinc}(\theta_i(t)-\theta_j(t)) > 0$ on every edge.
• (44): for $P = P^T$, $\tfrac{d}{dt}D\theta = -L(w_{ij}(t))\theta$ in the rotating frame, so $\sum_i D_i\theta_i(t) = \sum_i D_i\theta_i(0) + \bar\omega t\sum_i D_i$.
• Theorem V.10 1) on its own.

Significance

Theorem V.10 gives phase synchronization of non-uniform oscillators from every initial configuration inside an open half-circle, for directed coupling graphs that need only a globally reachable node, and identifies the limit phase: anywhere in the initial range in general, the damping-weighted mean under symmetric coupling. The rate $\lambda_{\mathrm{ps}}$ separates three effects: graph connectivity ($\lambda_2$), initial phase cohesiveness ($\mathrm{sinc}(\gamma)$) and non-uniformity of the damping ($\cos\angle(D\mathbf 1,\mathbf 1)$, $D_{\max}$). For the classic Kuramoto model the statements reduce to known results ([19] and [33] of the paper).

The result is proved in the paper; no machine-checked proof is known. A formal proof needs a Dini-derivative argument for the maximum of finitely many smooth functions along an ODE, a contraction theorem for linear time-varying consensus with positive, bounded weights on a graph with a globally reachable node, and a Courant–Fischer bound on the weighted disagreement. Each is a reusable piece of infrastructure for networked control and synchronization.

Difficulty

The obvious argument linearizes: $\sin(\theta_i - \theta_j) \approx \theta_i - \theta_j$, so the model looks like linear consensus. That only holds near the diagonal; from an arc of length close to $\pi$ the effective weights $\mathrm{sinc}(\theta_i - \theta_j)$ come close to zero and depend on the state, so neither a fixed Laplacian nor a local linearization covers the claimed region. One needs first that the arc never grows, which is a statement about a non-smooth function (the arc length) along trajectories, and then exponential convergence of a consensus protocol whose weights are time-varying and, for directed graphs, non-symmetric, so there is no quadratic Lyapunov function in general. For the explicit rate the disagreement vector is orthogonal to $D\mathbf 1$, not to $\mathbf 1$, and the angle between the two hyperplanes has to be controlled.

Formalization scope

Angles are represented by real lifts $\theta \in \mathbb R^n$; the vector field is $2\pi$-periodic in each coordinate, so solutions on the torus are projections of solutions on $\mathbb R^n$. $\theta \in \bar\Delta(\gamma)$ is encoded as $\theta_i - \theta_j \le \gamma$ for all $i,j$ on the initial lift, and every conclusion (invariance, $\theta_{\min}(0)$, $\theta_{\max}(0)$, the weighted mean, the convergence bounds) refers to that lift and the same continuous lifted trajectory, which is at least as strong as the statement on the torus. A solution is a curve on $[0,\infty)$ whose derivative within $[0,\infty)$ equals the vector field at every $t \ge 0$; statements quantify over every such curve. The model is stated with the paper's general field $D_i^{-1}(\omega_i - \sum_j P_{ij}\sin(\theta_i - \theta_j + \varphi_{ij}))$ and hypotheses $\varphi_{ij} = 0$, $\omega_i/D_i = \bar\omega$ ($\bar\omega$ arbitrary, not set to zero).

Conventions and corrections, each disclosed in the item's Formalization Note:

• $n \ge 2$, $P_{ii} = 0$ and $P_{ij} \ge 0$ for $i \ne j$ are the paper's conventions (§V allows zero and non-symmetric weights).
• (42) is printed with a leading minus sign; the rate is a decay exponent and is stated as the positive number above.
• The page says only "a rate no worse than $\lambda_{\mathrm{ps}}$"; the stated bound uses the constant $\sqrt{D_{\max}/D_{\min}}$ that the referenced proof of Theorem V.1 2) produces (p. 19).
• "$a_{ij}(t)$ is strictly positive" in (43) is stated on the edges, $P_{ij} > 0$; elsewhere $a_{ij} = 0$.
• $\|\cdot\|_2$ is the Euclidean norm, written out (Mathlib's default norm on functions $\mathrm{Fin}\,n \to \mathbb R$ is the sup norm); $\lambda_2$ is the platform definition AlonMilman.PropertyT.lambda1 applied to the symmetric Laplacian.

Exponential convergence requires a rate $\lambda > 0$; a statement with $\lambda = 0$ or with a solution predicate that no curve satisfies would be trivial. The synchronous trajectory $\theta_i(t) = c + \bar\omega t$ satisfies the solution predicate, so the hypotheses are satisfiable.

Infrastructure welcome beyond this mission: upper Dini derivatives and comparison lemmas for maxima of differentiable functions, exponential stability of linear time-varying consensus (Moreau's theorem), and Laplacian ordering $\mathrm{sinc}(\gamma)L(P) \preceq L(w)$ for entrywise-dominated weights.

Selected references

• F. Dörfler and F. Bullo, Synchronization and Transient Stability in Power Networks and Nonuniform Kuramoto Oscillators, SIAM J. Control Optim. 50(3), 2012; this mission cites the arXiv preprint v4. https://arxiv.org/abs/0910.5673v4, https://doi.org/10.1137/110851584
• Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, Lecture Notes in Physics 39, Springer, 1975. https://doi.org/10.1007/BFb0013365
• L. Moreau, Stability of continuous-time distributed consensus algorithms, IEEE Conference on Decision and Control, 2004. https://arxiv.org/abs/math/0409010
• Z. Lin, B. Francis and M. Maggiore, State agreement for continuous-time coupled nonlinear systems, SIAM J. Control Optim. 46(1), 2007. https://doi.org/10.1137/050626405
• N. Chopra and M. W. Spong, On exponential synchronization of Kuramoto oscillators, IEEE Trans. Automatic Control 54(2), 2009. https://doi.org/10.1109/TAC.2008.2007884

Motivation

A differential variational inequality (DVI) couples an ordinary differential equation with a finite-dimensional variational inequality whose solution acts as an algebraic input to the dynamics. The class contains linear complementarity systems, mechanical systems with unilateral contact and friction, dynamic Nash games and many hybrid engineering models. Pang and Stewart's paper (Math. Program. 113 (2008); author's version hal-01366027v1) set up a unified theory of such systems: existence of solutions, uniqueness, and convergence of numerical time-stepping schemes.

The numerical question is the practical one. A DVI is solved by discretizing time and solving one finite-dimensional variational inequality per step. Whether the discrete trajectories approximate a true solution as the step size shrinks is not automatic: the algebraic variable $u$ can switch discontinuously, so the discrete $u$'s need not converge pointwise. Section 7 of the paper answers this for a semi-implicit Euler scheme applied to the initial-value problem, and the answer is also an existence proof for the DVI.

Setting

Fix $T>0$ and write $\Omega=[0,T]\times\mathbb R^n$. Let $K\subseteq\mathbb R^m$ be a nonempty closed convex set, let $f:\Omega\to\mathbb R^n$, $B:\Omega\to\mathbb R^{n\times m}$, $G:\Omega\to\mathbb R^m$ and $F:\mathbb R^m\to\mathbb R^m$. For $\Phi:\mathbb R^m\to\mathbb R^m$, $\mathrm{SOL}(K,\Phi)$ is the set of $u\in K$ with $(u'-u)^{\mathsf T}\Phi(u)\ge0$ for all $u'\in K$. The initial-value DVI (6.2) is

$$\dot x=f(t,x)+B(t,x)u,\qquad x(0)=x^0,\qquad u\in\mathrm{SOL}(K,G(t,x)+F).$$

The standing assumptions are (A): $f$, $B$ and $G$ are Lipschitz continuous on $\Omega$ (for the metric $|t-t'|+\|x-x'\|$), and (B): $\sigma_B=\sup_\Omega\|B(t,x)\|<\infty$.

A weak solution is a pair $(x,u)$ with $x(0)=x^0$, $u$ integrable on $[0,T]$ with $u(t)\in K$ for almost every $t$, the integral equation $x(t)-x(s)=\int_s^t[f(\tau,x(\tau))+B(\tau,x(\tau))u(\tau)]\,d\tau$ for $0\le s\le t\le T$, and the integral form of the VI,

$$\int_0^T(\tilde u(t)-u(t))^{\mathsf T}[G(t,x(t))+F(u(t))]\,dt\ge0\quad\text{for every continuous }\tilde u:[0,T]\to K.$$

The time-stepping scheme (7.2) uses a step $h=T/N$, grid points $t_{h,i}=ih$, a parameter $\theta\in[0,1]$, and $x^{h,0}=x^0$:

$$x^{h,i+1}=x^{h,i}+h\big[f(t_{h,i+1},\theta x^{h,i}+(1-\theta)x^{h,i+1})+B(t_{h,i},x^{h,i})u^{h,i+1}\big],\quad u^{h,i+1}\in\mathrm{SOL}(K,G(t_{h,i+1},x^{h,i+1})+F),$$

for $i=0,\dots,N-1$. The iterates are turned into functions of time: $\hat x^h$ is the continuous piecewise linear interpolant of the $x^{h,i}$, and $\hat u^h$ is the piecewise constant interpolant equal to $u^{h,i+1}$ on $(t_{h,i},t_{h,i+1}]$.

Formalization targets

Goal: Theorem 7.1 (p. 44)

Suppose the iterates satisfy the uniform bounds (7.5),

$$\|x^{h,i+1}\|\le c_{0,x}+c_{1,x}\|x^0\|,\qquad\|u^{h,i+1}\|\le c_{0,u}+c_{1,u}\|x^0\|,$$

for all small $h$. Then along some $h_\nu\downarrow0$,

$$\hat x^{h_\nu}\to\hat x\ \text{uniformly on }[0,T],\qquad\hat u^{h_\nu}\rightharpoonup\hat u\ \text{weakly in }L^2(0,T).$$

If moreover (a) $F=\Psi\circ E$ with $\Psi$ Lipschitz and $\|Eu^{h,i+1}-Eu^{h,i}\|\le hc_{2,u}$ (7.6), or (b) $F(u)=Du$ with $D$ positive semidefinite, then every such limit pair $(\hat x,\hat u)$ is a weak solution of (6.2).

Milestones

The proof decomposes into the paper's own steps: Lemma 7.1 (the implicit Euler step is uniquely solvable, with the bounds (7.4)); the uniform step bound $\|x^{h,i+1}-x^{h,i}\|\le Lh$, which makes $\{\hat x^h\}$ equicontinuous; the fact that a weak $L^2$ limit of $K$-valued functions is $K$-valued almost everywhere; the integral VI (7.7) under (a); the integral equation (7.8); and the weak lower semicontinuity (7.9) of $u\mapsto\int_0^Tu^{\mathsf T}Du$ for positive semidefinite $D$, which handles (b).

Companion results

Lemma 7.2 (a discrete Gronwall bound giving (7.5) from one-step growth conditions), Lemma 7.3 ((7.5) from linear growth of the VI solutions) and Proposition 7.1 (existence of the iterates for small $h$) are included as further targets: together with Theorem 7.1 they reduce the hypotheses on the iterates to conditions on $(K,F)$.

Significance

Theorem 7.1 is simultaneously a convergence theorem for a practical numerical method and an existence theorem: any scheme with bounded iterates produces, in the limit, a weak solution of the DVI. Combined with Proposition 7.1 and Lemma 7.3 it yields existence for DVIs whose VI part has linearly growing solution sets, including linear complementarity systems with positive semidefinite $D$, without the upper-semicontinuity machinery of differential inclusions. Case (b) covers monotone linear complementarity systems, the class arising in passive electrical networks and frictionless contact.

The results are proved in the paper. None of them is machine-checked. A formalization would produce a reusable Lean account of Euler polygons, equicontinuity estimates for discrete schemes, and weak $L^2$ limits of constrained sequences, all of which recur in the numerical analysis of nonsmooth dynamical systems.

Difficulty

The obvious argument would pass to the limit in each step of the scheme pointwise. This fails because $\hat u^h$ has no pointwise or strong limit in general: only weak $L^2$ compactness is available. The nonlinear term $F(\hat u^h)$ does not commute with weak limits, so the variational inequality cannot be passed to the limit directly. Each of the two cases supplies the missing compactness or convexity: in (a) the bound (7.6) makes $E\hat u^h$ converge uniformly, and in (b) positive semidefiniteness makes the quadratic term weakly lower semicontinuous. Membership $\hat u(t)\in K$ also does not follow from pointwise reasoning and needs convexity of $K$.

Formalization scope

Vectors are EuclideanSpace ℝ (Fin k), inner products $u^{\mathsf T}v$ are ⟪u, v⟫_ℝ, and matrices are continuous linear maps with the operator norm. $\mathrm{SOL}(K,\Phi)$ is the published definition SolodovSvaiterVI.Alg21.viSol Φ K. Positive semidefinite matrices are not assumed symmetric.

Step sizes are $h=T/N$ with an integer $N\ge1$, so that the grid ends exactly at $T$. "For all $h\in(0,\bar h]$" becomes "for all $N\ge\bar N$", and "$h_\nu\downarrow0$" becomes a strictly increasing sequence of $N$'s. The iterates are a family indexed by $N$ and are hypotheses of the goal; their existence is the separate Proposition 7.1. Weak convergence in $L^2(0,T;\mathbb R^m)$ is tested against every $L^2$ function, and the uniform convergence of $\hat x^h$ is on all of $[0,T]$. Condition (7.6) is imposed for $1\le i\le N-1$, because the paper's convention $u^{h,0}\equiv u^{h,1}$ makes $i=0$ vacuous.

A weak solution carries integrability of both integrands as explicit conditions. Without them a non-integrable pair would satisfy the integral conditions vacuously, since a Lean integral of a non-integrable function is zero. The following readings are also excluded: convergence of the scheme only at grid points; a scheme in which $u^{h,i+1}$ solves the VI at $x^{h,i}$ instead of $x^{h,i+1}$; and a conclusion about one particular subsequence instead of every subsequence along which both limits exist.

A complete development needs: compactness in $C([0,T])$ (Arzelà–Ascoli), weak sequential compactness of bounded sets in $L^2$, Mazur's lemma on convex combinations, and estimates for piecewise linear and piecewise constant interpolants. The last two items are reusable for other time-stepping schemes. Proofs of any milestone, and of auxiliary interpolation lemmas, are welcome.

Out of scope: the weaker Carathéodory assumption (A′), the implicit scheme (7.3), the cited fixed-point theorems 7.2–7.3, and Theorem 7.4 (which combines Theorem 7.1 with the existence results of §6).

Selected references

• J.-S. Pang and D. E. Stewart, Differential variational inequalities, Mathematical Programming 113(2), 345–424, 2008. https://doi.org/10.1007/s10107-006-0052-x ; author's version hal-01366027v1, https://hal.science/hal-01366027v1
• F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, 2003. https://doi.org/10.1007/b97543
• D. E. Stewart, Rigid-body dynamics with friction and impact, SIAM Review 42(1), 3–39, 2000. https://doi.org/10.1137/S0036144599360110