arXiv CSopen access

Teacher Forcing as Generalized Bayes: Optimization Geometry Mismatch in Switching Surrogates for Chaotic Dynamics

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Teacher Forcing as Generalized Bayes: Optimization Geometry Mismatch in Switching Surrogates for Chaotic Dynamics Andre Herz1,∗

Daniel Durstewitz1,2

Georgia Koppe1,3,4

1 Interdisciplinary Center for Scientific Computing, Heidelberg University, Germany Dept. of Theoretical Neuroscience, Central Institute of Mental Health (CIMH), Mannheim, Germany 3 Hector Institute for AI in Psychiatry and Dept. of Psychiatry and Psychotherapy, CIMH, Mannheim, Germany 4 Hertie Institute for AI in Brain Health, University of Tübingen, Germany ∗ Correspondence to: [email protected]

arXiv:2604.25904v1 [cs.LG] 28 Apr 2026

Abstract Identity teacher forcing (ITF) enables stable training of deterministic recurrent surrogates for chaotic dynamical systems and has been highly effective for dynamical systems reconstruction (DSR) with recurrent neural networks (RNNs), including interpretable almostlinear RNNs (AL-RNNs). However, as an intervention-based prediction loss (and thus a generalized Bayes update), teacher forcing need not match the free-running model’s marginal likelihood geometry. We compare the objective-induced curvatures of ITF and marginal likelihood in a probabilistic switching augmentation of AL-RNNs, estimating ambiguity-aware observed information via Louis’ identity. In the switching setting studied here, conditioning on a single forced regime path (as ITF does) inflates curvature, while marginal likelihood curvature is reduced by a missing-information correction when multiple switching explanations remain plausible. In Lorenz–63 experiments, windowed evidence fine-tuning improves held-out evidence but can degrade dynamical quantities of interest (QoIs) relative to ITF-pretrained models.

Introduction

In dynamical systems reconstruction (DSR), one aims to learn surrogate models of unknown dynamical systems from time series data. The central goal is to recover the long-term behavior of the underlying system, beyond short-horizon prediction, commonly assessed via dynamical invariants and statistics such as Lyapunov exponents (LEs), the physical invariant measure This work was presented at the Workshop on Optimisation and Post-Bayesian Inference in Machine Learning at AISTATS 2026.

governing time-averages, and topological and geometric properties of the attractor (Durstewitz et al., 2023; Göring et al., 2024). In practice, such models are typically trained using short-horizon prediction objectives. In chaotic systems, however, short-horizon predictive fit can be statistically decoupled from long-run invariant structure (Park et al., 2024). Positive LEs cause small state errors to grow exponentially, rendering trajectory-level gradients propagated through backpropagation through time (BPTT) ill-conditioned due to exploding gradients (Mikhaeil et al., 2022). As a result, direct long-rollout training of recurrent neural network (RNN) surrogates is numerically unstable. This creates a training paradox: we seek long-run invariant behavior, yet optimization breaks down at long horizons. A common workaround relies on sparse teacher forcing (STF) (Mikhaeil et al., 2022). In STF, the latent state of a deterministic model is periodically fed a control signal derived from observations, anchoring rollouts to the data and thereby stabilizing gradient propagation (Mikhaeil et al., 2022). This enables effective learning even when the autonomous (free-running) model would quickly diverge. Empirically, models trained under this principle have demonstrated striking success in reconstructing chaotic attractors, including almost-linear RNNs (AL-RNNs) with minimal nonlinear structure (Hess et al., 2023; Brenner et al., 2024). Nevertheless, STF optimizes prediction error under an explicit intervention that alters the rollout distribution, meaning it generally does not correspond to maximumlikelihood training of the unforced generative model. From a statistical viewpoint, STF instead defines a task loss that can be interpreted as a generalized Bayes (Gibbs) posterior, inducing its own local geometry in parameter space (Catoni, 2007; Bissiri et al., 2016). Under this generalized framework, this local geometry is fundamental: it underpins second-order optimization,

Laplace uncertainty quantification, local identifiability analysis, and optimal experimental design. By contrast, stochastic state-space models (SSMs) start from an explicit autonomous generative model and define a marginal likelihood by integrating over latent variables (Durbin and Koopman, 2001). The local geometry of the resulting posterior is summarized by the observed information. For latent-variable SSMs, this can be computed via Louis’ identity (Louis, 1982). As a modeling stance consistent with representing probabilistic beliefs, stochastic SSMs provide a principled default for real-world time series, since partial observability, measurement noise, and unmodeled perturbations are ubiquitous. Yet deterministic STF still dominates many chaotic DSR benchmarks (Hess et al., 2023). We compare the local geometries induced by STF and marginal likelihood-based learning via ambiguity-aware observed information in a probabilistic switching SSM. AL-RNNs provide an interpretable testbed with explicit switching codes and tractable (closed-form) Fisherinformation expressions aligned with ITF training. We then test whether windowed marginal likelihood optimization aligns with long-horizon dynamical invariants. Related Work. Chaotic DSR surrogates are commonly trained with sparse/generalized teacher forcing to stabilize BPTT (Mikhaeil et al., 2022; Hess et al., 2023; Brenner et al., 2024). Probabilistic switching SSMs (e.g., SLDS/rSLDS) instead learn by marginal likelihood objectives with approximate inference/EM/SMC (Ghahramani and Hinton, 2000; Murphy, 1998; Linderman et al., 2017). We connect these views by (i) treating STF as a generalized Bayes update (Catoni, 2007; Bissiri et al., 2016) and (ii) estimating ambiguity-aware likelihood curvature via Louis’ identity (Louis, 1982).

Two Posteriors and Local Curvatures

2.1

AL-RNNs and Switching Dynamics

AL-RNNs combine linear recurrent structure with a small number of gated nonlinear units, resulting in dynamics that switch between linear regimes (Brenner et al., 2024). Let zt ∈ RM be the latent state at time t. The AL-RNN transition model is zt+1 = Fθ (zt ) := Azt + W ϕ∗ (zt ) + h,

(1)

where A, W ∈ RM ×M are diagonal and full, respectively, and h ∈ RM is a bias term. ϕ∗ is identity on the first M − P coordinates and ReLU on the last P . The architecture induces switching dynamics. Define the binary switching code ct ∈ {0, 1}P by ct,j = 1{zt,M −P +j > 0}, which specifies the active linear regime at time t. Let D(ct ) denote the diagonal

matrix that applies the corresponding gates: D(ct ) = diag 1, . . . , 1, ct,1 , . . . , ct,P . | {z }

(2)

M −P

Then the dynamics are piecewise affine and can be written as zt+1 = A + W D(ct ) zt + h. (3) The induced switching code (gate) sequence {ct } defines a symbolic partition of state space and determines which linear regime governs local evolution (Brenner et al., 2024). Near switching boundaries, multiple gate configurations may be locally plausible, giving rise to latent uncertainty about which linear regime generated the data. We will empirically show that whether or not this uncertainty is properly accounted for has direct consequences for the local curvature of the training objective. 2.2

ITF Loss Posterior and Local Curvature

