Extrapolating the emergence of Hamiltonian chaos with random-feature Hamiltonian neural networks

Near-integrable Hamiltonian systems exhibit one of the most thoroughly studied routes to chaos in nonlinear dynamics. As a perturbation or control parameter grows, Kolmogorov–Arnold–Moser (KAM) tori are progressively destroyed, resonance layers widen and overlap, and an initially thin chaotic layer grows into a sea that dominates the energy shell [7, 22]. The resulting mixed phase space, in which islands of regular motion coexist with chaotic regions, is generic for conservative systems ranging from galactic potentials to molecular vibrations [14, 2, 21]. On a Poincaré section, the progression appears as the gradual breakup of invariant curves and the spread of scattered points over the accessible area, as exemplified by the Hénon–Heiles system [14]. A one-parameter family of Hamiltonians therefore spans the full range of behavior from predominantly regular to strongly chaotic motion.

Reconstructing such dynamics from data becomes necessary when the governing Hamiltonian is not known in explicit form. Conventional neural networks handle this task poorly, as a model trained directly on trajectory data does not inherit the conservation laws of the underlying flow, and its long-time predictions acquire spurious dissipation or energy drift even when training data are abundant [12]. Hamiltonian neural networks (HNNs) address this failure by representing a scalar function $H_{\theta}$ and obtaining the vector field from Hamilton’s equations, so that the learned flow conserves the learned energy and phase-space volume by construction [12]. Symplectic architectures take a complementary route and learn the discrete-time flow map itself, either by composing elementary symplectic modules or through parameterizations that are symplectic by construction [17, 6]. Structure-preserving models of this kind have delivered stable long-time predictions for systems ranging from simple oscillators to many-body problems and have conserved energy over horizons where conventional networks drift [12, 17, 6]. Trained across the order-to-chaos transition of the Hénon–Heiles system, HNNs also reproduce both regular and chaotic behavior, at conditions covered by the training data [10]. The generalization has been pushed even further by adaptable HNNs, which take the control parameter as an additional input and have reconstructed the Hamiltonian dynamics of a family at control values absent from the training set [13]. In these demonstrations, however, the training data already contained chaotic orbits, so predictions at new control values, including values outside the training set, reproduced dynamical regimes that were represented in training rather than a qualitatively new one.

Extrapolation beyond the training range is a far more demanding task [11], but a series of recent results in reservoir computing has shown that learned models can cross dynamical regimes [18, 20, 8]. In this approach, a fixed, randomly generated recurrent network is driven by the observed time series, and its internal states serve as random features on which only a linear readout is trained [15, 9]. Parameter-aware reservoir computers that receive the control as an additional input have reconstructed attractors at controls not sampled in training [18], predicted critical transitions and system collapse from pre-transition data [20, 24], and extrapolated tipping points of non-stationary systems when internal hyperparameters such as the spectral radius are optimized [19]. Regime-crossing results continue to accumulate in other settings: parameter-dependent neural ordinary differential equations extrapolate across bifurcations of dissipative systems [28], and an evolution-operator emulator trained only on chaotic dynamics predicts an unseen laminar regime [27]. Extrapolation of Hamiltonian dynamics, however, remains largely unexplored. Reservoir computing has been applied to Hamiltonian systems, reconstructing KAM diagrams from trajectories recorded at a few parameter labels [29]. In the mixed-regime demonstrations, however, chaotic motions were part of the training data, and the authors note that quasi-periodic training motions alone do not suffice to replicate the chaotic dynamics. The framework also offers no mechanism that enforces Hamiltonian structure by construction. Yet the random-feature strategy underlying reservoir computing is not tied to recurrent networks. Random-feature Hamiltonian models replace the reservoir with a static set of nonlinear features, drawn at random or sampled from data, and use the readout to represent the Hamiltonian itself, so that the fit reduces to a ridge regression solved in closed form [16, 5, 26, 25]. Models of this type retain the conservation structure of an HNN while avoiding back-propagation entirely, and for a single training system, data-driven sampling schemes such as SWIM have reached accuracies comparable to those of back-propagated networks at a fraction of the training cost [5, 26, 25].

