Symplectic Recurrent Neural Networks
Year: 2020

Abstract:
We propose Symplectic Recurrent Neural Networks (SRNNs) as learning algorithms that capture the dynamics of physical systems from observed trajectories. An SRNN models the Hamiltonian function of the system by a neural network and furthermore leverages symplectic integration, multiplestep training and initial state optimization to address the challenging numerical issues associated with Hamiltonian systems. We show SRNNs succeed reliably on complex and noisy Hamiltonian systems. We also show how to augment the SRNN integration scheme in order to handle stiff dynamical systems such as bouncing billiards.

heading: 1 Introduction
text: Can machines learn physical laws from data? A recent paper (Greydanus et al., 2019), Hamiltonian Neural Networks (HNN), proposes to do so representing the Hamiltonian function H(q, p) as a multilayer neural network. The partial derivatives of this network are then trained to match the time derivatives ṗ and q̇ observed along the trajectories in state space. The ordinary differential equations (ODEs) that express Hamiltonian dynamics are famous for both their mathematical elegance and their challenges to numerical integration techniques. Except maybe for the simplest Hamiltonian systems, discretization errors and measurement noise lead to quickly diverging trajectories. In other words, Hamiltonian systems can often be stiff, a concept that usually refers to differential equations where we have to take very small time-steps of integration so that the numerical solution remain stable (Lambert, 1991). A plethora of numerical integration methods, symplectic integrators, have been developed to respect the conserved quantities in Hamiltonian systems, thereby usually being more stable and structure-preserving than non-symplectic ones (Hairer et al., 2002). For example, the simplest symplectic integrator is the well-known leapfrog method, also known as the Stömer-Verlet integrator (Leimkuhler and Reich, 2005). However, even the best integrators remain severely challenged by phenomena as intuitive as a mechanical rebound or a slingshot effect, which are more severe forms of stiffness. Such numerical issues are almost doomed to conflict with the inherently approximate nature of a learning algorithm. In the first part of this paper, we propose Symplectic Recurrent Neural Networks (SRNNs), where (i) the partial derivatives of the neural-network-parametrized Hamiltonian are integrated with the leapfrog integrator and where (ii) the loss is back-propagated through the ODE integration over multiple time steps. We find that in the presence of observation noise, SRNN are far more usable than HNNs. Further improvements are achieved by simultaneously optimizing the initial state and the Hamiltionian network, presenting an interesting contrast to previous literature on the hardness of general initial state optimization (Peifer and Timmer, 2007). The optimization can be motivated from a maximum likelihood estimation perspective, and we provide heuristic arguments for why the initial state optimization is likely convex given the symplecticness of the system. Furthermore, experiments in the three-body problem show that the SRNN-trained Hamiltonian compensates for discretization errors and can even outperform numerically solving the ODE using the true Hamiltonian and the same time-step size. This could be of particular interest to researchers who study the application of machine learning to numerically solving differential equations. The second part of this paper focuses on perfect rebound as an example of the more severe form of stiffness. When a point mass rebounds without loss of energy on a perfectly rigid
obstacle, the motion of the point mass is changed in ways that can be interpreted as an infinite force applied during an infinitesimal time. The precise timing of this event affects the trajectory of the point mass in ways that essentially make it impossible to merely simulate the Hamiltonian system on a predefined grid of time points. In order to address such events in learning, we augment the leapfrog integrator used in our SRNN with an additional trainable operator that models the rebound events and relates their occurrence to visual hints. Training such an augmented SRNN on observed trajectories not only learns the point mass dynamics but also learns the visual appearance of the obstacles.

heading: 2 Related work
text: Learning physics with neural networks A popular category of methods attempts to replicate the intuitive ways in which humans perceive simple physical interactions, identifying objects and learning how they relate to each other (Battaglia et al., 2016; Chang et al., 2016). Since such methods cannot be used for more general physical systems, another category of methods seeks to learn which differential equations govern the evolution of a physical system on the basis of observed trajectories. Brunton et al. (2016) assemble a small number of predefined primitives in order to find an algebraically simple solution. Lutter et al. (2019) use a neural network to model the Lagrangian function of a robotic system. Most closely related to ours, Greydanus et al. (2019) use a neural network to learn the Hamiltonian of the dynamical system in such a way that its partial derivatives match the time derivatives of the position and momentum variables, which are both assumed to be observed. Although the authors show success on a simple pendulum system, this approach does not perform well on a more complex system such as a three-body problem. ODE-based learning and recurrent neural networks (RNNs) To learn an ODE that underlies some observed time series data, Chen et al. (2018a) proposes to solve a neuralnetwork-parameterized ODE numerically and minimize the distance between the generated time series with the observed data. To save memory, they propose to use the adjoint ODE instead of back-propagating through the ODE solver. Using stability analysis of ODEs, Chang et al. (2019) propose the AntisymmetricRNN with better trainability. Niu et al. (2019) establish a correspondence between RNNs and ODEs, and propose an RNN architecture inspired by a universal quantum computation scheme. Summary of our main contributions: In this paper, we propose SRNN, which
• learns Hamiltonian dynamics directly from position and momentum time series • performs well on noisy and complex systems such as a spring-chain system and a
three-body system, and is compatible with initial state optimzation • is augmented to handle perfect rebound, an example of very stiff Hamiltonian dy-
namics

