Panda: A pretrained forecast model for chaotic dynamics
Authors: Jeffrey B. Lai
Anthony Bao
William Gilpin
Venue: ICLR 2026 Poster (2026)
Keywords: chaos
nonlinear dynamics
forecasting
physics
scientific machine learning
dynamical systems
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We introduce a novel large chaotic systems dataset and use it to pretrain a patch-based foundation model for chaos.

Abstract: Chaotic systems are intrinsically sensitive to small errors, challenging efforts to construct predictive data-driven models of real-world dynamical systems such as fluid flows or neuronal activity. Prior efforts comprise either specialized models trained on individual time series, or foundation models trained on vast time series databases with little underlying dynamical structure. Motivated by dynamical systems theory, we present Panda, Patched Attention for Nonlinear Dynamics. We train Panda on a novel synthetic, extensible dataset of 20,000 chaotic dynamical systems that we discover using an evolutionary algorithm. Trained purely on simulated data, Panda exhibits emergent properties: zero-shot forecasting of unseen chaotic systems preserving both short-term accuracy and distributional measures, nonlinear resonance patterns in attention heads, and effective prediction of real-world experimental time series. Despite having been trained only on low-dimensional ordinary differential equations, Panda spontaneously develops the ability to predict partial differential equations without retraining. We also demonstrate a neural scaling law for differential equations, underscoring the potential of pretrained models for probing abstract mathematical domains like nonlinear dynamics.

Published as a conference paper at ICLR 2026
PANDA: A PRETRAINED FORECAST MODEL FOR CHAOTIC DYNAMICS
Jeffrey Lai" The Oden Institute
Anthony Bao* ECE Department UT Austin abao wetoouttexa.coomm edu
William Gilpin'
Department of Physics UT Austin
UT Austin
Ó jlai Putexas..ddi
gilpinechass. www.aaasssuuexa...o.. edu
ABSTRACT
Chaotic systems are intrinsically sensitive to small errors, challenging efforts to con- struct predictive data-driven models of real-world dynamical systems such fluid flows or neuronal activity. Prior efforts comprise either specialized models trained individual time series, foundation models trained vast time series databases with little underlying dynamical structure. Motivated by dynamical systems theory, we present Panda, Patched Attention for Nonlinear DynAmics. We train Panda on a novel synthetic, extensible dataset of 2 x 10 chaotic dynamical systems that we discover using an evolutionary algorithm. Trained purely on simulated data, Panda exhibits emergent properties: zero-shot forecasting of unseen chaotic systems preserving both short-term accuracy and distributional measures, nonlinear resonance patterns in attention heads, and effective prediction of real-world experi- mental time series. Despite having been trained only on oow-dimensional ordinary differential equations, Panda spontancously develops the ability to predict partial differential equations without retraining. We also demonstrate a neural scaling law for differential equations, underscoring the potential of pretrained models for probing abstract mathematical domains like nonlinear dynamics.
INTRODUCTION
Nonlinear dynamical systems test the limits of scientific machine learning (SciML). When an approximate model is constructed a chaotic nonlinear system, any small error grows exponentially over time, precluding long-ter ter-termml forecasting. This intrinsic property underscores the practical difficulty of accurately forecasting systems like weather fronts, neural activity, or economic markets (Li et al., 2022; Mikhaeil et al., 2022; Price al., 2025).
Recent empirical studies show surprising progress on the classical problem of forecasting chaos, including the ability to predict these systems well-beyond the classical predictability timescale for nonlinear systems (Gilpin, 2021; 2023; Pathak et al., 2018). These approaches construct local forecast models trained on past observations of a single dynamical system, and then forecast future, unseen states of the same system. For dynamical systems, this represents an in-domain generalization task, because future timepoints are drawn from the same underlying differential equations. This problem thus reduces to learning the numerical propagator for the true underlying governing equations.
However, a frontier in SciML. oft-o-ndomain generalization (Göring al., 2024; Wang et al., 2022):
Can a dynamics model effectively forecast unseen dynamical systems?
This task requires a global forecast model, which combines training on a large body of background knowledge with local adaptation to generate meaningful forecasts of unseen systems (Sen al., 2019). Moreover, what kind of data is required to train a forecasting model for dynamical systems in order achieve generalization? A global nonlinear forecast model has intrinsic theoretical interest in SciML, which has long questioned the degree to which complexity can be "transformed out" i.e. whether the predictability system is determined by its intrinsic properties or by the choice of measurement coordinates (Brunton et al., 2022; Mezić, 2013).
Equual contribution. 'Corresponiig author.

---

