1 Introduction↩︎

The Hamiltonian formalism is standard in mathematical physics for expressing dynamics in many physical systems. Beyond yielding a systematic framework for deriving the dynamics of the system, there are deep structural, geometric, topological, and analytic properties of Hamiltonian functions that can be used to understand the physics of the system and its dynamics [1], [2]. Hamiltonian systems are described by a single scalar function in the phase space \(\mathcal{H}: \Omega \subseteq \mathbb{R}^{2d} \rightarrow \mathbb{R}\) and the corresponding system for \(y \in \Omega\) in Hamilton’s equations [3], [4].

\[\label{eq:Hamilt95ODE} \begin{align} \dot{y} = {J}^{-1}\nabla \mathcal{H}(y), \quad {J} = \begin{pmatrix} 0 & 1\\ -1 & 0 \end{pmatrix}, \end{align}\tag{1}\] where \(\mathcal{J}\) is the canonical symplectic matrix, which we will refer to repeatedly in subsequent sections. The solutions of the system are trajectories \(y(t)\). If the Hamiltonian does not have an explicit time-dependence, then it is the Noether quantity associated with time-translational symmetry and is conserved along the flow, i.e., \(\mathcal{H}(y(t))\) remains constant along trajectories. This implies many important conservation properties across physical systems of interest. For a set of canonical coordinates {\(p_i\), \(q_i\)} \(\in \mathbb{R}^{2d}\), an integrator, or numerical simulation method, is said to be symplectic if it preserves the canonical symplectic form, \(\omega = \sum_{i=1}^{n} dp_i \wedge dq^i\); see subsequent sections and also [5].

Geometric integrators preserve geometric invariants of the flow, and the use of such integrators has been instrumental in understanding the long-term qualitative properties of many important physical systems [6]–[10]. As such, by learning the Hamiltonian itself rather than seeking to learn trajectories, we are able to directly learn the qualitative properties of the system. In contrast, non-geometric numerical methods such as the Euler and explicit Runge–Kutta methods do not conserve the Hamiltonian [7]. A result in [11] shows that it is not possible for a fixed-timestep integrator to simultaneously preserve the energy and the symplectic structure unless it samples the exact solution, and this result inspired the development of energy-preserving and symplectic integration schemes.

In recent years, the field of Physics-Informed Machine Learning has brought significant advances in architecture and learning techniques. By enforcing physical inductive biases while learning the dynamics through appropriate structural modifications to standard NN training, the learning process requires much fewer samples to accurately fit the data and achieves better out-of-distribution accuracy. This includes the recent development of Hamiltonian Neural Networks [12]. By representing the Hamiltonian as a Neural Network, a more physical representation of the system is available. While incorporating symplectic integrators yields more qualitatively accurate approximations of the Hamiltonian dynamics, it poses technical challenges, namely (i) reconstructing a learned Hamiltonian function from noisy trajectory observations, (ii) solving the implicit equations required by symplectic integration methods, and (iii) ensuring that the learned Hamiltonian model generalizes to out-of-distribution data. Motivated by these technical challenges, the remainder of this paper is structured as follows:

In Section 2, we provide background on two complementary formulations of dynamical systems: the Lagrangian and Hamiltonian perspectives. Section 3 reviews symplectic integrators and their application to solving general Hamiltonian systems. We also discuss briefly backward error analysis in the context of symplectic methods. In Section 4, we describe the architecture of our Hamiltonian Neural Network and introduce an implicit symplectic integrator that is implemented with a predictor-corrector method in the forward pass. Section 4.2 describes the backward pass and introduces the adjoint sensitivity equations for Hamiltonian Learning. Finally, Section 5 presents numerical experiments.

1.1 Contributions↩︎

The paper makes the following contributions:

• Showing that a neural framework can be constructed which learns generalized Hamiltonians up to some error tolerance from noisy observations of trajectories without having any implicit structural bias (no assumptions on the separability of the Hamiltonian).

• Showing that the Hamiltonian Neural Net generalizes to out-of-distribution data and is able to learn the governing Hamiltonians for systems with chaotic dynamics.

• Using the fact that the adjoint sensitivity of a symplectic discretization of a Hamiltonian ODE can be computed using the same symplectic discretization applied to the adjoint system, we are able to compute exact gradients of the loss with respect to neural network parameters without the need for backpropagation.

• We also describe in brief how backward error analysis applied to our Learned Hamiltonian results in a higher order reconstruction of the true Hamiltonian.

2 Hamiltonian Dynamics: A Background↩︎

Before discussing the Hamiltonian formulations, we recall the definition of a manifold in the context of mechanics.

There are two equivalent and foundational approaches to formulating classical mechanics, rooted in geometric and variational principles: the Lagrangian and Hamiltonian frameworks. These frameworks describe the evolution of physical systems using generalized coordinates and velocities or momenta, rather than relying directly on Newton’s second law. The Lagrangian formalism is based on the principle of stationary action and encodes the dynamics in terms of the Lagrangian, which is difference between kinetic and potential energy. The Hamiltonian formalism reformulates the problem as a system of first-order differential equations on phase space, where the state of the system is expressed through canonical coordinates and conjugate momenta. These approaches systematically reveal conserved quantities through symmetries, as formalized by Noether’s theorem, and reveals the mathematical structure underlying the dynamics, including the preservation of the symplectic geometric structure in the Hamiltonian setting. Beyond their classical role, these formulations have become central in the analysis of dynamical systems, numerical methods, and emerging applications in optimization and machine learning.

2.1 The Lagrangian formulation↩︎

Lagrangian mechanics describes the motion in a mechanical system by means of the configuration space. The configuration space of a mechanical system has the structure of a differentiable manifold, on which its group of diffeomorphisms acts. The basic ideas and theorems of Lagrangian mechanics are invariant under the right action of this group [13], which reflect the coordinate independence of this formulation. A Lagrangian mechanical system is given by a configuration manifold and a Lagrangian function on its tangent bundle. Every one-parameter group of diffeomorphisms of the configuration space induces a tangent lifted action on the tangent bundle. If the Lagrangian is invariant under this tangent lifted action, then there is an associated Noether quantity that is conserved. A Newtonian potential system is a particular case of a Lagrangian system where the configuration space is Euclidean and the Lagrangian function is the difference between kinetic and potential energies. The equations of motion are derived via the variational principle by extremizing the action functional 2 whose domain is the infinite-dimensional space of curves and the dynamics is defined by the curve that extremizes this action¹.

Let \(Q\) be a \(d\)-dimensional smooth manifold representing the configuration space of a mechanical system. The tangent bundle \(TQ\) of \(Q\) consists of all pairs \((q, \dot{q})\), where \(q \in Q\) denotes a configuration and \(\dot{q} \in T_q Q\) is a tangent vector representing the velocity at \(q\). Given local coordinates \((q_1, \ldots, q_d)\) on \(Q\), the tangent bundle has local coordinates \((q_i, \dot{q}_i)\) for \(i = 1, \ldots, d\). Lagrangian mechanics is formulated on \(TQ\) with the action functional \(S: C^2([t_0,t_1],Q) \to \mathbb{R}\) is defined as: \[\label{eq:action} \begin{align} S[q] = \int_{t_0}^{t_1} L(t, q(t), \dot{q}(t)) dt. \end{align}\tag{2}\] Then, the dynamics of the system is governed by Hamilton’s principle: \[\delta S[q] = 0,\] for all variations \(\delta q(t)\) of \(q(t)\) that vanish at the endpoints, i.e., \(\delta q(t_0) = \delta q(t_1) = 0\). Computing the variation, integrating by parts, and evoking the fundamental theorem of the calculus of variations yields the Euler–Lagrange equations: \[\begin{align} \label{eq:lagrange} \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = 0, \quad \text{for } i = 1, \ldots, d. \end{align}\tag{3}\]

2.2 The Hamiltonian formulation↩︎

