Do learned stochastic AI simulators respond correctly to forcing?

Why forecast skill cannot certify forcing response

Machine-learning emulators of chaotic and stochastic systems are usually validated on forecast skill and long-run statistics. The paper argues that neither certifies the property projection, attribution and tipping-point studies depend on: whether the emulator responds correctly to forcing. A model can reproduce the invariant measure closely and still misrepresent how that measure deforms under forcing.

Linear response theory makes this testable. For a large class of chaotic and stochastic systems, the forced response follows from unperturbed correlations through a generalized fluctuation-dissipation relation, and it decomposes over the stochastic Ruelle-Pollicott resonances of the Koopman generator and their eigenfunctions. The Koopmanism Response framework computes these objects from a single unperturbed trajectory with extended dynamic mode decomposition (EDMD), with ResDMD residuals guarding against spurious modes. "Did the model learn the physics?" then becomes falsifiable: does a surrogate trained only on unperturbed data carry the correct resonances, eigenfunctions and static susceptibility?

A calibrated, mode-resolved test per rollout

The paper turns response theory into a calibrated evaluation protocol for learned stochastic simulators. Each surrogate rollout passes or fails each check, and failure rates are compared with the false-alarm rate of independent realizations of the true system. Three families of checks are applied to every rollout: invariant statistics, the static susceptibility χ(0), and a mode-resolved eigenfunction test.

The testbed is stochastic Lorenz-63, a three-variable toy model of convection with no spatial structure, integrated by Euler-Maruyama at δt = 0.001 with (σ, ρ, β) = (10, 28, 8/3) and isotropic noise ς = 2. The x-equation is forced and x is tracked, because x is the most strongly coupled variable, driving both y and z through the nonlinear terms. From one unperturbed trajectory of T = 10000 time units sampled at τ = 0.1, a Fourier-EDMD approximation of the generator is built with dictionary size (2Kmax + 1)^3 = 729 and TSVD relative threshold 10^-4. The reference response is the instrument reconstruction on the held-out true realization (seed 101).

The authors are explicit that this is a controlled testbed with a known true response, not an operational climate emulator. Its purpose is to map how network design and training formulation shape response fidelity.

Invariant statistics and response fidelity dissociate

The central result is that the two kinds of fidelity come apart in both directions. The reservoir computer reproduces the invariant measure closely, with W1 inside its floor in 76 of 78 rollouts, and 26 of its rollouts lie inside all four invariant-statistics floors and the χ(0) band. Of those 26, 20 fail the slow-mode test, against 3 of the 48 truths that pass the same checks (Fisher p = 4.9×10^-10). The failure is not a disguised switching error: RC fails the slow modes in 70% of rollouts inside the switch-ratio floor and 80% outside it (p = 0.424). Spectrally, RC places its two leading oscillatory resonances within 0.133% of the reference frequencies, but its second slow mode decays faster than in any true realization in 77% of rollouts (median Reλ = -0.639 against the true range [-0.611, -0.483]), and its slow eigenfunctions pass in under a third of rollouts.

The neural ODE and neural SDE show the reverse pattern. Both over-switch in every rollout (median switch ratio 1.05 and 1.05) and fail invariant statistics in 93% and 100% of rollouts, yet they fail the slow modes in only 24% and 23%. An exploratory causal transformer fine-tuned on short rollouts, the architecture family of current weather emulators, shows the same split.

The MLP baselines fail every rung. The drift MLP learns the O(∆t) secant drift rather than F, and its Euler rollout collapses onto one wing in 33 of 80 rollouts. Fitted as a flow map, the same network restores two-wing switching (0 of 80 rollouts one-wing) but not the measure (median f+ = 0.474); its slowest real mode decays too slowly (-0.146 against the true -0.389), so only 33% of its χ(0) has accrued by t = 1, against a median 89% for the truths, and χ(0) overshoots.

Why χ(0) cannot separate these cases

For SINDy, RC, NODE and NSDE the χ(0) failure rate never exceeds 9%, while their slow-mode failure rates range from 0% to 76%. The paper explains this as expected for a time-integrated quantity: for forcing and observing x, χ(0) is set mostly by fast relaxation. Across the true realizations a median 89% of it has accrued by t = 1, before the slow modes (e-folding times 1.97 and 2.73) have relaxed.

This is why a scalar response summary such as χ(0) does not detect the dissociation here, while the mode-resolved test does. The paper notes that forcings or observables coupled to lobe occupation may load on exactly these slow modes; the resonances and eigenfunctions, unlike χ(0), do not depend on that choice.

