arXiv CSopen access

Dirac-Frenkel dynamics with inertia for nonlinearly parametrized solutions of evolution problems

Dirac–Frenkel dynamics with inertia for nonlinearly parametrized solutions of evolution problems

Matteo Raviola

arXiv:2606.24769v1 [math.NA] 23 Jun 2026

Scientific Computing and Uncertainty Quantification - CADMOS Chair, EPFL

Benjamin Peherstorfer Courant Institute of Mathematical Sciences, New York University

Abstract Even when Dirac–Frenkel dynamics determine a well-defined evolution in function space, the corresponding parameter dynamics can be non-unique or ill-conditioned for redundant nonlinear parametrizations such as neural networks or mixture models. We propose to add inertia to the Dirac-Frenkel dynamics and show that this allows useful parameter velocity information to persist from the past trajectory in directions that are weakly informed, while well-informed parameter velocity directions continue to follow the Dirac–Frenkel dynamics. We prove that the inertial formulation yields well-posed parameter dynamics and provide a posteriori error bounds. After time discretization, the method requires the solution of the same type of regularized linear least-squares problem as standard Dirac–Frenkel dynamics, but with the previous velocity appearing as an anchor. Numerical experiments demonstrate the increased robustness obtained with inertia.

1 Introduction Let u̇(t) = F (u(t)),

u(0) = u0 ,

(1)

be an evolution problem on a Hilbert space H. We approximate the solution by a nonlinear parametrization û(t) = Φ(θ(t)), θ(t) ∈ Θ, (2) where Θ is a p-dimensional Hilbert space and Φ : Θ → H. For example, Φ could be a neural network and θ(t) the weights of the neural network. The Dirac–Frenkel variational principle determines the parameter velocity θ̇(t) by asking that the tangent vector J(θ(t))θ̇(t), with J(θ(t))w := DΦ(θ(t))[w], be the best instantaneous approximation of the vector field F (Φ(θ(t))) in H. This point of view is classical in quantum dynamics and dynamical lowrank approximation [11, 16, 26, 22, 17] and it has become increasingly relevant for nonlinear parametrizations given by neural networks such as Neural Galerkin schemes [6] and other techniques [13, 2, 31, 19, 15, 14, 27, 1, 20, 32, 10, 5, 18].

1

A central difficulty is that the Dirac–Frenkel principle determines the function-space tangent vector J(θ)θ̇(t), but the corresponding parameter velocity θ̇(t) need not be unique or well conditioned. For redundant parametrizations, such as neural networks or Gaussian mixtures, the Jacobian J(θ(t)) may have a nontrivial kernel, so different parameter velocities θ̇(t) produce the same function-space velocity J(θ(t))θ̇(t). Even if the Jacobian is mathematically full rank, small singular values of J(θ(t)) make the instantaneous problem of fitting F (Φ(θ(t))) by a tangent vector J(θ(t))θ̇(t) ill-conditioned. Components of the desired tangent vector that can be matched only through small singular directions of J(θ(t)) then require large, unstable coefficients in the parameter velocity. Thus, while the evolution in function space may be well determined, the induced dynamics in parameter space are non-unique or ill-conditioned. This phenomenon is called tangent space collapse [32] or matrix singularity issue [21] and is widely recognized; see, e.g., [25, 4, 15, 14, 9]. One remedy is to regularize the instantaneous Dirac–Frenkel problem so that a unique parameter velocity is selected. The work [14] provides a detailed analysis of Tikhonovregularized Dirac–Frenkel dynamics. Tikhonov regularization and related truncated singularvalue rules select a velocity by damping or removing directions that are weakly informed by the Jacobian J(θ(t)) at current parameter θ(t). Another line of work builds on randomization to use partial or compressed information from the current Jacobian. For example, one may evaluate the Jacobian action J(θ(t))v only at randomly sampled spatial points, or apply a random sketch matrix before solving for the velocity [4, 24, 12, 7]. These randomized methods can reduce the cost of forming and solving the least-squares problem, and in some cases improve conditioning. Other works propose to re-train the neural network [15] or use different dynamics than Dirac-Frenkel dynamics [23, 8, 32, 9]. In most of these approaches, the velocity is still based on the current time step only, through a regularizer, truncation, randomization rule, or modified local problem. Closest to our work is the Dirac–Frenkel–Onsager approach of [29], which interprets the non-uniqueness as gauge freedom. That method preserves the instantaneous Dirac–Frenkel residual minimization and uses an additional Onsager-type variable to select how the parameters move in nullspace directions. The present paper takes a different route. We introduce Dirac–Frenkel dynamics with inertia (DFI), which replace the instantaneous selection of a parameter velocity by an evolution equation for the velocity itself. The velocity is driven toward satisfying the current Dirac– Frenkel condition, but it is not forced to forget its past at every time step. This distinguishes DFI from the Dirac–Frenkel–Onsager approach of [29], which preserves the instantaneous Dirac–Frenkel residual minimization and uses history only to resolve gauge freedom in nullspace directions. DFI instead applies inertia to the full parameter velocity. As a result, DFI behaves like the usual Dirac–Frenkel dynamics in directions that are weakly informed by the current Jacobian, directions whose instantaneous correction is strongly regularized, or directions seen only through sketched information. We show that, once the initial velocity is fixed, DFI yields a well-defined evolution in parameter space and improves robustness by using past trajectory information when the current least-squares problem is weakly informative. After time discretization, a DFI time step leads to a previous-velocity-anchored least-squares solve, which requires the same type of regularized least-squares solve as standard regularized Dirac–Frenkel methods, with the previous velocity appearing as an anchor. The theory and

2

experiments below show that the inertial memory improves robustness precisely when the instantaneous least-squares problem is ill-conditioned or otherwise weakly informative.

2 Preliminaries and problem formulation We briefly recap the Dirac-Frenkel variational principle and discuss that parameter dynamics can be under-determined by it. 2.1 Dirac–Frenkel variational principle We work with a real Hilbert state space H and a p-dimensional real Hilbert parameter space Θ. Let Φ : Θ → H be of class C 2 , let F : H → H be the vector field of (1), and write û(t) = Φ(θ(t)) along a trajectory. In the following, when discussing the instantaneous Dirac–Frenkel problem at a fixed state, we suppress the time argument and, for θ ∈ Θ, set û = Φ(θ),

J(θ)w := DΦ(θ)[w] (w ∈ Θ),

f (θ) := F (û) ∈ H.

(3)

Note that since Θ is a finite-dimensional Hilbert vector space, we identify each tangent space with Θ. Thus parameter velocities are in Θ. At fixed θ, define the function-space Dirac–Frenkel defect by ˙ := 1 ∥û˙ − F (û)∥2 . E(û) H 2

(4)

The Dirac–Frenkel principle selects a velocity by minimizing (4) over all velocities û˙ in the range Im(DΦ(θ)), which is the tangent space at the representative θ. At the parameter level, using the chain rule û˙ = J(θ)θ̇, this corresponds to selecting the parameter velocity θ̇ θ̇ ∈ arg min Eθ (v) v∈Θ

(5)

by minimizing the parameter-space Dirac–Frenkel defect 1 Eθ (v) := E(J(θ)v) = ∥J(θ)v − f (θ)∥2H . 2

(6)

2.2 Non-unique or ill-conditioned parameter dynamics We now discuss that the Dirac–Frenkel dynamics impose dynamics in function space that can lead to under-determined parameter dynamics. 2.2.1 Non-unique parameter dynamics For θ ∈ Θ, let S0 (θ) := arg min Eθ (v). v∈Θ

(7)

All elements of S0 (θ) are (unregularized) Dirac–Frenkel velocities at θ that minimize the defect (6). If v̄(θ) ∈ S0 (θ) is any fixed minimizer, then
S0 (θ) = v̄(θ) + ker J(θ) = v̄(θ) + ker J(θ)∗ J(θ) . (8) 3

Thus, the ambiguity lies in directions of the velocity that do not change the tangent vector in H. 2.2.2 Tikhonov-regularized Dirac–Frenkel dynamics One way to select a unique velocity is regularization. The work [14] introduces the Tikhonovregularized Dirac–Frenkel functional 1 ε2 Eθε (v) := ∥J(θ)v − f (θ)∥2H + ∥v∥2Θ , 2 2

v ∈ Θ,

(9)

with regularization parameter ε > 0, and set θ̇ = v̄ε (θ) with v̄ε (θ) given by v̄ε (θ) = arg min Eθε (v).

(10)

v∈Θ

This yields the regularized Dirac–Frenkel dynamics with θ̇ = v̄ε (θ) ,

v̄ε (θ) = Mε (θ)−1 g(θ),

(11)

g(θ) := J(θ)∗ f (θ) = J(θ)∗ F (Φ(θ)) ,

(12)

where Mε (θ) := J(θ)∗ J(θ) + ε2 I,

with the adjoint J(θ)∗ of J(θ) with respect to the inner products of H and Θ. 2.2.3 Tikhonov regularization is a shrinkage rule To see how Tikhonov regularization acts on parameter directions that are weakly or not at all informed by the Jacobian, freeze θ and let qi ∈ Θ be right singular vectors of J(θ) with singular values σi ≥ 0. Let us now consider the velocity in these singular directions by setting v̄iε := ⟨v̄ε (θ), qi ⟩Θ ,

gi := ⟨g(θ), qi ⟩Θ .

The Tikhonov solution (11) decouples in this singular-vector basis as (σi2 + ε2 )v̄iε = gi .

(13)

For σi > 0, writing J(θ)qi = σi ui , this gives v̄iε =