Another approach is the Hamiltonian formulation of dynamics, which provides an alternative (but equivalent) description of the system on phase space when the Lagrangian is hyperregular. In this picture, one works on the cotangent bundle \(T^*Q\) (the space of pairs \((q,p)\) of generalized coordinates and conjugate momenta) instead of the tangent bundle \(TQ\) consisting of \((q,\dot{q})\) used in Lagrangian mechanics. The cotangent bundle \(T^*Q\) carries a natural symplectic structure, which is central to Hamiltonian mechanics. Specifically, let \(\mathcal{M} = T^*Q\) be a \(2d\)-dimensional smooth manifold. Then the symplectic form \(\omega\) is a closed, non-degenerate 2-form on \(\mathcal{M}\), and in canonical (Darboux) coordinates \((q_1, \ldots, q_d, p_1, \ldots, p_d)\) on \(T^*Q\), it is given by: \[\begin{align} \label{eq:symplecticform} \omega = \sum_{i=1}^d dq_i \wedge dp_i. \end{align}\tag{4}\]

To transition from the Lagrangian \(L: TQ \to \mathbb{R}\) to the Hamiltonian formalism, we first define the conjugate momenta by taking partial derivatives of \(L\) with respect to the generalized velocities. In local coordinates \(q_i\)² on \(Q\) , the \(i\)-th component of momentum is defined as:

\[\begin{align} p_i = \frac{\partial L}{\partial \dot{q_i}}, \quad \text{for } i=1,\dots,d \end{align}\] This definition induces a map known as the Legendre transform, often denoted \(\mathbb{F}L: TQ \to T^*Q\). The Legendre transform sends a point \((q,\dot{q})\) in the tangent bundle to a corresponding point \((q,p)\) in the cotangent bundle by pairing velocities with momenta. In coordinates, \(\mathbb{F}L(q,\dot{q}) = (q,p)\), where \(p_i = \partial L/\partial \dot{q}_i\)³.

Given this correspondence, we define the Hamiltonian \(H: T^*Q \to \mathbb{R}\) as the Legendre transform of the Lagrangian, i.e., the function whose value is the “energy” obtained by trading the velocity dependence of \(L\) for momentum dependence. It is given by the formula, \[H(q, p) = \sum_{i=1}^d p_i \dot{q}_i - L(q, \dot{q}),\label{H95in95terms95of95L}\tag{5}\] where \(\dot{q}_i\) is expressed as a function of \((q, p)\) by inverting the Legendre transform. Equivalently, the dependence on the velocities on the right-hand side can be eliminated by extremizing with respect to the velocities, which is analogous to the approach adopted in the Pontryagin maximum principle. Given \(H\), we can define a unique vector field \(X_H\) on \(T^*Q\), the Hamiltonian vector field, by the condition \[\begin{align} \label{Hamilton95eq95abstract} dH = \omega(X_H,\cdot). \end{align}\tag{6}\] The Hamiltonian vector field \(X_H\) takes the form \[\begin{align} X_H = \sum_{i=1}^{d} \frac{\partial H}{\partial p_i}\frac{\partial }{\partial q_i} - \frac{\partial H}{\partial q_i} \frac{\partial}{\partial p_i}. \end{align}\] The integral curves of \(X_H\) are the solutions to Hamilton’s equations \[\dot{p}_i = -\frac{\partial H}{\partial q_i}, \quad \dot{q}_i = \frac{\partial H}{\partial p_i}, \quad \text{for } i = 1, \ldots, d.\] While we have derived Hamilton’s equations using local canonical coordinates, 6 defines a global vector field on \(\mathcal{M}\). These equations describe the flow of the system in phase space and are equivalent to the Euler–Lagrange equations 3 derived from Hamilton’s principle if the Legendre transformation is globally invertible and the Lagrangian \(L\) is related to the Hamiltonian \(H\) by 5 . For completeness, the Hamilton–Jacobi partial differential equation, \[\partial_t S + H(q, \partial_q S) = 0,\] describes the generating function \(S\) for the canonical transformation that maps \((q(t_0),p(t_0))\) to \((q(t_1),p(t_1))\). This provides an alternative method for solving Hamilton’s equations [14], and Jacobi’s solution to the Hamilton–Jacobi equation is given in terms of the action functional evaluated along the two-point boundary value solution of the Euler–Lagrange equations.

3 Symplectic Integrators↩︎

This section presents symplectic integration in the context of the current Hamiltonian learning literature. Then, we will review some properties of symplectic flow maps which motivate their incorporation into the learning process for dynamical systems. Finally, we present some standard methods for constructing symplectic integrators that we later use in this work.

3.1 Symplectic Integration in NNs: Previous Work↩︎

Since the original proposal by [15] and the concurrent work by [12], HNNs have generated much scientific interest. This has spawned generative [16], recurrent [17], and constrained [18] versions, as well as Lagrangian Neural Networks [19] have been proposed. In contrast to standard HNNs, the approach introduced in [17] directly optimizes the actual states observed at each timestep for a given initialization by integrating the Hamiltonian vector field using a symplectic integrator (leapfrog algorithm) and backpropagating each squared error through time. The state at the next timestep is predicted using the symplectic integrator; in this way the entire time series is predicted, which is then compared with the observed data. Note that they make the assumption that the Hamiltonian is separable, which is significant as the leapfrog algorithm is generally implicit, but if the Hamiltonian is separable, then the algorithm becomes explicit. Notably, the work in [20] uses an implicit midpoint scheme to train a Poisson Neural Net, which learns the Poisson structure of a dynamical system, and the network is trained using a standard backpropagation technique.

Most of these works considered either explicit integration schemes, some semi-implicit schemes, or used a separable ansatz for the neural nets. In [21], [22], additive separability biases were introduced in the HNN architecture/training, allowing the network to learn Hamiltonians of the form \(H(q,p)=T(p)+V(q)+\text{const}\). This yields better performance by making the Störmer–Verlet (leapfrog) integrators fully explicit and easier to train, at least when the true Hamiltonian is well-approximated by a separable Hamiltonian. The error bounds for such symplectic HNNs have been analyzed for the noiseless case in recent works (e.g., showing that energy error grows linearly in time under a symplectic integrator) [23], [24], and it was established that the learned \(H\) is reasonable but is slightly perturbed from the true \(H\). These limitations highlight that while semi-implicit methods make training feasible, they come at the cost of a small modeling bias as the integrator inherently assumes that the Hamiltonian is separable in order to simplify the computation. More details can be found in [25], [26].

A generalized HNN framework is proposed in [27] proposes, which can be used for non-separable Hamiltonians where they approximate the original non-separable Hamiltonian by an augmented one that was introduced in [28], with an extended phase space and a tunable parameter \(\omega\) which controls the binding between the two copies of the Hamiltonian and model them using NNs. The approach in [28] yields an explicit symplectic integration scheme for non-separable Hamiltonians, at the expense of increased complexity associated with two coupled Hamiltonian systems, and the need to tune the coupling parameter. As we mentioned earlier, all such Symplectic HNN architectures learn a valid Hamiltonian which is a perturbation of the true Hamiltonian, due to the approximation error in replacing the exact flow map with a numerical integrator. Establishing the error bounds is necessary, and there have been a few classic works which derive these error bounds for symplectic integrators, notably [24], [25].

Most existing works have either relied on explicit methods or employed imperfect implicit schemes, often limited to noiseless settings. This is largely due to the computational cost of the nonlinear solves, which can be addressed by using a predictor-corrector approach. The predictor is an explicit one-step method with the same order as the implicit symplectic integrator, and the corrector involves a few fixed-point iterations of the symplectic integrator. Another challenge is posed by backpropagation through implicit ODE solvers, which scales in memory with the simulation length and number of parameters. In [29], a backpropagation-free framework using sampled neural networks is proposed for learning Hamiltonians, which is effective with rich trajectory data but assumes access to the true time derivatives, which is unrealistic in practice. Notably another approach that avoids solving a forward integral is in form of SympNets[30] which learn the Neural NEt itself as a symplectic map but this comes with restrictions on freedom of choosing an arbitrary stepsize during simulation. In contrast, the adjoint method offers greater generality and handles sparse supervision.