As training protocol, we use identity TF (ITF), a variant of STF where the first N latent coordinates (the observed subspace) are overwritten at forcing times Tτ (Mikhaeil et al., 2022; Brenner et al., 2022). Let B = [IN 0] ∈ RN ×M denote the projection onto observed coordinates such that x̂t = Bzt , and let Tτ = {t : t ≡ 0 (mod τ ), t > 0} denote the forcing times. The ITF-modified rollout can be written schematically as ( Fθ (z̃t ), t ∈ Tτ , T z̃t = zt + B xt − Bzt , zt+1 = Fθ (zt ), else. (4) STF can also be generalized using a forcing strength parameter to linearly interpolate between zt and z̃t (Hess et al., 2023). Here, we focus on “strict” forcing. Training based on BPTT then minimizes the one-step prediction loss along an ITF-modified trajectory: LITF (θ; x1:T ) =

T −1 1 X 2 Bzt+1 − xt+1 2 . T − 1 t=1

(5)

This loss evaluates prediction error under an intervention that alters the system’s dynamics during training. It therefore does not correspond to the maximum likelihood estimator for the autonomous generative model one uses at test time (i.e., free-running dynamics). A useful framework for interpreting this objective statistically is generalized Bayes (Catoni, 2007; Bissiri et al., 2016). In chaotic DSR, global likelihoods are brittle, and any surrogate model is necessarily misspecified relative to the underlying physical system. Rather than assuming a single canonical likelihood, it is therefore

natural to treat training objectives as scoring rules that induce generalized posteriors and, consequently, distinct local parameter geometries. Interpreting LITF as such a scoring rule, ITF corresponds to the Gibbs posterior πITF (θ | x) ∝ π0 (θ) exp − β LITF (θ; x) , (6) where π0 is a prior and β > 0 is an inverse-temperature (loss-scale) parameter. We define the ITF curvature proxy as the generalized per-step Gauss–Newton/Fisher matrix aligned to ITF: T −1 1 X T −1 I ITF = J Λ Jt , T − 1 t=1 t

∂(Bzt+1 ) Jt := , (7) ∂θ

where Λ is a chosen loss weighting. At forcing times Tτ , the sensitivity recursion is modified to reflect the overwrite of the observed coordinates (Appendix A1.3). Crucially, Jt is evaluated along the forced rollout’s induced switching code sequence {ct }. Thus, the local geometry is conditioned on a single code path and does not marginalize over switching ambiguity (alternative plausible regime paths). All curvature diagnostics we report are data curvature only (no prior Hessian added). Since regularization is fixed across objectives, comparisons isolate objectiveinduced curvature. 2.3

Marginal Likelihood Curvature and Observed Information via PAL-RNN

To obtain a marginal likelihood-based reference geometry, we probabilistically augment the AL-RNN into a hybrid switching SSM (PAL-RNN). Instead of treating the switching codes ct as deterministic functions of zt , we model them as latent random variables governed by probabilistic gating. Starting from a switching linearGaussian SSM: zt+1 | zt , ct ∼ N A + W D(ct ) zt + h, Q , (8) xt | zt ∼ N Bzt , R , we introduce probit gate noise: zt,M −P +j ct,j | zt ∼ Bernoulli Φ , σg

(9)

where Φ is the standard normal CDF and σg > 0 is a hyperparameter that controls gate stochasticity. Q and R are process/observation-noise covariances. In 2 our experiments we use isotropic noise, i.e., Q = σproc I 2 and R = σobs I on the relevant state/observation dimensions. We further interpret the squared error in (5) as a Gaussian negative log-pseudo-likelihood with covariance R, so R = Λ, the loss weight in the ITF

curvature (7). Under the PAL-RNN formulation, the local curvature of the marginal likelihood pθ (x1:T ) of an observed time series is given by the observed information, which, by Louis’ identity (Louis, 1982), decomposes into a complete-data term minus a missing-information correction arising from latent uncertainty: I obs (θ) = E[I comp (θ) | x] − Cov(∇θ log p(x, z, c | θ) | x), | {z }

(10)

I miss (θ)⪰0

where I comp (θ) := −∇2θ log pθ (x, z, c) is the completedata information. Missing information is large when multiple switching code explanations remain plausible a posteriori. Under a Laplace approximation at θ̂, I obs (θ̂) provides the marginal likelihood contribution to posterior precision. Thus, objective curvature dictates local Gaussian uncertainty under both likelihood and loss-based frameworks. Next, we compare the ITFinduced pseudo-likelihood geometry with the marginal likelihood-based observed-information geometry and examine their implications for local identifiability and long-horizon DSR.

Results

Missing Information vs. Posterior Entropy To isolate the Louis mechanism in a setting where switching ambiguity can be varied without changing the underlying linear regimes, we study a probit-gated switching AR(1) toy model (Appendix A1.1 and Appendix A1.2), and sweep the gate-noise scale σg . For each run, we relate the (time-averaged) mean posterior gate entropy to the missing-information ratio (MIR := 1 − tr(I obs )/tr(E[I comp | x1:T ])), and the observed curvature proxy log10 tr(I obs ). Figure 1(a) shows that increasing entropy increases MIR and decreases observed curvature. This validates that posterior switching ambiguity flattens the marginal likelihood geometry through the Louis missing-information term. ITF-aligned Curvature Exceeds AmbiguityAware Observed Information. To examine whether ITF yields sharper local curvature than ambiguity-aware observed information under switching ambiguity, we compare ITF-aligned curvature to ambiguity-aware observed information. We test this on AL-RNNs trained under ITF on Lorenz–63 data (Lorenz, 1963). For each model with drift parameters θ̂ = {A, W , h}, we construct a probabilistic counterpart by applying the PAL-RNN formulation (8, 9) and adding isotropic process/observation and probit gate noise (Sec. 2). We compute the ITF-aligned curvature

Figure 1: Summary of results. (a) In a probit-gated switching AR(1) model, increasing gate noise σg leads to higher posterior switching ambiguity. This increases both mean posterior gate entropy and missing-information ratio (MIR), while decreasing observed-curvature proxy log10 tr(I obs ). Faint points = individual runs; markers = mean ± SEM across 20 seeds. (b) Curvature gap gQ (measuring how much ITF curvature exceeds ambiguity-aware observed information) vs. filtering switching code entropy Hc evaluated at the same drift parameters. The gap indicates that conditioning on a single forced switching path inflates curvature under switching ambiguity. Points are colored by observation-noise level σobs . (c) Particle-SAEM fine-tuning: held-out windowed conditional log evidence (mean over windows; normalized per step and per observed dimension) and two hard-gated QoIs: Dstsp b1 − λref (negative = overly contracting dynamics). calib updates process/observation noise and signed LE error λ 1 only (QoIs unchanged); full SAEM updates drift+noise, improves evidence the most, but degrades both QoIs. I ITF (θ̂) along the forced rollout (4). Using a Rao– Blackwellized particle smoother, we estimate I obs (θ̂) via Louis’ identity (10) while sweeping over noise levels (see Appendix A1.2 for details). We quantify the mismatch by the local curvature gap gQ = log10 (T tr(I ITF ) / tr(I obs )) ,

(11)

which estimates curvature inflation under ITF for a fixed process-noise scale Q. Relating gQ to (timeaveraged) filtering switching code entropy Hc , Figure 1(b) shows gaps ranging from ≈ 10 up to ≈ 103 that tend to increase with Hc . Appendix A4.3 shows that this curvature gap is not a uniform matrix dominance: most generalized-eigenvalue directions are not sharper under ITF, while trace summaries and leadingsubspace overlaps still indicate an anisotropic, subspacedependent geometry mismatch. Evidence–QoI Misalignment To assess whether higher evidence implies better dynamics, we finetune PAL-RNNs initialized from ITF-trained ALRNN checkpoints by maximizing windowed conditional log evidence with particle stochastic approximation expectation–maximization (particle-SAEM; Estep via a Rao–Blackwellized particle smoother; Appendix A2.3). We compare (i) baseline (no finetuning), (ii) calib (process/observation noise only), and (iii) full SAEM (noise + drift). Across 12 ITFpretrained initializations and five window lengths (detailed settings in Appendix A2), we evaluate held-out evidence and QoIs: overlap in state distribution Dstsp

b1 − λref (QoI details in Apand largest LE error λ 1 pendix A3). Hard-gated QoIs use deterministic rollouts with thresholded gates and zero noise. Figure 1(c) reveals an evidence–QoI misalignment across models and window lengths. Although both calib and full SAEM increase held-out evidence (with calib leaving the hard-gated QoIs unchanged by construction), updating drift parameters (full SAEM) substantially degrades invariant reconstruction. In particular, it drives the signed largest Lyapunov error b1 − λref from near zero to strongly negative, indicatλ 1 ing that windowed evidence trades chaos for stability while distorting long-run invariant structure.

Implications for Post-Bayesian DSR

Local geometry matters in chaotic DSR: it drives local optimization, uncertainty quantification, identifiability, and experimental design. In switching surrogates, it is strongly objective-dependent: teacher-forced curvature conditions on a single forced switching path, whereas ambiguity-aware observed information flattens curvature when multiple regime explanations remain plausible. Windowed marginal likelihood fine-tuning can also be misaligned with long-horizon goals: evidence may improve while dynamical invariants degrade when drift parameters are updated. Together, these findings suggest there is no canonical posterior for chaotic DSR. The update and geometry should be chosen for the target QoI, motivating QoI-aware post-Bayesian learning and geometry-aware active data collection.

Acknowledgments This work was funded by the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments, by the Federal Ministry of Research, Technology and Space (BMFTR) under the neuroAI initiative (01GQ2509B), by the German Research Foundation (Du 354/15-1), and the Hector II foundation. References Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory. Meccanica, 15(1):9–20. Bissiri, P. G., Holmes, C., and Walker, S. (2016). A general framework for updating belief distributions. Bayesian Analysis, 11(3):677–698. Brenner, M., Hemmer, C. J., Monfared, Z., and Durstewitz, D. (2024). Almost-linear RNNs yield highly interpretable symbolic codes in dynamical systems reconstruction. In Advances in Neural Information Processing Systems (NeurIPS). arXiv:2410.14240. Brenner, M., Hess, F., Mikhaeil, J. M., Bereska, L. F., Monfared, Z., Kuo, P.-C., and Durstewitz, D. (2022). Tractable dendritic RNNs for reconstructing nonlinear dynamical systems. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research. PMLR. Brenner, M. and Koppe, G. (2026). Uncovering the computational roles of nonlinearity in sequence modeling using almost-linear RNNs. Transactions on Machine Learning Research. Catoni, O. (2007). PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning, volume 56 of Institute of Mathematical Statistics Lecture Notes–Monograph Series. Institute of Mathematical Statistics. Chen, R. and Liu, J. S. (2000). Mixture Kalman filters. Journal of the Royal Statistical Society: Series B, 62(3):493–508. Cheng, X., He, Y., Yang, Y., Xue, X., Cheng, S., Giles, D., Tang, X., and Hu, Y. (2025). Learning chaos in a linear way. In International Conference on Learning Representations. Delyon, B., Lavielle, M., and Moulines, E. (1999). Convergence of a stochastic approximation version of the EM algorithm. The Annals of Statistics, 27(1):94– 128. Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the

EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1):1–38. Doucet, A., de Freitas, N., Murphy, K., and Russell, S. (2000a). Rao–blackwellised particle filtering for dynamic bayesian networks. In Proceedings of the 16th Conference on Uncertainty in Artificial Intelligence (UAI). Doucet, A., Godsill, S., and Andrieu, C. (2000b). On sequential monte carlo sampling methods for bayesian filtering. Statistics and computing, 10(3):197–208. Durbin, J. and Koopman, S. (2001). Time series analysis by state space methods. Oxford University Press, Oxford; New York. Durstewitz, D., Koppe, G., and Thurm, M. I. (2023). Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24(11):693–710. Eisenmann, L., Brändle, A., Monfared, Z., and Durstewitz, D. (2026). Detecting invariant manifolds in relu-based rnns. arXiv:2510.03814. Eisenmann, L., Monfared, Z., Göring, N., and Durstewitz, D. (2023). Bifurcations and loss jumps in RNN training. In Advances in Neural Information Processing Systems, volume 36. Curran Associates, Inc. Emonds, N., Herberg, E., Gerchen, M. F., Pritsch, M., Rocha, J., Zamoscik, V., Kirsch, P., Herzog, R., and Koppe, G. (2025). A data-driven closed-loop control approach to drive neural state transitions for mechanistic insight. bioRxiv. Ghahramani, Z. and Hinton, G. E. (2000). Variational learning for switching state-space models. Neural Computation, 12(4):831–864. Göring, N. A., Hess, F., Brenner, M., Monfared, Z., and Durstewitz, D. (2024). Out-of-domain generalization in dynamical systems reconstruction. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 16071–16114. PMLR. Hess, F., Monfared, Z., Brenner, M., and Durstewitz, D. (2023). Generalized teacher forcing for learning chaotic dynamics. In Proceedings of the 40th International Conference on Machine Learning (ICML), volume 202 of Proceedings of Machine Learning Research, pages 13017–13049. arXiv:2306.04406. Kutta, W. (1901). Beitrag zur näherungsweisen integration totaler differentialgleichungen. Zeitschrift für Mathematik und Physik, 46:434–453. Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., and Paninski, L. (2017). Bayesian Learning and Inference in Recurrent Switching Linear Dynamical Systems. In Singh, A. and Zhu, J., editors, Proceedings of the 20th International Conference on

