Source-linked AI summary
Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics
Conrad Ainslie, Pedram Hassanzadeh, Michael W. Mahoney, Ashesh Chattopadhyay
TL;DR
Long-rollout error in autoregressive neural emulators is poorly understood because training focuses on one-step accuracy without constraining propagated error. This paper analyzes the state Jacobian of the learned update map and finds that direct-step models are unstable, whereas integration-constrained models achieve near-neutral stability and linear error scaling.
Problem
Autoregressive emulators are trained for one-step error but evaluated over long rollouts, while propagated error and stability lack comparable a priori analysis.
Method
The framework linearizes rollout error dynamics using the learned update map’s state Jacobian and analyzes its spectral radius.
Results
Direct-step models admit eigenvalues above one, whereas integration-constrained models’ eigenvalues collapse near the unit circle, yielding neutral stability and linear error scaling.
Takeaways & Limitations
The largest Jacobian eigenvalue provides a rollout-free diagnostic of error amplification and supports a priori stability analysis across emulator designs.
Takeaways & Limitations
Near-neutral spectra do not guarantee long-term stability because accumulated one-step error and eigenvector alignment also affect realized amplification.
Abstract
from arXiv · showhide
Neural autoregressive models have rapidly emerged as powerful emulators of high-dimensional chaotic systems, yet their long-term instability and error growth remain poorly understood, leading to ad-hoc solutions. Here, we develop an eigenanalysis framework that reveals the dynamical origin of this error growth. By analyzing the Jacobian of the learned one-step update map with respect to the state, we show how inference-time error growth, and thus model stability, is governed by its spectral radius. Direct-step architectures (models that predict the next state from the previous one) generically admit unstable eigenvalues with magnitudes exceeding one, explaining the rapid divergence of these widely used models. In contrast, integration-constrained models (where the time derivative is estimated and integrated with a higher-order integrator) collapse their eigenspectrum onto the unit circle, yielding neutral stability and a universal linear error-scaling law. The largest eigenvalue of this Jacobian provides an architecture-agnostic, a priori diagnostic of short-term skill, long-term stability, and spectral bias, without requiring an expensive rollout. Leveraging this theory, we introduce a stability-promoting loss that explicitly regularizes Jacobian-driven error amplification, improving both forecast accuracy and dynamical robustness. Demonstrated across $29$ models spanning two architectures, several explicit and implicit integrators, and multiple loss functions on the Kuramoto-Sivashinsky system, our results establish a theoretical foundation for the design and evaluation of neural emulators of chaotic multi-scale dynamics. More broadly, our framework is a step toward the kind of a priori stability analysis that numerical analysis provides for discretizations of differential equations and that scientific machine learning currently lacks.
Significance Statement
Neural networks can emulate chaotic, multi-scale physical systems faster than traditional simulations, but may drift or become unstable during long rollouts; the paper introduces an eigenanalysis theory to address this challenge.
- Motivation: Neural networks can rival traditional simulations for chaotic, multi-scale flows including the atmosphere, climate, and ocean.The passage identifies these systems as application domains for neural-network emulators.
- Motivation: These neural emulators can run far faster than traditional numerical methods.
- Motivation: Long rollouts can cause neural emulators to drift or blow up, while remedies are typically found by trial and error through replacement simulations.The passage frames current remedy evaluation as requiring the simulation the network was designed to replace.
- Contribution: The paper develops an eigenanalysis theory intended to identify a governing quantity for these networks.The supplied passage introduces this theory but does not specify the quantity’s name or role.
1 Introduction
Neural autoregressive emulators perform well on high-dimensional chaotic systems but often lose long-term stability without a principled explanation. This work addresses that gap by using the state Jacobian’s eigenanalysis to diagnose rollout error amplification independently of architecture, integration scheme, and loss function.
- Motivation: Long-term emulation remains difficult because models can become unstable or excessively diffusive, while existing remedies rely on tuning or problem-specific fixes.Even carefully tuned models may still blow up, leaving unclear when and why stability is lost [Guan et al., 2025, Pedersen et al., 2025, Sambamurthy and Chattopadhyay, 20
- Motivation: Spectral bias may cause single-step models to miss high-wavenumber dynamics, whose errors grow through small–large-scale coupling and corrupt resolved dynamics.This explanation has been proposed for multi-scale systems such as turbulent flows with decaying energy spectra [Chattopadhyay et al., 2023b, Yu et al., 2024].
- Contribution: The paper analyzes the Jacobian of the learned one-step map with respect to the system state, using its spectrum to diagnose rollout error amplification without expensive long rollouts.This derivative governs discrete-map stability and differs from the parameter Jacobian commonly studied through backpropagation.
- Contribution: The Jacobian-based diagnostic is independent of architecture, integration scheme, and loss function, while combining its error-amplification measure with one-step error addresses short-term accuracy and spectral-bias growth.The framework is intended to answer whether integration constraints control rollout error growth and improve short-term accuracy using a pre-rollout quantity.
- Contribution: Direct-step models generically admit eigenvalues with magnitude well above one and lack a mechanism to control them, explaining their instability during autoregressive rollout.The theory is developed and evaluated on the Kuramoto–Sivashinsky system across two architectures, several explicit and implicit integration schemes, and losses with and without spectral regularization.
2 Linear stability analysis
The analysis separates one-step generalization error from Jacobian-mediated propagated error, identifying the largest Jacobian eigenvalue as the local stability control. Integration-constrained models keep eigenvalues near unity by design, whereas direct-step models lack this spectral constraint and are prone to instability.
- Error dynamics: The Jacobian separates generalization error from propagated error, with its eigenvalue spectrum determining whether rollout errors are amplified or damped.Training and validation directly control the one-step generalization error, while products of state Jacobians govern dynamical stability during rollout.
- Eigenvalue criterion: The largest eigenvalue magnitude sets local stability: |λmax| > 1 produces geometric error growth, whereas |λmax| = 1 gives neutral stability.This criterion follows from the eigendecomposition of the learned update-map Jacobian.
- Architecture dependence: Integration-constrained models have near-unit spectra because their Jacobians differ from the identity by a small timestep-weighted learned tangent operator, while direct-step models lack this control.The implicit formulation is J = (I − ∆t ∇H[N])^-1 and places an eigenvalue outside the unit circle only when the learned tangent operator has an expanding direction.
- Architecture dependence: For the Kuramoto–Sivashinsky resolution analyzed, the exact tangent operator would make explicit Euler unstable by three orders of magnitude, whereas learned operators satisfy ∥∆t ∇H[N]∥∼10^-3.The resulting near-unit spectrum reflects smooth learned high-wavenumber time derivatives rather than the algebra of I + ∆t ∇H alone.
- Error-growth law: When p∆t ≪1 and eigenvalues are tightly clustered near unity, rollout error grows linearly with step count, with measured slopes 2, 11, 101, and 1001 at lead times 1, 10, 100, and 1000.Models separate within these lines through small departures of |λmax| from one, with larger |λmax| producing faster error growth.
- Validity regime: The linear accumulation law is guaranteed only for p ≪ 1/∆t; longer rollouts exceed this regime, so continued slow accumulation there is empirical rather than predicted by the derivation.At ∆t = 10^-3, the guaranteed window is p ≪ 1000.
3 Results
Direct-step models rapidly become unstable, whereas integration-constrained models remain physically plausible with near-linear error growth and orders-of-magnitude smaller errors. Jacobian eigenvalues explain these outcomes and provide an a priori stability diagnostic, while spectral regularization and implicit integration further reduce error growth.
- 3 Results: Figure 2 shows direct-step models becoming unstable within a few steps, while integration-constrained models remain stable with physically plausible states and orders-of-magnitude smaller, step-proportional errors.Direct-step errors saturate within tens to hundreds of steps.
- 3 Results: The direct-step models have |λmax| > 1, whereas integration-constrained eigenvalues collapse near the unit circle, explaining rapid divergence versus near-neutral error accumulation.The measured Euler departures from unity, 1.75 × 10−4 and 2.01 × 10−4, agree within a factor of 1.5 with the linear-operator estimate 2.5 × 10−4.
- 3 Results: Figure 3 confirms the predicted linear law for integration-constrained models when p∆t ≪1, with deviations emerging near p∆t ≈0.1 and breakdown at p∆t ≈1.Direct-step MLP and FNO models depart after about 10 and 100 steps, respectively, then saturate at the climatological decorrelation level.
- 3 Results: Integration constraints make linear stability analysis usable by fixing |λmax| ≈1, whereas direct-step spectra span a wide range with |λmax| > 1 and lack a comparable rollout-wide linear law.Realized error amplification also depends on error direction through the Jacobian eigenbasis, not solely on eigenvalue magnitudes.
- 3 Results: Implicit integration constraints produce |λmax| closer to unity and lower error-growth rates than corresponding explicit constraints at every tested ∆t, with eigenvalue ordering matching measured RMSE ordering.This establishes the largest eigenvalue as a rollout-free diagnostic: larger |λmax| predicts faster error growth, poorer short-term accuracy, and weaker long-term stability.
- 3 Results: With identical PEC4 constraints, spectral regularization visibly reduces early-rollout error even though both FNO models have |λmax| = 1.0002 to legend precision.The regularizer suppresses high-wavenumber one-step error, and accumulated generalization error is the leading early-rollout contribution.
4 Conclusion
The paper addresses the disconnect between one-step training and long-horizon inference by using the Jacobian as a stability-analysis object. It distinguishes error injection from error amplification, explains architectural differences, and identifies important limits on the diagnostic.
- 4 Conclusion: The Jacobian framework separates one-step generalization error, which injects error, from the propagator J, which amplifies it during rollout.This explains why models with comparable one-step accuracy can diverge at different rates.
- 4 Conclusion: Direct-step and integration-constrained emulators differ fundamentally because only the latter builds the update through a derivative network and numerical integrator that constrain amplification.Direct-step Jacobians have no time-step dependence or spectral constraint, whereas integration-constrained operators exhibit near-unit spectral behavior.
- 4 Conclusion: A near-neutral spectrum supports slow error growth but cannot guarantee long-term stability when accumulated one-step error is large or eigenvector orientation increases realized amplification.The diagnostic describes a local tangent map, and eigenvalue magnitudes alone do not determine amplification.
- 4 Conclusion: The theory was developed and tested on the Kuramoto–Sivashinsky equation, with Jacobians evaluated at initial conditions and robustness reported across 100 random initial conditions.Its scaling law is restricted to small time steps, while the diagnostic does not track tangent-map evolution along the attractor.
- 4 Conclusion: The framework aims to provide scientific machine learning with architecture- and loss-agnostic stability analysis analogous to numerical analysis for differential-equation discretizations.It derives inference-relevant quantities from the trained network’s dynamical operator rather than inspecting internal representations.
6 Methods and Systems
The study evaluates autoregressive neural emulators of the Kuramoto–Sivashinsky system using direct-step and integration-constrained architectures, explicit and implicit integrators, and spectral or stability-promoting losses. Models are trained on DNS data and assessed through long autoregressive emulations, with Jacobian-informed regularization targeting next-step error growth.
- Experimental system and data: Experiments use Kuramoto–Sivashinsky DNS with 1024 spatial points, Δt_DNS = 10^-3, 150000 training samples, and 100000-step evaluation rollouts from a new initial condition.Depending on the model, emulation becomes unstable, unphysical, or remains stable and physically consistent.
- Autoregressive architectures: The model suite compares direct-step networks with integration-constrained models that learn the time derivative and apply differentiable higher-order numerical integration.The integration-based hard constraint parameterizes F with a neural network while H performs the numerical integration.
- Time integration: Integration constraints include explicit Euler, RK4, and PEC4 schemes, while implicit integration uses a fixed-point implicit layer and time steps ranging from Δt_DNS to 100Δt_DNS.For explicit schemes, training and inference use Δt = Δt_DNS except for the large-step implicit comparison.
- Neural architectures: The neural operators are implemented with a six-layer, 2000-neuron ReLU MLP or a six-block FNO retaining 512 modes with width 32.Hyperparameter trials selected these configurations by optimizing short-term accuracy.
- Loss functions: To address spectral bias, the training framework supplements standard RMSE with a spectral regularizer whose Lagrange multiplier is set to γ = 0.10.The regularizer targets high-wavenumber dynamics that conventional L2 training may fail to learn, contributing to instability and physical inconsistency.
- Stability-promoting loss: The proposed stability-promoting loss minimizes the linearized next-step error ||e(t + Δt)||2 by projecting the predicted state through the model Jacobian.A more aggressive largest-eigenvector projection is discussed, but the proposed loss uses the Jacobian-based next-step error approximation.
Supporting Information Text
The supporting information provides four result sets covering the full 29-model spectral analysis, implicit integration constraints, explicit-scheme comparisons across time steps, and Fourier-space spectral-bias analysis.
- Section 1: The supporting results comprise a full 29-model suite with measured Jacobian spectra and a time-normalized growth rate for comparing models trained at different time steps.
- Section 2: They derive the Jacobian of an implicit integration-based update map and compare implicit and explicit schemes across a range of ∆t.
- Section 3: A separate analysis examines spectral bias in Fourier space.
1 The full model suite
The full suite contains 29 models, and growth-rate normalization separates direct-step from integration-constrained models by more than two orders of magnitude. Because σ derives from per-step |λmax|, comparisons across time steps require this normalization and should respect its precision limits.
- Model suite: 29 models comprise 27 integration-constrained and 2 direct-step models, with 24 explicit integration-constrained models summarized jointly by RMSE growth and Jacobian eigenspectra.All explicit models use ∆t = ∆tDNS = 10^-3, except two large-∆t implicit models.
- Normalization: σ = ln |λmax|/∆t converts the dominant per-step multiplier into a growth rate comparable across time steps and carrying Lyapunov-exponent units.The relation |λmax| = exp(σ ∆t) interprets the learned update map as a time-∆t flow map.
- Growth-rate separation: 140–180 times larger growth rates distinguish the two direct-step models from the integration-constrained models: σ = 48.1 and 60.2 versus σ ≤0.34.The majority of integration-constrained models have σ between 0.01 and 0.25; this growth-rate separation, rather than per-step eigenvalue differences, explains the rollout contrast.
- Precision and scope: At ∆t = 10^-3, Table S1 values resolve σ to only one or two significant figures, so small within-class differences should not be over-interpreted.Equation (S1) describes the modulus of the dominant multiplier, not the full complex spectrum.
2 Implicit integration-based hard constraints
Implicit integration-based hard constraints implement higher-order implicit time steppers through fixed-point implicit layers, producing Jacobians with larger stability regions than explicit constraints. Across tested step sizes, they reduce error growth while preserving physically consistent long-term emulation, with validity limits at large time steps and lead times.
- Implementation: Implicit schemes are implemented by solving the higher-order time-stepper equation with a fixed-point iteration inside an implicit layer at each training epoch.The method initializes y* = u(x,t) and iterates toward u(x,t+Δt), with experiments spanning ΔtDNS through 100 ΔtDNS.
- Validity and limitations: The implicit Jacobian matches the explicit integration-constrained Jacobian to first order in Δt, with differences beginning at O(Δt^2), which matters in large-Δt experiments.At Δt = 10^-1, neglected O(Δt^2) terms make the linear theory’s quantitative predictions less precise.
- Spectral mechanism: Implicit Jacobians place eigenvalues outside the unit circle only for expanding learned directions, whereas explicit Jacobians can expand mildly contracting directions, yielding smaller |λmax|.Thus, at fixed Δt, the implicit spectrum lies closer to and inside the unit disk, reflecting a larger stability region.
- Empirical comparison: Implicit constraints yield lower error-growth rates than explicit constraints at every tested Δt: σ = 0.100, 0.004, and 0.187 versus 0.200, 0.155, and 0.243.The comparison uses growth rates per unit time because |λmax| is a per-step multiplier and differs across step sizes.
- Empirical comparison: An implicit Euler FNO with |λmax| = 1.0189 remains stable and physically consistent, while its σ = 0.187 is roughly 260–320 times below two direct-step models.The per-unit-time comparison gives 0.187 for the implicit model versus 48.1 and 60.2 for the direct-step models.
- Validity and limitations: At the longest lead time, models 26 and 27 reach pΔt = 50 and 100, exceeding the pΔt ≪ 1 validity condition by two orders of magnitude.Their departure from the reference scaling line is therefore predicted by the stated validity window rather than evidence against the scaling law.
3 Spectral bias and the Fourier spectrum of the emulation
The section connects Jacobian eigenvalues to Fourier-spectrum fidelity and spectral bias in neural autoregressive models. It finds that |λmax| reflects implicit diffusivity and complements one-step error, while its association with high-wavenumber fidelity is established only empirically for one matched model pair.
- 3 Spectral bias and the Fourier spectrum of the emulation: The Fourier-space analysis links model-Jacobian eigenvalues to the fidelity of predicted spectra and their time derivatives.The section specifically examines whether |λmax| explains differing spectral-bias values across neural autoregressive architectures.
- 3 Spectral bias and the Fourier spectrum of the emulation: A spectral regularizer can reduce spectral bias and improve stability, but does not eliminate bias in the state or its derivative for the decaying KS spectrum.This limitation is shown in Fig. S3(b,c) and follows earlier work [Chattopadhyay et al., 2023b, Guan et al., 2025].
- 3 Spectral bias and the Fourier spectrum of the emulation: |λmax| reflects a model’s implicit diffusivity: smaller values indicate greater diffusion and often correlate with reduced spectral bias.Neural operators have also been shown to better capture spatial Fourier spectra of chaotic systems [Azizzadenesheli et al., 2024, Oommen et al., 2025].
- 3 Spectral bias and the Fourier spectrum of the emulation: |λmax| complements one-step error as an a priori indicator by setting amplification scale, while one-step error ε sets error injection.For the matched model pair, it also tracks implicit diffusivity and high-wavenumber spectral fidelity.
- 3 Spectral bias and the Fourier spectrum of the emulation: The correspondence between |λmax|, implicit diffusivity, and high-wavenumber fidelity is an empirical association established for a single matched model pair, not a general law.The main text’s Section 3.3 provides the more general statement.
4 Near-normality of the integration-constrained Jacobian and validity of the eigenvalue approximation
For integration-constrained models, the Jacobian is near-normal because its departure from normality is O(∆t^2), making eigenvalue-based error-growth estimates valid at leading order. Direction-dependent effects remain O(∆t) corrections, while the guarantee weakens at larger time steps and does not extend to direct-step models.
- Near-normality: The constrained Jacobian is normal to first order in ∆t, because its self-commutator is O(∆t^2), regardless of the network Jacobian’s normality.Writing J = I + ∆t A causes the identity to commute and the terms linear in ∆t to cancel.
- Eigenvalue approximation: At leading order, replacing ||J e|| with |λmax| ||e|| is valid because all singular values and |λmax| equal 1 + O(∆t).For any error-vector orientation, the induced amplification differs from the eigenvalue estimate only by corrections of the same O(∆t) size already neglected.
- Finite-step corrections: At ∆t = ∆tDNS = 10−3, the approximation residue is small, consistent with the fourth-decimal spread of measured spectra in Fig. 2(d) and Table S1.The remaining directional effects arise from error-eigenvector projections and accumulated one-step error within the O(∆t) band around unity.
- Limitation: Direct-step models lack ∆t scaling in their Jacobians, so their departure from normality is uncontrolled and eigenvalue-based amplification estimates are not guaranteed.Clustered eigenvalues alone cannot constrain transient amplification for strongly non-normal matrices.
- Validity range: The near-normality guarantee weakens at ∆t = 10−1, where the O(∆t^2) departure is four orders of magnitude larger than at ∆t = ∆tDNS.The same leading-order expansion, J = I + ∆t A + O(∆t^2), applies to implicit models.