cs.AIAug 17, 2026
SaveEigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics
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 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 (Section 3.4 ).
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 steps. The two direct-step models (Direct MLP, Direct FNO) grow rapidly and saturate near the climatological decorrelation level, where the RMSE is in these normalized units. The two integration-constrained Euler models (Euler MLP, Euler FNO), whose eigenvalue structure is representative of all 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 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 at the initial condition on the Argand plane for all four models, with in the legend. The Direct MLP ( ) and Direct FNO ( ) admit eigenvalues with ( and ) spread well inside and beyond the unit circle, whereas the Euler MLP ( ) and Euler FNO ( ) eigenvalues collapse onto the unit circle (dotted) and cluster tightly at ( ), 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 for the Euler MLP and 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 for every model.
Figure 3: A single linear scaling law governs error growth across the full model suite. Norm of the accumulated error plotted against for all models, where is the largest eigenvalue of the model Jacobian evaluated at the initial condition and is the Euclidean ( ) norm of the error vector. Here , so is the first-step (generalization) error and the lead time is measured from the first step; the reference slope at lead time is therefore the step count of Eq. ( 16 ). Panels (a)–(d) show the successive lead times , 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 in (a) to slopes , , and in (b)–(d). The panels correspond to , , , and , so they span the validity window of Eq. ( 16 ) from well inside it to its boundary. The integration-constrained models ( ) 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 and the derivation no longer applies. The two direct-step models sit at much larger and rise above the reference line at the lead time of (b), reflecting the faster-than-linear growth driven by their ; 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 . This faster-than-linear growth arises from the higher-order terms in the Taylor expansion of , 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 and therefore exceed the 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 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 of Eq. ( 36 ). (a) RMSE over the autoregressive rollout on log–log axes. For each model the solid line is the total error and the black dashed line is the cumulative generalization error , 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 for the FNO ( ) and Spectral FNO ( ), with the unit circle dotted; to the precision of the legend the two models share . Each marker is shaded by the cosine similarity between its eigenvector and the error vector (grayscale bar). (c, d) The linearized error law of Eq. ( 20 ), evaluated with from the Jacobian at the initial condition, tested after rollout steps ( ) 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 . (c) Total error, with the measured one-step error. (d) The same law restricted to the high-wavenumber band ( , with ), with the injection term replaced by , the high- 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 , versus , 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 (gray), confirming that the new loss preserves the linear scaling law, and for every integrator the stability-promoting loss yields a smaller and correspondingly smaller 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 , the slope of a linear fit of against of the timestep over the first rollout steps (the exponent of ). For every integrator, the stability-promoting loss produces a smaller , 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 | RMSE ( ) | Spec. bias | |||
| 1 | MLP | Euler | RMSE | |||||
| 2 | MLP | Euler | Spectral | |||||
| 3 | MLP | RK4 | RMSE | |||||
| 4 | MLP | RK4 | Spectral | |||||
| 5 | MLP | PEC4 | RMSE | |||||
| 6 | MLP | PEC4 | Spectral |
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 ; the two large- implicit models (26, 27) are the exception. is the largest-magnitude eigenvalue of the model Jacobian at the initial condition. is the corresponding growth rate per unit time of Eq. ( S1 ), which is comparable across ; its precision is limited by the five decimal places reported for . RMSE is the short-term error at steps. The spectral-bias metric is 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 , 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 , rather than in slope, consistent with the near-neutral spectra shared across the suite. The window shown satisfies , so it lies inside the validity range of that scaling law. (b) Eigenvalues of the model Jacobian , evaluated at the initial condition , on the Argand plane, with the unit circle dotted. For every model the eigenvalues collapse onto the unit circle and cluster tightly at , with the spread of each spectrum confined to the third decimal place and the models separated in the fourth (Table S1 reports 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.
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, in blue, in vermilion, 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 by the FNO model constrained with first-order implicit Euler at , stable and physically plausible over steps of length . (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 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 , closer to unity for the implicit scheme in each pair. Within each pair the ordering of matches the ordering of the RMSE curves, so acts as an a priori diagnostic of the error-growth rate. Comparison across pairs requires the normalized rate of Eq. ( S1 ), reported in Section 2.3 . (d) Eigenvalues on the Argand plane for both schemes at and , with the unit circle dotted. At each 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, ) and the FNO-based model (blue open stars, ), 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 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 of the predicted state after time steps and (c) the spectrum of its time derivative . Up to both models track the reference numerical data (dotted black). In the shaded high-wavenumber band ( ), 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, thus tracks the implicit diffusivity and the high- noise floor, with the smaller corresponding to greater implicit diffusion and more faithful capture of the high-wavenumber tail of the spectrum.