Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics

Conrad Ainslie; Pedram Hassanzadeh; Michael W. Mahoney; Ashesh Chattopadhyay

cs.AIAug 17, 2026

Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics

Authors: Conrad Ainslie, Pedram Hassanzadeh, Michael W. Mahoney, Ashesh Chattopadhyay

Organizations: Department of Applied Mathematics, University of California, Santa Cruz, Santa Cruz, 95064, CA · Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637 · Committee on Computational and Applied Mathematics, University of Chicago, Chicago, IL 60637 · International Computer Science Institute, Berkeley, CA 94720 · Lawrence Berkeley National Laboratory, Berkeley, CA 94720 · Department of Statistics, University of California, Berkeley, CA 94720

Abstract

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.

Figures & tables

Figure 1: Schematic representation of the main contributions of this paper. An eigenvalue decomposition of the Jacobian of the learned one-step update map, taken with respect to the system state and not with respect to the network parameters, yields a linear theory for the role the Jacobian plays in short-term error growth, the rate of error accumulation over rollout, and the spectral fidelity of the emulation. Direct-step models, which predict the next state in one shot, admit eigenvalues outside the unit circle and diverge during rollout, whereas integration-constrained models, in which the network estimates a time derivative that is advanced by a numerical integrator, collapse their eigenspectrum onto the unit circle and remain neutrally stable, so that a perturbation is neither amplified nor damped from one step to the next (Section 2.3 ). The resulting stability-promoting loss regularizes the Jacobian-driven amplification of error, reducing the rollout error-growth rate: for some models this appears directly as eigenvalues pulled toward the unit circle (illustrated), while more generally it suppresses the realized, direction-dependent amplification of the error without necessarily changing $|\lambda_{\max}|$ (Section 3.4 ).

Panel of Figure 2 — Algorithm 1 Implicit time integration scheme

Figure 2: Short-term error, rollout error accumulation, and eigenanalysis of the models. Throughout, direct-step models are shown in vermilion and integration-constrained models in blue, with MLP as solid lines and FNO as dotted lines. (a) RMSE over the first $100$ steps. The two direct-step models (Direct MLP, Direct FNO) grow rapidly and saturate near the climatological decorrelation level, where the RMSE is $\mathcal{O}(1)$ in these normalized units. The two integration-constrained Euler models (Euler MLP, Euler FNO), whose eigenvalue structure is representative of all $27$ integration-constrained models examined here, remain stable over this window, with error orders of magnitude smaller and indistinguishable from zero on the linear axis. (b) The same RMSE on log–log axes over the first $10^{4}$ steps. The direct-step models saturate within a few tens of steps (MLP) to a few hundred (FNO), whereas the Euler models grow slowly and in proportion to the step count, remaining orders of magnitude more accurate over most of the window and reaching the decorrelation level only at its right edge ; over the full emulation several MLP-based models eventually drift as their accumulated generalization error grows. (c) Eigenvalues of the model Jacobian $\mathbf{J}$ at the initial condition $t_{0}$ on the Argand plane for all four models, with $|\lambda_{\max}|$ in the legend. The Direct MLP ( $\circ$ ) and Direct FNO ( $\square$ ) admit eigenvalues with $|\lambda_{\max}|>1$ ( $1.0493$ and $1.0620$ ) spread well inside and beyond the unit circle, whereas the Euler MLP ( $\star$ ) and Euler FNO ( $\triangle$ ) eigenvalues collapse onto the unit circle (dotted) and cluster tightly at $\lambda=1$ ( $|\lambda_{\max}|=1.0002$ ), appearing as a single dense cluster. (d) Zoomed-in view of the integration-constrained eigenvalues from (c), unit circle again dotted. Each spectrum spreads only in the third decimal place, and the two models separate in the fourth , with $|\lambda_{\max}|=1.000175$ for the Euler MLP and $1.000201$ for the Euler FNO. This neutral spectrum, consistent across architecture, explains the slow accumulation in proportion to the step count in (a) and (b) through the Jacobian structure of Eq. ( 10 ). SI, Fig. S1 shows the equivalent curves and spectra for the remaining explicit models, and Table S1 reports $|\lambda_{\max}|$ for every model.