heading: 3 Framework

heading: 3.1 Hamiltonian systems
text: A Hamiltonian system of dimension d is described by two vectors p, q ∈ Rd. Typically, they correspond to the momentum and position variables, respectively. The evolution of the system is determined by the Hamiltonian function H : (p, q, t) ∈ R2d+1 7→ H(p, q, t) ∈ R through a system of ordinary differential equations called Hamilton’s equations,
ṗ = −∂H ∂q , q̇ = + ∂H ∂p , (1)
where we use the dot notation to compactly represent derivatives with respect to the time variable t. We are focusing in this work on Hamiltonians that are conservative,1 that is, they do not depend on the time variable t, and separable,2 that is, they can be written as a
1In our opinion, extending this work to non-conservative systems should not be done by adding a time dependency in the Hamiltonian, but by adding additional dissipation or intervention operators in the numerical integration schema, as illustrated in section 6.
2Extending this work to non-separable Hamiltonians can be achieved by rewriting the numerical integration schema using an extended phase space (Tao, 2016).
sum H(p, q) = K(p) + V (q). In this case, (1) becomes
ṗ = −V ′(q), q̇ = K ′(p) (2) With a proper choice of the p and q variables, the evolution of essentially all physical systems can be described with the Hamiltonian framework. In other words, Hamilton’s equations restrict the vast space of dynamical systems to the considerably smaller space of dynamical systems that are physically plausible. Therefore, instead of modeling the dynamics of a physical system with a neural network fθ(p, q) whose outputs are interpreted as estimates of the time derivatives ṗ and q̇, we can also use a neural network Hθ(p, q) = Kθ1(p) + Vθ2(q) with θ = [θ1, θ2], whose partial derivatives −V ′θ2(q) and K ′ θ1
(p) are interpreted as the time derivatives ṗ and q̇. We refer to the former as ODE neural networks (O-NET) and the latter approach as Hamiltonian neural networks (H-NET). In order to define a complete learning system, we need to explain how to determine the parameter θ of the neural networks on the basis of observed discrete trajectories. For instance, Greydanus et al. (2019) trains H-NET in a fully supervised manner using the observed tuples (p, q, ṗ, q̇).

heading: 3.2 From ODEs to discrete trajectories
text: A numerical integrator (or ODE solver) approximates the true solution of an ODE of the form ż = f(z, t) at discrete time steps t0, t1. . .tT . For instance, the simplest integrator, Euler’s integrator, starts from the initial state z0 at time t0 and estimates the function z(t) at uniformly spaced time points tn = t0 + n∆t with the recursive expression
zn+1 = zn + ∆t f(zn, tn) (3)
In stiff ODE systems, however, using Euler’s method could easily lead to unstable solutions unless the time-step is chosen to be very small (Lambert, 1991). The development of efficient and accurate numerical integrators is the object of considerable research (Hairer et al., 2008; Hairer and Wanner, 2013). Symplectic integrators3 are particularly attractive for the integration of Hamilton’s equations (Leimkuhler and Reich, 2005). They are able to preserve quadratic invariants, and therefore usually have desired stability properties as well as being structure-preserving (McLachlan et al., 2004), even for certain non-Hamiltonian systems (Chen et al., 2018b). A simple and widely-used symplectic integrator is the leapfrog integrator. When the Hamiltonian is conservative and separable (2), it computes successive estimates (pn, qn) with
pn+1/2 = pn − 12∆t V ′(qn)
qn+1 = qn + ∆tK ′(pn+1/2)
pn+1 = pn+1/2 − 12∆tV ′(qn+1)
(4)
Repeatedly executing update equations (4) is called the leapfrog algorithm, which is as computationally efficient as Euler’s method yet considerably more accurate when the ODE belongs to a Hamiltonian system (Leimkuhler and Reich, 2005).