σi σi2 ⟨u , f (θ)⟩ = i H σi2 + ε2 σi2 + ε2



 1 ⟨ui , f (θ)⟩H . σi | {z }

(14)

v̄i0

Thus each positive singular direction is damped relative to the unregularized pseudoinverse coefficient v̄i0 by the factor σi2 /(σi2 + ε2 ). This damping is strongest for small singular values, precisely the directions that are only weakly informed by the current Jacobian. In exact null directions, qi ∈ ker J(θ) and hence σi = 0, the forcing coefficient satisfies gi = ⟨J(θ)∗ f (θ), qi ⟩Θ = ⟨f (θ), J(θ)qi ⟩H = 0, 4

so that the scalar Tikhonov equation (13) reduces to ε2 v̄iε = 0 ,

(15)

because ε > 0. Thus, the Tikhonov-selected velocity has zero component in ker J(θ). In this sense, Tikhonov regularization is a static, instantaneous selection rule: it damps weakly informed directions and removes exactly uninformed directions. A truncated pseudoinverse used in, e.g., [6] is even more abrupt. Directions whose singular values fall below the truncation threshold are simply ignored and their velocity components are set to zero. In summary, both mechanisms—Tikhonov regularization and truncated singular value decomposition (SVD)—resolve the parameter non-uniqueness by suppressing velocity directions that are not sufficiently visible to the current Jacobian. In particular, damping occurs through an instantaneous shrink-or-delete rule so that Tikhonov and truncated SVD regularization are an instantaneous shrinkage rules in parameter space.

3 Dirac–Frenkel dynamics with inertia We now introduce Dirac–Frenkel dynamics with inertia (DFI). Through inertia, the parameter velocity carries information from the past of the trajectory, so motion can persist in weakly informed directions even when these directions are singular or close to singular in the instantaneous Jacobian. 3.1 Dirac–Frenkel dynamics with inertia For fixed θ ∈ Θ and current velocity θ̇ ∈ Θ, we choose an acceleration θ̈ ∈ Θ by Onsager’s principle applied to the unregularized defect energy (6), namely θ̈ ∈ arg min R0θ,θ̇ (a) := a∈Θ

τ2 ∥a∥2Θ + ⟨∇Eθ (θ̇), a⟩Θ , 2

τ > 0,

(16)

where ∇Eθ denotes the gradient with respect to the Hilbert-space inner product on Θ. The above selects an acceleration θ̈ that minimizes the sum of two terms: a quadratic penalty on acceleration, which prevents rapid changes in the velocity, and the rate of change of the unregularized Dirac–Frenkel defect energy in the direction of the acceleration, which promotes minimization of the defect energy. The parameter τ plays the role of a mass in parameter space, which prevents the parameter velocity from changing direction abruptly. Because R0θ,θ̇ is strictly convex, the minimizer is unique and solves
τ 2 θ̈ + J(θ)∗ J(θ)θ̇ − f (θ) = 0.

(17)

We therefore obtain the DFI system θ̇ = v,

τ 2 v̇ + J(θ)∗ J(θ)v = J(θ)∗ f (θ).

5

(18)

3.2 Tikhonov-regularized Dirac–Frenkel dynamics with inertia The same inertial mechanism can be combined with the Tikhonov-regularized Dirac–Frenkel functional (9). For fixed θ ∈ Θ and current velocity θ̇ ∈ Θ, we choose an acceleration θ̈ ∈ Θ by Onsager’s principle applied to (9): θ̈ ∈ arg min Rεθ,θ̇ (a) := a∈Θ

τ2 ∥a∥2Θ + ⟨∇θ̇ Eθε (θ̇), a⟩Θ , 2

τ > 0.

(19)

Because Rεθ,θ̇ is strictly convex, the minimizer is unique and is characterized by τ 2 θ̈ + ∇θ̇ Eθε (θ̇) = 0.

(20)

∇Eθε (v) = J(θ)∗ J(θ)v − f (θ) + ε2 v

(21)

Using and the definitions in (12), we obtain the coupled DFI system θ̇ = v,

τ 2 v̇ + Mε (θ)v = g(θ) ,

(22)

which can be written in second-order form as τ 2 θ̈ + Mε (θ)θ̇ = g(θ).

(23)

For fixed θ, the v-equation in (22) is the gradient flow of Eθε in the velocity variable v, with mobility τ −2 IΘ . If ε = 0, the relaxation is toward the affine set of unregularized Dirac–Frenkel minimizers. In the frozen-θ dynamics, if ε > 0, the gradient flows relaxes toward the unique regularized Dirac–Frenkel velocity. We note that the second-order system (23) is analogous to heavy-ball dynamics used in optimization; see, e.g., [28, 3]. 3.3 Interpretation of DFI Let us now interpret DFI direction by direction, in direct analogy with the regularization rules discussed in Section 2.2. Freeze θ and choose an orthonormal basis {qi }pi=1 of Θ consisting of eigenvectors of J(θ)∗ J(θ), J(θ)∗ J(θ)qi = σi2 qi ,

σi ≥ 0.

with singular values σi ≥ 0. Set vi := ⟨v, qi ⟩Θ ,

gi := ⟨g(θ), qi ⟩Θ .

Since Mε (θ)qi = (σi2 + ε2 )qi , the qi -component of (22) is τ 2 v̇i + (σi2 + ε2 )vi = gi ,

6

(24)

which is in stark contrast to Tikhonov and other instantaneous regularizers that implement an instantaneous algebraic selection such as (13). We can further rewrite (24) as τ 2 v̇i + (σi2 + ε2 )(vi − v̄iε ) = 0 to make explicit the connection to v̄iε given by (13). Thus, for frozen θ, the Tikhonov coefficient v̄iε is not imposed immediately; it is approached via relaxation. The relaxation time in this direction is τ 2 /(σi2 + ε2 ). In directions where τ 2 ≪ σi2 + ε2 , the velocity rapidly tracks the Tikhonov velocity v̄ε given by (11), so DFI behaves like the Tikhonov-regularized DF dynamics. In directions where τ 2 ≫ σi2 + ε2 , the relaxation is slow and the velocity is mainly transported by inertia. For σi > 0, the Tikhonov rule leads to the damped velocity coefficient given in (14). In DFI, the damped velocity direction (14) is only the asymptotic target of the frozen-θ velocity dynamics. In particular, if the direction is weakly informed, the velocity need not be reset immediately to this small Tikhonov value; it can retain motion from the past trajectory controlled by τ . The contrast of DFI to Tikhonov and truncated regularization is most pronounced in exact null directions. If qi ∈ ker J(θ), then Tikhonov regularization selects v̄iε = 0 (see (15)), whereas DFI gives τ 2 v̇i + ε2 vi = 0. In the frozen-θ model, if ε = 0, then the velocity in nullspace directions is maintained and not changed by DFI. For ε > 0, it is not removed instantaneously but decays exponentially. Thus, in the frozen-θ interpretation, DFI permits velocity components in parameter directions that are weakly seen by J(θ), including exact kernel directions, to persist over a relaxation time scale rather than being removed instantaneously. This makes DFI history-aware because the current residual still corrects the velocity in directions resolved by the Jacobian, while components in weakly resolved directions are damped dynamically rather than eliminated by a pointwise algebraic rule. 3.4 Well-posedness of DFI We now turn to the well-posedness of the DFI system (22). We first establish local existence and uniqueness, and then global existence under a growth condition on the force map g defined in (12). Proposition 1 (Local existence and uniqueness of DFI solution). Recall that Φ ∈ C 2 (Θ; H) and assume that F : H → H is locally Lipschitz. Then, for τ > 0, ε ≥ 0, and initial datum (θ0 , v 0 ) ∈ Θ × Θ, there exists a time Tmax ∈ (0, ∞] and a unique solution (θ(t), v(t)) of the DFI system (22) on [0, Tmax ) with (θ(0), v(0)) = (θ0 , v 0 ). Moreover, this solution is maximal in the sense that it cannot be extended to any larger interval [0, T ′ ) with T ′ > Tmax . Proof. Since Θ is finite-dimensional, the DFI system is an ordinary differential equation on the finite-dimensional phase space Θ×Θ. Since Φ ∈ C 2 , the map θ 7→ J(θ) = DΦ(θ) is locally Lipschitz. Since F is locally Lipschitz and Φ is C 2 , the composition θ 7→ f (θ) = F (Φ(θ)) is

7

locally Lipschitz as well. Therefore g(θ) is locally Lipschitz, because

g(θ1 ) − g(θ2 ) = J(θ1 )∗ − J(θ2 )∗ f (θ1 ) + J(θ2 )∗ f (θ1 ) − f (θ2 ) . On each bounded neighborhood in Θ, the maps J, J ∗ , and f are locally bounded, and J and f are locally Lipschitz; hence the right-hand side is bounded by a constant times ∥θ1 − θ2 ∥Θ . Similarly, Mε (θ) = J(θ)∗ J(θ) + ε2 I is locally Lipschitz and thus the map (θ, v) 7→ Mε (θ)v is locally Lipschitz on Θ × Θ. Hence the phase-space vector field
Gτ,ε (θ, v) := v, τ −2 (J(θ)∗ f (θ) − Mε (θ)v) is locally Lipschitz on Θ × Θ. The classical Cauchy–Lipschitz theorem therefore yields a unique maximal local solution. Proposition 2 (Global existence of DFI solution). Assume the hypotheses of Proposition 1 and, in addition, that the map g defined in (12) satisfies ∥g(θ)∥Θ ≤ Cg (1 + ∥θ∥Θ ),

θ ∈ Θ,

(25)