Published as a conference paper at ICLR 2026
A
B
Masked Competone &C. TAR
10 Human-curated ODE
8
Input
Tokenize
P999
Mutate
Generate
(parameter pitte)
Mutivariate Tragectory TI Fasches TOP,
Recombine (skew product)
f(x(())
g(x(())
Foresssta (C, H)
Select
Rejected
(chaoticit)
Dynamic Mode Embedding
Cross-Chanmn Transformer
Novel System g ([e01
Aftertion
Patch
Rand, Fourerr
Patch Enibedding
Temporal Allintion
10' Discovered ODE
Polyromia
Figure 1: large-scale chaotic dynamics dataset and dynamics-informet forecast model. (A) Evolutionary creation of large dataset of chaotic ODEs through mutation and recombination of known systems.
(B) Patch model architecture with forecasting and masked completion output modes. (C) The dlyamiic--nformee time series embedding and attention modules.
To address these questions, we introduce Punda' Patched Attention for Nonlinear DynAmics. Our key contributions are as follows:
1. We introduce a framework for generating novel chaotic dynamical systems, allowing us create dataset of ~ 2 10 ODEs, algorithnicallyddscoeee recombination of 129 chaotic systems such as the Lorenz attractor, double pendulum, etc.
based evolutionary
2. We pretrain global forecast model for nonlinear dynamics purely on chaotic trajectories integrated from our dataset. Our model exhibits competitive zero-shot forecasts for real systems including mechanical motion of C. Elegans, electronic circuits, and turbulent flows.
3. We demonstrate the effectiveness of features motivated by dynamical systems theory: (a) masked pretraining for dynamical continuity (b) channel attention for dynamical coupling, (c) kernelized patch embeddings based dynamic mode decomposition.
4. Despite being trained only low-dimensional ODEs, Panda develops emergent ability to zero-shot forecast high-dimensional PDEs.
2 RELATED WORK
Machine learning for dynamical systems. Machine learning models for dynamical systems (MLDS) leverage as inductive biases the unique properties of dynamical systems, relative to traditional time series. These include: (1) Strong channel coupling: The evolution of system variables is governed by deterministic differential or difference equations, implying coupled functional dependencies among variables rather than statistically correlations. Several MLDS approaches perform large- scale multivariate dynamical modeling, or infer interactions networks among measurement channels (Bhaskar al., 2024; Brunton al., 2022; Chen et al., 2018; Li al., 2022). (2) Invariant statistical measures: Ergodic dynamical systems possess invariant probability measures supported on non- wandering sets, such as limit cycles or strange attractors, resulting in well-defined long-term statistical distributions for all observables. Recent works incorporate these properties inductive biases in modern methods in MLDS settings (Cheng et al., 2025; Koppe et al., 2019; Pedersen et al., 2025).
Cooe available: https:/ //g1 tuwithub.cco com abaol 999 /panda
2

---

Published as a conference paper at ICLR 2026
A
SMAPE
12
70
10-
Kinda
Chicneo 2004
SMAPE
MAE
Prediction Length
MAE
30
20
Punce
Chrunss 2014
0.2
Chronos 200M
10
Time MOE S0M
TimesFM 2009
Prediction Lenget
Figure Panda zero-shot forecasts unseen nonlinear dynamics. (A) Example zero-shot forecasts on novel chaotic skew-systems. (B) sMAPE and MAE of Panda compared to zero-shot time series models over a 128 timepoint prediction horizon. (C) Error versus forecast horizon. Error ranges correspond to median and semi-interquartille range across 9.3 103 held-out dynamical systems, 6 forecasts per system. Note: i indicates some NaNs present in forecasts (more examples in Appendix C; dataset description in Section 3). See Table 6 in Appendix D for statistical significance tests.
Discovering new dynamical systems. Small datasets of dynamical systems have previously been curated from the published literature (Gilpin, 2021; 2023; La Cava et al., 2021). Several pretrained models, particularly for partial differential equations (PDE), generate new equations for training by randomly-perturbing parameters or initial conditions from known systems (Chen et al., 2024; Herde et al., 2024; Sun et al., 2025; Tripura & Chakraborty, 2023). Others construct de novo systems by combining terms from a fixed function library (Ziegler et al., 2024), or leveraging language models create candidate symbolic expressions (d' Ascoli et al., 2023; Du al., 2024). However, these approaches do not address the harder task of sampling based on whether a system exhibits a unique dynamical attractor. Richer sampling requires post-hoc analysis of candidate dynamics, akin to intrinsically- -motiicallymooivateed discovery previously used in domains such as cellular automata and coupled oscillators (Crutchfield & Mitchell, 1995; Falk al., 2024; Kumar et al., 2024; Reinke et al., 2020). Some foundation models generate synthetic time series using stochastic dynamics like Gaussian processes (Ansari et al., 2024; Das al., 2024), or simulated physics environments (Lin al., 2025; Wang et al., 2024).
Pretrained models for SciMI.. Pretrained foundation models dynamics enable transfer learning and zero-shot inference. One study trains transformers across diverse PDEs to create a shared multiphysics embedding space (McCabe et al., 2024). Another study proposes supervised pretraining to enable out-of-domain generalization for scientific foundation models, and derives scaling laws for transfer learning on PDEs (Subramanian et al., 2023). Several recent studies evaluate the zero- shot performance of time series and language models in MLDS, and observe performance only comparable to standard time series tasks (Liu et al., 2024; Zhang & Gilpin, 2025b). Several studies apply pretrained transformers to control symbolic equation discovery tasks (Becker al., 2023; Assooli et al., 2023; Lee et al., 2023; Zhang et al., 2024). One work generates pretraining data by randomizing the parameters four named ODE (Song et al., 2024), similar the first step of our evolutionary algorithm described below, with small founder pool. Another work samples systems from a fixed function space, selecting based on total variation over time (Ziegler al., 2024), while another study uses latent ODE prior for zero-shot imputation (Seifner et al., 2025). contemporaneous work to our study, DynaMix, zero-shot dynamical systems reconstruction (Hemmer & Durstewitz, 2025), built from Almost-Linear RNN experts and trained on the founder pool for our dataset (Brenner et al., 2024). Our work is distinguished by (1) rich data generation process, which discovers novel chaotic flows with diverse properties, and (2) a multivariate patched-basee baaehed--bsedd architecture which demonstrates emergent forecasting capabilities like zero-shot PDE inference.
multivariate mixture-of-expertt model for
3

---

Published as a conference paper at ICLR 2026
3 DATASET
Evolutionary search. We discover 2 x 10s novel chaotic ODEs (schematic in Fig. 1A, example systems in Appendix A). 1. Founding population: We start from a human-curated dataset of 129 low-dimensional chaotic systems (Gilpin, 2021), consisting of curated ODEs from the literature (e.g. the Lorenz equations blinking vortex flow) of the form X - fx 1). The default parameters of each system 0 and initial conditions x(0) were hand-tuned to the chaotic regime, and the integration timescales were standardized based on calculations of invariant mathematical properties of the underlying equations, such as the largest Lyapunov exponent. 2. Mutation: We randomly sample pairs of systems fa. f.. For each ODE's default system parameters, we add random Gaussian noise, 0% ~ (Ma,,, 7). 8% ~ N(0,,o). 3. Recombination: We combine the mutated parents using an additive seeeuurrouuctee coupling:
x (x, t)
(1)
y ky, + No tx,
(2)
This coupling between flows is asymmetric, and thus we refer to f, as the driver and f. as the response. In general, skew-product coupling maps be symmetric and nonlinear, but may be harder to integrate as a result. This particular recombination scheme, for appropriate scale factors, preserves chaoticity because the response system either synchronizes to the chaotic driver or continues to exhibit chaotic dynamics (Gilpin, 2025; Pecora & Carroll, 1990). For the scale factors, we compute the inverse RMS norm = 1/ V/vVEl f t) for each individual flow over representative trajectory. Selection: We integrate trajectories from multiple initial conditions using 5th order implicit Runge- Kunge-Kutta integrator (see Appendix A), and use a suite of attractor tests to cull systems that fail to exhibit chaotic behavior. First, transient systems that converge to fixed point or diverge to infinity are filtered. Then, we apply the chaos 0-1 test, which distinguishes quasiperiodic dynamics from true chaos (Falconer et al., 2007). We also apply near-recureence test to reject limit cycles, a power spectrum test to reject trajectories with only few distinct sharp peaks, and the data-driven Rosenstein estimator (Rosenstein et al., 1993) to ensure a positive maximum Lyapunov exponent. Finally, filter for stationarity using the Kwiatkowski-Phiile Wawlllps-Sehmiit---e- et al., 1992) and augmented Dickey-Fuller (ADF) (Dickey & Fuller, 1979) statistical tests.
(KPSS) (Kwiatkowski
Augmentations. On top of the integrated trajectories, we expand the training data by applying ponamics-inspires augmentations that preserve the property that the transformed timeseries arise from closed nonlinear dynamical system. Our augmentations are: Random time-delay embedding x; (1) - 2r(l - Ti), Ti ~ L (1, dembeed) This augmentation produces dynamics diffeomorphic to the original trajectory due Takens' embedding theorem (Packard et al., 1980; Takens, 1981). Given X E RCxT and d ~ V (dmin, dmaz). Convex combinations take random linear combinations of coordinates with coefficients sampled from a Dirichlet distribution: X <- MX R M E Rdxc, Mr. ~ Dir(alc). Affine transforms implement X + AX +b, [A bj E IRd((et1), [A b]ay ~ o/.v. We set dmnn- =3 3, dmax = 10, and dembed 10 for our experiments.
Standardization. For all trajectories, we apply instance-normalizatio to standardize the scales per channel. For integration, we standardize the integration horizon and granularity based on the number of timepoints (4096) and dominant timescale; note, however, that the numerical integrator ultimately dictates the stepsizes (Gilpin, 2021; Rosenstein et al., 1993). We observe no decrease in the range of invariant properties (maximum Lyapunov exponents, fractal dimension) across generations, suggesting that the starting population is sufficiently large and diverse (see Appendix A more details).
Held-out systems. For zero-shot test metrics, we evaluate on an unseen set of 9.3 10 systems. We form the test set by holding out a random subset of 20 systems from the 129 founding system population and ensure that none of these systems or their descendants (systems where the parent is the driver or response) appear in the training set. We then evolve these systems into the test set and include all skew product systems descended from these held-out systems.
MODEL ARCHITECTURE
Dynamical systems differ from traditional time series, and we introduce novel architecture motivated by dynamical systems theory (Fig. 1B). Time series foundation models with causal decoders that tokenize time series per-observation basis tend parrot motifs from their context,
4

---

Published as a conference paper at ICLR 2026
leading to over-confident predictions out-of-domain tasks (Olsson et al., 2022; Zhang & Gilpin, 2025a;b). Parroting is a useful emergent inductive bias when modeling invariant properties long forecasts is prioritized over accuracy otherwise known forecasting the climate. However, opt for an encoder-only, fixed prediction horizon forecast model that maximizes short-term pointwise accuracy, known as predicting the weather in SciMI..
Panda generalizes PStchTST, a transformer for univariate forecasting that tokenizes timeseries on a per-patch basis (Nie et al., 2022). Section 5.2 shows that univariate-only architectures underperform on dynamical systems, motivating channel attention. Moreover, patching admits an inductive bias for dynamical systems due to Takens' theorem which states that time-delayed copies of a dimennenii measurement of a dynamical system result in a multivariate time series that preserves key topological features of the true attractor (Packard al., 1980; Takens, 1981).
Patching. We tokenize length T trajectory T E CxT by patching it into a token sequence of size P patches with stride S so that in general, Tr.s E [ +1)xXP We choose stride S = so that the token sequences are Tp E (CX(T/P)XP We choose P = 16, unless stated otherwise.
Dynamics Embedding. We lift the patched multivariate timeseries to higher-dimensionall em- bedding space (dmodel) by concatenating each patch token P E RCx with random polyno- mial and random Fourier features. For random polynomial features with degree d, we sam- ple a fixed index set I C {1..... P]d of [I] = Npuly (number of features) d-tuples such that for E I: (P) = II Pa.l, = Pe.,, ..... Pe.ls. The random Fourier features sam- ple parameters W PP(NNa)))) b E RNaa22 such that Wij,b; NN o2) and F(P) = [sin(PW + b) cos(PW + b]] RMr,, where added across channels (Rahimi & Recht, 2007). The overall patch embedding is E(P) == [P ((P) F(P)] RCx(P+NNot We use degrees d E 32 and choose Npoly and Nor such that dmodel = = Npply + Nut = 512. The use of polynomial and Fourier features as lifted dynamics embedding is motivated by prior approximations of the Koopman operator via extended dynamic mode decomposition (eDMD) (Kutz et al., 2016; reservoir computers, which use polynomial features to
Williams et al., 2015) and mext-generatie forecast chaotic systems (Gauthier et al., 2021). See Appendix B for our choices of hyperparameters.
Temporal Attention. We mix information over the temporal dimension by taking the channel dimension as a batch dimension and performing self-attention with p-RoPE. (Barbero et al., 2025) (a modification of rotary positional encoding, RoPE (Su et al., 2023)) over the T/ P univariate patches of dimension dmedat For all experiments, we use RoPE wavelength of 500 and p 75%.
Multivariate Attention. Several time series foundation models are univariate, and thus, channel- independent; they solely employ temporal attention for information mixing (Nie et al., 2022). However, chaotic dynamical systems exhibit strong channel coupling. We demonstrate this em- pirically for the electronic circuits dataset in Fig. 4D, where show the benefit of channel attention as the coupling strength increases. We interleave channel attention layers without po- sitional encoding after each temporal attention layer. Each layer transposes the token sequence, treating the token dimension as a batch dimension and the channels as a set before self-attention ChannelAttention (Tt) = SelfAttention(TT (TT)) T. E RT/PPC daaad Temporal attention is fol- lowed by feedforward residual network, GeLU activations (Hendrycks & Gimpel, 2016), and RMSNorm (Zhang& & Sennrich, 2019). In the prediction head, processed tokens are aggregated along the sequence dimension T/ P and mapped with linear layer into a length H channel-wise forecast. The architecture is further described in Appendix B.
5 RESULTS
5.1 Panda ZERO-SHOT FORECASTS UNSEEN NONLINEAR DYNAMICS
To evaluate the quality of the generated forecasts, we measure (1) short-term forecast accuracy via mean squared error (MSE), mean absolute error (MAE), symmetric mean absolute percentage error (sMAPE), and Spearman correlation, as well as (2) attractor reconstrution accuracy via correlation dimension, spectral Hellinger distance, and Kullback-Leible (KL) divergence from the ground-truth attractor. For brevity, we report only the sMAPE and MAE (short-term), and KL divergence and spectral Hellinger distance (global) in the main text; the other metrics show similar results and
5

---

Published as a conference paper at ICLR 2026
included in the Appendix D and Appendix C. We compute all metrics for forecasts generated from zeroshot (held-out) systems never seen during training. Specifically, these are Nest = 9.3 x 10 unique skew-pree skew-product dynamical systems found using the methodology described in section 3. We additionally include results for scaling up model size and training dataset size in Appendix K.
Comparison to baseline models. We train Panda with 2/M parameters and evaluate against several time series foundation models comparable or larger scale: Chronos 20M, a causal univariate model which was recently shown to produce competitive forecasts of chaos systems (Ansari al., 2024; Zhang & Gilpin, 2025b). Chromos 20M SFT: Chronos supervised-fineed
on our entire chaotic systems dataset (Section 3). Time MOE 50M: A 50M parameter univariate model based on sparse mixture of experts (Shi et al., 2024). TimesFM 200M: A aatch-based 200M parameter decoder-only univariate model (Das et al., 2024). DynaMix: A parameter-efficte recurrent multivariate forecasting model pretrained on 3D dynamical systems (Hemmer & Durstewitz, 2025). For univariate models, each dimension forecast independently.
Across 9.3 x 10 held-out systems, we find Panda outperforms the baselines across variety of prediction horizons and error metries (Fig. 2). While we train our model exclusively on dimensional dynamical systems, the evaluation set includes arbitrary dimension systems, indicating that channel attention enables multivariate generalization. Moreover, we use window autoregression to extend our evaluation forecast horizon well beyond the forecast horizon used during training. Our model maintains its performance advantage, indicating that learns an effective dynamical propagator independent of single timescale. In Appendix D, show that results robust to the choice of metric (see Fig. 15).
A blations. We also ablate several features of Panda, in order to verify the contributions of our dynamics- based architectural choices. These include (1) Chan- nel Attention, (2) Dynamics Embedding, and (3) Masked Language Modeling (MLM) Pretraining.
Lprea
Lored
52.5
We observe a significant improvement due to chan- nel attention and MLM pretraining (See Section for example zero-shot completions). However, the
275-
47.5
combined effect the MLM with the dynamics em- bedding appears to be more complex: with no MLM, the dynamics embedding helps, but with MLM, it reduces performance. Moreover, the dynamics bedding improves the error on autoregressive rollout, whereas MLM reduces performance rollout. We conclude that using the dynamics embedding with polynomial features (PolyEmbed) gives us the best model for long prediction horizons.
22.-
Multivariate
Univariate
Chattn MLM
Chattn RFF
Univar MLM RFF
Chattn MLM RFF
Univar (wider) + RFF
PolyEmbed
(deeper)
Chattn
Univar (wider)
We include additional forecast metrics in Fig. 16 in Appendix D. We continue the discussion and evalua- tion of MLM the completions task in Appendix E. In particular, we compare the correlation dimension of the completions against that of the ground truth trajectories (Fig. 18) and show strong match. Fur- thermore, we investigate the effect of patch size on
Figure 3: Ablations of key architectural fea- tures of Panda: MLM pretraining, chan- nel attention (Chattn), and components of the dynamics embedding (RFF denotes ran- dom Fourier features and PolyEmbed includes polynomial features).
Panda's 's performance in Appendix J.
5.2 Panda ZERO-SHOT FORECASTS EXPERIMENTAL DATA
We next show that Panda generalizes to experimental time series from real-world dynamical systems. These experimental datasets have nonstationarity, missing values, noise, and other complexities not seen during training. Following prior works, we select systems in which the experimental data is known to have an underlying dynamical process generating it: the positions and momenta of the tips of the two rods in experimental recording of double pendulum (Asseman al., 2018), the leading independent components of body posture from a light microscopy video of C. elegans worms crawling agar gel (Ahamed et al., 2021), and voltage recordings from networks of 28 randomly

---

Published as a conference paper at ICLR 2026
Double Pendulum
EigenWorms
Electronic Circuit
Log SMAPE Ratio
Figure 4: Zero-shot forecasts of experimental data from (a) Double Pendulum (Asseman et al.,
2018), (b) Eigenworms (Ahamed et al., 2021), and (c) Electronic Circuit (Vera-Avila et al.,
2020). (d) Relative change in forecast error for Panda compared to Chronos-SFT (as measured in log (sMAPErunaa/ SMAPEEhhonssSPP showing the advantage of our approach as the coupling strength between variables increases, for various prediction horizons.
connected electrical oscillators (Vera-Ávila al., 2020). In all cases, the zero-shot performance of Panda outperforms Chronos-SFT (Fig. 4a).
For the circuit dataset in particular, find that the experimental coupling strength increases, the relative advantage of Panda over Chronos-SFT increases (red regions), particularly at long visible Pareto front between the two models (Fig, 4b). This finding
prediction horizons-leadien underscores the importance of channel attention for capturing nonlinear couplings typical in real world dynamical systems.
5.3 Panda EXHIBITS A DYNAMICAL SYSTEMS SCALING LAW
We create eight independent pretraining datasets that are subsets of the 2 10a unique systems generated using our methodology in Section 3. Across these eight datasets, we maintain a constant num- ber of total timepoints while taking. at extreme, a single trajectory (one initial con- dition) from each unique system, and at the other extreme, several trajectories (multiple initial conditions) from only few unique systems. These datasets thus allow us to measure how dynamical diversity (unique systems versus initial conditions) affects generalization. We repeat our zero-shot evaluations on our set of 9.3 x 103 held- out systems for each model trained on the eight datasets.
Scaling Law Fits
[[SMAPE] 28.12 151 67-N
-
Fitted Curve
42.5
40.0
SMAPE
37.5
35.0
32.
30.0
R2 0 9896
2/0 211 gl
g 3/s Number of Systems
1 g/o
Number Systems
Lyved
Lyred
192 Larer = 320 Lares m 448
Figure 5: Scaling laws in zero-shot forecast error the number unique dynamical systems increases. The total amount of training timepoints is held constant.
In particular, let Niss sampled initial conditions and N the number of unique systems. Keeping Nice x Neys fixed, our eight dataset splits are constructed (Nays R 156, Nice = 128)
number
as {Neye & 2 x 10*, Nics = 1). (Nswe RS 104 Nice = 2). where each subsequent split uses striet subset of the systems of the previous split, but with double the number of sampled initial conditions Nics
We observe clear scaling of zero-shot performance unseen systems with the number of new dynamical systems encountered. We emphasize that this scaling law is distinct from traditional neural scaling laws for total training data, because hold the number of timepoints constant while
7

---

Published as a conference paper at ICLR 2026
controlling the diversity of the data (Kaplan et al., 2020). These results highlight the advantages of scaling with diverse synthetic data. This finding accords with classical nonlinear dynamics theory: additional on-attractor trajectories continuously produce new information about that particular attractor's measure (a result of Pesin's 's theorem), but beyond certain point they fail to provide new topological information about winding, voids, etc (Gilmore, 1998; Pesin, 1977). The distinction between these information types partly explains the gap between in-domain and out-of-domain generalization in MLDS (Göring al., 2024).
5.4
Panda EXHIBITS EMERGENT PDE FORECASTING CAPABILITY
Kuramoto Sivashinsky
von Kärmán Vortex Street
Kurameto Sivashinsky
100
40
Truthh
Predictions
Dino ponar
40
Prediction
von Kärmän Vortex Street
2:251
Eroor
Predictions
40
Time
equation. The time axis t 25 to t > 44 contains 768 timepoints (512 context 256 prediction): solid black line marks end of context window, and dashed gray line marks length 128 prediction horizon. (B) Zero-shot forecasts of the von Kármán vortext street. (C) The horizoned MAE (with standard errors bars) compared to baselines. We show point-wise MAE, due to sMAPE's saturation the upper bound. We include two baselines, Fourier Neural Operators and DeepONet (black traces) fully-trained the context (see Appendix H).
Figure 6: (A) Zero-shot forecasts of the
Partial differential equations (PDEs) are dynamical systems on continuous domains, with diverse applications in weather prediction materials science (Kochkov al., 2024). Conceptually, PDEs may be seen as coupled ordinary differential equations evolving in an infiniteddimensional space. We apply our trained model the problem of forecasting two weakly-turbulent PDEs representing standard SciML benchmarks: the Von-Karman vortex street (VKVS) describes the unsteady motion of flow past a cylinder, and the Kuramoto Kuramoto-Sivashinksy (KS) models a laminar flame front (Cvitanović et al., 2010). More details on the PDE evaluation setup can be found in Appendix H.
Surprisingly, Panda outperforms baselines in zero-shot forecasting these systems (Fig. 6), despite having never encountered PDE during training. Unlike baselines, our model predicts nonlinear phenomena like merging of flame fronts in the KS equation or vortex pinchoff in the VKVS. While prior works require apecially-trained models to forecast chaotic PDEs (Pathak et al., 2018), our zero- shot approach does not require extensive in-distribution training data, highlighting the advantages of cross-channel attention and multivariate training in generalization.
5.5 Panda DEVELOPS INTERPRETABLE INTERNAL REPRESENTATIONS OF COMPLEX DYNAMICS
To probe the role of channel attention in Panda, feed two-tone sinusoids into the model and measure the response. The frequencies fi, f2 are each swept out over the range [2i,, 5mr). Let Á denote the attention rollout (Abnar & Zuidema, 2020) of the temporal attention matrices. Since Á is the product of row-stochastic matrices, Á remains row-stochastic. Thus, we measure the response
8

---

Published as a conference paper at ICLR 2026
Figure 7: (A) Nonlinear resonance structure measured by average row-wise entropy of temporal attention rollout matrices. (B) Mean row-wise entropy of the final layer during training. (C) Lack of nonlinear resonance structure the temporal attention rollout entropy for the univariate ablation.
from "shaking" the model at frequency mixtures f1, fz by measuring the average of the rowwise entropies of A (c.f. Fig. 7A). The attention maps exhibit complex, multiscale structure indicating nonlinear resonance, a phenomenon in dynamical systems where a physical system (such as kicked rotor or forced pendulum) exhibits gain with nonlinear dependence on the input frequencies. As control, the frequency response of an equally trained univariate model does not exhibit the same nonlinear multiscale structure (Fig. 7C).
next analyze Panda's attention maps to probe its underlying forecast strategy. The attention maps largely concentrate mass away from the diagonals, which indicates that Panda effectively uses the context. In contrast, a model implementing purely local rule (like a numerical integrator) would exhibit predominant diagonal structure, indicating that Panda performs more complex operations than few-step integration. For example, some attention maps form recurrence maps, which encode large-scale attractor geometry in classical nonlinear dynamics (Donner al., 2010; Gilpin, 2025). Other layers show banding and circulant structure (Fig. 8), consistent with global integral transforms like Fourier series.
WI IN
Www.
Figure 8: Temporal attention maps from Panda on context from different chaotic systems, showing Toeplitz, block, selector, and hybrid/combined structures (left to right). Appendix G further diseusses spatiotemporal coupling and cross-channel maps (Fig. 22).
5.6 LONG-TERM BEHAVIOR OF PREDICTIONS
Divergence of error is inevitable when forecasting chaos with finite precision. Eventually, the prediction error grows to the point where a point forecast is useless, but invariant and geometric properties of the chaotic system can still be estimated in the long horizon regime. We quantify the utility long horizon forecasts by measuring the geometry of forecasts much farther than 4x the training prediction horizon. Specifically, we compute: the maximum Lyapunov exponents (Appendix F.6); the forward KL divergence (Table 11) between the attractor and the predictions (Hess et al., 2023); and the spectral Hellinger distance (Table 12), an f-divergence between power spectra of the attractor and predictions (Mikhaeil et al., 2022). We compare long-term forecast metrics in Table 1, report mean + std. dev. across all test systems, averaged over 5 context windows for

---

Published as a conference paper at ICLR 2026
prediction horizons Lered {512. 1024, 2048, 3072). We report additional distributional metrics and computed invariant quantities in Appendix F.
We find that Panda outperforms transformer- based TSFM like Chronos and TimesFM, and that its performance remains competitive even when forecasts during testing are extended up to 8x the prediction horizon used for training. We attribute this finding partly to the pretraining data imposing an inductive bias for sustained long-term dynamics, particularly because Chronos-SFT improves over Chronos across forecast horizons and metrics. However, most transformer-based TSFM trained using a loss tied to short-term forecasts, including Panda, eventually regress to the mean at sufficiently-lnng enough forecast horizons (Appendix M). We can quantify this failure mode by computing the distributional metrics on the tail forecasts in Appendix F.5, However, Chronos avoids this outcome because was trained using cross-entropy on tokenized time series, leading to parrot sequences from context at long horizons (Fig. 31, 33, 32). This performs well on long-term metrics because ergodic dynamical systems have frequent repeated motifs known as unstable periodic orbits (Auerbach et al., 1987). A contemporaneous model to our work, DynaMix uses mixture of recurrent neural networks, allowing it to better capture long term geometry than transformer- based models (Hemmer & Durstewitz, 2025). Because experiment setting and dataset structure differs from this model's original setting, we highlight here comparisons to transformer based TSFMs.
DKL (Ground Truth pred) Model Prediction (prea)
Model
Irrdd
512
Lpredd = 1024
Iprrdd 2048
Lpred 3072
Ipred 3584
Panda
3.93 1 3.51
4.72 3.64
5.63 1 3.71
6.14 : 3.68
6.39 3.90
Chronos 20M SFT
4.72 + 5.00
5.09 + 4.90
5.62 + 4.86
5.93 + 4.84
6.05 + 5.34
Chronos 20M
5.99 5.07
6.19 1 4.85
6.51 1 4.76
6.76 1 4.74
6.94 5.41
Chronos 200M
5.12 5.25
5.22 5.49
6.05 1 5.30
6.36 + 5.28
6.47 5.67
A% (1)
+16.7%
+7.3%
+0.0%
-3.5%
-5.6%
Average H²[
H Ground Truth (Zeea)
Model Prediction (Lyea))
Panda
0.25
0.14
0.25
0.12
0.25 1 0.11
0.25 +:: 0.11
0.26 0.12
Chronos 20M SFT
0.29 + 0.17
0.29 + 0.16
0.29 + 0.16
0.30 = 0.16
0.30
0.18
Chronos 20M
0.37 + 0.16
0.36 + 0.16
0.37 + 0.16
0.38 + 0.16
0.38
0.17
Chronos 200M
0.28 + 0.16
0.28 + 0.15
0.29 + 0.15
0.30 + 0.15
0.30
0.17
A% (t)
+10.7%
+10.7%
+13.8%
+16.7%
+13.3%
Table I: KL divergence (top) and Hellinger Distance (bottom) between ground truth and model predictions. A% denotes percent gain of Panda over the best bascline. See Table 13 for per-system differences, and Appendix F.1 for implementation details.
CONCLUSION AND FUTURE DIRECTIONS
Our work demonstrates the feasibility of pretrained models in discovering generalizable properties of dynamical systems, mathematical objects of intrinsic interest to the SciML and forecasting communi- ties. Our model's emergent ability to predict nigher-dimensional partial differential equations, and the scaling of its performance with the diversity of dynamical systems, show that its generalization signal stems from unique properties of dynamics relative to time series.
A limitation of our work stems from our focus on low-dimensional dynamical systems. We argue that www-dimensional dynamics are the building block for sigher-dimensional systems like weather fronts spiking neurons, because they capture essential properties like bifurcations that become more complex in extended systems. A future variant of our approach for high-dimensional dynamics could exploit the structure of coupling such sparsity blocks typical in these systems by allowing the channel attention layers to receive custom attention masks. Another limitation is the degradation of rollout performance from MLM pretraining. Future work will investigate the question of what pretraining task is most natural for modeling dynamical systems. We believe this is basic question that necessitates further progress SciML.
10

---

Published as a conference paper at ICLR 2026
7 ACKNOWLEDGEMENTS
JL was supported by the UT CSEM Fellowship. AB was supported by the UT PGEF Fellowship and the Basdall Gardner Memorial Fellowship. WG was supported by NSF DMS 2436233 and NSF CMMI 2440490. This project has been made possible in by Grant No. DAF2023-329528 from the Chan Zuckerberg Initiative DAE, an advised fund Silicon Valley Community Foundation. Financial support for this publication results from grant 1A-CCS--2026-077 from Research Corporation for Science Advancement. Computational analyses were performed using the Biomedical Research Computing Facility at UT Austin, Center for Biomedical Research Support (RRID: SCR 021979). The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas Austin for providing computational resources.
REPRODUCIBILITY STATEMENT
Python code is available at https: http:///ggghaa.cccaaaa////////// Weights for the original model are available tttpps//aaaggaaffe....////baaaaa and the larger forecast httee///nntttss///ooonoooo////// 72M and MLM https: http:///tprp/////ooooooooooo///a
REFERENCES
Samira Abnar and Willem Zuidema. Quantifying attention flow in transformers, 2020. URL
https://arxe org/ abs /oss22055 00928.
Tosif Ahamed, Antonio C Costa, and Greg J Stephens. Capturing the continuous complexity of behaviour caenorhabditis elegans. Nature Physics, 17(2):255-28) -237, 2021.
Abdul Fatir Ansari, Lorenzo Stella, Ali Caner Turkmen, Xiyuan Zhang, Pedro Mercado, Huibin Shen, Oleksandr Shchur, Syama Sundar Rangapuram, Sebastian Pineda Arango, Shubham Kapoor, et al. Chronos: Learning the language of time series. Transactions on Machine Learning Research, 2024.
Alexis Asseman, Tomasz Koruta, and Ahmet Ozcan. Learning beyond simulated physics. In Modeling and Decision-making in the Spatiotemporal Domain Workshop, 2018. URL https:
//operrevi net //tt did-Hyllllll jW3RF7.
Ditza Auerbach, Predrag Cvitanović, Jean-Pierre Eckmann, Gemunu Gunaratne, and Itamar Procaccia. Exploring chaotic motion through periodic orbits. Physical Review Letters, 58(23):3387, 1987.
Federico Barbero, Alex Vitvitskyi, Christos Perivolaropoulos, Razvan Pascanu, and Petar Velicković. Round and round we go! what makes rotary positional encodings useful? In The Thirteenth International Conference Learning Representations, 2025, URL https: openrevie.. net/ forum?ormmmmmm wtwprrmmiid-Gttaeee... G
Sören Becker, Michal Klein, Alexander Neitz, Giambattista Parascandolo, and Niki Kilbertus. Predicting ordinary differential equations with transformers. International conference on machine learning, pp.
PMLR, 2023.
Dhananjay Bhaskar, Daniel Sumner Magruder, Matheo Morales, Edward De Brouwer, Aarthi Venkat, Frederik Wenkel, Guy Wolf, and Smita Krishnaswamy. Inferring dynamic regulatory interaction graphs from time series data with perturbations. In Learning on Graphs Conference, pp. 22-1. PMLR. 2024.
Manuel Brenner, Christoph Jürgen Hemmer, Zahra Monfared, and Daniel Durstewitz. Almost-linear rans yield highly interpretable symbolic codes in dynamical systems reconstruction. Advances in Neural Information Processing Systems, 37:36829-3660 -68668, 2024.
Steven L Brunton, Marko Budisić, Eurika Kaiser, and Nathan Kutz. Modern koopman theory for dynamical systems. SIAM Review, 64(2):229-340, 2022.

---

Published as a conference paper at ICLR 2026
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Wuyang Chen, Jialin Song, Pu Ren, Shashank Subramanian, Dmitriy Morozov, and Michael W Mahoney. Data-efficient operator learning via unsupervised pretraining and in-context learning. Advances in Neural Information Processing Systems, 37:6213 177611-6244,, 2024.
Xiaoyuan Cheng, Yi He, Yiming Yang, Xiao Xue, Sibo Chen, Daniel Giles, Xiaohang Tang, and Yukun Hu. Learning chaos in a linear way. In The Thirteenth International Conference on Learning Representations, 2025.
Maximilian Christ, Nils Braun, Julius Neuffer, and Andreas W Kempa-Liehr, Time series feature extraction on basis of scalable hypothesis tests (tsfresh-a python package). Neurocomputing, 307: 72-77, 2018.
Aaron Clauset, Cosma Rohilla Shalizi, and M. E. J. Newman. Power-law distributions in empirical data. SIAM Review, 51(4):661-703, November 2009. ISSN 1095-7200. doi: 10. 10.1137/070700111 URL http: / /xx doi wwa.oi..rorooooooooooooooooo 10. wtt.://..rrr01111177000000000000
James P Crutchfield and Melanie Mitchell. The evolution of emergent computation. Proceedings of the National Academy of Sciences, 1077 1995.
Predrag Cvitanović, Ruslan Davidchack, and Evangelos Siminos. On the state space geometry of the kuramoto-sivashinsk. flow in periodic domain. SIAM Journal on Applied Dynamical Systems, 911):1-33, 2010.
2010.
Abhimanyu Das, Weihao Kong, Rajat Sen, and Yichen Zhou. decoder-only foundation model for time-series forecasting. In Forty-first International Conference on Machine Learning, 2024.
Stéphane d' Ascoli, Sören Becker, Alexander Mathis, Philippe Schwaller, and Niki Kilbertus. Odeformer: Symbolic regression of dynamical systems with transformers. arXiv preprint arXiv:2310.00553. 2023.
David A. Dickey and Wayne A. Fuller. Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366):42--311, 1979. ISSN http:/staeeeee222333
01621459, 1537274X. URL http:: /www. j jstor.
Reik V Donner, Yong Zou, Jonathan Donges, Norbert Marwan, and Jürgen Kurths. Recurrence networks-a novel paradigm for nonlinear time series analysis. New Journal of Physics, 12(3): 033025, 2010.
Mengge Du, Yuntian Chen, Zhongzheng Wang, Longfeng Nie, and Dongxiao Zhang. Large language models for automatic equation discovery of nonlinear dynamics. Physics of Fluids, 36(9), 2024.
Ian Falconer, Georg A Gottwald, Ian Melbourne, and Kjetil Wormnes. Application of the 0-1 test for chaos to experimental data. SIAM Journal Applied Dynamical Systems, 602):395-402, 2007.
Martin Falk, Finnegan D Roach, William Gilpin, and Arvind Murugan. Curiosity-driven search for novel nonequilibrium behaviors. Physical Review Research, 6(3):033052, 2024.
Daniel Gauthier, Erik Bollt, Aaron Griffith, and Wendson AS Barbosa. Next generation reservoir computing. Nature communications, 12(1):1-8, 2021.
Robert Gilmore. Topological analysis of chaotic dynamical systems. Reviews of Modern Physics, 70 (4):1455, 1998.
William Gilpin. Chaos as an interpretable benchmark for forecasting and data-driven modelling. NeurIPS, 34, 2021.
William Gilpin. Model scale versus domain knowledge in statistical forecasting of chaotic systems. Phys. Rew. Research, 544:04425 2023.
William Gilpin. Recurrences reveal shared causal drivers of complex time series. Physical Review X. :5 1005, 2025.
12

---

Published as a conference paper at ICLR 2026
Niclas Alexander Göring, Florian Hess, Manuel Brenner, Zahra Monfared, and Daniel Durstewitz. Out-of-domain generalization in dynamical systems reconstruction. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 16071-16114. PMLR, 21-27 Jul 2024. URL html.
ht pttts://proceede/ www.g.smm..........
mlr. wtt.e/eeee///////ee/eeee
Peter Grassberger and Itamar Procaccia. Characterization of strange attractors. Phys. Rev. Lett., 50:346 339, Jan 1983a. doi: 10.1 100.1000oyyyseeeeeet....... URL
org/
50.346.
Peter Grassberger and Itamar Procaccia. Estimation of the kolmogorov entropy from a chaotic signal. Phys. Rev. A, 28:2591 Oct 1983b. doi: 10.1 10s//hyyseee https: 2591.
//link. org /tog 11tp//...1///////////ee
Christoph Jürgen Hemmer and Daniel Durstewitz. True zero-shot inference dynamical systems preserving long-term statisties. Advances in Neural Information Processing Systems, 38, 2025.
Dan Hendrycks and Kevin Gimpel. Gaussian error linear units (gelus). arXiv preprint
arXiv:1606.08411 084411 2016.
Maximilian Herde, Bogdan Raonic, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes. Advances in Neural Information Processing Systems, 37:72525-72624, 2024.
Florian Hess, Zahra Monfared, Manuel Brenner, and Daniel Durstewitz. Generalized teacher forcing for learning chaotic dynamics. In Proceedings of the 40th International Conference Machine Learning, JCMI.223 JMLR.org, 2023.
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXiv:2001.0836/ 2020.
Dmitrii Kochkov, Janni Yuval, Ian Langmore, Peter Norgaard, Jamie Smith, Griffin Mooers, Milan Klöwer, James Lottes, Stephan Rasp. Peter Düben, et al. Neural general circulation models for weather and climate. Nature, 632(8027)
2024.
Georgia Koppe. Hazem Toutounji. Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fmni. PLoS computational biology, 15(()=e0007663, 2019.
Jean Kossaifi, Nikola Kovachki, Zongyi Li, Davit Pitt, Miguel Siu-Schiaffiii, Robert Joseph George, Boris Bonev, Kamyar Azizzadenesheli, Julius Berner, and Anima Anandkumar. A library for learning neural operators, 2024.
Nikola B. Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew M. Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. CoRR, abs/2108.08.. 084881, 2021.
Akarsh Kumar, Chris Lu, Louis Kirsch, Yujin Tang, Kenneth O Stanley, Phillip Isola, and David Ha. Automating the search for artificial life with foundation models. arXiv preprint arXiv:24/2. 17799, 2024.
Jose Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor. Dynamic Mode Decomposition : Data-Drits DDiven Modeling of Complex Systems. SIAM, 2016. ISBN 9781611974492 1611974 URL http://////// /ww.. dmdbook. com/.
Denis Kwiatkowski, Peter CB Phillips, Peter Schmidt. and Yongcheol Shin. Testing the null hypothesis of stationarity against the alternative unit root. Journal of Econometrics, 54(1-3): 159-178, 1992.
13

---

Published as a conference paper at ICLR 2026
William La Cava, Bogdan Burlacu, Marco Virgolin, Michael Kommenda, Patryk Orzechowski, Fabrício Olivetti de França, Ying Jin, and Jason H Moore. Contemporary symbolic regression methods and their relative performance. Advances in neural information processing systems, 2021 (DB1):1, 2021.
Jonathan Lee, Annie Xie, Aldo Pacchiano, Yash Chandak, Chelsea Finn, Ofir Nachum, and Emma Brunskill. In-context recision-making from supervised pretraining. In ICML, Workshop on New Frontiers in Learning, Control, and Dynamical Systems, 2023.
Zongyi Li, Miguel Sluu-Sciaffiii, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar, Learning dissipative dynamics in chaotic systems. In Proceedings of the 36th International Conference on Neural Information Processing Systems, pp. 16768-16781, 2022.
Yixin Lin, Jan Humplik, Sandy H Huang, Leonard Hasenclever, Francesco Romano, Stefano Saliceti, Daniel Zheng, Jose Enrique Chen, Catarina Barros, Adrian Collister, et al. Proc4gem: Foundation models for physical agency through procedural generation. arXiv preprint arXiv:2503.08593, 2025.
Toni JB Liu, Nicolas Boullé, Raphaël Sarfati, and Christopher J Earls. Llms learn governing principles of dynamical systems, revealing an in-context neural scaling law. arXiv preprint arXiv:4402.00795. 2024.
Edward Lorenz. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2): 1963.
Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. DeepXDE: A deep learning library for solving differential equations. SIAM Review, 63(1):208-228, 2021. doi: 10.1137/19M/274067.
Michael McCabe, Régaldo-Saint Blancard, Liam Parker, Ruben Ohana, Miles Cranmer, Alberto Bietti, Michael Eickenberg, Siavash Golkar, Geraud Krawezik, Francois Lanusse, et al. Multiple physics pretraining for spatiotemporal surrogate models. Advances in Neural Information Processing Systems, 37:1 17:1199011111 19935,, 2024.
Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger. Umap: Uniform manifold approximation and projection. Journal of Open Source Software, 3(29):861, 2018.
Igor Mezić. Analysis of fluid flows via spectral properties of the koopman operator. Annual review of fluid mechanics, 45(1):357-378, 2013.
Jonas M. Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with mns. In Proceedings of the 36th International Conference on Neural Informa- tion Processing Systems, NIPS '22, Red Hook, NY, USA, 2022. Curran Associates Inc. ISBN 9781713871088.
Yuqi Nic, Nam H Nguyen, Phanwadee Sinthong. and Jayant Kalagnanam. time series is worth 64 words: Long-tem forecasting with transformers. In The Eleventh International Conference Learning Representations, 2022.
Catherine Olsson, Nelson Elhage, Neel Nanda, Nicholas Joseph, Nova DasSarma, Tom Henighan, Ben Mann, Amanda Askell, Yuntao Bai, Anna Chen, et In-context learning and induction heads. arXiv preprint arXiv: 1895, 2022.
Edward Ott. Strange attractors and chaotic motions of dynamical systems. Rev. Mod. Phys., 53: 655-671, Oct 1981. doi: 10. 1 doi/10.
URL http//innkk..........epp org/
1tt//00.1003//ee/oodhyyss. RewModPhyMs 53.65...6.. 655.
Norman H Packard, James Crutchfield, Doyne Farmer, and Robert S Shaw. Geometry from time series. Physical review letters, 45(9):712, 1980.
Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Physical review letters, 12242) 2018
14

