Learning Energy Conserving Dynamics Efficiently with Hamiltonian Gaussian Processes
Authors: Magnus Ross
Markus Heinonen
Venue: Accepted by TMLR ()

Abstract: Hamiltonian mechanics is one of the cornerstones of the natural sciences. Recently there has been significant interest in learning Hamiltonian systems in a free-form way directly from trajectory data. Previous methods have tackled the problem of learning from many short, low-noise trajectories, but learning from a small number of long, noisy trajectories, whilst accounting for model uncertainty has not been addressed. In this work, we present a Gaussian process model for Hamiltonian systems with efficient decoupled parameterisation, and introduce an energy-conserving shooting method that allows robust inference from both short and long trajectories. We demonstrate the method's success in learning Hamiltonian systems in various data settings.

Published in Transactions on Machine Learning Research (03/2023)
Learning Energy Conserving Dynamics Efficiently with Hamiltonian Gaussian Processes
Magnus Ross University of Manchester
magnus. mmgpuu...msfpptotouuummmmmeeeeee
Markus Heinonen Aallo University
markus. o. h.s.ooonnnnaalloofffff
Reviewed on OpenReview: https: // openreview. net/ forum? id= DHEZukStat
net/forum?
Abstract
Hamiltonian mechanics is one of the cornerstones of the natural sciences. Recently there has been significant interest in lcarning Hamiltonian systems in a free-form way directly from trajectory data. Previous methods have tackled the problem of learning from many short, low-noise trajectories, but learning from a small number of long, noisy trajectories, whilst accounting for model uncertainty has not been addressed. In this work, we present Gaussian process model for Hamiltonian systems with efficient decoupled paramoterisation, and introduce energy-oonserving shooting method that allows robust inference from both short and long trajectories. We demonstrate the method's success in learning Hamiltonian systems in various data settings.
1 Introduction
Hamiltonian mechanics represent of the most important classos of dynamical systems, describing wide variety of natural phenomena from electromagnetism to the motion of planets (Salmon 1988 Taylor 2005). In this work we consider time-invariant systems characterised by Hamiltonian H p) R over position q(t) E R and momenta (t)) e R& over time & E R,, which system's 's configuration. The rules of evolution of a Hamiltonian system are defined by Hamilton's equations
can be thought of as the total energy of the
8H Öp
0H aq
p
(1)
We consider problem of learning free-form Hamiltonian H(.) entirely from observed system trajectories. The conventional mechanistic approach involves manually deriving the Hamiltonian H. and its evolution equations for a system of interest, and possibly estimating system cocfficients from data (Scinfeld 1970 Hernandez & Poznyak b020 However, for many systems the Hamiltomian is either unknown or too complex to derive from first principles (Schmidt Lipson 2009 Battaglia et al. 2018). Recently, numerous data-driven approaches have been introduced learn Hamiltonian systems with neural networks (Greydanus al. 2019 Cranmer et al. 2020: Zhong et al. 2019 Finzi et al. 2020). These methods return point solutions, and not able to characterise the uncertainty of the solution, which is important when data are limited.
Bayesian approaches, such as Gaussian processes (GPs), can be used to place a prior distribution over the derivative function of a dynamical system, allowing the posterior distributtion over the system dynamics to be computed in light of observations (Ridderbusch et al., 2021 Hegde et al. 2022). In this paper we place a GP prior over the Hamiltonian and infer its posterior directly from noisy trajectory data, shown in Figure I
Work partially completed whilst visiting Aalto University.

---

Published in Transactions on Machine Learning Research (03/2023)
True H(x)
Pred. H(x)
Pred. f(x)
Pred. x(t)
Inducing points x(t)
x(f) samples
Figure 1: The proposed model. We place a GP prior the Hamiltonian H(x) and, using a set of inducing points which lie on H. map function samples through Hamilton's 's equations to obtain system derivative f(x) samples, to which we apply ODE solver to obtain sample trajectories x(t). Shading represents model uncertainty.
We develop energy conserving variational multiple shooting scheme, which allows for efficient inference over long trajectorics, which usually present challenge for dynamical models due to the problem vanishing or exploding gradients (Ribeiro et al. 2020 Metz el al. 2021 p. Recently, two studies have introduced GP models that also aim to learn the Hamiltonian directly from trajectory data: the symplectic spectrum GP (SSGP) learns the system in a Fourier domain (Tanaka et al. 2022). while the structure preserving GP (SPGP) embeds the model within a numerical symplectic integrator (Ensinger et al.] 2022). Both models attempt to skirt the problem of trajectory length by using heuristic methods based learning from short. sub-sequences of the full data, which find works poorly for systems with complex behaviour.
In this work we present number of contributions, to both the inference methodology. and the experimental evaluation, of GP models for Hamiltonian systems. Together they can be summarised as follows:
We propose Hamiltonian Gaussian process parameterised by inducing variables with efficient functional sampling.
We adapt the variational multiple shooting method of Hogde et al. (2022) energy conserving Hamillonian systems, for highly accelerated optimisation and increased performance.
We provide an extensive experimental evaluation, and find that our method shows strong performance in a number of settings, whilst discussing the areas in which GP based methods are limited.
2 Primer on Hamiltonian mechanics
We consider a Hamiltonian dynamical system over the 2D-dimensional phase space of canonical positions q R (for example, coordinates angles) and their associated canonical momenta p R The system is characterised by the Hamiltonian energy function H(q. p, t) E R (Thornton & Marion 2004). In this work we will restrict the discussion to the common case of time-invariant Hamiltonian, i.e. Hqq. = H(q. p), in which case, the Hamiltonian energy is conserved over the system trajectories. The temporal evolution of the system is described by set of coupled first order differential equations, known as Hamilton's equations,
dq
0H 0p
dp 0H
q
p
di
dl
Qq
The structure of Hamilton's equations ensures that various constants of motion are conserved over the system trajectorics. Chief among these constants of motion is the energy of the system, but they can also include linear momentum, angular momentum, and other quantities. The energy conservation is shown by
A dHl 0H.
07 0H 0H
0H 0H
+
p
0.
(3)
dt oq
0p
oq 0p
Öp 0q
2

---

Published in Transactions on Machine Learning Research (03/2023)
By using Hamilton's equations we obtain system of 2D first order differential equations,
dx dt = f(x) =
OH
q E ER
Öp
(4)
where we have concatenated q and p into the 2D dimensional phase space state vector x, and f R2D RR22 is the system time differential. We obtain system trajectories by forward integration
x(t) = x(t; xo) = X0 S. f(x(T))dr,
(5)
where x(t) the state of the system at time t, Xy is the initial state of the system, and T E [0, is an integration time variable. Typically for mechanical systems the Hamiltonian is determined by identifying the kinetic and potential energies of the constituent parts, which form the total energy H. with the process becoming increasingly diffieult for more complex systems. In the this work, aim to forgo this process and to learn the Hamiltonian directly from trajectory data, with no assumptions on its functional form.
3 Hamiltonian Gaussian processes
In order to infer the Hamiltonian of the system, we assume it follows a Gaussian process (GP) prior (For revicw, see Williams & Rasmussen 2006b)
H(x) ~ ((x))~P((O,A) (xx,xx)
(6)
which models distribution over energy surfaces with zero mean E[H((x) 0 and kernelised covariance,
cov[H((x), H(x')] (x(xx x'), kH R2D
R 2DD + R
In a GP the probability of a function at any finite subset of evaluations follows a multivariate Gaussian (H(x))...., H((x)) H((x)) ~ AO0,, Kx ()),
where KH (X, X) E NNN is positive definite kernel matrix with elements [>x X)];; = lxy (x,,,)) In order to perform inference whilst using GP representation of the Hamiltonian, we adapt the method described by Hegde et al. 2022), which allows for inference of ODE systems with GP based derivative function. This is possible because GP prior over the Hamiltonian implies GP prior over the derivative function, as we will see in the following section.
3.1 The time derivative
Given the prior over H it is necessary to define the system time derivative x, in order to compute trajectories for inference and sampling. Equation {4] can be rewritten (Rath et al. 2021) as
0
0 -0,)
x = f(xxx) = CH(x), with L =
Vx
Poisson tensor
Gaussian processes are closed under linear operators Williams & Rasmussen 2006
Agrell
2019
and hence
the Hamiltonian and its vector field follow a zero-mean joint GP
(H(xx) ~ GP (".(k (x)(xx') x!) F(x)
kHr(x, x') Kr(x,x')
(10)
(xx(xx)))
with covariances induced by the Hamilton's equation,
cov [7((x), f(x') = kxr((xxx)) =
Op
ki kx x')
222x1
(11)
-
02
-02
cov [f(x), = kf(x,x))= =
Kx (x(xx))
E 2D222
(12)
82
82
qp
3

---

Published in Transactions on Machine Learning Research (03/2023)
3.2
Inducing points
In order to allow tractable inference and sampling, we introduce a set of inducing points (Snelson & Ghahra- mani 2006) to obtain finite, parametric representation of the infinite joint GP (H,f). We condition the Hamil- tonian with M inducing energies u = (M1,..., uM) E R at phase space locations Z = (Zi...., ZM) E R xxD that encode energy pseuuooobservations u H(z)). The vector field conditioned on these pseudo observaiossrvtsonns is again Gaussian process,
(H()) = =u) gP( kHr(x, Z)) Z) u, Kr(x,x') - (Pp(x, 2), (104))('pp) KHH (e, x)):
(13)
mean
covariance
where the set inputs refer to expanding the corresponding function over them. By varying the inducing parameters (u, Z), we can represent approximately arbitrary Hamiltonians and their unique vector fields, which can then be forward integrated to obtain simulated trajectory solutions. In later sections these parameters are the main variables be learnt. Our approach differs from Hegde al. (2022) by placing the inducing points on the Hamillonian, not the derivative function. We then require set of scalar inducing points, instead of vector-valued 2D-dimensional inducing points.
3.3 Sampling
In order to simulate trajectories x(0). x(T) we must be able to sample a persistent derivative function or vector field = f GP from the GP (3), and evaluate it along the trajectory x(t). Standard methods for sampling functions from GP are based kernel matrix decompositions, which have prohibitive complexity of O(N5) for N evaluation points. We follow Hegde et al. (2022), and bypass this problem by 'decoupled'
paramcterisation [Wilson et al. 2020). sampling from the Hamitonian instead of the derivative function,
H Hw.u.z(xx) = Σ wx +
Σ (x,
(14)
where wi NNO, 1) are random weights of the S Fourier basis functions dx(x) cos(a, , x frequencies a sampled from the spectral density of k, and Bg (0, 2m) (Rahimi & Recht 2007 and v = k2 (u fw) where +- <(Z) E M represents the evaluation of bases over inducing inputs. We denote all Fourier-related parameters by w = {w, a. B). By sampling and fixing w, we obtain a deterministic function sample that can be evaluated anywhere. Finally, we transform the energy surface H. to a derivative vector field via fx,,,z(x)= = Cxw,u.z(x), which allows trajectory rollout using numerical integration.
B,) with
3.4
The probabilistic model
We aim to infer H. from one or more noisy rcalisations of trajectorios from the true system, with arbitrary length and time irregularity, each of which have different unknown initial conditions. We present simplified notation for a single observed trajectory (See Appendix A.3 for multiple trajectory derivations).
Let the trajectory obscrvation be denoted Y (y1....yN) y Nx20 where y; (t.) + E RD is the i'th noisy state of the system al time 4 E (41,..., (N),, and N is the number of time steps observed. We begin by assuming Gaussian prior on the inducing energies and the initial state Xo,
p(u) N (u/o, k(Z,
(15)
p(xp) (x, I)
(16)
where k(Z, Z) E MRM is covariance matrix for the inducing stales, with [k(Z, Z) = (H 2z). We follow the convention in sparse Gaussian processes to treat the inducing locations as hyperparameters that are merely optimised, instcad of inferring their posterior distribution (Hensman et al. 2015a) (For
in-depth discussion, see Hensman et al. (2015b); Rossi et al. (0221) We omit all hyperparameters from the probabilistie model notation below for simplicity.
4

---

Published in Transactions on Machine Learning Research (03/2023)
Z
YN
W
H
X0
x(4)..
x(LN)
f
Hamiltonian model
Plate diagram over the model (blue) and ODE system (red). We denote observations {y:} with shaded nodes, random variables (w, u, xp) as white nodes, hyperparameters (Z) with rectangles, and functions (H, x) without shapes. Both sides contain the same
Figure 2:
The joint distribution over data Y, sampling parameters w, energy variables u and initial state xo is
p(Y, w, u, xo) = pl Y u, (.))))())((()ppxx
(17)
N
I
u. u.)
(18)
where the likelihood is
ply: W. u. xo) = -N)xx|x((())
x(() = Xw,uxx (1)
(19)
where first compute the Hamiltonian H. with Equation
14 then its derivative f using Equation
and
finally integrate forward Xo - x(() with Equation 5 The form of the joint distribution here mirrors Hegde
et al. 2022p. Our goal is to infer the intractable posterior
p(u, Xo Y) = p(u, Xo Y, W)p(W)dw.
(20)
over initial state Xg and inducing variables u, where assume the tractable sampling variables to be marginalised. Figure 2 shows the dependency structure between the variables in the model in graphical form. We turn to variational inference to approximate the posterior.
3.5
Variational inference
We follow Hegde et al. (2022) and the framework of stochastic VI (Hofman al. 2013 Hensman et al. 2013 to infer the posterior of Equation 20 We assume factorised postcrior approximation
(u, xo) (x)))(Xo)
(21)
(u) NNulmm Q)
(22)
q (Xx)) A (ml Qo).
(23)
with variational free parameters 0 (m, Q, m.. Qu). Variational inference seeks to find
arg ming KT. (u, ( and the true posterior, which is equivalent to maximising the evidence lower bound Blei et al. 2017
x Y)], the Kullback-Leibler divergence between the approximate posterior
F(0,2 Z) [g(u)g(xo)p(W))
Σ log pys/ xo) KL [a(u)][p(u)) - KL [a(x))]lp(xx)]
(24)
5

