Symplectic Recurrent Neural Networks
Authors: Zhengdao Chen
Jianyu Zhang
Martin Arjovsky
Léon Bottou
Venue:  (2020)
Keywords: Hamiltonian systems
learning physical laws
symplectic integrators
recurrent neural networks
inverse problems

Abstract: We propose Symplectic Recurrent Neural Networks (SRNNs) as learning algorithms that capture the dynamics of physical systems from observed trajectories. SRNNs model the Hamiltonian function of the system by a neural networks, and leverage symplectic integration, multiple-step 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. Finally, we show how to augment the SRNN integration scheme in order to handle stiff dynamical systems such as bouncing billiards.

Published as a conference paper at ICLR 2020
SYMPLECTIC RECURRENT NEURAL NETWORKS
Zhengdao C Jianyu Zhang
Martin Arjovsky", Lêon Bottou
New York University, New York, USA
h Tianjin University, Tianjin, China
c Facebooo AI Research, New York, USA
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, multiple- step training and initial state optimization to address the challenging nu- merical issues associated with Hamiltonian systems. We show SRNNs suc- ceed 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.
INTRODUCTION
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 func- tion H(g.)) as multilayer neural network. The partial derivatives of this network are then trained to match the time derivatives ý 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 tech- niques. Except maybe for the simplest Hamiltonian systems, discretization errors and mea- surement 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 (Lam- bert, 1991). A plethora numerical inlegration methods, symplectic inteyrators, 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- integrator (Leimkuhler and Reich, 2005). However, even the best inte- grators remain severely challenged by phenomena intuitive as mechanical rebound slingshot effect, which are more severe forms of stiffness. Such numorical issucs 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- paramerriee Hamiltonian are in- legrated 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 observa- tion noise, SRNN are far more usable than HNNs. Further improvements are achieved by simultaneously optimizing the initial state and the Hamiltionian network, presenting an in- teresting contrast to previous literature on the hardness of general initial state optimization (Peifer and Timmer, 2007). The optimization can be motivated from a maximum likeli- hood estimation perspective, and we provide heuristic arguments for why the initial state optimization is likely convex given the symplectieness of the system. Furthermore, experi- ments 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 stifness. When a point mass rebounds without loss of energy on perfectly rigid

---

Published as a conference paper at ICLR 2020
obstacle, the motion of the point mass is changed in ways that be interpreted 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 the Hamiltonian system on a predefined grid of time points. In order to address such events in learning, augment the leapfrog integrator used in SRNN with 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.
merely simulate
2 RELATED WORK
Learning physics with neural networks A popular category of methods attempts to replicate the intuitive ways in which humans perceive simple physical interactions, iden- tifying objects and learning how they relate each other (Battaglia et al., 2016; Chang et al., 2016). Since such methods cannot be used for more general physical systems, an- other category 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 neural network to model the Lagrangian function of robotic system. Most closely related to ours, Greydamus et al. (2019) use a neural net to learn the Hamiltonian of the dynamical system in such 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 simple pendulum system, this approach docs not perform well on more complex system such as 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 neural- network-parameterize ODE numerically and minimize the distance between the generated time scries with the observed data. To ODE instead of back-propagating through the ODE solver. Using stability analysis of ODEs, Chang el. al. (2019) propose the ANiisymme 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.
memory, they propose to the adjoint
Summary of our main contribulions: In this paper, propose SRNN, which
learns Hamiltonian dynamics directly from position and momentum time series
performs well on noisy and complex systems such spring-chain system and
three-body system, and compatible with initial state optimzation
is augmented to handle perfect rebound, an example of very stiff Hamiltonian dy- namics
FRAMEWORK
3.1 HAMILTONIAN SYSTEMS
A Hamiltonian system of dimension d is described by two vectors p.q E Rd. Typically, they correspond to the momentum and position variables, respectively. The evolution of the system is determined by the Hamiltonian function H (p,g, () E R2d+1 + H(p, e R through system of ordinary differential equations called Hamillon's equations,
HH 0g
+ 0H
q
(1)
0p
where 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 the time variable f, and separable," that is, they can be written
"I opinion, extending this work to non-conservative systems should not be done by adding time dependency in the Hamiltonian, but by adding additional dissipation or intervention operators in the numerical integration schema, as illustrated in section "Extending this work to non-separable Hamiltonians can be achieved by rewriting the numerical integration schema using an extended phase space (Tao, 2016).
2

---

Published as a conference paper at ICLR 2020
sum H(p.)) (p) + V(g). In this case, (1) becomes
p= -==) i = k)))
(2)
With proper choice of the p and 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 fe(p, q) whose outputs are interpreted as estimates of the time derivatives and á. we can also use a neural network Ho(p. q) = Ks, (p) + Vo,(g) with [01. 02). whose partial derivatives - -V (q) and K%, (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 I 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, Q).
3.2 FROM ODES TO DISCRETE TRAJECTORIES
A numerical integrator (or ODE solver) approximates the true solution of an ODE of the form & = f(z,t) t) at discrete time steps tor tT. For instance, the simplest integrator, Euler's integrator, starts from the initial state 20 at time to and estimates the function =(t)
at uniformly spaced time points t = to + n At with the recursive expression
2n+1 2n + At f( tx)
(3)
In stiff ODE systems, however, using Euler's '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 the object of considerable research (Hairer et al., 2008; Hairer and Wanner, 2013). Symplectie integrators" 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 well as being structure- perectuurrpreervingg (McLachlan et al., 2004), even for certain non-
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 (Pa, with
Pn+1/2 Pn At (g.))
dut1 9n At K' (pn-1/2)
Px+1 Pn+1/2 1 (
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 Hamiltonian system (Leimkuhler and Reich, 2005).
3.3 LEARNING ODEs FROM DISCRETE TRAJECTORIES
Following Chen el. al. (2018a), let the right hand side of the ODE be parametric function fe(z, t) and let 2T. be an observed trajectory measured at uniformly spaced time points tt. We can estimate the parameter 0 that best represents the dynamics of the observed trajectory by minimizing the mean squared error Σ || 2;(6)ll2 between the observed trajectory {=){oo and the trajectory (0)) generated with our integrator of choice,
(2)(0) Integrator fer (ta)To))
For instance, this minimization can be achieved using stochastic gradient descent after back-propaggting through the steps of our numerical integration algorithm of choice and then through each call to the functions fe. This can be done when fe(2) is a neural network (O-NET), or the concatenation (-V,, (g), K2, (p)) of the partial derivatives of an H-NET
3 An integrator is symplectic if applied to Hamiltonian systems, its flow maps symplectic for short enough time-steps. For details, Hairer al. (2002) and Leimkuhler and Reich (2005).
3