In this study, we apply this random-feature construction to the extrapolation of a Hamiltonian transition. Instead of drawing the features once and using them as given, we optimize internal hyperparameters such as the feature scales, as in reservoir computing, and the fast training provided by the single closed-form ridge fit is what makes this optimization affordable. We show that the resulting parameter-aware random-feature HNN (RF-HNN; Fig. 1), trained only at control values where regular motion dominates, extrapolates the emergence of widespread chaos beyond the training range. The extrapolation regime contributes nothing to training or model selection, and the test is carried out on four two-degree-of-freedom families, namely two Hénon–Heiles families, a confined Barbanis-type model, and a bilinearly coupled Morse pair [14, 2, 23], against conventionally trained parameter-aware HNNs [12, 13], both as a single network and as a 20-network ensemble. Because every model in the comparison is Hamiltonian by construction, differences in the extrapolated dynamics can be attributed to the fitted continuation in the control rather than to the presence of conservation structure. In all four families, the RF-HNN tracks the growth of the chaotic sea quantitatively, whereas the conventionally trained models, including the ensemble, fail to track that growth and remain too regular at the largest extrapolation controls. To our knowledge, no learned Hamiltonian has previously crossed this transition from a regular-dominated training regime, and the comparison identifies parameter continuation, rather than conservation structure alone, as the decisive factor.

We consider one-parameter families of two-degree-of-freedom Hamiltonians $H_{c}(z)$ with $z=(x,y,p_{x},p_{y})$, in which increasing the control $c$ drives the standard near-integrable route to chaos. We write $H(z,c)\equiv H_{c}(z)$ when the control appears as an explicit argument. Along this route, invariant tori break up, resonance layers overlap, and a chaotic sea spreads over the energy shell. Figure 1(a) illustrates the extrapolation problem posed in this work for the canonical Hénon–Heiles family of Sec. II.1. At the observed control, the Poincaré section is organized by invariant curves and the orbit-wise finite-time Lyapunov exponents remain uniformly small. At the unseen target control, the same sampling yields a section dominated by a chaotic sea. The parameter axis beneath the sections summarizes the data constraints of the task. The learner receives vector-field samples $f_{i}=J\nabla_{z}H(z_{i},c_{i})$, where $J$ is the canonical symplectic matrix, on a fixed energy shell $E_{0}$, drawn only at the filled training controls, at which the section dynamics is dominated by invariant tori. Model selection, including all hyperparameter tuning, is likewise restricted to controls at or inside the training interval. The trained model is then integrated autonomously at the open test control, which lies strictly above the training interval and at which the true family develops a broad mixed phase space. Success is judged not by short-time trajectory error but by the long-time dynamical content of the learned flow, namely the occupied area and geometry of its Poincaré sections and its orbit-wise finite-time Lyapunov exponents, compared with truth by identical procedures (Secs. III.1 and III.2). No data from the extrapolation regime are used for training or model selection at any stage.

The two panels of Fig. 1(a) make clear that this task differs from interpolation in kind rather than in degree. Within the training interval, a faithful model needs only to connect regular behavior of the kind it has already observed. Beyond this interval, it must produce a qualitatively new form of long-time behavior from a fitted Hamiltonian whose continuation in $c$ was never supervised. This is the transition from the left panel of Fig. 1(a) to the right. Every model considered below is Hamiltonian by construction, but nothing in that constraint determines how the fitted Hamiltonian behaves at controls beyond the training interval. A model can conserve its learned energy exactly and still predict regular motion where the true family has become chaotic.

We instantiate the task on four families, summarized in Table 1. Two Hénon–Heiles families carry the principal claims, and two auxiliary families test whether the results generalize. We report a dedicated interpolation test for the bounded family. Every control designated for extrapolation lies strictly above its family’s training interval.

Let $c$ denote the control parameter of a given family. Our model, the parameter-aware RF-HNN shown in Fig. 1(b), represents the Hamiltonian of the family as

where $\sigma$ is an elementwise activation function and the state weights $W_{z}$, the control weights $w_{c}$, the second-layer matrix $W$, and the biases $b_{1}$ and $b_{2}$ are sampled once and kept fixed. Only the readout $\beta$ is trained. For a training state $z_{i}$ with target field $f_{i}=J\nabla_{z}H(z_{i},c_{i})$,