---

Published as a conference paper at ICLR 2026
Louis M Pecora and Thomas L Carroll. Synchronization in chaotic systems. Physical review letters, 64(8):821, 1990.
Chris Pedersen, Laure Zanna, and Joan Bruma. Thermalizer: Stable autoregressive neural emulation of spatiotemporal chaos. arXiv preprint arXiv:2503./8731, 2025.
Ya B Pesin, Characteristic Iyapunov exponents and smooth ergodic theory. Russian Mathematical Surveys, 32(4):55, 1977.
Ilan Price, Alvaro Sanchez-Gonzale Ferran Alet, Tom Andersson, Andrew El-Kadi, Dominic Masters, Timo Ewalds, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, et al. Probabilistic weather forecasting with machine learning. Nature, 637(8044):84-90, 2025.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In J. Platt, D. Koller, Y. Singer, and S. Roweis (eds.), Advances in Neural In- formation Processing Systems, volume 20. Curran Associates, Inc., 2007. URL https: httsss/ppocceeddgg neurips
httppperfiieeeppppp//0////////// pdf.
Chris Reinke, Mayalen Etcheverry, and Pierre-Yves Oudeyer. Intrinsically motivated discovery of di- verse patterns in self-organizing systems. In International Conference on Learning Representations, 2020.
Michael T, Rosenstein, James J. Collins, and Carlo De Luca. A practical method for calculating largest Iyapunov exponents from small data sets. Phys. D, 65(1-2):117-33,, May 1993. ISSN 0167-2789. doi: 10. 10110160166-227990039900000
URL /ttp::s///oo..oo/1..110.0.0
0167-2789 (93) 90009-P.
David Ruelle and Floris Takens. On the nature of turbulence. Communications in Mathematical Physics, 20(3):167 192, 1971.
Patrick Seifner, Kostadin Cvejoski, Antonia Körner, and Ramses J Sanchez. Zero-shot imputation with foundation inference models for dynamical systems. In The Thirteenth International Conference on Learning Representations, 2025.
Rajat Sen, Hsiang-Fu Yu, and Inderjit Dhillon. Think globally, act locally: deep neural network approach to high-dimennional time series forecasting. Advances in neural information processing systems, 32, 2019.
Xiaoming Shi, Shiyu Wang, Yuqi Nie, Dianqi Li, Zhou Ye, Qingsong Wen, and Ming Jin. Time- moe: Billion-scale time series foundation models with mixture of experts, 2024. URL https: //arxiv. org/ abs 2409. 16400
Zezheng Song, Jiaxin Yuan, and Haizhao Yang. Fmint: Bridging human designed and data pretrained models for differential equation foundation model. arXiv preprint arXiv:2404. 14688, 2024.
Steven H Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press, 2018.
Jianlin Su, Yu Lu, Shengfeng Pan, Ahmed Murtadha, Bo Wen, and Yunfeng Liu. Roformer: Enhanced transformer with rotary position embedding, 2023. URL http::::s////xxvx.rr///eea 09864.
Shashank Subramanian, Peter Harrington, Kurt Keutzer, Wahid Bhimji, Dmitriy Morozov, Michael W Mahoney, and Amir Gholami. Towards foundation models for seientific machine learning: Charac- terizing scaling and transfer behavior. Advances in Neural Information Processing Systems, 36: 71242-71262, 2023.
Jingmin Sun, Yuxuan Liu, Zecheng Zhang, and Hayden Schaeffer. Towards foundation model for partial differential equations: Multioperator learning and extrapolation. Physical Review E, 111(3): 035304, 2025.
Floris Takens. Dynamical systems and turbulence. Warwick, 1980, pp. 366-381, 1981.
15