---

Published as a conference paper at ICLR 2020
Ho(p. = Kø, (p) + Vo, (9), where the partial derivatives can be expressed using the same parameters 0 as the Hamiltonian (p, q), for instance using automatic differentiation. We can then predict trajectories at testing time using the trained fo- and initial state atsst
T
Integrator( a0 for (tel=m
Note that neither the integrator, nor the mumber of steps, nor the stop size, need to be the same at training and testing.
3.4 SYMPLECTIC RECURRENT NEURAL NETWORK
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 integralor can be any explieit integrators. In our experiments, we only focus on Euler's s integrator and the leapfrog integrator.
The training trajectories can consist of 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 recurrent training, because prok-prrogggting through multiple steps of the training integrator is comparable to paok-propggg through time in recurrent networks.
The testing integrator can also be chosen freely and can use a different time-step size it does not involve poee-propagation.
In order to save space while describing the possibly different integrators used for training and lesting, 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.
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, "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 al., 2002), we conjecture that vanishing and exploding gradients' issues in backpropagating through entire state sequences are ameliorated (Ar- jovsky et al., 2015). Finally, because the leapfrog integrator reversible in time, there 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 topic of future work.
SRNN CAN LEARN COMPLEX AND NOISY HAMILTONIAN DYNAMICS
As an example of a complex Hamiltonian system, we first present experiments performed on the spring-chain system: 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 gravily 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 7 when performing recurrent training. When performing single-step training, each training trajectory of length 10 is instcad considered as consecutive trajectories of length 10 (T-9) are turned into 9000 sample trajectories of length 2 (T-1), allowing for 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 PyTorc implementation can be found at https: //tgtgggttttttt//////ee//t com/zhengdao- saataaao-hheee/Seee chen/SRNN git.
In this way, 1000 sample trajectories of length
4

---

Published as a conference paper at ICLR 2020
10
Single-sepp
(HNN) H-NET
(HNN) Single E-E NET
Single-step F-I H-NET
Sinli-stapp EL D-NET
Single step LLL H-NET
13 erroro
Singiestepp
stope L-L NET
position
Observations
Single-3lpp L. H-NT
10
50
time step
me step
Figure 1: Testing results in the noiseless case by single-step methods. Left: Prediction error of each method over time, measured by the L2 distance between the true and predicted positions of the 20 masses. Right: Each curve represents the position of one of the masses (number 5) as a function of time prodicted by the three singlo-step-trained H-NET models. Plots of the other masses' positions are provided in Appendix D.1.
4.1 GOING SYMPLECTIC RESCUING HNN WITH THE LEAPFROG INTEGRATOR
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 a single-step E-E H-NET deviates from the ground truth quickly and is unable to capture the periodic motion. By comparison, single-step E-E O- NET yields predictions that is qualitatively reasonable. This shows that using Hamiltonian models without paying attention to the integration scheme may not be good idea.
We then replace Euler's integrator used during testing by leapfrog integrator, yielding Single-step F-T. H-NET. Figure shows that this helps the H-NET produce predictions that remain stable and periodic over longer period of time. Since the training process remains the same, this implies that part of the instability and degeneration of H-NET'ss 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 modified Hamiltonian that corrects some of the errors caused by the discretization scheme.
4.2 GOING RECURRENT USING MULTI-STEP TRAINING WHEN NOISE IS PRESENT
Since noise is prevalent in real-world observations, we also test our models on noisy trajec- tories. Independent and identically distributed Gaussian noise is added to both the position and the momentum variables at cach time step. Applying the single-step methods described above yield considerably worse prodictions, as shown in Figure (left).
This phenomenon can be controlled by training type of recurrent neural network: if noise is added independently at cach time-step, then having data from multiple consccutive time steps may allow to discern the actual noiscless 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-F: H-NET, which means that recurrent training does not help if one uses a naïve integrator.
multiple steps, effectively arriving at
4.3 INITIAL STATE OPTIMIZATION (ISO)
However, one issue remains to be addressed: in the framework that we have adopted so far, the initial states Po and go are treated the actual initial states from which the system
5

---

Published as a conference paper at ICLR 2020
10
Single-sepp
Recurient O-NET
(HNN) Siople-sipp E-E NET
Recurrent EE H-NET
Singil-Slpp L- O-NET
Recurrent L-L O-NET
errorr
Singlestepp
Recurrent
Recurrent E:E: O-NET
Recurrent EE O-NET I50
Recurent E-E H-NET
Recurrent E-E H-NET IS0
L2
(SRNN) Returrent O-NET
(SANN-1S)) Recurrent H-NET w/ 150
Vanilla RNN
Vanilla ANK w/I0
LSTH IS0
10
50
29 time step
30
time step
Figure 2: Prediction error of all methods in the noisy case measured by L2 distance, pre- sented in two plots due to the large number of methods. Included in the left plot are the single-ttep-trained methods, recurrently trained methods, vanilla RNN and LSTM. Included in the right plot are the (same) recurrently trained methods, the recurrently trained methods with initial state optimization (ISO), as well as vanilla RNN and LSTM with ISO.
HNN
SRNN
position Mmm
M
mmw
Mmm
SANN-I5O
Noiseless dynamic
position
position
Noisy observations
-2
40
time step
time step
time step
(a) HNN
(b) SRNN
(c) SRNN-ISO
Figure 3: Predictions made by three methods in the noisy case. The Y-axis corresponds to the position of one of the masses (number 5) on the chain.
begins lo evolve despite the added noise in observation. With noise added to the observation of po and 2o; 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 roduce two new parameter vectors for each sample, Pn and 2o: interpreted as our estimate of the actual initial states, and we let our dynamical models evolve starting from them instead of the observed po and qo- Treating Po and do as parameters,
can optimize them based on the loss function while fixing the model's paramoters, 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 observa- tion 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 Po and qo parameters for every training trajectory. At testing time, the model is given the noisy values of and 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 fo before advancing in time to predict the entire trajectory.
As seen in Figure (right), SRNN-ISO (i.c. 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 T-BFGG-B for ISO.
In summary, have proposed three extensions to learning complex and noisy dynamics with H-NET and demonstrated the improvements they lead to: a) using the leapfrog inte- grator instead of Fuler's integrator; b) using recurrent instead of single-step training; and c)

---

Published as a conference paper at ICLR 2020
Table 1: Testing results of predicting the dynamics of the spring-chain system by methods based on fixed Po, Qo (i.c., not optimizing Po. go as paramcters). The error is defined as the discrepancy between the (noisy) ground truth and the predictions at cach time stop averaged over the first 200 time steps, where the discrepancy is measured by the L2 distance between the true and predicted positions of the 20 masses in the chain, both of which considered as 20-dimensional vectors. The mean and standard devial are computed based on 32 testing samples, each starting from a random configuration of the chain.
Model
Integrator (tr)
Integrator (te)
Error mean
Error
single-step
O-NET
Euler
Euler
6.93
1.22
Euler
Leapfrog
5.87
1.04
Leaptrog
Leaptrog
7.28
1.48
H-NET
Euler
Euler
7.24
0.64
Euler
Leapfrog
3.32
0.89
Leaplrog
Leaptrog
3.36
0.67
recurrent
O-NET
Euler
Euler
2.88
0.45
Euler
Leapfrog
4.12
0.41
Leapfrog
Leapfrog
3.34
0.86
H-NET
Enler
Euler
7.58
0.63
Euler
Leapfrog
5.26
0.63
Leaptrog
Leapirog
2.37
0.87
Vanilla RNN
N/A
N/A
4.80
0.82
LSTM
N/A
N/A
5.95
1.05
Table Testing results of predicting the dynamics of the spring-chain system by methods that optimize on and starting from their observed (noisy) values using L-BFGS-B, as explained in the text. The definition of the errors the same in the above table.
Model O-NET
Integrator (Ir)
Integrator (te)
Error
Error std
Euler
Euler
2.13
0.37
Euler
Leapfrog
3.59
0.50
Leuplrog
Leapfrog
2.27
0.60
H-NET
Fuler
Euler
6.26
0.60
Euler
Leapfrog
3.00
0.63
Leaplrog
Leapirog
1.45
0.32
Vanilla RNN
N/A
N/A
4.72
0.94
LSTM
N/A
N/A
5.81
0.98
optimizing the initial states of each Irajectory as parameters when data 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.
SRNN CAN LEARN THE DYNAMICS OF A THREE-BODY SYSTEM
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 dif- ferent 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 qualita- tively 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. the noiseless three-body system
Here, we test SRNN together with other baselines 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 single- they even outperform the
step L-L H-NET. Interestingly, and perhaps counter-intuitieel time-step size. How is this baseline method of simulating the correct equation with the same possible? In short, our explanation is that the error introduced by numerical discretization could be learned and therefore compensated for by the models train. More concretely, once again using the concept of modified equations mentioned Section 4.1, we argue that the ODE-based learning models, including both H-NET and O-NET models, could learn not

---

Published as a conference paper at ICLR 2020
Table Prediction error results for the three-body system with time-step At 1. The last row corresponds to numerically solving the correct underlying equations using the leapfrog integrator with time-step At 1. The other rows correspond to the different learning- based methods, same as in the spring-chain experiments.
Model
Integrator (tx)
Integrator (te)
Error mean
Error std
single-step
O-NET
Euler
Buler
0.65
0.16
Euler
Leaplrog
1.36
0.18
Leapfrog
Leaplrog
1.33
0.20
H-NET
Euler
Euler
1.64
0.25
Euler
Leaptrog
0.88
0.33
Leapfrog
Leapfrog
0.35
0.09
recurrent
O-NET
Euler
Euler
0.51
0.11
Euler
Leapirog
1.27
0.18
Leapfrog
Leapfrog
0.49
0.10
H-NET
Euler
Euler
0.79
0.17
Euler
Leaptrog
1.76
0.62
Leapfrog
Leaptrog
0.26
0.07
simulation
true equs.
(no training)
Leapfrog
0.47
0.18
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.
LEARNING PERFECT REBOUND WITH AN AUGMENTED SRNN
We focus in this section on the perfect rebound problem as a prototypical example of stiff ODE in physical system. We consider heavy billiard, subject to gravitational forces pointing downwards, and bouncing around a two-dimensional square domain delimited by impenctrable walls. Whenever it hits wall, the billiard rebounds without loss of energy, by reversing the component its momentum orthogonal to the wall surface. Microscopically, when the billiard hits the wall, the atomic structure deformation produces strong electro- magnetic forces that reverse the momentum during a very brief timescale. Simulating this microscopic phenomenon with a Hamiltonian ODE would not only be computationally ex- pensive, 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 event. with zero duration infinite forces. Although this approximation is convenient for high-school level derivations, the singularity makes 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 possible rebound event,
poott t pf 2(r .n)
(5)
where pt is the pre-rebound momentum vector and 1t 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 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, also need a scalar a E [0, 1] that preciscly places the rebound event at time + aAt between the successive time steps and l + At. The augmented leapfrog schema then becomes
W are detacel lcopfrog
[pt, 4t]
(6)
aAt
Pttat Pttast 2(p n)n
(7)
8

---

Published as a conference paper at ICLR 2020
Actual Preditted
Attal Prodicted
10-
05-
0.-
-0.5
-0.
-10
-1.0
-0.5
0.
05
10
-1.9
-05
0.0
0.5
10
Actual Predicted
Adual Predicted
05-
05-
-05
-05
-10-
-10
0.0
0.5
-1.0
-0.5
0.0
1.0
Figure 4: Actual versus predicted trajectorics the heavy billiard with perfect rebound. The predictions are obtained by an SRNN plus the rebound module described in section 6.
leapfrog
Pe+At: T+41)]
Wetaser prast
(1 (1-a)At
where equations (6) and (8) represent ordinary leapfrog updates (4) for time steps of re- spective durations a&t and (1t More precisely, we first compute a tentative position l+t and momentum PitAr assuming no rebound,
(et leapfrog
[pe, qel
(9)
At
then compute both and a parametric functions of the tentative position t+ as well as the current position Gt, 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 a with neural network taking ]e+ar as the input would be very inefficient because we would need to train with 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 and magnitude E [0, 1], and we want the latter to take value close to I when perfect rebound actually occurs between and + At and close 0 otherwise. Both and a are modcled as MLPs that take as input two 10x10 neighborhoods of the background image, centered at positions ge and Tea respectively. In contrast, ? is modeled as an MLP that takes as input smaller 2x2 neighborhood centered al π++ and is trained with an additional not needed.
regularization term 1/// in order to switch the rebound module off when
Training is achieved by back-propagating through the successive copies of the augmented leapfrog scheme, through the models of n, a, and and also through the computation of the tentative Pi+Ae and qi+Ae. 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 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 rebound module that does not learn a. 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 in this case.

---

Published as a conference paper at ICLR 2020
7 CONCLUSION
We propose the Symplectie Recurrent Neural Network, which learns the dynamics of Hamil- tonian systems from data. Thanks to symplectic integration, milti-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. can perporf 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 handle stiff systems using ODE-based learning algorithms.
ACKNOWLEDGMENTS
The authors acknowledge stimulating discussions with Dan Roberts, Marylou Gabrić, Anna Klimovskaia, Yann Ollivier and Joan Bruna.
REFERENCES
Arjovsky, M., Shah, A., and Bengio, Y. (2015). Unitary evolution recurrent neural networks. CoRR, abs/ 1bs/1511.000 11141.06444
Battaglia, P. W., Pascanu, R., Lai, M., Rezende, D. J., and Kavukcuoglu, K. (2016). Inter- action networks for learning about objects, relations and physics. CoRR, abs 1612.00222.
Brunton, S. L.., Proctor, J. L., and Kutz, J. N. (2016). Discovering governing equations from data by sparse identification nonlincar dynamical systems. Proceedings of the National Academy of Sciences, 113:33)33332 3937.
Chang, B., Chen, M., Haber, E., and Chi, E. H. (2019). AntisymmetricRNN: A dynami- cal system view on recurrent neural networks. In International Conference on Learning Representations.
Chang, M. B., Ullman, T., Torralba, A., and Tenenbaum, J. B. (2016). compositional object-basee based approach to learning physical dynamics. CoRR, abs/ 1612.00341.
Chen, T. Q., Rubanova, Y., Bettencourl, J., and Duvenaud, D. K. (2018a). Neural ordinary differential equations. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa- Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems 91, pages 6571-6583. Curran Associates, Inc.
numerical integrators for Greydamus, S., Dzamba, M., and Yosinski, J. (2019). Hamiltonian neural networks. arXiv preprint arXiv: 01563.
Chen, Z., Raman, B., and Stern, A. (2018b). hodgkin--huxley-type systems.
Hairer, F. (1994). Backward analysis of numerical integrators and symplectic methods. Annals of Numerical Mathematics, 1:107-132. ID: un4ge:
Hairer, E., Lubich, C., and Wanner, G. (2002). Geomelric Numerical Integration: Structure- Preserving Algoritlems for Ordinary Differential Equalions. Springer series computa- tional mathematics. Springer.
Hairer, E., Nørsett, S. P., and Wanner, G. (2008). Solving Ordinary Differential Equations I: Nonstiff Problems. Springer Series Computational Matherattccs. Springer Berlin Heidelberg.
Hairer, E. and Wanner, G. (2013). Solving Ordinary Differential Equations II: Stiff and Dif- ferential - Algebracc Problems. Springer Scries Computational Mathematics. Springer Berlin Heidelberg.
Hochreiter, S. and Schmidhuber, J. (1997). Long short-term memory. Neural Computation, 7735
Kingma, D. P. and Ba, J. (2014). Adam: method for stochastic optimization. arXiv preprint 42
10

---

Published as a conference paper at ICLR 2020
Lambert, J. D. (1991). Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. John Wiley & Sons, Inc., New York, NY, USA.
Leimkuhler, B. and Reich, S. (2005). Simmlaaing Hamiltonian Dynamics. Cambridge Mono- graphs on Applied and Computational Mathematics. Cambridge University Press.
Lutter, M., Ritter, C., and Peters, J. (2019). Deep lagrangian networks: Using physics as model prior for deep learning. In International Conference on Learning Representations.
McLachlan, R. I., Perlmutter, M., and Quispel, G. R. W. (2004). On the nonlinear stability of symplectic integrators. BIT Numericol Mothematics, 44(1):99-117,
Niu, M. Y., Horesh, L.., and Chuang, (2019). Recurrent neural networks in the eye of differential equations. Xiu preprint urXiv:1904.1200
Peifer, M. and Timmer, (2007). Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shooting. The Institution of En- gineering and Technology, Systems Biology-
Stapor, P., Fröhlich, F., and Hasenauer, (2018). Optimization and profile calculation of ODE models using second order adjoint sensitivity analysis. Bioinformatics, 11(13) i159.
Tao, M. (2016). Explicit symplectic approximation of nonseparable hamiltonians: Algorithm and long time performance. Phys. Rev. E, 94:043303.
Zhu, C., Byrd, H., Lu, P., and Nocedal, (1997). Algorithm 778: L-bffgsb: Fortran
subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw., 23(4):550 560.
11

---

Published as a conference paper at ICLR 2020
A EXPERIMENT SETUP
A.1 THE SPRING-CHAIN EXPERIMENT
We set At = 0.1. The ground truth trajectories in both training and testing are simulated by the leapfrog integrator using At' = 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 fo(p, q) is a one-hidden-layer MLP with 40 input units, 2048 hidden units and 40 output units. The H-NET that represents Hp = Ko, (p) + Vo, (g) consists of two one- mee-hidden-layer MLPs, one for Ke, and the other for Va,, 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 Reduce LRRn Plateau scheduler" with patience 15 and factor 0.7.
A.2 THE THREE- THREE-BODY BODY EXPERIMENT
The ground truth trajectories are simulated by SciPy's solve_ivp adaptive solver with method RK45. We coarse-grain the simulated ground truth trajectories into time-steps of Al --=1 1, so that the models developed in section 4 are numerically integrated with time-step Al = in both training and lesting. 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 At = 1. In addition, the training data consist of 100 sample trajectories of length 10-Ar 10, which are then turned into 900 trajectorics 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 go(p,g) is three-hiddeenllyerr MLP with 12 input units, 512 hidden units in each hidden layer and 12 output units. The H-NET that represents He(p, q) Ke, (p) + Vo, (q) consists of two turee-hidden--lyee MLPs, one for Kg, and the other for Vo2- Each of the MLPs have 6 input units, 512 hidden units in each hidden layer and 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 15 and factor 0.7.
ROnPlateau scheduler with patience
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 fo(p, q) is one-hidden--lyyr MLP with 4 input units, 32 hidden units and 4 output. units. The H-NET that. represents Ho 4) = Kø, (p) + Vo, (q) consists of two one- aaedhiddenr MLPs, one for Ke, and the other for Ve. Each of the MLPs have 2 input units, 32 hidden units and 1 output unit. The vanilla RŇN model also has hidden states of size 32. For the rebound module, two-hidden-layer MLP, with 200 input units, 128 units in the first hidden layer, 32 units in the second hidden layer and output units. a is also computed using two-hidden- layer MLP, sharing the first hidden layer units with the MLP for n, and having 32 units in the second hidden layer and output unit. y is computed by passing through sigmoid the output of a two- two-hidden-lage and output unit. All of the activation functions are tanh except for the hidden-to-outpe activation in the MLP for a, where ReLU is used. Implemented in Py Torch, the models are trained over 1500 epochs with the Adam optimizer with initial learning rate 0.005 and using the ExponentialLft scheduler with decay factor 0.99 until the learning rate reaches 0.0001. We set At = 0.1, and use 5000 trajectories of length 10 At = 1 as training data, and 32 trajectories of length 60 At = 6 as testing data.
is computed as the normalized output of
MLP, with 4 input units, 16 units in each hidden layer
htw.storch.. www.rrchh.piii...... acc.duleccerrr.... Reducal.RnnPlaooo
5 "attps://boca.sccig htw.://gww...eaae///c/////ttenne htw.aateeennnengatteeeeeeeeeeeee gegggaaateaaaapppp wttpatte//cccyy.tttrtttee.e wtw.gnato.oooooooooooooo. wwwva_ivp..taaeee "nttpp::/pytorcch.. http:///pwooorc.ooo/////ssaaaaa/ htmi wwwrrcchh.ppin..... ar_achadulerr Exponentiall.
12

---

Published as a conference paper at ICLR 2020
B THE MAXIMUM LIKELIHOOD ESTIMATION PERSPECTIVE
In the presence of noise, 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:
(0) Integrator &o 20, fo. (=))
((=;0))
Ix 1=-18(011//2222)
/2x)/o22
(10)
P([ {[ 1=) II ((==;)) (V2To² nd
II =
111.-8,(0) (,00) 12/(22*)
and [(0 P((=)) is the likelihood function given the time-series data { [=) Note that. this model assumes independence between z; and 2j for i i once fixed.
If are to perform maximum likelihood inference, we arrive at the following:
max log [(01<=]]=)) = nd
log log(2mo²) 2
21 (0))
(11)
20
which is equivalent to
min Σ || - (0) 1=1
(12)
This provides motivation for using the I.2 loss, as we did the experiments.
So far, consider therefore the only argument of the likelihood function, while ån is fixed to be the observed initial state 2g. As generalization, we can consider strictly larger family of models by allowing 20 to vary as well. In this way, we treat both 0 and as the parameters in the model and therefore arguments of the likelihood function that we optimize on. In other words, the model becomes
the only parameter of the model defined by equations 10, and
2,0 20) Integrator 80. for (t.
a1z;0, 0, 20)
1
2, (0,50) (1)
V (2m))dodd
(13)
P([=))20 10 II (z;;0) (V nd II - 1| 8.000//20)))2200
and the optimization problem becomes
min 0,50
Σ 1 1 20, 1=1
(14)
which justifies optimizing over the initial states po, qo in addition to the neural network parameters 8 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 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 called the initial value approach. However, despite the success demonstrale in section 4.2, two difficulties of this approach have been pointed out: 1) The optimization could converge 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 suffer from bag local minima, what about optimizing the initial value &o?
C SYMPLECTICNESS AND INTTAA--TTATEE.OOTTMAAAAO The success of optimizing on the initial state of the system in addition to the recurrent H-NET and O-NET models described in section 4.2 raises the following question: If we
CONVEXITY
13

---

Published as a conference paper at ICLR 2020
already have relatively well-trained H-NET or O-NET, is the optimization the initial values convex? We formalize the question below and provide a heuristic answer.
For simplicity, we restrict our attention autonomous ODEs, which means that the function f in da = f()) does not depend on t. Assuming existence and uniqueness of solutions, there di
exists a function o, that maps each initial state 2. to the state of the system after evolving from &n for time i, (zn)) This function is usually called the flow map. Flow maps have also been defined for numerical solutions of ODEs, by letting (E)) = Integrator(as J: {1, [ with to O. 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 (<()) to be the state of the system after letting the system evolve from initial state for time t, for suitable choices of t. For example, for O-NET, we have (20) = Integrator for (t
Suppose we impose an L2 loss on Φe())), (²
|10e(20) zell2. where 2 corresponds to the observed data at time t. The question is, is (2)2) (perhaps locally) convex function of 2, for what functions and numorical integrators? To understand convexity, we compute the
gradient and the Hessian as follow.
①
[((2)= 2((x(()) TT)T.
(15)
Oz
8 (x) = 2()x(z) 21) (3) + Fe(2),
(16)
022
where F((2) is the Jacobian matrix of the flow map, defined F((2);) à (dr (=)()) and Ge(z) is a third-order tensor contains the second order derivatives of the flow map, defined as Gz(z)ezk 0z:2, (6))))). We use .(3) to denote the dot product in the third dimension.
āe,
F((2)* (=()) is symmetric positive semidefinite for any matrix F.(:). Φe corresponds to either the exact flow map of Hamiltonian system or the flow of symplectic integrator, such as the leapfrog integrator, applied to Hamiltonian system, then F²(z) is a symplectic = 1, which further implies
matrix, implying that det(F)(=)) 1. Hence, (F..() F+(=)) that F((z)T. F2(2 is positive definite matrix. Therefore, non-rigorously, when (2) 2 ll small and so the first term on the right hand side of equation 16 is negligible compared
to the least cigenvalue of F((z)T . (2),,, the entire Hessian matrix az2 (2(2 is also positive definite, implying strong convexity of the optimization problem.
8a
If 6, is the exact flow map, then (x)x()) zell being small means that noise in the data small. If de is the flow map of a learned model, then it mcans that we have model close to the true underlying system in addition to not having too much noisc 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 Hamiltonian syslem, and trained to be close enough to the true underlying system, and the noise in the data is small enough, then the optimization problem the initial state is strongly convex.
14

---

Published as a conference paper at ICLR 2020
D ADDITIONAL PLOTS OF THE SPRING-CHAIN EXPERIMENTS
D.1 NOISELESS DATA (SECTION 4.1)
(HNN) Single-step E-E H-NET
Single-step E-L H-NTT
Single-step L-L H-NET
Observations
position
position
20
40
60
80
100
20
40
60
80
100
mass
mass
2-
2-
position
position
0.
80
mass #5
mass #6
position
position
20
40
80
100
20
40
80
100
2.
position
position
20
40
60
80
100
20
40
60
80
100
mass
mass
4.
2
2.
position
position
0.
0.
-2
-2
160 -4
20
40
80
100
20
40
80
time step
time
step
Figure 5: Extension of Figure 1 to 10 masses and 10th being in the center).
the chain (Ist being the closest to one end,
15

---

Published as a conference paper at ICLR 2020
D.2 NOISY DATA (SECTION 4.2)
Single-step E-E O-NET Noiseless dynamic Observations
4
M
2
position
0
20
40
80
100
time step
Figure 6: Single-step E-E O-NET
(HNN) Single-step H-NET
Noiseless dynamic
Observations
2
position
0.
-2
-4
20
40
60
80
100
time step
Figure 7: Single-step E-E H-NET
Single-step L-L O-NET Noiseless dynamic
4
Observations
2
position
-2
0
20
40
60
80
100
time step
Figure 8: Single-step L-L O-NET
16

---

Published as a conference paper at ICLR 2020
Single-step L-L H-NET
4
Noiseless dynamic Observations
M
position
-2-
20
40
60
80
100
time step
Figure 9: Single-step L-I. H-NET
Recurrent E-E O-NET
Noiseless dynamic
Observations
position
-2
-4
20
40
60
80
100
time step
Figure 10: Recurrent E-E O-NET
Recurrent E-E H-NET
Noiseless dynamic
MM
Observations
position
-2
-4
20
40
60
80
100
time step
Figure 11: Recurrent E-E H-NET
17

---

Published as a conference paper at ICLR 2020
Recurrent O-NET Noiseless dynamic Observations
position
-2
20
80
100
time step
Figure 12: Recurrent I.-. O-NET
(SRNN) Recurrent L-L H-NET
4
Noiseless dynamic
Observations
2
position
0
-2
-4
20
40
60
80
100
time step
Figure 13: (SRNN) Recurrent L-L H-NET
Recurrent E-E O-NET w/ ISO Noiseless dynamic Observations
4
M ww
2
position
0
-2
-4
0
20
40
60
80
100
time step
Figure 14: Recurrent E-E O-NET w/ ISO
18

---

Published as a conference paper at ICLR 2020
Recurrent E-E O-NET w/ IS0 Noiseless dynamic
Mmm
Observations
position
0
40
60
80
100
time step
Figure 15: Recurrent E-E H-NET w/ ISO
Recurrent L-L O-NET ISO
Noiseless dynamic
Observations
position
40
80
100
time step
Figure 16: Recurrent L-L O-NET w/ ISO
(SRNN-ISO) Recurrent L-L H-NET w/ Noiseless dynamic
Observations
position
0
-2
20
40
60
80
100
time step
Figure 17: (SRNN-ISO) Recurrent L-L H-NET w/ ISO
19

---

Published as a conference paper at ICLR 2020
Vanilla RNN Noiseless dynamic Observations
4
M
position
0
-2
-4
20
40
60
80
100
time step
Figure 18: Vanilla RNN
LSTM
4
Noiseless dynamic
Observations
position
0
-4
20
40
60
80
100
time step
Figure 19: LSTM
20

---

Published as a conference paper at ICLR 2020
E ADDITIONAL PLOTS OF THE THREE-BODY EXPERIMENTS
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.
H-NET
18
-3.
20
15
15
O-NET
4:
Euler
Euler-Leapfrop
Leapfrog-Leapfrog
Figure 20: Actual versus predicted trajectories of the three-body system by the various suggeestep--rraned methods.
20
H-NET
O-NET
44
-1
30
28
13
Euler-Euler
Euler-Leapfrog
Ceapfrog-Teaplree
Figure 21: Actual versus predicted trajectories of the three-body system by the various recurrently trained methods.
-1
Figure 22: Actual trajectory versus the trajectory simulated by the leapfrog integrator with time-step (left) and 0.1 (right).
21

---

Published as a conference paper at ICLR 2020
F
ADDITIONAL PLOTS OF THE HEAVY BILLIARD EXPERIMENT
10-
10
00
00
-05
-10
-10
Actual Preditted
Attal Predicted
05-
05
-05
-05
-10
-10
-05
0.0
0.5
-1.0
-05
0.
10
Figure 23: SRNN with a rebound module that does not learn a (and effectively treats 1).
Actual Preditted
Attal hreditted
10
05-
05-
-05
-0.
-1.0
-0.
0.
0.5
100
-1.9
-05
0.
1.0
Aduall
1.0
00
0.-
-05
-10
-3.0
Figure 24: SRNN without the rebound module.
22

---

Published as a conference paper at ICLR 2020
Adual
10
00
0.-
-10
-10
Actual Predicted
Atual Tredeted
05-
05
-05
-05
-10
0.0
0.5
-1.
-05
0.0
1.0
Figure 25: Recurrent L-L. O-NET with the rebound module.
Adual
10
00
0.-
-10
-1.0
Actual
Atal
Preditted
hedicted
05
05
-05
-05
-10
-10
-05
08
05
-1.0
-0.5
80
10
Figure 26: Vanilla RNN.
23


Reviews:
[{"review_id": "Hyg1X9VTqS", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\n\nThe paper improves upon Hamiltonian Neural Networks to model physical systems from observed trajectories. Specifically, the authors propose to use (i) better integrator, (ii) multi-step learning, and (iii) initial state optimization. Authors experimentally show that all these improvements are beneficial in modelling complex and noisy Hamiltonian systems.\n\n\nMy Comments:\n\nI come from RNNs background and I am not an expert on physical systems. But the paper is extremely well written that I can easily follow. Experiments are convincing.\n\nI do not have any issues with the claims.\n\nAre the authors willing to release their code to reproduce the results?\n\n========================\n\nPost rebuttal: I stand by my decision.\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "H1lAbHqTYH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper introduces SRNN to model the hamiltonian of a dynamical system. Authors break down the design choices made in the algorithm and validate each one through experiments. \n\nI like the motivation of the paper of solving general physics problems using neural networks. The paper is well-written and the ideas are communicated properly.\n\nHowever, I have some concerns regarding the experimental results:\n\nIn Figure 1, even though the fifth mass seems to follow the exact pattern of observations for single-step L-L H-NET, the L2 error is increasing after each time-step suggesting that the fifth mass might not be the best one to consider for comparison between different algorithms here. Would be nice if a comparison between all trajectories was presented and discussed (perhaps in the Appendix)\n\nThe idea tested in section 4.2 seems not novel since previous works have already used recurrence to mitigate problems faced using one-step training. (https://arxiv.org/abs/1902.09689) \n\nI like the optimization method proposed in Section 4.3 and the results in Table 1 and 2 seems to justify its effect. Does this approach also handle degenerate cases in which different initial states can lead to the same trajectories after some amount of time? To me, it seems like a boundary on the noise variance should be assumed. Otherwise, no amount of optimization would be actually able to retrieve the correct initial state. Is this true or am I missing something?\n\n\n\n-minor comments:\n\n1. The conclusion section is missing from the paper. It would be nice to recap your findings there.\n2.  In section 2, the \"universality property\" of Niu's recurrent model should be explained or referenced. \n3. In Figure 1, the second \"Left:\" should be changed to \"Right:\".\n4. In Table 2, the Error std. for O-NET (E-L) should be bold not H-NET (L-L) unless it was intended to highlight the values for the model with the lowest mean error. (please clarify this in the paper.)\n5. Section 6: \"We focus on this section\" should be changed to \"We focus in this section\".\n\n \nOverall I think the work is interesting but it lacks some justifications regarding the claims made (as mentioned above), and although it is generalizable to other tasks and systems, it does not have sufficient novelty in its algorithm and approach.\n\nAs of now, I am recommending a rejection, but I am willing to reconsider my score should the authors address the above concerns.\n\n\n----------------------------------------------------------------------------------------------------------------------------\nUpdates:\n\nThank you for your response and addressing my concerns in the revised version of the paper. I also see new updates to the text which has improved the readability considerably. Thank you for your work.\nI have updated my score accordingly. \n\n\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "BJlTuVBRFr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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. \n\nOverall, the work is solid contribution and a reasonable improvement over the recent work on HNNs is demonstrated. Therefore I recommend acceptance of the paper. However, I have some fundamental doubts on the motivation on this line of works. This might be, because I'm not too familiar with the subject, and I'm willing to increase my score if the doubts are cleared.\n\nIn the shown examples, to best of my knowledge, the \"exact\" Hamiltonians describing the physics of the system are well-known. Therefore, I'm unsure what is the advantage of trying to learn physical laws, that are already well understood. The paper argues that the learned Hamiltonian will correct for errors in the discretization, but one could instead use a better integration scheme or a finer time-discretization, based on classical theory which has been developed over the last 50 years which comes with strong convergence guarantees and error bounds. I would have liked to see a stronger motivation, why it is interesting to learn an Hamiltonian of a system, where the exact Hamiltonian is already known. It would also be enlightening to see some plots, which illustrate how \"far\" the learned Hamiltonian is from the analytical one.\n\nOf course, one might argue that the ultimate goal is to have a learning based approach discover physical laws so far unknown to humans,  just from observations. But it is unclear why the inductive bias that the observations are generated by a Hamiltonian might be reasonable. It could very well be, that the law cannot be described by a Hamiltonian system. \n\nFrom a high-level point of view, one might even argue that it is not too surprising that one can fit a parametrized Hamiltonian to observations generated by a Hamiltonian system better than a general purpose function approximator without such an inductive bias or better than a system based on a naive/unsuitable non-conservative integrator.\n\nAs a remark, often the exact Hamiltonian is known to be (strictly) convex. I'm wondering whether convex function approximators such as convex neural networks could provide an even stronger inductive bias. But it might be that a general purpose RNN can account better for the discretization errors. ", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I made a quick assessment of this paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a novel architecture for learning Hamiltonian dynamics from data. The model outperforms the existing state of the art Hamiltonian Neural Networks on challenging physical datasets. It also goes further by proposing a way to deal with observation noise and a way to model stiff dynamical systems, like bouncing balls. The paper is well written, the model works well and the experimental evaluation is solid. All reviewers agree that this is an excellent contribution to the field, hence I am happy to recommend acceptance as an oral.
Decision: Accept (Spotlight)