with the feature Jacobian $\partial_{z}\phi_{2}$ evaluated in closed form through the fixed layers, and stacking the $A_{i}$ and $f_{i}$ reduces the fit to a ridge regression,

with regularization parameter $\alpha$. For fixed features and hyperparameters, the readout fit is a convex least-squares problem, and Eq. (5) is its unique solution, with no iterative optimization involved. Hamiltonian values are not required for training, and the additive constant of $\widehat{H}_{\beta}$ left undetermined by field matching does not affect the vector field. Hamilton’s equations applied to $\widehat{H}_{\beta}$ then define the autonomous flow used for all rollouts and diagnostics.

The scales of the random features act as internal hyperparameters and, as in reservoir computing, are optimized rather than fixed in advance. The state weights $W_{z}$ and control weights $w_{c}$ are independent standard normal draws multiplied by $s_{z}$ and $s_{c}$, respectively. The second-layer matrix $W$ is a standard normal draw multiplied by $g_{2}/\sqrt{N}$, where $N$ is the common width of the two feature layers. The factor $1/\sqrt{N}$ offsets the $N$-term sums in $W\phi_{1}$, so that $g_{2}$ sets the second-layer scale independently of the width. Every bias component is drawn independently and uniformly from $[-1,1]$. For each of the three candidate activations ($\tanh$, $\mathrm{GELU}$, and $\mathrm{softplus}$), 100 Bayesian-optimization trials independently optimize $N$, $s_{z}$, $s_{c}$, $g_{2}$, and the regularization parameter $\alpha$ [1]. Each trial fits the readout at a designated subset of the training controls and is scored by the root-mean-square error (RMSE) of the predicted vector field on an independently sampled shell at one validation control inside the training interval. The activation and all remaining hyperparameters are then selected by that validation RMSE. No extrapolation control is used in model selection. Because each trial costs only one closed-form ridge fit, the 300 trials per family remain affordable. After selection, the readout is refitted at all training controls using 8000 independently sampled shell states per control, and reported results average three random-feature seeds. Full details of the search and the selected configurations are given in Appendix B.

Training states are sampled independently rather than extracted from trajectories. We draw $(x,y)$ uniformly from a family-specific rectangular proposal box and reject any point with potential energy $V(x,y;c)\geq E_{0}$, for which the shell leaves no kinetic energy. At each accepted point, the remaining energy $E_{0}-V$ fixes the momentum magnitude and leaves only its direction free: an angle $\theta$ drawn uniformly on $[0,2\pi)$ gives $(p_{x},p_{y})=[2(E_{0}-V)]^{1/2}(\cos\theta,\sin\theta)$, which places the state on the designated energy shell. Appendix A gives the proposal boxes.

For comparison, we train a conventional HNN in which every weight is trained end to end by back-propagation on all training controls, hereafter the MLP-HNN: two hidden $\tanh$ layers of width 200 map $(z,c)$ to a scalar $H_{\theta}$, following the standard parameter-aware design [12, 13]. The baseline receives the same field-matching objective and states from the same designated energy shell as the RF-HNN, and its architecture and training budget are fixed without consulting any extrapolation control. We evaluate both single networks and a 20-network ensemble, whose prediction is the gradient of the mean of the learned Hamiltonians. This is a controlled comparison of fitting methods on data drawn from the same shells and controls. The training details are given in Appendix A.

We first test the canonical family, whose potential is unmodified, and then follow the transition deeper in the bounded family, whose confining wall permits extrapolation beyond the escape limit of the canonical potential. Section IV.4 then extends the comparison to all four families. All reported occupancy, Chamfer, and Lyapunov statistics are computed from the 28-orbit sample of Sec. III.1. Truth and every model are advanced by the same classical fourth-order Runge–Kutta implementation with the family time steps $\Delta t$ of Table 1, so that model-specific time steppers cannot confound the comparison and the finite-step integration error is common to all flows.