---

Published as a conference paper at ICLR 2026
Tapas Tripura and Souvik Chakraborty. A foundational neural operator that continuously learns without forgetting. arXiv preprint arXiv:2310. 2023.
VP Vera-Avila, Ricardo Sevilla- Eailla--scobza,, AA
RR Rivera-Durde and Javier M
Buldú. Experimental datasets networks nonlinear oscillators: Structure and dynamics during the path to synchronization. Data in brief, 28: 105012, 2020.
Rui Wang, Robin Walters, and Rose Yu. Data augmentation vs. equivariant networks: theory of generalization on dynamics forecasting, arXiv preprint arXiv:2206. 09450, 2022.
Yufei Wang, Zhou Xian, Feng Chen, Tsun-Hsuan Wang, Yian Wang, Katerina Fragkiadaki, Zackory Erickson, David Held, and Chuang Gan. Robogen: Towards unleashing infinite data for automated robot learning via generative simulation. In International Conference Machine Learning, pp. 51936-51 51983. PMLR, 2024.
Matthew O Williams, loannis G Kevrekidis, and Clarence W Rowley. A data-driven approximation of the koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science, 15:1307-1346, 2015.
Biao Zhang and Rico Sennrich. Root mean square layer normalization. In H. Wal- lach, H. Larochelle, Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (eds.), Ad- vances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL httpss:/prooeedinggs newrpip.s....nn............... http:ppperfffllppppppeeeee file/le88 9029222222222 2aw.a2222222eeeeeeeeeeeeeeeeeee2
Paper,,.ee Paper.pdf.
Xiangyuan Zhang, Weichao Mao, Haoran Qiu, and Tamer Başar. Decision transformer foundation model for partially observable continuous control. arXiv preprint WxXiv2244.002077 2024.
Yuanzhao Zhang and William Gilpin. Context parroting: A simple but tough-to-beat baseline for foundation models in scientific machine learning. arXiv preprint arXiv:2505.11 11349, 2025a.
Yuanzhao Zhang and William Gilpin. Zero-shot forecasting of chaotic systems. In The Thirteenth International Conference on Learning Representations, 2025b. URL pesttps//arxiv..oo/ abs/ /abs/2405
Martin Ziegler, Andres Felipe Dosada-Moreno, Friedrich Solowjow, and Sebastian Trimpe. On foun- dation models for dynamical systems from purely synthetic data. arXiv preprint arXiv:2412.003355 2024.
16

---

Published as a conference paper at ICLR 2026
A GENERATION OF A NOVEL CHAOTIC SYSTEMS DATASET
A.1 SKEW-PRODULE SYSTEMS
We algorithmically discover skew-produs systems following the methodology described in Section 3.
Here, we present a subset 30 of these systems, out total of 2 x 10 in our training set.
Figure 9: Examples of novel chaotic skew-product systems discovered via evolutionary search. Shaded regions correspond to two-dimensional projections onto the corresponding axes.
Our starting point is hand-curated, crowdsourced public dataset of 129 chaotic low-dimensional dynamical systems from the nonlinear dynamics literature (Gilpin, 2021; 2023; Zhang & Gilpin, 2025b). Each entry comprises set of coupled ordinary differential equations with dimensionality between three and ten. The parameters and initial conditions for each system have been hand-tuned into the chaotic regime, based values used in meviously-pubbiste studies.
17

---

Published as a conference paper at ICLR 2026
A.2 MUTATION OF BASE SYSTEMS
We also generate new instances of the base 129 chaotic systems by perturbing the ODE parameters.
Figure 10: Examples of parameter perturbations of base systems. Top row Unperturbed original systems. Botton rows Parameter perturbations of the top row systems.
A.3 NUMERICAL INTEGRATION
For all ODEs, we integrate trajectories of 4096 timesteps with an integration time-span dictated by system's characteristic timescale based off the dominant modes the power spectrum. We thus call this timescale the period of the system and set the integration upan to be [0. N2 y) where Np is the number periods, and yo is the "period" measured from integrating test trajectories of the base system; for skew systems we take the period to be the maximum period between the driver and response systems. For all experiments in the main text, we take Np = 40, but use larger mixed-period dataset in our scaled up experiments (Appendix K).
The numerical integration timestep controlled via adaptive step- seep-sizing from the Radau solver, a 5th order implicit Runge Kutta Scheme. For high quality trajectory data, integrate using relative tolerance x 10-9 and an absolute tolerance systems are obtained by integrating test trajectory at lower tolerance (rtol 1e 6, atol 1e - 7) and sampling a point from the coarse trajectory which approximates starting at a point on attractor.
10-10, The initial conditions for discovered
A.4 ATTRACTOR SELECTION
The only general way to identify properties about chaotic dynamical systems is to integrate them. This fundamental fact makes the system discovery process described in Section 3 very expensive at scale. To effectively reduce the number of incoming candidates chaoticity selection and validation, we employ callhacks during integration that will immediately kill the process and prune that system candidate. Specifically, we terminate integration whenever the step size falls below 10-10, any bounded (non-driving dimension) coordinate exceeds 10a in value, and whenever the integration time exceeds minutes. The surviving systems will finish integration and move on to the chaoticity selection phase (see the overview of our selection for chaoticity in Section 3).
18

---