for some constant Cg > 0. Then every maximal solution of (22) is global. Proof. Let (θ(t), v(t)) be a maximal local solution and define X(t) := ∥θ(t)∥2Θ + τ 2 ∥v(t)∥2Θ ,

(26)

which is differentiable because by the Cauchy–Lipschitz theorem, the local maximal solution is continuously differentiable. Along the solution we compute Ẋ(t) = 2⟨θ(t), v(t)⟩Θ + 2⟨g(θ(t)), v(t)⟩Θ − 2∥J(θ(t))v(t)∥2H − 2ε2 ∥v(t)∥2Θ ≤ 2⟨θ(t), v(t)⟩Θ + 2⟨g(θ(t)), v(t)⟩Θ Now use Cauchy-Schwarz to obtain ⟨θ(t), v(t)⟩Θ ≤ ∥θ(t)∥Θ ∥v(t)∥Θ ,

⟨g(θ(t)), v(t)⟩Θ ≤ ∥g(θ(t))∥Θ ∥v(t)∥Θ .

Using Young’s inequality 2ab ≤ a2 + b2 , we obtain Ẋ(t) ≤ ∥θ(t)∥2Θ + 2∥v(t)∥2Θ + ∥g(θ(t))∥2Θ . Recall the growth condition (25) to obtain ∥g(θ(t))∥2Θ ≤ Cg2 (1 + ∥θ(t)∥Θ )2 ≤ 2Cg2 + 2Cg2 ∥θ(t)∥2Θ ,

8

where we used (a + b)2 ≤ 2a2 + 2b2 in the last step, and bound Ẋ(t) as Ẋ(t) ≤ 2Cg2 + (1 + 2Cg2 )∥θ(t)∥2Θ + 2∥v(t)∥2Θ .

(27)

The terms involving J(θ)∗ J(θ) and ε2 I are dissipative in this estimate; therefore no growth assumption on J(θ) is needed for this particular global bound. To write the right-hand side of (27) in terms of X(t) notice that both terms in (26) are non-negative so that ∥θ(t)∥2Θ ≤ X(t) , τ 2 ∥v(t)∥2Θ ≤ X(t) . and thus

 2 Ẋ(t) ≤ 2Cg2 + 1 + 2Cg2 + 2 X(t). τ Gronwall’s lemma therefore yields a bound on X(t) on every finite time interval. Thus (θ(t), v(t)) remains bounded on every finite time interval. Since the vector field is locally Lipschitz on the finite-dimensional phase space Θ × Θ, the standard continuation theorem for ODEs implies that the maximal existence time is infinite [30, Theorem 2.17]. Proposition 2 shows that the nonuniqueness in the parameters given by (unregularized) Dirac–Frenkel dynamics (5) is removed by DFI. Instead of selecting, independently at each θ, one element of the affine set S0 (θ), DFI treats the parameter velocity as part of the state. Once the initial phase point (θ0 , v 0 ) is fixed, Proposition 1 yields a unique parameter trajectory even for ε = 0, and Proposition 2 shows that this trajectory exists for all times under the stated growth condition. 3.5 A posteriori error analysis of DFI in continuous time Throughout this section we fix τ > 0 and ε > 0 and let (θ(t), v(t)) be a sufficiently smooth solution of (22). For each θ, recall the velocity v̄ε (t) given by (10) and define the corresponding ū˙ ε (t) := J(θ(t))v̄ε (t). We consider the projection defect as δε (t)2 := ∥ū˙ ε (t) − f (θ(t))∥2H + ε2 ∥v̄ε (t)∥2Θ ,

(28)

which is also used in [14, Section 3.1]. Along the DFI trajectory, define the relaxation lag by
˙ − ū˙ ε (t). r(t) := J(θ(t)) v(t) − v̄ε (t) = û(t) This motivates introducing the relaxation defect as ρε (t)2 := ∥r(t)∥2H + ε2 ∥v(t) − v̄ε (t)∥2Θ .

(29)

The following proposition shows that the instantaneous total defect can be decomposed into the projection defect δε and the relaxation defect ρε . Proposition 3 (Instantaneous defect decomposition). For every t in the interval of existence, ˙ − F (û(t))∥2 + ε2 ∥v(t)∥2 = δε (t)2 + ρε (t)2 . ∥û(t) H Θ 9

(30)

Proof. Decompose v(t) as v(t) = v̄ε (t) + (v(t) − v̄ε (t)) = v̄ε (t) + w(t) , to write ˙ − F (û(t))∥2 + ε2 ∥v(t)∥2 = ∥J(θ(t))v(t) − f (θ(t))∥2 + ε2 ∥v(t)∥2 ∥û(t) H Θ H Θ = ∥J(θ(t))v̄ε (t) − f (θ(t)) + J(θ(t))w(t)∥2H + ε2 ∥v̄ε (t) + w(t)∥2Θ = (∥J(θ(t))v̄ε (t) − f (θ(t))∥2H + ε2 ∥v̄ε (t)∥2Θ ) + (∥J(θ(t))w(t)∥2H + ε2 ∥w(t)∥2Θ ) + 2⟨J(θ(t))v̄ε (t) − f (θ(t)), J(θ(t))w(t)⟩H + 2ε2 ⟨v̄ε (t), w(t)⟩Θ = δε (t)2 + ρε (t)2 + 2⟨J(θ(t))v̄ε (t) − f (θ(t)), J(θ(t))w(t)⟩H + 2ε2 ⟨v̄ε (t), w(t)⟩Θ . We now show that the cross terms vanish, which then leads to the decomposition (30). Consider ⟨J(θ(t))v̄ε (t) − f (θ(t)), J(θ(t))w(t)⟩H = ⟨J(θ(t))∗ (J(θ(t))v̄ε (t) − f (θ(t))), w(t)⟩Θ and thus ⟨J(θ(t))v̄ε (t) − f (θ(t)), J(θ(t))w(t)⟩H + ε2 ⟨v̄ε (t), w(t)⟩Θ = ⟨J(θ(t))∗ (J(θ(t))v̄ε (t) − f (θ(t))) + ε2 v̄ε (t), w(t)⟩Θ . (31) The left argument of the inner product in (31) are the first-order optimality conditions of the objective (9), which is minimized by v̄ε (t) and thus the left argument of the inner product of (31) is zero and the cross terms vanish. Theorem 3.1 (Continuous a posteriori error bound). We assume that the vector field F satisfies the one-sided Lipschitz estimate ⟨u − u e, F (u) − F (e u)⟩H ≤ ℓ∥u − u e∥2H ,

u, u e ∈ H,

(32)

for some ℓ ∈ R. Let u(t) be a sufficiently smooth solution of (1) on [0, T ], and let û(t) = Φ(θ(t)) be the DFI approximation. Then, for every t ∈ [0, T ], Z t
1/2 ℓt ∥û(t) − u(t)∥H ≤ e ∥û(0) − u(0)∥H + eℓ(t−s) δε (s)2 + ρε (s)2 ds. (33) 0

In particular, if û(0) = u(0), then Z t ∥û(t) − u(t)∥H ≤

eℓ(t−s) δε (s)2 + ρε (s)2

1/2

0

Proof. Define the error e(t) := û(t) − u(t). Because u solves (1),
˙ − F (û(t)) . ė(t) = F (û(t)) − F (u(t)) + û(t) 10

ds.

(34)

Now consider 1d ˙ − F (û(t))⟩H . ∥e(t)∥2H = ⟨e(t), ė(t)⟩H = ⟨û(t) − u(t), F (û(t)) − F (u(t))⟩H + ⟨e(t), û(t) 2 dt Using (32), we obtain 1d ˙ − F (û(t))⟩H . ∥e(t)∥2H ≤ ℓ∥e(t)∥2H + ⟨e(t), û(t) 2 dt Applying Cauchy-Schwarz to the second term and using (30) to obtain ˙ − F (û(t))∥2 ≤ δε (t)2 + ρε (t)2 , ∥û(t) H leads to

1/2 1d ∥e(t)∥2H ≤ ℓ∥e(t)∥2H + ∥e(t)∥H δε (t)2 + ρε (t)2 . 2 dt Whenever ∥e(t)∥H ̸= 0, division by ∥e(t)∥H yields
1/2 d ∥e(t)∥H ≤ ℓ∥e(t)∥H + δε (t)2 + ρε (t)2 , dt and when ∥e(t)∥H = 0 then we restart the same argument from time t when ∥e(t)∥H ̸= 0. By continuity this differential inequality extends to all t ∈ [0, T ], and Gronwall’s lemma gives (33). The second estimate is the specialization to exact initial data. Compared with the Tikhonov-regularized Dirac–Frenkel estimate [14, Section 3.1], this bound separates the instantaneous defect along the DFI trajectory into the regularized projection defect δε and the relaxation defect ρε ; however, it should not be read as the same estimate with an extra nonnegative term added along the same path. The quantities δε and ρε are evaluated along the DFI trajectory θ(t), whereas the Tikhonov-regularized Dirac-Frenkel defect is evaluated along the trajectory generated by θ̇ = v̄ε (θ). These trajectories may visit different regions of the parameter space Θ, even when their function-space approximations are close. In particular, the inertial dynamics can move through kernel and near-kernel directions and thereby in principle can sample representatives for which the instantaneous least-squares problem is better conditioned. Thus, δε (t) may be smaller than the corresponding regularized Dirac-Frenkel defect along its own trajectory, but this is a trajectory-dependent effect rather than a pointwise comparison of the two bounds.

4 Euler time discretization of DFI We now turn to the time discretization of DFI. We derive a time-discrete scheme based on a semi-implicit Euler discretization of (22) that treats the velocity relaxation implicitly and the parameter update explicitly. This choice yields a velocity update through an implicit solve, while keeping the parameter update explicit; equivalently, each step becomes a previousvelocity-anchored least-squares problem that exposes the inertial memory mechanism.