Figure 3: A single linear scaling law governs error growth across the full model suite. Norm of the accumulated error $\|\mathbf{e}(t+k\Delta t)\|_{2}$ plotted against $|\lambda_{\max}|\,\|\mathbf{e}(t)\|_{2}$ for all $29$ models, where $\lambda_{\max}$ is the largest eigenvalue of the model Jacobian $\mathbf{J}$ evaluated at the initial condition and $\|\cdot\|_{2}$ is the Euclidean ( $L_{2}$ ) norm of the error vector. Here $t=t_{1}$ , so $\mathbf{e}(t)=\mathbf{e}(t_{1})=\boldsymbol{\epsilon}(t_{0})$ is the first-step (generalization) error and the lead time $k$ is measured from the first step; the reference slope at lead time $k$ is therefore the step count $p=k+1$ of Eq. ( 16 ). Panels (a)–(d) show the successive lead times $k=1,10,100,1000$ , one per panel, on log–log axes. Marker color denotes model class (direct-step in vermilion, integration-constrained in blue) and fill denotes architecture (MLP filled, FNO open), while marker shape indicates the integration scheme (Euler, RK4, PEC4, implicit Euler). The gray reference line in each panel is the linear scaling law of Eq. ( 16 ) evaluated at that lead time; on log–log axes it is a unit-slope line whose vertical offset grows with the step count, from slope $2$ in (a) to slopes $11$ , $101$ , and $1001$ in (b)–(d). The panels correspond to $p\,\Delta t=0.002$ , $0.011$ , $0.101$ , and $1.001$ , so they span the validity window of Eq. ( 16 ) from well inside it to its boundary. The integration-constrained models ( $|\lambda_{\max}|\approx 1$ ) lie on the predicted line in (a) and (b), scatter modestly about it in (c), and depart from it in both directions in (d), where $p\,\Delta t\approx 1$ and the derivation no longer applies. The two direct-step models sit at much larger $|\lambda_{\max}|\,\|\mathbf{e}(t)\|_{2}$ and rise above the reference line at the lead time of (b), reflecting the faster-than-linear growth driven by their $|\lambda_{\max}|>1$ ; by (c) and (d) they have reached the climatological decorrelation level and stopped growing, so they fall below a reference line that continues to grow in proportion to $p$ . This faster-than-linear growth arises from the $\mathcal{O}(\|\mathbf{e}(t)\|_{2}^{2})$ higher-order terms in the Taylor expansion of $\tilde{\mathbf{G}}$ , which are neglected in the linear analysis but become non-negligible once the accumulated error grows large. Each model is identified in Table S1. Two of the implicit-Euler models are trained at $\Delta t\gg\Delta t_{\text{DNS}}$ and therefore exceed the $p\,\Delta t\ll 1$ window by a larger margin at the longest lead time (SI, Section 2.4).

Figure 4: Error growth is set by the accumulated generalization error and the error–eigenvector projection, not by $|\lambda_{\max}|$ alone, and is predicted by the linearized error law. (a, b) Two FNO models constrained with the same PEC4 integrator, trained without (FNO, blue) and with (Spectral FNO, green) the spectral regularizer $\mu(\theta)$ of Eq. ( 36 ). (a) RMSE over the autoregressive rollout on log–log axes. For each model the solid line is the total error $\mathbf{e}(u(t))$ and the black dashed line is the cumulative generalization error $\sum_{i=1}^{t}\boldsymbol{\epsilon}(u(i))$ , the leading term of Eq. ( 14 ), with the dash length distinguishing the two models. The spectrally regularized model has the smaller error over the early rollout. (b) Eigenvalues of the model Jacobian near $\lambda=1$ for the FNO ( $\circ$ ) and Spectral FNO ( $\star$ ), with the unit circle dotted; to the precision of the legend the two models share $|\lambda_{\max}|=1.0002$ . Each marker is shaded by the cosine similarity $|\cos(\mathbf{e},\mathbf{v}_{i})|$ between its eigenvector $\mathbf{v}_{i}$ and the error vector $\mathbf{e}(t)$ (grayscale bar). (c, d) The linearized error law of Eq. ( 20 ), evaluated with $a=|\lambda_{\max}|$ from the Jacobian at the initial condition, tested after $k=10$ rollout steps ( $t_{k}=t_{0}+k\Delta t$ ) for every explicit integration-constrained model (marker shape: integration scheme, Euler / RK4 / PEC4; fill: MLP filled, FNO open); each point is one model and the gray line is $y=x$ . (c) Total error, with $\|\epsilon\|$ the measured one-step error. (d) The same law restricted to the high-wavenumber band ( $\omega\geq\omega_{c}$ , with $\omega_{c}=100$ ), with the injection term replaced by $\beta$ , the high- $\omega$ content of the one-step error. See Section 3.3 for the analysis. The two models in (a, b) are models 11 and 12 of Table S1, and the models in (c, d) are models 1 to 24.

