Artifacts of Numerical Integration in Learning Dynamical Systems

1. Problem Statement

The paper addresses a critical issue in dynamical system learning: the process of inferring the underlying vector field $f(y)$ of an autonomous system $\dot{y} = f(y)$ from discrete trajectory observations $\{(t_n, y(t_n))\}$.

In practice, this learning problem is formulated as an optimization task where a candidate function $\hat{f}$ is selected to minimize the mismatch between observed data and the numerical solution of the candidate system. Crucially, the numerical solution is generated using a numerical integrator (e.g., Runge-Kutta, Euler methods) with a specific step size $h$.

The core problem identified is that the choice of numerical integrator introduces systematic artifacts that distort the learned dynamics. Even if the learned model fits the training data perfectly, the resulting continuous-time system $\dot{y} = \hat{f}(y)$ may exhibit qualitatively incorrect behaviors, such as:

• Anti-damping: Identifying a dissipative (decaying) system as expansive (growing).
• Reversed oscillation: Incorrectly predicting the direction of rotation in oscillatory systems.
• Instability: Learning a system that is unstable for initial conditions outside the training distribution.

The authors challenge the common assumption that reducing the step size $h$ or increasing the order of an explicit integrator will automatically resolve these issues.

2. Methodology

The authors employ a rigorous mathematical analysis combining numerical analysis (stability theory) and inverse problem theory.

• Scalar Reduction: To isolate the effects of the integrator, the authors reduce the general nonlinear problem to a scalar linear model $\dot{z} = \lambda z$, where $\lambda$ is an eigenvalue of the true system. They analyze the learning of the parameter $\hat{\lambda}$ (the learned eigenvalue) using least-squares fitting.
• Stability Region Analysis: The core of the methodology relies on the stability function $p(\xi)$ of the numerical integrator (where $\xi = \lambda h$). The learned parameter $\hat{\lambda}_h$ is constrained by the equation $p(\hat{\lambda}_h h) = e^{\lambda h}$.
• Geometric Characterization: Instead of relying on asymptotic error expansions (which assume $h \to 0$), the authors analyze the finite-$h$ geometry of the integrator's stability region in the complex plane. They demonstrate that the learned $\hat{\lambda}_h$ must lie within the closure of the integrator's stability region.
• Perturbation Analysis: The paper extends the analysis to noisy data, quantifying how observational noise and the number of trajectories ($M$) affect the learned parameter.
• Nonlinear Extension: The scalar results are extended to nonlinear systems (e.g., Lotka-Volterra, pendulums) via local linearization (Jacobian analysis), showing that the artifacts observed in the scalar case propagate to the local behavior of the learned vector field.

3. Key Contributions

A. The "Stability Region" Artifact

The paper establishes that the geometry of the integrator's stability region dictates the qualitative nature of the learned dynamics, not just the accuracy of the fit.

• The Artifact: If an integrator's stability region extends into the right half of the complex plane (where $\text{Re}(z) > 0$), it is possible to learn a parameter $\hat{\lambda}$ with a positive real part even when the true system is dissipative ($\text{Re}(\lambda) < 0$) or conservative ($\text{Re}(\lambda) = 0$).
• Counter-Intuitive Finding: Higher-order explicit methods (like RK4) have stability regions that extend further into the right half-plane than lower-order methods (like Euler). Therefore, increasing the order or decreasing the step size does not necessarily fix the artifact; it can sometimes worsen it.

B. Resolution of Non-Uniqueness in Linear Multistep Methods

For Linear Multistep Methods (LMMs), the paper resolves a known non-uniqueness issue in learning. It shows that while the learned parameter might not lie in the absolute stability region, the optimization process effectively selects a specific root of the characteristic polynomial. However, this selection can lead to "spurious roots" that are excited by noise, causing the learned system to diverge for initial conditions not present in the training data.

C. Noise Propagation and Data Aggregation

The authors provide a perturbation analysis showing that the error in the learned parameter scales with the noise level $\sigma$. Crucially, they prove that aggregating multiple independent trajectories ($M$) reduces the variance of the learned parameter by a factor of $M$, effectively suppressing the excitation of spurious unstable modes.

D. The Prescriptive Recommendation: Implicit Midpoint Rule

The paper identifies the Implicit Midpoint Rule as the unique principled choice for learning autonomous systems when the structure (conservative vs. dissipative) is unknown.

• Its stability region coincides exactly with the left half-plane (including the imaginary axis).
• It preserves conservative modes exactly ($\text{Re}(\hat{\lambda}) = 0$ if $\text{Re}(\lambda) = 0$).
• It correctly captures dissipative behavior ($\text{Re}(\hat{\lambda}) < 0$ if $\text{Re}(\lambda) < 0$).
• Unlike the Implicit Euler method (which biases towards expansion) or Explicit Euler (which biases towards over-damping), the Implicit Midpoint Rule avoids sign errors in the real part of the eigenvalues.

4. Key Results and Numerical Experiments

The theoretical findings are validated through several numerical experiments:

1. Damped Nonlinear Pendulum:

   • Observation: Using RK4 (a high-order explicit method), the learned system correctly fits the sampled data points but, when simulated with high resolution, exhibits expansive behavior (growing amplitude) instead of damping.
   • Cause: The RK4 stability region allows $\hat{\lambda}$ to have a positive real part.
   • Solution: The Implicit Midpoint rule correctly recovers the damping.

2. Convection-Diffusion Equation:

   • Observation: Using the Implicit Euler method, the learned diffusion coefficient $\hat{\epsilon}$ can become negative, implying anti-diffusion (instability), even though the true system is diffusive.
   • Cause: The stability region of Implicit Euler includes parts of the right half-plane, allowing $\text{Re}(\hat{\lambda}) > 0$.

3. Lotka-Volterra System (Neural ODEs):

   • Observation: Neural ODEs trained with Implicit Euler produce trajectories that spiral outward (expansive), failing to capture the closed orbits of the true conservative system.
   • Solution: Neural ODEs trained with the Implicit Midpoint rule successfully preserve the closed orbit structure and limit cycles.

4. Relaxation Oscillator:

   • Observation: Explicit Euler learning results in excessive damping, causing trajectories to fail to converge to the true limit cycle.
   • Mechanism: The stability region geometry forces the learned eigenvalues to have a larger negative real part than the truth.

5. Significance and Implications

• Redefining Stability in Learning: The paper argues that traditional numerical stability (ensuring the solver doesn't blow up for a fixed equation) is insufficient for learning dynamical systems. The "stability" of the learning process depends on the geometric alignment of the integrator's stability region with the physical properties of the system being learned.
• Impact on Deep Learning: For Neural ODEs and other data-driven dynamical system models, the choice of the internal integrator is not just a computational detail but a structural prior. Using standard high-order explicit solvers (like RK4) can inadvertently introduce physics-violating artifacts (e.g., energy creation in dissipative systems).
• Practical Guidance: The authors provide a clear, actionable recommendation: Use the Implicit Midpoint Rule as the default integrator for learning autonomous dynamical systems from data. It offers a robust balance, preserving both conservative and dissipative structures without requiring prior knowledge of the system's specific nature.

In conclusion, this work reveals that numerical artifacts in dynamical system learning are not merely convergence errors but fundamental structural distortions caused by the geometry of numerical stability regions. Correcting these artifacts requires a principled choice of integrator rather than simply refining the step size.