11

4.1 Euler time discretization Let hk > 0 be the step size and tk+1 = tk + hk . Given (θk , vk ), we set ûk := Φ(θk ),

Jk := J(θk ),

fk := f (θk ) = F (ûk ).

(35)

We discretize the DFI system (22) with implicit Euler in v and explicit Euler in θ, θk+1 =θk + hk vk+1 , v − vk k+1 τ2 = − Jk∗ (Jk vk+1 − fk ) − ε2 vk+1 . hk

(36)

It is convenient to further introduce ηk2 := ε2 +

τ2 ∈ (0, ∞), hk

βk :=

τ2 ∈ (0, 1], ε2 hk + τ 2

Mη,k := Jk∗ Jk + ηk2 I.

(37)

Note that 0 < βk < 1 when ε > 0, while βk = 1 when ε = 0. A direct computation shows that the discretized system (36) can be written as θk+1 = θk + hk vk+1 ,
−1 ∗ vk+1 = βk vk + Mη,k Jk fk − βk Jk vk .

(38)

The formulation (38) mirrors the continuous direction-by-direction interpretation of DFI. If qk,i is a right singular vector of Jk with singular value σk,i , and gk,i := ⟨Jk∗ fk , qk,i ⟩Θ ,

vk,i := ⟨vk , qk,i ⟩Θ , then ⟨vk+1 , qk,i ⟩Θ =

gk,i βk ηk2 + 2 + η2 2 + η 2 vk,i . σk,i σk,i k k

Thus, for frozen θk , the discrete update has the same type of direction-dependent balance 2 ≫ η 2 , the memory factor as the continuous DFI dynamics in the sense that when σk,i k 2 2 2 βk ηk /(σk,i + ηk ) is small, so the update is dominated by the current least-squares information. 2 ≪ η 2 , this factor is close to β , so the update retains most of the damped previous When σk,i k k velocity coefficient. On ker Jk , this reduces to ⟨vk+1 , qk,i ⟩Θ = βk ⟨vk , qk,i ⟩Θ , so null-space motion is conserved when ε = 0 and damped otherwise. 4.2 Variational characterization of time-discrete DFI Let us now give a variational characterization of the DFI step (38), which will be useful for the further interpretation of the time-discrete DFI scheme and error analysis. Proposition 4 (Discrete variational characterization). For fixed k, the velocity update vk+1 defined by the update (38) is the unique minimizer of 1 Jk (w) := ∥Jk w − fk ∥2H + ηk2 ∥w − βk vk ∥2Θ . 2

(39)

vk+1 = arg min Jk (w).

(40)

That is, w∈Θ

12

Proof. Since ηk2 > 0, the functional Jk is strongly convex and therefore has a unique minimizer. Its first-order optimality condition is Jk∗ (Jk w − fk ) + ηk2 (w − βk vk ) = 0. For w = vk+1 this becomes (Jk∗ Jk + ηk2 I)vk+1 = Jk∗ fk + βk ηk2 vk . Using βk ηk2 = τ 2 /hk , this becomes 

τ Jk∗ Jk + ε2 IΘ +

2

hk



IΘ vk+1 = Jk∗ fk +

τ2 vk , hk

which is precisely the velocity equation in (36). Thus the unique minimizer of Jk is the velocity given by the semi-implicit Euler step. The variational characterization shows that the DFI step is a (damped-)previous-velocityanchored least-squares problem. The new velocity is chosen to reduce the instantaneous Dirac–Frenkel residual while remaining close to the damped previous velocity βk vk . Thus ηk2 controls the strength of the anchoring, whereas βk controls how much of the previous velocity is retained in the anchor. 4.3 Algorithm With the change of variables z = w − βk vk , the same minimization problem can be written as zk+1 = arg min∥Jk z − (fk − βk Jk vk )∥2H + ηk2 ∥z∥2Θ , z∈Θ

(41)

and the velocity is recovered by vk+1 = βk vk + zk+1 . The shifted form (41) is convenient for implementation. It requires one solve of a standard Tikhonov-regularized least-squares problem for the correction zk+1 , while the retained part βk vk carries the inertial memory of the method. The resulting semi-implicit Euler DFI algorithm is summarized in Algorithm 1.

5 A posteriori error analysis of Euler-discretized DFI We now provide an a posteriori analysis of the Euler-discretized DFI scheme (38). Throughout this section we fix τ > 0 and ε > 0. Note that the restriction ε > 0 is used in the velocity estimates below, where 1 − βk > 0 is required. The case ε = 0 has βk = 1 and is therefore not covered by the estimates involving (1 − βk )−1 .

13

Algorithm 1 Semi-implicit Euler Dirac–Frenkel dynamics with inertia (DFI) Require: Initial data θ0 ∈ Θ, v0 ∈ Θ, step sizes hk > 0, parameters τ > 0 and ε ≥ 0 Ensure: Iterates (ûk )k≥0 1: for k = 0, 1, 2, . . . do 2: Set ûk = Φ(θk ), Jk = J(θk ), and fk = F (ûk ). 3: Set τ2 τ2 ηk2 = ε2 + , . βk = 2 hk ε hk + τ 2 4:

Solve the regularized least-squares problem zk+1 = arg min∥Jk z − (fk − βk Jk vk )∥2H + ηk2 ∥z∥2Θ . z∈Θ

Set vk+1 = βk vk + zk+1 . 6: Set θk+1 = θk + hk vk+1 . 7: end for 5:

5.1 Local error bound We start by making stronger assumptions on the problem than in previous sections. We denote the flow of (1) as φt : H → H,

u(t) = φt (u(0)).

(42)

Assumption 5.1 (Flow stability and regularity). Fix T > 0. Assume that there exist constants ℓ ∈ R, CΦ ≥ 0, cΦ ≥ 0, and Ca ≥ 0 such that the following hold for all relevant states. (i) The exact flow satisfies ∥φt (u) − φt (e u)∥H ≤ eℓt ∥u − u e∥H ,

0 ≤ t ≤ T.

(43)

(ii) The second derivative of Φ is bounded, relative to the frozen Jacobian J(θ)w := DΦ(θ)[w] at the base point, as follows for all relevant θ, ζ, ξ: ∥D2 Φ(θ + sζ)[ξ, ξ]∥H ≤ CΦ ∥J(θ)ξ∥2H + cΦ ∥ξ∥2Θ ,

0 ≤ s ≤ 1.

(44)

(iii) For every relevant initial state w and the corresponding exact trajectory y(s) = φs (w), one has ∥ÿ(s)∥H ≤ Ca for all times s for which the trajectory is used below. We denote the discrete-time defect at time step k as ∆2k := ∥Jk vk+1 − fk ∥2H + ηk2 ∥vk+1 − βk vk ∥2Θ

(45)

This is twice the normalized anchored least-squares objective Jk if the latter is written with the conventional factor 1/2. 14

Lemma 1 (One-step local error). Under Assumption 5.1, for every time step k, the local error of one Euler step û+ := Φ(θk + hk vk+1 ) . (46) satisfies CΦ βk2 ηk2 2 cΦ Ca 2 hk ∥vk ∥2Θ + h2k ∥vk+1 ∥2Θ + h . (47) 4 2 2 k

∥û+ − φhk (ûk )∥H ≤ hk ∆k + CΦ h2k ∥fk ∥2H +

Proof. Let yk (s) := φs (ûk ) for 0 ≤ s ≤ hk . Then yk (0) = ûk and, since yk solves the exact evolution equation, ẏk (0) = F (ûk ) = fk . Taylor’s formula with integral remainder gives Z hk φhk (ûk ) = yk (hk ) = yk (0) + hk ẏk (0) + (hk − s)ÿk (s) ds. (48) 0

Hence φhk (ûk ) = ûk + hk fk + rkflow , where

(49)

Z hk

rkflow :=

(hk − s)ÿk (s) ds.

(50)

0

By Assumption 5.1(iii), ∥rkflow ∥H ≤

Z hk (hk − s)∥ÿk (s)∥H ds 0

Z hk ≤ Ca

(hk − s) ds = 0

Ca 2 h . 2 k

(51)

Next, Taylor’s formula for Φ applied to the curve s 7→ θk + shk vk+1 , which parametrizes the straight line segment in Θ from θk to θk+1 = θk + hk vk+1 , gives û+ = Φ(θk + hk vk+1 ) = Φ(θk ) + hk DΦ(θk )[vk+1 ] + h2k

Z 1

(1 − s)D2 Φ(θk + shk vk+1 )[vk+1 , vk+1 ] ds.

(52)

0

Since ûk = Φ(θk ) and Jk vk+1 = DΦ(θk )[vk+1 ] (see (3)), this becomes

where rkΦ := h2k

Z 1

û+ = ûk + hk Jk vk+1 + rkΦ ,

(53)

(1 − s)D2 Φ(θk + shk vk+1 )[vk+1 , vk+1 ] ds.

(54)

0

Applying Assumption 5.1(ii) with θ = θk , ζ = hk vk+1 , ξ = vk+1 yields Z 1 Φ 2 ∥rk ∥H ≤ hk (1 − s)∥D2 Φ(θk + shk vk+1 )[vk+1 , vk+1 ]∥H ds 0 Z 1
≤ h2k (1 − s) CΦ ∥Jk vk+1 ∥2H + cΦ ∥vk+1 ∥2Θ ds 0

h2 = k CΦ ∥Jk vk+1 ∥2H + cΦ ∥vk+1 ∥2Θ . 2 15