Figure 5: The stability-promoting loss lowers error growth. FNO-based models constrained with the Euler, PEC4, and RK4 integrators, each trained with the baseline RMSE loss and compared against its counterpart trained with the proposed Jacobian-based stability-promoting loss of Eq. ( 38 ). (a) The linear scaling law of Fig. 3 , $\|\mathbf{e}(t+\Delta t)\|$ versus $|\lambda_{\max}|\,\|\mathbf{e}(t)\|$ , with the RMSE baselines shown as circles and the stability-promoting loss as stars. All points lie on the predicted unit-slope line of offset $2$ (gray), confirming that the new loss preserves the linear scaling law, and for every integrator the stability-promoting loss yields a smaller $|\lambda_{\max}|\,\|\mathbf{e}(t)\|$ and correspondingly smaller $\|\mathbf{e}(t+\Delta t)\|$ than the baseline. (b) RMSE over the autoregressive rollout for the same models, with the RMSE baselines shown as solid lines and the stability-promoting loss as dotted lines. The legend reports $\alpha$ , the slope of a linear fit of $\log_{10}(\mathrm{RMSE})$ against $\log_{10}$ of the timestep over the first $100$ rollout steps (the exponent of $\mathrm{RMSE}\propto t^{\alpha}$ ). For every integrator, the stability-promoting loss produces a smaller $\alpha$ , and hence lower error growth, than the corresponding RMSE baseline, because it constrains the Jacobian-driven amplification of error, which the RMSE objective alone leaves uncontrolled. The six models shown are models 7, 9, 11, 19, 21, and 23 of Table S1.

#	Arch.	Scheme	Loss	$\Delta t$	$\|\lambda_{\max}\|(t_{0})$	$\sigma$	RMSE ( $k=10$ )	Spec. bias
1	MLP	Euler	RMSE	$10^{-3}$	$1.000\,18$	$0.180$	$1.289\text{\times}{10}^{-03}$	$12.46$
2	MLP	Euler	Spectral	$10^{-3}$	$1.000\,02$	$0.020$	$3.850\text{\times}{10}^{-04}$	$10.46$
3	MLP	RK4	RMSE	$10^{-3}$	$1.000\,21$	$0.210$	$1.310\text{\times}{10}^{-03}$	$12.47$
4	MLP	RK4	Spectral	$10^{-3}$	$1.000\,07$	$0.070$	$3.970\text{\times}{10}^{-04}$	$10.44$
5	MLP	PEC4	RMSE	$10^{-3}$	$1.000\,18$	$0.180$	$1.273\text{\times}{10}^{-03}$	$12.49$
6	MLP	PEC4	Spectral	$10^{-3}$	$1.000\,30$	$0.300$	$3.577\text{\times}{10}^{-04}$	$10.35$

Table S1: The full model suite. All 29 models evaluated in this work, comprising 27 integration-constrained models and 2 direct-step models. Every explicit model is trained and run at

\Delta t=\Delta t_{\mathrm{DNS}}=10^{-3}

; the two large-

\Delta t

implicit models (26, 27) are the exception.

|\lambda_{\max}|(t_{0})

is the largest-magnitude eigenvalue of the model Jacobian at the initial condition.

\sigma=\ln|\lambda_{\max}|/\Delta t

is the corresponding growth rate per unit time of Eq. ( S1 ), which is comparable across

\Delta t

; its precision is limited by the five decimal places reported for

|\lambda_{\max}|

. RMSE is the short-term error at

k=10

steps. The spectral-bias metric is

\langle\log_{10}(|\hat{u}_{\mathrm{pred}}|/|\hat{u}_{\mathrm{true}}|)\rangle_{\omega\geq 100}

at step 100, with positive values indicating excess high-wavenumber energy relative to truth.