Artificial Intelligence and Statistics, volume 54 of Proceedings of Machine Learning Research, pages 914–922. PMLR. Liu, L., Jiang, H., He, P., Chen, W., Liu, X., Gao, J., and Han, J. (2020). On the variance of the adaptive learning rate and beyond. In International Conference on Learning Representations. Lorenz, E. N. (1963). Deterministic nonperiodic flow. Journal of the Atmospheric Sciences, 20(2):130–141. Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society: Series B, 44(2):226–233. Mikhaeil, J. M., Monfared, Z., and Durstewitz, D. (2022). On the difficulty of learning chaotic dynamics with rnns. In Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS ’22, Red Hook, NY, USA. Curran Associates Inc. Murphy, K. P. (1998). Switching kalman filters. Technical Report 98-10, Compaq Cambridge Research Laboratory. Park, J., Yang, N. T., and Chandramoorthy, N. (2024). When are dynamical systems learned from time series data statistically accurate? In Advances in Neural Information Processing Systems, volume 37, pages 43975–44008. Curran Associates, Inc. Schiff, Y., Wan, Z. Y., Parker, J. B., Hoyer, S., Kuleshov, V., Sha, F., and Zepeda-Núñez, L. (2024). DySLIM: Dynamics stable learning by invariant measure for chaotic systems. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 43649–43684. PMLR. Sparrow, C. (1982). The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors, volume 41 of Applied Mathematical Sciences. Springer New York, New York, NY.

Appendix for Teacher Forcing as Generalized Bayes: Optimization Geometry Mismatch in Switching Surrogates for Chaotic Dynamics

Contents Appendix A0: Dynamical systems reconstruction (DSR) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Appendix A1: Methodological details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Appendix A2: Datasets and experimental settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Appendix A3: Metric definitions and computation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Appendix A4: Additional analyses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Appendix A5: Limitations and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Dynamical systems reconstruction (DSR)

This section provides additional details on the field of DSR, its dynamical quantities of interest (QoIs), and the challenges of training surrogates for chaotic systems. It also introduces the AL-RNN surrogate model used in this work and motivates why objective geometry matters in DSR. A0.1

DSR overview

The goal of DSR is to learn a generative surrogate model of an unknown time-evolving process from observed time series. Unlike conventional forecasting, the scientific target is typically not only to remain aligned with one particular trajectory, but to recover the underlying dynamical mechanisms well enough to reproduce long-run behavior and support downstream analysis via interpretable model architectures. A0.1.1

A minimal formalization

Throughout, we focus on discrete time, t ∈ Z. A deterministic DS is then expressed by an iterated map zt+1 = F ⋆ (zt ),

zt ∈ Z ⊆ RM ,

(A1)

where Z is the (compact, measurable) state space of the system and observations are generated through an (often non-invertible) measurement operator g, xt = g(zt ) + ηt , (A2) with observation noise ηt . DSR trains a generative model to approximate a surrogate map Fθ (and possibly a surrogate observation operator) so that the induced dynamics match the true system in an appropriate DS sense. An idealized target is topological conjugacy: there exists a homeomorphism h such that h(F ⋆ (z)) = Fθ (h(z)). In practice this is often too strong globally. Instead, one typically aims for agreement restricted to the attractor (and its basin of attraction) that is supported by the observed data, and assesses it through dynamical QoIs (Durstewitz et al., 2023; Göring et al., 2024). An attractor is a forward-invariant set that attracts a neighborhood of initial conditions in state space. In applications, we only observe trajectories drawn from (or quickly converging to) this attracting set, so reconstruction is evaluated on the portion of state space that carries mass under the corresponding invariant measure of this set. Conventions used throughout this work. We follow the state-space notation zt ∈ RM (latent state) and denote observations by xt ∈ RN . We choose the measurement operator g to be a linear projection xt = Bzt + ηt ,

ηt ∼ N (0, R).

(A3)

For identity teacher forcing (ITF), we use B = [IN 0] ∈ RN ×M to select the first N coordinates of the state as observations. At forcing times Tτ = {t : t ≡ 0 (mod τ ), t > 0}, ITF applies the forcing intervention z̃t = zt + B T (xt − Bzt ),

(A4)

which overwrites the observed coordinates with the data while leaving the unobserved coordinates unchanged. A0.2

Measuring success: dynamical QoIs

The success of DSR is typically judged by whether the free-running surrogate reproduces the system’s typical long-run behavior on the attractor, rather than by trajectory-level alignment. Under standard ergodicity/mixing assumptions, long-run statistics (time averages of observables) can be expressed as expectations under the system’s invariant measure. Concretely, an invariant measure µ is a stationary distribution over states that satisfies µ(F ⋆−1 (Z)) = µ(Z) for all measurable sets Z ⊆ Z. Let µ denote the invariant measure of the true system on the attractor of interest, and µθ the invariant measure induced by the surrogate. Matching long-run properties can then be expressed as agreement of expectations of relevant observables ψ, Ez∼µ [ψ(z)] ≈ Ez∼µθ [ψ(z)].

(A5)

In practice, agreement in Eq. (A5) is assessed by comparing long observed sequences to long free-running surrogate rollouts. Examples of ψ include (i) overlap in long-run state occupancy (a proxy for agreement between invariant measures), (ii) temporal summary statistics (e.g., power spectrum), and (iii) Lyapunov exponents (LEs). Typically, one discards an initial transient (“burn-in”) from rollouts to reduce dependence on initialization.1 A0.3

Challenges of DSR for chaotic systems

Chaotic systems exhibit sensitive dependence on initial conditions: first-order perturbations typically grow like δzT +1 ≈

T Y

Jt δz1 ,

where Jt := ∂F ⋆ (zt )/∂zt .

(A6)

t=1

The largest LE of the system (“maximal LE”) is then given by the exponential growth rate of these Jacobian products, T Y 1 λmax = lim sup log ∥ Jt ∥. (A7) T →∞ T t=1 Chaotic systems have a positive maximal LE, λmax > 0, so small initial errors grow exponentially, making long-term trajectory alignment impossible in practice. This implies a finite effective forecast horizon even for a correct model. This has two important consequences. First, trajectory-alignment losses (e.g. MSE over long horizons) become uninformative for DSR on chaotic systems. Second, the same Jacobian products that drive chaotic expansion also appear in backpropagation through time (BPTT) and real time recurrent learning, creating exploding gradients that make naive long-horizon training unstable (Mikhaeil et al., 2022). A common rough forecast-horizon scale is the Lyapunov time tλ := 1/λmax , i.e. a characteristic time over which small errors amplify by a factor e. A0.4

Teacher forcing as intervention-based training

To alleviate the ill-conditioning of BPTT in chaotic DSR, modern DSR training methods introduce controlled interventions to stabilize gradients for recurrent neural networks (RNNs) while still training the model to generate plausible long rollouts. Sparse teacher forcing (STF) and generalized teacher forcing (GTF) are two prominent examples (Mikhaeil et al., 2022; Hess et al., 2023). In GTF, one interpolates between model state and a data-inferred state ẑt = (1 − αt )zt + αt z̃t , z̃t ≈ E(xt ), (A8) 1

Early in a rollout, the state distribution can still depend strongly on the initial condition (before the trajectory has mixed onto the attracting set). Discarding an initial “burn-in” reduces this initialization bias so that empirical averages better approximate invariant-measure expectations.

and advances the dynamics using this interpolated state (e.g., zt+1 = Fθ (ẑt )), with αt ∈ [0, 1] controlling forcing strength (typically, αt ≡ α and forcing is applied at every step t). Here E denotes an encoder mapping observations to a surrogate latent state. STF (including ITF) can be viewed as a sparse special case where the intervention is applied only at forcing times, together with a partial overwrite of observed coordinates. Both schemes can be understood as optimizing a loss under an intervention policy that modifies the autonomous dynamics during training. Crucially, the geometry (curvature) of the resulting training objective need not match the geometry of any likelihood for the free-running model. A0.5 A0.5.1

AL-RNN surrogates for DSR Piecewise-linear RNNs for DSR

DSR models are typically not only prediction tools but also objects of downstream analysis. Piecewise-linear RNNs have been shown to be an attractive choice because they can represent complex nonlinear dynamics while retaining locally linear structure. This structure leads to closed-form Jacobians and dynamical objects (e.g., fixed points, cycles), thus enabling mechanistic interpretation of the learned dynamics (Eisenmann et al., 2023, 2026). The interpretability of these models has been exploited to analyze the dynamics of trained surrogates and to relate them to known mechanisms in the target system, e.g. in neuroscience applications (Durstewitz et al., 2023; Emonds et al., 2025). A0.5.2

AL-RNNs

Almost-linear RNNs (AL-RNNs) particularly emphasize interpretability by concentrating the nonlinearity into a small block of P < M ReLU-gated coordinates, yielding a parsimonious representation of nonlinear dynamics with 2P linear subregions and an explicit symbolic regime code (Brenner et al., 2024; Brenner and Koppe, 2026). Formally, the AL-RNN models the transition map by zt+1 = Fθ (zt ) := Azt + W ϕ∗ (zt ) + h,

(A9)

using a linear part A ∈ RM ×M (diagonal), a gated part W ∈ RM ×M (full), and an affine bias h ∈ RM . The activation ϕ∗ applies ReLU only on the last P coordinates and is identity elsewhere. While external inputs can also be included into the state transition, we focus on autonomous dynamics throughout this work. Let ct ∈ {0, 1}P be the switching code that records which of the last P coordinates are active (positive) at time t. This allows us to rewrite the transition as a switching affine system, zt+1 = A + W D(ct ) zt + h, D(ct ) = diag 1, . . . , 1, ct,1 , . . . , ct,P . (A10) | {z } M −P

Hence conditional on ct , the dynamics are linear. The switching code sequence {ct }t=1:T can be analyzed as a symbolic description of how the model partitions state space (e.g., regime occupancy, switching rate, regime-specific linearization and stability). Learned linear initialization embedding. When initializing an AL-RNN rollout from the first observation x1 ∈ RN (rather than from a full latent initial state z1 ∈ RM ), we embed x1 into the full state space using a learned linear map E ∈ RM ×N , z1 ← Ex1 ,

z1 ← z1 + B T (x1 − Bz1 ),

(A11)

so that the observed coordinates match the data (Bz1 = x1 ) while the unobserved coordinates are provided by the initializer. The initializer E is trained jointly with the AL-RNN parameters under the same training objective (i.e., it is a linear map that produces an initial latent state from the first observation). A0.6

Why objective geometry matters in DSR