(55)

Let us now bound the state velocity Jk vk+1 . Using the definitions −1 ∗ v̄η,k := Mη,k Jk fk ,

Mη,k := Jk∗ Jk + ηk2 I,

(56)

where v̄η,k is the instantaneous Tikhonov-regularized Dirac–Frenkel velocity at θk with effective regularization parameter ηk , we obtain with (38), −1 Jk vk+1 = Jk v̄η,k + βk ηk2 Jk Mη,k vk .

(57)

First, we bound the ηk -regularized instantaneous Dirac–Frenkel Jk v̄η,k of Jk vk+1 . By the push-through identity Jk (Jk∗ Jk + ηk2 IΘ )−1 Jk∗ = Jk Jk∗ (Jk Jk∗ + ηk2 IH )−1 , we have Jk v̄η,k = Jk Jk∗ (Jk Jk∗ + ηk2 I)−1 fk .

(58)

The self-adjoint nonnegative operator Jk Jk∗ (Jk Jk∗ + ηk2 I)−1 has spectral values λ , λ + ηk2

λ ≥ 0,

which all lie in [0, 1]. Therefore, ∥Jk v̄η,k ∥H ≤ ∥fk ∥H .

(59)

Second, we bound the inertial-memory part of Jk vk+1 . Spectral calculus for the self-adjoint nonnegative operator Jk∗ Jk gives −1 ∥ηk2 Jk Mη,k vk ∥2H = ⟨ηk4 Jk∗ Jk (Jk∗ Jk + ηk2 I)−2 vk , vk ⟩Θ .

(60)

For every spectral value λ ≥ 0 of Jk∗ Jk , ηk4

ηk2 λ/ηk2 λ 2 = η ≤ , k 4 (λ + ηk2 )2 (1 + λ/ηk2 )2

(61)

because x/(1 + x)2 ≤ 1/4 for x ≥ 0. Hence −1 vk ∥2H ≤ ∥ηk2 Jk Mη,k

ηk2 ∥vk ∥2Θ . 4

(62)

Using (57) and the elementary inequality ∥a + b∥2H ≤ 2∥a∥2H + 2∥b∥2H , we get −1 ∥Jk vk+1 ∥2H ≤ 2∥Jk v̄η,k ∥2H + 2βk2 ∥ηk2 Jk Mη,k vk ∥2H ≤ 2∥fk ∥2H +

βk2 ηk2 ∥vk ∥2Θ . 2

(63)

We now subtract the exact-flow expansion (49) from the parametric expansion (53) and obtain with the triangle inequality, ∥ûk+1 − φhk (ûk )∥H ≤ hk ∥Jk vk+1 − fk ∥H + ∥rkΦ ∥H + ∥rkflow ∥H .

(64)

Inserting the bounds (55), (51), and (63), we obtain ∥ûk+1 − φhk (ûk )∥H ≤ hk ∥Jk vk+1 − fk ∥H + CΦ h2k ∥fk ∥2H + +

cΦ 2 Ca 2 hk ∥vk+1 ∥2Θ + h . 2 2 k

Finally, by using (45), we obtain the bound (47). 16

CΦ βk2 ηk2 2 hk ∥vk ∥2Θ 4 (65)

The bound (47) in Lemma 1 still depends on the norm of the velocity vk and vk+1 . The following lemma bounds these. Lemma 2 (Squared velocity recursion). Under the setup of Lemma 1, we have ηk ∥vk+1 − βk vk ∥Θ ≤ ∆k .

(66)

Consequently, ∥vk+1 ∥2Θ ≤ βk ∥vk ∥2Θ +

∆2k . (1 − βk )ηk2

Iterating this recursion yields     k k k Y X Y  ∥vk+1 ∥2Θ ≤  βj  ∥v0 ∥2Θ + βj  m=0

j=0

j=m+1

(67)

∆2m . 2 (1 − βm )ηm

(68)

Proof. Because ε > 0, we have 0 < βk < 1 and thus (66) follows directly from (45). Next define the velocity correction ζk := vk+1 − βk vk . (69) Then vk+1 = βk vk + ζk and taking the squared Θ-norm gives ∥vk+1 ∥2Θ = βk2 ∥vk ∥2Θ + 2βk ⟨vk , ζk ⟩Θ + ∥ζk ∥2Θ .

(70)

We now estimate the cross term. Since 0 < βk < 1, we have 1 − βk > 0. By Cauchy–Schwarz and Young’s inequality, s ! p  βk 2βk ⟨vk , ζk ⟩Θ ≤ 2βk ∥vk ∥Θ ∥ζk ∥Θ = 2 βk (1 − βk ) ∥vk ∥Θ ∥ζk ∥Θ 1 − βk ≤ βk (1 − βk )∥vk ∥2Θ +

βk ∥ζk ∥2Θ . 1 − βk

(71)

Substituting (71) into (70), we obtain ∥vk+1 ∥2Θ ≤ βk ∥vk ∥2Θ +

1 ∥ζk ∥2Θ , 1 − βk

(72)

where we used βk2 + βk (1 − βk ) = βk2 + βk − βk2 = βk ,

βk βk 1 − βk 1 +1= + = . 1 − βk 1 − βk 1 − βk 1 − βk

From (66) we obtain ∥ζk ∥2Θ ≤

∆2k . ηk2

and thus ∥vk+1 ∥2Θ ≤ βk ∥vk ∥2Θ + This proves (67). Iterating (67) gives (68). 17

∆2k . (1 − βk )ηk2

(73)

(74)

The local estimate in Lemma 1 still contains the velocity norms ∥vk ∥Θ and ∥vk+1 ∥Θ . These velocities are controlled by the squared velocity recursion (67), whose right-hand side contains the term ∆2k /((1 − βk )ηk2 ). One way to control this term would be to impose a step-size restriction of the form hk ∆k ≤ c∆ ε2 using the ε coming from the Tikhonov regularizer in (9). Indeed, since (1 − βk )ηk2 = ε2 , this condition implies ∆2k c∆ ≤ ∆k . 2 hk (1 − βk )ηk This is the same type of restriction as in Tikhonov-regularized Dirac–Frenkel dynamics, where the relevant regularization scale is ε2 ; see [14]. In the DFI scheme, however, the regularization scale at step k is not only ε2 , but rather ηk2 = ε2 + τ 2 /hk , as shown in (39) and, equivalently, in the shifted least-squares formulation (41). We therefore impose the DFI defect-control condition hk ∆k ≤ c∆ ηk2 .

(75)

This condition uses the effective regularization of the DFI velocity update, which contains both the Tikhonov contribution ε2 and the inertial contribution τ 2 /hk . Under (75), the velocity recursion term satisfies ∆2k ∆k ∆k c∆ = ≤ ∆k . 1 − βk ηk2 (1 − βk )hk (1 − βk )ηk2 Thus DFI allows us to use the larger effective regularization scale ηk2 , and therefore the weaker step-size restriction (75). The price is the explicit factor (1 − βk )−1 in the resulting velocity estimates. Since βk = τ 2 /(ε2 hk + τ 2 ), the regime βk ≈ 1 corresponds to weak damping of the previous velocity. So from a stability perspective, DFI is most useful in a moderate-memory regime where βk is large enough that the method benefits from inertial transport of the velocity, but not so close to one that the damping encoded by 1 − βk becomes ineffective. The following proposition show the resulting local error bound. Proposition 5 (Local error under defect control). Under the setup of Lemma 1, suppose that there exists a constant c∆ ≥ 0 such that (75) holds for all k. Define the squared velocity envelope by B02 := ∥v0 ∥2Θ , (76) and, for k ≥ 0, by  2 Bk+1 := 

k Y

 βj  ∥v0 ∥2Θ + c∆

k X

k Y

 m=0

j=0



 βj 

j=m+1

∆m . (1 − βm )hm

(77)

Then ∥vk ∥2Θ ≤ Bk2

for k ≥ 0,

(78)

CΦ βk2 ηk2 2 2 cΦ 2 2 Ca 2 hk Bk + hk Bk+1 + h . 4 2 2 k

(79)

and the local error of û+ given in (46) satisfies ∥û+ − φhk (ûk )∥H ≤ hk ∆k + CΦ h2k ∥fk ∥2H + 18

Proof. We start from the squared velocity growth estimate (68) and defect-control condition 2 and (75). Since hm > 0, ηm > 0, and ∆m ≥ 0, we may divide (75) by hm ηm ∆2m ∆m . ≤ c∆ 2 (1 − βm )ηm (1 − βm )hm

(80)

2 Substituting (80) into (68) and by the definition of Bk+1 given in (77), we obtain 2 ∥vk+1 ∥2Θ ≤ Bk+1 .

(81)

The corresponding estimate for vk holds because for k = 0, we have by definition of B02 = ∥v0 ∥2Θ that ∥v0 ∥2Θ ≤ B02 holds with equality, and for k ≥ 1, we use (81) with k − 1 in place of k to obtain ∥vk ∥2Θ ≤ Bk2 . Hence (78) holds. The bound (79) follows now from (47). Corollary 5.2 (Uniform-step local error under defect control). Assume the hypotheses of Proposition 5. In addition, suppose that hk = h for all k. Then ηk2 = η 2 = ε2 +

τ2 , h

βk = β =

τ2 . ε2 h + τ 2

(82)

Define ¯ k := max ∆j . ∆

(83)

0≤j≤k

Then ∥vk+1 ∥2Θ ≤ β k+1 ∥v0 ∥2Θ +

c∆ ¯ k. ∆ (1 − β)2 h

(84)