Figure S1: Short-term error and eigenspectra of the full explicit integration-constrained suite. All 24 explicit integration-constrained models of Table S1 , spanning the Euler, RK4, and PEC4 schemes, the MLP and FNO architectures, and the RMSE, spectral, Jacobian, and spectral+Jacobian losses. Every model is trained and emulated at $\Delta t=\Delta t_{\text{DNS}}=10^{-3}$ , with MLP models drawn as solid lines and FNO models as dotted lines, and color and marker distinguishing individual models as indicated in the legend. (a) RMSE over the first 100 rollout steps on log–log axes. Every model exhibits the slow, near-linear error accumulation predicted by Eq. ( 16 ) of the main text; the curves separate in vertical offset, set by the one-step generalization error $\boldsymbol{\epsilon}$ , rather than in slope, consistent with the near-neutral spectra shared across the suite. The window shown satisfies $p\,\Delta t\leq 0.1$ , so it lies inside the validity range of that scaling law. (b) Eigenvalues of the model Jacobian $\mathbf{J}$ , evaluated at the initial condition $t_{0}$ , on the Argand plane, with the unit circle dotted. For every model the eigenvalues collapse onto the unit circle and cluster tightly at $\lambda=1$ , with the spread of each spectrum confined to the third decimal place and the models separated in the fourth (Table S1 reports $|\lambda_{\max}|(t_{0})$ for each). This confirms that the qualitative spectrum is insensitive to the architecture, the explicit integration scheme, and the loss function, as expected from the Jacobian structure of Eq. ( 10 ) of the main text.

Panel of Figure 9 — Algorithm S1 Implicit time integration scheme

Figure S2: Short-term and long-term performance of neural autoregressive models constrained with implicit time integration. Color encodes the time step, $\Delta t=0.1$ in blue, $0.05$ in vermilion, $0.001$ in black. Implicit schemes are solid lines and filled pentagons, explicit schemes dotted lines and open stars. (a) Long-term emulation of the KS state $u(x,t)$ by the FNO model constrained with first-order implicit Euler at $\Delta t=0.05=50\,\Delta t_{\text{DNS}}$ , stable and physically plausible over $10^{5}$ steps of length $\Delta t_{\text{DNS}}$ . (b) Pointwise difference between the emulation in (a) and the reference numerical solution. The error stays small and develops slowly, with structure inherited from the flow rather than diverging, which quantifies the physical consistency. (c) RMSE growth on linear axes for the FNO models with implicit and explicit Euler at the three time steps. At each $\Delta t$ the implicit scheme gives lower error growth, and the linear scale makes the separation visible late in the rollout, where accumulated error dominates. The legend reports $|\lambda_{\max}|$ , closer to unity for the implicit scheme in each pair. Within each $\Delta t$ pair the ordering of $|\lambda_{\max}|$ matches the ordering of the RMSE curves, so $|\lambda_{\max}|$ acts as an a priori diagnostic of the error-growth rate. Comparison across pairs requires the normalized rate $\sigma$ of Eq. ( S1 ), reported in Section 2.3 . (d) Eigenvalues on the Argand plane for both schemes at $\Delta t=0.1$ and $0.05$ , with the unit circle dotted. At each $\Delta t$ the implicit spectrum has a smaller maximum and a more diffusive left tail than the explicit one, consistent with its lower error growth. Equation ( S6 ) accounts for both. The implicit map places an eigenvalue outside the unit circle only for expanding directions of the learned tangent operator, the explicit map even for mildly contracting ones.

Figure S3: Spectral bias and its connection to the eigenvalues of the model Jacobian. (a) Eigenvalues on the Argand plane for the MLP-based model (vermilion plus signs, $|\lambda_{\max}|=1.0003$ ) and the FNO-based model (blue open stars, $|\lambda_{\max}|=1.0002$ ), both constrained with the PEC4 integrator and trained with the spectral regularizer of Eq. ( 36 ) of the main text, with the unit circle dotted. The FNO-based model has the smaller $|\lambda_{\max}|$ and a longer left tail, that is, eigenvalues further inside the unit circle, indicating greater implicit diffusivity. This is reflected in (b) the spatial Fourier spectrum $|\hat{u}(\omega)|$ of the predicted state after $100$ time steps and (c) the spectrum of its time derivative $d|\hat{u}|/dt$ . Up to $\omega\approx 100$ both models track the reference numerical data (dotted black). In the shaded high-wavenumber band ( $\omega\gtrsim 100$ ), where the true spectrum continues to decay, both learned models instead saturate at a flat noise floor set by the smallest scale each can represent rather than by the true dynamics. The floor is markedly lower for the FNO (blue) than for the MLP (vermilion) in both the state and its derivative, so the FNO more faithfully captures the high-wavenumber tail. For this matched pair, $|\lambda_{\max}|$ thus tracks the implicit diffusivity and the high- $\omega$ noise floor, with the smaller $|\lambda_{\max}|$ corresponding to greater implicit diffusion and more faithful capture of the high-wavenumber tail of the spectrum.