---

Published in Transactions on Machine Learning Research (03/2023)
The expectations integrate the inducing energies u, initial state Xo estimates, and Fourier function sampling determined by w. The likelihood requires solving the trajectory Xw.u,xo (1) numerically. The expectations are Monte Carlo averaged, while the KL terms have closed-form solutions. The bound is maximised by gradient ascent wrt and Z. See appendix A.1 for details.
4 Energy-conserving shooting parallelisation
Optimisation of the variational bound 24 for long sequences is challenging in practice due to the problem of unstable gradients (Ribeiro et al., 2020 Metz et al. 2021). For long trajectories, small perturbations at early times can compound into large effects at later times, resulting in gradients of the bound vanishing or exploding (Haber & Ruthotto 2017). Although this problem has primarily been discussed in the context of neural networks (Kim et al. 2021 Choromanski et al., 2020), Hegde et al. (2022) introduces probabilistic shooting
solution for GP-ODEs. Shooting is an ODE optimisation technique where the 'long' solution x(0) + x(T) is split into consecutive segments xx(t), (((+1)) that are solved in parallel, while ensuring that the neighboring
segments match (Hemker
1974 Bock & Plitt 1984h (For review, see Diehl & Gros (2000
4.1 Energy-conserving shooting model
We adapt the shooting formulation of Hegde et al. (2022) to Hamiltonian GP ODEs, in order to stabilise gradients and allow learning over longer soquences. We begin by augmenting the system with a set of L < shooting variables = {= E IR4x2D, which represent the state of the system at times lt E (41 and split the continuous state solution x(); into L distinct segments from initial states ss with segment solutions
x(t: st) = / f(x(T)))T,, for t [ti, t(1).
(25)
That is, to solve particular time f, need to find its interval [t,, 4+1] and solve from st. The splitting leads to more stable gradients since each segment has less non-linear solution map (Diehl & Gros 2020 Furthermore, the segments can be solved in parallel. The key problem of shooting is to ensure continuity x(1;; st-1) = St at segment boundaries 4,, otherwise the combined solution will be discontinuous.
We formulate this with an error model
x(tran))
x(t((1)a(2)
x ((((()):(0)
S x(f(; E f o(
(26)
where & 2R22D is the between-segment tolerance of posi- tion and momenta. To conserve energy, also introduce energy Volerances to encode the permissible energy change between segments,
So
81
82
Shooting ODE
Figure The shooting system.
H(s:)) H(x(tr:s)))) + +X+ XI ~(0,0)
(27)
Together these constraints lead us to product prior
p(s||st-1, w. u) = (((ss)((tssss))) (H)
H()) = H..u(-)
(28)
where H. is the Hamiltonian function conditioned by the w, u, and x(11;81-1) depends CH via Equations
(4h and (5) We place Gaussian prior on the initial state, p(so) = 00 I). Let Sig) denote the last model is
shooting variable before time t, The joint distribution of the
P(Y W. S)
I w,y u,
I p(s///-1,, s1-1, w, w) p(W)p(u)p(so).
(29)
shooting likellhood
tolerance prior

---

Published in Transactions on Machine Learning Research (03/2023)
4.2 Variational shooting inference
The inference of the posterior S|Y) is again intractable. We define a variational posterior approximation g(u, q(u) II ((si) with independent Gaussians q(s:) (sss Ex) on the shooting states. This results in new evidence lower bound
N Be(e [log(yy/ww u, s()) + Σ (() [logp(s) 81-1, w, u)]
F(0,2) Z) (v(w)((w))
(30)
L
Σ H[[gs:)] - (L lI p(s)) - KL(g(u) || p(u)],
where 0 = (m, Q, {ar, EtH). The first expectation is the 'short' likelihood of each observation solved only from previons shooting state. The second term ensures the quality of the match between the shooting segments. The remaining terms regularise the model. This bound is similar to that of Hegde et al. (2022), but places the inducing points on the Hamiltonian, not. the derivative function, and includes an additional energy based matching term for the shooting states. See appendix A.2 for details.
5
Related work
Mothodologically our work is founded on Hegde al. (2022), who present learning of non-Hamiltonian ODEs with GPs using variational multiple shooting for ellicient long trajectory inference. We extend with Hamiltonian structure on both the system and the shooting method to control energy conservation.
Hamiltonian Gaussian processes. The most similar existing work to ours are the symplectic spectrum GP (SSGP) model of Tanaka et al. (2022), , published concurrently, and the structure-peesrringg GP (SPGP) model [Ensinger et al. [2022), both of which place a GP prior over the Hamiltonian function H. The SSGP model closely related lo ours, but dillers inference. They use staudard RFTs for sampling, which suffers from variance starvation and poor uncertainty representation (Wilson et al. 2020 and place a variational distribution on the RFF weights instead of inducing energies. The SPGP proposes embed a Hamiltonian GP specifically within symplectic integrator to ensure numerical volume preservation. Neither model supports initial state estimation and neither use shooting approximations, but instead resort to heurist.c minibatching of subsequences of the trajectories.
Hamiltonian maps with Gaussian processes. A number of works aim to learn the Hamiltonian flow map x(finit) + x(tānat); or related mappings, from initial and final conditions. These methods typically require large number of low-noise trajectories (10s or 100s). Rath et al. (2021) use GPs to learn the flow map of Hamiltonian systems, providing implicit learning scheme for non-separable systems, and efficient explicit scheme for separable systems. Offen & Ober-Blöbaum (0222) develop scheme known as shadow symplectic integration (SSI), which allows for compensation for the error incurred by forward numerical integration, leading to more accurate preservation of symplectic structure. Bertalan et al. (2019) learn the Hamiltonian directly in phase space, which requires an approximation of the time derivatives of the coordinates using first differences, introducing additional error, especially for noisy data.
Hamiltonian neural networks. There has been a significant amount of interest in learning Hamiltonian dynamics with neural networks (NNs). Greydanus al. (2019) introduce the Hamiltonian NN (HNN), which trains an NN to predict the dynamics with auxiliary loss based on Hamilton's equations. The original formulation of the HNN requires computation of the time derivatives of the coordinates, and cannot learn directly from trajectory data. Zhong et al. (2019) introduce version of the HNN which computes trajectory rollouts using an ODE solver and can learn directly from trajectories, in addition to incorporating control signals. Finzi et al. (2020) introduce constrained HNNs for mechanical systems, in which the structure of the system is encoded as constraints, and the Hamiltonian is learned Cartesian coordinates. Gruver et al. (2021) investigate the effect of the inductive biases in Hamiltonian and Lagrangian NNs, and find that for realistic systems, baseline models without energy conservation often perform better.
7

---

Published in Transactions on Machine Learning Research (03/2023)
6 Experiments
In this section we provide experimental evaluation of our method for variety of Hamiltonian systems. The code for our implementation of the model is available at httts:///iithut http://wwww....////guuuuu///////
6.1
Experimental Setup
Tasks. Aside from an initial illustrative example, we study two distinct. tasks:
Task 1: Trajectory forecasting. In this task, the model is provided with a single noisy trajectory of length T with unknown initial condition. The model learns the system from [0, 2J] interval, and is lasked to forecast the trajectory forward for [T,2T]. 2T].
Task 2: Initial condition extrapolation. In this task, the model is provided with K noisy trajectories of fixed length from different, unknown initial conditions. The model is evaluated on its ability to forecast trajectories from new set of noise-free initial conditions.
Most prior work learning Hamiltonians deal with variations of task 2, for example Tanaka al. (0022;
; Rath
et al. (022)). Experiments for task 1 are provided by Ensinger al. (2022) and Hegde et. al. (2022) among others, but prior work has provided a joint evaluation of both tasks. For both tasks report the trajectory and Hamiltonian energy root mean square error (RMSE), and the mean negative log likelihood (MNLL). We repeat each experiment 10 times with different initial conditions, and report median and interquartile range on all tables and plots due to the high variability the results for all models. We give tables of means and standard errors in appendix C
Systems. We evaluate our model on three true Hamiltonian systems: the fixed pendulum (FP), the spring pendulum (SP) (Lynch] (200),, and the Henon-Heiles (HH) system (Henon & Heiles) 1964). The FP is two-dimensional with predictable periodic dynamics, and serves as simple test case. The SP and HH are both four-dimensional, and exhibit more complex, chaotic dynamics. The system definitions are
(fixed pendulum)
HFP(q.P) p) = mgr(1 - cos q)
2mr
Har P1, P2) = 1
1
(spring pendulum)
pi + (g) + r)²
+
kqp - mgr cos Q2
2m
2
[Henon-Heiles)
(2.P1, Pz) = (9 + 2 +
r4 + p2)
9229
_02
where m is the pendulum mass, r is its resting length, is the gravitational field strength, k is the spring constant, and is a parameter controlling the magnitude of the HH potential. For more details on each system, see the appendix B.1 We add Gaussian noise to the training data for each task, with variance set to 5% of signal variance, unless otherwise mentioned. For cach experiment we sample different initial conditions from the phase space of the system generate data for each repeat, using the same data across each model test.
Our HGP setup. We M 48 inducing points for task and M 128 for task 2. Throughout use S 256 basis functions, and run optimisation for 2500 iterations using the Adam optimiser (Kingma & Ba 2015) with lcarning rate 3 During training we use single sample from the model to estimate the intractable expectations, when making predictions we use 32 samples. Unless otherwise stated, we use the shooting approximation with single shooting state per 4 data points, iee. L [N/4]. We use the torchdiffeg package (Chen 2018) with implicit dopri5 solver. For more details on the setup of the HGP and all the baseline models see appendix B.3
Hamiltonian aware initialisation. We optimise the bounds in Equations 24 and 30 using gradient ascent. the variational parameters, and the set of model hyperparameters. This optimisation problem is
8

---

Published in Transactions on Machine Learning Research (03/2023)
True vector field
Inferred vector field, } cycle
Posterior samples, 1 sycl
Energy MSE, A cycle
15-
0
0
1 cycle
Posterior samples, cydle
Energy cycle
Training data
True trjeetory HGP
GPODE
2.-
1.5
0.0
Figure The HGP can accurately model energy-preserving Hamiltonians. Top row: the true pendulum system (FP) followed by learnt vector field of HGP (orange) and non- Non-Hamillonian GPODE baseline (green) from half a cycle of data, followed by sample trajectories and the error predicted energy along the trajectory. Bottom row: full cycle estimates, where the baseline is able to fit the vector field, but is oblivious to the energy.
challenging, such good initialisations very important, particularly for the variational parameters of the Hamiltonian GP. We can obtain a good initialisation by observing that the derivative function f and H are jointly a GP. Using the approximate numerical derivatives of the trajectory data, we can form some data for f, which can then be conditioned on, using equations (10) to (2), to obtain an estimate for the mean function of the Hamiltonian. We use this estimate as the initial value for the variational mean m. See appendix B.2 for details.
NN baselines. We compare against two NN models: a standard (non-Hamiltonial) Neural ODE (NODE) CChen al. 2018)) with a neural time derivative (x),, and HNN, based on the Unstructured SymODEN model of Zhong et al. (2019) with neural Hamiltonian Ho(x), Both methods perform poorly when trained raw trajectories due to vanishing/explodiing gradients. To allow for fair comparison, split the trajectories into smaller sub-trajectories by running sliding window of fixed length over the data, and form these into batches of training data.
GP baselines. In order to evaluate the effect of our Hamiltonian prior, compare against the non- Hamiltonian shooting GPODE baseline of Hegde et al. 2022).. We use the same hyperparameter settings and initialisations as in the HGP model, where applicable. To initialise the inducing variable distribution means, we follow the same procedure described by Hegde al. (2022). We are unable to provide results for the SPGP (Ensinger et al. 2022) as the authors did not release their code publicly, and would not make an implementation available on request. A public code exists for the SSGP, however it does not support systems with D 1. In order to provide some comparison re-implement one of the key points of difference between the models, the sub-sequence training method, and test its performance. We also provide a comparison with the SSGP for the toy FP system in Appendix C.3
6.2
Toy task
As an initial test of the effect of the Hamiltonian prior we train both HGP and GPODE, without shooting, the FP system with both half cycle and full cycle of data. Figure shows the inferred vector fields,

---

Published in Transactions on Machine Learning Research (03/2023)
the posterior trajectory samples, and the error in the energy of these trajectories. We see the HGP producing plausible trajectories and appropriate vector field estimate even far from the observations when given half a cycle of training data, whilst the GPODE does not, and produces samples with a large energy violation. With full cycle of data both models fit the data well, however the GPODE still has high energy violations. This loy lask illustrates the ability the HGP to learn accurate dynamics from a relatively small amount of data by its Hamiltonian inductive bias.
6.3 Task 1: trajectory forecasting
State RMSE (4)
State MNLL (L)
Energy RMSE (L)
Method
FP
НH
SP
FP
HH
SP
FP
HH
SP
NODE
0.12 (0:12)
0.42 (0.25)
1.13 (0.82)
0.30 (0.28)
0.02 (0.01)
1.02 (0.41)
INN
0.20 (0.16)
1.56 (0.26)
1.41 (0.30)
0.34 (0.33)
0.04 (0.10)
1.12 (0.25)
GPODE
0.23 (0.26)
0.54 (0.71)
0.91 (0.36)
-0.03 (1.21)
1.20 (4.19)
(2.77)
0.38 (0.34)
0.02 (0.02)
0.78 (0.50)
HGP
0.18 (0.19)
0.32 (0:30)
0.64 (0.50)
(0.67)
0.31 (1.89)
(2:20)
(0.42)
0.01 (0.01)
0.81 (0.36)
Table Performance comparison of different mothods on each system on the trajectory forecasting task.
Table compares the performance of HGP to baseline methods for each of system on the forecasting task. We use train trajectory lengths of 40, 8 and 16 seconds for the HH, FP and SP systems respectively! For the more complex systems (SP, HH) the HGP provides strong performance particularly when compared to the equivalent GP baseline, the GPODE, providing significantly better point predictions and uncertainty quantification. For the simple FP system the NODE provides a similar result to the HGP, for the SP and HH systems the HGP provides the best performance. In this and subsequent experiments, we found that the NODE unexpectedly outperformed the HNN, this issue is discussed in Section 7 Table [ additionally shows the energy RMSE, computed by applying the true Hamiltonian to the predicted trajectories, showing that the HGP is able to recover the energy conserving dynamics well. Figure 5 shows the cumulative errors accrued over the test period. Appendix C contains plots of the model predictions for sample trajectory from each system These animations are illuminating because they show that, in addition to providing strong predictive performance, the GP models produce trajectories that are more physically plausible than the NN models.
Spring
Henon-Heiles
NODE
200
GPODE
200
30
150
150
20-
10
50-
1.00
1.25
1.00
1.35
1.00
1.25
1.50
7 [rel))
T (rel.
T (rel.)
Figure 5: The HGP with energy shooting has lowest cumulative trajectory error for the complex systems. Cumullative error over test period (T,77) for each system in the forecasting task. The horizontal axis shows the time relative to the start the training period.
Effect of trajectory length. Figure Ga shows the effect of increasing trajectory length model performance for task 1 on the HH system. The point-wise predictive performance of each model stays approximately constant with increasing trajectory length, with performances degrading slightly for the longest trajectorics,
We aim give trajectory lengths that similar amount of phase space for cach system for fair comparison, with the different times reflecting the different properties of each system.
2 A nimations the trajectories be found at hattps:/magnusrosss github. setthab.io/tamiiitoiiinnni
10

---

Published in Transactions on Machine Learning Research (03/2023)
HKN
3.-
HNN HGGPES
HGP-ES
1.5
NODE GRODE
GPODE
2.0
MPALL
RIISE
30
18
10.0
25.0
1.0
13.0
25.0
(())
(()
o3 Rel (%)
(%)
(a) Trajectory length
(b) Relative noise variance
Figure 6: The HGP performance robust to observation length and noise. Analysis effect of data for the HH system on task 1.
due to the difficulty the optimisation problem. The HGP gives consistently good performance over each trajectory length relative the bascline models.
Effect of noise. Figure 6b shows the effect of incrcasing noise on performance, with the horizontal axis showing the noise variance as a percentage of total signal variance. The results show that, in terms of RMSE, both the NODE and HGP largely robust increasing noise, although the is some drop, with the GPODE degrading more severely.
Initialisation
RMSE (4)
MNLL (4)
Random
1.03 (0.06)
1.49 (0.09)
Hamiltonian
0.36 (0.19)
0.42 (0.83)
Table 2: Effect of inducing initialisation.
Effect of initialisation. Table illustrates the effect of the Hamiltonian aware initialisation scheme discussed in section 6.1 We ran HGP task 1 for the HH system, with randomly initialised inducing mcan, and our proposed initialisation. The models initialised randomly fail to fit the data well, and produce solutions that extrapolate poorly.
Comparing inference schemes. Figure provides comparison between different inference methods for task 1 on the HH system with varying trajectory lengths. We compare the HGP with: energy conserving shooting (HGP-ES), standard shooting with no additional energy constraint (HGP-S), no shooting (HGP), and with inference on sub-sequences formed into batches (HGP-Batched).. With the HGP-Batched model, we aim emulate the inference method used for the SSGP and SPGP, see appendix B.3 for details. The results show that the HGP-ES provides the best equivalent performance over each trajectory length, with
HGP-ES
1.0
150
12
HGP-5
Runtine (s/step)
0.8
HGP
HGP Hatched
1.00
0.6
0.50
0.2
0.00
18
30
42
54
30
42
54
18
30
42
54
(9)
(a) Model performance
(b) Runtime
Figure 7: Energy shooting inference performs best across trajectory lengths, whilst speeding up inference 15x for long trajectories. Comparison of inference methods for task on the HH system, 7a shows performance whist varying trajectory lengths, 7b compares the time taken to compute the evidence lower bound for the shooting and standard methods for varying trajectory lengths.
11

---

Published in Transactions on Machine Learning Research (03/2023)
performance relative to alternative methods improving with increased trajectory length. The HGP-Batched method performs poorly for all trajectory lengths, with performance degrading for longer trajectories. Figure 7b shows the difference in runtimes for a single computation of the evidence lower for different trajectory lengths for the shooting and non-shooting models. The shooting approximation provides a 15x specdup for the longest (54s) trajectory and 6.5x speedup for the shortest (18s) trajectory, whilst matching or improving on the standard HGP in terms of RMSE and MNLL.
6.4
Task 2: Initial condition extrapolation
State RMSE (L)
State MNLL (1)
Energy RMSE (4)
Method
FP
HH
SP
FP
HH
SP
HH
NODE
0.48 (0.15)
1.15 (0.19)
0.92 (0.38)
0.26 (0.06)
0.03 (0.00)
1.35 (0.31)
HNN
0.91 (0.47)
1.35 (0.18)
1.13 (0.26)
0.28 (0.15)
0.05 (0.03)
0.96 (1.79)
GPODE
2.41 (2.28)
3.58 (0.79)
3.53 (1.14)
7.11
(8:22)
19.11 (10.51)
11.48 (10.73)
6.17 (10.30)
1.52 (0.95)
1.9x10 (2.5 x 10")
HGP
(0.46)
1.30 (0.28)
(0.50)
4.68 (2.95)
7.60 (3.03)
6.13 (0.81)
0.55 (0.60)
(0.01)
3.32 (8.86)
Table 3: Performance comparison of different methods on the initial condition extrapolation task, with K trajectories.
Table 3 shows the results for each model each systems for task 2. For this task each model was given set of K = 8 noisy trajectories, of lengths 12, 4 and 6 seconds for the HH, FP and SP systems respectively, with initial conditions sampled randomly in phase space. The test set consists of 25 trajectories sampled from phase space using the same procodure, with length triple that of the training trajectorics. In terms of state RMSE, the best performing model is the NODE, which provides siguificantly better predictions than the HGP for the simplest system (FP) and marginally better predictions for the more complex systems. The HGP performs better than the GPODE on all metrics, again illustrating advantages of the energy conserving prior. The HNN and NODE perform best in terms of energy RMSE, with HGP performing slightly worse, and the GPODE performing poorly.
40
Number of trajectorics Figure shows the per- formance of each model on the HH system for task 2 as the number of training trajectories K is increased. We can see that whilst the HGP performs better than the GPODE on both metrics for all K, for large K the NN models provide better point predictions. We believe the poor relative performance of the GP based models in task 2 can be attributed to the choice of prior, which we discuss in section
IKPP ES
35
30
NODF
30-
GPODE
25
MHLL
2.0
1.5
10-
1.0
0.5
32
16
Figure 8: The HGP outperforms GPODE for all K, but NNs are best for large K.
7 Discussion and limitations
Performance of HNN NODE. We found that in both tasks 1 and 2, the HNN model performed significantly worse than the NODE, which is unexpected and differs from previous results in the literature. We offer two possible explanations for this, both of which likely play a role. Firstly we use data that is both noisier, and use longer, more sparsely sampled trajectories than previous studies. Additionally in order to provide a fair comparison across methods, we learn Hx)x), so that there no restrictions on the set of Hamiltonians. For mechanical systems, the Hamiltonian can be rewritien in terms of the potential and kinetic energies, and NN can be used to represent each term. Zhong et al. (2019) find that learning the mechanical Hamiltonian in this way leads to improved performance. We believe it should be possible to modify our framework to learn the potential and kinetic terms separately using HGP, and this is an avenue we would like to pursue in future work.
Suitability of the prior. The performance of the HGP on task 2, somewhat underwhelming especially for larger K, relative to NODE. We believe some of this effect be attributed to our choice of GP prior.
12

---

Published in Transactions on Machine Learning Research (03/2023)
We use stationary GP prior, which is likely non-optimal for most Hamiltonian systems, which are typically nonstationary. The assumption of non-stationarity likely leads to the poor generalisation al new phase space points, since the GP prior will revert to zero mean functions there, which is not the case for the NN models. One option to improve the suitability of the GP prior is via the use of specifically designed kernels to represent symmelries the system Ridderbusch et al. 2021). or kernel structure learning (Kim & Teh, 2018). Another is extend our framework and instead represent H with deep GP prior (Damianou & Lawrence, 2013). Deep GPs are adept modelling complex, nonstationary functions, and would be well suited to the task. Using the SVI scheme proposed by (Salimbeni & Deisenroth (2017) would make integration of deep GPs into our framework relatively straight forward, and would be interesting to undertake as part of future work.
Control. The present method has significant potential at improving Bayesian online (Diiserrthh & Rasmussen 2011) or policy-based RL (Yildiz et al. 2021) by incorporating Hamiltonian inductive biases. This requires expanding the model towards Hamiltonian systems with external forces, which relax the energy conservation assumption.
8 Conclusion
In this work presented Gaussian process model to learn Hamiltonian dynamical systems from trajectory observations. We proposed a parameterisation that combines inducing points and Fourier bases, and introduced novel energy-conserving shooting method to allow reliable inference from long data. Our experiments show strong and stable performance under various learning settings.
Acknowledgements
We would like to thank Pashupati Hegde for useful conversations and advice, and in particular for help with the code for the GPODE. We would also like to thank Mauricio Álvarez and Tom McDonald for fcodback on earlier drafts of this work. This work was supported by MR's research visit to Aalto University through the ELLIS network. The calculations were performed using resources within the Aalto University Science-IT project. This work has been supported by the Academy of Finland (grant 334600).
References
Christian Agrell. Gaussian processes with linear operator inequality constraints. Journal of Machine Learning Research, 20(135):1 36, 2019.
Peter Battaglia, Jessica Hamrick, Victor Bapst, Alvaro Sanchez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Caglar Gulcehre, Francis Song, Andy Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey Allen, Charles Nash, Victoria Jayne Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet Kohli, Matt Botvinick, Oriol Vinyalsand Yujia Li, and Razvan Pascanu. Relational inductive biases, deep learning, and graph networks. arXiv,, 2018.
Tom Bertalan, Felix Dietrich, Igor Mezić, and Joannis G Kevrekidis. On learning Hamiltonian systems from data. Chaos: An Intendisciplinary Journal of Nonlinear Science, 29(12), 2019.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe, Variational inference: review for statisticians, Journal of the American statistical Association, (2(1)))) 877, 2017.
Hans Georg Bock and Karl-Josef Plitt. A multiple shooting algorithm for direct solution of optimal control problems. IFAC Procendings Volumes, 17(2):1603 1608, 1984.
Ricky T. Q. Chen. torchdilfee,, 2018. URL
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In NewrIPS, 2018
Krzysztof M Choromanski, Jared Quincy Davis, Valerii Likhosherstov, Xingyou Song, Jean-Jacques Slotine, Jacob Varley, Honglak Lee, Adrian Weller, and Vikas Sindhwani. Ode to an ODE. In NewrIPS, 2020.
13

---

Published in Transactions on Machine Learning Research (03/2023)
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Batlaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXio, 2020.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In AISTATS, 2013.
Marc Deisenroth and Carl E Rasmussen. PILCO:: A model-based and data-efficient approach to policy search. In
JCML, 2011.
Moritz Diehl and Sébastien Gros. Numerical optimal control University of Freiburg, 2020.
John R Dormand and Peter J Prince. A lamily of embedded runge-kutla formulaae. Journal of compalationa and applied mathematics, 6(1):19-26, 1980.
Katharina Ensinger, Friedrich Solowjow, Sebastian Ziesche, Michael Tiemann, and Sebastian Trimpe. Structure- preserving Gaussian process dynamics. In ECML/PKDD, 2022.
Mare Finzi, Ke Alexander Wang, and Andrew Wilson, Simplifying Hamiltonian and Lagrangian neural networks via explicit constraints. In NewrIPS, 2020.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Newr/PS, 2019.
Nate Gruver, Mare Anton Finzi, Samuel Don Stanton, and Andrew Gordon Wilson. Deconstructing the inductive biases of Hamiltonian neural networks. In ICLR, 2021.
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse problems, 34(1), 2017.
Charles R. Harris, K. Jarrod Millman, Stéfan van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchanl, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis Oliphant. Array programming with NumPyy. Nature, (55()82)):357 362, 2020.
Pashupati Hegde, Cagatay Yildiz, Harri Lähdesmäki, Samuel Kaski, and Markus Heinonen. Variational multiple shooting for Bayesian ODEs with Gaussian processes. In UAI, 2022.
PW Hemker. Nonlinear parameter estimation in initial value problems, 1974.
Michel Henon and Carl Heiles. The Applicability of the Third Integral of Motion: Some Numerical Experiments. Astronomical Journal, 69:73-79, 1964.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In UAT, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramann.. Scalable variational Gaussian process classification. In AISTATS, 2015a.
James Hensman, Alexander G Matthews, Maurizio Filippone, and Zoubin Ghahramani. MCMC for variationally sparse Gaussian processes. In NewrIPS, 2015b.
Alejandra Hernandez and Alexander S Poznyak. Nonlinear parametric estimation of Hamiltonian systems: Identification as stabilization. Automation and Remote Control, (1(9):1611 1626, 2020.
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(5), 2013.
Hyunjik Kim and Yee Whye Teh. Scaling up the automatic statistician: Scalable structure discovery using Gaussian processes. Tn A ATSTATS, 2018.
Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9), 2021.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Peter Lynch. The swinging spring: A simple model of atmospheric balance, 2000.
Lake Metz, Daniel Freeman, Samuel Schoenholz, and Tal Kachman. Gradients are not all you need. arXin, 2021.
14

---

Published in Transactions on Machine Learning Research (03/2023)
C. Offen and S. Ober-Blöbaum. Symplectic integration of learned Hamillonian systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(1), 2022.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoît Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, bigh-ppeforrance deep learning library. NeurIPS, 2019.
Ali Rahimi Benjamin Recht. Random features for large-scale kernel machines. NeurIPS, 2007.
Katharina Rath, Christopher G. Albert, Bernd Bischl, and Udo von Toussaint, Symplectic Gaussian process regression of maps in Hamiltonian systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, (1)) 2021.
Antônio H Ribeiro, Koen Tiels, Jack Umenberger, Thomas B Schön, and Luis A Aguirre. On the smoothness of nonlinear system identification. Automatica, 121, 2020.
Steffen Ridderbusch, Christian Offen, Sina Ober-Blöbaum, and Paul Goulart. Learning ODE models with qualitative structure using Ganssian processes. JEEE Conference on Decision and Control, pp. 2896-2896, 2021.
Simone Rossi, Markus Heinonen, Edwin Bonilla, Zheyang Shen, and Maurizio Filippone. Sparse Gaussian processes revisiled: Bayesian approaches inducing-variable approximations. In AISTATS, 2021.
Hugh Salimbeni and Marc Peter Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. NewrIPS, 2017.
Rick Salmon. Hamiltonian fluid mechanics. Annuel review of fluid mechanics, 20(1):225-2555 1988.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324(5923):81-85, 2009.
John H Seinfeld. Nonlinear estimation theory. Industrial & Engincering Chemistry, 62(1):32 42, 1970.
Edward Snelson and Zoubin Ghahramani. Sparse Ganssian processes using psendo-inputs. In NeurIPS, 2006.
Yusuke Tanaka, Tomoharu Iwala, and Naonori Ueda. Symplectic spectrum Gaussian processes: Learuing Hamiltonians from noisy and sparse data. In NewrIPS, 2022.
John R Taylor. Classical Mechanics. University Science Books, 2005.
S.T. Thornton and J.B. Marion. Classical Dynamics of Particles and Systems. Brooks/Cole 2004.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. MIT press, Cambridge, MA, 2006.
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In ICML, 2020.
Cagatay Yildiz, Markus Heinonen, and Harri Lähdesmäki. Contimuous-time model- based reinforcement learning. In ICML, 2021.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ODE-Net: Learning Hamiltonian dynamics with control. JCLR, 2019.
15

---

Published in Transactions on Machine Learning Research (03/2023)
A Variational bound derivations
A.1 Standard bound derivation
We wish to minimise the KL arg min, KL (g)(u. xp) Y)] which corresponds maximising the evidence lower bound lY) 2 F(0, Z), with
Z) // q(u, Xo) log p Y. u, xo)
dudxo
(31)
quu xo)
// x()) log Y xu) p(u) p(xo) dudxo
(32)
(u) (xxo)
// (() Xo) log p( Y X )duxxo +
p(u) q(u)
((xx) log p(xx)
(u) log
du
du.
(33)
g(xo)
Fo
Fxo
The ELBO decomposes into three terms,
F F
(34)
The latter terms are KL divergences,
Fu (T.(gu))
Fx. KL[g(xo)]]pe
which be computed in closed form, since Gaussian priors p and variational approximations g for both variables. Using the fact that
pl Y)) = p(Y) u, Xo)p(w)dw
(35)
= Ep(w) ip(Y u, xo)]
(36)
can write
Fx // q(u, xo) log p( Y|u, xo)dudxo
(37)
// q(u, xo) log Eptm [p( Y ww u, Xo)] dudxo
(38)
2 Ep(w)) [ʃ q(u, xo) log Y ww u, x duudxo
(39)
Ep(w) Eg(())((e)) log W, u. Xp)
(40)
= (p(())((u)((xo))
u,)),
(41)
where applied Jensen's 's inequality in line
39 We compute this term using Monte Carlo averaging.
16