Consequently,    c∆ CΦ β 2 η 2 cΦ ¯ + ∆k ∥û+ − φh (ûk )∥H ≤ h ∆k + (1 − β)2 4 2     Ca CΦ β 2 η 2 k cΦ k+1 2 2 2 + h CΦ ∥fk ∥H + + β + β ∥v0 ∥Θ . 2 4 2 

(85)

Proof. For constant β, (77) gives 2 Bk+1 ≤ β k+1 ∥v0 ∥2Θ +

k ¯k X c∆ ∆ β k−m . (1 − β)h m=0

P Since km=0 β k−m ≤ (1 − β)−1 , this yields (84). Substituting the corresponding bounds for 2 Bk2 and Bk+1 into (79) gives (85).

19

5.2 Global error The following global bound shows that the semi-implicit Euler DFI trajectory is controlled by the accumulated local DFI defects along the computed path. In the uniform-step case hk = h, the estimate takes the form
¯ n−1 + h , ∥ûn − u(tn )∥H ≤ C ∆ where the constant C depends on the flow stability constant in (43), the curvature constants of Φ, the final time, and the DFI parameters. Thus the discretization error contribution is of the expected first order in time. Our bound is of the a posteriori type, it bounds the global error in terms of the observed defects ∆j along the computed trajectory; however, it does not by itself prove that these defects vanish under smaller time-step size h → 0. Proposition 6 (A posteriori global error bound). Under the hypotheses of Proposition 5, let u(t) denote the exact solution of (1) with u(0) = û0 . Then for every k ≥ 1 with tk ≤ T , " # k−1 X CΦ βj2 ηj2 2 2 cΦ 2 2 Ca 2 ℓ(tk −tj+1 ) 2 2 ∥ûk − u(tk )∥H ≤ e hj ∆j + CΦ hj ∥fj ∥H + hj Bj + hj Bj+1 + h , 4 2 2 j j=0

(86) where ℓ is the constant from the flow stability in Assumption 5.1(i). If, in addition, hk = h for all k, so that (82) holds, then there exists a constant CT > 0, depending only on T , ℓ, and universal numerical constants, such that for every k with tk ≤ T , " ∥ûk − u(tk )∥H ≤ CT

1+


c∆ 2 ¯ k−1 + ∆ c + C η Φ Φ (1 − β)2 # h

2 + Ca + CΦ F̄k−1

cΦ + CΦ η

2

∥v0 ∥2Θ

, (87)

where F̄k−1 :=

max ∥fj ∥H .

0≤j≤k−1

(88)

Proof. The dependence on the DFI memory parameter is explicit through the factor (1−β)−2 . Fix k ≥ 1 with tk ≤ T . The index k denotes the final time at which we want to estimate the global error ∥ûk − u(tk )∥H . We use j = 0, . . . , k − 1 to index the individual time steps whose local errors are propagated to the final time tk . Since u(0) = û0 , we have u(tk ) = φtk (û0 ). Now we introduce the intermediate propagated states j = 0, . . . , k . Wj := φtk −tj (ûj ) , At the end points, we have W0 = φtk (û0 ) = u(tk ) and Wk = φ0 (ûk ) = ûk . Therefore ûk − u(tk ) = Wk − W0 and telescoping (Lady Windermere fan) gives us ûk − u(tk ) = Wk − W0 =

k−1 X j=0

20

Wj+1 − Wj .

(89)

Now use the semi-group property φtk −tj = φtk −tj+1 ◦ φhj of the exact flow φ to obtain Wj+1 − Wj = φtk −tj+1 (ûj+1 ) − φtk −tj+1 (φhj (ûj )). Therefore, by the triangle inequality and flow stability, ∥ûk − u(tk )∥H ≤

k−1 X

eℓ(tk −tj+1 ) ∥ûj+1 − φhj (ûj )∥H ,

(90)

j=0

where we took the H-norm, applied the triangle inequality, and used the flow stability (43). Now notice that the terms ∥ûj+1 − φhj (ûj )∥H denote local errors in (90), so we can plug the bound (79) from Proposition 5 and Corollary 5.2 into (90) to obtain (86). It remains to prove the simplified estimate (87). Let ℓ+ := max{ℓ, 0}. Since 0 ≤ tk −tj+1 ≤ T , we have k−1 k−1 X X ℓ(tk −tj+1 ) ℓ+ T e h≤e h = eℓ+ T tk ≤ eℓ+ T T. (91) j=0

j=0

For 0 ≤ j ≤ k − 1, we have ¯ k−1 , ∆j ≤ ∆

¯j ≤ ∆ ¯ k−1 , ∆

∥fj ∥H ≤ F̄k−1 .

Also, since 0 < β < 1, β j ≤ 1, Hence

β j+1 ≤ 1,

β 2 ≤ 1.

CΦ β 2 η 2 cΦ + ≤ C CΦ η 2 + c Φ , 4 2

and

CΦ β 2 η 2 j cΦ j+1 β + β ≤ C CΦ η 2 + c Φ , 4 2 where C > 0 is a universal numerical constant. Applying these estimates to (86), and using (91), gives (87).

6 Numerical experiments We demonstrate the DFI scheme on examples with the Allen–Cahn and Fokker-Planck equations. 6.1 Allen–Cahn equation 6.1.1 Setup We consider the one-dimensional Allen–Cahn equation on the periodic domain [0, 2π), ∂t u(t, x) = ν∂xx u(t, x) + u(t, x) − u(t, x)3 ,

21

(t, x) ∈ [0, T ] × [0, 2π),

(92)

with periodic boundary conditions. Recall that (92) is an L2 -gradient flow because if we set Z 2π V (u) = 0

2 1 ν u(x)2 − 1 dx, |∂x u(x)|2 + 2 4

then ∂t u = −∇u V (u). We set the viscosity parameter to ν = 0.2 and consider the end time T = 15. The initial condition is u0 (x) =

π 2 1 2 2 tanh(2 sin x) − e−23.5(x− 2 ) + e−27(x−4.2) + e−38(x−5.4) . 3

6.1.2 Nonlinear parametrization with neural network We parametrize the approximate solution by ûk (x) = Φ(θk )(x), where Φ(θ) is a periodic fully connected feedforward neural network. The periodicity is imposed at the input level. Namely, instead of feeding x directly into the network, we first define the periodic feature vector  
2πx 16 ϕ(x) = ϕ1 (x), . . . , ϕ16 (x) ∈ R , ϕj (x) = cos + sj , L = 2π, L where the shifts sj are fixed. Since L = 2π, these features satisfy ϕj (x + L) = ϕj (x), and therefore every function obtained by composing a standard feedforward network with ϕ(x) is L-periodic. The network has five hidden layers, each of width 15, and uses the swish activation applied componentwise. The trainable parameters are the entries of the weight matrices, biases, output weights, and output bias. Their total number is p = 1231. Hence the parameter space is Θ = R1231 . 6.1.3 Setup of DFI scheme We use the L2 (0, 2π) inner product in space. In the numerical experiments this inner product is approximated on the uniform periodic grid x0 , . . . , xNx −1 with Nx = 450 points. We use the trapezoidal rule for approximating norms and inner products on it: for a function w ∈ L2 (0, 2π), we use Nx −1 2π X ∥w∥22,Nx = |w(xi )|2 . Nx i=0

In the empirical least-squares solves, both Jk and fk are evaluated on this grid using JAX’s automatic differentiation to compute all derivatives involved. The semi-implicit Euler DFI update is executed as described in Algorithm 1. In particular, the initialization parameter θ0 is computed by fitting the initial condition with the Adam optimizer run for 5 × 104 full batch iterations on the same grid, while the initial momentum is set to v0 = 0. We compare DFI to Tikhonov-regularized Dirac-Frenkel dynamics, which corresponds to the same update with βk = 0 for all k. For DFI and Tikhonov-regularized Dirac-Frenkel we keep β = βk and η = ηk fixed over all time steps k. Furthermore, we use a fixed time step h so that tk = kh,

k = 0, . . . , Nt ,

22

Nt =

T . h

10 2

0.8

time step size h 0.001 0.002 0.005 0.01

opt

Integrated relative error

10 1

1.0

Tikhonov-DF DFI

0.6 0.4 0.2

10 4 10 3 10 2 10 1 regularization parameter 2

10 3 10 2 regularization parameter 2

(a) integrated relative error

(b) optimal memory parameter βopt

Figure 1: DFI uses memory to compensate for loss of information in the instantaneous regularized least-squares solve. Plot (a) shows that as η 2 increases, Tikhonov-DF deteriorates because the velocity is computed from an increasingly regularized local problem alone, while DFI remains accurate by retaining more of the previous velocity. Plot (b) provides more evidence of this mechanism, showing that the error-minimizing βopt increases with η 2 , meaning that more previous-velocity information is used when the instantaneous least-squares signal is more strongly regularized

We also consider the left-sketched version of DFI and Tikhonov-regularized Dirac-Frenkel dynamics, where the correction is obtained by solving (s)

zk+1 = arg minp ∥Sk (Jk z − (fk − βJk vk ))∥22 + η 2 ∥z∥2Rp ,

(93)

(s)

(94)

z∈R

(s)

vk+1 = βvk + zk+1 ,

with a left-sketch (sub-sampling) matrix Sk ∈ Rs×Nx , with s ≤ Nx , where the quadrature weights from the empirical L2 norm are absorbed into Sk . The reference solution, denoted by uk , is computed to near machine precision with a Fourier spectral discretization in space and fourth-order Runge–Kutta in time with step size href = 2.5 · 10−4 . We regard this trajectory as the true solution for reporting error diagnostics. We report the pointwise relative error given by ek =