If the Hamiltonian is separable, it has the form, \[\label{eq:sep95H} \begin{align} \mathcal{H}(q,p) = T(p) + V(q) \end{align}\tag{7}\] and the Hamiltonian vector field can be written as \[\label{eq:seperable95Hamiltonian95dynamics} \begin{align} \dot{q} &= \frac{\partial T(p)}{\partial p}\\ \dot{p} &= -\frac{\partial V(q)}{\partial q} \end{align}\tag{8}\] For separable Hamiltonians, one can then use the standard Störmer-Verlet method, whose numerical update is given by \[\label{eq:sv95int} \begin{align} p_{n+\frac{1}{2}} &= p_n - \frac{h}{2}\nabla V(q_n), \\ q_{n+1}&=q_n + h\nabla T(q_{n+\frac{1}{2}}),\\ p_{n+1}&=p_{n+\frac{1}{2}} - \frac{h}{2}V(q_{n+1}). \end{align}\tag{9}\]

The method is second-order accurate but is only defined when the Hamiltonian is separable and the Hamiltonian vector field is given by 8 . In the next part, we will describe a well-known second-order symplectic method that works for general (non-separable) Hamiltonians.

3.2 Symplectic Integrators for Nonseparable Hamiltonians↩︎

When simulating physical systems, it is desirable to consider geometric integrators that preserve the geometric properties of the flow. More specifically, for Hamiltonian systems, this property is the preservation of the symplectic \(2\)-form described in 4 . An automatic consequence of preserving the symplectic form is that the numerical integrator preserves the phase space volume on the cotangent bundle \(T^*Q\).

There are various kinds of symplectic methods, the simplest of which are the first-order symplectic Euler schemes, which we will not discuss here. Consider the canonical Hamiltonian dynamics where we do not assume separability, \[\label{eq:generic95Hamiltonian95dyn} \begin{align} \dot{q} &= \frac{\partial \mathcal{H}}{\partial p},\\ \dot{p} &= -\frac{\partial \mathcal{H}}{\partial q}. \end{align}\tag{10}\] Observe that we cannot use methods of the form 9 . For systems with stiff dynamics or strongly nonlinear Hamiltonians, the semi-implicit methods may suffer from poor accuracy or numerical instability. In such cases, a fully implicit method is often desirable. The implicit midpoint method is often used to integrate such stiff dynamics, and it is symplectic and thus automatically phase space volume-preserving. Additionally, we observe that the map is symmetric, i.e., \(\Phi_h^{-1} = \Phi_{-h}\). For the Hamiltonian vector field, the implicit midpoint method is written as \[\label{eq:im95itegrator} \begin{align} q_{n+1} &= q_n + h \nabla_p \mathcal{H}\Big(\frac{q_n + q_{n+1}}{2}, \frac{p_n + p_{n+1}}{2}\Big), \\ p_{n+1} &= p_n - h \nabla_q \mathcal{H}\Big(\frac{q_n + q_{n+1}}{2},\frac{p_n + p_{n+1}}{2} \Big), \end{align}\tag{11}\] This update scheme is second-order accuracy and it is fully implicit in both variables \(q\) and \(p\). Hence, it requires a nonlinear solver.

3.3 The Numerical Advantage of Using Symplectic methods↩︎

Consider an initial value problem, \[\begin{align} y' = f(y), \quad y(0) = y_0, \end{align}\] and a numerical method \(\Phi_h\) which approximates the trajectory with a sequence of points \(y_0, y_1,\dots, y_n\). Then, the forward error is concerned with the extent to which the numerical method approximates the exact flow, \[\begin{align} \epsilon = |y(nh) - \Phi_{nh}(y_0)|. \end{align}\] In contrast, backward error analysis is concerned with constructing a modified differential equation \(\tilde{y}' = f_h(\tilde{y})\) whose exact flow agrees with the numerical method. This relationship is expressed in Figure 1, where \(\varphi_t\) is the exact flow map for the ODE. It is in this context that the true power of symplectic integrators for simulating Hamiltonian dynamics is revealed. In general, a one-step method cannot be exactly expressed as the time-\(h\) flow of a differential equation, but we can construct an asymptotic expansion for this modified equation, \[\label{eq:modified95eq} \begin{align} \dot{\tilde{y}} = f_h(\tilde{y}) = f(\tilde{y}) + hf_1(\tilde{y}) + h^2f_2(\tilde{y})\dots ,\qquad \tilde{y}(0) = y_0 \end{align}\tag{12}\] Where \(f_i\) denotes terms involving up to \(i\)-th derivatives of \(f\). If the numerical method is \(p\)-th order accurate, then the expansion \(f_h\) agrees with \(f\) up to the \(h^{p-1}\) term, or that the first nontrivial correction term is \(h^p\).

This backward error analysis can be applied to any numerical method, but when applied to a symplectic integrator, the modified vector field is Hamiltonian, which allows us to view a symplectic integrator as being very well approximated by the exact flow of a modified Hamiltonian. The precise theorem is as follows:

The theorem above implies that the numerical flow generated by the numerical map \(\Phi\) is well approximated by the exact flow the modified Hamiltonian: \[\label{eq:modified95Hamiltonian} \begin{align} \tilde{\mathcal{H}} = \mathcal{H}+ h^{p}\mathcal{H}_{p}+\dots \end{align}\tag{13}\] where \(p\) denotes the order of accuracy of the numerical integrator and \(\mathcal{H}_p\) satisfies \(f_j=J^{-1}\mathcal{H}_j\) in the modified vector field. This is significant as symplectic integrators can be viewed as the flow of a modified Hamiltonian whereas the modified vector field associated with a non-symplectic integrator cannot be associated with a modified Hamiltonian flow. This ensures that the numerical reconstruction more physical by preserving the total energy. Moreover, we will show in Section 5.4.1 an additional benefit this presents for learning physical systems. In essence, it allows us to use a lower-order symplectic integrator to learn a lower-order approximation of the true Hamiltonian, and use backward error analysis to post-process the learned Hamiltonian to obtain a higher-order approximation of the true Hamiltonian. For the implicit midpoint method 11 that we use in our analysis, the modified Hamiltonian is given by \[\label{eq:modified95Hamiltonian95IM} \begin{align} \tilde{\mathcal{H}} = \mathcal{H}- \frac{h^2}{24} \nabla^2\mathcal{H}(J^{-1}\nabla \mathcal{H}, J^{-1}\nabla \mathcal{H}) + \mathcal{O}(h^4). \end{align}\tag{14}\] Note that there terms involving odd powers of \(h\) vanish, which is a consequence of the fact that the method is symmetric. The expression \(\nabla^2\mathcal{H}(J\nabla \mathcal{H}, J\nabla \mathcal{H})\) represents a bilinear form that is explicitly given by \[\begin{align} \nabla^2\mathcal{H}(J^{-1}\nabla \mathcal{H}, J^{-1}\nabla \mathcal{H}) = (J^{-1}\nabla \mathcal{H})^\top\nabla^2\mathcal{H}(J^{-1}\nabla \mathcal{H}). \end{align}\] More about the modified Hamiltonian will be presented in the numerical section.

3.4 Constructing Symplectic integrators↩︎

3.4.0.1 Partitioned Runge–Kutta Methods

Partitioned Runge–Kutta (PRK) methods are a class of numerical integrators that are particularly well-suited for systems where the state can be naturally split into multiple components, such as position and momentum in Hamiltonian systems as given by: \[\label{eq:prk95dyn} \dfrac{dq}{dt} = \frac{\partial\mathcal{H}(q,p,t)}{\partial p}, \quad \dfrac{dp}{dt} = -\frac{\partial\mathcal{H}(q,p,t)}{\partial q}.\tag{15}\] Now generally, when we are solving for systems, we can combine these equations in a single vector of state \(y(t)\) and iteratively solve for \(y_n\) using RK methods, however, for Hamiltonian systems these variables play different roles and may evolve at different rates hence these methods treat these two set of variables differently. Now, equation 15 can be integrated using a partitioned Runge–Kutta scheme: \[\begin{align} q_{n+1} &= q_n + h_n \sum_{i=1}^{s} b_i k_{n,i},\\ \qquad p_{n+1} &= p_n + h_n \sum_{i=1}^s B_i l_{n,i}, \end{align}\] where \[\begin{align} k_{n,i} = f(Q_{n,i}, P_{n, i}, t_n + c_ih_n), \\ \\ \qquad l_{n,i} = g(Q_{n,i}, P_{n,i}, t_n + C_ih_n). \end{align}\] which are evaluated at the internal stages, \[\begin{align} Q_{n,i} = q_n + h_n \sum_{j=1}^{s} a_{ij}k_{n,j}, \qquad P_{n,i} = p_n + h_n \sum_{j=1}^s A_{ij} l_{nj}. \end{align}\] A partitioned Runge–Kutta scheme is symplectic if the following conditions hold: \[\label{eq:symp95condition} \begin{align} c_i = C_i, b_i = B_i, \qquad i &= 1, ..., s; \\ b_i A_{ij} + B_j a_{ji} - b_i B_j = 0, \quad i,j&=1,...,s. \\ \end{align}\tag{16}\] All the symplectic integrators discussed above, including the one used in our subsequent analysis, are examples of symplectic PRK schemes.