For the goals of DSR, there are multiple reasonable learning criteria (teacher-forced prediction losses, likelihoodbased objectives with latent states, and hybrid criteria), and they can disagree strongly in which parameter directions they treat as sensitive. Curvature provides a local summary of this geometry, quantifying “which directions matter”. In particular, under chaos, small changes in parameters lead to large changes in trajectories,

so the curvature of the objective can have a strong influence on which local minima are found and how well they match the true system’s invariant structure. In switching latent-variable models, Louis’ identity implies that latent ambiguity reduces marginal-likelihood curvature through a missing-information correction term. This is the mechanism behind our “geometry mismatch” narrative: training-time forcing can induce sharp curvature (overconfident geometry) when the marginal likelihood, after properly accounting for switching uncertainty, is comparatively flat. QoI-aware objectives. Related approaches in data-driven surrogate modeling emphasize matching long-horizon statistics (e.g., invariant measures and distributional summaries) via losses defined on long rollouts or on learned operators acting on observables (Schiff et al., 2024; Cheng et al., 2025). These objectives are complementary to the local-geometry perspective used here: we focus on how training-time forcing and latent switching uncertainty shape the curvature of the fitted criterion, and how this geometry can disagree with ambiguity-aware likelihood curvature even when the resulting rollouts appear qualitatively plausible.

Methodological details

A1.1

Louis missing-information mechanism for two-state gates

This section states the two-state Louis mechanism used for intuition in the main paper, based on Louis’ identity (Louis, 1982; Dempster et al., 1977). It also specifies the controlled-ambiguity toy model used to visualize the mechanism. For completeness, we restate Louis’ identity in standard latent-variable notation. Let z denote latent variables and define the complete-data log density ℓcomp (θ) := log pθ (x, z). Then the observed information for the marginal likelihood pθ (x) satisfies I obs (θ) = −∇2θ log pθ (x) = E − ∇2θ ℓcomp (θ) | x − Cov ∇θ ℓcomp (θ) | x , (A12) where the conditional moments are taken under the posterior pθ (z | x). Proposition A1 (Louis missing-information term under state-dependent gating). Let C be a latent gate variable and let x denote the observed quantity. For brevity, we write pθ (x, c) := pθ (x, C=c),

pθ (x | c) := pθ (x | C=c),

pθ (c | x) := Pθ (C=c | x).

(A13)

Assume the joint model can be written as pθ (x, c) = π(c | x) pθ (x | c),

c ∈ {0, 1},

(A14)

where π(c | x) is a (possibly x-dependent) gating term that is independent of θ for the parameter block of interest. Define, for each c ∈ {0, 1}, sc (θ) := ∇θ log pθ (x | c), Let p := pθ (1 | x) =

I c (θ) := −∇2θ log pθ (x | c).

(A15)

π(1 | x) pθ (x | 1) . π(0 | x) pθ (x | 0) + π(1 | x) pθ (x | 1)

(A16)

Then the observed information for the marginal likelihood pθ (x) satisfies I obs (θ) = p I 1 (θ) + (1 − p) I 0 (θ) − I miss (θ).

(A17)

T I miss (θ) := p(1 − p) s1 (θ) − s0 (θ) s1 (θ) − s0 (θ) ⪰ 0.

(A18)

In particular, the missing-information term I miss (θ) is large when (i) the posterior is ambiguous (p(1 − p) is maximized at p = 1/2) and (ii) the two explanations disagree in score (large ∥s1 (θ) − s0 (θ)∥). Remark A1 (When gating depends on θ). If the gating distribution depends on θ, then the complete-data score for (x, C) includes an additional term ∇θ log πθ (C | x), so both terms in Louis’ identity acquire extra contributions. In that case, the simple rank-one form in Eq. (A17) need not hold exactly, but the qualitative mechanism remains: greater posterior ambiguity and larger between-component score gaps reduce observed information.

A1.1.1

Proof

Proof of Proposition A1. Because π(c | x) is θ-independent, we have log pθ (x, c) = log π(c | x) + log pθ (x | c).

(A19)

Therefore the complete-data score and complete-data negative Hessian for the parameter block θ are ∇θ log pθ (x, c) = sc (θ),

−∇2θ log pθ (x, c) = I c (θ).

(A20)

Louis’ identity for this two-point latent gate can be written as X I obs (θ) = pθ (c | x) I c (θ) c∈{0,1}

−

pθ (c | x) sc (θ)sc (θ) −

c∈{0,1}

! X T pθ (c | x) sc (θ) pθ (c | x) sc (θ) .

c∈{0,1}

(A21)

c∈{0,1}

Writing p := pθ (1 | x) and simplifying yields T I obs (θ) = p I 1 (θ) + (1 − p) I 0 (θ) − p(1 − p) s1 (θ) − s0 (θ) s1 (θ) − s0 (θ) ,

(A22)

which is Eq. (A17) with the missing-information term shown in Eq. (A18). A1.1.2

Toy model settings (controlled switching ambiguity)

The mechanism illustrated in the main text (Fig. 1a) was computed using a two-regime, scalar switching AR(1) toy model with a probit gate. Let xt ∈ R denote the (observed) scalar state and let ct ∈ {0, 1} denote the latent regime indicator at time t. We generate trajectories from ct | xt ∼ Bernoulli(Φ(xt /σg )) , xt+1 | xt , ct ∼ N act xt , σ 2 , (A23) where Φ is the standard normal CDF and σg controls gate stochasticity (and therefore posterior switching ambiguity). We fix a0 = 0.90, a1 = 0.60, σ = 0.15, and simulate length-T trajectories with T = 600 for 20 random seeds. All time indices in this toy model use t = 1, . . . , T − 1. To vary ambiguity, we sweep σg over 25 logarithmically spaced values between 0.03 and 0.9. For each simulated pair (xt , xt+1 ), the gate posterior is available in closed form: letting πt := Φ(xt /σg ), we have pt := P(ct = 1 | xt , xt+1 ) =

πt N (xt+1 ; a1 xt , σ 2 ) . (1 − πt ) N (xt+1 ; a0 xt , σ 2 ) + πt N (xt+1 ; a1 xt , σ 2 )

(A24)

The switching-uncertainty summary reported in the main paper (Fig. 1a) is the time-average of the Bernoulli entropy (in bits), T −1 1 X − pt log2 pt − (1 − pt ) log2 (1 − pt ) . (A25) H(pt ) = T − 1 t=1 For the same simulated trajectories we compute the toy-model observed information I obs for the drift-parameter block θ = (a0 , a1 ) by applying the two-state Louis decomposition in Eq. (A17) at each time step using the closed-form posteriors pt , and then summing over t. A1.2

Observed information via RBPF and Louis’ identity

This section describes how we estimate observed information for the probabilistic AL-RNN (PAL-RNN), which augments the AL-RNN to a switching linear-Gaussian state-space model (SSM) with latent gate variables (Durbin and Koopman, 2001). The key point is that, conditional on a gate path c1:T −1 , the continuous state admits closed-form Kalman filtering/smoothing and closed-form complete-data score and information contributions. Thus, we approximate the required smoothing distribution with a Rao–Blackwellized particle filter/smoother (RBPF) over gate paths and then apply Louis’ identity to estimate I obs . See e.g. Doucet et al. (2000a) for RBPF background.

Figure 2: Directed graphical model of the probabilistic AL-RNN (PAL-RNN), shown as a conditional Bayesian network over (z1:T , c1:T , x1:T ) given shared drift parameters θ = (A, W , h). The continuous latent state zt emits observations xt , and the discrete switching code ct ∈ {0, 1}P is drawn from a state-dependent probit model p(ct | zt ) and indexes the transition p(zt+1 | zt , ct ; θ). Dashed arrows indicate that the same parameters θ are shared across time. Noise (hyper-)parameters (e.g., Q, R, σg ) are treated as fixed and omitted for visual clarity. A1.2.1

PAL-RNN

The PAL-RNN is a probabilistic augmentation of the AL-RNN (Eq. (A10)) that retains the same switching structure. However, it treats the switching code as a latent random variable via a probit gate model and adds Gaussian process/observation noise, so that conditional on a gate path it becomes a linear-Gaussian SSM. Figure 2 shows the corresponding directed graphical model. We use the same notation as in Eq. (A10). In particular, ct ∈ {0, 1}P encodes which of the last P coordinates are gated (active) and D(ct ) is the corresponding diagonal gate matrix. For a length-T segment, gates index transitions for t = 1, . . . , T − 1. The conditional dynamics and observation model are zt+1 = A + W D(ct ) zt + h + εt , εt ∼ N (0, Q), xt = Bzt + ηt ,

ηt ∼ N (0, R).

(A26)

We use the linear observation model in Eq. (A3) with the identity projection B = [IN 0]. In our experiments, we 2 2 use isotropic noise Q = σproc IM and R = σobs IN . The switching law is state dependent but has no explicit dependence on the drift-parameter block of interest θ = (A, W , h): zt,M −P +j P(ct,j = 1 | zt ) = Φ , j = 1, . . . , P, (A27) σg where Φ is the standard normal CDF and σg is a fixed gate-noise scale. For the curvature diagnostics in this section, we hold (σproc , σobs , σg ) fixed as hyperparameters (Table 1). In the particle-SAEM experiments later in the supplement (Section A2.3), we update (σproc , σobs ) in the calib and full SAEM configurations, while keeping σg fixed. We also treat the learned linear initialization embedding E (Eq. (A11)) as fixed when computing curvature. Under these assumptions, the switching enters Louis’ identity only through posterior uncertainty over gate paths, matching the setting of Proposition A1. Relation to (r)SLDS. The PAL-RNN is closely related to switching linear dynamical systems (SLDS) and recurrent SLDS variants used for hybrid time-series modeling (Ghahramani and Hinton, 2000; Murphy, 1998; Linderman et al., 2017). A key difference is that we do not introduce a separate Markov transition model over discrete regimes: instead, the latent gate bits ct ∈ {0, 1}P are state dependent and (given zt ) conditionally independent across time, reflecting the AL-RNN’s per-unit gating structure. Moreover, SLDS/rSLDS typically use a categorical regime variable with regime-specific linear dynamics, whereas PAL-RNN uses a factorial code that deterministically masks the same shared drift parameters via D(ct ). A1.2.2

RBPF filtering, approximate smoothing, and Louis identity

The RBPF samples gate paths and Rao–Blackwellizes the continuous state using a Kalman filter, yielding a (particle) mixture of Kalman filters (Chen and Liu, 2000). At each time t, particle i carries a sampled gate prefix

