1 Introduction↩︎

An observed time series can often be viewed as a projection of an underlying dynamical system, though the generative process may be deterministic, stochastic, or a mixture of both. The source of predictive structure shapes which modeling strategies are appropriate. Deterministic chaotic systems are governed by an underlying flow that supports geometry-aware representations, even when long-horizon forecasting is limited by sensitive dependence on initial conditions. Purely stochastic processes are better described through statistical regularities. Both regimes can share low-order moments such as autocorrelation and power spectrum, making the distinction difficult to draw from observations alone. A signal may reflect low-dimensional deterministic dynamics, a stochastic process, or an intermediate regime that retains some structure while losing others.

Classical tools approach this distinction within a single system. Surrogate hypothesis tests [1], [2] ask whether a signal departs from a chosen stochastic null. Lyapunov exponent estimators [3], [4] quantify instability under an assumed deterministic model. Complexity statistics such as permutation entropy [5] and sample entropy [6] compress a signal into a scalar intended to separate deterministic from stochastic behavior. Machine learning approaches have extended work on predictive structure in dynamical systems, including neural differential equation models [7], neural operators [8], [9], and pretrained models for chaotic forecasting [10]. These methods characterize structure within a system or family, and a new signal typically requires new calibration, new surrogates, or a fresh thresholding decision. A detailed comparison is provided in Appendix 7.

This paper addresses a specific problem: given a trajectory window, assign it a position on an ordered continuum of predictability regimes. We use predictability to mean the degree to which the underlying dynamics constrain future trajectory evolution, a property of the governing flow rather than of any particular model or observation process. We instantiate this continuum as a five-level predictability ladder spanning structured deterministic dynamics, weakly chaotic regimes, strongly chaotic regimes, structured stochastic surrogates, and unstructured noise. Standard deterministic-versus-stochastic discrimination corresponds to one boundary in this richer ordinal problem. The central question is whether a learned ordinal score transfers across systems, rather than merely separating classes within a single one.

The ordinal setting introduces a structural ambiguity. If a scalar score correctly orders the ladder levels, any strictly increasing reparameterization preserves the same ordinal decisions after a corresponding threshold adjustment. Correct ordering alone does not determine the numerical coordinate of the score, an identifiability property of ordered regression models [11], leaving scores across systems non-comparable even when both induce identical ordinal predictions. We refer to this as the gauge freedom of ordinal scoring and identify it as the structural barrier to cross-system transfer; post-hoc per-system recalibration cannot substitute, because the intended use case is zero-shot transfer where no target-system data is available to fit a calibration mapping.

To address this, we introduce the Gauge-Fixed Ordinal Network (GON), a temporal convolutional model trained with an anchor-based ordinal objective. The anchor fixes target locations for level-wise score means and penalizes excessive within-level spread, producing a scalar score with a stable coordinate convention. GON operates on 2-jet features (position, velocity, and acceleration), which expose local trajectory geometry [12] preserved by smooth flows and disrupted by surrogate procedures; this representation has been applied to learned physical dynamics in [13]. We train on twelve dissipative chaotic ODE systems, chosen for their bounded attractors and varied geometric structure, and evaluate zero-shot transfer and few-shot adaptation on five fully held-out systems.

The main contributions are as follows.

1. Problem formulation. We formalize ordinal predictability estimation and introduce the predictability ladder, identifying gauge freedom as the structural source of cross-system score ambiguity.

2. Method. We propose GON, which combines 2-jet trajectory features with an anchor-based objective to learn a shared score coordinate across systems.

3. Empirical results. We show that gauge fixing enables cross-system ordinal transfer: pretrained GON outperforms training from scratch across all held-out systems, zero-shot scores retain ordinal coherence at the surrogate boundary, and adaptation depth reflects geometric proximity to the training family.

2 Problem Formulation↩︎

We formalize ordinal predictability estimation: given a trajectory window, assign it a position on an ordered predictability spectrum that remains meaningful across systems, subsuming binary deterministic-versus-stochastic discrimination as a special case.

2.1 The Predictability Ladder↩︎

Let \(x = (x_t)_{t=1}^T \in \mathbb{R}^{d \times T}\) denote a multivariate time series generated by an unknown process. We define a five-level ordered taxonomy \[\mathcal{R} = \{L_0, L_1, L_2, L_3, L_4\}\] called the predictability ladder, arranged from more predictable to less predictable regimes.

This ladder reflects two axes of limited predictability. Along \(L_0\)–\(L_2\), regimes remain deterministic, with predictability degrading as Lyapunov instability shortens the forecast horizon. The \(L_2\!\to\!L_3\) transition instead destroys deterministic geometry: surrogate procedures preserve linear statistics while removing nonlinear flow structure. Level \(L_4\) removes all temporal structure beyond the marginal distribution. Finer granularity can be introduced within adjacent levels by progressively increasing instability (\(L_0\)–\(L_2\)), disrupting deterministic geometry (\(L_2\)–\(L_3\)), or removing residual temporal dependence (\(L_3\)–\(L_4\)). This ordering parallels predictability-horizon arguments in nonlinear dynamical systems, where increasing trajectory divergence reduces forecastability.

2.2 Ordinal Detectability↩︎

Let \(\mu_k = \mathbb{E}[E_\theta(x)\mid y=k]\). To quantify global ordinal structure beyond pairwise separation, we use the monotonicity slope \[\beta = \frac{\sum_k (k-\bar{k})(\mu_k-\bar{\mu})}{\sum_k (k-\bar{k})^2},\] the least-squares slope of level means on ladder indices. A method may achieve a strong AUROC at a single boundary while producing \(\beta \approx 0\) if the full ladder is not coherently ordered. We report \(\beta\) both in-distribution and under transfer to systems withheld entirely from training.

2.3 Cross-System Score Ambiguity↩︎