3.4.0.2 Hamiltonian Variational Integrators

Discrete variational analysis is another way to construct symplectic integrators. These apply a discrete variational principle, which automatically guarantees symplecticity. In fact, discrete Lagrangian mechanics is expressed in terms of a discrete Lagrangian \(L(q_k, q_{k+1})\) that is a Type I generating function and can be used to generate symplectic maps. However, if the Hamiltonian is not hyperregular and the Legendre transformation is not invertible, as is the case for \(\mathcal{H}= qp\), it is not possible to transform the Hamiltonian into an equivalent Lagrangian and apply the discrete Lagrangian construction. Instead, it was showed in [32] that a discrete Hamiltonian \(H^+_d(q_k, p_{k+1})\) that is a Type II generating function can be constructed directly from the continuous Hamiltonian, thereby avoiding the need to evoke the Lagrangian perspective. As with the Lagrangian theory, the discrete Hamiltonian perspective extremizes a discrete action sum, yielding a Hamiltonian variational integrator that is automatically symplectic.

4 Fully Symplectic Hamiltonian Neural Network↩︎

Now, we are at the point where we can apply the insights from symplectic integration theory to the data-driven setting. Consider the problem of learning a Hamiltonian \(\mathcal{H}\) from data, with noisy observations \((\tilde{\mathbf{q}}, \tilde{\mathbf{p}}) \in \mathbb{R}^{2d}\). At this point let us set up a few assumptions

• We have a set of noisy trajectories \(N_{train}\) in the training set and \(N_{val}\) in the validation set.

• Our integrator is represented by a flow map \(\Phi\) such that \(y_T = \Phi(y, f(y), 0, T, h)\), where \(y\) is the initial state, \(f(y)\), the generator of dynamics, \(T\) is the final time and \(h\) is stepsize of the integrator.

• The Hamiltonian Neural Network is represented by \(\mathcal{H}_{\theta}(\mathbf{q}, \mathbf{p})\) where \(\theta\in \mathbb{R}^p\), \((\mathbf{q}, \mathbf{p})\in \mathbb{R}^{2d}\) is a combined vector of phase space coordinates and are also the inputs to our neural network.

Now, the problem boils down to learning the functional form of Hamiltonian \(\mathcal{H}_{\theta^*}(\mathbf{q}, \mathbf{p})\) where inference in the network, for given trained \(\boldsymbol{\theta}^*\), amounts to computing the value of Hamiltonian function \(\mathcal{H}_{\theta^*}(\mathbf{q}, \mathbf{p})\) over phase space at a point \((\mathbf{q}, \mathbf{p}) \in \mathbb{R}^{2d}\). We start at a set of initial phase space coordinates \((\mathbf{q}, \mathbf{p})\) and simulate the trajectory forward using a symplectic integrator to compute up to final time \(T\) to get \(\{\mathbf{q}_T, \mathbf{p}_T\}\). At this stage, we construct a loss function \[\label{eq:loss95function} \begin{align} Loss = \frac{1}{n}\Big(\sum_{i=1}^N\|\mathbf{q}^i_T - \tilde{\mathbf{q}}^i_T\|^2 + \|\mathbf{p}^i_T - \tilde{\mathbf{p}}^i_T\|^2\Big). \end{align}\tag{17}\] Here, \(i\) represents a data sample/trajectory in phase space and \(n\) is the batch size since we compute a gradient estimate over a minbatch. The vector \((\tilde{\mathbf{q}},\tilde{\mathbf{p}})\) represents the noisy ground truth values.

4.1 Numerical Forward Pass with Implicit-Midpoint method↩︎

In the forward pass, as represented by lower block in schematic 2, we solve the Hamiltonian ODEs in equation 18 to generate the forward trajectory up to final time \(T\).

\[\label{eq:hamilton95odes} \begin{align} \dfrac{dq_i}{dt} = \frac{\partial \mathcal{H}}{\partial p_i}, \quad \dfrac{dp_i}{dt} = -\frac{\partial \mathcal{H}}{\partial q_i}. \end{align}\tag{18}\]

As a choice of integrator, we choose the second-order implicit midpoint method with updates given by equation 19 , \[\label{eq:im95int} \begin{align} q_{i+1} &= q_i + h\nabla_p\mathcal{H}_{\theta}\left(\frac{q_{i} + q_{i+1}}{2},\frac{p_{i} + p_{i+1}}{2}\right),\\ p_{i+1} &= p_i - h\nabla_q\mathcal{H}_{\theta}\left(\frac{q_{i} + q_{i+1}}{2},\frac{p_{i} + p_{i+1}}{2}\right). \end{align}\tag{19}\] The above equations are implicit in the unknowns and thus requires a fixed-point iteration scheme to solve. The details of the integrator with fixed-point iteration is discussed in Algorithm 4.