(i)

c1:t−1 and a Gaussian filtering belief (i)

(i)

p(zt | x1:t , c1:t−1 ) = N (mt , Pt ).

(A28)

The RBPF represents the gate-prefix posterior p(c1:t−1 | x1:t ) with finitely many particles (Monte Carlo approximation). Conditional on a fixed gate prefix, the Kalman update is exact. Gate bits are proposed with a factorized integrated probit rule, i.e. we analytically average the probit link under the particle’s Gaussian belief. Concretely, if a gated coordinate has marginal z ∼ N (µ, ν) under the particle’s current belief, then  P(c=1 | µ, ν) = Φ q

 µ ν + σg2

.

(A29)

(i)

Given sampled ct , the conditional dynamics are linear-Gaussian and the particle performs a Kalman pre(i) (i) dict/update to obtain (mt+1 , Pt+1 ), with importance weights updated by the one-step predictive likelihood. Resampling is triggered when the effective sample size (ESS) falls below a fraction of Np . With normalized (i) (i) weights wt ∝ exp(log wt ), the ESS is 1 ESSt := PNp (i) , (A30) 2 i=1 (wt ) and when ESSt < τESS Np , where τESS ∈ (0, 1) is the resampling threshold (here τESS =0.5), we resample particle indices. This is done via i.i.d. sampling from the categorical distribution with probabilities wt using multinomial resampling with replacement. Afterwards, we reset weights to 1/Np , and continue. The RBPF used for the experiments employs a bootstrap-style proposal: gate bits are sampled from the integrated gate prior described above, and the weight update uses the Kalman log-likelihood increment log p(xt+1 | (i) x1:t , c1:t ; θ). To approximate Louis’ identity on a length-T segment, we generate S approximate smoothed (c, z) trajectories from the RBPF by tracing ancestors and sampling backward in time under the conditional Gaussian dynamics. For each sampled trajectory, we compute per-step complete-data score contributions st (θ) and complete-data information contributions I comp,t (θ) for the drift-parameter block θ = (A, W , h). The observed information estimate then is given by "T −1 # ! T −1 X X b d b obs = E I I comp,t − Cov st t=1