Published as a conference paper at ICLR 2026
A.5 DATASET PROPERTIES
We verify that our integrated trajectories exhibit chaotic dynamics by measuring the average number of Lyapunov times in various prediction horizons. A chaotic flow separates nearby initial conditions according to ((() R (xppA11))6(0 where d the time dependent separation and À the maximum Lyapunov exponent. Thus, Lyapunov time is defined to be Tyyap 1/A1. Given arbitrary time series with timestep öt, the Lyapunov times per N timepoints is then A x åt N. Since we rely on implicit integrator with adaptive step-sizing, timespan and estimate the maximum Lyapumov exponent using the Rosenstein estimator (Rosenstein al., 1993) compute the distribution of Lyapunov times per horizon length in Fig. 12a where it is clear that we are predicting in the chaotic regime most of the time.
compute the average timestep
the integration
To ensure consistency between the founder population and offspring, we featurize all pretraining time series using the same procedure as previous works reporting chaotic systems datasets (Gilpin, 2021; 2023). For each channel of a D-dimensional multivariate time scries (4096 timepoints, 100 points per dominant Fourier period) we compute vector of 749 standard time-shift invariant time series features like wavelet modes, signal power, reversion rate, etc. using the t tsfresh library (Christ al., 2018).
We average the D feature vectors for each system to produce a single channel- permuttttoooooooonnnn invariant feature vector for each skew-product system. We then project all 10 pretraining skew-product systems into 2D using UMAP, nonlinear embedding algorithm that preserves the local neighbors of each point from the high dimensional space (distances, however, are not necessarily preserved) (McInnes et al., 2018). We next featurize and embed the 135 parent systems from the founder population into the same space. We observe broad dispersion the parent systems across the child population, implying the absence of mode collapse or strong distribution shift between the parents and offspring (Fig. 11). We interpret this result as the absence of strong founder or bottleneck effects in the offspring generation.
10
Offspring
Parent
Lorenz
UMAP2
6
10
Figure 11: A low-dimensional embedding of the 2 x 10 skew-product systems (gray), as well as the 135 founding parent systems (magenta) from which these offspring systems are evolved. The well-known chaotic Lorenz attractor is starred on the plot.
Additionally, for all trajectories in the test set, we measure the empirical stiffness score defined as at time
S log10 (max, /A//E, Ar)), where A, is the finite difference (forward or t. Fig. 12b shows that most test systems have at least an order of magnitude scaling between the largest observation jump compared to the average change per channel. This distribution suggests that the dataset generation algorithm generates stiff systems and reinforces the fact that the integrated trajectories exhibit non-trivial dynamics.
19

---

Published as a conference paper at ICLR 2026
Lyapunov Times
Stiffness Score
w28 (72%
Lyw cvve
20r
Count
300
307
207
22
22
0.5
2.0
(a)
(b)
Figure 12: Dynamical properties of systems in the test set. (a) Distribution of Lyapunov times within pped timepoints; annotated with the proportion of systems which exceed Lyapunov time in the horizon. (b) Distribution of sitffness scores (log-ratio of largest delta compared to the average delta).
Base Systems Systems
System
GP Dim (mean + std)
Base Systems
2.09 0.27
Skew Systems
2.11 + 0.23
Table 2:
Correlation dimension for Base
2.5
3.5
(Founder) systems and Skew (Children) systems.
Figure 13: Distribution of correlation dimension Craasberger-Procaceee of skew systems and their founder (base) systems.
As shown in Fig. 13, our skew-product generation preserves the distribution of the correlation dimension, invariant quantity used as proxy for fractal dimension. This suggests that dataset does not suffer from a "founder effect" that would kill off diversity. Table 3 further presents a comparison of the Kullback Leibler divergence between the invariant measures of the attractors. for skew systems: with the same parents; with different parents; with one parent shared; and between parent (Base) and child (Skew) systems.
Metric
std
Neombos
Dri (Skew
Response)
5.35 + 5.81
10,000
DeL (Skew
Driver)
8.46 + 6.58
10,000
DKL DKI(Sken (Skew
Non-Parent)
9.01 + 6.64
10,000
De Dkı(Skew Intra)
3.39 + 5.11
10,000
Dwu Inter)
7.07 + 6.58
10,000
DKL (Base Intra)
2.54 4: 4.13
6,000
DR (Base Inter)
8.24 + 6.48
6,000
Table 3: DEL. between skew systems and: 1) param perts of response; param perts of driver; 3) param perts of non-parent system in the founder pool. (Skew Intra) DxL. between param perts of skew systems with the same parents. (Skew Inter) DEL between param perts of skew systems with different parents. (Base Intra) Drz between parameter perturbations of the same founder system. (Base Inter): DEL. between parameter perturbations of different founder systems.
20

---

Published as a conference paper at ICLR 2026
Table 4: Model Architecture
Table 5: Model Architecture (Continued)
Parameter
Value
Parameter
Value
Context length
512
Patch length stride
16/ 16
Prediction length
128
Rope percent
0.75
Hidden layers
8
Max wavelength
500
dioood
512
Poly features
120
FFN dimension
512
Poly degrees
2
Attention heads
RFF count
256
Activation
gelu
RFF scale
1.0
Pre-norm
True
Normalization
RMSNorm
Init std
0.02
B TRAINING DETAILS
technical difficulty of training a global multivariate model is forming batches of trajectories with mixed channel dimensions. We look to dynamical systems theory and note that it is well known that least 3 coupled dynamical variables are necessary for a system to exhibit deterministic chaos in continuous-time (Strogatz, 2018). To this end, we fix the dimensions of each input trajectory to 3 only during training by randomly sampling 3 channels from each multivariate trajectory to enable efficient batching. During inference time, we process the full multivariate trajectories. For the 21.3M parameter Panda model we use dmooell = dm = 512, Nueuts = Nuavrr = 8, Npoly = 120 with degree 2, and Not 256. For the 41.5M parameter model (Appendix K), use dmmdel dm 640, Nhends Nuayrr 10, Npoly = 156 with degree 2, and NoF 312. And for the 71.5M parameter Panda-722 (Appendix K), we use dmnstel- dmn 768, Nhends Nuwyers 12, Npoly 188 with degree 2, and Nur = 376. Additionally, data augmentations (Section 3) are uniformly randomly applied to 20% of the training trajectories.
We use a patch size (and patch stride) of 16. All models are trained with context length of 512, which corresponds to 32 patches, and use non-causal transformer encoder with 8 layers, each with dmodel = 512 and 8 heads. Each attention block maps (batch size, channels, patches, hidden) sized hidden state via:
H H RopeTTemporalAAtention RMSNorm(H) (H)
H H + ChannelAtttention o RMSNorm(H))
H +- H + FFN O RMSNorm(H)
Where the transpose is applied to the channel and patch (sequence) dimension.
For models optimized with masked language modeling (MLM) style pretraining (masking and reconstructing intermediate patch tokens), a linear head is used to infill masked patches. For the forecasting model, a prediction head aggregates the encoder hidden states via a mean along the sequence (patch) dimension and a linear layer maps this representation to fixed-length 128 forecast for all models. All models are trained with MSE loss and the Adam W optimizer with maximum learning rate 1 x 10 on cosine schedule with 10% warmup. Additionally, we train with gradient norm clipping at a value of 1.0. See Tables 4, for comprehensive details about model architecture.
The 20M Panda MLM models are trained for 200K iterations (~52 wallclock hours across 4 GPUs or 208 GPU hours) with a batch size of 1024 with 50% of tokens randomly masked out each batch.
The forecast models are trained for 100K iterations with batch size of 1024 and are optionally initialized with a pretrained encoder from an MLM model (Section 5.1). The 20M parameter forecasting checkpoints are trained for ~26 wallclock hours or 104 GPU hours. The Chronos-SFT models use considerably more memory during training permitting a batch size of 160 for 300K training iterations which required ~48 wallclock hours 192 GPU hours.
21

---

Published as a conference paper at ICLR 2026
C FORECASTS
www.n hammt bomann Vmmann Mnn Mann pmmm
w
vornenlas werual Munt wehraa
www.mmp
W wwN.wm W M WWMMM
nnaos Aukuya Ma mm my
WWW. WAAAAAAAAAA Maddehhaosshaaa
/ m b - -
Figure 14: Examples of zero-shot forecasts (Tepedd 256) on held-out chaotic dynamical systems.
For additional forecasts, see Appendix I. (Figs. 28, 29, and 30). All forecasts plotted are with prediction length Lped = 256.
D ADDITIONAL FORECAST METRICS
Table 6 shows statistical significance testing of Panda metrics against other baselines. Note that we do not report results for TimeMoE due to the presence of NaNs, and instead test against a 200M Chronos baseline in greedy decoding and probabilistic mode. Panda clearly achieves lower error across the board; the gap closes with Chronos 20M SFT but still remains statistically significant.
In Fig. 2 we presented the sMAPE and MAE comparison for Panda versus our baseline models. We now present more zero-shot forecast metrics, using the probabilistic forecasting mode for Chronos and Chronos-SPP SFT. When finetuning Chronos 20M the default top-k and top-p (nucleus sampling) and temperature settings. We use these same settings, top-k 50, top-p 1.0, temperature 1.0 for the Chronos probabilistic forecasting, aggregating our metrics over 10 sample forecasts per context window per system.
dataset (i.e. Chronos 20M SFT), used
Table 7 highlights the efficiency of Panda against Chronos. This 60x speedup is largely due the fact that Chronos represents time series tokens as quantized individual points, whereas Panda relies on patches.
22

---

Published as a conference paper at ICLR 2026
Table 6: Wilcoxon Signed Ranked test for Panda errors vs. baseline errors (Holm-Sidák adjusted p-values)
Model
Prediction
MSE
MAE
sMAPE
Horizon
p-value
statistic
p-value
statistic
7-value
statistic
L, 128
5.96 10 48 13-484
9.01 x 10
12 153
6.35 x 10
494
Chronos 20M SFT
I = 512
8.24
10
50 12:275
3.07
54 10432
3.34 x 10
51 252
L 128
5.81 x 10-49
12778
6.88 x 10
9182
1.20 x 10
23 770
Chronos 200M
L 512
5.33 10
10-48 48
13233
1.45 x 10-56
9323
1.36 x 10
310
L 128
6.93 10-41
16.851
3.09
10 33
21 198
1.30 x 10
44 464
Chronos 20M SFT Prob
L 512
3.83 10-41
16313
1,18 x 10 10-35
19
3.26 x 10
54 730
128
2.91 10
12577
1.22 10 55
9290
x 31
22 353
Chronos 200M Prob
L 512
2.37 x 10-38
18024
1.83 X 10-49
12541
7.73 x 10
37 624
L 128
4.05
10-55
9863
1.67 X 10 10-69
3863
4.89 X 10
972
TimesFM 200M
512
2.50
10-38
103
1.59 10-59
7949
2.66 x 10
3588
L 128
1.03
10-41
14499
9.66 X 10-58
7320
1.68 x 10-76
170
Time MOE 50M
512
2.58
10-29
152
3.14
10-52
9639
3.13 x 10
589
128
1.46
10-23 23
724
3.03
35
1.86 x 10 10
16 533
DynaMix
L 512
4.54 10-8
39 970
7.00 x 10 22
27 602
1.83 x 10-27
23 592
MSE
SMAPE
1 Spearman
36-
6..10
63-00
Panda
Chronos 20M SFT
Chronas 200M
TimesFM 200M
50M
Figure 15: Zero-shot forecast metrics for baselines, using probabilistic (10 samples) forecasts for the Chronos models. Dash-dotted lines indicate presence of NaNs for some systems.
Time per Forecast (s)
Model
Time (mean std)
Panda
0,031 0.001
TimeMOE 50M
0.336 0.060
TimesFM 200M
0.605 + 0.032
Chronos 20M
1.880 1 0.041
Chronos 200M
4.233 1 0.121
Table Inference time per forecast (T-pred 512), computed over N 1000 calls to each model, on single H100 GPU. Each model call uses context length $12 timesteps. from our multivariate data, which has variable number of channels (at least 3). The univariate models (Chronos, TimesFM, TimeMOE) treat the channels as batch dimension, for each call.
also provide median forecast metries with IQR for the metrics in Fig. 2 over multiple prediction horizons for the best baselines in Tables 9.
23

---

Published as a conference paper at ICLR 2026
sMAPE Median (P25, P75J
Model
Lyred 128
Lyrod 256
Lyree 512
Panda
27.6 [18.5, 39.3]
36.7 [26.2, 47.6]
46.3 [37.0, 57.0]
Chronos 20M SFT
30.3 [21.9, 40.0]
40.1 [30.3, 48.3]
48.8 [37.3, 56.8]
Chronos 200M
36.0 [26.4, 44.6]
44.6 [34.6, 52.9]
53.5 [42.8, 60.8]
Chronos 20M SFT Probabilistic
29.7 [21.3, 40.7]
39.4 [29.3, 48.7]
48.3 [37.8, 57.1]
Chronos 200M Probabilistic
36.4 [26.7. 44.7]
45.0 [34.2, 53.1]
53.8 [42.8, 60.6]
DynaMix
47.1 [37.9, 57.91
54.7 [46.8, 62.6] 60.8 [53.9, 65.8)
Table Median sMAPE and interquartile range [P25, P75]. Note: while we include DynaMix for completeness, it recurrent model (not transformer) that benefits from longer context lengths than the short context settings of our experiments. Its performance improves with additional context. In contrast, most transformer-basee models, including Panda, degrade when tested on longer context than they were trained on.
MAE Median [P25, P75]
Model
Ipred 128
Lyrod 256
Lxpred 512
Panda
0.35 [0.22, 0.54]
0.49 [0.35, 0.70]
0.65 [0.48, 0.84]
Chronos 20M SFT
0.56 [0.34, 0.80]
0.85 [0.52, 1.26]
1.25 [0.72, 2.14]
Chronos 200M
0.61 [0.41, 0.80]
0.86 [0.58, 1.18]
1,07 [0.75, 1.83]
Chronos 20M SFT Probabilistic
0,46 [0.29, 0.69]
0.65 [0.43, 0,93]
0.85 [0.57, 1,34]
Chronos 200M Probabilistic
0.55 [0.38, 0.73]
0.71 [0.51, 0.91]
0.84 [0.63, 1.12]
DynaMix
0.79 [0.60, 0.94]
0.94 [0.76, 1.06]
1.02 [0.88, 1.15]
Table Median MAE and interquartile range [P25, P75J. Note that the DynaMix comparison does not match the model's original experiment setting, see the note in Table 8
D.0.1 ADDITIONAL. METRICS FOR MODEI. ABLATIONS
In Fig. we presented a sMAPE comparison for several key ablations of our model. Here, we provide additional zero-shot forecast metries for these ablations, supporting our conclusion that our dynamics embedding with polynomial features (PolyEmbed) is best for long-horizon forecasting via rollouts.
242
3.5
67
Spearman
M6SE
SMAPE
3.25
0.89
322
6.8
022
Multivariate
Chattn MLM
Univariate Univar MLM
Univar (wider) Univar (deeper)
Chattn MLM
Chattn + PolyEmbed
Chatte
Chattn
RFF
Univar (wider)
Figure 16: Zero-shot forecast metrics for our ablation experiments.
24

---

Published as a conference paper at ICLR 2026
E MLM COMPLETIONS
wmw mwma
WWw. wWW.koo0000000
acmn
A www an www.e
Jymn taant
www.mm mmmmmm
A www.mmmmmmmmm
moo anl Momoomy boreooolooo womming born
MMw
MMM
M
M NNW M m
ww wwwmmm
Figure 17: Examples of zero-shot completions on held-out chaotic dynamical systems. Each com- pletion plotted was with a context length of 512 time points, with half the patches (patch length 16) randomly masked out in a channel-inconsiste parameter checkpoint, completing the masked-ou trajectories ie. 256 time points. See Appendix N for more examples from Panda MLM and from our scaled-up model Panda MLM-66M.
manner. These plots show Panda MLM, our 20M
We present examples of Panda MLM completions on our held-out test set Fig. 17. For more examples of zero-shot completions, see Appendix N.
For the completions task, we randomly mask out half of patches for each coordinate dimension separately i.e. channel-independen masking. We trained Panda MLM with patch length 16 and context length 512, so each context window has 32 patches on the time axis. But we can generate completions with any context length. We refer to the masked-out portions of the trajectory as the erasures. We seek to measure how the model learns the cross-channel coupling relationships and statistical dependencies.
In future work, we hope to investigate more sophisticated masking strategics, such as masking out contiguous blocks of patches and investigating channel-dependent masking. which is closer forecasting task. Recall from our discussion of Fig. that MLM pretraining reduces performance on autoregressive rollout (c.f. Fig. 16). Determining the optimal MLM pretraining objective for long horizon forecasts on autoregressive rollout remains an intriguing area to investigate.
To quantify the performance of our MLM checkpoint the completion task, we compute the correlation dimension (Fig. 18) of completions versus ground truth trajectories using the Grassberger- Procaccia algorithm (Grassberger Procaccia, 1983a;b). This algorithm was developed quantify the strangeness (Lorenz, 1963; Ott, 1981; Ruelle & Takens, 1971) of chaotic attractors via com- putable metric related the fractal (Hausdorff) dimension and information entropy. For the result in Fig. 18, we take the entire length 4096 trajectory for each of our 9.3 100 held-out systems and we randomly mask out (erase) half of the patches (patch length 16) in a channel-inconsistent manner.
25

---

Published as a conference paper at ICLR 2026
Correlation Dimension (Panda MLM-66M)
Let {2t C CR be time series of T points in D dimensions. First, we compute pairwise Euclidean distances (excluding j):
2.8
2.6
R = rig ||= 21l122 15i1 S T, +
2.4
Completions
Next, we select the scaling region. Let ((6%) and ((50%) denote the empirical 5th and 50th percentiles of R. Then truncate to:
2.0
1.0
R² = { E R ((5%) < < 7(50%)
1.6
Now denote : |Ra| and Tmin min r.
y 0.84x + 0.41 (R 0.91)
Following Clauset. Shalizi, and Newman (Clauset et al., 2009), we identify a power law fit using maximum likelihood estimation (MLE). Assume for r 2 mii that the distances follow p(r) = Zr where the normalizing constant, Then,
1.4
2.8
Ground Truth
Figure 18: Correlation dimension comparison held-out systems. Computed for ground truth and completions from Panda MI.M666 using the Grassberger-Procee Procaccia method, av- eraged over independent masks for each trajectory, using context length 4096 with half the patches (patch length 16) randomly masked out in a Stannel-incoosister manner.
ā 1+
Σ FER
In
In the Grassherger-Procaceii method one examines correlation integral with unbiased estimator:
2 C(r) T(T - 1)
Σ i2j
H(r - C(r) ~ r (r 0),
that D2 d'I(((r), d'p Fitting C(r) P2 is equivalent fitting the distribution of pairwise distances to power law, yielding D2 R as the estimated correlation dimension.
Panda MLM,, with 20M parameters, shows promise in recovering the correlation dimension, statistical invariant of the attractor, even when given much longer context (length 4096) than seen during traning (recall the context length for training was 512), and with half of the timesteps masked out (in patches) per dimension. We also trained scaled-up checkpoint, Panda MLM-66M (with details in Appendix K), which demonstrates improved performance (Fig. 18).
In 10, we present the comparison against interpolation baselines. For polynomial in- terpolation, we use numpy polyfit
fit polynomial the unmasked timesteps. and evaluate with numpy. polyval at masked positions. For linear interpolation, we use scipy. interpolater inteiintee interpid, with extrapolation for timesteps out- side the range of known values. For the piecewise cubic spline baseline, we use with k = 3 (cubic spline).
w www.inineeeettteeee make backe__nterp_spp__
Comparison with Baselines for Completions Task
Method
Loontezt = 4096, with 50% Erasure (in patches)
Panda MLM-66M
0.91
Panda MLM
0.78
Piccewise Cubic Spline
0.71
Linear Interpolation
0.61
Polynomial Interpolation (Deg 3)
0.21
Table 10: Coefficient of Determination (R²) between the correlation dimension (via Grassberger- Procaccia) computed the completions versus the full length 4096 ground truth trajectories. For each of our 9347 held-out test systems, we average across random seeds, which determine the timestep masks for each trajectory (in patches of length 16). See Appendix N for further discussion.
26

---

Published as a conference paper at ICLR 2026
F ADDITIONAL DISTRIBUTIONAL METRICS AND INVARIANT QUANTITIES
We compare long-term forecast metrics in Tables 11 and 12, we report mean + std. dev. across all test systems, averaged over 5 context windows for prediction horizons Lpred E {512, 1024, 2048, 3072). Only one context window is available for evaluating Lpred = 3584, as our dataset contains trajectories of length 4096. Here, we include comparisons DynaMix due its strength in maintaining long term statistics during long rollouts, but note that differences in inference settings and our focus on shorter context lengths and rollouts compared to DynaMix, a fair comparison not possible. We thus intend to primarily compare with other transformer-based foundation models.
DeL. (Ground Truth [ proa) Model Prediction T ))
Model
Lpred 512
Lpred 1024
Lpred 2048
Lpred 3072
Lpred 3584
Panda
3.93 + 3.51
4.72 + 3.64
5.63 + 3.71
6.14 + 3.68
6.39 3.90
Chronos 20M SFT
4.72 5.00
5.09 + 4.90
5.62 + 4.86
5.93 4.84
6.05 5.34
Chronos 20M
5.99 + 5.07
6.19 +
6.51 + 4.76
6.76 4.74
6.94 5.41
Chronos 200M
5.12 5.25
5.49 + 5.22
6,05 + 5.30
6,36 + 5.28
6.47 5.67
DynaMix
4.75 1: 5.70
4.90 1: 5.65
5.22 +: 5.72
5,40 + 5.70
5.51 de 6.13
Table 11: KL divergence between ground truth and model predictions. See Table 13 for per-system differences, and Appendix F.1 for implementation details.
Average 2 (Srronnd (Lyest) || Model Prediction (Lyed))
Model
Lpred 512
Lpred 1024
Lpred 2048
4pred 3072
Lpred 3584
Panda
0.25 + 0.14
0.25 + 0.12
0.25 + 0.11
0.25 0.11
0.26 0.12
Chronos 20M SFT
0.29 0.17
0.29 0.16
0.29 + 0.16
0.30 0.16
0.30 0.18
Chronos 20M
0.37 + 0.16
0.36 0.16
0.37 0.16
0.38 0.16
0.38 0.17
Chronos 200M
0.28 + 0.16
0.28 + 0.15
0.29 + 0.15
0.30 + 0.15
0.30 0.17
DynaMix
0.36 + 0.19
0.31 0.19
0.33 + 0.19
0.33 + 0.19
0.32 + 0.21
Table 12: Average per-dimension spectral Hellinger distance between ground truth and model predictions. See Table 14 for per-system differences. We use Welch's method for estimating the PSD. Note that the DynaMix comparison does not match that model's original experiment setting, see the note in Table 8
27

---

Published as a conference paper at ICLR 2026
F.1 IMPLEMENTATION DETAILS: KL DIVERGENCE VIA GMMS
Algorithm presents our implementation of the Kullback Leibler divergence between ground truth and model predictions. This is the implementation we for main evaluations (Tables 11 and 13), although in Subsection F.3 we also present results using an alternative implementation found in the literature. In particular, we construct GMMs by fitting Gaussians each point, with local scale parameter determined by the simplex neighbors algorithm.
Algorithm 1 KL Divergence Estimation via Gaussian Mixture Models
Require: Ground Truth X {x+ J Generated Predictions Y ( Number of Monte Carlo Samples small > 0 and tol > 0
1: Function DIIIMATEEKLDDVRREENNCH
Y, n), Step 1: Local bandwidth (scale) estimation
2: a < SIMPLEXNEIGHIBRRS/ X, k = 10)
3: O < SIMPLEXNEIGHBORS Y, k = 10)
// Step 2: Construct Gaussian Mixture Models
4: < GAUSSIANMIXTURE
X, covariances = σx
5: q <- MAUSSIANMIXTURE/MSOO
Y, covariances )
// Step 3: Monte Carlo KLD Estimate
6: {z.) T + p.SAMPLE(1+)
7: for i 1 to ng do
8:
Ps + p(z,)
9:
4 + a(z,)
q 6 max(qi,)) 8)
11:
<- log(p:/0g /1))
12: end for
13: return KLD
T Di=17
14: Function SIMPLEXNEIGHRORS(Z,,.
15: Let Z = ( with z; E Rd
Build (k++))-narrst neighbor search structure Z (e.g., using Euclidean distance)
17: for = 1 to n do
Query (k++1) nearest neighbors of z,, including itself
19: Discard the self-neighbor to obtain neighbors {Zi, 1/ with distances di,j
20: a; + FINDSIGMA ((d........d..).
Estimate local scale parameter
Let P: + miny d.j
22: end for
23: return (a:)
24: Function FINDSIGMA(d,
(di....,, de) are distances to k nearest neighbors
25: p f miny dj
26: Define A f max(d; 0) for j- 1,...,k
b ReI.U on shifted distances
27: Define objective
( = exp(-Aj/(o - (ogy(k)
0 > 0
28: Minimize (()) using ID optimization (e.g., root-finding in log with initial guess p).
29: Let a² be the resulting positive solution
30: return 0"
28

---

Published as a conference paper at ICLR 2026
F.2 PER-SYSTEM DIFFERENCES DISTRIBUTIONAL METRICS
We report mean + std. dev. of per-system differences across all test systems, averaged over context windows for prediction horizons Irred E {512, 1024, 2048, 3072).
Per-system Ditference in Dxz (Ground Truth (pred) Model Prediction (pra)) )) between Baselines
Comparison
Lirad 512
Lpeal 1024
Lpoat 2048
Lyeat 3072
Lpmot 3584
Chronos 20M SFT Panda DynaMix- Panda
0.79 + 4.60
0.36 1: 5.03
-0.01 5.14
-0.22 1: 5.09
-0.33 + 5.84
0.81 + 5.12
0.18 + 5.58
-0.42 5.59
-0.75 + 5.52
-0.88 + 6.07
Chronos 20M - Chronos 20M SET
1.27 = 4.64
1.09 +: 5.39
0.89 + 5.50
0.84 + 5.53
0.89 + 6.62
Chronos 200M Chronos 20M SFT
0.39 = 4.35
0.41 +: 5.36
0.43 # 5.61
0.44 = 5.73
0.42 = 6.59
Table 13: The (mean + std) of per-system diff. in KI. divergence between models, a fine- fine-grained view of Table 11. DynaMix and Chronos 20M SFT outperform Panda on very long prediction horizons.
Per Per-System Difference in Average (SGround Truth (Lnnd)
S| Model Prediction (Lyna)) between Baselines
Comparison
Ixpead 512
Lipeed
1024
Liped 2048
Lpeod 3072
Lpres 3584
Chronos 20M SFT Panda
0.04 + 0.19
0.04 + 0.18
0.04 1: 0.18
0.04 0.18
0.04 0.20
DynaMix- Panda
0.11 = 0.22
0.09 : 0.22
0.08 + 0.21
0.07 0.21
0.07 0.24
Chronos 20M Chronos 20M SFT
0.08 + 0.20
0.08 + 0.20
0.08 + 0.20
0.08 + 0.20
0.08 + 0.23
Chronos 200M Chronos 20M SFT
0.00 = 0.20
0.00 + 0.20
0.00 + 0.20
0.00 + 0.20
0.00 + 0.23
Table 14: The (mean + std) of per-system differences in average spectral Hellinger distance between models, a mne-grained view of Table 12 showing that Panda outperforms the baselines.
F.3 AN ALTERNATIVE KL DIVERGENCE IMPLEMENTATION (GEOMETRIC MISALIGNMENT)
In addition our GMM-based KL divergence implementation (Tables 11 and 13), also use the implementation of (Hemmer & Durstewitz, 2025) based on geometric misalignment.
DEL (Ground Truth (Lpre)) Model Prediction(4prea) )) via Geometric Misalignment
Model
Lpred 512
Lpred 1024
Lpred 2048
Lpreat 3072
Lprest 3584
Panda
2.82 +: 2.67
3.29 + 2.79
3.88 + 2.85
1.26 + 2.88
4.44 3.14
Chronos 20M SFT
2.52 + 2.63
2.81 + 2.94
3.09 + 3.16
3.25 + 3.28
3.34 : 3.72
Chronos 20M
4.33 + 3.20
4.67 + 3.53
5.03 + 3.76
5.24 + 3.88
5.37 + 4.40
Chronos 200M
2.96 + 2.86
+ +.19+3.155
3.47
3.37
3.64 +.49
3.73 + 4.01
DynaMix
3.06 1: 4.07
3.15 1 4.42
3.24 1 4.68
3.30 4.81
3.37 +: 5.40
A% (t)
11.9%
-17.1%
-25.6%
-31.1%
-32.9%
Table 15: KL divergence between the ground truth and model predictions. A% denotes percentage gain of Panda over the best baseline. See Table 16 for per-system differences. Here, the implementation of (Hess et al., 2023). However, note that our experiment setting does not match the DynaMix's 's original experiment setting, see the note in Table 8
Per-system Difference DEL (Ground Truth( Model Prediction (ppretd between Baselines
Compurison
Lyeet 512
Lyss
1024
Lyoi 2048
Lyeat 3072
Lyeat 3584
Chronos 20M SFT - Panda
-0.30 1: 3.50
-0.49 1 3.75
-0.79 3.95
-1.00 +: 4.01
-1.09 + 4.54
DynaMix Panda
0.23 + 4.71
-0.14 + 5.06
-0.64 + 5.25
-0.96 + 5.38
-1.07 + 6.01
Chronos 20M
Chronos 20M SFT
1.82 = 3.60
1.86 + 3.93
1.94
4.25
1.99 + 4.33
2.03 2.033:5.11
Chronos 200M Chronos 20M SFT
0.44 + 3.33
0.39 + 3.66
0.38 + 4.01
0.39 + 4.10
0.38 + 4.88
Table 16: The (mean + std) of per-system diff. in KL divergence between models, a fine-grained view of Table 15. DynaMix and Chronoy 20M SFT outperform Panda on very long prediction horizons.
29

---

Published as a conference paper at ICLR 2026
F.4 VISUALIZATION OF METRICS DISTRIBUTION ACROSS ALL TEST SYSTEMS
We visualize distributional metrics of different models in Figure 19 and 20, which quantify how well the forecasts align with the global invariant statistics of the underlying dynamical systems. We emphasize that a fair comparison to Dynamix is not entirely feasible due to the fact that it requires a much longer context length during inference and is only trained on 3 dimensional dynamical systems. We hypothesize that the latter fact explains the anomalous mode in Figure 19.
KL Divergence (Lane 1024)
KL Divergence (Lund 2048)
KL Divergence (L 3584)
Panda Chronos 20M SFT DynaMix
Panda
Panda
Chrones 20M SFT Dyna
Chrones 20M SFT DynaMix
301
Count
Count
Count
300
101
320
00
2.5
5.0
75
36.0
125
25.0
37.
00
25
5.0
75
100
32.5
1500
37.5
0.
25
500
75
1000
125
15.0
315
20.3
Figure 19: KL divergence (via geometric misalignment) between ground truth (Iep and model predictions (Lpred). See Table 15 for aggregate values. Note that Lpred 3584 is 28x the prediction length used in training Panda.
Avg Hellinger (1 1024)
Avg Hellinger (lpres
2048)
Avg Hellinger (Lared 3584)
300
Panda
Pands
Pando
Chronos 20M SFT
Chronos 20M SFT
Chrenos 20M SFT
600
Dynelix
DynaMiz
532
DynaMix
500
500
402
400-
Count
Count
922
300
200
200
200
252
300
1.00
102
8.0
02
3.A
3.5
3.0
0.
0.2
0.4
0.5
0.1
3.3
0.0
02
0.4
0.6
0.8
Figure 20: Average spectral Hellinger distance, between the power spectra of ground truth (Toped) and model predictions (Lpred). See Table 12 for aggregate values.
F.5 QUANTIFYING MEAN REGRESSION ON VERY LONG HORIZONS
Mean regression is a common failure mode for TSFMs on very long prediction horizons. To quantify this failure mode, we compute the distributional metrics at Leea 3584, which is the longest possible horizon for evaluation, since our dataset contains trajectories of 4096 timepoints. However, we cut off the first Newetf 1536 timepoints of model predictions and ground truth, to compute the metries on the last 2048 timepoints solidly within the mean regression regime.
30

---

Published as a conference paper at ICLR 2026
Metrics pred [1536, 3584] (Cut Off First 1536 Timepoints)
Model
KL Divergence (Geometric Misalignment) Spectral Hellinger Distance
Panda
15.25 + 2.16
0.49 + 0.11
Chronos 20M SFT
7.00 + 5.63
0.34 + 0.18
Chronos 20M
9.68 + 6.12
0,48 + 0.18
Chronos 200M
7.37 + 6.05
0.36 + 0.19
DynaMix
3.50 + 5.44
0.35 + 0.22
Table 17: Metries between ground truth and model predictions after cutting off the first 1536 time- points Lpedd = 3584 (keeping only the last 2048). We present (mean E std) across test systems.
F.6 MAXIMUM LYAPUNOV EXPONENT
Anas Distribution (Land 1024)
Distribution (Uur 2048)
Area Distribution (Crre: 3200)
107
Pancaa
Panda
Panda
DynaMix
Ground Truth
Truth
707
Count
Count
Count
10r
20"
1.1
100
1
2
2
Aras
Pess
Figure 21: Distributional comparison of the maximum Lyapunov exponents estimated from the ground truth (Eerea) versus estimates from the model predictions (Iered) of Panda and of DynaMix. Note that the y-axis on a log scale.
We compute the maximum Lyapunov exponents for long prediction horizons, using the data-driven Rosenstein estimator. Fig. 21, we compare the distribution of estimated (Amax) for Panda versus that of DynaMix at prediction lengths Lped 1024, 2048, 3200. Note that DynaMix was trained pointwise autoregressively for prediction horizon of 550 points. Despite Panda being trained non-autoregressively for 4x shorter prediction horizon, it is capable of producing forecasts which maintain the characteristic Lyapunov exponent out to 25 x the prediction horizon was trained on.
We do observe that Panda struggles to capture systems with Amar > 6 in Fig. 21. This is likely due to the failure mode of mean regression over long enough prediction horizons. As a recurrent model, DynaMix avoids this limitation at long horizons.
G IMPLICIT SPATIO-TEMPORAL COUPLING
Temporal attention and channel attention layers independently mix information along the patch and channel dimensions. For system like the Lorenz attractor with coupled phase coordinates {z, y, 2). we would ideally want information to mix across space and time. We will show that by composing temporal and channel attention in sequence, Panda implicitly performs tpatio-temporal coupling.
Let Wo, WK, Wy denote the learned projections for temporal attention and W Wo, WK K, Wv for channel attention. For simplicity, we will focus on the linear attention setting without the row-wise softmax. Let P E N be a stack of N, dimae--iimenionnal patch embeddings with C channels. and pj attention output is (PW W PWW. In vector form,
(c) Rddas individual patch embedding for patch i and channel The linear
T (TA): : -0000-- Σ
[ (WTp. 9 WWK,,") WTD,
(p,)) P ATA,P, WTD,
(3)
1=1
31

---

Published as a conference paper at ICLR 2026
(CA) 0i
(k) W q%i W
Wkd(")
W
Σ
(k) of
Aci,
W Wro!!
(4)
E=1
Mha
Where TA denotes temporal attention and CA channel attention, and ATA WoWk and ACA W q W K- Looking at an element of the 3-tensor M we see that:
T
Mat
(k) Pi
ATAP;
(1)) WVDj (6) AcA
10 (p/ ( WTp, (2)
(5)
j'=1
((p!") Aup, WTp, (k)
(p!)) AmDp))) Aca WT.,
(0)
(6)
3/'=1
= (p!*) AnAP) (p. (,00 (WVAcaW pj
T
(7)
=/=1
Ācs
Where Aca prescribes how patches from different channels attend to each other. In matrix form,
AcAp!") :
(p AcADy")
(p!²) ATAP) )
M (PTx) ) T
(b)
PTA
PTA
(pt Acspi
(pf)) (k) AcADT
(p! (TAP!) )
Cross-Chamen Channel Mixing Map
JM Madan
Figure 22: Cross-channel mixing maps across patches for different channels from different held-out systems. Each mixing map is max-scaled to the range [-1,1].
32

---

Published as a conference paper at ICLR 2026
H FORECASTS ON PDES
Other than foundation models, we include a Fourier neural operator (FNO) baseline that was trained using the neuraloperaton framework (Kossaifi et al., 2024; Kovachki et al., 2021), and DeepONet baseline trained with the deepxde framework (Lu et al., 2021). Unless otherwise specified, experiment parameters follow the default values in these libraries. For both operator learning baselines, we tune the parameters for each PDE and use the best training checkpoint according the validation loss. Both operator learning baselines are trained (on a single MI100X GPU) for one-step-ahead prediction on length 512 context windows and rolled out for 512 prediction steps for each evaluation window in Fig. 6c
For the Kuramoto-Shivan Shivashinsky (KS) PDE, we integrate the equations pseudospectrally with 64 Fourier modes and the spatial length parameter L 100. We use an explicit eighth-order Dormand- Prince scheme (DOP853) to integrate the discretized PDE with a relative and absolute tolerance of 1e-8 from = 0 to 100 save the trajectory at 4096 uniformly spaced timepoints. We sample 40 initial conditions from U, =0, window starting at the 1024-th timepoint for training and the following 512 for prediction/rollovt for each sample to produce the error bars in Fig. 6c. See Tables 18 for comprehensive details.
where we choose
0.1 and use the length 512 context
For the KarmKarman vortex street (VKVS) data, we use 4600 timepoints of velocity field data in the domain Ω [0, 2] x [0, 1] on 256 X 128 grid simulated via Lattice Boltzmann simula- tion at a Reynolds number of 450. We then compute the vorticity field via second-order finite difference and reduce the dimensionality by keeping the top 512 principal components. For eval- uation in Fig. 6c, we train on length 512 training context windows starting at the time indices {0, 1024, 2048, 3072] and cross-validate on the length 512 prediction windows starting at the time indices {512, 1536, 2560, 3584) (avoiding train-set leakage) to produce the error bars in Fig. 6c. See Tables 19 for comprehensive details on the operator learning baselines.
Component
Specification
Component
Specification
Model
Fourier Neural Operator (FNO)
Model
DeepONet
Modes
256
Branch Net
128, 256 x 6]
Hidden Channels
256
Trunk Net
[1, 256 x 6]
Layers
Activation
tanh
Activation
GELU
Initializer
He normal AdamW
Optimizer
AdamW
Optimizer
Learning Rate
x 10
Learning Rate
1 1x 10-3
LR Scheduler
Cosine decay
IRR Scheduler
Cosine decay
Epochs
5000
Iterations
2 x 106
Batch Size
512
Batch Size
512
Loss Function
Lg
Metric
Mean relative L error
(a) FNO configuration
(b) DeepONet configuration
Table 18: Kuramoto-Shivashinsky PDE Operator Learning Configurations.
Component
Specification
Component
Specification
Model
Fourier Neural Operator (FNO)
Model
DeepONet
Modes
512
Branch Net
[512, 512 x 5]
Hidden Channels
256
Trunk Net
[1. 512 x 5]
Layers
Activation
tanh
Activation
GELU
Initializer
He normal
Optimizer
AdamW
Optimizer
AdamW
Learning Rate
x 1
Learning Rate
1 10
LR Scheduler
Cosine decay
LRR Scheduler
Cosine decay
Epochs
5000
Iterations
1 x 10c
Batch Size
512
Batch Size
512
Loss Function
La
Metric
Mean relative I.2
(a) FNO configuration
(b) DeepONet configuration
Table 19: Von-Karman PDE Operator Learning Configurations.
33

---

Published as a conference paper at ICLR 2026
Compared to the foundation models, the operator learning baselines under-perfoni mostly since they are limited to a context and prediction length of 1 for one-step-ahert prediction in contrast to the much larger context and prediction lengths of foundation models. We do not claim that foundation models are superior operator learning methods, but merely aim to provide a baseline for the PDE problems. The dash-dotted lines in Fig. 6 indicate that these methods are not directly comparable.
1 COMPUTING AND HARDWARE REQUIREMENTS
All training runs were conducted on 4x AMD MI100X GPUs, cach with 32 GB of HBM2 memory. Inference was performed on a single AMD MI100X GPU.
J EFFECT OF PATCH LENGTH
To investigate the effect of patch length which we train a version of our model with various patch lengths. To isolate the effect of patch size, we remove our dynamics embedding for these ablations. This is because each patch gets embedded to dimension dmodet, making the dynamics embedding incomparable between models using different patch lengths. Keeping fixed compute budget, also must halve the batch size every time halve the patch length, as a tradeoff patch length x 1 //btcch size exists between the two quantities: half the patch length implies twice as many patches, all embedded to dimension dmodel- In the tables below, A% denotes percentage improvement of the best ablation over the next closest.
model's performance, we conduct an ablation study in
MAAPE Median [P25, P75]
Model
Ixpred
128
p 256
Ipred 512
Patch 4
26.6 [17.5, 37.4]
[2..,,
47.55
48.1 [39.0, 58.4]
Patch
28.6 40.2)
(27.77
[39.9, 58.7]
Patch 12
28.9 [19.4, 40.6]
37.8
[399 59.0]
Patch 16
29.1 [19.7,41.1 41.1]
37.6[27.8,
49.44
47.7 [38.5, 58.7]
Patch 24
30.1 [20.6, 41.4]
49.5] 37.9127.9,
47.0138.7, 58.3]
Patch 32
30.1 [20.3, 42.2)
37.7
(28.2,,9
46.6 57.4]
A% (1)
+7.0%
+3.7%
+0.9%
Table 20: Median sMAPE and interquartile range [P25, P75] for various patch lengths.
MAE Median [P25, P75]
Model
Lpred 128
Lpred 256
Lpred 512
Patch 4
0.321 10.1 0.509]
|0.32, 0.685
(0.41, 0.824]
Patch
0.35910.222 0.543]
0.498 [0.359, 0.710]
0.668 [0.502, 0.860]
Patch 12
0.364 (0.234, 0.559]
0.501 [0.366, 0.71
0.662 [0.490, 0.859]
Patch 16
0.380 (0.246, 0.561]
0.507 [0.365. 0.724]
0.656 [0.486, 0.846]
Patch 24
0.382 (0.245, 0.578]
10.37, 0.729]
0.668 [0.495, 0.857]
Patch 32
0.382 (0.248, 0.584]
0.519 10.376, 0.722
0.664 [0.492, 0.852]
A% (t)
+10.7%
+5.2%
+2.1%
Table 21: Median MAE and interquartile range [P25, P75] for various patch lengths.
1 (Spearman distance) Median [P25, P75]
Model
Lpred = 128
Lpred 256
Lired 512
Patch 4
(..1110.117, 0.349]
0.322 [0. 197, 0.464]
0.480 0.480 [0.326, 0.612]
Patch
0.244 10.138, 0.376]
0.35010.223 0.485]
01.495 J0.350, 0.637]
Patch 12
0.246|0.139 0.384]
0.357 [0.225, 0.507]
0.496 J0.338, 0.640]
Patch 16
0.261 [3, 0.412]
0.360 10.221, 0.504]
0.486 [0.341, 0.644]
Patch 24
0.264 (0.150. 0.408]
0.361 [0.235, 0.512]
0.492 [0.348, 0.643]
Patch 32
0.268 (0.150, 0.413]
0.364 [0.230, 0.513]
0.485 [0.342, 0.644]
A% (t)
+10.2%
+8.0%
+1.0%
Table 22: Median P and interquartile range [P25, P75] for various patch lengths.
34

---

Published as a conference paper at ICLR 2026
K SCALING UP
We scale up our model parameters and training to investigate the improvement in performance. Our 21M, 42M, and 72M parameter models have values of (Theads Tlayers mmddet (8, 8, 512), (10, 10, 640), and (12, 12, 768) respectively. For the scaled-up training, we had (Naeers: batch size per device, number of GPUs) set (400K, 1024, 4), (400K, 512, 6), and (800K, 384, respectively.
We also scaled up our MLM checkpoint to create Panda MLMM666 with (7hends Maayers dmodel) set to (12, 12, 768), trained for 800K iterations, with batch size 192, and on 6 GPUs.
Lyre
0/5
Panda 72M
Panda
Panda 21M improved
Panda 21M
Figure 23: Zero-shot metrics for scaled-up checkpoints of Panda with increasing number of parame- ters. Here, Panda 21M is our original model presented in the main body. For Panda 21M improved, we trained for 400K iterations (compared to the 100K iterations for our original Panda 2 21M), and on an improved dataset with = 20% more systems, which we also use for the Panda 42M training. For Panda 72M, we trained on a larger version of our improved dataset with 8 initial conditions per system and with mixed periods. For presentation, bars show semi-interquarti range (40th to 60th percentile); for numerical values medians and interquartile ranges, see Tables 23, 24, 25.
K.1 SCALED-UP BASELINES
We also scaled up the model parameters for Chronos well as the training for the Chronos SFT baseline. Hardware limitations prevented us from fine-tuning Chronos 200M and larger model classes. We observe that our model continues to beat the baselines.
For our scaled-up training of Chronos 46M SFT, we used (Naers = 400K, batch size per device 100, number of GPUs 6).
For Panda 72M and for Panda MLM-66M, we trained on a larger dataset with 8 initial conditions per system and with mixed periods.
MSE
MAE
SMAPE
1 Spearman
08
03
06-
03
Prediction Langth
Prediction Langth
Prediction Eangth
Prediction Length
Panda 72M
Panda 42M
Chronos 46M SFT
Chronos 200M
Chronos 46M
Figure 26: Zero-shot forecast metrics for scaled-up baselines, using probabilistic forecasts for Chronos and Chronos SFT. Dash-dotted lines indicate presence NaNs for some systems (4% of systems for Spearman).
35

---

Published as a conference paper at ICLR 2026
Panda 21M Forecasts
tin the ih la / y Mel mmm
Wyyhl Mo Apl po pampuie whd Mhhhulann
the
gomm mmm
mMMVVV
Scaled up training mixed period multiple IC
Panda 72M Forecasts
m my veel y 40m
My My Mw poomea inmar imman wrm mmm mmm
mat wall
Figure 24: Comparison of sample zero-shot forecasts between the Panda 21M model (8 heads, layers), and the Panda 72M model (12 heads, layers), with the latter trained on a larger dataset with 8 initial conditions and mixed periods. As reflected in the metrics of Fig. 23, the scaled-up model forecasts appear to decrease error and capture higher-frequens details.
43-
MAE
SMAPE
14
3MAPE
10
35
62
32
Length
Length
Panda 72M
Panda 42M
Chronos 46M SFT
Chronos 200M
Chronos 46M
Figure 25: Zero-shot forecast metrics for scaled-up baselines, using deterministic forecasts for
Chronos and Chronos SFT. For Panda 72M, we trained on larger dataset with 8 initial conditions per system and with mixed periods. Spearman correlation is not shown because of the high proportion of NaNs for the Chronos deterministic forecasts, which we attribute to mean regression.
36

---

Published as a conference paper at ICLR 2026
We present the metries shown in Fig. 25 and Fig. 26 in tabular form in Tables 23, 24, 25.
sMAPE Median [P25, P75]
Model
Lpeed 128
Lpred 256
Lpred 512
Panda 72M
20.6 [12.4, 32.91
29.4[20.0, 41,4]
39.9 [30.9, 50.6]
Panda 42M
26.2 [16.6, 37.8]
33.9 45.7]
43.5 [35.0, 54.9]
Chronos 46M SFT
27.0 [18.2, 36.9]
45.5] 36.0
44.4 [32.7. 54.3]
Chronos 200M
36.0 [26.3, 44.7]
44.6[34.6, 52.9]
53.5 [42.7, 60.8)
Chronos 46M
(28.8, 47.1]
[38.6 54.6]
56.0 [46.9, 61.8]
Chronos 46M SET Probabilistic
[18.0, 37.3]
46.0] 35.8
45.1 [33.8, 54.1]
Chronos 200M Probabilistic
36.4 (26.7, 44.7]
45.0 53.1]
53.8 (42.8, 60.6]
Chronos 46M Probabilistic
38.3 [28.6, 46.6]
47.5 (37.9, 54.7]
55.9 [47.2, 61.6]
Table 23: Median SMAPE and interquartile range [P25, P75] for scaled-up models.
MAE Median [P25, P75]
Model
Ipedd = 128
Lprdd 256
Lpred = 512
Panda 72M
0.26 [0.14, 0,45]
0.41 [0.26, 0.60]
0.57 [0.41, 0.78]
Panda 42M
0.34 [0.20, 0.55]
0.48 [0.33, 0.71]
0.64 [0.46, 0.85]
Chronos 46M SFT
0.49 [0.28, 0.72]
0.78 1.17]
1.15 [0.61, 1,95]
Chronos 200M
0.61 [0.41, 0.80]
0.86 [0.58, 1.18]
1.07 [0.75, 1,83]
Chronos 46M
0.66 [0.46, 0.85]
0.91 [0.63, 1.23]
1.15 [0.81, 1,95]
Chronos 46M SFT Probabilistic
0.41 [0.25, 0.64]
0.60 [0.377 0.89]
0.77 [0.50, 1.21]
Chronos 200M Probabilistic
0.55 [0.38, 0.73]
0.71 [0.51, 0.91]
0.84 [0.63, 1.12]
Chronos 46M Probabilistic
0.58 [0.41, 0.75]
0.74 0.93]
0.86 (0.69, 1.17]
Table 24: Median MAE and interquartile range [P25, P75] for scaled-up models.
p (Spearman distance) Median [P25, P75]
Model
Tpredd = 128
Tppead 256
Txpred = 512
Panda 72M
0.16 [0.07, 0.29]
0.26 [0.14, 0.40]
0.40 [0.26, 0.56]
Panda 42M
0.23 [0.11, 0.39]
0.33 [0.20, 0.48]
0.47 [0.32, 0.62]
Chronos 46M SFT Probabilistic
0.25 [0.14, 0.40]
(.405, 0.54]
0.54 [0.36, 0.68]
Chronos 200M Probabilistic
0.41 [0.25, 0.55]
0.53 [0.35, 0.67]
0.65 [0.48, 0.79]
Chronos 46M Probabilistic
0.45 [0.28, 0.60]
0.57 [0.40, 0.71]
0.70 (0.55, 0.82]
Table 25: Median 1 - and interquartile range [P25, P75] for scaled-up models.
K.2 DATASET WITH MULTIPLE INITIAL CONDITIONS AND MIXED PERIODS
For our scaled-up training, we used larger dataset of multiple initial conditions and mixed periods. We present a sample of this dataset in Fig. 27. We vary the number of periods (on Fourier timescale), from 20 to 100 to produce multiple periods, and carry out the numerical integration with up to 16 different initial conditions (although we only use 8 initial conditions per system for training, due to compute budget restrictions). We integrate the same set of 2 x 10 systems used in our training set. The scaled-ut up training thus allows us to assess the effect of varying the timescales present in our training data.
37

---

Published as a conference paper at ICLR 2026
Thenstop
Theestop
2005
3008
Theestee
Timestpp
Figure 27: Examples of systems from our mixed period multi-IC training dataset. Fach subplot shows
multiple (4) initial conditions for single system (integrated with different timescale).
38

---

Published as a conference paper at ICLR 2026
L ADDITIONAL FORECASTS
In Appendix C, we presented sample forecasts from Panda on our held-out test set. Here, we provide more forecasts. As done previously, keep the prediction length fixed Lpred 256 for consistency and clearer visibility. Our model was trained with Lpred 128, so these forecasts include an autoregressive rollout. Fig. presents more forecasts, and Appendix Section presents failure modes and comparison against baseline models.
www.mm mmn
www
MMmM hrach M
do pomtooon momet 2y mag Mo pon Aap
MMm
Awmmmm
Figure 28: Examples of zero-shot forecasts (Trpred 256) on held-out chaotic dynamical systems.
MM MM M tn MM MMMM MM "n
M0 www. foreeasoooroeeoo Matheys /e 1 Moppoly
Figure Examples of zero-shot forecasts (Tored = 256) on held-out base systems (from the founder pool, parents of the skew-produc systems).
39

---

Published as a conference paper at ICLR 2026
w mongemanaaa manhong amaayun Ma Aub maturrn Lab
M W woo h 11100000000000000000000000000000
www.mmm MOANA WMMM
100 107"~~00 wwwpanyna
m www.M MMMay Monas X
~ moymmmmmnn rr wam wwmmmmom MMM
Man War WWwAMm wwm.eemmamme www.mmmmmmm
wwm.moo
Figure 30: Examples of zero-shot forecasts (Lored 256) on held-out chaotic dynamical systems.
40

---

Published as a conference paper at ICLR 2026
M COMPARISON WITH BASELINE MODEL FORECASTS AND FAILURE MODES
We compare long-term (Trrdd = 512) forecasts between Panda and the Chronos SFT and Chronos baselines. The following plots highlight some failure modes of each model, and also emphasize the advantage of our multivariate approach. Clearly, a univariate model can do well on single channel (dimension) but fail respect the attractor geometry. The coupling between channels encodes important information.
museramann boutentessen "onson Min monemmonn sww pommamam
Mane Woll woth Mynth www.wma wat.oao
porots -
10
M
Mancon Mine MANGOAAN www. petmaiyeaa
www Whi
Mag Worf Myrthhg www.mommm
Figure 31: Long-term zero-shot forecasts (oyred 512) held-out chaotic dynamical systems. Comparison between Panda (Red), Chronos SFT (Blue), and Chronos (Purple).
Figure 32: Comparison (T.ped 512) between Panda (Red), Chronos SFT (Blue), and Chronos (Purple). An illustrative example of held-out system where Chronos appears to parrot (limit cycle), Chronos SFT does not respect the attractor geometry, and Panda mean regresses.
41

---

Published as a conference paper at ICLR 2026
P0 10 WoN
Figure 33: Long-term zero-shot forecasts (Lpred 512) on held-out chaotic dynamical systems. Comparison between Panda (Red), Chronos SET (Blue), and Chronos (Purple).
42

---

Published as a conference paper at ICLR 2026
N ADDITIONAL COMPLETIONS
In Appendix E, we presented a sample completions from our Panda MLM checkpoint on our held-out test Here, we provide more completions from Panda MLM (Fig. 34) and long-context completions from our scaled-up checkpoint Panda MLM-66M (Fig. 36). We also provide qualitative comparison between Panda MLM completions and piecewise cubic spline interpolation (Fig. 35) to further demonstrate the advantage of our method.
Wohoos/ MMMMMM wwwmmmm
wwmamm m mwmaamm mmmnnnn
Ma war
ا /
www mmm
Aummannn w
wwmmmmm
Ae NAMA
wwmmmm
V Mmmammo me u
www.mmannm yar
www ww
www.
Figure 34: Examples of zero-shot completions on held-out chaotic dynamical systems. Each com- pletion plotted was with a context length of 512 time points, with half the patches (patch length 16) randomly masked out in a channel-inconsiste manner. These plots show Panda MLM, our 20M parameter checkpoint, completing the masked-out trajectories i.c. 256 time points.
43

---

Published as a conference paper at ICLR 2026
Panda MLM
60m
Piecewise Cubic Spline
Meh
Au
10
Figure 35: Qualitative comparison between completions generated by Panda MLM (with 20M parameters) and by piecewise cubic splines. First (leftmost) panel provides example with context length 512 for clearer presentation; all other panels show context length 2048. Shaded red regions show the difference from the ground truth. Piecewise cubic spline interpolation is the most successful naive baseline, and although it achieves near competitive performance on preserving the correlation dimension (Table 10), is not competitive with respect pointwise error or attractor geometry.
44

---

Published as a conference paper at ICLR 2026
MHr Myffh wareere " Jyln Amaa
Figure 36: Examples of zero-shot completions on held-out chaotic dynamical systems. Each com- pletion plotted was with a context length of 4096 time points, with half the patches (patch length 16) randomly masked out in channel-inconsisten manner. These plots show Panda MLM-66M, completing the masked-out trajectories i.e. 2048 time points, despite only being trained on context length 512.
45


Reviews:
[{"review_id": "6W4ALviGiQ", "reviewer": "Reviewer_t2gr", "summary": "This paper introduces Panda, which is a large foundation model for chaotic time series data. To generate the training dataset, the authors introduce a novel, genetic algorithm-like method for discovering new chaotic systems. The model architecture is an extension of the PatchTST model, utilizing patching along the time dimension, feature augmentation with random polynomial and Fourier features, and the use of both temporal attention and channel attention. Comparing the trained model to previous time series foundation models of similar size, the authors find that Pandas gives more accurate predictions, both in the short-term exact trajectory prediction and in the long term predictions of the statistical properties of the attractor. Furthermore, Pandas was able to perform zero-shot forecasts for both experimental data (which contains imperfections not presented in the training dataset) and for PDE systems (which were not seen by the model since it is only trained on ODEs). The paper contains lots of additional results, such as the scaling law with respect to the diversity of the training dataset, and analyses of the attention maps learned by the model.", "questions": "1.\tI am curious about the additive skew-product coupling used in the recombination step for the data generation process. From equations 1 and 2, it seems that recombination will lead to a ODE of larger dimensionality than its parent, since given ODEs for x and y, recombination gives a ODE for the concatenated state [x, y]. If so, I would imagine that the dimensions of the ODEs in the dataset will grow over iterations due to the repeated recombination process. Is this indeed the case?\n\n2. The authors mention that the integration horizon and the granularity (which I assume indicates the sampling period of the trajectories) is standardized using the dominant timescale. How was this timescale determined? Are the authors looking at the dominant frequency in the frequency spectrum, or using the inverse of the maximum Lyapunov exponent?\n\n3. Looking the generated chaotic systems in A.1 and A.2, it seems that loosely there are families of related systems (single, double, triple scrolls, etc). I wonder if it would be possible to do a low dimensional embedding of the generated dataset (or if it is too expensive, some subset of it) to create a 2D plot that shows how different systems in the dataset are positioned (in the chaotic behavior space) with respect to one another, as done in Figure 1 of [1].\n\n4. How long does it take to train the Panda model? How does this compare to the competitors considered such as Chronos and Time MOE?\n\n5. On a similar note, how long does the data generation process take? \n\n[1] W. Gilpin. Chaos as an interpretable benchmark for forecasting and data-driven modelling. NeurIPS (2021).", "limitations": "", "rating": 10, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1.\tThe proposed model shows excellent performance in both short term trajectory predictions and long term statistical property predictions. It is carefully designed to predict chaotic systems in particular and contain lots of architectural considerations that will be useful to researchers working on similar problems.\n\n2. The automatized dataset generation process is also novel, and the suite of criteria used to automatically sift for chaotic systems seems quite useful. The final dataset of 20000 new chaotic ODEs is also, to my knowledge, first of its kind. Such large scale dataset will also help the further development of foundational models for chaos.\n\n3. Lots of experiments were performed to demonstrate the model's performance. The out-of-distribution results for experimental data and the PDE is quite impressive too, and suggests that indeed the model is able to properly forecast chaos, instead of just having memorized motifs in the training dataset.", "weaknesses": "1.\tThis is a very well written paper, and there seem to be no significant weaknesses.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gKiHggO9zO", "reviewer": "Reviewer_tpp7", "summary": "This work presents Panda, an attention-based foundation model (FM) dedicated to zero-shot forecasting chaotic dynamics. Using a novel synthetic dataset constructed from a library of known chaotic systems using principles from dynamical systems theory (DST), Panda is trained to forecast trajectories based on limited dynamical context. The model demonstrates robust generalization, successfully providing zero-shot forecasts for entirely novel chaotic systems, diverse real-world datasets, and even PDE problems. Panda is compared to other existing time series FMs in the field and demonstrates superior performance.", "questions": "1. What is the rationale for comparing context + ground truth to context + prediction for several experiments and measures throughout the paper? Why not simply ground truth vs. prediction? Doesn’t inclusion of the context heavily skew the outcome of empirical estimation of the Lyapunov exponent, for example? Especially in cases where $L_{pred} = 128$, where the estimated value is dominated by the context which is $4 \\times L_{pred} =  512$?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. Panda shows great scaling behavior both in the amount of data as well as number of parameters.\n2. Great to see DST motivated architectural design and the focus on setting a suiting inductive bias, which other FMs in the field lack.\n3. The manuscript is well written and easy to follow.\n4. Generating new chaotic systems using skew-product coupling and evolutionary search is a great idea and powerful method to generate vast amounts of valid data for this class of FMs. I think this greatly benefits research in the entire field.", "weaknesses": "My two main concerns are the following:\n1. The paper misses relevant literature which introduces another DS FM model: DynaMix [1]. Both methods train on a synthetic dataset comprised of low-d chaotic DS, both methods have the aim of zero-shot forecasting and both question the efficacy of existing TS foundation models. I think it would be highly valuable to the SciML community if Panda is compared to such a similarly specialized model, which seems to perform much better than context parroting architectures such as Chronos.\n2. “Panda exhibits emergent properties: zero-shot forecasting of unseen chaotic systems preserving both short-term accuracy and long-term statistics.” I find the claim that Panda preserves long-term statistics not sufficiently backed up by evidence. The main metrics employed to verify this statement are comparison of the fairly short prediction lengths (up to 512) w.r.t. the ground truth attractor (4096 time steps) in Section 5.4. On the same note, the authors identify regression-to-the-mean (or convergence to fixed point dynamics) as a failure mode of Panda. Proper evaluation of long-term dynamics would consist in performing *longer autoregressive roll-outs* in which the model, ideally, converges to the “climate” hinted by the context (see e.g. [1]). In this case, according to the argumentation of the authors that Chronos exhibits the more suiting inductive bias for climate forecasting (parroting), Chronos should actually perform better in both Hellinger distance as well as KL measure. The way these measures are currently used only strengthens the point that Panda is superior in preserving e.g. geometric features only in the *short-term*, also indicated by the consistently rising mean measure values with prediction length (see Table 1, 2, 10). A SciML motivated architecture should in the end aim to preserve the physics of the underlying system in the long-term, i.e. $T \\rightarrow \\infty$.\n\nOther concerns:\n1. Since the emergent zero-shot PDE inference is highlighted as a major contribution of the paper, the authors should include comparisons to baselines given by domain specific models trained on the context window (e.g. Neural Operator based architectures such as DeepONets [2] or FNOs [3])\n2. (Appendix E / Figure 15) ground-truth + completions to ground truth using the correlation dimension lacks baselines: As a reader I can not quantify how well naive imputation methods would perform in estimating the correlation dimension.\n\nMinor comments:\n- Fig. 2 misses labels A, B, C mentioned in the caption\n- Abbrv. “MLM” is first used in Section 5.1 but hasn’t been introduced in the main text before\n- Double citation in References: “Jonas Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with rnns. Advances in neural information processing systems, 35:11297–11312, 2022a.”  and “Jonas M. Mikhaeil, Zahra Monfared, and Daniel Durstewitz. 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, 2022b. Curran Associates Inc. ISBN 9781713871088.”\n\n**References**:\n\n[1] Hemmer, Christoph Jürgen, and Daniel Durstewitz. \"True zero-shot inference of dynamical systems preserving long-term statistics.\" arXiv preprint arXiv:2505.13192 (2025).\n\n[2] Lu, Lu, et al. \"Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators.\" Nature Machine Intelligence, vol. 3, no. 3, Mar. 2021. https://doi.org/10.1038/s42256-021-00302-5\n\n[3] Liu-Schiaffini, Miguel, et al. \"Neural Operators with Localized Integral and Differential Kernels.\" International Conference on Machine Learning. PMLR, 2024.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PBV6wMlEyh", "reviewer": "Reviewer_SotH", "summary": "This paper introduces Panda, a pretrained transformer model for forecasting chaotic dynamical systems. The key contributions are: (1) an evolutionary algorithm that discovers ~20,000 novel chaotic ODE systems through mutation and recombination of 129 base systems; (2) a multivariate patch-based architecture incorporating dynamics-inspired features. in particular, channel attention; (3) demonstration of zero-shot forecasting capabilities on held-out chaotic systems, real-world experimental data, and emergent PDE forecasting; (4) a scaling law showing that diversity of dynamical systems, rather than data volume, drives generalization performance.", "questions": "1. Please explicitly describe the inference procedure for univariate baselines on multivariate data. Do they process each channel independently and concatenate forecasts? How does Chronos-SFT handle multivariate data during fine-tuning?\n2. Why does MLM pretraining help at training horizon but hurt rollout? Do you have some intuition here?\n3. How can mean regression achieve better distributional metrics than \"parroting\"?\n4. What are the training costs?\n5. Have you tried computing Lyapunov spectra?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1. **Synthetic dataset generation**: The evolutionary discovery framework is interesting and well-executed. I think this is a valuable contribution.\n2. **Scaling law for dynamics**: The finding that diversity of unique systems matters more than total data volume is very interesting scientifically.\n3. **Architectural choices**: Channel attention for coupled systems and dynamics embeddings are good.\n4. **Experimental results**: The evaluations are comprehensive and the results strong compared to baselines.\n5. **Emergent PDE forecasting**: The zero-shot ability to forecast PDEs after training only on low-dimensional ODEs is impressive.\n6. **Computational efficiency**: Table 8 shows significant speedup over Chronos, the nearest baseline.", "weaknesses": "1. **Clarity**: While the paper is generally very well presented, there are a few places where clarity could be improved. For example, the fact that training uses only 3 randomly sampled channels while evaluating on arbitrary dimensions is buried in Appendix B, but for me this an important detail. In addition, the paper doesn't explicitly state how the univariate baselines are applied to multivariate data—do they process channels independently and concatenate forecasts? How does Chronos-SFT handle multivariate data during fine-tuning?\n2. **Mean regression problem**: Section 5.4 reveals that Panda regresses to the mean on long rollouts, yet the paper claims superior distributional metrics compared to Chronos, which \"parrots\" the input. How can mean regression preserve attractor statistics at all? I think this issues deserves closer study. The claims about attractor statistics would be strengthened if, for example, the model was able to reproduce Lyapunov spectra.\n3. **Training details**: While Table 8 reports inference time showing favorable comparison to Chronos, the paper provides no analysis of training setup or timing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "P9Gt8SsjGn", "reviewer": "Reviewer_bbv2", "summary": "This paper introduces  PANDA (Patched Attention for Nonlinear DynAmics), a pretrained transformer model designed to forecast chaotic dynamical systems. The model provides the results of the model to forecast unseen chaotic systems and experiment data snippets in the zero-shot manner; the paper explores the instrinsic properties of Panda in Scaling law, long term statistical property, and some ability in PDE forecasting.", "questions": "1. As the authors claimed 'a framework to create the dataset of 2e4 ODEs' in the contribution 1, the framework would be a great contribution if the authors can clarify a) what algorithmic evolutionary method is applied; b) how to demonstrate the recombined thing is equivalent to physic solvers and simulators? c) What the phyics meanings of the recombined systems and how the community can access and benefit?\n2. The emergent PDE performance is baseline-limited, how the performance is comparable specific surrogate models?\n3. The emergent PDE performance is case-limited, e.g., how the performance is in classic benchmarks of 2D/3D Kolmogrov Flows?\n4. Why Panda uses this set of attention design? What can be expected if other attention types (e.g Fourier Attention, fast attention, attention as gating filters) are applied?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The manuscript shows the following strengths worth credits:\n\n1. the idea of new chaotic trajectories dataset creation is novel;\n2. the emprical visualization of the scaling law Panda Performance;\n3. key indicators in long-term properties preserving and patch-attention maps are included", "weaknesses": "1. Buzz words cause overclaims, such as 'a universal model ...'. The model does not specify whether the model learns from one modality in the system and is able to forecast other modalities in the real-world systems and PDEs.\n2. The emergent PDE performance is less solid, how is the long term statistical property evaluated?\n3. Though the results presentation is extensive, the model layers, param and inference computing details are not found.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)