heading: 3.3 Learning ODEs from discrete trajectories
text: Following Chen et al. (2018a), let the right hand side of the ODE be a parametric function fθ(z, t) and let z0 . . . zT be an observed trajectory measured at uniformly spaced time points t0 . . . tT . We can estimate the parameter θ that best represents the dynamics of the observed trajectory by minimizing the mean squared error ∑T i=1 ‖zi − ẑi(θ)‖2 between the observed trajectory {zi}Ti=0 and the trajectory {ẑi(θ)}Ti=0 generated with our integrator of choice,
{ẑi(θ)}Ti=0 = Integrator(z0, fθ, {ti}Ti=0) . For instance, this minimization can be achieved using stochastic gradient descent after back-propagating through the steps of our numerical integration algorithm of choice and then through each call to the functions fθ. This can be done when fθ(z) is a neural network (O-NET), or is the concatenation [−V ′θ2(q),K ′ θ1 (p)] of the partial derivatives of an H-NET
3An integrator is symplectic if applied to Hamiltonian systems, its flow maps are symplectic for short enough time-steps. For details, see Hairer et al. (2002) and Leimkuhler and Reich (2005).
Hθ(p, q) = Kθ1(p) + Vθ2(q), where the partial derivatives can be expressed using the same parameters θ as the Hamiltonian Hθ(p, q), for instance using automatic differentiation. We can then predict trajectories at testing time using the trained fθ∗ and initial state ztest0 ,
{ẑtesti }T test i=0 = Integrator(z test 0 , fθ∗ , {ti}T test i=0 ) .
Note that neither the integrator, nor the number of steps, nor the step size, need to be the same at training and testing.

heading: 3.4 Symplectic Recurrent Neural Network
text: This framework provides a number of nearly orthogonal design options for the construction of algorithms that model dynamical systems using trajectories:
• The time derivative model could be an O-NET or H-NET. • The training integrator can be any explicit integrators. In our experiments, we only
focus on Euler’s integrator and the leapfrog integrator. • The training trajectories can consist of a single step, T=1, or multiple steps, T>1.
We refer to the first case as single-step and the second case as multi-step or recurrent training, because back-propagating through multiple steps of the training integrator is comparable to back-propagating through time in recurrent networks. • The testing integrator can also be chosen freely and can use a different time-step size as it does not involve back-propagation.
In order to save space while describing the possibly different integrators used for training and testing, we use the labels “E-E”, “E-L”, and “L-L”, where the first letter tells which integrator was used for training —“E” for Euler and “L” for leapfrog— and the second letter indicates which integrator was used as testing time. For instance, with our terminology, the HNN model of Greydanus et al. (2019) is a “single-step E-E H-NET” with the additional subtlety that they supervise the training with actual derivatives instead of relying on finite differences between successive steps of the observed trajectories. A Symplectic Recurrent Neural Network (SRNN) is a recurrent H-NET that relies on a symplectic integrator for both training and testing, such as, for instance, a "recurrent L-L H-NET". As shown in the rest of this paper, SRNNs are far more usable and robust than the alternatives, especially when the Hamiltonian gets complex and the data gets noisy. We believe that SRNNs may also have other potential benefits: because leapfrog preserves volumes in the state space (Hairer et al., 2002), we conjecture that vanishing and exploding gradients’ issues in backpropagating through entire state sequences are ameliorated (Arjovsky et al., 2015). Finally, because the leapfrog integrator is reversible in time, there is no need to store states during the forward pass as they can be recomputed exactly during the backward pass. We leave studying these other computational and optimization benefits as a topic of future work.

heading: 4 SRNN can learn complex and noisy Hamiltonian dynamics
text: As an example of a complex Hamiltonian system, we first present experiments performed on the spring-chain system: a chain of 20 masses with neighbors connected via springs. Each of the two masses on the ends are connected to fixed ground via another spring. The chain can be assumed to lay horizontally and the masses move vertically but no gravity is assumed. The 20 masses and the 21 spring constants are chosen randomly and independently. The training data consist of 1000 trajectories of the same chain, each of which starts from a random initial state of positions and momenta of the masses and is 10-time-step long (including the initial state). We thus take T = 9 when performing recurrent training. When performing single-step training, each training trajectory of length 10 is instead considered as 9 consecutive trajectories of length 2. In this way, 1000 sample trajectories of length 10 (T=9) are turned into 9000 sample trajectories of length 2 (T=1), allowing for a fair comparison between single-step training and recurrent training. During testing, the trained model is given 32 random initial states in order to predict 32 trajectories of length 100. Detailed experiment setups and model architectures are provided in Appendix A.1, and a PyTorch implementation can be found at https://github.com/zhengdao-chen/SRNN.git.