---

Published in Transactions on Machine Learning Research (03/2023)
A.2 Shooting bound derivation
We wish to minimise the KI. arg ming KL [uo(u, S) (u, S|Y)] which corresponds to maximising the evidence lower bound logp(Y) Y) 2 F(0, Z), with
F(0, = // q(u, log
p(Y, u.
dudS
(42)
q(u. S)
/ʃ q(u, log Bpw)
N
II
ply. w. u, s(())
II 1=1
p(s S1- w, w. u) dudS
(43)
Fva
-f (S) log II (s) + ((5)) log
p(so) dso + g((u) log
p(u) du,
(44)
g(so)
g(u)
Fug
Fu
where in the second line we have used the fact that,
/ p(Y, u. S. W)d
p(Y, S)
(45)
plys s(() II (( w. w, ())(W))dp((uup(s))
(46)
1=1
L II p(ss/ w.u) p(u)p(so).
E(w)
I(y. u. s(())
(47)
The latter terms in the bound are given by,
Fg = Σ H[g(s))
(48)
1=1
F., KTl(g(s)) [[))]
(49)
Fu = KL[g(u) [|p())
(50)
where IH represents the entropy, and each can be computed closed form, since we again use Gaussian priors p and variational approximations q. We can dccomposc the term Fys further,
Fya /ʃ (u, log Ep(w)
I.
Il ply W. s(()
II p(s wr, u) duds
(51)
E=1
[] q(u, log II plyd w. u s(a)
NV
L
2 Ep(w))
II p( w. u)dudS
(52)
N
/f
=F ())
si(s) )dudsu())))))
(53)
|] log logr(si(st_1, w. u)duds;dsi-1
Ep(w)
(54)
(p()))
Σ Egtsao) [lgg(y,/ u, s(a))) +
E a) [logp(s/|s:-1) w. u)]
(55)
where we applied Jenson's 's incquality in line Together the bound is given by,
52
We can compute this term using Monte Carlo averaging.
F= Jvs
Fx Fso Fu.
(56)
17