The instrument itself is checked. With Γ = -∂x ln ρ0 and forcing along x, the value of the Green's function at zero lag is fixed analytically at Gx(0) = 1, and nothing in the pipeline normalizes Gx, so the measured deviation from unity is a free per-series estimate of dictionary-truncation error: |Gx(0) - 1| is at most 0.83% over the truths and at most 1.92% over all scored rollouts. The invariant-statistics floors are the maxima of each distance over the 50 independent true realizations scored against the held-out reference realization (seed 101) rather than the master trajectory, which is every model's training data; no floor fails its leave-one-out null test.

Training formulation decides which properties are learned

For a fixed network, the training formulation decides which of these properties it gets right. The paper builds a formulation ladder sharing architecture and data: one-step drift regression, flow map, multi-step neural ODE, and neural SDE with learned diffusion. This isolates the effects of the training target, the multi-step loss and learned noise.

A one-step drift regressed at a coarse step builds the O(∆t) secant bias into the model and breaks the attractor. One-step flow maps can match the attractor, as RC does, yet misplace the slow, weakly damped modes, plausibly because these change little over one step and barely enter a one-step loss. Training on multi-step windows through the integrator (NODE, NSDE) recovers the slow modes in three quarters of rollouts but over-switches.

A discretization control separates learning error from integrator error. Rolling the true drift through the same RK4 path at h = ∆tfit = 0.01 with oracle noise, 10 draws each, the switch ratio falls below one in every draw (median 0.953 and 0.948, maximum 0.975 and 0.968), whereas NODE and NSDE lie above one in 80 of 80 and 80 of 80 rollouts. The over-switching is therefore learned rather than numerical. The step shifts the median χ(0) by 0.0443 and 0.0549 from the truth median (0.221) and leaves the slow modes largely intact (slow-mode failure rate at most 10%).

What the positive control and thresholds reveal

SINDy with a quadratic library containing the true equations serves as the positive control. Trained on the full record, it recovers the structural coefficients to within 0.116% and fails at rates consistent with the truth on every rung, with no mode pass decided by the 0.99 cap. The instrument therefore certifies a correct model rather than only rejecting wrong ones.

Threshold sensitivity matters for reading the failure rates. The thresholds of the ω ≈ 8.42 and 16.68 pairs sit at the 0.99 cap (uncapped 0.9984 and 0.9960), so passes on these modes can be decided by the cap, as for 75 of RC's 77 passes at 8.42 and 63 and 71 of the 80 and 80 NODE and NSDE passes there, against 0 for SINDy. Rebuilding the two slow-mode thresholds from the 50 fresh truths (leave-one-out minimum: 0.9257 and 0.8693 instead of 0.9450 and 0.9216) lowers RC's slow-mode failure rate from 76% to 55% and NODE's from 24% to 7%, and the RC exhibit becomes 14 of 26 rollouts against 1 of 48 truths (Fisher p = 2.6×10^-7).

Per-mode overlaps also show that RC, NODE and NSDE miss the individual eigenvectors of the near-degenerate oscillatory pair in most rollouts while recovering their joint subspace; the MLPs recover neither, and only SINDy resolves the split.

Where the probe applies and where it stops

The test is calibrated, not exhaustive: one low-dimensional system, one forcing direction and one observable. Thresholds come from 16 calibration truths, and rebuilding the slow-mode thresholds from the 50 fresh truths changes the reported failure rates, as does the 0.99 cap on oscillatory passes. Except for the NSDE, rollouts use the true noise intensity, an oracle that favours the surrogates; as a robustness check, SINDy fitted to only 100 time units acquires symmetry-breaking terms.

Within those bounds, the paper's practical priority is clear. Neither invariant statistics nor χ(0) certifies the slow dynamics: a quarter of RC's rollouts pass both with the wrong slow modes. Multi-step losses, now common in weather emulators, appear to carry the slow dynamics, and could be paired with invariant-measure constraints or spectral priors rather than with more flexible architectures. Because slow modes govern transitions between the quasi-invariant sets of a metastable system and, for forcings that couple to them, its response, a mode-resolved check needing one unperturbed trajectory per model should accompany invariant statistics in validating learned simulators. The natural next step is to test whether these patterns persist in spatially extended systems and for state-of-the-art architectures of current emulators.

Sources used in this answer