Explore similar work

Apr 22, 2026

A Hybridizable Neural Time Integrator for Stable Autoregressive Forecasting

For autoregressive modeling of chaotic dynamical systems over long time horizons, the stability of both training and inference is a major challenge in building scientific foundation models. We present a hybrid technique in which an autoregressive transformer is embedded within a novel shooting-based mixed finite element scheme, exposing topological structure that enables provable stability. For forward problems, we prove preservation of discrete energies, while for training we prove uniform bounds on gradients, provably avoiding the exploding gradient problem. Combined with a vision transformer, this yields latent tokens admitting structure-preserving dynamics. We outperform modern foundation models with a

65\times

reduction in model parameters and long-horizon forecasting of chaotic systems. A "mini-foundation" model of a fusion component shows that 12 simulations suffice to train a real-time surrogate, achieving a

9{,}000\times

speedup over particle-in-cell simulation.

Brooks Kinch, Xiaozhe Hu, Yilong Huang +6

Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA, USA · Department of Mathematics, Tufts University, Medford, MA, USA · Department of Mathematics and Cybernetics, SINTEF Digital, Oslo, Norway +3

Neural Surrogate Modeling Time Series Transformers Chaotic Time-Series Prediction Autoregressive Modeling

Ask PDF

Aug 22, 2026

Symbolic Neural ODEs: Learning interpretable models from time-series data

We present a machine learning framework for identifying sparse, interpretable models of dynamical systems directly from time-series data. Our approach parameterizes the underlying vector field using a neural architecture and trains it by minimizing a multi-step prediction loss over a finite horizon. To ensure numerical tractability, we optimize a mean absolute error objective averaged across prediction steps, and progressively increase the horizon during training. A key feature of this formulation is that it enforces consistency under repeated composition of the learned dynamics. As a result, the identified models exhibit significantly improved stability compared with approaches based on one-step regression of the vector field. When combined with sparsity-promoting regularization, this leads to parsimonious models that generalize beyond the training data. We demonstrate accurate recovery of systems exhibiting a wide range of behaviors, including stable and unstable fixed points, periodic orbits, and chaotic attractors. For chaotic systems, while long-term trajectory prediction is inherently limited by sensitivity to initial conditions, we show that multi-step training yields models with accurate short-term dynamics and strong agreement in long-time statistical properties, including mean, variance, and Lyapunov exponents. Moreover, we establish theoretical bounds linking trajectory error to statistical accuracy, providing a step toward a principled explanation for this behavior.

Nibodh Boddupalli, Jeff Moehlis

Dept. of Mechanical Engineering, University of California, Santa Barbara, CA 93106

Nonlinear System Identification System Identification Neural ODEs Sparse Identification of Nonlinear Dynamics Chaotic Time-Series Prediction

Ask PDF

Apr 22, 2026

Learning to Emulate Chaos: Adversarial Optimal Transport Regularization

Chaos arises in many complex dynamical systems, from weather to power grids, but is difficult to accurately model with data-driven methods such as machine learning emulators. While emulators are promising tools for accelerating simulations and solving inverse problems, they still struggle to learn chaotic dynamics, where sensitivity to initial conditions renders exact long-term forecasts infeasible, especially given noisy data. Recent work instead trains emulators to match the statistical properties of chaotic attractors, but these approaches often rely on handcrafted summary statistics or large, diverse multi-environment datasets. In this work, we propose a family of adversarial optimal transport objectives that can jointly learn high-quality summary statistics and a physically consistent emulator from a single noisy trajectory. We theoretically analyze and experimentally validate a Sinkhorn divergence formulation (2-Wasserstein) and a WGAN-style dual formulation (1-Wasserstein) of our approach. Numerical experiments across a variety of chaotic systems, including ones with high-dimensional spatiotemporal chaos, show that emulators trained using our proposed objectives have significantly improved long-term statistical fidelity.

Gabriel Melo, Leonardo Santiago, Peter Y. Lu

LTCI, Télécom Paris, Institut Polytechnique de Paris, Palaiseau · Department of Mechanical and Aerospace Engineering, North Carolina State University, Raleigh, NC · University of Campinas +1

Neural Surrogate Modeling Chaotic Dynamical Systems Optimal Transport Chaotic Time-Series Prediction

Ask PDF