heading: 4.1 Going symplectic - rescuing HNN with the leapfrog integrator
text: First, we consider the noiseless case, where the training data consist of exact values of the positions (q) and momenta (p) of the masses on the chain at each discrete time point. As shown in figure 1, the prediction of a single-step E-E H-NET deviates from the ground truth quickly and is unable to capture the periodic motion. By comparison, a single-step E-E ONET yields predictions that is qualitatively reasonable. This shows that using Hamiltonian models without paying attention to the integration scheme may not be a good idea. We then replace Euler’s integrator used during testing by a leapfrog integrator, yielding a Single-step E-L H-NET. Figure 1 shows that this helps the H-NET produce predictions that remain stable and periodic over a longer period of time. Since the training process remains the same, this implies that part of the instability and degeneration of H-NET’s predictions comes from the nature of Euler’s integrator rather than the lack of proper training. In contrast, using a leapfrog integrator for both training and testing substantially improve the performance, as also shown again in figure 1. This improvement shows the importance of consistency between the integrators used in training and predicting modes. This can be understood with the concept of modified equations (Hairer, 1994): when we use a numerical integrator to solve an ODE, the numerical solution usually does not strictly follow the original equation due to discretization, but can be regarded as a solution to a modified version of the original equation that depends on the integrator and the time-step size. Therefore, training and testing with the same numerical integrator and time-step size could allow the system to learn a modified Hamiltonian that corrects some of the errors caused by the discretization scheme.

heading: 4.2 Going recurrent - using multi-step training when noise is present
text: Since noise is prevalent in real-world observations, we also test our models on noisy trajectories. Independent and identically distributed Gaussian noise is added to both the position and the momentum variables at each time step. Applying the single-step methods described above yield considerably worse predictions, as shown in Figure 2 (left). This phenomenon can be controlled by training on multiple steps, effectively arriving at a type of recurrent neural network: if noise is added independently at each time-step, then having data from multiple consecutive time steps may allow us to discern the actual noiseless trajectory, analogous to performing linear regression on multiple (more than 2) noisy data points. As we see in Figure 2 (left), recurrent training consistently improves the predictions except for E-E H-NET. The best performing model is the SRNN (recurrent L-L H-Net) which improves substantially over the single-step L-L H-NET. Interestingly, the recurrent E-E H-NET does not improve over the single-step E-E H-NET, which means that recurrent training does not help if one uses a naïve integrator.

heading: 4.3 Initial state optimization (ISO)
text: However, one issue remains to be addressed: in the framework that we have adopted so far, the initial states p0 and q0 are treated as the actual initial states from which the system
begins to evolve despite the added noise in observation. With noise added to the observation of p0 and q0, our dynamical models will start from these noisy states and remain biased as we advance in time in both the training and the testing mode. To mitigate this issue, we propose to introduce two new parameter vectors for each sample, p̂0 and q̂0, interpreted as our estimate of the actual initial states, and we let our dynamical models evolve starting from them instead of the observed p0 and q0. Treating p̂0 and q̂0 as parameters, we can optimize them based on the loss function while fixing the model’s parameters, a process that we call initial state optimization (ISO). When the model is good enough, we hope that this will guide us towards the true initial states without observation noise. In Appendix B, we motivate the use of ISO from the perspective of maximum likelihood inference. In actual training, we first train the neural network parameters for 100 epochs as usual, and starting from the 101st, after every epoch we perform ISO with the L-BFGS-B algorithm (Zhu et al., 1997) on the p̂0 and q̂0 parameters for every training trajectory. At testing time, the model is given the noisy values of p and q for the first 10 time steps and must complete the trajectory for the next 200 steps. These 10 initial time steps allow us to perform the same L-BFGS-B optimization to determine the initial state p̂0 before advancing in time to predict the entire trajectory. As seen in Figure 2 (right), SRNN-ISO (i.e. SRNN equipped with ISO) clearly yields the best prediction among all the methods. Figure 3 shows the predictions of HNN, SRNN and SRNN-ISO on one test sample, and we clearly see the qualitative improvements thanks to recurrent training and ISO. O-NET also benefits from ISO while vanilla RNN and LSTM do not seem to, likely because the initial state optimization only works when we already have a reasonable model of the system. In Appendix C, we give a heuristic argument for the convexity of ISO, which helps to explain the success of using L-BFGS-B for ISO. In summary, we have proposed three extensions to learning complex and noisy dynamics with H-NET and demonstrated the improvements they lead to: a) using the leapfrog integrator instead of Euler’s integrator; b) using recurrent instead of single-step training; and c)
optimizing the initial states of each trajectory as parameters when data are noisy. Thorough comparisons of test errors are given in Tables 1 and 2, where we highlight that the SRNN (recurrent L-L H-NET) models achieve the lowest errors.