∥ûk − uk ∥2,Nx . ∥uk ∥2,Nx

(95)

and the integrated relative error, PNt E=

2 k=0 h ∥ûk − uk ∥2,Nx PNt 2 k=0 h ∥uk ∥2,Nx

23

!1/2 .

time step size h 0.001 0.002 0.005 0.01

10 7 10 9 10 4 10 3 10 2 10 1 regularization parameter 2

pointwise relative error

final energy

10 5

Tikhonov-DF DFI

Tikhonov-DF DFI

10 1 10 2 10 3 10 4 10 5

0

5

time

10

15

(b) point-wise error for η 2 = 2 × 10−2

(a) final energy

Figure 2: Plot (a) shows that DFI reaches a smaller final energy gap than Tikhonov-DF across regularization strengths, indicating that the inertial dynamics better follow the long-time energy decay in this example. Plot (b) shows that this improvement is not only a final-time effect. After the initial transient, DFI also gives smaller pointwise errors along the trajectory.

6.1.4 Benefit of inertia: robustness when the local Jacobian information is limited The DFI scheme remains accurate even when the instantaneous least-squares problem built from the instantaneous Jacobian Jk provides only limited reliable information about the next velocity. As we discussed earlier, this can happen because Jk is ill-conditioned or nearly rank-deficient, so some parameter directions are only weakly visible in the tangent vector Jk v. Additionally, the least-squares problem may be too strongly regularized locally, which suppresses the correction obtained from the instantaneous Jacobian. We also note a third option, which is that the residual may be evaluated only through a left sketch Sk Jk as in (93), so the update uses a compressed approximation of least-squares problem, which can be beneficial for speedups [4, 12]. In all three cases, the instantaneous local solve contains less usable information about the velocity that should be taken in parameter space. In TikhonovDF, weakly informed velocity directions are shrunk by the regularized instantaneous solve. By contrast, DFI combines the information available from the instantaneous Dirac–Frenkel residual with the history given by the inertia of the parameter velocity. This mechanism explains the behavior we see in Figure 1. The left plot shows that DFI is more robust than Tikhonov-DF as the regularization strength η 2 increases. For small and moderate values of η 2 , the two methods give comparable integrated errors. When η 2 becomes large, the correction obtained from the instantaneous least-squares problem is strongly penalized. In this regime, the Tikhonov-DF error grows rapidly because Tikhonov-DF recomputes a shrunk velocity from scratch at each step. DFI, on the other hand, remains accurate over a wider range of η 2 because the transported velocity βvk supplies information that is not obtained from the instantaneous regularized correction alone. The right plot in Figure 1 provides further evidence that this is the active mechanism. The optimal value of β to

24

integrated relative error

10 1

reg. param. 2 0.0007 0.007 0.07 101

Tikhonov-DF DFI 2 × 101

sketch size

3 × 101

4 × 101

Figure 3: As the sketch size s decreases, the local solve uses less information from the residual, so TikhonovDF deteriorates because it relies entirely on the instantaneous sketched leastsquares problem. DFI remains more accurate for small sketch sizes because the transported velocity supplements the information missing from the instantaneous sketched solve.

minimize the integrated error E increases with η 2 . Thus, when the correction from the instantaneous Jacobian is more strongly regularized, the best DFI trajectory compensates by retaining more of the previous velocity. This matches the direction-wise interpretation from Section 3.3. Directions that are well resolved by Jk are updated using the instantaneous Dirac–Frenkel residual, while weakly resolved or strongly damped directions can be carried forward by inertia instead of being instantaneously suppressed. The gradient-flow diagnostics in Figure 2 show the same behavior from the perspective of the energy. Across different regularization strengths, DFI reaches a smaller final energy gap than Tikhonov-DF, suggesting that the inertial dynamics help preserve the long-time relaxation structure of the Allen–Cahn flow. The pointwise-in-time comparison in Figure 2 shows that the improvement is not only a final-time effect: DFI reduces both the pointwise error and the least-squares defect after the initial transient. 6.1.5 DFI is more robust under sketching A smaller sketch size s reduces the cost of the least-squares solve, but it also means that the method uses less information from the residual. Since Tikhonov-DF relies entirely on this instantaneous sketched solve, its accuracy is more sensitive to a small sketch size. DFI is less sensitive because it also uses the velocity from the past. Thus, DFI can maintain accuracy even when the local least-squares problem is made less informative in order to improve stability or reduce cost. The results in Figure 3 demonstrate this effect. For small sketch sizes, DFI achieves lower integrated error than Tikhonov-DF. Thus the inertial memory is not merely a qualitative difference in the dynamics; it leads to a computational benefit, allowing one to use cheaper sketched Dirac–Frenkel updates with less loss of accuracy.

25

integrated relative error

Figure 4: As β approaches one, equivalently as 1 − β becomes small, the transported velocity dominates and the instantaneous least-squares correction has too little influence. The growth of the integrated error in this regime matches the (1 − β)−1 deterioration predicted by the error estimates in Section 5.

time step size h 0.001 0.002 0.005 0.01

10 1

10 2 10 3 1

10 1

6.1.6 Inertia helps, but excessive memory hurts. The results in Figure 4 show the expected tradeoff in the memory parameter β. Increasing β allows DFI to retain more information from the previous velocity, which is the main mechanism behind its robustness. However, taking β too close to one makes the method insufficiently responsive to the instantaneous Dirac–Frenkel residual. Then inaccurate or outdated velocity components can persist for too long, and the integrated error grows. This agrees the error analysis in Section 5. The bounds contain factors involving (1 − β)−1 , which deteriorate as β → 1. 6.2 High-dimensional Fokker–Planck equation 6.2.1 Setup We consider a Fokker–Planck equation for a probability density u(t, ·) on Rd , with d = 10, x ∈ Rd ,

∂t u(t, x) = −∇x · (u(t, x)b(t, x)) + D∆x u(t, x),

t ∈ [0, T ].

(96)

Here b(t, x) ∈ Rd is the drift field and D > 0 is the diffusion coefficient. The equation describes the evolution of the probability density of the stochastic differential equation √ dXi (t) = bi (t, X(t)) dt + 2D dWi (t), i = 1, . . . , d. (97) We set D = 10−2 and integrate until T = 2. The drift describes an interacting anharmonic trap [6]. Its components are d

bi (t, x) = (a(t) − xi )3 + α(x̄ − xi ),

x̄ =

1X xj , d

i = 1, . . . , d,

j=1

where a(t) = 1.25(sin(πt) + 1.5),

α = −0.5.

The cubic term gives a nonlinear restoring force centered at a(t), while the mean-field term couples each coordinate to the empirical mean x̄. The initial density is a Gaussian density.

26

On the computational box, we evaluate it using the wrapped displacement   L L δi (x, µ) := xi − µi + mod L − , L = 4, 2 2 and set 0

2 −d/2

u (x) = (2πσ )

  1 0 2 exp − 2 ∥δ(x, µ )∥2 , 2σ

σ 2 = 0.1.

The initial mean is placed along a line in the coordinate index, µ0i = 1.46 + 1.27

i−1 , d−1

i = 1, . . . , d.

Although the equation is posed on Rd , the computation is carried out on the box [0, 4)d . The box is chosen large enough to contain essentially all probability mass over the time interval considered. 6.2.2 Nonlinear parametrization with neural network Since the unknown is a probability density, we enforce positivity directly in the parametrization by setting Φ(θ)(x) = exp(−φ(θ)(x)). Here φ(θ) is a periodic fully connected feedforward neural network; see Section 6.1.2. This makes the neural representation periodic on the computational box [0, 4)d . The network φ(θ) has four hidden layers, each of width 64, and uses the swish activation applied componentwise. The trainable parameters are the entries of the weight matrices, biases, output weights, and output bias. In dimension d = 10, this architecture gives p = 18625 trainable parameters. 6.2.3 Setup of the DFI scheme The empirical least-squares problem is formed from importance samples on the computational box. At each time step, we draw d s Xk = {xℓ }N ℓ=1 ⊂ [0, 4) ,

Ns = 2000,

from an adaptive Gaussian mixture proposal. The proposal consists of 1500 samples from a Gaussian fitted to the current approximation ûk and 500 samples from the uniform distribution on [0, 4)d . Let qk denote this proposal density. The empirical norm used in the least-squares problem is the importance-sampling approximation of the L2 ([0, 4)d ) norm, N

∥w∥22,Xk =

s 1 X |w(xℓ )|2 . Ns qk (xℓ )

(98)

ℓ=1

The semi-implicit Euler DFI update is then conducted as in Algorithm 1 with the empirical norm (98) and the Gaussian fitted to the current solution ûk at each time step. We use a time step size of 10−3 and keep β and η fixed over all time steps. We initialize the parameter 27

variable to θ0 computed by fitting the initial condition with the Adam optimizer for 2 × 105 full batch iterations. To ensure our results are not polluted by the initial error, we do this employing a finer sample of size 2 × 105 , half being drawn from the Gaussian initial condition itself and the other half from the uniform measure, and using importance sampling to estimate the L2 -error. The momentum variable is initialized to v0 = 0. We compare DFI and Tikhonov-DF using moment diagnostics. The reference moments ref 5 µref k and Σk are computed from the 10 reference particles from the SDE (97) with the b k are computed by Euler–Maruyama method and time step 10−4 . The moments µ bk and Σ importance sampling from the same adaptive proposal used for evaluation. The weights are self-normalized using the ratio between the approximate density ûk and the proposal density qk . Because the computation is carried out on a periodic box, the mean is estimated coordinate-wise by circular averaging after mapping each coordinate from [0, 4) to the unit circle. The covariance is then computed from wrapped displacements around this circular mean. We report the pointwise relative errors eµ,k =