A scalar ordinal score is useful across systems only if its numerical values retain a stable interpretation; ranking alone does not guarantee this. Given a score \(E:\mathcal{X}\to\mathbb{R}\) and thresholds \(\tau_0 < \tau_1 < \tau_2 < \tau_3\), the induced predictor is \(h_E(x) = \sum_{m=0}^{3}\mathbf{1}\{E(x)>\tau_m\}\), mapping \(x\) to one of five ordered levels.

Ordinal supervision therefore determines score ordering but not its numerical coordinate: the score is identifiable only up to a monotone reparameterization. While this ambiguity is harmless within a single system, it becomes consequential under transfer, where identical score values may correspond to different predictability regimes. Cross-system comparison thus requires a shared coordinate convention, which the anchor objective provides by fixing level-wise score locations globally.

3 Method↩︎

GON has three components: 2-jet preprocessing, a temporal convolutional encoder, and a gauge-fixed ordinal objective.

3.1 2-Jet Preprocessing↩︎

Raw state coordinates carry system-specific scaling, orientation, and units. Representing each window by its 2-jet exposes local trajectory geometry independent of these system-specific factors. For a smooth trajectory \(x(t)\), the 2-jet is the tuple \[(x(t),\dot{x}(t),\ddot{x}(t)),\] consisting of position, velocity, and acceleration. These quantities summarize local behavior up to second order and provide direct access to geometric features such as tangent direction and curvature. For trajectories generated by smooth ODEs, these components are linked by the governing flow. Surrogate procedures, by contrast, preserve selected linear statistics while disrupting the original nonlinear relation among position, velocity, and acceleration. This makes the 2-jet a natural representation for distinguishing deterministic trajectories from structured stochastic surrogates.

3.1.0.1 Implementation.

Given a window \(x\in\mathbb{R}^{d\times T}\), we first apply a Savitzky–Golay filter (window length \(7\), polynomial order \(3\); clamped to signal length and kept odd if shorter) to obtain a smoothed trajectory \(\hat{x}\). Central finite differences are then used to estimate derivatives, as they provide symmetric second-order accurate approximations while minimizing phase bias in the local trajectory geometry: \[\begin{align} \dot{x}_t &= \frac{\hat{x}_{t+1}-\hat{x}_{t-1}}{2\Delta t}, \\ \ddot{x}_t &= \frac{\hat{x}_{t+1}-2\hat{x}_t+\hat{x}_{t-1}}{\Delta t^2}. \end{align}\] One timestep is removed from each boundary to avoid edge effects. The resulting representation is \[J^2(x) = [\hat{x},\,\dot{x},\,\ddot{x}] \in \mathbb{R}^{3d \times (T-2)},\] which serves as the network input.

3.2 Architecture↩︎