b comp ] − Cov d = E[I

t=1 T −1 X

! st

(A31)

t=1

b and Cov d denote the empirical mean and covariance over the S sampled smoothed trajectories. Hence, Here E this is Louis’ identity with Monte Carlo estimates of the complete-data expectation and the conditional score covariance under the RBPF smoothing distribution. Pseudo-code for the estimator is given in Algorithm A1. A1.2.3

Filtering switching code entropy and ESS

We summarize switching uncertainty using the time-averaged RBPF filtering (online) switching code entropy Hc (in bits), computed from the RBPF filtering distribution at each transition time t = 1, . . . , T − 1. We use filtering rather than smoothing entropy because it is the most directly aligned with ITF, which is a forward (i) (online) intervention. Let wt denote the normalized particle weights after the weight update at transition t (i) (i.e., after incorporating xt+1 via the one-step predictive likelihood) and let ct ∈ {0, 1}P be the sampled gate bits. Equivalently, the induced full-code distribution p̂t (c) defined below estimates the online gate uncertainty at transition t: p̂t (c) ≈ pθ (ct = c | x1:t+1 ).

(A32)

Algorithm A1 RBPF and Louis estimator of observed information 1: Input: observations x1:T ; parameters θ = (A, W , h); hyperparameters (σproc , σobs , σg ); number of particles Np ; number of smoothing trajectories S. (i) (i) (i) 2: Init: for each particle i, set w1 ← 1/Np and initialize the filter (m1 , P1 ). 3: for t = 1, . . . , T − 1 do 4: for i = 1, . . . , Np do (i) (i) (i) 5: Sample gate bits ct from the factorized integrated-probit proposal under N (mt , Pt ). (i) (i) (i) 6: Run Kalman predict/update to obtain (mt+1 , Pt+1 ) given ct . (i) (i) (i) 7: Update weight: wt+1 ∝ wt p(xt+1 | x1:t , c1:t ; θ). 8: end for 9: Normalize weights. 10: If ESS < τESS Np (default τESS = 0.5), resample and reset weights; store ancestry and sampled codes. 11: end for 12: for s = 1, . . . , S do (s) 13: Sample a gate path c1:T −1 via ancestry tracing / backward sampling. (s)

(s)

Sample a continuous trajectory z1:T from the conditional Gaussian smoother (given x1:T and c1:T −1 ). P (s) P (s) 15: Accumulate complete-data score t st and complete-data information t I comp,t for θ = (A, W , h). 16: end for b comp ] − Cov( d P st ) and diagnostics (ESS, entropies). b obs = E[I 17: Return I t 14:

Full-code entropy. Since the size of the full codebook is manageable in our experiments, we can form an empirical filtering posterior over full codes c ∈ C, where C := {0, 1}P , via p̂t (c) :=

Np X

(i)

1{c(i) t = c},

i=1

Ht,code := −

p̂t (c) log p̂t (c).

(A33)

c∈C

When 2P is too large to tabulate the full code distribution, one can instead use a bit-marginal entropy sum as a cheaper proxy. Time-averaged entropy in bits. The entropy values reported in the curvature-gap results (Fig. 1b) are the time average of the full-code entropy, converted to bits: Hc :=

T −1 X 1 Ht,code . (T − 1) log 2 t=1

(A34)

Missing-information ratio (MIR). To summarize the magnitude of the Louis missing-information correction relative to the complete-data curvature, we report the missing-information ratio b obs ) tr(I b tr(E[I comp | x1:T ]) b comp | x1:T ] − I b obs ) tr(E[I

MIR := 1 − =

b comp | x1:T ]) tr(E[I

(A35)

b comp | x1:T ] is the Monte Carlo estimate of the complete-data information term from the same smoothing where E[I b obs . Values closer to 1 indicate larger curvature reductions due to latent uncertainty. samples used for I A1.2.4

Settings used for curvature-gap experiments

In the curvature-gap experiment (Fig. 1b), the RBPF/Louis diagnostics are computed on fixed-length segments. The hyperparameter settings are summarized in Table 1.

Table 1: Settings used for the curvature-gap analysis in Fig. 1b. Setting AL-RNN training settings Dataset Training trajectory length Ttrain Noise regimes Forcing interval τ Latent dimension M Number of gated units P Batch size B BPTT sequence length LBPTT Epochs Batches per epoch RBPF/Louis settings Sequence length for curvature diagnostics T Particles Np Smoothing samples S Resampling threshold τESS Gate noise σg

A1.3

Value Lorenz-63 8 · 104 σproc ∈ {0.1, 0.3, 0.5}, σobs ∈ {0.1, 0.3, 0.5} {4, 8, 16, 32, 64} 30 10 16 200 2000 50 200 64 8 0.5 0.1

ITF-aligned curvature proxy

This section gives the closed-form computation of the ITF-aligned generalized Gauss–Newton/Fisher curvature proxy used throughout the experiments. Under identity teacher forcing (ITF), forcing is applied at times Tτ = {t : t ≡ 0 (mod τ ), t > 0}, as described in Appendix A0 (Eq. (A4)). We compute the per-step generalized Gauss–Newton/Fisher matrix aligned with this intervention-based training objective, denoted I ITF , by differentiating the forced rollout and accumulating the resulting Jacobians. We use tr(I ITF ) as a curvature proxy for ITF training, and as a comparator to ambiguity-aware observed information (Section A1.2). Gaussian pseudo-likelihood aligned with ITF. Let {xt }Tt=1 be an observed trajectory with xt ∈ RN and let x̂t := Bzt be the one-step prediction of the observed coordinates. The ITF loss used for training is the MSE along the forced rollout, T −1 1 X 2 Bzt+1 − xt+1 2 . (A36) LITF (θ; x1:T ) = T − 1 t=1 For curvature diagnostics, we interpret this as a Gaussian negative log pseudo-likelihood with an isotropic weighting 2 Λ = σobs IN , so that the generalized Gauss–Newton matrix equals the Fisher matrix for this pseudo-likelihood. ITF intervention and sensitivity reset. Let zt ∈ RM be the latent state and let B = [IN 0] be the projection from Eq. (A3). At forcing times, ITF overwrites the observed coordinates, z̃t = zt + B T (xt − Bzt ) = M zt + B T x t ,

M := I − B T B.

(A37)

Since xt does not depend on parameters, this implies the sensitivity boundary condition ∂ z̃t ∂zt =M . ∂θ ∂θ

(A38)

Thus, at forcing times, the sensitivity ∂zt /∂θ is reset by multiplication with M , which zeros out the first N rows (i.e. the observed coordinates) and leaves the remaining M − N rows unchanged. Piecewise-linear derivatives along the forced rollout. The forced rollout evolves according to the AL-RNN transition map evaluated at the forced or unforced state, ( z̃t , t ∈ Tτ , ∗ zt+1 = Fθ (z̄t ) = Az̄t + W ϕ (z̄t ) + h, z̄t := (A39) zt , else.

Let ct ∈ {0, 1}P denote the induced switching code at time t along this rollout, i.e. ct,j = 1{z̄t,M −P +j > 0}. With D(ct ) as in Eq. (A10), the Jacobian of the one-step map with respect to the state is Jt :=

∂Fθ (z̄t ) = A + W D(ct ), ∂z

(A40)

where we adopted the subgradient convention of setting derivatives at the ReLU-kink to zero. Closed-form sensitivity recursion for ∂zt /∂θ. We evaluate curvature for the drift-parameter block θ = (A, W , h) and treat the initial condition as fixed. (In our implementation, z1 is obtained from a learned linear embedding of the first observation, which we keep fixed for curvature diagnostics, see Eq. (A11)). For the sensitivity recursion, we parameterize A by its diagonal entries a := diag(A) ∈ RM and vectorize   a θ := vec(W ) ∈ Rp , p = M + M 2 + M, (A41) h using row-major vectorization for vec(W ). Define the sensitivity matrix St := ∂zt /∂θ ∈ RM ×p and initialize S1 = 0. For each step t = 1, . . . , T − 1 we compute ( M St , t ∈ T τ , St+1 = Jt S̄t + Vt , S̄t := (A42) St , else, where Vt := ∂Fθ /∂θ has the following closed form at state z̄t :

(a)

Thus Vt = [Vt

(W )

(h)

(a)

(a block)

(W block)

(h block)

= diag(z̄t ) ∈ RM ×M ,

(W ) (h)

= IM ⊗ ϕ∗ (z̄t )T ∈ RM ×M ,

= IM ∈ RM ×M .

(A43)

Information accumulation and per-step normalization. The Jacobian of the one-step prediction x̂t+1 = Bzt+1 with respect to parameters is BSt+1 ∈ RN ×p . For the isotropic pseudo-likelihood weighting Λ = σ 2 IN , the ITF-aligned Fisher is I ITF =

T −1 1 X T (BSt+1 ) Λ−1 (BSt+1 ) T − 1 t=1

T −1 X 1 T (BSt+1 ) (BSt+1 ) , (T − 1) σ 2 t=1

(A44)

which is positive semidefinite by construction (up to floating-point roundoff). For notational simplicity, we write I ITF for this per-step-normalized ITF-aligned Fisher in the main text and the remainder of this supplement.

Datasets and experimental settings

A2.1

Lorenz-63 dataset

All AL-RNN checkpoints analyzed in this work (both for curvature-gap diagnostics and as initializations in the SAEM experiments) are trained on trajectories generated from the Lorenz-63 system (Lorenz, 1963) across a sweep of observation- and process-noise regimes. The continuous-time Lorenz-63 dynamics are dz1 = σ (z2 − z1 ), dt

dz2 = z1 (ρ − z3 ) − z2 , dt

dz3 = z1 z2 − βz3 , dt

(A45)

with the standard chaotic parameter setting σ = 10, ρ = 28, and β = 8/3. We denote the corresponding discrete-time latent state by zt ∈ R3 (with components corresponding to (z1 , z2 , z3 )). We discretize with step size

∆t = 0.01. In the deterministic regime (σproc = 0), we integrate the ODE using the fourth order Runge-Kutta method (RK4) (Kutta, 1901). In the stochastic regime (σproc > 0), we use an additive-noise SDE discretization implemented as an RK4 drift step plus Euler–Maruyama diffusion, √ zt+1 = RK4∆t (zt ) + σproc ∆t diag(sref )εt , εt ∼ N (0, I3 ), (A46) where zt denotes the latent Lorenz state and sref is the per-coordinate standard deviation of the deterministic reference trajectory. Thus σproc is reported in standardized-coordinate units, although the perturbation is applied before standardization. Observation noise is added in raw coordinates as xt = zt + σobs diag(strain )ηt ,

ηt ∼ N (0, I3 ),

(A47)

where strain is the per-coordinate standard deviation used by the training standardizer. After standardization, σobs is therefore also on a standardized-coordinate noise scale. A2.2

AL-RNN training settings

All ITF-trained model checkpoints analyzed in this work (both for curvature-gap diagnostics and as initializations in the SAEM experiments) are trained on the Lorenz-63 time series generated by the procedure above. For training and evaluation, we standardize observations per dimension to zero mean and unit variance using the mean and standard deviation computed on the training trajectory. We denote the interleaving forcing interval by τ . The loss minimized is MSE on the next-step prediction of the observed coordinates along the interleaved rollout, as described in Appendix A0 (Eq. (A36)). Parameter initialization. For each random initialization seed, we initialize the drift parameters by setting h = 0, drawing entries of W i.i.d. from N (0, 0.01), and taking the diagonal entries of A from the diagonal of a random symmetric positive definite matrix normalized to have largest-magnitude eigenvalue √ √1. The learned initialization embedding matrix E ∈ RM ×N is initialized by drawing i.i.d. from U (−1/ N , 1/ N ). BPTT training details. Training is performed by BPTT on mini-batches of contiguous subsequences sampled from the standardized training trajectory x1:Ttrain . Concretely, each training update uses a batch of B sequences (here B = 16) of length LBPTT (here LBPTT = 200) drawn uniformly at random (with replacement) from the training trajectory, and we run a fixed number of batches per epoch (here 50). We optimize the ITF loss with the RAdam optimizer (Liu et al., 2020) using default momentum parameters β1 = 0.9, β2 = 0.999 and ϵ = 10−8 , and an exponential learning rate decay schedule from 10−3 to 10−5 over the course of training. The training settings and the subset of regimes used in the curvature-gap figure are summarized in Table 1. A2.3

Particle-SAEM fine-tuning for windowed evidence

This section documents the particle-SAEM procedure used for the analysis of windowed evidence and dynamical QoIs in Fig. 1c. The results in Fig. 1c were aggregated over 12 ITF-pretrained AL-RNN initializations, one per training noise, each evaluated at five window lengths. The settings used in the experiments are summarized in Table 2. A2.3.1

Dataset used in SAEM experiments

For the SAEM fine-tuning experiments we simulate deterministic Lorenz-63 data (σproc = 0) for T = 5000 steps, add i.i.d. observation noise with standard deviation σobs = 0.085, and stack nseq = 4 independent sequences. We remove potential transients by discarding the first 1000 time steps and standardize each observed coordinate to zero mean and unit variance using statistics pooled over all sequences and time indices. A2.3.2

Windowed evidence objective

Given a sequence x1:T (i.e. one of the nseq Lorenz-63 sequences), we partition it into non-overlapping windows of length L ∈ N with stride L. Concretely, for window starts tk := 1 + (k − 1)L,

k = 1, . . . , nwin ,

with tk + L ≤ T,

(A48)

we define the kth window as the observation block xtk +1:tk +L together with its left boundary condition xtk . For a PAL-RNN parameter setting θ (drift parameters plus any noise parameters being optimized in the chosen configuration), the windowed evidence objective is the conditional log-evidence log pθ xt+1:t+L | xt , (A49) where the conditioning on xt fixes the window boundary condition through the learned initialization described in Eq. (A11). In practice we approximate this quantity with the RBPF estimate of the particle filter normalizing constant, bθ (t, L), log pθ xt+1:t+L | xt ≈ log Z (A50) (i)

bθ (t, L) is obtained from the same RBPF forward pass used for smoothing (Section A1.2). Let w where Z t+s−1 denote the normalized particle weights before incorporating xt+s , and let (i) (i) ℓt+s := pθ xt+s | xt , xt+1:t+s−1 , ct:t+s−1 ,

(A51) (i)

denote the one-step predictive likelihood under particle i (available in closed form because, conditional on ct:t+s−1 , the model is linear-Gaussian). These are the same increments used in the RBPF weight update. The corresponding particle-filter normalizing-constant estimator over the window is Zbθ (t, L) :=

L Y

pbθ xt+s | xt , xt+1:t+s−1 ,

s=1 N

p X (i) (i) pbθ xt+s | xt , xt+1:t+s−1 := wt+s−1 ℓt+s .

(A52)

i=1

This sequential Monte Carlo estimator targets the desired conditional marginal likelihood and becomes accurate as particle budgets increase (Doucet et al., 2000b). We therefore use log Zbθ (t, L) as an approximation to the windowed conditional log-evidence. bθ (t, L) on held-out windows, averaged over windows and normalized per step and per In Fig. 1c we report log Z observed dimension, n win X 1 log Zbθ (tk , L), (A53) Evidence(θ) := nwin L N k=1

computed on a fixed set of nwin held-out windows (Table 2), where N is the observation dimension. A2.3.3

SAEM updates and configurations

For each initialization and each window length L, we consider three parameter-update configurations. We write the full parameter vector as θ := (θdrift , σproc , σobs ), θdrift := (A, W , h),

(A54)

and keep the gate-noise scale fixed at σg = 0.7 throughout. The three configurations differ only in which blocks of θ are updated: • baseline: no updates, the objective is evaluated at initialization. • calib: update only (σproc , σobs ) and keep θdrift fixed. • full SAEM: update both θdrift and (σproc , σobs ). All configurations share the same PAL-RNN switching structure with P gated coordinates (Section A1.2); only the updated parameter blocks differ. The underlying AL-RNN settings and the additional SAEM/RBPF hyperparameters are listed in Table 2.

Algorithm A2 Particle-SAEM update on windowed evidence (j)

seq 1: Input: sequences {x1:T }j=1 ; window length L, initial parameters θ (1) , configuration (baseline / calib /

full SAEM), fixed gate noise σg .

2: Hyperparameters: SAEM iterations R, windows per iteration B, RBPF particles Np , smoothing samples

S, resampling threshold τESS , ridge regularization weight λR , blend factor αM .

3: for r = 1, . . . , R do 4: Sample a mini-batch of B windows (sequence index j and start time tk ) from the training windows. 5: for each sampled window (j, tk ) do (j)

(j)

Run RBPF filtering on xtk :tk +L (conditioning on xtk ) and generate S approximate smoothed trajectories by backward sampling. From the smoothing samples, form Monte Carlo estimates of the complete-data sufficient statistics needed 7: for the M-step (restricted to the parameter blocks updated by the chosen configuration). 8: end for 9: Aggregate sufficient statistics over the mini-batch (empirical mean over windows and smoothing samples). 10: Compute a ridge-regularized M-step solution θ̃ (r) using ridge parameter λR . 11: Apply the damped update θ (r+1) ← (1 − αM )θ (r) + αM θ̃ (r) , respecting the chosen configuration. 12: end for 13: Return: final parameters θ (R+1) . 6:

train train Note that the training noise scales σproc and σobs refer to the data-generation regime used when fitting the initialization checkpoint, and are distinct from the SAEM/RBPF noise parameters (σproc , σobs ) optimized during fine-tuning. For numerical stability, whenever (σproc , σobs ) are updated we impose a small positive floor 2 σmin = 10−4 in the M-step (equivalently, estimated residual variances are clamped below by σmin ).

For each initialization we run SAEM for all window sizes L ∈ {16, 32, 64, 128, 200}, yielding 60 initialization-window settings per configuration. We optimize the windowed evidence objective via stochastic approximation expectation maximization (SAEM), using a Monte Carlo E-step based on samples from the RBPF smoothing distribution (Delyon et al., 1999). Operationally, we implement SAEM as damped EM updates on random mini-batches of windows. At each iteration we (i) approximate the E-step expectations under the window-wise smoothing distribution via RBPF backward sampling, (ii) compute the corresponding complete-data sufficient statistics (restricted to the parameter blocks updated in the chosen configuration), (iii) compute a ridge-regularized M-step solution to obtain a provisional update, and (iv) blend this provisional update into the running parameters. We use a fixed blend factor αM (“blend step” in Table 2) and update rule θ (r+1) := (1 − αM ) θ (r) + αM θ̃ (r) ,

(A55)

where θ̃ (r) denotes the ridge-regularized M-step solution for the chosen parameter blocks based on the Monte Carlo E-step at iteration r. Pseudo-code is given in Algorithm A2.

Metric definitions and computation details

This section defines the scalar metrics reported in the main-text figures and fixes their normalizations. Unless otherwise stated, all quantities are computed on standardized observation coordinates. A3.1

Toy-model diagnostics (Fig. 1a)

b obs The toy model in Fig. 1a reports two curvature/ambiguity summaries derived from Louis’ identity. Let I b denote the observed information estimate for the toy-model parameter block of interest, and let E[I comp | x] denote the corresponding complete-data information term (both computed for the same length-T trajectory). In this toy model the gate posterior pt is available in closed form (Section A1.1), so these quantities are computed by directly evaluating Eq. (A17) (with θ = (a0 , a1 )) and aggregating over time, rather than via RBPF.

Table 2: Settings used in the SAEM fine-tuning experiments in Fig. 1c. Quantity Value AL-RNN initialization (ITF pretraining) Number of models 12 AL-RNN initializations Dataset Lorenz-63 Pretraining objective ITF next-step MSE on standardized observations (Eq. (A36)) Training trajectory length Ttrain 8 · 104 Test trajectory length Ttest 2 · 104 train train Noise regimes σproc ∈ {0.1, 0.3, 0.5}, σobs ∈ {0, 0.1, 0.3, 0.5} Forcing interval τ 16 Latent dimension M 30 Number of gated units P 10 Batch size B 16 BPTT sequence length LBPTT 200 Epochs 2000 Batches per epoch 50 RBPF/Louis settings Sequence length for curvature diagnostics T 200 Particles Np 256 Smoothing samples S 8 Resampling threshold τESS 0.5 Gate noise σg 0.7 SAEM settings SAEM iterations R 8 Windows per iteration 80 Ridge λR 10−2 Blend step αM 0.25 Held-out windows nwin 120 Window sizes L {16, 32, 64, 128, 200}

Missing-information ratio (MIR).

We report the missing-information ratio MIR := 1 −

b obs ) tr(I b comp | x] tr E[I

b comp | x] − I b obs tr E[I , = b comp | x] tr E[I

(A56)

which takes values in [0, 1] and increases with posterior latent uncertainty (Section A1.2, “Missing-information ratio (MIR)”). b obs ). In the toy-model Observed-curvature proxy. As a scalar curvature proxy we also report log10 tr(I experiment, this quantity is used only for relative comparison across gate-noise settings σg at fixed trajectory length T . A3.2

Curvature gap / sharpness mismatch (Fig.1b)

For a trained AL-RNN checkpoint, with (i) the ITF-aligned curvature proxy I ITF (Section A1.3) and (ii) an b obs computed on a trajectory of length T (Section A1.2), we define RBPF/Louis observed information estimate I the per-step observed information trace as b obs,step ) := tr(I

b obs ) tr(I . T

(A57)

The curvature-gap statistic reported in Fig. 1b is gQ := log10

tr(I ITF ) b obs,step ) tr(I

! .

(A58)

Both numerator and denominator are interpreted on a per-step scale (in our implementation I ITF is stored with the per-step normalization from Eq. (A36)). A3.3

Windowed evidence (held-out; RBPF estimate) (Fig. 1c)

For a fixed window length L, the windowed conditional log-evidence is approximated by the RBPF normalizingconstant estimate log Zbθ (t, L) (Section A2.3). The quantity shown on the y-axis in Fig. 1c is the held-out b over a fixed set of held-out windowed evidence score Evidence(θ) (Eq. (A53)), defined as the mean of log Z windows, normalized per step and per observed dimension by L N (Section A2.3). We emphasize that this is a windowed and conditional evidence estimate, not the marginal likelihood of an entire long free-running trajectory. A3.4

Dynamical QoI metrics (Fig. 1c)

For the QoI evaluations in Fig. 1c we use long free-running rollouts of length T = 104 time steps (discarding a fixed burn-in of 1000 steps), and compute both Dstsp and Lyapunov-based quantities on the resulting generated and reference trajectories. T true T State-space divergence Dstsp . Given generated and reference trajectories {xgen }t=1 , we t }t=1 and {xt construct a joint histogram over a fixed bounding box defined by the minimum/maximum of the reference trajectory, using nbins bins per dimension (nbins = 30 throughout this work). Let p̂ and q̂ be the resulting discrete distributions after Laplace smoothing with αsmooth = 10−5 . Then

Dstsp := KL(p̂ ∥ q̂) =

X b

p̂(b) log

p̂(b) . q̂(b)

(A59)

In the SAEM rollouts, the model is simulated in the standardized coordinates used by the objective, and we then apply the inverse standardization map before computing Dstsp so that the histogram comparison is carried out in the physical (raw) Lorenz coordinates. Lyapunov exponent and error. We evaluate the largest LE λmax of the learned discrete-time dynamics using a standard method based on QR re-orthonormalization (Benettin et al., 1980), applied to a sequence of closed-form one-step model Jacobians along a rollout. Rollouts are initialized from the first observation x1 ∈ RN using the learned linear embedding E (Eq. (A11)). Along the rollout, let ct ∈ {0, 1}P denote the induced switching code at time t (Eq. (A10)) and let D(ct ) be the corresponding diagonal gate matrix. The local Jacobian Jt of the one-step map is given by Eq. (A40). In the SAEM rollout evaluation we use hard-gated deterministic rollouts for the results shown in Fig. 1c, i.e. we evolve the state with the deterministic AL-RNN transition and set the gate bits by the ReLU sign pattern at each time step. For additional analyses (Section A4), we also consider a stochastic probit-gated rollout that samples Bernoulli gates with probabilities Φ(zt,j /σg ) on the gated coordinates. Given Jacobians {Jt }Tt=1 of the model, the discrete-time Lyapunov spectrum is estimated by QR reorthonormalization (Benettin et al., 1980): T X bdisc = 1 λ log |Rt,kk | , k T t=1

where Jt Qt−1 = Qt Rt .

(A60)

b1 := λ bmax as the largest LE of Continuous-time exponents are then obtained through dividing by ∆t. We report λ the model after discarding an initial 1000 time steps of the rollout. For the Lorenz-63 reference value we use the literature value λref max = 0.9056 (Sparrow, 1982). For σproc > 0 we interpret this as a reference for the underlying noiseless drift dynamics; a noise-matched comparison would require a random-DS/SDE notion of LEs and a separate stochastic reference estimate, which we do not compute here. b1 − λref . We report the signed Lyapunov error λ max

Table 3: Fixed-Q stratification of the entropy–gap association (stochastic checkpoints with σobs > 0). Within each fixed-σproc stratum (n = 60 checkpoints), we report the Spearman rank correlation between RBPF filtering code entropy Hc and curvature gap gQ (Fig. 1b). “Partial r” controls for σobs via rank-residualization. 95% CIs are bootstrap percentiles (2000 resamples). σproc 0.1 0.3 0.5

A3.5

n Hmin Hmax H span Spearman r Spearman p Partial r Spearman 95% CI 60 2.173 4.352 2.180 0.630 2.0 × 10−7 0.016 [0.426, 0.773] 60 4.314 5.929 1.615 0.628 7.8 × 10−8 -0.039 [0.438, 0.751] 60 5.488 6.443 0.955 0.585 9.3 × 10−7 -0.310 [0.410, 0.717]

Switching code diagnostics

In the RBPF/Louis recomputation we report Hc , the time-averaged RBPF filtering switching code entropy (in bits) defined in Section A1.2.

Additional analyses

This section collects additional analyses not included in the main text. We (i) probe robustness of the entropy–gap association under stratifications, (ii) sanity-check RBPF computations with time-resolved diagnostics, (iii) provide qualitative visualizations that help interpret the quantitative QoI metrics and (iv) explore additional curvature diagnostics that complement the trace-based curvature-gap statistic in Fig. 1b. We also briefly discuss the computational cost of the RBPF/Louis computations and the SAEM fine-tuning procedure. A4.1

Fixed-Q stratification of entropy–gap association

To check whether the entropy–gap association reported in the main paper (Fig. 1b; gQ versus Hc at fixed Q) depends strongly on the fixed process-noise level σproc (i.e., fixed-Q strata), we re-analyzed the fixed-Q stochastic checkpoints (σobs > 0) separately for each σproc ∈ {0.1, 0.3, 0.5}. Each stratum contained n = 60 checkpoints. Across strata, the rank association between RBPF filtering switching code entropy Hc and curvature gap gQ remained positive. However, after rank-residualizing both variables against σobs , the partial correlations are near zero or negative. Thus, the positive within-stratum rank association is largely explained by variation in observation-noise level inside each process-noise stratum, rather than by a strong residual model-to-model entropy effect at fixed noise settings. Table 3 summarizes the per-stratum results. Note that Hc is computed from the RBPF filtering distribution (Section A1.2), and therefore reflects an online ambiguity diagnostic (ITF-aligned). Smoothing-based entropies will typically be smaller. A4.2

RBPF diagnostics

The main text focuses on scalar summaries of switching ambiguity and curvature (e.g., Hc and trace-based curvature proxies). As a basic sanity check for the RBPF computations underlying these summaries (Section A1.2), we optionally inspect time series diagnostics that probe particle degeneracy and switching uncertainty. Effective sample size (ESS). We monitor the effective sample size ESSt defined in Eq. (A30). Resampling is triggered when ESSt < τESS Np , where the threshold τESS is reported in Table 2. Resampling events. We record the time indices at which resampling occurs. Frequent resampling is expected in highly informative (low-noise) regimes but can also indicate proposal mismatch. Conversely, persistently low ESS without resampling indicates degeneracy. Time-resolved switching-code entropy. At each transition time, we compute the filtering-distribution switching uncertainty via the full-code entropy Ht,code (Section A1.2). We report this trace in bits, i.e. we plot (bits)

Ht,code := Ht,code / log 2. This provides a time-resolved view of ambiguity that complements the time-averaged summary Hc .

(A61)

Figure 3: RBPF diagnostics on a representative SAEM evaluation window. Panels: (a) normalized ESS (ESSt /Np ) with resampling events (triggered when ESSt < τESS Np ); (b) full-code filtering entropy in bits, (bits) Ht,code := Ht,code / log 2. These diagnostics help confirm that particle degeneracy is controlled under the chosen RBPF budget and that the switching-uncertainty summaries (e.g., Hc ) reflect switching ambiguity rather than particle collapse. Figure 3 illustrates these optional diagnostics on a representative SAEM evaluation window (window length L = 32). All RBPF settings for this run match Table 2. This visualization serves to exemplify that the reported entropy and curvature summaries are not dominated by particle collapse artifacts. A4.3

Matrix-aware mismatch diagnostics (anisotropy and eigenspace misalignment)

The main text focuses on scalar (trace-based) curvature summaries because they are robust and easy to interpret. Figure 4 complements this view by probing whether the mismatch between the ITF-aligned curvature proxy I ITF b obs (Section A1.2) is primarily a (Section A1.3) and the ambiguity-aware Louis observed information estimate I uniform rescaling, or whether it is anisotropic and/or accompanied by eigenspace misalignment. Throughout, we compare matrices on the same segment and on a per-step scale. Let b obs,step := I b obs /T. I

(A62)

Recall that I ITF already includes the per-step normalization from the ITF loss (Section A1.3). To stabilize the Monte Carlo Louis estimator at finite Np , S, we apply ridge stabilization and symmetrize the respective matrices before computing the following diagnostics. Let p be the parameter dimension of the block under study (here, the drift parameters (A, W , h) with p = M 2 + 2M ) and define µ := ε

b obs,step ) + tr(I ITF ) tr(I , 2p

ε = 10−6 .

(A63)

We then set 1 b b (µ) bT I obs,step := 2 (I obs,step + I obs,step ) + µI, (µ) I ITF := 12 I ITF + I T ITF + µI.

(A64)

Log-det gap ∆log det .

We report the ridge-stabilized log-determinant difference in decades, (µ) b (µ) log det I ITF − log det I obs,step ∆log det := . log 10

(A65)

This is a global (volume) mismatch summary that is sensitive to anisotropy, not only to trace scaling. γ Direction-wise curvature ratios via generalized eigenvalues q0.9 . To characterize curvature ratios in a direction-dependent way, we consider the generalized eigenproblem (µ)

(µ)

b I ITF ui = γi I obs,step ui ,

i = 1, . . . , p,

(A66)

−1 (µ) b (µ) I ITF ). Here γi compares the ITF curvature to the ambiguity-aware equivalently {γi }pi=1 = eig((I obs,step ) observed information along the generalized-eigenvector direction ui ; γi > 1 indicates a direction that is sharper under ITF than under the marginal likelihood geometry. We summarize the spectrum by quantiles of log10 γi and define the shorthand qαγ := Quantileα {log10 γi }pi=1 , α ∈ (0, 1). (A67) γ In Fig. 4(b) we report the 90th percentile summary q0.9 , which captures the upper end of the curvature ratio distribution and is sensitive to the presence of directions that are much sharper under ITF than under the γ marginal likelihood geometry. Note that q0.9 > 0 indicates that at least 10% of the generalized eigenvalues γi exceed 1, i.e. there are nontrivial directions along which the ITF curvature is sharper than the ambiguity-aware γ observed information. Conversely, q0.9 < 0 indicates that at least 90% of the γi are below 1, i.e. most directions are sharper under the marginal likelihood geometry than under ITF. (µ) b (µ) Top-k sensitive-subspace overlap ovk . Finally, we compare the leading eigenspaces of I ITF and I obs,step . (µ) (µ) p×k b Let UITF , Uobs ∈ R collect the top-k eigenvectors of I and I , respectively. We define ITF

ovk :=

obs,step

1 T ∥UITF Uobs ∥2F ∈ [0, 1], k

(A68)

the mean squared cosine of principal angles between the two k-dimensional subspaces (since the eigenvectors can T be chosen orthonormal, the singular values of UITF Uobs equal cos θi and ∥ · ∥2F sums cos2 θi ). In Fig. 4(c) we use k = 50 and denote the reported overlap by ov50 . Interpretation (Fig. 4). Across the same checkpoints and fixed-σproc segments used in Fig. 1b, the matrixγ aware summaries show that ∆log det is negative and q0.9 < 0, meaning that at least 90% of generalized-eigenvalue directions have γi < 1. Thus, most directions are not sharper under ITF. The mismatch is instead anisotropic and concentrated: trace-level gaps can coexist with broad direction-wise shrinkage, while the leading-subspace overlap remains low and decreases with switching entropy. We interpret these diagnostics as evidence for subspace-dependent geometry mismatch instead of broad direction-wise overconfidence of ITF. A4.4

Qualitative attractor reconstruction visualization

As a qualitative complement to the scalar QoI metrics in Section A3, we visualize attractor reconstruction for a representative AL-RNN initialization checkpoint and its SAEM-fine-tuned PAL-RNN counterpart (“full SAEM” in Fig. 1c). Figure 5 complements the scalar QoI metrics in Fig. 1c and illustrates that improving windowed evidence need not preserve the qualitative attractor geometry. A4.5

Computational complexity considerations

Let T be the sequence length, M the latent dimension, N the observation dimension, and p the parameter dimension of the drift block (here p = M 2 + 2M for A, W , h). Computing the ITF Fisher proxy in closed form (Section A1.3) requires propagating sensitivities St ∈ RM ×p and accumulating a p × p outer product. In a straightforward implementation, this costs on the order of O T M p + T N p + T p2 = O(T M 4 ) (for fixed N ) (A69)

Figure 4: Matrix-aware mismatch diagnostics comparing ITF curvature I ITF to ambiguity-aware observed b obs on the same segments, stratified by fixed process noise σproc . Points are colored by observationinformation I noise level σobs . (a) Ridge-stabilized log-det gap in decades, ∆log det . (b) Upper-quantile direction-wise ratio, (µ) b (µ) q γ , where γi are generalized eigenvalues of (I )−1 I (Eq. (A66)). (c) Leading-subspace overlap, ov50 , obs,step (µ) (µ) b between the top-50 eigenspaces of I ITF and I obs,step . 0.9

ITF

time, with memory O(M p + p2 ) = O(M 4 ) if the full matrix is stored. Thus the cost is dominated by the number of latent dimensions M of the model. b obs via RBPF/Louis (Section A1.2), with Np particles and S approximate Considering the estimation of I smoothing trajectories, the cost scales roughly as O (Np + S) T KF(M ) , (A70) where KF(M ) denotes the per-step cost of the conditional linear-Gaussian update (O(M 3 ) in the worst case, i.e. for dense covariances).

Limitations and future work

All curvature quantities used here (I ITF and I obs ) are inherently local summaries around a parameter value. Interpreting them as posterior precision (e.g., via Laplace/Bernstein–von Mises-type arguments) requires additional regularity and, at least locally, approximately Gaussian behavior of the objective. In non-convex switching models, the global geometry can differ substantially due to multimodality and flat directions, so curvature is best read as a diagnostic of local sensitivity rather than a complete description of posterior shape. Further, our observed-information and fine-tuning results rely on (i) Monte Carlo RBPF/Louis approximations for the PAL-RNN, (ii) a windowed conditional evidence objective, and (iii) experiments on Lorenz–63 initialized from ITF-trained AL-RNN checkpoints. These choices bound the scope of the conclusions. First, the evidence objective in Fig. 1c is computed on fixed-length windows as log pθ (xt+1:t+L | xt ) (Section A2.3), rather than as a full-sequence marginal likelihood for a single long free-running trajectory. Consequently, improving windowed evidence is not guaranteed to improve long-horizon dynamical QoIs, and can favor parameter changes that improve short-window predictability while altering invariant structure. This misalignment is visible empirically when drift parameters are updated (full SAEM). We view windowing as a pragmatic stabilization of likelihood estimation for chaotic, long-horizon sequences, not as a substitute for QoI-aligned training. Second, the RBPF/Louis observed-information estimates are Monte Carlo approximations whose accuracy depends on particle/smoothing budgets, resampling behavior, and the chosen PAL-RNN augmentation (Section A1.2). In more nonlinear regimes or at longer horizons, these estimators can suffer from particle degeneracy or high variance; stability under increased budgets and alternative proposals is therefore an important practical consideration (Fig. 3). Third, the PAL-RNN uses simplified noise structure (isotropic Gaussian process and observation noise, and fixed gate-noise scale σg ). Extending the analysis to richer noise models (e.g., anisotropic or state-dependent noise), to learning or calibrating σg , and to partial-observation settings would improve realism and help identify which aspects of the geometry mismatch persist.

Figure 5: Qualitative attractor reconstruction visualization for a representative AL-RNN initialization checkpoint and its SAEM-fine-tuned PAL-RNN counterpart (“full SAEM” in Fig. 1c). Rollouts are length T = 5 · 104 time steps and initialized at the first ground-truth state. (a) Attractor reconstruction. Baseline closely matches the ground truth (GT) attractor geometry, while the deterministic rollout of the SAEM-fine-tuned model collapses to a strongly contracted periodic orbit. (b) Zoomed-in view of the trajectories, with GT, baseline, and SAEMfine-tuned rollouts initialized at the same point to highlight sensitive dependence on initial conditions: while the baseline trajectory (necessarily) diverges in phase from the GT reference, it still captures the overall attractor geometry. Fourth, the dynamical QoIs reported for particle-SAEM fine-tuning are evaluated on deterministic, hard-gated rollouts (noise set to zero after thresholding the learned gate probabilities). This isolates drift dynamics for comparability, but it does not evaluate invariants of the full stochastic PAL-RNN implied by the fitted noise parameters. A natural next step is to estimate QoIs directly from sampled free-running PAL-RNN trajectories (and their induced stationary distributions), and to compare drift-only invariants to their distributional/stochastic counterparts. On the inference/geometry side, it would be valuable to develop scalable ambiguity-aware curvature summaries beyond trace proxies (e.g., leading-eigenspace estimates or structured approximations) and to connect them more directly to downstream goals such as uncertainty quantification/calibration, robust fine-tuning, and optimal experimental design. Extending the empirical comparison beyond strict ITF to GTF (e.g., sweeping forcing strength α in Eq. (A8)) would further clarify how intervention policies interpolate between stable training geometry and ambiguity-aware likelihood geometry. Finally, the ITF-aligned Fisher proxy I ITF exploits AL-RNN structure to yield a closed-form, training-aligned curvature summary under a specific intervention policy. Extending such intervention-aligned curvature proxies to other surrogate classes and downstream applications is a promising direction.

Record · ID 141467 · SHA-256 8fa236cc576da79c
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.