The main result of this study is that a parameter-aware RF-HNN can extrapolate the emergence of a broad chaotic regime from data dominated by regular motion. No data from the extrapolation regime enter training or model selection, yet the learned Hamiltonian reproduces the breakup of invariant structures and the growth of chaotic components at unseen controls across four Hamiltonian families. Because the flow of the learned scalar Hamiltonian is integrated autonomously, the prediction is not a direct classification of chaos or a short-time extension of observed trajectories. The out-of-range continuation of the learned Hamiltonian generates the long-time phase-space organization and instability of the new regime.

Evaluating this prediction requires measures that survive the rapid decorrelation of chaotic trajectories. We therefore use Poincaré-section occupancy, section geometry, and finite-time Lyapunov exponents, rather than long-horizon pointwise error, to probe the extent, location, and instability of the learned flow. Their combined agreement is important because the transition appears differently across the systems, with some showing a strong change in occupied area and others a clearer increase in instability. The consistent RF-HNN behavior across different potentials, symmetries, and couplings indicates that the result is not tied to one Hénon–Heiles parameterization.

The contrast with the MLP-HNN baselines exposes an ambiguity specific to parameter extrapolation. Matching the field $J\nabla_{z}H$ on energy shells at a finite set of controls constrains the vector field where data are available, but it does not uniquely determine its continuation in $c$. Distinct Hamiltonian models can therefore reproduce the observed regular regime while generating different phase-space organizations outside it. This is what the comparison shows. The MLP-HNNs reproduce the regular in-band behavior (inside the training interval), while their continuations delay or suppress the breakup of invariant structures and the associated growth of instability. The ensemble does not repair this discrepancy because averaging reduces member-to-member variance, not a bias shared by the learned continuations. The relevant distinction is therefore not only between Hamiltonian and non-Hamiltonian learning, but among Hamiltonian continuations that are nearly indistinguishable on the observed interval.

The mechanism by which the RF-HNN selects the more faithful continuation is not yet isolated, but the results suggest a specific role for its learning procedure. Once the random features are fixed, fitting is a regularized linear inverse problem, and the feature scales set how rapidly the representation can vary in state and control. In-band validation over these scales may therefore favor a lower-complexity continuation in $c$ without any access to the target regime, and the convex readout makes this selection computationally practical. The exploratory tests in Appendix C disfavor the simplest alternatives: increasing the MLP-HNN width, changing its activation, or refitting its final layer by ridge regression does not systematically recover the transition. These tests do not identify a unique mechanism, however. Matched ablations that equalize in-band approximation error and effective capacity will be needed to separate the roles of the fixed representation, the control-channel scaling, the regularized readout, and the hyperparameter selection.

A second implication is that the predictability of a dynamical regime can survive the unpredictability of its individual trajectories. Beyond the transition, pointwise rollouts decorrelate rapidly, yet the RF-HNN reproduces the occupied region of the section, the surviving regular structures, and the finite-time instability of the extrapolated flow. Vector-field data from a regular-dominated regime can thus contain enough information to extrapolate the parameter dependence that governs later resonance overlap, provided that the learning procedure selects an appropriate continuation. The claim is consequently about long-time phase-space organization, not about pointwise forecasting or global recovery of the Hamiltonian.

The present study establishes these results in a controlled setting: synthetic two-degree-of-freedom families, one smoothly varying control, clean vector-field data on a single energy shell, and moderate extrapolation distances. Extending the analysis to noisy and sparsely sampled trajectories, multiple energies or controls, and higher-dimensional Hamiltonian systems will determine how broadly the result persists. A further challenge is to decide, from in-band information alone, when an extrapolated phase portrait should be trusted, because a small in-band field error does not by itself certify the extrapolated dynamics.

The central conclusion is that matching the Hamiltonian vector field at the observed controls and reproducing the phase-space transition of the underlying family beyond them are distinct requirements. The former can be optimized and validated inside the observed interval, whereas the latter depends on how the fitted Hamiltonian continues in the control beyond it. In the four families studied here, the RF-HNN meets both requirements, generating the breakup of invariant structures and the growth of a chaotic sea absent from training and model selection.

The code and data that support the findings of this study are available from the corresponding author upon reasonable request.