Starting from an initial point \(\mathbf{y}\equiv (\mathbf{q}, \mathbf{p})\), a preliminary estimate for the next step \(\mathbf{y}'\) is first obtained using an explicit predictor, followed by refinement through repeated application of the implicit midpoint update via fixed-point iteration. Alternatively, when access to the full trajectory \((\mathbf{q}_i, \mathbf{p}_i)\) is available, as is the case in training against noisy trajectories, one may optionally bypass the predictor step and directly use the noisy values as an initial guess.

In our implementation, the simulation starts from an initial point \(\mathbf{y}\equiv (\mathbf{q}, \mathbf{p})\). Algorithm 4 is first applied with an explicit update to obtain a predicted value \(\mathbf{y}'\), which is then used as the initial guess for the fixed-point iteration. For simulating up to final time \(T\), each solve runs for \(\tau\) fixed-point iterations.

4.2 Backward pass with gradient computation via the Adjoint Sensitivity method↩︎

To perform an optimization over parameters \(\theta\) to obtain the final \(\theta^*\), we need to compute the sensitivity of the loss function w.r.t the parameters. As discussed in the previous section, pairing the ODE solver with fixed-point iterations incurs multiple function evaluations per timestep. To compute gradients for optimization, the automatic differentiation engine must retain the entire sequence of intermediate operations, to backpropagate gradients through the solver and evaluate the sensitivity of the loss with respect to \(\boldsymbol{\theta}\). This standard approach, often implemented via reverse-mode automatic differentiation (backpropagation), is memory-intensive, as it requires storing all intermediate values of \((\mathbf{q}_t, \mathbf{p}_t)\) and internal solver states across the entire integration window. This memory overhead grows linearly with the number of timesteps and becomes prohibitive for long-horizon simulations or high-dimensional dynamical systems.

To overcome this bottleneck, one can instead formulate the gradient computation as a boundary-value problem via the adjoint sensitivity method. This method derives a backward-in-time differential equation for the adjoint variables (or co-states), which represent the sensitivity of the loss functional with respect to the trajectory. Crucially, the backward pass is decoupled from the original forward pass: the adjoint equations can be solved without retaining the entire forward trajectory. This reformulation allows memory-efficient gradient computation, particularly when the number of scalar outputs (loss terms) is much smaller than the number of parameters.

In this section, we merely give the form of adjoint equations that we solve to get the gradient of the loss function w.r.t network parameters. In the next section, we give a detailed discussion of this approach. Given a Neural Network parameter estimation model subject to Hamiltonian ODEs

|l|
 {\theta}{\Jc(\boldsymbol{\theta}, q, p) = \frac{1}{n}\sum \limits_{i=0}^n \|q^i_T - \tilde{q}^i_T\|^2_2 + \|p^i_T - \tilde{p}^i_T \|^2_2}{}{}
\addConstraint{\dot{q}(t) - \frac{\partial\mathcal{H}_{\theta}(q, p)}{\partial p}}{=0, \hspace{1cm} q(0) = q_0}
\addConstraint{\dot{p}(t) + \frac{\partial\mathcal{H}_{\theta}(q, p)}{\partial q}}{=0, \hspace{1cm} p(0) = p_0},
\label{eq:opt95prob951}

where \((q_T, p_T)\) is the predicted value of states through the symplectic map \(\Phi_h\) at final tim \(T\), \(h\) is the stepsize, and \((\tilde{q}_T, \tilde{p}_T)\) is the true value of those variables. For simplicity, here we consider \(q\) and \(p\) as scalar functions. Now, considering \(\Phi_h\) as our symplectic map 19 and using a shorthand notation \(y_t\) for \((q_t, p_t)\), we can write \[\begin{align} \mathbf{y}_t = \Phi_h(\mathbf{y}_t, \mathbf{y}_{t-1}). \end{align}\] In the usual cases where a neural network is used as a function approximator, the gradients are computed using backpropagation which is a standard way to optimize the loss function by adjusting the network weights. In our case, the forward pass includes a symplectic ODE solver which means we have to backpropagate through this solver which will drastically increase the memory requirements, which can be understood if we consider our loss function \[\begin{align} \mathcal{J} = \frac{1}{n}\sum\limits_{i=1}^n (y_T - \tilde{y}_T)^2, \end{align}\] where each \(y_T\) is obtained via forward propagation through a symplectic ODE solver. It is clear that if we are solving for larger simulation times, the number of evaluations in the backward pass will scale with the number of timesteps as shown below: \[\begin{align} \frac{d\mathcal{J}}{d\boldsymbol{\theta}} = \sum\limits_{j=1}^t\Big( \prod\limits_{k=j+1}^t\frac{\partial y_k}{\partial y_{k-1}}\Big)\frac{\partial y_j}{\partial \theta}, \end{align}\] where each \(y_t\) is a function of \(y_{t}, y_{t-1}, y_{t-2},\dots,y_{1}\). In contrast to this, if we use adjoint sensitivity analysis, we first have to derive the adjoint equations which is a system of ODEs in the adjoint variable \(\lambda(t)\) which is a time dependent version of Lagrange multiplier which get introduced while solving equation [eq:opt95prob951] using variational calculus. Given Hamiltonian ODEs in matrix notation \[\label{eq:Hamiltonian95system} \begin{align} \frac{d}{dt} \begin{pmatrix} q(t) \\ p(t) \end{pmatrix} = \begin{pmatrix} \nabla_p \mathcal{H}(q, p, \boldsymbol{\theta}) \\ -\nabla_q \mathcal{H}(q, p, \boldsymbol{\theta}) \end{pmatrix}, \end{align}\tag{20}\] the system of adjoint ODEs is then given as: \[\label{eq:adj95eqns} \begin{align} \frac{d}{dt} \begin{pmatrix} \lambda^q \\ \lambda^p \\ \end{pmatrix} &= -\begin{pmatrix} \nabla_{qp} \mathcal{H}& -\nabla_{qq} \mathcal{H}\\ \nabla_{pp} \mathcal{H}& -\nabla_{pq} \mathcal{H} \end{pmatrix} \begin{pmatrix} \lambda^q \\ \lambda^p \\ \end{pmatrix}. \end{align}\tag{21}\]

Solving these equations backward in time subject to terminal conditions: \[\label{eq:adj95terminal95condns} \begin{align} \lambda^q(T) = -\frac{\partial \mathcal{J}}{\partial q}_{(t=T)}\text{ and } \quad \lambda^p(T) = -\frac{\partial \mathcal{J}}{\partial p}_{(t=T)}, \end{align}\tag{22}\] gives the adjoint state which is in turn used to calculate the loss gradients by computing the following integral: \[\label{eq:lossgrad} \frac{d\mathcal{J}}{d\theta} = \frac{\partial \mathcal{J}}{\partial \theta}-\int_0^T \lambda^T \frac{\partial f_{\rm aug}}{\partial \boldsymbol{\theta}} dt.\tag{23}\] Introducing the adjoint sensitivity method for gradient calculation will amount to solving the adjoint equations backward in time which is a constant memory task where we only need to store the current variable and its partial derivatives in the memory at any particular instant. The complete derivation of HNN adjoint state and gradients is provided in the subsequent section.

Notice that to get the adjoint state, we have to solve equation 21 subject to the terminal conditions 22 \(t=T\) to \(t=0\) and store the results \(\lambda(t)\) in memory to later compute the integral. Now compare it to backpropagation, where the memory requirements scale with simulation length \(T\), where we have to store all intermediate sensitivities \(\frac{\partial \mathbf{y}_i}{\partial\mathbf{y}_j}\), \(0\leq j\leq i\), but here, we only need to solve the two equations for each pair of states, hence a constant memory operation⁴.

4.3 Derivation of Adjoint sensitivity for Hamiltonian system from variational perspective↩︎

Here, we will present a proof of adjoint sensitivity rooted in variational analysis. For that consider a continuous version of the problem that we are solving:

|l|
 {y}{\Jc(\theta, y) = \int \limits_0^T f(\theta, y, t)dt}{}{}
\addConstraint{\dot{y} - f(y, \theta)}{=0}
\addConstraint{y(0)}{=y_0}.
\label{eq:opt95prob95adj}

where \(\theta\in \mathbb{R}^p\) is a vector of unknown parameters, \(y\in \mathbb{R}^{2d}\) is a vector of state variables and \(h:\mathbb{R}^{2d}\rightarrow\mathbb{R}^{2d}\) is the forward dynamics.

In our case, we use a forward map \(\Phi_h\) to generate the flow \(y_1, y_2\dots y_T\) and evaluate the loss/mismatch at the terminal point \(y_T\).

The loss over a batch of trajectories can be written as: \[\begin{align} Loss = \frac{1}{N} \sum_{i=1}^N \|y_i(T) - \tilde{y}_i(T)\|^2. \end{align}\] In the continuous limit, this can be written as \[\begin{align} \mathcal{J}(y) = (y(T) - \tilde{y}(T))^2. \end{align}\]

Writing the Lagrangian for this, we get

\[\label{eq:lagrangian} \begin{align} \mathcal{L} = \mathcal{J}(y(T)) + \int_0^T \lambda^T(\dot{y} - f(y,\theta)) dt + \mu^T(y(0) - y_0). \end{align}\tag{24}\] Taking variation w.r.t \(y\) of the Lagrangian and following a similar procedure as above after integration by parts gives us

\[\begin{align} \delta \mathcal{L} &= \frac{\partial \mathcal{J}(y(T))}{\partial y}\delta y - \int_0^T\Big( \dot{\lambda}^T y + \frac{\partial f}{\partial y}\lambda \Big)\delta y dt + \Big[\lambda(t)\frac{\partial y}{\partial y}\delta y\Big]_0^T. \end{align}\] Solving further, and collecting terms, we get \[\begin{align} \delta \mathcal{L} &= \Big(\frac{\partial \mathcal{J}(y(T))}{\partial y} + \lambda(T)^T \Big)\delta y(T) - \int_0^T\Big( \dot{\lambda}^T y + \frac{\partial f}{\partial y}\lambda \Big)\delta y dt -\lambda(0)\delta y(0). \end{align}\] Now, by the fundamental theorem of the calculus of variations, the coefficients of \(\delta y\) must vanish, which yields the adjoint equation \[\label{eq:adjoint95dynamics952} \begin{align} \dot{\lambda}(t) = -\frac{\partial f}{\partial y}\lambda(t), \qquad \text{s.t.} \quad \lambda(T) =- \frac{\partial \mathcal{J}(T)}{\partial y(T)}. \end{align}\tag{25}\]

4.3.0.1 Gradient w.r.t parameters:

During the optimization, we require that the forward ODE and the boundary value constraint is always satisfied, formally: \[\begin{align} \dot{y} - f(y,\theta) = 0, \quad y(0) - y_0 = 0. \end{align}\] Thus, the gradient of the loss w.r.t. the parameters is just the gradient of Lagrangian as the other two terms are identically zero, \[\begin{align} \frac{d \mathcal{J}}{d\theta} = \frac{d\mathcal{L}}{d\theta}. \end{align}\] Differentiating the Lagrangian again in 24 w.r.t. \(\theta\) (network parameters): \[\begin{align} \frac{d\mathcal{L}}{d\theta} &= \frac{\partial \mathcal{J}}{\partial \theta} + \frac{\partial \mathcal{J}}{\partial y} \frac{\partial y}{\partial \theta} +\int_0^T \left[ \left( \frac{\partial \dot{y}}{\partial \theta} - \frac{\partial f}{\partial y} \frac{\partial y}{\partial \theta} - \frac{\partial f}{\partial \theta} \right)\lambda \right] dt. \end{align}\] We require that \(\mathcal{J}\) only depends on \(y(T)\), hence the second term vanishes. For the term inside the integral, using integration by parts on \(\int_0^T \frac{\partial \dot{y}}{\partial \theta}\lambda\), we obtain \[\begin{align} \frac{d\mathcal{L}}{d\theta} = \frac{\partial \mathcal{J}}{\partial \theta} + \int_0^T \Bigg[ - \lambda^T \frac{\partial f}{\partial \theta} + \left( - \dot{\lambda} - \frac{\partial f}{\partial y}\lambda\right)\frac{\partial y}{\partial \theta}\Bigg] dt + \Big[\lambda^T\frac{\partial y}{\partial \theta} \Big]^T_0. \end{align}\] Now, the coefficient to the \(\frac{\partial y}{\partial \theta}\) term is exactly our adjoint equation, so when restricted to solutions \(\lambda(t)\) of the adjoint equation, it vanishes identically. Also, we can choose \(\lambda\) such that \(\lambda(T)=0\), both of these choices help us in the sense that now we do not have to compute and store the derivatives of each intermediate state \(y(t)\) w.r.t. \(\theta\). This is a huge gain as tracking \(\frac{\partial y(t)}{\partial \theta}\) for each time is a memory-intensive task.

So, by solving \(\lambda(t)\) backward in time from the initial data \(\lambda(T)=0\) along the adjoint equation, we obtain the following expression, \[\begin{align} \frac{d\mathcal{L}}{d \theta} = \frac{\partial \mathcal{J}}{\partial \theta} - \int_0^T\lambda^T\frac{\partial f}{\partial \theta}dt, \end{align}\] which is also the gradient of the loss function w.r.t. parameters \(\theta\) as reasoned above. Therefore, \[\begin{align} \frac{d\mathcal{J}}{d\theta} = \frac{\partial \mathcal{J}}{\partial \theta}- \int_0^T \lambda^T \frac{\partial f}{\partial \theta} dt. \end{align} \label{grad95case2}\tag{26}\]

4.3.0.2 Symplecticity Analysis

If we look at the form of the adjoint dynamics for this case \[\label{eq:definitions} \begin{align} y \equiv \begin{pmatrix} q \\ p \end{pmatrix}, \qquad f\equiv \begin{pmatrix} \nabla_p\mathcal{H}\\ -\nabla_q \mathcal{H} \end{pmatrix}, \qquad \frac{\partial f}{\partial y} \equiv \begin{pmatrix} \nabla_{qp}\mathcal{H}& \nabla_{pp}\mathcal{H}\\ -\nabla_{qq}\mathcal{H}&-\nabla_{pq}\mathcal{H} \end{pmatrix}, \end{align}\tag{27}\] using these definitions, we see that the adjoint equation can be written as : \[\begin{align} \dot{\lambda} = - J^T \nabla^2 H \lambda \qquad \text{where }\quad J = \begin{pmatrix} 0 & 1\\ -1 & 0 \end{pmatrix}. \end{align}\] We see that the adjoint system has a symplectic structure.

It should be noted that more generally, the adjoint system of any ODE is symplectic. Given an ODE \(\dot{(}y)=f(y)\), the adjoint system is given by the Hamiltonian vector field associated with the formal Hamiltonian \(H(y,p)=p^T f(y)\). In particular, the formal Hamiltonian can be derived from the Pontryagin maximum principle by considering the case of trivial controls.

5 Numerical Results↩︎

In this section, we present the performance of our method on three representative Hamiltonian systems: the double well potential, the coupled Harmonic Oscillator, the Hénon–Heiles system, and the 2-body Kepler problem. These systems are chosen to cover a range of problems, including separable and generalized non-separable Hamiltonians. For training, the data is generated by sampling initial conditions from a uniform random distribution within a specified bounded domain. In contrast, for evaluation, the model is tested on test data sampled from 3 separate distributions different from training data in order to evaluate the ability of the learned dynamics to generalize. All the computations were performed on the RCI SLURM cluster (2 × Intel Xeon Gold 6146, 46 cores @ 3.2 GHz, 384 GB RAM, 2 TB NVMe).

5.1 Implementation details↩︎

Below, we provide the main implementation details, including data generation and hyperparameters for our model.

5.1.1 Data generation↩︎

For the first three experiments, we generated 16,384 initial conditions for the training set and 8,192 for the validation set. We used a high-order symplectic integrator to ensure that trajectories evolve on a symplectic manifold. All trajectories were simulated with a time discretization of \(\Delta t = 0.01\) using a 4th order integrator. Additive Gaussian noise \(\mathcal{N}(0, 1)\) was introduced at each timestep, scaled by a noise level coefficient of \(0.01\). The details on the distribution of data and sampling are plotted in the respective figures below.

5.1.2 Model architecture details:↩︎

We parametrize the Hamiltonian ansatz with a neural network \(\mathcal{H}_{\theta}(\mathbf{q}, \mathbf{p})\) where \(\mathbf{q}\) and \(\mathbf{p}\) are one-dimensional vectors of canonical phase space coordinates. Our neural network model has configurable input and hidden layers, which are passed as parameters. The network can learn a Hamiltonian of any number of particles in any number of dimensions as long as we have \(dim(\mathbf{q}) = dim(\mathbf{p})\) as shown in Figure 5.

5.1.3 Hyperparameter details↩︎

A minibatch size of 512 was used for both the training and validation datasets. The ansatz neural network has a configurable number of layers and dimensions passed by the user as command line arguments. For the first three experiments, we used a single hidden-layer NN model with dimension 32. The input layer had dimension \(2d\), which again is configurable, where \(2d\) is the dimension of the phase space of system under study. Model parameters were optimized using the Adam optimizer with an initial learning rate of \(10^{-3}\) . A learning rate scheduler based on the ReduceLROnPlateau strategy was employed to adapt the learning rate during training in response to stagnation in validation loss. The HNN was trained for 10 epochs. For each batch, the starting point was sampled at a random timestep along the trajectory, allowing the model to observe a wide variety of initial conditions and improve generalization across the phase space. In the forward pass, we simulated trajectories for 2 timesteps and evaluated the loss according to 28 . Gradients were computed using the adjoint state method, employing the same integrator as used in the forward pass, which leverages the fact that a symplectic Runge–Kutta integrator is invariant under cotangent lift. The loss function is \[\label{eq:final95objective} Loss = \frac{1}{n}\sum\limits_{i=1}^{n} \|p^i_{T} - \tilde{p}^i_T\|_2^2 + \|q^i_T - \tilde{q}^i_T\|_2^2\tag{28}\] where the state variables \((q_t, p_t)\) are evolved using’ a symplectic integrator.

5.1.4 Evaluation Criteria↩︎

We aim to learn the functional form of the Hamiltonian \(\mathcal{H}(q, p)\) beyond simple trajectory matching. Unlike prior works such as Hamiltonian Neural Networks (HNN) and Neural ODEs, which primarily assess the learned models by comparing predicted and ground-truth trajectories, we explicitly evaluate the learned Hamiltonian function \(\mathcal{H}_{NN}\) across the broader phase space. To do so, we sample test points \((q^i, p^i)\) from three distinct distributions over the phase space:

• Random Uniform: Points are sampled independently from a uniform distribution over a fixed bounding box that encompasses the training trajectories. This probes generalization to randomly scattered unseen states.

• Uniform Square Grid: A structured grid of evenly spaced points is generated within the same bounding box. This enables a systematic and resolution-controlled evaluation of \(\mathcal{H}_{NN}\) over the phase space.

• Gaussian: Points are sampled from a multivariate Gaussian distribution centered around typical states observed during training.

The results for the Hamiltonian prediction error \(\|\mathcal{H}_{true} - \mathcal{H}_{NN}\|_1\) over the test domain are plotted in Figures 11, 12 and 13 for the three Hamiltonian systems under evaluation respectively. Table 1 provides a direct comparison of our method against two of the chosen baselines in terms of mean Hamiltonian prediction error over the test domain.

5.2 Hamiltonian Systems and Results↩︎

System 1: Double Well potential
The particle in a double-well potential is another commonly studied system in classical and quantum mechanics where the system has 2 stable fixed points. In our case, we consider a symmetrical double potential well with Hamiltonian and governing equations given by \[\label{eq:state95dw} \begin{align} \mathcal{H}&= \frac{p^2}{2} + \frac{q^4}{4} - \frac{q^2}{2}, \\ \dot{q}&= p,\quad \dot{p}= q - q^3. \end{align}\tag{29}\] The plots in Figure 6 show the input data distribution, training and validation loss, and the predicted dynamics for the double well system.

System 2: Coupled Harmonic Oscillator
A coupled harmonic oscillator is a simple 1-D non-separable Hamiltonian system with Hamiltonian and governing dynamics given by \[\label{eq:state95coupled95ho} \begin{align} \mathcal{H}=& \frac{p^2}{2} + \frac{q^2}{2} + \alpha pq,\\ \quad \dot{q}=& p + \alpha q, \quad \dot{p}= -(q + \alpha p). \end{align}\tag{30}\] The plots in Figure 7 show the input data distribution, training and validation loss, and the predicted dynamics for the coupled oscillator system.

System 3: Henon-Hieles Potential
We now explore higher-dimensional systems where chaos can emerge. Chaotic systems are deterministic yet unpredictable over long timescales due to exponential error growth, governed by the Lyapunov exponent. However, since these systems follow well-defined Hamiltonians, their dynamics can still be learned from limited observations. A key example is the Hénon–Heiles (HH) system, a Hamiltonian system describing a star’s planar motion around a galactic center. While the system exhibits chaotic behavior, stable regions exist, aiding in learning its governing dynamics [33]. The Hamiltonian \(\mathcal{H}\) and the corresponding equations of motion are given by \[\label{eq:state95hh} \begin{array}{ll} \mathcal{H} = \frac{p_x^2 + p_y^2}{2} + \frac{q_x^2 + q_y^2}{2} + q_x^2q_y -\frac{q_y^3}{3},\\ \dot{q_x} = p_x,\quad \dot{q_y} = p_y, \\ \dot{p_x} = -q_x - 2q_xq_y,\quad \dot{p_y} = -q_y - q_x^2 + q_y^2.\\ \end{array}\tag{31}\]

The plots in Figure 8 show the input data distribution, training and validation loss, and the predicted dynamics for Hénon–Heiles system.

System 4: Kepler’s potential

Kepler’s potential is a spherically-symmetric central potential given by \(V(r) = -\frac{k}{r}\), where \(r\) is the radial distance between interacting particles. The problem has a singularity at \(r=0\) but this is not a fundamental impediment to NN learning and symplectic integrators [34]. However, the trajectories exhibit finer scales and we have to significantly increase the sampling resolution to compensate for this. It represents the potential energy of a particle under gravitational or electrostatic attraction, resulting in elliptical, parabolic, or hyperbolic orbits. It describes Kepler’s laws for planetary motion. These paths, described by Kepler’s Laws, range from bound, closed elliptical orbits to unbounded, open parabolic or hyperbolic trajectories, depending on the system’s total energy as shown in Figure 9. The total energy of a particle is given by \[\begin{align} E_{tot} &= \frac{1}{2}m\dot{r}^2 + V_{eff}, \qquad \text{where} \qquad V_{eff} = \frac{J^2}{2mr^2} -\frac{k}{r}, \end{align}\] where \(J\) is the angular momentum and is a constant of motion.
We analyze the Hamiltonian of two interacting particles which, in Cartesian coordinates of position and momenta, is given by \[\label{eq:kepler95295body} \begin{align} \mathcal{H}&= \frac{1}{2m_1}(p_{x_1}^2 + p_{y_1}^2) + \frac{1}{2m_2}(p_{x_2}^2 + p_{y_2}^2) - \frac{G m_1m_2}{\sqrt{(q_{x_1} - q_{x_2})^2 + (q_{y_1} - q_{y_2})^2}},\\ \dot{q}_{x_1} &= \frac{p_{x_1}}{m_1},\quad \dot{q}_{y_1} = \frac{p_{y_1}}{m_1}, \quad \dot{q}_{x_2} = \frac{p_{x_2}}{m_2}, \quad \dot{q}_{y_2} = \frac{p_{y_2}}{m_2},\\ \dot{p}_{x_1}&=-\frac{Gm_1m_2(q_{x_1} - q_{x_2})}{r^3_{12}}, \quad \dot{p}_{y_1}=-\frac{Gm_1m_2(q_{y_1} - q_{y_2})}{r^3_{12}},\\ \dot{p}_{x_2}&=\frac{Gm_1m_2(q_{x_1} - q_{x_2})}{r^3_{12}}, \quad \dot{p}_{y_2}=\frac{Gm_1m_2(q_{y_1} - q_{y_2})}{r^3_{12}}. \end{align}\tag{32}\] In the equations above, \(q_{x_1}, p_{x_1}, q_{y_1}, p_{y_1}\) represent the phase space coordinates of particle 1 and similarly, \(q_{x_2}, p_{x_2}, q_{y_2}, p_{y_2}\) represent the phase space coordinates for particle 2.
Data generation Due to the dimensionality and functional form of the problem, we generated a larger set of training data for Kepler’s problem, where the training set consists of 32768 trajectories and the validation set consists of 8192 trajectories. To avoid sampling at the singularity, we sampled initial conditions such that \(r\sim\mathcal{N}(\mu=4, \sigma=2.5)\) where r is the radial distance between the bodies and \(\mathcal{N}(\mu, \sigma)\) represents a random normal distribution. The data was generated with a timestep of 0.01 sec with 100 total timesteps. After data generation, the trajectories were corrupted at each timestep with Gaussian noise sampled from \(\mathcal{N}(\mu=0, \sigma=10^{-3})\).

The results for the system are presented in Figure 10.

5.3 Runtime and Memory performance↩︎

In our study, we systematically analyzed the runtime and memory consumption of the adjoint method and backpropagation across increasing simulation lengths, ranging from 4 to 32 timesteps. The evaluation was conducted on a benchmark problem involving the training of a coupled harmonic oscillator system with a single batch of size 512. The results as shown in Figure 15 show a stark contrast in memory scalability between the two approaches: while the adjoint method maintains a constant memory footprint irrespective of the simulation length, the memory usage in backpropagation exhibits a linear growth pattern. This discrepancy arises due to the fundamental difference in how gradients are computed. Backpropagation explicitly stores intermediate states for every timestep, whereas the adjoint method reconstructs gradients via a reverse-time integration of the system dynamics, circumventing the need for extensive memory allocation. Interestingly, our runtime analysis also favors the adjoint method for smaller-scale problems, where it demonstrates superior computational efficiency, surpassing backpropagation in execution speed.

4pt

5.4 Modified Hamiltonian and correction of learned Hamiltonian↩︎

In this section, we will validate the claims in Section 3.3. Let us establish some terminology.

• \(\mathbf{\mathcal{H}_{true}}\): The true Hamiltonian that we wish to learn. Ideally, it should represent a true physical system from which our training trajectories \((q, p)\) are generated.

• \({\Phi_h}\): The one-step numerical map which we use in the forward pass to generate the trajectories.

• \(\mathbf{\mathcal{H}_{NN}}\): The learned Hamiltonian, represented by a neural network.

• \(\mathbf{\tilde{\mathcal{H}}_{NN}}\) : The modified Hamiltonian of the learned Hamiltonian, the exact flow of which agrees with the numerical flow of the learned Hamiltonian, i.e., \(\Phi_h^{\mathcal{H}_{NN}}=\varphi_h^{\tilde{\mathcal{H}}_{NN}}\).

If we use a neural network to train over the data \((q_i, p_i)\) to learn a Hamiltonian \(\mathcal{H}_{true}\). Let us assume that the data is noiseless and that the learning is perfect, then the symplectic integrator applied to the learned Hamiltonian \(\mathcal{H}_{NN}\) agrees with the exact flow of the true Hamiltonian \(\mathcal{H}_{true}\). By definition, the modified Hamiltonian \(\tilde{\mathcal{H}}_{NN}\) of the learned Hamiltonian \(\mathcal{H}_{NN}\) is such that the exact flow of the modified Hamiltonian \(\tilde{\mathcal{H}}_{NN}\) agrees with the symplectic integrator applied to the learned Hamiltonian \(\mathcal{H}_{NN}\). From this, we conclude that the modified Hamiltonian of the learned Hamiltonian is the true Hamiltonian that we set out to learn in the first place, as seen in Figure 16. Our neural network model aims to learn the Hamiltonian which, when passed through our one-step method \(\Phi_h\), generates the true flow.

Consider a simple Hamiltonian system of a harmonic oscillator given by \[\label{eq:harmonic95oscillator} \begin{align} \mathcal{H}_{true} = \frac{1}{2}(p^2 + q^2). \end{align}\tag{33}\] We use a high-order symplectic integrator (4th order in this case) to generate the trajectory samples \((\hat{q}, \hat{p})\). Now, for our learned Hamiltonian \(\mathcal{H}_{NN}\), we consider a very simple 2nd degree polynomial ansatz given by \[\begin{align} \mathcal{H}_{NN} = \theta_0 + \theta_1q + \theta_2p + \theta_3q^2 + \theta_4qp + \theta_5p^2. \end{align}\] We wish to learn the coefficients \(\theta_i\). We use a second-order symplectic map \(\Phi_h\) in the forward pass of our neural net, which generates the flow \((q_h, p_h)\). We use a stepsize \(h=0.1\) for the forward map. Now we formulate the loss function over the true and predicted flow and optimize according to the scheme described above.

5.4.1 Modified Hamiltonian↩︎

For the implicit midpoint rule, the modified Hamiltonian for the learned Hamiltonian \(\mathcal{H}_{NN}\) is given using 14 by \[\label{eq:NN95modified95Hamiltonian} \begin{align} \tilde{\mathcal{H}}_{NN} = \mathcal{H}_{NN} - \frac{h^2}{24} \nabla^2\mathcal{H}_{NN}(J^{-1}\nabla \mathcal{H}_{NN}, J^{-1}\nabla \mathcal{H}_{NN}) + \mathcal{O}(h^4), \end{align}\tag{34}\] where the term \(\nabla\mathcal{H}_{NN}\) is a vector of gradients w.r.t. \(q\) and \(p\). Now, from the training, the learned Hamiltonian is \[\begin{align} \mathcal{H}_{NN} = 0.15410 + (1.1111\times10^{-7})q -(6.1064\times10^{-8})p + 0.500420q^2 -(1.5463\times 10^{-8})qp + 0.500415p^2. \end{align}\] We can immediately observe that the coefficients of \(q^2\) and \(p^2\) are not the same as in 33 , also we note that there is a non-zero offset which is due to the fact that while training, our symplectic map 19 only sees the derivative of \(\mathcal{H}\) w.r.t. \(q\) and \(p\) and hence any Hamiltonian function plus a constant also produces the same dynamics. Neglecting the numerically insignificant terms and correcting for the offset, we can write this as \[\label{eq:learned95H} \begin{align} \mathcal{H}_{NN} = 0.500420q^2 + 0.500415p^2. \end{align}\tag{35}\] Now, we plug this into the expression for the modified Hamiltonian 34 , \[\begin{align} \tilde{\mathcal{H}} = 0.500420q^2 + 0.500415p^2 - \frac{(0.1)^2}{24} (1.00084q^2 + 1.00083q^2) + \mathcal{O}(h^4). \end{align}\] Neglecting the fourth-order terms, we obtain the modified Hamiltonian for our learned Hamiltonian, \[\label{learned95modified95H} \begin{align} \tilde{\mathcal{H}}_{NN} \approx 0.500004q^2 + 0.500001p^2, \end{align}\tag{36}\] where we have neglected the fourth-order terms and we observe that it agrees better with the true Hamiltonian, as seen in Figure 17. For the other examples we considered, it is possible to obtain a better approximation of the true Hamiltonian by post-processing the learned Hamiltonian \(\mathcal{H}_{NN}(q, p)\) using the modified Hamiltonian approach to obtain an improved approximation.

6 Conclusion↩︎

The adjoint approach normally results in gradients that differ from backpropagation, unless the adjoint system is computed using the cotangent lift of the numerical integrator used in the forward propagation, in which case the adjoint approach yields gradients that coincide with backpropagation. When the forward flow is integrated using a Runge–Kutta method, the cotangent lift is a symplectic partitioned Runge–Kutta method, where the original Runge–Kutta integrator to the primal variables, and the adjoint variables are integrated by the symplectic adjoint of the original Runge–Kutta integrator, defined by 16 .

When the forward flow is Hamiltonian, as is the case for Hamiltonian Neural Networks, it is natural to discretize the forward flow using a symplectic integrator. If the forward integrator is a symplectic Runge–Kutta method, then the symplectic adjoint is itself. As such, using the same symplectic integrator on the adjoint variables will lead to a discretization of the adjoint system that yields gradients that also coincide with backpropagation, leading to an efficient method for training Hamiltonian Neural Networks.

In our work, we adopt such an approach, using implicit symplectic partitioned Runge–Kutta methods. Symplectic methods are generally implicit for non-separable Hamiltonians, unless one artificially doubles the number of variables in an augmented formulation. However, contrary to conventional wisdom, implicit SPRK methods can be very efficiently implemented by using an explicit RK method of the same order as a predictor, and using a few fixed-point iterations of the SPRK method as a corrector. Therefore, Hamiltonian Neural Networks with a non-separable Hamiltonian ansatz can be efficiently trained using implicit SPRK discretization by applying the adjoint method combined with the predictor-corrector fixed-point iteration.

Finally, we observe that the learned Hamiltonian obtained using a lower-order symplectic integrator can be post-processed using backward error analysis to obtain a modified Hamiltonian that better agrees with the true Hamiltonian, without the need to use higher-order integrators in the learning process.

Acknowledgements↩︎

• HC and VK acknowledge support from the Czech National Science Foundation under Project 24-11664S

• HC acknowledges the computational resources provided by RCI project.

Disclaimer↩︎

This paper was prepared for information purposes and is not a product of HSBC Bank Plc. or its affiliates. Neither HSBC Bank Plc. nor any of its affiliates make any explicit or implied representation or warranty and none of them accept any liability in connection with this paper, including, but not limited to, the completeness, accuracy, reliability of information contained herein and the potential legal, compliance, tax or accounting effects thereof. Copyright HSBC Group 2025.

References↩︎

---

1. The action functional is extremized when the first variation vanishes; it may correspond to a minimum, maximum, or saddle point.↩︎

2. Strictly speaking, \(q_i\) should be written as \(q^i\), since it denotes a component of a vector, while \(p_i\) is a component of a covector. However, for notational simplicity and because we work primarily with individual components when formulating our loss functionals, we use subscripts for both.↩︎

3. We assume \(L\) is regular, meaning the Hessian matrix \(\big(\partial^2 L/\partial \dot{q}_i \partial \dot{q}_j\big)\) is nonsingular, so that this transformation is locally invertible. If this induces a global diffeomorphism \(TQ \cong T^*Q\), then the Lagrangian is said to be hyperregular. Under this assumption, for each \((q,p)\), there is a unique \(\dot{q}\) such that \(p_i = \partial L/\partial \dot{q}_i\).↩︎

4. An important point to note is that the continuous adjoint computes gradients of the underlying ODE, whereas backpropagation corresponds to the discrete adjoint (cotangent lift) of the numerical integrator; these coincide only in the limit of vanishing stepsize or when the adjoint discretization is the symplectic (cotangent) lift of the forward scheme.↩︎