The encoder is a temporal convolutional network (TCN) of dilated residual blocks that maps the 2-jet input to a hidden sequence \(h \in \mathbb{R}^{H\times T'}\). After normalization, multi-scale average pooling aggregates \(h\) into a vector \(z \in \mathbb{R}^{4H}\), which a two-layer MLP maps to a scalar score: \[E_\theta(x) = \mathbf{W}_2\,\mathrm{GELU}\!\left(\mathrm{LN}(\mathbf{W}_1 z + b_1)\right) + b_2.\] Outputs are smoothly clipped via \(E_\theta \leftarrow c\,\tanh(E_\theta/c)\), \(c=5\), to keep scores aligned with the anchor range. Full details appear in Appendix 9.

3.3 Gauge-Fixed Ordinal Objective↩︎

Ranking-based supervision does not pin a unique score coordinate (Section 2.3). We fix one by assigning target locations \[\{t_k\}_{k=0}^4 = \{-4,-2,0,+2,+4\}\] and training the model so that each level concentrates around its anchor. These numeric targets are arbitrarily ordered anchors chosen to fix a shared coordinate; any monotonic rescaling yields an equivalent ordinal representation. Equal spacing is chosen for simplicity; the ordinal structure is invariant to any monotone rescaling of the targets, so the specific values affect only the absolute score scale and not relative ordering. The objective is \[\mathcal{L}(\theta) = \lambda_a \sum_{k=0}^{4} (\mu_k - t_k)^2 + \lambda_v \sum_{k=0}^{4} [\sigma_k - \sigma^*]_+^2, \label{eq:loss}\tag{1}\] where \(\mu_k\) and \(\sigma_k\) are the mini-batch mean and standard deviation of scores for level \(k\), \([\cdot]_+=\max(0,\cdot)\), and \(\sigma^*=0.3\). The anchor term fixes level-wise score locations, selecting a canonical representative from the class of monotone-equivalent solutions (Appendix 11); the variance term limits within-level dispersion, together defining a shared affine coordinate across systems.

3.3.0.1 Relation to margin-based objectives.

A natural alternative is margin loss \(\mathcal{L}_{\text{margin}} = \sum_{k=0}^{3}\max(0,\,\gamma-(\mu_{k+1}-\mu_k))\), which enforces adjacent-level separation but leaves the absolute score axis unfixed. Scores therefore remain non-comparable across systems even when levels are well-separated; Section 4.4 tests this directly.

3.4 Training Details↩︎

We set \(\lambda_a=1.0\) with linear warmup over the first five epochs and \(\lambda_v=0.01\). Optimization uses AdamW with learning rate \(2\times10^{-3}\), weight decay \(10^{-4}\), cosine decay over \(30\) epochs, gradient clipping at \(5.0\), and EMA decay \(0.995\). Training uses a batch size of \(64\). Weighted sampling balances ladder levels. The full experimental suite required approximately 40–50 GPU-hours on a single A100.

3.4.0.1 Data Augmentation.

Augmentation is applied to raw trajectory windows \(x \in \mathbb{R}^{d \times T}\) before 2-jet preprocessing and includes additive Gaussian noise (\(\sigma_{\text{rel}}=0.05\)), multiplicative scaling (\(\times[0.9,1.1]\)), temporal masking (\(5\%\)), linear drift (\(\pm 0.05\)), and causal time shifts (\(\leq 5\) steps). Each operation is applied independently with probability \(0.5\); the full pipeline activates with probability \(0.7\). Augmentation precedes smoothing, so the filter attenuates discontinuities introduced by masking or time shifts before finite differences are computed.

4 Experiments↩︎

The experiments address four questions:

1. Does the gauge-fixed coordinate remain calibrated on the source distribution, and does it transfer in zero-shot to held-out systems?

2. Does pretraining provide a better initialization than scratch for adapting to unseen systems under limited labeled windows?

3. Which objective component is responsible for cross-system coordinate stability?

4. Does coherent 2-jet structure contribute beyond dimensionality-matched controls?

4.1 Experimental Setup↩︎

4.1.0.1 Source and target systems.

Training uses 12 three-dimensional dissipative chaotic ODE systems spanning a range of attractor geometries: Chen, Chua, Duffing, Finance, Genesio–Tesi, Hastings–Powell, Lorenz-63, Lorenz-84, Rössler, Rucklidge, Shimizu–Morioka, and Halvorsen. Each system is simulated to produce trajectories of length \(N=4096\), with system-specific timesteps chosen to resolve dominant dynamical timescales; \(\Delta t\) is determined from the observed sampling rate and requires no system-level knowledge.

Ladder levels follow Definition 1: deterministic regimes (\(L_0\)–\(L_2\)) are assigned using \(\lambda_1\), \(L_3\) is generated by applying surrogate procedures (AAFT, IAAFT, WLS, phase-shuffle) to \(L_2\) trajectories, and \(L_4\) is i.i.d.noise. Full simulation and labeling details appear in Appendix 10.

Cross-system transfer is evaluated on five held-out systems: Thomas, Tigan, Newton–Leipnik, Rabinovich–Fabrikant, and the forced pendulum.

4.1.0.2 Windowing and preprocessing.

Trajectories are segmented into windows of length \(T=256\) with stride \(128\). Inputs are normalized via per-trajectory, per-channel z-score standardization, requiring no source-distribution statistics and applying equally to any held-out system. After 2-jet preprocessing, each example is a \(9\times 254\) tensor. Windows are the sole input unit; system-level information is unknown throughout. The adaptation experiment measures how many windows suffice to recalibrate a pretrained model to a new system under limited data.

4.1.0.3 Metrics.

The monotonicity slope \(\beta\) (Definition 2) measures global ordinal coherence. Cross-method zero-shot comparisons use the scale-normalized \(\beta_{\text{norm}}\); raw \(\beta\) is used within the GON family, where the score scale is fixed by construction. Adjacent-pair AUROC at each boundary \((L_0,L_1)\) through \((L_3,L_4)\) serves as a diagnostic, since a method can separate one transition well while failing global coherence. For in-distribution calibration, we also report anchor drift \[\delta_{\text{anchor}}=\frac{1}{5}\sum_{k=0}^4 |\mu_k-t_k|,\] measuring how closely the learned level means match the anchor targets.

4.2 Zero-Shot Evaluation↩︎

Models are trained on the source systems and applied directly to the held-out systems.

4.2.0.1 In-distribution coordinate calibration.

Before assessing transfer, it is useful to verify that the gauge-fixed score remains well aligned with its intended coordinate system on held-out source-domain data. Table 1 reports level-wise score statistics on the test split of the training systems. GON closely matches the anchor targets, with mean anchor drift \(\delta_{\text{anchor}}=0.0097\). The within-level standard deviations (\(0.08\)–\(0.32\)) remain well below the 2-unit spacing between adjacent anchors.

4.2.0.2 Cross-system zero-shot transfer.

Figure 2 compares GON with neural baselines on the held-out systems. All supervised baselines cluster tightly at \(\beta_{\text{norm}} \approx 0.46\)–\(0.48\) despite using different objectives. GON achieves the highest value, \(\beta_{\text{norm}}=0.559\). Full per-boundary AUROC and per-system neural baseline results are reported in Appendix 13. The transfer pattern is selective rather than uniform. The strongest zero-shot separation occurs at the stochastic boundary (\(L_2 \to L_3\) and \(L_3 \to L_4\)), while deterministic-side boundaries remain weaker: \(\text{AUROC}_{01}\) is near chance and \(\text{AUROC}_{12}\) is only moderate (Appendix 13). The gauge-fixed score thus retains meaningful ordinal structure under distribution shift, with transfer concentrated where surrogate procedures most strongly disrupt nonlinear geometry.

4.3 Adaptation from Pretrained Initialization↩︎

We evaluate whether GON pretraining provides a better initialization than training from scratch across \(k \in \{5,10,\ldots,100\}\) labeled windows per held-out system, averaged over five random seeds. Three strategies are compared: scratch (random initialization), pre_head (frozen encoder, readout fine-tuned), and pre_all (full model fine-tuned from the pretrained checkpoint). Table 2 reports averages over all five systems; Figure 3 and Appendix 14 show per-system results and fine-tuning details.

4.3.0.1 Pretrained initialization is faster and less variable than scratch.

At \(k=5\), pre_all leads scratch by \(+0.776\) in \(\beta\) (\(1.319\) vs.\(0.543\)) and \(+0.191\) in \(\text{AUROC}_{23}\) (\(0.758\) vs.\(0.567\)), and this gap does not close at \(k=100\) (Table 2), indicating that the pretrained encoder encodes geometric structure that random initialization cannot recover within this window budget. pre_head reliably exceeds scratch on \(\text{AUROC}_{23}\) but captures little of pre_all’s gain in \(\beta\), confirming that effective adaptation requires updating the encoder.

4.3.0.2 Transfer depth scales with geometric proximity.

Systems with polynomial or exponential vector fields represented in the source distribution (Rabinovich–Fabrikant, Newton–Leipnik, Tigan) adapt immediately or are solved at zero-shot. Rabinovich–Fabrikant is the clearest case: scratch yields \(\beta\approx0.01\)–\(0.54\) across all \(k\), never recovering ordinal structure, whereas pre_all exceeds zero-shot from \(k=10\). Thomas and the forced pendulum are harder: both contain state-dependent sinusoidal vector-field terms absent from training (e.g.,\(\sin\theta\), \(\sin y\)), producing unseen oscillatory curvature in the 2-jet. In both cases, pre_all recovers global ordinal structure but fails to resolve the \(L_2\to L_3\) boundary, identifying sinusoidal vector-field geometry as the primary gap in source-distribution coverage.

4.4 Ablations↩︎

The ablations ask which parts of the method matter for transfer.

4.4.0.1 Objective ablation: why the anchor matters.

Table 3 shows that coordinate stability, not ordinal coherence alone, is critical for transfer. Replacing the anchor with a margin objective leaves in-distribution AUROC unchanged by at most \(0.010\) across adjacent pairs and produces comparable scale-normalized transfer (\(\beta_\text{norm}=0.554\) vs.\(0.559\)), yet anchor drift rises from \(0.005\) to \(2.556\). A model can therefore retain ordinal coherence while losing coordinate stability, confirming that relative level separation is insufficient for cross-system comparability: without a fixed coordinate, scores on unseen systems carry no consistent numerical meaning. The variance term plays a secondary but consistent role: removing it increases within-level spread from \(0.195\) to \(0.230\) and lowers transfer \(\beta\) from \(1.118\) to \(1.055\).

4.4.0.2 Input ablation: why the 2-jet helps.

Table 4 shows that coherent derivative structure improves transfer beyond raw states. Adding first derivatives improves \(\beta\) from \(1.132\) to \(1.336\); the full 2-jet improves further to \(\beta=1.504\) and \(\text{AUROC}_{23}=0.846\). Shuffling derivatives across time drops \(\beta\) to \(0.983\), confirming that temporal coherence among position, velocity, and acceleration drives the gain. A control replacing derivatives with per-pass Gaussian noise remains near the raw-state baseline, ruling out channel count.

4.4.0.3 Noise robustness.

GON remains robust under moderate corruption, with \(\text{AUROC}_{23}\) staying above \(0.89\) at \({\geq}\,20\) dB SNR and \(30\%\) temporal dropout. Performance degrades gracefully in \(\beta\), indicating reduced transfer strength before loss of discriminability. Full results appear in Appendix 15.

5 Discussion↩︎

5.0.0.1 A fixed score coordinate enables transferable predictability estimates.

Replacing the anchor term with a margin loss leaves in-distribution discrimination nearly unchanged but increases anchor drift and reduces transfer \(\beta\) by \(47\%\) (Table 3). Source-domain separability alone, therefore, does not ensure that a scalar score retains a stable meaning on unseen systems. Gauge fixing instead learns a shared ordinal coordinate, enabling predictability comparisons that are transferred across heterogeneous dynamics.

5.0.0.2 Transfer concentrates where predictability is structurally limited.

Zero-shot transfer is strongest at the stochastic boundary (\(L_2 \to L_3\) and \(L_3 \to L_4\)); deterministic-side boundaries remain weaker. The same asymmetry holds under perturbation: \(\text{AUROC}_{23}\) is more robust than \(\text{AUROC}_{01}\) to noise and dropout. GON captures geometric structure disrupted by surrogate generation, whereas finer distinctions within deterministic chaos are more system-specific.

5.0.0.3 Geometric representations support cross-system generalization.

The 2-jet representation improves transfer by exposing local trajectory geometry. The benefit depends on coherent relationships among position, velocity, and acceleration; disrupting this alignment reduces the transfer even at a matched dimensionality. This supports the view that transferable predictability cues arise from shared geometric structure rather than system-specific statistics.

5.0.0.4 Pretraining helps, but transfer depth is system-dependent.

Systems whose vector fields resemble those seen in training adapt immediately or are solved at zero-shot. In contrast, systems with unseen structure (e.g., sinusoidal forcing in Thomas and the forced pendulum) require more data and may fail to resolve finer boundaries. Performance degrades smoothly with geometric distance from the source distribution.

5.0.0.5 Scope, limitations, and broader implications.

The benchmark covers 12 source and 5 target systems, all low-dimensional and dissipative; extensions to higher-dimensional systems, conservative dynamics, or real-world signals, remain unestablished. This work takes one of the earliest steps toward a general theory of transferable predictability estimation, establishing gauge-fixing as a necessary condition for cross-system ordinal comparability. Predictive structure in real systems is often limited by chaos, stochastic forcing, or their interaction, arising across weather and climate dynamics, turbulent flows, biological oscillations, financial processes, and engineered systems. By learning a transferable ordinal predictability coordinate, the framework enables cross-system comparison without system-specific calibration, informing model selection, early-warning diagnostics, and the design of resilient coupled human–natural systems. The method operates exclusively on synthetic ODE signals, involves no personal data, and has no foreseeable negative societal impacts.

6 Conclusion↩︎

This paper identifies gauge-fixing as the missing ingredient for transferable ordinal predictability scores. The central challenge is not ranking regimes within a single system, but assigning scores whose numerical values retain a consistent interpretation across systems. The loss ablation makes the point cleanly: replacing the anchor with a margin objective leaves the source-domain discrimination nearly unchanged while reducing transfer \(\beta\) by \(47\%\). Ordinal coherence is robust to 4-bit quantization and Gaussian noise down to 20 dB SNR, confirming the learned coordinate survives common observational degradations. Pairwise discrimination and globally coherent ordinal scoring are distinct properties that can dissociate, and a stable coordinate is necessary to bridge them under transfer. The benchmark is limited to low-dimensional dissipative systems; this bound does not invalidate the core finding, since the ablation isolates coordinate convention as the operative variable within a fixed system class. Extension to higher-dimensional attractors, partial observability, and real-world signals in weather, climate, and physiology remain the most consequential direction for future work.

Appendix↩︎

7 Related Work↩︎

7.0.0.1 Determinism versus stochasticity.

Geometric diagnostics [14], [15] and prediction-based tests study rejection of the stochastic null hypothesis within a single system. Permutation entropy [5] and sample entropy [6] provide scalar complexity measures widely used as binary indicators of dynamical complexity. The present work places this binary separation within a broader ordinal hierarchy, where pairwise discriminability at a single boundary and global ordinal coherence across levels need not coincide.

7.0.0.2 Chaos detection and surrogate testing.

The classical toolbox for the binary problem is surrogate hypothesis testing [1], [2], which constructs a null distribution from the signal and tests whether its nonlinear statistics exceed those of matched stochastic alternatives. Complementary approaches estimate the maximal Lyapunov exponent directly [3], [4]. Both are designed for single-system analysis: each new signal requires new surrogates or freshly tuned estimators, and both address a binary question. The present work considers an ordinal formulation in which predictability is represented on a continuum, with scores intended to remain comparable across systems without per-signal calibration.

7.0.0.3 Early warning signals.

Critical slowing down theory [16] and spectral precursor methods [17] study statistical signatures of approaching bifurcations. The predictability ladder subsumes early-warning detection as a special case: \(L_1\) is precisely the transition regime where such signatures emerge. The present benchmark does not directly evaluate early-warning performance. Transfer is strongest at the stochastic boundary rather than the deterministic-side transition, and calibrated cross-system early-warning detection remains a natural extension.

7.0.0.4 Ordinal learning.

Cumulative link models [18] and neural extensions such as CORAL [19] enforce a consistent rank structure across ordered labels. Ordinal contrastive methods [20] encode ordering in representation space. Proxy-based metric learning [21] uses class-level anchor prototypes to stabilize embedding geometry. These formulations study ordering within a single system or dataset, while the setting considered here involves scores that are intended to remain comparable across heterogeneous systems.

7.0.0.5 Physics-informed and geometric representations.

Neural differential equation models [7] and Lagrangian networks [13] encode geometric and physics-aware inductive biases in learned representations. Reservoir computing has been applied to chaotic forecasting [22]. The 2-jet representation used here is grounded in differential geometry: position, velocity, and acceleration jointly encode local trajectory curvature in a way that is invariant to system-specific scaling and is disrupted by surrogate generation procedures.

7.0.0.6 Transfer and meta-learning for dynamical systems.

Gradient-based meta-learning [23] has been applied to dynamical systems forecasting [24]. Neural operators [8], [9] learn solution operators across parameter families. [10] trained a transformer on a large synthetic corpus of chaotic systems and reported zero-shot forecasting on unseen ODEs and PDEs. These approaches study predictive generalization across systems under supervised forecasting objectives. The setting considered here focuses on learning representations of predictability regimes, where the goal is to assign comparable ordinal scores rather than to predict future states.

8 Proof of Proposition 1↩︎

Let \(f:\mathbb{R}\to\mathbb{R}\) be strictly increasing. For any threshold \(\tau_m\) and input \(x\), \[E(x) > \tau_m \iff f(E(x)) > f(\tau_m),\] since strictly increasing functions preserve the direction of all inequalities. The induced ordinal predictor \(h_E(x) = \sum_{m=0}^{3}\mathbf{1}\{E(x) > \tau_m\}\) depends only on these threshold comparisons. Under the transformed score \(\tilde{E} = f \circ E\) with transformed thresholds \(\tilde{\tau}_m = f(\tau_m)\), each comparison \(\tilde{E}(x) > \tilde{\tau}_m\) is equivalent to \(E(x) > \tau_m\). Therefore, \(h_{\tilde{E}}(x) = h_E(x)\) for all \(x\), and the induced ordinal prediction is unchanged. 0◻

9 Architecture Details↩︎

9.1 Temporal Convolutional Encoder↩︎

The encoder is a temporal convolutional network (TCN) composed of \(B = 6\) residual blocks with exponentially increasing dilation factors \([1, 2, 4, 8, 16, 32]\). Each block consists of two 1D dilated convolutions with kernel size \(k = 3\), followed by Group Normalization and GELU activation. Residual connections are applied at the block level.

The hidden dimension is fixed at \(H = 128\) across all layers.

9.2 Receptive Field↩︎

The effective receptive field of the network is \[\mathcal{R} = 1 + 2(k - 1)\sum_{i=0}^{B-1} 2^i = 253,\] which covers nearly the full temporal extent of the input after 2-jet boundary trimming (\(T' = T - 2\)). This allows the model to capture long-range temporal dependencies relevant for distinguishing predictability regimes.

9.3 Readout and Multi-Scale Aggregation↩︎

The encoder produces a feature map \(h \in \mathbb{R}^{H \times T'}\). RMS normalization is applied before aggregation. Average pooling is performed at scales \(s \in \{4, 16, 64, 256\}\) to capture structure at multiple temporal resolutions; for \(s > T'\) the kernel is clamped to \(k = \min(s, T')\), so at scale \(256\) this reduces to a global average over all \(254\) timesteps. Pooling uses no padding; non-integer-multiple tails are truncated by F.avg_pool1d with default padding=0. The pooled representations are concatenated to form a feature vector \(z \in \mathbb{R}^{4H}\).

9.4 MLP Head and Output Scaling↩︎

The pooled feature vector \(z \in \mathbb{R}^{4H}\) (dimension \(512\)) is mapped to a scalar score via a two-layer MLP with hidden dimension \(H=128\): \[E_\theta(x) = \mathbf{W}_2\,\mathrm{GELU}\!\left(\mathrm{LN}(\mathbf{W}_1 z + b_1)\right) + b_2,\] where \(\mathbf{W}_1 \in \mathbb{R}^{128 \times 512}\) and \(\mathbf{W}_2 \in \mathbb{R}^{1 \times 128}\). To prevent unbounded growth and maintain consistency with the fixed anchor targets, the output is smoothly clipped: \[E_\theta \leftarrow c \tanh(E_\theta / c), \quad c = 5.\] The full model contains \(616{,}233\) trainable parameters.

10 Data Generation and Labeling Details↩︎

10.1 Numerical Integration↩︎

Each system is simulated to produce \(N = 4096\) uniformly sampled time steps. The timestep \(\Delta t\) is selected to resolve the dominant dynamical timescale of each system. We estimate a characteristic timescale \(\tau\) from the mean inter-peak interval of the trajectory, and set \(\Delta t\) such that approximately 64 dominant cycles are covered, with \(\Delta t\) clamped to \([10^{-4}, 10^{-1}]\).

Integration is performed using either a fixed-step fourth-order Runge–Kutta method or the adaptive Dormand–Prince (DOPRI5) solver, depending on system stiffness.

10.2 Predictability Ladder Construction↩︎

Deterministic regimes follow Definition 1, with \(\lambda_1\) estimated via [4]. \(L_0\) corresponds to \(\lambda_1 < 0\); \(L_1\) and \(L_2\) are distinguished by the relative magnitude of \(\lambda_1\), cross-checked by visual inspection of attractor phase portraits. Because attractor geometries differ substantially across systems, no single numerical threshold is applied uniformly; regime boundaries are confirmed by qualitative separation in both the exponent estimate and the phase portrait. These boundaries are fixed at data generation time and held constant throughout all experiments; the qualitative inspection step affects only this initial label assignment and does not enter any learned component of the method. \(L_3\) is generated from \(L_2\) trajectories using AAFT, IAAFT, WLS, and phase randomization [1], [2], [25], [26]. \(L_4\) is i.i.d.Gaussian noise matched in length and dimensionality.

10.3 Data Splits↩︎

Train/validation/test splits are performed at the trajectory level using an 80/10/10 partition. All normalization statistics are computed on the training split and held fixed for validation and testing.

11 Gauge Fixing and Ordinal Identifiability↩︎

Ordinal regression determines only the ordering of latent scores. For any strictly increasing \(f\), the transformed score \(\tilde{E}(x)=f(E(x))\) induces identical ordinal decisions under the corresponding threshold transformation, so \(E\) is identifiable only up to a monotone reparameterization [11]. Within a single system, this ambiguity is harmless: thresholds are learned jointly with the score, and ranking suffices. Under cross-system transfer, it becomes consequential because identical ordinal predictions on two systems need not assign the same numeric value to the same predictability regime.

The anchor objective removes this degree of freedom by fixing the level-wise score means, \[\mathbb{E}[E(x)\mid y=k] = t_k,\] standardizing the macroscopic location and scale of the score distribution, breaking the unbounded gauge freedom of ordinal regression. The variance term further restricts within-level dispersion. Together, they establish a shared coordinate convention across systems without claiming strict uniqueness: infinitely many nonlinear strictly increasing functions could, in principle, satisfy these moment constraints, but the anchor fixes the operationally relevant degrees of freedom for cross-system comparison.

This differs from post-hoc calibration, which rescales scores after training on a per-system basis and, therefore, cannot produce a universal numeric coordinate. Gauge fixing constrains the representation during learning, so the coordinate generalizes by construction.

12 Benchmark System Equations and Parameters↩︎

For each system, we report the governing equations, fixed parameter values, the single parameter varied to traverse the predictability ladder, and the parameter range assigned to each regime. Regime boundaries were determined by the protocol described in Appendix 10: the sign and relative magnitude of the leading Lyapunov exponent \(\lambda_1\) cross-checked by visual inspection of attractor phase portraits. Initial conditions were sampled from multiple distributions per system (Gaussian perturbations about a base point, uniform on a sphere, and system-specific defaults); all distributions produced qualitatively consistent regime behavior.

The stable range corresponds to L0, the transition range to L1 (weakly chaotic), and the chaotic range to L2 (strongly chaotic).

Unless noted otherwise, initial conditions were sampled from three distributions: Gaussian perturbations about a base point, uniform on a sphere, and system-specific defaults. All distributions produced qualitatively consistent regime behavior. Two exceptions are noted explicitly: Hastings–Powell, which requires ecologically structured non-negative initial conditions, and Chua, which additionally uses double-scroll sampling centered on the two attractor scrolls.

Source Systems (Training Set)↩︎

Chen [27]↩︎

\[\begin{align} \dot{x} &= a(y - x) \\ \dot{y} &= (c - a)x - xz + cy \\ \dot{z} &= xy - bz \end{align}\]

Fixed parameters: \(a = 35\), \(b = 3\). Varied parameter: \(c\).

Chua [28]↩︎

\[\begin{align} \dot{x} &= \alpha\bigl(y - x - f(x)\bigr) \\ \dot{y} &= x - y + z \\ \dot{z} &= -\beta y \end{align}\]

where the piecewise-linear Chua diode is \[f(x) = m_1 x + \tfrac{1}{2}(m_0 - m_1)\bigl(|x+1| - |x-1|\bigr).\]

Fixed parameters: \(\beta = 28.0\), \(m_0 = -1.143\), \(m_1 = -0.714\). Varied parameter: \(\alpha\).

Chua additionally uses double-scroll sampling (\(\sigma = 0.3\) about \((\pm 1.5,\;0,\;0)\)) and Gaussian perturbations (\(\sigma = 0.01\)) about the base point \((0.7,\;0,\;0)\).

Duffing [29]↩︎

The Duffing oscillator is non-autonomous; the cosine forcing term renders the effective phase space three-dimensional.

\[\begin{align} \dot{x} &= y \\ \dot{y} &= -\delta y - \alpha x - \beta x^3 + \gamma\cos(\omega t) \end{align}\]

Fixed parameters: \(\delta = 0.1\), \(\alpha = -1.0\), \(\beta = 1.0\), \(\omega = 1.0\). Varied parameter: \(\gamma\) (forcing amplitude).

Finance [30]↩︎

\[\begin{align} \dot{x} &= \left(\tfrac{1}{b} - a\right)x + z + xy \\ \dot{y} &= -by - x^2 \\ \dot{z} &= -x - cz \end{align}\]

Fixed parameters: \(b = 0.2\), \(c = 1.1\). Varied parameter: \(a\).

Genesio–Tesi [31]↩︎

\[\begin{align} \dot{x} &= y \\ \dot{y} &= z \\ \dot{z} &= -cx - by - az + x^2 \end{align}\]

Fixed parameters: \(b = 1.1\), \(c = 1.0\). Varied parameter: \(a\).

Halvorsen [32]↩︎

\[\begin{align} \dot{x} &= -ax - by - bz - y^2 \\ \dot{y} &= -ay - bz - bx - z^2 \\ \dot{z} &= -az - bx - by - x^2 \end{align}\]

Fixed parameter: \(b = 4.0\). Varied parameter: \(a\).

Hastings–Powell [33]↩︎

\[\begin{align} \dot{x} &= x(1-x) - \frac{a_1 xy}{1 + b_1 x} \\ \dot{y} &= \frac{a_1 xy}{1 + b_1 x} - \frac{a_2 yz}{1 + b_2 y} - d_1 y \\ \dot{z} &= \frac{a_2 yz}{1 + b_2 y} - d_2 z \end{align}\]

Fixed parameters: \(a_2 = 0.1\), \(b_1 = 3.0\), \(b_2 = 2.0\), \(d_1 = 0.4\), \(d_2 = 0.01\). Varied parameter: \(a_1\) (prey–predator interaction rate).

Initial conditions were absolute-valued to ensure non-negative populations, with ecological structured sampling (prey \(\in [0.5,1.0]\), predator \(\in [0.1,0.4]\), superpredator \(\in [0.2,0.8]\)) used in place of sphere sampling. Integration used an adaptive timescale with an \(8\times\) ecological multiplier applied to the estimated dominant period \(\tau\), reflecting the longer characteristic timescales of trophic oscillations.

Lorenz-63 [34]↩︎

\[\begin{align} \dot{x} &= \sigma(y - x) \\ \dot{y} &= x(\rho - z) - y \\ \dot{z} &= xy - \beta z \end{align}\]

Fixed parameters: \(\sigma = 10.0\), \(\beta = 8/3\). Varied parameter: \(\rho\).

Lorenz-84 [35]↩︎

\[\begin{align} \dot{x} &= -ax - y^2 - z^2 + aF \\ \dot{y} &= -y + xy - bxz + G \\ \dot{z} &= -z + bxy + xz \end{align}\]

Fixed parameters: \(a = 0.25\), \(b = 4.0\), \(G = 1.0\). Varied parameter: \(F\) (forcing).

Rössler [36]↩︎

\[\begin{align} \dot{x} &= -y - z \\ \dot{y} &= x + ay \\ \dot{z} &= b + z(x - c) \end{align}\]

Fixed parameters: \(a = 0.2\), \(b = 0.2\). Varied parameter: \(c\).

Rucklidge [37]↩︎

\[\begin{align} \dot{x} &= -ax + by - yz \\ \dot{y} &= x \\ \dot{z} &= -z + y^2 \end{align}\]

Fixed parameter: \(a = 2.0\). Varied parameter: \(b\).

Shimizu–Morioka [38]↩︎

\[\begin{align} \dot{x} &= y \\ \dot{y} &= x - ay - xz \\ \dot{z} &= -bz + x^2 \end{align}\]

Fixed parameter: \(b = 0.5\). Varied parameter: \(a\).

Target Systems (Held-Out Set)↩︎

Forced Pendulum [39]↩︎

Non-autonomous; cosine forcing renders the effective phase space three-dimensional.

\[\begin{align} \dot{\theta} &= v \\ \dot{v} &= -\gamma v - \sin\theta + A\cos(\omega t) \end{align}\]

Fixed parameters: \(\gamma = 0.2\), \(\omega = 2/3\). Varied parameter: \(A\) (forcing amplitude).

Newton–Leipnik [40]↩︎

\[\begin{align} \dot{x} &= -ax + y + 10yz \\ \dot{y} &= -x - 0.4y + 5xz \\ \dot{z} &= bz - 5xy \end{align}\]

Fixed parameter: \(a = 0.4\). Varied parameter: \(b\) (z-axis damping).

Rabinovich–Fabrikant [41]↩︎

\[\begin{align} \dot{x} &= y(z - 1 + x^2) + \gamma x \\ \dot{y} &= x(3z + 1 - x^2) + \gamma y \\ \dot{z} &= -2z(\alpha + xy) \end{align}\]

Fixed parameter: \(\alpha = 0.1\). Varied parameter: \(\gamma\).

Trajectories exceeding \(\|s\|_\infty = 25\) were truncated to the bounded portion before divergence; trajectories with fewer than 300 retained points were discarded.

Thomas [42]↩︎

\[\begin{align} \dot{x} &= \sin y - bx \\ \dot{y} &= \sin z - by \\ \dot{z} &= \sin x - bz \end{align}\]

No fixed parameters. Varied parameter: \(b\) (dissipation).

Tigan [43]↩︎

\[\begin{align} \dot{x} &= a(y - x) \\ \dot{y} &= (c - a)x - axz \\ \dot{z} &= -bz + xy \end{align}\]

Fixed parameters: \(a = 2.1\), \(b = 0.6\). Varied parameter: \(c\).

13 Neural Baseline Comparison↩︎

We compare GON with three neural baselines trained on the same 12 source systems using identical data, preprocessing, encoders, and optimization protocols. The only change is the training objective.

Regression minimizes mean squared error to integer labels \(\{0,1,2,3,4\}\). Classification minimizes cross-entropy over five classes. CORAL [19] uses cumulative link constraints for ordinal supervision.

Each method is converted to a scalar score for evaluation: the raw output for regression, the expected label \(\sum_k k\,p_k\) for classification, and the cumulative-logit score for CORAL.

13.0.0.1 In-distribution.

Table 5 reports test performance on the 12 source systems. Regression, classification, CORAL, and GON all achieve near-perfect adjacent-pair discrimination in-distribution. To make cross-method comparisons meaningful despite different score scales, we report the scale-normalized monotonicity slope \(\beta_{\text{norm}}\) rather than raw \(\beta\). Under this normalization, regression, classification, CORAL, and GON all show nearly perfect global ordering on the source distribution. This pattern indicates that source-domain performance alone is insufficient to reveal whether a learned score will transfer with a stable numerical meaning.

13.0.0.2 Zero-shot transfer.

Table 6 reports zero-shot performance on the 5 held-out target systems. All supervised baselines cluster at similar values of \(\beta_{\text{norm}} \approx 0.46\)–\(0.48\) despite their different objectives, suggesting a shared limitation: none learns a score with a stable cross-system coordinate.

GON achieves the highest scale-normalized slope, \(\beta_{\text{norm}}=0.559\), while also remaining competitive on the strongest transferable boundary \(L_2 \rightarrow L_3\). This is the relevant comparison for cross-system ordinal structure, since raw \(\beta\) is not directly comparable across methods with different score scales.

14 Per-System Adaptation Results↩︎

Tables 7–16 report monotonicity slope \(\beta\) and all four adjacent-pair AUROC values for each of the five held-out systems, across \(k \in \{5,10,\ldots,100\}\) labeled windows. All entries are mean \(\pm\)  std over 5 random seeds. Zero-shot (ZS) performance is listed in the caption of each system block for reference. Columns correspond to the three adaptation strategies: scratch (random initialization), pre_head (frozen encoder, readout fine-tuned), and pre_all (full model fine-tuned from the pretrained checkpoint).

14.0.0.1 Fine-tuning protocol.

All three adaptation strategies use AdamW with learning rate \(2\times10^{-3}\), weight decay \(10^{-4}\), cosine annealing over 30 epochs, and gradient clipping at 5.0; identical to pretraining. Data augmentation and EMA are disabled during fine-tuning. For pre_head, all encoder parameters are frozen including GroupNorm weight/bias; only the readout, normalization layer, and scale/bias embeddings are updated. The trajectory-level train/val/test split is fixed across all runs (seed 42). For each \((k, \text{seed})\) pair, \(k\) windows are drawn without replacement from the fixed training pool, with at least one window per ladder level before filling the remaining budget uniformly at random.

Newton–Leipnik↩︎

Zero-shot: \(\beta=0.737\), \(\text{AUROC}_{01}=0.539\), \(\text{AUROC}_{12}=0.514\), \(\text{AUROC}_{23}=0.856\), \(\text{AUROC}_{34}=1.000\).

Forced Pendulum↩︎

Zero-shot: \(\beta=0.489\), \(\text{AUROC}_{01}=0.295\), \(\text{AUROC}_{12}=0.390\), \(\text{AUROC}_{23}=0.604\), \(\text{AUROC}_{34}=1.000\).

Rabinovich–Fabrikant↩︎

Zero-shot: \(\beta=1.270\), \(\text{AUROC}_{01}=0.483\), \(\text{AUROC}_{12}=0.576\), \(\text{AUROC}_{23}=0.725\), \(\text{AUROC}_{34}=0.872\).

Thomas↩︎

Zero-shot: \(\beta=0.958\), \(\text{AUROC}_{01}=0.508\), \(\text{AUROC}_{12}=0.452\), \(\text{AUROC}_{23}=0.391\), \(\text{AUROC}_{34}=0.995\).

Tigan↩︎

Zero-shot: \(\beta=2.105\), \(\text{AUROC}_{01}=0.645\), \(\text{AUROC}_{12}=0.985\), \(\text{AUROC}_{23}=0.997\), \(\text{AUROC}_{34}=1.000\).

15 Full Noise Robustness Results↩︎

Table 17 reports the complete inference-time perturbation sweep on the 12 source systems (in-distribution). No retraining is performed; the pretrained GON checkpoint is applied directly.

15.0.0.1 Perturbation details.

All perturbations are applied at inference only under torch.no_grad(). Gaussian noise: SNR is defined on a mean-square power basis, \(\mathrm{SNR}_\mathrm{dB} = 10\log_{10}(P_s/P_n)\) where \(P_s = \mathrm{mean}(x^2)\) over all channels and timesteps jointly per sample; a single noise power scalar is shared across channels. Temporal dropout: selected timesteps are zero-masked in-place (sequence length preserved), with the same mask applied across all channels simultaneously; each sample draws an independent mask. Quantization: post-training uniform asymmetric quantization applied per channel per sample (min-max over the time axis) in pure PyTorch; zero is a representable level (uniform mid-tread).

15.0.0.2 Results.

(1) Ordinal coherence is preserved down to 20 dB SNR and falls below \(1.0\) between 30 and 20 dB; at 5 dB, \(\beta=0.216\), but the score remains positively monotone on average. (2) Degradation is boundary-asymmetric throughout: \(\text{AUROC}_{34}\) holds at \(1.000\) through 30 dB Gaussian noise and all dropout fractions up to 30%, while \(\text{AUROC}_{01}\) falls below \(0.5\) at 20 dB and at 20% dropout; \(\text{AUROC}_{12}\) is intermediate, tracking \(\text{AUROC}_{01}\) in shape but degrading more slowly. (3) Quantization is negligible at all tested precisions: no metric changes by more than \(0.002\) between 16-bit and 4-bit, and \(\beta\) is unchanged to three decimal places down to 6-bit.

References↩︎