---

Published in Transactions on Machine Learning Research (03/2023)
A.3 Multiple trajectory variational bounds
To extend the bounds given in Equations [24 and 30 consider a set of K trajectories, Y = (Y,...... Yx), for ease of notation we assume each trajectory has the same number of points N, and the same observation times i, E (f1...., initial conditions, so we aim to infer the set Xo = (x1.0...., XK,0). We assume the variational posterior over each initial state factorises, so
but this is not a limitation of the model itself. Each trajectory will have distinct
K II (xx,0)
((x))
II N Mx mk,o: Q6,0)- k=1
k=1
Given this factorisation, for the standard HGP model without the shooting approximation, we obtain the bound
K
N
K
F(0, 2) = (oop))(
ΣΣ
log Oyg,(y u, u, x4.0) - -KL [a(u))l(())] - Σ KL [a(xx,))ll(xxx))] (57) k=1
where 0 = (m, Q, (ma.o: Qu, To extented the the shooting HGP to multiple trajectories, we require a set of shooting states for each trajectory, and aim to infer = .S,,..,, Sk). where we assume that each trajectory has the same number of shooting states L. We again assume that the variational posterior over the shooting states factorises, so
K ( = (Sa) = II IIII (ss.//aa. Ek.))
K
k=110
which leads the bound
K N
F(0, Z) = Eg(u)p(w))
E (( [logply w. u. sk)(
(58)
K
(o(s 1) pl. w. w)]
(59)
k=
1=0
K
[[g(s))) KLlg(s..))lI KL[g(u) [I|g())|lpu\))
where 0 = (m, Q, (ax,1, Ex,1),, and $k,(()) represents the last shooting state before time t, for trajectory k.
A.4 Model Complexities
The dominant complexity comes from the sampling of inducing points, which is O(M3) due to the requirement to compute the Cholesky decomposition of the covariance matrix. The number of evaluations of the derivative function is approximately proportional to the length of the training period T,, meaning that for D dimensional system the complexity is O(DT), due to the decoupled sampling of GPs being linear with respect to the number of inputs. For the shooting model, integration over T can be parallelised, leading a lower runtime. Large systems are likely be more complex, and require more inducing points. Therefore, the complexity with respect to the number of inducing points provides some upper limit on the system size for which we obtain good performance, although it is typically possible to use 1000s of inducing points, which would correspond to very complex Hamiltonians.
B Experiments
B.1 Hamiltonian systems
This section gives further details on the systems used for the experiments. For all systems, we scale the data so the training points have zero mean and unit variance, and compute metrics on the scaled data, so they are comparable across systems.
18

---

Published in Transactions on Machine Learning Research (03/2023)
Fixed Pendulum The fixed pendulum Hamiltonian represents the motion of a frictionless pendulum rotating about single axis, under the influence of gravity, with the coordinate y representing the angle of the pendulum to the vertical, and the coordinate representing the corresponding angular momentum. We set m = 1kg. r = Im and g 9.8lm/s. We sample initial conditions by drawing p. q ~ U(1,,1) and discard those with energy > mgr to avoid the pendulum swinging over its pivot. We use sampling rate of 8Hz for the training data and 15Hz for the test data.
Spring pendulum The spring pendulum system (Lynch 2000) is similar to the fixed pendulum, but additionally models the pendulum shaft Hookian spring, which can extend and contract as forces act upon it. In this system, Qi represents the angle of the pendulum to the vertical, Ø2 the extension of the shaft relative to its resting position, and P1 and P2 are the corresponding angular and linear momentum respectively. We set m = 1kg, r r= 3m and g m.81m/s, we set the spring constant to k = 10N/m. We sample initial conditions as q1. 92, PI. P2 ~ U(-0.25 0.25). We use a sampling rate of 6Hz for the training data, and 10Hz for the test data.
Henon-Heiles system The Henon-Heiles Hamiltonian represents the motion of star around galactic centre, with 91 and 92 representing spatial coordinates in the plane and p1 and p2 the corresponding linear momenta (Henon & Heiles) 1964). We set μ 0.8m-1 We sample initial conditions by first sampling 92. Pl,P2 P2 ~ U(-1,1) and then discarding initial conditions with E' > हह to ensure trajectories stay in the region of phase space with connected level sets, that is to say the trajectories stay localised after long times, and do not tend to infinity (Offen & Blob-Blöbaum 2022 We use a sampling rate of 4Hz for the training data, and 10Hz for the test data.
B.2 Hamiltonian aware initialisation
We form approximate data for the system's derivative function by computing the numerical deriva- tives of the trajectory data. For example using the first differences approximation we obtain Y IN Yy (Y==Y: Y2-Yz -N- NN2D We treat Y as observed data for the derivative function at input 89-th 1a-t5
locations Y E Nx2D. We wish to obtain the mean of the Hamiltonian implied by this data, at the inducing input locations Z, which we can use to initialise the variational distributions for the inducing variables. Let vec(Y) E RRNDD denote a 'flattened' Y, then this mean is given by
m = kkr (Z, YY,, Y) vea(Y)
(60)
where (Z, Y) E MX2ND the matrix from by evaluating the covariance in equation 11 for each pair of input points, and similarly Kv(Y, Y) RR2DD2ND is the covariance equation (12) evaluated at each
pair of input points.
B.3 Models
This section provides further details on the model sctup.
HGP For the HGP model we use the ARD kernel. We use a whitened representation of the inducing variables for optimisation. We optimisc the variational bound with respect to the variational paramcters 6, and the various hyperparameters of the model: o the noise variance, the kernel lengthscales and signal variance, and the inducing imput locations. We fix the shooting constraint variance as of = 1 x 10-6, and the energy constraint variance to o² 2.5 -3
For the GPODE we mirror the setup of the HGP where possible. We follow Hegde et al. (0022 and place independent GP priors over the each component of the derivative function, each with an ARD kernel and a distinct. set of inducing variables. We optimise the variational bound with respect to the same set of variational parameters and model hyperparameters. We use the same initialisation process for the inducing variational mcans as described by Hegde al. (2022)
GPODE
In practice we estimate derivatives using the gradient function, from the numpy package Harris
2020
19

---

Published in Transactions on Machine Learning Research (03/2023)
NN models For both the HNN and the NODE we use 3 hidden layers of size 256, with tanh activation. To train the NN models we form a set of sub-sequences by sliding a window of length 6 samples over the data, moving the window 1 sample at a time over the training range, e.g. for a single trajectory of length 100, we obtain 100- - 6 = 94 sub-sequences. We take these sub-sequences, shuffic them, and form them into batches. We train the model by rolling out the NN parameterised derivative function over each sub-sequence in the batch, using the first value as the initial condition, and optimising the L1 loss with respect to the training data. We use the Adam optimiser with learning rate 3 x 10-3 Because the model is trained on short sub-sequences, and is unaware that continuity is required between sub-sequences, it is liable to overfit and produce a solution that performs poorly when integrated over the entire long trajectory. To avoid this, every 10 epochs we integrate over the entire trajectory from the first point in the training data and compute the loss, we select the model with the best full trajectory loss after training has proceeded for a fixed number of iterations. For the experiments on task 1 we use a batch size of 16 and for those on task 2 we use a batch size of 32.
HGG- in relation to SSGP/SPGP To compare the our energy conserving shooting method with the training method used in the SSGP and SPGP we re-implement the training method based on short sub-sequences for our HGP model. The SSGP and SPGP use slightly different schemes, but both use a method similar to that described in the proceeding paragraph. Short sub-sequences are formed from the full, long trajectories by sliding a fixed size window over the data. We again shuffle the sub-sequences and form roll out the model and compute the bound over the short sub sequences
them into batches. In training, in each batch. This avoids the problem of vanishing/explooinng gradients, but means are optimising the model on a related but different task. A model can perform well on the short sub trajectories, but poorly when rolled out over the full trajectory. To avoid overfitting we roll out the model over the full trajectory after a set number of epochs, and compute the MNLL, selecting the model with the lowest value after a fixed number of training iterations. For the HGP-Batched model we used a sub sequence size of 6, and a batch size of 16 for all experiments. Additionally we follow (Ensinger et al.) 2022) and do not learn initial conditions, but instead use the first. data point in each sequence as its initial condition. It should be noted that Ensinger et al. (2022) use batch size of 1 and different sequence lengths (10-50) for each problem under consideration. We found that larger batch size and a shorter sequence length performed better for our experiments. The procedure for the SSGP differs slightly because Tanaka et al. (2022) only consider task 2, that is to say learning from multiple trajectories. Instead of rolling out over the full training trajectory for validation, the authors keep a set of hold out trajectories, which the use to avoid overfitting, this is of course not applicable in the case of task 1. Additionally Tanaka al. (2022) a sub sequence length equivalent 1s of data, which means different number of points for different systems depending on the sampling rate. For our implementation of the HGP-Batchhed model, we attempted to capture the key components of the methods proposed in the SSGP and the SPGP, and in order to make the comparison as fair as possible chose settings that provided the best performance on our experiments.
Symplectic integration Symplectic integrators specifically designed for Hamiltonian systems respect the conservation of the Hamiltonian during forward integration (Ensinger et al. 2022 Rath et al. 2021). This can lead to improved integration accuracy over long time scales. In our experiments we chose to use the non-symplectic Runge-Kutta 4/5 method (Dormand & Prince 1980), due to the implementation being easily available in Torch.. Our framework does not make any restrictions on the type of solver, and one could chose to use a symplectic solver if an implementation were available. It should be noted that for the HGP with shooting the choice of solver has very little effect on the solution Icarned in training, since integration is only happening over the very short shooting segments, which means the energy loss within segment is negligible, and so a symplectic solver would be of little benefit to the learning algorithm.
B.4 Implementation and hardware
All experiments were run on a MacBook Pro (14-inch, 2021) laptop with M1 Pro chip and 32 GB memory, using the CPU. We implement the HGP, and baseline models in Python, with the PyTorch framework (Paszke et al. 2019). We base our implementation on that of Hegde et al. (2022), which is available at.
20

---

Published in Transactions on Machine Learning Research (03/2023)
15
10-
-41-
-104
Figure 9: A sample trajectory showing training and test data, and model predictions for each model for single repeal of lask 1 on the FP system. For the GPODE and HGP, 32 samples from the model are shown.
https://github.
processs-o--oooss-odeeee We will on releasing the code along with the camera ready version of the paper if accepted.
C Additional results
C.1
Trajectory plots
10 and 11 show the data and model predictions for the FP, SP and HH systems respectively, for one repeat of the results shown in Table [ These plots are additionally shown in animation form in the supplementary material. Note that the animations show the noise free data for the training period, not the noisy data shown in the plots.
Figures
C.2 Mean results tables
Tables and show the equivalent of tables [ and 3 but with means and standard errors, as opposed to the median and interquartile range shown in the main lext.
State RMSE (L)
State MNLL (1)
Energy RMSE (L)
Method
FP
HH
SP
HH
SP
FP
NODE
0.16 (0.03)
0.48 (0.08)
1.04 (0.19)
0.33 (0.08)
0.02 (0.00)
1.67 (0.72)
HNN
0.38 (0.14)
1.49 (0.11)
1.40 (0.06)
0.39 (0.09)
0.07 (0.02)
1.05 (0.11)
GPODE
0.32
(0.08)
0.65 (0.13)
1.06 (0.22)
0.77 (0.69)
2.86 (1.11)
4.31 (1.54)
0.37 (0.09)
0.02 (0.00)
3.34 (1.78)
HGP
0.23 (0.04)
0.50 (0.11)
0.56 (0.09)
0.04 (0.20)
1.63 (0.89)
1.64 (0.55)
0.31 (0.08)
0.01 (0.00)
0.80 (0.08)
Table 4: Performance comparison of different methods on each system on the trajectory forecasting task.
Table shows mean and standard error over 10 repeats, for the same data as Table
State RMSE (4)
State MNLI (L)
Energy RMSE (1)
Method
FP
IIII
SP
FP
IIII
SP
FP
III1
SP
NODE
0.47 (0.03)
1.18 (0.01)
1.05 (0.13)
0.27 (0.03)
0.03 (0.00)
1.26 (0.18)
IINN
0.88 (0.12)
1.42 (0.12)
1.20 (0.07)
0.43 (0.19)
0.02 (0.01)
1.53 (0.73)
GPODE
2.48
(0.41)
3.50 (0.26)
3.38 (0.27)
8.56
(1.76
18.38 (2.55)
11.07
(1.89
4.67
(2.21)
1.42 (0.44)
3.1x
(3.0x107))
HCP
1.33
(0.13)
1.39 (0.13)
1.49
(0.11)
5.83
(1.23)
7.96 (0.69)
6.15 (0.31)
0.80 (0.41)
0.02 (0.00)
6.06 (2.52)
Table 5: Performance comparison different methods on the initial condition extrapolation task, with K trajectories. Table shows mean and standard error over 10 repeats, for the same data as Table
21

---

Published in Transactions on Machine Learning Research (03/2023)
1M 10000000000
100000000000000
MM- MMM
Figure 10: A sample trajectory showing training and test data, and model predictions for each model for a single repeat of task 1 on the SP system. For the GPODE and HGP, 32 samples from the model are shown.
MAMM
V 1
M M
MMMMMM
1
Mo0000000000
1M000100000
10000000000
1001100000
1000
1000000000
M
Figure 11: A sample trajectory showing training and test data, and model predictions for each model for a single repeat of task on the HH system. For the GPODE and HGP, 32 samples from the model shown.
22

---

Published in Transactions on Machine Learning Research (03/2023)
C.3 Comparison with SSGP
In this section we provide comparsion with the SSGP model of Tanaka al. (2022). The implementation of the SSGP provide by the authors only works for systems with D = 1, and is non-trivial to reimplement,
we are only able to obtain comparisons on the toy FP system. For the SSGP we used the same settings as
Tanaka et al. (0222) throughout.
Method
State RMSE (L)
State MNLL (L)
Energy RMSE (L)
NODE
0.12 (0.12)
0.30 (0.28)
HNN
0.20 (0.16)
0.34 (0.33)
CPODE
0.23 (0:26)
-0.03 (1.21)
0.38 (0.34)
HGP-ES
0.18 (0.19)
-0.21 (0.67)
0.25 (0.42)
SSGP
1.49 (0.39)
1.50 (0.20)
0.23 (0.32)
Table 6: Performance comparison with the SSGP on the trajectory forecasting task for the FP system.
Table shows the results of the SSGP on task 1. As expected, the SSGP model performs poorly on the trajectory forecasting lask, which mirrors our findings in our comparison with the HGP-Batched model, discussed in section 6.3 In terms of both RMSE and MNLL performance the HGP produces considerable better results. On this task the best model is the NODE, which gives slightly better results than the HGP.
Method
State RMSE (4)
State MNLL (1) Energy RMSE (1)
NODE
0.48 (0.15)
0.26 (0.06)
HNN
0.91 (0:17)
0.25 (0.15)
GPODE
2.41 (2.28)
7.11 (8:22)
6.17 (10.30)
HGP-ES ES
1.41 (0.46)
4.68 (2.95)
0.55 (0.60)
HGG--Btccedd
0.51 (0.31)
0.66 (0.51)
0.10 (0.03)
SSGP
0.75 (0.50)
(0.37)
0.13 (0.03)
Table 7: Performance comparison with the SSGP on the initial condition extrapolation task with K trajectories for the FP system.
Table shows the results of the SSGP on task 2. The SSGP model performs well for task 2 in comparison to the HGP-ES, which is perhaps not surprising given that this is the task that SSGP focused on (D = 1, multiple trajectories). The SSGP performs considerably worse than NODE in terms of state RMSE, although in terms of energy RMSE provides good result. In order to determine the cause of the performance difference between SSGP and HGP-ES we also ran the HGP model with batching. We found the HGP-Batthh provided better performance than the SSGP and HGP-ES on this task. This indicates that the performance difference between the HGP-ES and the SSGP on task can be attributed to the use of batching opposed to shooting inference, and implies that the other differences between the models (inference via inducing points vs RFFs, inferring the initial conditions, Hamiltonian aware initialisation) have a positive effect on performance, since the HGP-batched provides better results than the SSGP
Overall these results show that for systems with ==1, the HGP with shooting provides significantly better trajectory forecasting performance, whereas the batching- based inference provides better extrapolation to initial conditions, from given set of initial conditions. The results for task 2 show that shooting based inference is likely non-optimal for the problem of trajectory extrapolation with multiple trajectories. Fitting multiple trajectories jointly is a more difficult optimisation problem, since the inducing points must be effectively distributed over a larger area of phase space. Adding the shooting objectives for this problem narrows down the space of feasible solutions and makes the optimisation problem more difficult, making it harder find a good solution. We plan investigate this effect further as part of future work.
C.4 Investigating FP performance
On the FP system the HGP performs worse than the neural network baseline on both tasks 1 and 2, while this does not hold for systems HH and SP. To investigate this phenomenon, extended experiments of
23

---

Published in Transactions on Machine Learning Research (03/2023)
HGP-ES
HINN
GPODE
0.4
12
Denity
(a) Data density
I HGP-ES GPODE
0.5
HIP--ES CPOCE
04
0.4
10-
0.3
SNILL
0.3
S/NLL
0.0
-055
0.5
(b) Number of inducing points
(c) Relative noise variance
Figure 12: Effect of different variables on model performance for the FP system on task 1.
the FP system to see if the performance of the HGP is improved. We vary (i) data density, (ii) number of inducing points, and (iii) number of random Fourier basis functions.
The results for task are shown in Figure 12 We
that increasing the data density 12a provides
small improvement in the results for both the NODE and HGP, and seems to have little effect on performance for the HNN and GPODE. Increasing the number of inducing points (1bb) and basis functions 12c does not effect performance significantly for the HGP.
The results for task are shown Figure Figure 13a shows that increasing the data density for task has small positive effect on RMSE performance for the GP based models, and no significant effect for the NN based models. The effect on MNLL performance not significant. Increasing the number of inducing points and basis functions again does not have noticeable effect the performance for the HGP.
These results do not suggest the density, number of basis functions, number of inducing points is the cause for the poor relative performance of the GP models on lask 2. We believe that the poor performance on task caused by the shooting bound being non optimal in the case of multiple trajectory inference as discussed in Scction CR
24

---

Published in Transactions on Machine Learning Research (03/2023)
35
HINN
NODE
2.5
MALL
2.0
1.0
Dematy (He)
Demitty (Ho)
(a) Data density
2000
.-E HGP-FS
HGP-ES GPOOL
35
RASSE
MRLU
25
MALL
0-
15-
(b) Number of inducing points
(c) Relative noise variance
Figure 13: Effect of different variables on model performance for the FP system on task 2.
25


Reviews:
[{"review_id": "udR6lXn0L3", "reviewer": "Reviewer_xxKa", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The overall scheme is similar to the proposal of Hegde et al, where learning potentially non-conservative vector fields is considered. The two primary differences here are that (1) the authors propose considering inducing points on the Hamiltonian rather than the force field (recall, there's no Hamiltonian in the other work), (2) add energy-conservation terms at the end points of the segments in the shooting method for inference (again, only makes sense to conserve energy for conservative fields). Although the rest of the recipe is similar, the work demonstrates impressive gains in comparison to GPODE, and perhaps less convincingly with respect to other comparators.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Sincerest apologies to the editors and the authors for handing in this review late. I hope that the review is still helpful. If so, I would like to request that the authors be granted some extra time to respond to the review. I would also like to add that I am reviewing the revised version of the manuscript.\n\nThis work consider the task of learning Hamiltonian systems from a few long-horizon trajectories via Gaussian processes. It infers inducing points from observational data by variational inference; ELBO is maximized by using a shooting method to avoid blow up of gradients on long trajectories. \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Direct comparisons to the experiments in Hegde et al are difficult since the experiments therein involve non-conservative dynamics, i.e. there is no potential function. However, in the present context, it would helpful to know which of the two primary modifications (direct prior over Hamiltonian, or energy constraints in shooting) is responsible in greater part for the empirical gains overs GPODE. \n\nDisentangling contributions -- neural net vs GP, energy-conserving, imposing Hamiltonian -- in empirical terms would certainly add value to the present work.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "HD5G3oUPN2", "reviewer": "Reviewer_GDmG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n- The paper addresses a problem that might be of broad interest in machine learning and applied mathematics.\n- In general, the paper is well-written with appealing visualizations.\n- The technical developments of the work are rigorous and intuitive. Background information and relevant literature are adequately covered.\n\n\n### Weaknesses\n\n1. In my opinion:\n    - The paper overstates the significance of its technical contributions to inference methodology.\n    - The current presentation of the paper doesn’t highlight well enough its connection to the closest reference, Hegde et al. (2022). This might give the readers a false sense of novelty.\n    \n    As pointed out by the authors (but very late into the manuscript), in terms of methodology, the work depends heavily on Hegde et al. (2022), which introduces variational multiple shooting for general ODEs using Gaussian processes. As far as I can tell, the technical constructions of the two papers are the same, and the only major algorithmic adjustment of the paper is to hardcode the conservation of energy into the distribution (Equation 27-28) to specifically target Hamiltonian systems.\n    \n    While this idea is intuitive and reasonable, given the close relation of the two frameworks, I suggest that the authors:\n    \n    - remove their claims of significance to the inference methodology\n    - restructure the framing of the work to explicitly acknowledge Hegde et al. (2022) (early, in the Methodological sections, not Related works) as the foundation of the framework. I also suggest that the authors highlight in each part of the framework the deviations of the methods from Hegde et al. (2022), and possibly explain why such changes are made/necessary.\n\n\n2.  Experiments: I find the experimental evaluations of the work unconvincing:\n    - Since the paper follows the framework of Hegde et al. (2022) with an additional hardcoding of the law of conservation of energy, it is expected that the method will outperform Hegde et al. (2022) for Hamiltonian systems. The main question for the evaluations is whether the method is comparable with pre-existing Hamiltonian Gaussian Process baselines.\n    - It is thus disappointing that such comparisons are not available:\n        - “A public code exists for the SSGP, however at the time of writing it does not support experiments with a single trajectory, or systems with D > 1.” Why don’t we test and compare the performance on multiple trajectories, or for D=1?\n        - “Code for the SPGP is not publicly available”: It seems that the paper SPGP stated that code will be published upon request. Have the authors reached out to SPGP’s authors for access?\n        \n        I’m aware that the paper attempt to use HGP-Batched to emulate SSGP/SPGP, but comparisons with a proxy is less convincing and should only be done if there are no alternatives.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The manuscript presents a Gaussian process model for Hamiltonian systems with decoupled parameterization and introduces an energy-conserving shooting method that allows inference from dynamics trajectories. The paper then compares the performances of the approach in learning Hamiltonian systems with existing methods in various data settings.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- [Optional] Inducing variables and inducing points (Section 3.2) are mentioned as an established concept. I suggest adding some relevant background to make it accessible to a general audience with less exposure to the Gaussian process.\n- [Typo] Page 2: In this work\n- Figure 1: There shouldn’t be an abbreviation in the legend. What are inducing energies?", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6WasTmX0la", "reviewer": "Reviewer_14sk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "S1. To my knowledge, this is the first example of modeling a Hamiltonian using decoupled sampling.\n\nS2. The variational inference algorithm with multiple shooting is practically useful.\n\nS3. The evaluation was conducted for two tasks, i.e., forecasting and simulation.\n\nS4. The manuscript is readable.\n\nW1. Experiments have not shown the benefits of using decoupled sampling.\n\nW2. There are several unclear points regarding the experimental setup.\n\nW3. The advantage of the proposed method is unclear.\n\nW4. There is insufficient discussion of the advantages and adaptability of the proposed method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper extends [Hegde et al., 2022] to handle Hamiltonian dynamics. The Hamiltonian is modeled by the decoupled representation of Gaussian processes, which allows one to efficiently sample the Hamiltonian vector fields from the approximated GP posterior. The model learning is based on variational inference with multiple shooting schemes. This scheme can be used for effectively learning dynamics from long trajectories. The experiments on several physical systems show that the proposed method outperform the baselines (i.e., NODE, HNN, GPODE) for the HH and SP systems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Compared to standard RFF, decoupling sampling has the advantage of being able to properly evaluate predictive variance. Can this be demonstrated in an experiment? Is it possible to visualize and analyze the variance in phase space, for example?\n\n- HNN uses a learning method based on finite differences, so a comparison with SymODEN, which can be trained using ODE solver, would be appropriate. Would it be possible to add a comparison with SymODEN?\n\n- I could not find the reason for the poor performance of HGP in the experiments on Pendulum. Would it be possible to analyze the results of the proposed method from the following perspectives.\nFor example,\n  - What happens to the performance when the temporal resolution of the trajectory is low?\n  - Is the number of inducing points and basis functions sufficient?\n\n- Can the proposed method be applied to large systems? It would be good to add a discussion on the computational cost in that case.\n\n- A learning method using symplectic integrators has been proposed to make energy conservation more rigorous, and I think it would be good to mention that and discuss whether the proposed method can take advantage of it.\n\n- It is also good to mention whether it can be applied to systems with friction.\n\n- I am having trouble understanding the relationship between the variables s,y,x. Please clarify.\n\n- Please explain in detail how to sample the initial conditions when generating data in an experiment. Do you limit the energy range, etc.?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors have already done a good job in the initial submission and in the responses tot he referees. They are welcome to works those responses into the final version of the paper.
Decision: Accept (poster)