heading: 5 SRNN can learn the dynamics of a three-body system
text: Next, we test SRNN with the three-body system, which is a well-known example of a chaotic system, meaning that a small difference in the initial condition could lead to drastically different evolution trajectories, even without noise added. As a result, even when the exact equations are known, simulating it with different time-step sizes could also lead to qualitatively different solutions. Moreover, Greydanus et al. (2019) mentions that HNN does not outperform a baseline method using O-NET in learning the three-body system’s evolution. Here, we test our SRNN together with other baselines on the noiseless three-body system with the same configurations as Greydanus et al. (2019). The detailed experimental setup and model architectures are provided in Appendix A.2. As we see in Table 3, the best-performing model is SRNN and the second-best is the singlestep L-L H-NET. Interestingly, and perhaps counter-intuitively, they even outperform the baseline method of simulating the correct equation with the same time-step size. How is this possible? In short, our explanation is that the error introduced by numerical discretization could be learned and therefore compensated for by the models we train. More concretely, once again using the concept of modified equations mentioned in Section 4.1, we argue that the ODE-based learning models, including both H-NET and O-NET models, could learn not
the correct underlying equation, but rather the equation whose modified equation associated with our choice of numerical integrator and time-step size is the original equation. Hence, when the time-step size is large and the error of numerical discretization is not negligible, it is possible that the learned equation could yield better predictions than the correct one. In addition, we also see that the recurrently trained models outperform the corresponding single-step-trained models. Plots of the predicted trajectories are provided in Appendix E.

heading: 6 Learning perfect rebound with an augmented SRNN
text: We focus in this section on the perfect rebound problem as a prototypical example of stiff ODE in a physical system. We consider a heavy billiard, subject to gravitational forces pointing downwards, and bouncing around a two-dimensional square domain delimited by impenetrable walls. Whenever it hits a wall, the billiard rebounds without loss of energy, by reversing the component of its momentum orthogonal to the wall surface. Microscopically, when the billiard hits the wall, the atomic structure deformation produces strong electromagnetic forces that reverse the momentum during a very brief timescale. Simulating this microscopic phenomenon with a Hamiltonian ODE would not only be computationally expensive, but also require a detailed knowledge of the atomic structures of the billiard and the walls. The perfect rebound is a macroscopic approximation that treats the billiard as a point mass and the rebound as an event with zero duration infinite forces. Although this approximation is convenient for high-school level derivations, the singularity makes it hard to simulate using Hamiltonian dynamics. We propose to approach this problem by augmenting each time step of a leapfrog-based SRNN with an additional operation that models a possible rebound event,
ppostt ← p pre t − 2(p pre t · n)n , (5)
where ppret is the pre-rebound momentum vector and p post t is the post-rebound momentum vector. When the vector n is zero, this operation does not change the momentum in any way. When n is a unit vector orthogonal to a wall, this operation computes the momentum reversal that is characteristic of a perfect rebound. Vectors n of smaller length could also be used to model energy dissipation in manner that is reminiscent of the famous LSTM forget gate (Hochreiter and Schmidhuber, 1997). Because the billiard trajectory depends on the exact timing of the rebound event, we also need a scalar α ∈ [0, 1] that precisely places the rebound event at time t+α∆t between the successive time steps t and t+ ∆t. The augmented leapfrog schema then becomes
[ppret+α∆t, q pre t+α∆t] leapfrog←−−−−−− α∆t [pt, qt] (6)
ppostt+α∆t = p pre t+α∆t − 2(p pre t+α∆t · n)n (7)
[pt+∆t, qt+∆t] leapfrog←−−−−−− (1−α)∆t [ppostt+α∆t, q post t+α∆t] (8)
where equations (6) and (8) represent ordinary leapfrog updates (4) for time steps of respective durations α∆t and (1−α)∆t. More precisely, we first compute a tentative position q̃t+∆t and momentum p̃t+∆t assuming no rebound,
[p̃t+∆t, q̃t+α∆t] leapfrog←−−−−−−
∆t [pt, qt] , (9)
then compute both n and α as parametric functions of the tentative position q̃t+∆t as well as the current position qt, and finally apply the forward model (6–8). Note that the final state is equal to the tentative state when no rebound occurs, that is, when n = 0. Directly modeling n and α with a neural network taking q̃t+∆t as the input would be very inefficient because we would need to train with a lot of rebound events to precisely reveal the location of the walls. We chose instead to use visual cues in the form of a background image representing the walls. We model n as the product of a direction vector n̄ and a magnitude γ ∈ [0, 1], and we want the latter to take value close to 1 when perfect rebound actually occurs between t and t + ∆t and close to 0 otherwise. Both n̄ and α are modeled as MLPs that take as input two 10x10 neighborhoods of the background image, centered at positions qt and q̃t+∆t, respectively. In contrast, γ is modeled as an MLP that takes as input a smaller 2x2 neighborhood centered at q̃t+∆t and is trained with an additional regularization term ‖γ‖1 in order to switch the rebound module off when it is not needed. Training is achieved by back-propagating through the successive copies of the augmented leapfrog scheme, through the models of n̄, α, and γ, and also through the computation of the tentative p̃t+∆t and q̃t+∆t. We use 5000 training trajectories of length 10 starting from a randomly-sampled initial positions and velocities. Similarly, we use 32 testing trajectories of length 60. Detailed exprimental setup is included in Appendix A.3. Figure 5 plots some predicted and actual testing trajectories. Appendix F compares these results with the inferior results obtained with several baseline methods, including SRNN without the rebound module, and SRNN with a rebound module that does not learn α. One limitation of our method, however, results from the assumption that there is at most one rebound event per time step. Although this assumption fails when the billiard rebounds twice near a corner, as shown in the bottom right plot in Figure 5, our method still outperforms the baseline methods even in this case.