∥b µk − µref k ∥2 , ref ∥µk ∥2

eΣ,k =

b k ) − diag(Σref )∥2 ∥diag(Σ k . )∥ ∥diag(Σref 2 k

(99)

and the integrated moment errors, PNt Eµ =

2 µk − µref k=0 h ∥b k ∥2 PNt ref 2 k=0 h ∥µk ∥2

!1/2 ,

and PNt EΣ =

ref 2 b k=0 h ∥diag(Σk ) − diag(Σk )∥2 PNt ref 2 k=0 h ∥diag(Σk )∥2

!1/2 .

6.2.4 Results We use this high-dimensional Fokker–Planck example to demonstrate that the inertial memory in DFI improves robustness in a more challenging nonlinear parametrization. For both DFI and Tikhonov-DF, the regularization parameters are tuned for best performance via grid search. In DFI, this means tuning both the memory parameter β and the regularization strength η 2 . In Tikhonov-DF, we tune the corresponding Tikhonov regularization parameter. At each time step, the residual is evaluated at Ns = 2000 sample points from the adaptive mixture proposal. Thus the instantaneous Dirac–Frenkel signal is randomized. The results in Figure 5 show that DFI gives a more robust time evolution of the moment diagnostics in this regime. For both the mean and the diagonal covariance, Tikhonov-DF initially follows the reference dynamics, but the error grows sharply after a short time. DFI avoids this error growth and maintains substantially smaller moment errors over the time interval. Furthermore, Figure 6 confirms the improved stability of the DFI least-squares problem with respect to local information, compared with DF, consistently with the behavior observed in the Allen-Cahn example. In particular, the mean and covariance errors indicate that DFI is more robust under sketching with very small sketching dimensions. Again, this robustness can be attributed to the fact that the corresponding least-squares formulation depends 28

relative error

relative error

mean 10 1 10 2 0.0

0.5

covariance

102 101 100 10 1

1.0 1.5 2.0 time Tikhonov-DF

0.0

0.5

1.0 1.5 time

2.0

DFI

mean Tikhonov-DF DFI 10 2 101

sketch size

covariance relative error

relative error

Figure 5: Fokker-Planck in 10D: DFI is more robust in the sample-based high-dimensional Fokker–Planck experiment. The instantaneous least-squares problems are formed from Ns = 2000 only and the regularization parameter is set to η 2 = 106 . TikhonovDF develops large pointwise errors in the mean and diagonal covariance after a short time, whereas DFI maintains smaller errors over the time interval. Curves are averaged over five independent runs (same initialization).

102

4 × 10 1 3 × 10 1

reg. param. 2 10 1e+03 1e+04

2 × 10 1

10 1 101

sketch size

102

Figure 6: Fokker-Planck in 10D: DFI is more robust to smaller sketch size due to the memory term supplementing information from past solves. Results are averaged over five independent runs (same initialization).

less strongly on localized information because of the memory term. Overall, these results corroborate the trends previously observed for the Allen-Cahn problem.

7 Conclusions By adding inertia to the Dirac–Frenkel dynamics, the proposed DFI scheme allows parameter velocity information from the past trajectory to persist in directions that are weakly informed

29

by the instantaneous Jacobian. This mechanism is useful precisely in the regimes where the current Jacobian does not provide enough reliable information about all parameter velocity directions. Such situations occur when the parametrization is redundant, when the Jacobian is ill-conditioned or nearly rank deficient, when the least-squares problem is strongly regularized, or when the Jacobian is sketched. In these cases, Tikhonov-regularized Dirac–Frenkel dynamics suppress weakly informed directions instantaneously. DFI instead allows velocity in such directions to persist, while still using the current Dirac–Frenkel residual to correct directions that are well resolved by the Jacobian. We established well-posedness of the continuous DFI system and derived a posteriori error bounds that separate the usual Dirac–Frenkel projection defect from the additional relaxation defect introduced by inertia. We also introduced a semi-implicit Euler discretization, which requires one regularized least-squares solve per time step, with the previous velocity appearing as an anchor.

References [1] J. Aghili, J. Z. Atokple, M. Billaud-Friess, G. Garnier, O. Mula, and N. Tognon. A dynamical neural galerkin scheme for filtering problems. ESAIM: ProcS, 81:2–15, 2025. [2] W. Anderson and M. Farazmand. Evolution of nonlinear reduced-order solutions for PDEs with conserved quantities. SIAM J. Sci. Comput., 44(1):A176–A197, 2022. [3] H. Attouch, X. Goudou, and P. Redont. The heavy ball with friction method, i. the continuous dynamical system: Global exploration of the local minima of a real-valued function by asymptotic analysis of a dissipative dynamical system. Communications in Contemporary Mathematics, 02(01):1–34, 2000. [4] J. Berman and B. Peherstorfer. Randomized sparse Neural Galerkin schemes for solving evolution equations with deep networks. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 4097–4114, New Orleans, Louisiana, USA, 2023. Curran Associates, Inc. [5] D. Bon, B. Caris, and O. Mula. Stable nonlinear dynamical approximation with dynamical sampling. arXiv, 2505.11938, 2025. [6] J. Bruna, B. Peherstorfer, and E. Vanden-Eijnden. Neural Galerkin schemes with active learning for high-dimensional evolution equations. Journal of Computational Physics, 496:112588, Jan. 2024. [7] B. Carrel. Randomized methods for dynamical low-rank approximation. Journal of Computational Physics, 544:114421, 2026. [8] H. Chen, R. Wu, E. Grinspun, C. Zheng, and P. Y. Chen. Implicit neural spatial representations for time-dependent PDEs. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 5162–5177. PMLR, 23–29 Jul 2023. 30

[9] Z. Chen, J. Mccarran, E. Vizcaino, M. Soljacic, and D. Luo. TENG: Time-evolving natural gradient for solving PDEs with deep neural nets toward machine precision. In R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp, editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 7143–7162, Vienna, Austria, 21–27 Jul 2024. PMLR. [10] W. Dahmen, W. Li, Y. Teng, and Z. Wang. Expansive natural neural gradient flows for energy minimization, 2025. [11] P. A. M. Dirac. Note on exchange phenomena in the thomas atom. Mathematical Proceedings of the Cambridge Philosophical Society, 26(3):376–385, 1930. [12] Y. Dong, P. Schwerdtner, and B. Peherstorfer. Randomized time stepping of nonlinearly parametrized solutions of evolution problems. arXiv, 2512.19009, 2025. [13] Y. Du and T. A. Zaki. Evolutional deep neural network. Phys. Rev. E, 104:045303, Oct 2021. [14] M. Feischl, C. Lasser, C. Lubich, and J. Nick. Regularized dynamical parametric approximation. arXiv, 2403.19234, 2024. [15] M. A. Finzi, A. Potapczynski, M. Choptuik, and A. G. Wilson. A stable and scalable method for solving initial value PDEs with neural networks. In The Eleventh International Conference on Learning Representations, 2023. [16] J. Frenkel. Wave Mechanics, Advanced General Theory. Clarendon Press, Oxford, 1934. [17] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pizorn, H. Verschelde, and F. Verstraete. Time-dependent variational principle for quantum lattices. Physical Review Letters, 107:070601, 2011. [18] J. S. Hesthaven, B. Peherstorfer, and B. Unger. Nonlinear model reduction for transportdominated problems. Acta Numerica, 35:173–272, 2026. [19] Z. Hu, C. Liu, Y. Wang, and Z. Xu. Energetic variational neural network discretizations of gradient flows. SIAM Journal on Scientific Computing, 46(4):A2528–A2556, 2024. [20] M. Kast and J. S. Hesthaven. Positional embeddings for solving PDEs with evolutional deep neural networks. Journal of Computational Physics, 508:112986, 2024. [21] K. G. Kay. The matrix singularity problem in the time-dependent variational method. Chem. Phys., 137(1):165–175, 1989. [22] O. Koch and C. Lubich. Dynamical low-rank approximation. SIAM Journal on Matrix Analysis and Applications, 29(2):434–454, 2007. [23] S. Kvaal, C. Lasser, T. B. Pedersen, and L. Adamowicz. No need for a grid: Adaptive fullyflexible gaussians for the time-dependent Schrödinger equation. arXiv, 2207(00271):1–8, 2023. 31

[24] H. Y. Lam, G. Ceruti, and D. Kressner. Randomized low-rank Runge–Kutta methods. SIAM J. Matrix Anal. Appl., 46(2):1587–1615, 2025. [25] C. Lubich. On variational approximations in quantum molecular dynamics. Mathematics of Computation, 74(250):765–779, 2005. [26] C. Lubich. From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis. EMS Press, Berlin, Germany, 2008. [27] C. Lubich and J. Nick. Regularized dynamical parametric approximation of stiff evolution problems. arXiv, 2501.12118, 2025. [28] B. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964. [29] M. Raviola and B. Peherstorfer. A Dirac-Frenkel-Onsager principle: Instantaneous residual minimization with gauge momentum for nonlinear parametrizations of PDE solutions. In International Conference on Machine Learning (ICML), 2026. [30] G. Teschl. Ordinary Differential Equations and Dynamical Systems, volume 140 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 2012. [31] Y. Wang, J. Chen, C. Liu, and L. Kang. Particle-based energetic variational inference. Statistics and Computing, 31:1–17, 2021. [32] H. Zhang, Y. Chen, E. Vanden-Eijnden, and B. Peherstorfer. Sequential-in-time training of nonlinear parametrizations for solving time-dependent partial differential equations. SIAM Review, 2025. (accepted).

32