heading: 7 Conclusion
text: We propose the Symplectic Recurrent Neural Network, which learns the dynamics of Hamiltonian systems from data. Thanks to symplectic integration, multi-step training and initial state optimization, it outperforms previous methods in predicting the evolution of complex and noisy Hamiltonian systems, such as the spring-chain and the three-body systems. It can even outperform simulating with the exact equations, likely by learning to compensate for numerical discretization error. We further augment it to learn perfect rebound from data, opening up the possibility to handle stiff systems using ODE-based learning algorithms.

heading: Acknowledgments
text: The authors acknowledge stimulating discussions with Dan Roberts, Marylou Gabrié, Anna Klimovskaia, Yann Ollivier and Joan Bruna.

heading: A Experiment setup
text: A.1 The spring-chain experiment
We set ∆t = 0.1. The ground truth trajectories in both training and testing are simulated by the leapfrog integrator using ∆t′ = 0.001 and coarsened into time-grids of 0.1 with a factor of 100, since simulating with a much smaller time-step leads to much more accurate solution, which we will treat as the ground truth solution. The O-NET that represents fθ(p, q) is a one-hidden-layer MLP with 40 input units, 2048 hidden units and 40 output units. The H-NET that represents Hθ(p, q) = Kθ1(p) + Vθ2(q) consists of two one-hidden-layer MLPs, one for Kθ1 and the other for Vθ2 . Each of the MLPs have 20 input units, 2048 hidden units and 1 output unit. The vanilla RNN and LSTM models also have hidden states of size 2048. Implemented in PyTorch, the models are trained over 1000 epochs with the Adam optimizer (Kingma and Ba, 2014) with initial learning rate 0.001 and using the ReduceLROnPlateau scheduler4 with patience 15 and factor 0.7.
A.2 The three-body experiment
The ground truth trajectories are simulated by SciPy’s solve_ivp adaptive solver5 with method RK45. We coarse-grain the simulated ground truth trajectories into time-steps of ∆t = 1, so that the models developed in section 4 are numerically integrated with time-step ∆t = 1 in both training and testing. We intentionally set the time-step to be relatively large, so that it becomes interesting to compare these models with a baseline method of simulating the true equations with time-step ∆t = 1. In addition, the training data consist of 100 sample trajectories of length 10·∆t = 10, which are then turned into 900 trajectories of length 2 and 600 trajectories of length 5, respectively for single-step and recurrent training, in the same way as for the spring-chain experiments above. The O-NET that represents fθ(p, q) is a three-hidden-layer MLP with 12 input units, 512 hidden units in each hidden layer and 12 output units. The H-NET that represents Hθ(p, q) = Kθ1(p) + Vθ2(q) consists of two three-hidden-layer MLPs, one for Kθ1 and the other for Vθ2 . Each of the MLPs have 6 input units, 512 hidden units in each hidden layer and 1 output unit. The vanilla RNN and LSTM models also have hidden states of size 512. Implemented in PyTorch, the models are trained over 1000 epochs with the Adam optimizer with initial learning rate 0.0003 and using the ReduceLROnPlateau scheduler with patience 15 and factor 0.7.
A.3 The heavy billiard experiment
The full image has size 128x128 pixels. The thickness of the wall is 12 pixels on each of the four sides, which leaves the free space of size 104x104 pixels in the middle for the billiard to move within. The billiard has size 3x3 pixels. The O-NET that represents fθ(p, q) is a one-hidden-layer MLP with 4 input units, 32 hidden units and 4 output units. The H-NET that represents Hθ(p, q) = Kθ1(p) + Vθ2(q) consists of two one-hidden-layer MLPs, one for Kθ1 and the other for Vθ2 . Each of the MLPs have 2 input units, 32 hidden units and 1 output unit. The vanilla RNN model also has hidden states of size 32. For the rebound module, n̄ is computed as the normalized output of a two-hidden-layer MLP, with 200 input units, 128 units in the first hidden layer, 32 units in the second hidden layer and 2 output units. α is also computed using a two-hiddenlayer MLP, sharing the first hidden layer units with the MLP for n̄, and having 32 units in the second hidden layer and 1 output unit. γ is computed by passing through sigmoid the output of a two-hidden-layer MLP, with 4 input units, 16 units in each hidden layer and 1 output unit. All of the activation functions are tanh except for the hidden-to-output activation in the MLP for α, where ReLU is used. Implemented in PyTorch, the models are trained over 1500 epochs with the Adam optimizer with initial learning rate 0.005 and using the ExponentialLR scheduler6 with decay factor 0.99 until the learning rate reaches 0.0001. We set ∆t = 0.1, and use 5000 trajectories of length 10 ·∆t = 1 as training data, and 32 trajectories of length 60 ·∆t = 6 as testing data.
4https://pytorch.org/docs/stable/optim.html#torch.optim.lr_scheduler.ReduceLROnPlateau 5https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html 6https://pytorch.org/docs/stable/optim.html#torch.optim.lr_scheduler.ExponentialLR

heading: B The Maximum Likelihood Estimation perspective
text: In the presence of noise, we can interpret the learning problem described in section 3.3 above from the perspective of maximum likelihood inference, which also provides justification for treating the initial states as trainable parameters. We define models as follow:
ẑi(θ) = Integrator(ẑ0 = z0, fθ, {ti}Ti=0)
q(zi; θ) = 1√
(2π)dσ2d e−‖zi−ẑi(θ)‖ 2 2/(2σ 2)
P ({zi}ni=1|θ) = n∏ i=1 q(zi; θ) = 1 ( √ 2πσ2)nd n∏ i=1 e−‖zi−ẑi(θ)‖ 2 2/(2σ 2),
(10)
and L(θ|{zi}ni=1) = P ({zi}ni=1|θ) is the likelihood function given the time-series data {zi}ni=1. Note that this model assumes independence between zi and zj for i 6= j once θ is fixed. If we are to perform maximum likelihood inference, we arrive at the following:
max θ
: logL(θ|{zi}ni=1) = − nd 2 log ( 2πσ2 ) − 1 2σ2 n∑ i=1 ‖zi − ẑi(θ)‖22, (11)
which is equivalent to
min θ n∑ i=1 ‖zi − ẑi(θ)‖22 (12)
This provides a motivation for using the L2 loss, as we did in the experiments. So far, we consider θ as the only parameter of the model defined by equations 10, and therefore the only argument of the likelihood function, while ẑ0 is fixed to be the observed initial state z0. As a generalization, we can consider a strictly larger family of models by allowing z0 to vary as well. In this way, we treat both θ and z0 as the parameters in the model and therefore arguments of the likelihood function that we optimize on. In other words, the model becomes
ẑi(θ, ẑ0) = Integrator(ẑ0, fθ, {ti}Ti=0)
q(zi; θ, ẑ0) = 1√
(2π)dσ2d e−‖zi−ẑi(θ,ẑ0)‖ 2 2/(2σ 2)
P ({zi}ni=1|θ, ẑ0) = n∏ i=1 q(zi; θ) = 1 ( √ 2πσ2)nd n∏ i=1 e−‖zi−ẑi(θ)‖ 2 2/(2σ 2),
(13)
and the optimization problem becomes
min θ,ẑ0 n∑ i=1 ‖zi − ẑi(θ, ẑ0)‖22, (14)
which justifies optimizing over the initial states p0, q0 in addition to the neural network parameters θ as described in the previous section. Such an interpretation is similar to approaches for parameter estimation in the literature of inverse problems and systems biology, though in those cases the parameters of interest appear directly in ODEs instead of via neural networks (Peifer and Timmer, 2007; Stapor et al., 2018). In particular, jointly optimizing the parameters in the model as well as the initial value is called the initial value approach. However, despite the success we demonstrate in section 4.2, two difficulties of this approach have been pointed out: 1) The optimization could converge to local minima; 2) The numerical solution of the ODE can be unstable (Peifer and Timmer, 2007). As explained in section 4.1, using HNN together with the leapfrog integrator mitigates the second issue. But what about the first issue? In particular, even if we assume that the optimization of the neural network can work “magically” well and do not suffer from bag local minima, what about optimizing the initial value ẑ0?

heading: C Symplecticness and initial-state-optimization convexity
text: The success of optimizing on the initial state of the system in addition to the recurrent H-NET and O-NET models as described in section 4.2 raises the following question: If we
already have a relatively well-trained H-NET or O-NET, is the optimization on the initial values convex? We formalize the question below and provide a heuristic answer. For simplicity, we restrict our attention to autonomous ODEs, which means that the function f in dzdt = f(z) does not depend on t. Assuming existence and uniqueness of solutions, there exists a function φt that maps each initial state ẑ0 to the state of the system after evolving from ẑ0 for time t, φt(ẑ0). This function is usually called the flow map. Flow maps have also been defined for numerical solutions of ODEs, by letting φt(ẑ0) = Integrator(z0, f, {ti}Ti=0) with t0 = 0. We can extend this definition to all the trainable models we have considered, including the models based on O-NET and H-NET by defining φt(ẑ0) to be the state of the system after letting the system evolve from initial state ẑ0 for time t, for suitable choices of t. For example, for O-NET, we have φt(ẑ0) = Integrator(z0, fθ, {ti}Ti=0). Suppose we impose an L2 loss on φt(ẑ0), et(ẑ0) = ‖φt(ẑ0) − zt‖22, where zt corresponds to the observed data at time t. The question is, is et(z) a (perhaps locally) convex function of z, for what functions and numerical integrators? To understand convexity, we compute the gradient and the Hessian as follow.
∂
∂z et(z) = 2(φt(z)− zt)ᵀ · Ft(z) (15)
∂2
∂z2 et(z) = 2(φt(z)− zt) ·(3) Gt(z) + Ft(z)ᵀ · Ft(z), (16)
where Ft(z) is the Jacobian matrix of the flow map, defined as Ft(z)ij = ∂∂zj (φt(z)i), and Gt(z) is a third-order tensor contains the second order derivatives of the flow map, defined as Gt(z)ijk = ∂ 2
∂zizj (φt(z)k). We use ·(3) to denote the dot product in the third dimension.
Ft(z) ᵀ · Ft(z) is symmetric positive semidefinite for any matrix Ft(z). If φt corresponds to either the exact flow map of a Hamiltonian system or the flow of a symplectic integrator, such as the leapfrog integrator, applied to a Hamiltonian system, then Ft(z) is a symplectic matrix, implying that det(Ft(z)) = 1. Hence, det(Ft(z)ᵀ · Ft(z)) = 1, which further implies that Ft(z)ᵀ ·Ft(z) is a positive definite matrix. Therefore, non-rigorously, when ‖φt(z)−zt‖2 is small and so the first term on the right hand side of equation 16 is negligible compared to the least eigenvalue of Ft(z)ᵀ · Ft(z), the entire Hessian matrix ∂ 2
∂z2 et(z) is also positive definite, implying strong convexity of the optimization problem. If φt is the exact flow map, then ‖φt(z) − zt‖ being small means that noise in the data is small. If φt is the flow map of a learned model, then it means that we have a model close to the true underlying system in addition to not having too much noise in the data. Translating back to the learning problem, we see that, heuristically, when the model we use is close to symplectic, which is likely if the underlying system is a Hamiltonian system, and trained to be close enough to the true underlying system, and the noise in the data is small enough, then the optimization problem on the initial state is strongly convex.

heading: D Additional plots of the spring-chain experiments
text: D.1 Noiseless data (section 4.1)
D.2 Noisy data (section 4.2)

heading: E Additional plots of the three-body experiments
text: In each of the plots below, the three dashdot curves represent the ground truth trajectories of the three masses, and the three sequences of dots are the predictions made by each method.

heading: F Additional plots of the heavy billiard experiment

Review summaries:
This paper proposes to represent a Hamiltonian model of a physical system by a neural network. The parameters are then adjusted, so that the observations are considered maximally likely under a probabilistic model. The novelty is to consider a symplectic Leapfrog integration scheme for the Hamiltonian system, which is known to conserve important quantities such as volume in the state space. The proposed approach is shown to outperform the recent work "Hamiltonian Neural Networks" by a large margin on mass-spring chain dynamics and three body systems. The approach can even handle stiff dynamical systems such as bouncing billiards.