Learning interacting dynamical systems with latent Gaussian process ODEs
Authors: Cagatay Yildiz
Melih Kandemir
Barbara Rakitsch
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian Processes
Ordinary Differential Equations
Interacting Dynamical Systems

Abstract: We study uncertainty-aware modeling of continuous-time dynamics of interacting objects. We introduce a new model that decomposes independent dynamics of single objects accurately from their interactions. By employing latent Gaussian process ordinary differential equations, our model infers both independent dynamics and their interactions with reliable uncertainty estimates. In our formulation, each object is represented as a graph node and interactions are modeled by accumulating the messages coming from neighboring objects. We show that efficient inference of such a complex network of variables is possible with modern variational sparse Gaussian process inference techniques. We empirically demonstrate that our model improves the reliability of long-term predictions over neural network based alternatives and it successfully handles missing dynamic or static information. Furthermore, we observe that only our model can successfully encapsulate independent dynamics and interaction information in distinct functions and show the benefit from this disentanglement in extrapolation scenarios.

Learning interacting dynamical systems with latent Gaussian process ODEs
Çagatay Yıldız"
Melih Kandemir
University of Tübingen cagatay. yildizluni www.:yldiezu--uuueeeeee...... tuebingen.de
University of Southern Denmark kandemir@imada sdu. dk
Barbara Rakitsch Bosch Center for Artificial Intelligence barbara..rkiitschheee bosch. com
Abstract
We study uncertainty-aware modeling of continuous-time dynamics of interacting objects. We introduce single objects accurately from their interactions. By employing latent Gaussian pro- cess ordinary differential equations. our model infers both independent dynamics and their interactions with reliable uncertainty estimates. In our formulation, each object represented as graph node and interactions are modeled by accumulating the messages coming from neighboring objects. We show that efficient infer- ence of such a complex network of variables is possible with modern variational sparse Gaussian process inference techniques. We empirically demonstrate that our model improves the reliability of long-term predictions over neural network based alternatives and it successfully handles missing dynamic or static informa- tion. Furthermore, we observe that only our model can successfully encapsulate independent dynamics and interaction information in distinct functions and show the benefit from this disentanglement in extrapolation scenarios.
model that decomposes independent dynamics of
Introduction
broad spectrum of dynamical systems consists of multiple interacting objects. Since their interplay is typically priori unknown, learning interaction dynamics of objects from data has become an
emerging field in dynamical systems [1, 2, 3]. The ever veerggoowing interest in interaction modeling is due to the diversity of real-world applications such as autonomous driving [4], physical simulators [5], and human roooot interactions [6]. Standard time-series algorithms or deep learning approaches (e.g. recurrent neural networks), that have been designed for single-object systems, do not scale to large number interacting objects since they do not exploit the structural information the data.
In recent years, graph neural networks (GNNs) have emerged promising tool for interactive systems, where objects are represented as graph nodes. State-of-the-at methods learn interactions by sending messages between objects in form of multi-layer perceptrons [1] or attention modules [7]. These methods yield highly flexible function approximators that achieve accurate predictions when trained large-scale datasets. However, their predictions come without calibrated uncertainties, hindering their reliable implementation for uncertainty-an -awar aware applications.
In contrast, Gaussian processes (GPs) are well-known for providing calibrated uncertainty estimates. They have been successfully employed on discrete time-series data [8, 9, 10] and, more recently,
Work performed during internship at Bosch Center for Artificial Intelligence
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Initial m observations
Encoders (arnortized posteriors)
Integration (forward pass)
Predictions FEM
Glabal Bucure
oe sylo 1222
Kioya)
(((()))))))
Σ (t),,)
Gauseian process postetiions
JeN.
f. ~ plt U.))
4 ~ p(& Up)
Interaction differential fonction
d (+) (1) + Σ (b200,b''
Figure I: An overview of our predictive model. (Top-left) An input bouncing balls sequence with four balls, which move independently other than collision (interaction) times. An encoder is used to extract initial values and global latents. (Bottom-left) The differential function formed by sampling from the GP posteriors the independent kinematics and interaction functions. (Middle Given the samples, predicted trajectories are computed using the forward integration of the differential function.
continuous-time generalizations of these methods [11, 12, 13]. Importantly, none of the these works adresses dynamical models for interacting systems. While it is possible to study each object in isolation, ignoring the interaction effects might lead to inaccurate predictions.
In this work, we address the shortcomings of both model families by presenting an uncertainty-aware continuous-time dynamical model for interacting objects. Our formulation decomposes the dynamics into independent (autonomous) and interaction dynamics. While the former governs the motion of an object isolation, the latter describes the effects that result from interactions with neighboring objects. For successful uncertainty characterization. we propose to infer the unknown independent and interaction dynamics by two distinct GPs. We demonstrate that having function-level GP prior on the individual dynamics components the key to successfully disentangling these dynamics, which in turn allows for interpretable predictions and leads improved extrapolation behavior.
We employ latent Gaussian process ordinary differential equations (GP-ODEs) for dynamics learning. allowing to learn complex relationships between interacting objects without the need of having access to fully observed systems. Thanks to recently proposed decoupled sampling scheme [14], the computational complexity of our model scales linearly with the number of time points at which the ODE system is evaluated. As a result, our algorithm scales gracefully to datasets with thousands of sequences. To demonstrate the benefits of our framework, exhaustively test our method range of scenarios varying in function complexity, signal-to-noise ratio, and system observability. Our model consistently outperforms non-interacting dynamical systems and alternative function approximators such as determinissic//aaessian neural networks.
wide
2 Background
In this section, give background continuous-tims systems and Gaussian processes. Both together form the backbone of our uncertainty- framework for interactive dynamical systems.
2.1 Continuous-time Dynamical Systems
Contimuous-time dynamical systems are often expressed using differential functions x(1) = 4x1)= f(x(t)), where x(t) E D represents the state of an ODE system at time l and and f: R RD the time differential function governing the dynamics evolution. The state solution x(1) at an arbitrary time & is characterized by the initial value at time point to and the differential function:
x(t1) x((t) 1/ f(x(T)) dr.
Existing work aims to approximate the unknown differential by Gaussian processes [15, 11] or neural networks [16]. These methods have shown to accurately capture the dynamics and outperform

---

their discrete-timg counterparts in a wide range of applications such as time series forecasting [17], classification [18] or reinforcement learning [19]. Furthermore, ODE models allow easily inject domain knowledge into the system, enabling interpretable and flexible hybrid models [20, 21, 22].
2.2 Gaussian Processes
Gaussian processes (GPs) define priors over functions [23]]
f(x) ~ (P(((x), k(x, x'),
where f :R'" + R maps D-dimensional inputs into one-dimensional outputs. GPs are fully specified terms of their mean and their covariance:
[[](x)) = (xxx,
cov covl [[(x) = f(x,x'),
where μ IRD + is the mean and A RD x RD + R is the kernel function. GPs can be treated as an extension of multivariate normal distribution infinitely many dimensions, where any fixed set of inputs X E IR ND follows the Gaussian distribution
pff = Wf( Kxx),
(I)
where the mean function 14 evaluated inputs X, and Kxx the kernel function evaluated all input pairs in X. While GPs provide a natural mechanism to handle uncertainties, their computational complexity grows cubically with the number of inputs. This problem is often tackled by sparse
GPs, which rely on augmenting the GP with M inducing inputs Z = [zi,... 2 2 E IR J x and
corresponding output variables u = [u,,.... "M]" E Mx1 with Uym = f(mm) [24, 25]. Assuming the commonly used zero-mean prior, the conditional distribution over S(X) follows the GP:
p(f u) = (f KxzKzz Kxx KxzK2i
(2)
where Kzz the covariance between all inducing points Z. and Kxz between the input points X and the inducing points Z. The inducing points can thereby be interpreted as compressed version of the training data in which the number of inducing points M acts as trade-off parameter between the goodness of the approximation and scalability.
In this work, we employ the squared exponential kernel k(xx, x') o² exp (- ΣD
where za denotes the d-th entry of the input x, x² is the output variance and La is the dimension-win -wise lengthscale parameter.
3 Interacting Dynamical Systems with Latent Gaussian Process ODEs
In Sec. 3.1, we describe our continuous-time formulation for systems of interacting objects. decomposes the dynamics into independent kinematics and an interaction component that takes the interactions to neighboring objects into account. Placing a GP prior over the individual components essential in order to arrive at (i) calibrated uncertainty estimates and (ii) disentangled representations as we later on also verify our experiments. In Sec. 3.2, we embed the GP dynamics into a latent space that can accomodate missing static or dynamic information. Both together allows the application of our continuous-tie interpretable dynamics. We conclude this section by our variational inference framework (Sec 3.3) based on sampling functions from GP posteriors.
formulation to a wide range of scenarios and allows for learning
3.1 Interacting Dynamical Systems
We assume dataset of P sequences y (YI...., YpJ. where each sequence Y = YI.N = y 1:4 composed of measurements of objects at time points T tt Without loss of generality, we assume that the measurement (t,,,) y R is related to the physical properties of object a, such as position and velocity, which can routinely be measured by standard sensors. The dynamic state of object at any arbitrary time is denoted by a latent vector ht) E IR' which does not necessarily live in the same space as the observations. We furthermore assume that each object is associated with a global feature vector E RC, which corresponds to the static attributes that remain constant over time. Finally, we denote the concatenation of all states by H(t) (t), (t)] [2 and all globals by C (c..... A E I8 AxC.
3

---

In the following, we propose to disentangle the complex continuous timee dynamics into independent kinematics and interaction differentials. More concretely, we introduce the following dynamics:
dt d H()) dt
d h h0))
h (4) ..... dt
(3)
(t) = f. (t) E fo (h"(t, (t), (0).).c))) a'EN.
(4)
dt
where N₂ denotes the set of neighbors of object in given graph. The first function f. RD4C Rº_ models the independent (autonomous) effects, which specifies how the object would behave without any interactions. The second function fo 222 + _D models the interactions by accumulating messages coming from all neighboring objects. Since message accumulation is the de-facto choice in interaction modeling [2, 26], the additive form of the differential function is very generic inductive bias.
Our formulation models interactions between pairs of objects explicitly via the differential equation (Eq. (5)). Higher-order interactions are taken into account via the continuous-time formulation that allows information to propagate through the complete graph via local messages over time such that the state of the object h,; can also depend on objects that are not directly connected in the graph. In contrast discrete formulations, for which the message passing speed is limited by the sampling rate, our continuous- time formulation enjoys instant propagation of information across objects. Finally, please see Section for an investigation of our interaction component under kernel perspective.
Remark-1 In Sec ?? we demonstrate two straightforward extensions of our formulation with non-linear message accumulation, which empirically show to have gain our formulation.
Remark-2 Previous GP-based ODE methods [15, 11, 12] assume black-box approximation the unknown system Ht) (()())) whereas our state representation gracefully scales a varying number of objects and also include global features.
3.2 Probabilistic Generative Model
Real-world data of interactive systems necessitates embedding the dynamics into a latent space in order allow for missing information; the observations may contain only partial estimates of the states and the globals C might not be observed at all. We account for those circumstances by treating the states H and the globals C as latent variables, leading following generative model (see Figure 1):
h ~ N(O,I), c² ~NO(0, fb(-)~ GP(O,Ao(:.)) h = hf + L f (T)()))
f(-) (P(O,,a(.:))
Σ fb (1"(7). (T) dr, d'EN.
(T),c"
(5)
yn ~ ply h(),
where we introduced standard Gaussian prior over the initial latent state, and assumed that the data likelihood decomposes across time and objects. furthermore model unknown functions f, and f. under independent vector-valued GP priors.
In our experiments. we further set p(y* from the latent to the observational space and o² E RQ is the noise variance. We further fix B = [I. 0] where 1 0 E R in order to arrive at an interpretable latent space in which the first dimensions correspond to the observables. This assumption is fairly standard in the GP state-space model literature since more complex emission models can be subsumed in the transition model without reducing the model complexity [27].
-y Bhg diag(o?)), where B Ox maps
Modeling partially observed systems often leads to non- con-identifiability issues that hamper optimization and ultimately lead to deteriorated generalization performance. One way counteract this behavior is to inject prior physical knowledge into the system by decomposing the state space of each
object he (t) = (t), (t)] into position (t) and velocity (t) components [28]. Using elementary physics, the differential function has then the form of d di (t) [ (t), de [t)] with
di vl) = f, (h (1)()) (1),c9) + EwEx. to (1)(1), ha (1), c² )

---

+1
Ball 42
e3
+1 e2 e3
y-g0s
Time (5)
Time
Time
Time
Time
Time
Test predictions for I-GPODE (left) and I-NODE (right) on bouncing balls dataset. The solid curves are the groundtruth trajectories and the shaded regions denote the predicted 95% confidence intervals. I-GPODE (ours) yields better calibrated long-term predictions than I-NODE. Additional results can be found in Figure ??.
Figure
Remark Unlike previous work, our formulation incorporates global features that modulate the dynamics. In many applications such as control engineering and reinforcement learning, the dynamics modulated by external control signals [29] which also be incorporated into framework.
3.3 Variational Inference
Next, we derive an efficient approximate inference scheme that provides high level of accuracy and model flexibility. In the above formulation, the model unknowns are the initial values H2 h global variables C I:A and the differentials f, and f.. Since the exact posterior p(f,, f.. H, C in non-linear ODE models is intractable, opt for stochastic variational inference [30]. We first describe the form of the approximate posterior and then discuss how to optimize its parameters.
Variational family Similarly to previous work [16], we resort to amortized inference for initial values and global variables, H, ~ 1s (H and, C ~ 4o (,(C) (C 1:N), where Φ is the parameters of neural network encoder that outputs Gaussian distribution with diagonal covariance. For the unknown time differentials f = {fs, fb ). we follow the standard sparse GP approximation in which each output dimension E DJ has its own independent set of inducing values Ua.di Us.d E R and kernel output variances o² o,. RR,. Inducing locations {Z.. Za) and lengthscales {La, 6 shared across output dimensions. Denoting the collection of inducing points by U {U,, Ub}. we arrive at the mean-fice field variational posterior:
(UU) = II N(U..a m.,d. S.,a) mi,d, St.a). d=1
where the means (m.,d, M and the covariances (S..d, S are free variational parame-
ters. Putting everything together, variational approximation becomes as follows [31]:
(H, C, U) = qH,, (f)((( U.)p(f. Ut),(U),
where p(f, and pff, U,) follow Eq. (2). Our variational family makes two assumptions that fairly standard in the (deep) GP literature (c.g. [32]): (i) we apply the same independence assumptions in the approximate posterior as in the prior resulting mean-field solution, and (ti) we assume that the inducing outputs U capture the sufficient statistics of the training data allowing the use of the prior pff U) in the approximate posterior.
Variational bound We then seek to optimize the parameters of the approximate posterior by maximizing lower bound to the evidence [30]:
logp())) > / qHH,, C, C, log P(Y, H,, C,C.U)
H C, dH, dCdfdU.
In the following, we detail its computation for single data instance YI,N, omitting its generalization to multiple sequences for the sake of better readability,
log p( (I,n)= 2 F. [logp(( Y IN H,, C, U)] (I.(((H,))l((H)) ((C(HH))))r((((C))) kI.[(())]p(()) ([U))))))))))))
5

---

where KL denotes the Bellback--eibler divergence.
Likelihood computation via decoupled sampling from GP posteriors Computing the condi- tional log-likelihood Deg((Y1s) 1:N HH, C, f, U) entails a forward pass in time (Eq. (5) which can | be done with any standard ODE solver. The difficulty lies in marginalizing over the approximate posterior of the initial latent states (H,,), global variables (C), and the GP functions q(f. U). Each marginalization step alone is already analytically intractable, let alone their combination. We therefore opt for Monte Carlo integration which gives us an unbiased estimate of the expected log-likelihood. We start by drawing L samples from the approximate posteriors
H g0 (H, (HHHYYLNN) IIN), CO go (C)(Y)):)
U ~ ((U), {))
(7)
where 1 denotes the sample index and (.) is a function drawn from the sparse GP posterior. Sampling from the GP posterior naively scales cubically with the number of data points. Moreover, since we do not know a-priori on which points the ODE solver evaluates the function, have to sequentially draw points from the posterior. While this can still be done cubically in time by performing low-rank updates, it often leads numerical instabilities for small step sizes. To overcome this challenge, we resort to the decoupled sampling scheme proposed in [14], where we first draw the prior samples from set of random Fourier features and then update them using Matheron's rule to obtain posterior samples. After having sampled the quadruple via Eq. (7), we can compute the trajectory H log-likelihood becomes
would
1&w deterministically by forward integrating Eq. (5). Monte Carlo estimate of the
1 Σ log logp(y>
By [log p Y 1:N H2, C,f,U)
i,na
where the log-likelihood term decomposes between objects and between time points, enabling doubly stochastic variational inference [33].
KL regularizer The prior distributions over the inducing variables follow the Gaussian process p(U) = p(U,)p(UU,) with p(U.) and p( Un) following the equivalent of Eq. (1). Since we further- more assumed standard Gaussian prior p(H,) dimensions for the initial values, all KI. terms can be computed in closed form.
NO, I) and p p(C1) (O, with suitable
Computational complexity Evaluating the differential function costs ((M* + 1)) for graph with A objects, / interactions and M inducing points. The term (MM2 stems from the independent kinematics, the term ((M² from the interaction effects. As a consequence, forward integration using RK4 solver takes (A + I)) time, where T is the number of time points.
3.4 Comparison to Standard GPODEs
Our approach enhances the capabilities of the GPODE model family in the following three aspects:
Explicit modeling of interactions Standard GPODEs model interactions by allowing the time differential to take the whole state vector H(t) as input and to learn one independent GP for each latent dimension. This shared latent space assumption entails three major drawbacks: (i) obligation to fix the object count as a model hyperparameter, (ii) dependency of the learned model on predefined ordering of the objects in the scene, (iii) inevitable growth of the latent dimensionality proportional to the object count. The latter sets a severe bottleneck especially for GP modeling as the performance of many kernel functions in widespread use are highly sensitive to input dimensionality. (For example, on the bouncing ball dataset with N, balls and D latent states per object, GPODE needs to learn a latent function with N. --dimension inputs and outputs.) In contrast, I-GPODE needs only to learn two functions, the independent kinematics f. and the interaction function f, whose input sizes scale independentl of N.. Our Table indicates that learning interaction dynamics without the strong inductive bias of our model is difficult and the GPODE model chooses to stay at the prior instead leading to deteriorated MSEs and ELLs.
Disentangled representation We infer object-specific latent variables that modulate the dynamics, which allows model to disentangle the dynamics from static object properties (e.g. charge information). The interpretability and likely physical correspondence of the disentangled factors have the potential facilitate the use of our approach in transfer learning and explainable Al applications.
Inference of latent state dynamics We perform the learning in a latent space, where the initial value of a trajectory is given by an encoder, leading to the following two advantages: First, the

---

Table 1: A taxonomy of the state-of-the-are ODE systems with NN/GP dynamics, with/without interactions, and with/without latent dynamics. These methods also form the baselines we compare with [neural ODE (NODE), Bayesian NODE (BNODE), latent dynamics (L), interacting (I)).
Function uncertainty
Interacting
Latents
Reference idea
Abbreviation
X
Neural ODE [16]
NODE
X
X
/
Latent NODE [17]
NODE-L
X
GDE [34]
I-NODE
/
X
V
LG-ODE [35]
I-NODE-L
X
X
Bayesian ODE [11]
GPODE
X
GP-SDE [13]
GPODE-L
X
V
ODE [28]
BNODE-L
V
Our work
I-GPODE
/
V V
Our work
I-GPODE-L
dynamical system and the data points may live in different spaces, which facilitates learning from high-dimensional sequences. Second, Bayesian modeling of state dynamics on a latent space enables reliable quantification and principled treatment of sources of uncertainty, such imprecision of modeling assumptions, approximation error, and measurement noise.
4 Related Work
GPs for ODEs GPs for modeling ODE systems have been studied in a number of publications (e.g. [11, 12, 13, 15, 36)). With the notable exception of [13], they only consider systems in which the dynamics are defined in the data space. The work that is closest to ours from a technical perspective is [11, 12] that also employ decoupled sampling in order compute consistent trajectories during inference. are not aware that interacting dynamical systems under a GP prior have been studied previously, either in the continuous or in the discrete time setting.
Neural ODEs for dynamical data Since the debut of Neural ODEs (NODEs) [16], much progress has been made on how to model sequential data in the continuous time domain using neural networks. These works [17, 37] assume a black-box approximation to unknown ODE systems d H(t) f (H))). where f a deterministic neural network, and the latent space separated from the observational space using neural network decoder. A subtle difference in our approach is that we decided use linear mapping instead. However, when the outputs are nigh-dimensional. e.g. video data, this can be easily changed. Few works have integrated function uncertainty into NODEs by putting a prior over the neural network weights [28, 38). To the best of our knowledge, none of these works addressed interacting systems. However, we still compare against an interactive adaptation of these methods.
Modeling interacting dynamics Interacting dynamical systems have first been considered for discrete time-step models and for the deterministic setting [1, 26, 39]. Many of these discrete formulations can be transferred to the continuous case as shown in [3. 34, 40]. This also holds true for our approach, for which the discretized version of the dynamics (Eq. (3), (4)) can be easily cast into one of the existing frameworks (e.g [39]). These works have also been extended to the probabilistic context using a variational auto-encoder [2, 35]. The hidden variables are thereby used to either encode the initial latent states or static information. In contrast to our work, none of these approaches allow for function uncertainty in the dynamics. Finally, [41] proposes a symbolic physics framework, differing from our method in its search-based fixed grammar describing the dynamics.
Finally, provide summary of related techniques and derived comparison partners for experiments in Table
5 Experiments
We compare our approach against state-of-the-art methods a large number of scenarios that differ in function complexity, signal-to-nois. ratio, and observability of the system. The empirical findings
7

---

suggest that our model leads to improved reliability of long-term predictions while being able to successfully encapsulate autonomous dynamics from interactions effects. In all our experiments, we use the RK4 ODE solver and ACA library [42] for integration and gradient computation. Due to the space limit, we refer to the Supplementary Material for more detailed information about the experimental setup and comparison methods (see also Table for an overview). Our Py PyTorch [43] implementation can be found in https: /tte:///tttbb com/ wwwwooeeeewrcciee....ooeee (GNU AGPL v.3.0 license).
5.1 Experiment Details
Datasets We illustrate the performance of our model two benchmark datasets: bouncing balls [44] and charges [2]. These datasets involve freely moving N₂ = 3 balls that collide with each other and N. = 5 particles that carry randomly assigned positive negative charges leading to attraction or repulsion, respectively. All simulations are performed in frictionless, square boxes in 2D. We generate 100 bouncing balls training sequences with different levels of Gaussian perturbations to simulate measurement noise (see Supplementary Section ?? for details). Since the charges dataset requires inferring the charge information, we use 10k train sequences without observation noise as in [2] and similarly use 500 training sequences when velocity information is missing.
Partial observability Prior works typically evaluate their methods on datasets with position and velocity observations [1, 2]. However, having access to the full state of an object is for many real- world problems unrealistic, e-g. many tracking devices can measure positions, but cannot measure velocity or acceleration. To test the model performance in such scenarios, we consider an additional bouncing balls dataset in which the velocities are not observed. On the charges dataset, we assume that position and velocity are observed, but treat the charge information as missing.
Reported metrics We quantify the model fit by computing the expected log- logllikel (ELL) of the data points under the predictive distribution. Further, we report the mean squared error (MSE) between ground truth and predictions over all predicted time steps and over all objects (see Supplementary Section ?? for the exact definitions). Each experiment is repeated five times and we report the mean and standard deviation of both metrics on test sequences.
5.2 Empirical Findings
Due to our latent variable construction, ODE state dimensionality D could be arbitrary even though the observations are -ou--diensional. To choose an appropriate state dimensionality, we study if the model performance can be improved by augmenting the state space with auxiliary dimensions (D >> 4). Table ?? shows that increasing the model flcxibility beyond need leads to overfitting as we observe lower training but significantly higher test error. Consequently, if not stated otherwise, we use a four-dimensional latent space for each object that corresponds to position and velocity in x and y coordinates and observations consist of their noisy versions. Next, we discuss the main findings.
5.2.1 Interacting dynamics superior over standard formulation
We consider three bouncing ball datasets with varying noise levels to reflect dif- ferent levels of problem difficulties. To demonstrate the merits of our decom- posed formulation in Eq. (4), we com- pare against a standard, GP-based non- interacting dynamical model (GPODE) in which the time differential takes the whole state vector Htt as input [11]. As shown in Table 2, our interaction model consistently outperforms its stan- dard counterpart irrespective of the noise level and function approximator. We also
Table 2: Comparison of standard and interacting GP-ODE systems on bouncing ball datasets.
NOISE LEVEL
MODEL
MSE +
ELL t
No NOISE
I-GPODE GPODE
17.3 + 1.2
-25.9 +: 3.6
23.9 + 0.9
-597.6 + 28.5
LOW NOISF
I-GPODE GPODE
17.8 + 0.9
-26.5 + 4.5
23.6 + 0.7
5800.8 + 19.8
HIGH NOISE
I-GPODE GPODE
18.4 0.9
29.44
+ 27.5 note that the results are consistent when we replace the GP approximation with deterministic and Bayesian neural networks (NNs), which indicates the robustness of the inductive bias. (see Table ??).
24.5 + 1.0
-569.3
8

---

5.2.2 GP approximation yields more calibrated uncertainties
Next, we compare our GP interaction model (I-GPODE) against neural net- work variants with/without uncertainty (I-NODE, I-BNODE), obtained by re- placing GPs with NNs and Bayesian NNs. The findings in Table strongly indicate that I-GPODE learns better cal- ibrated uncertainties, measured by ELL., than I-BNODE and I-NODE (see Fig- ure 2 and ?? for visual illustrations). terms of MSE, observe that I-NODE tends to attain the smallest values. How- ever, increasing noise levels deteriorates the performance of I-NODE and eventu- ally the results become comparable I-GPODE. The results stay consistent
Table Comparison of our GP approximation against alternative function approximators.
NOISE LEVEL
MODEL
MSE 4
ELL
I-GPODE I-NODE
17.3
2 3.6
No NOISE
8.5
-61.9 8.6
I-BNODE
2.6
-37.8 = 4.2
I-GPODE
1..
-6..5
LOW NOISE
I-NODE
13.3
-96.2 12.5
I-BNODE
2..0
-36.4 2.8
I-GPODE I-NODE I-BNODE
0.94
-91.4
HIGH NOISE
15.6
-128.9 + 23.2
20.8 0.7
-47.6 10.5
when we repeat the experiment with two balls (see Table ??).
5.2.3 GP approximation enables disentangled function learning
In the following, study whether the estimated functions can disentangle inde- pendent kinematics from interaction ef- fects. We train I-GPODE, I-NODE and I-BNODE dataset with three balls and evaluate the trained independent dy- namies function f, test dataset with one ball (since the test dataset incorpo- single object, the dynamics do not involve the interaction function f,). Three test sequences as well as the inde- pendent dynamics function predictions are illustrated Figure 3 (see Table ??
GPODE BNODE
Disentanglement between independent kine- matics and interactions. Our model variants are trained on a bouncing ball dataset with three balls and tested with single ball using the independent dynamics function f, only. We show the mean predictions with circles darkened time lapses.
Figure 3:
for quantitative comparison and Figure ?? for additional illustrations). We observe that I-NODE predictions tend to deviate from the test trajectory more quickly compared to I-GPODE predictions. We conjecture that this behaviour is because neural networks are overflexible and thus the learned functions may not necessarily decompose independent kinematics from interaction effects whereas the function-level regularization of I-GPODE helps with disentanglement.
5.2.4 Structured dynamics improve latent dynamics learning
In our last bouncing ball experiment, move to setting in which the velocities are no longer observed. First, keep the velocities la- tent states and contrast two variants of our model, i.e. with structured la- tent space (I-GPODE-L-S) and with unstructured latent space (I-GPODE- L). As shown in Table 4, inject- ing strong prior knowledge helps in this challenging setting in which the first order model clearly fails (see also Table ?? for training results).
Table 4: Comparison of I-GPODE variants on bouncing ball dataset without velocity observations. The suffices "-L" and "-S" stand for latent variable model and structured state space.
NOISE
MODEL
MSE
ELL 1
I-GPODE
24.1 + 1.2
-300 +: 23
No NOISE
I-GPODE-L.
23.3 + 0.8
-176 +: 301
I-GPODE-L-S
18.2 + 1.0
-33 2
I-GPODE I-GPODE-L
23.3 0.8
-239 +: 11
Low NOISE
47.0
-761 + 922
I-GPODE-L-S
2.9
-76 + 22
Finally, we compare I-GPODE-L- S with I-GPODE, which drops the velocity component from the latent space (hence learning the dynamics
I-GPODE I-GPODE-L I-GPODE-L-S
23.8 0.7
-254 + 16
HIGH NOISE
119.3 94.3
-90 982
21.8 $1.80.7 0.7
-115

---

in the data space). As demonstrated in Table 4, I-GPODE is clearly outperformed by I-GPODE-L-S. can thus be suggested that our latent variable construction necessary in presence of missing states.
5.2.5 Global latent variables boost performance
In the final part of the experiments, we consider more challenging dataset of charged particles. Since the dynamies modulated by unknown charges, we turn to our global latent variable formula- tion. In other words, our learning task be- comes simultaneously inferring the func- tions f, and f,, the initial latent states H1. as well as a latent variable ca R associated with each observed trajectory yf.N- To form the upper and lower per- formance bounds, include two base- lines in which the charges are either ob- served or completely dropped from the model description. The results are shown in Table 5. We notice that the structured state space formulation boosts the performance of I-GPODE. However, the effect is less pronounced compared the previous setting in which dynamic information is missing. See Table ?? for more results with different global variable encoders.
Table 5: Comparison of interacting GPODE and NODE
systems on the charges dataset (no noise).
MODEL
MSE +
ELL f
WITHOUT GLOBAL LATENTS
I-GPODE
19.2 + 0.9
-171 + 7 -460 = 36
I-NODE
14.1 + 0.8
I-GPODE I-NODE
15.4 + 2.3
172 = 88
WITH GLOBAL LATENTS
9.9 19::00. 0.6
282 -
I-GPODE-S
12.9 + 2.1
177 = 60
I-NODE-S
9.7 = 0.2
282 + 8
I-GPODE
10.7 + 1.1
99 9
GLOBALS OBSERVED
I-NODE
7.5 1.2
-148 - 27
6 Discussion
We have presented the first uncertainty-aware model for continuous-time interacting dynamical systems. By embedding the dynamics into latent space, formulation yields a flexible model family that can be applied to a variety different scenarios. In our experiments, we found that our approach leads to improved disentanglement and superior calibration of long-term predictions.
Exploring useful applications of our disentangled representation is also an interesting direction for future research. Accurate identification of independent kinematics and interaction effects could enable useful downstream funetionalities. For instance, one can perform algorithmic recourse [45] by counterfactual interventions at the object or interaction level.
Modeling limitations The capacity of our GP is limited by the choice of the kernel, e.g., the RBF kernel assumes that the dynamics are stationary. While our model formulation can be combined with arbitrary kernel functions and is possible to increase the kernel expressiveness, e.g by building com- posite kernels or coming up with hand-crafted features, these approaches are often time-consuming and lead highly parameterized kernels that are difficult to learn.
Approximation errors Our posterior inference scheme is inaccurate due to our variational frame- work, approximation errors that accumulate in time during future prediction, and numerical errors caused by numerical integration of ideally continuous dynamics and its solution.
Broader impact In our work, we propose a methodological contribution which is blind to specific data distributions. Its potential and unforeseeable side-effects in fairness-sensitive or safety-critical applications need to be investigated in a dedicated study.
Acknowledgements
The Bosch group is carbon neutral. Administration, manufacturing and research actvities do not longer leave carbon footprint. This also includes GPU clusters which the experiments have been performed. Cagatay Yildiz is funded by the Deutsche Forschungsgemeinscaatt (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-Number 2064/1 Project number 390727645. We would like to thank Jakob Lindinger and Michael Herman for discussions and proofreading.
10

---

References
[1] Nicholas Watters, Daniel Zoran, Theophane Weber, Peter Battaglia, Razvan Pascanu, and Andrea Tacchetti. Visual interaction networks: Learning physics simulator from video. In Advances in Neural Information Processing Systems, 2017.
[2] Thomas Kipf, Ethan Fetaya, Kuan-Chich Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In International Conference on Machine Learning, 2018.
[3] Louis-Pascal Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. In International Conference on Machine Learning, 2020.
[4] Frederik Diehl, Thomas Brunner, Michael Truong Le, and Alois Knoll. Graph neural networks for modelling traffic participant interaction. In IEEE Intelligent Vehicles Symposium, 2019.
[5] Alvaro Sanchez-Gonzallez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia, Learning to simulate complex physics with graph networks. In International Conference on Machine Learning. 2020.
[6] Edward Schmerling. Karen Leung, Wolf Vollprecht, and Marco Pavone. Multimodal probabilis- tic model-basee planning for human-robot interaction. In IEEE International Conference Robotics and Automation, 2018.
[7] Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. In International Conference Learning Representations, 2017.
[8] Jack M Wang, David J Fleet, and Aaron Hertzmann. Gaussian process dynamical models. In Advances Neural Information Processing Systems, 2005.
(9] Roger Frigola and Carl Edward Rasmussen. Integrated pre-processing for Bayesian nonlinear system identification with Gaussian processes. In IEEE Conference Decision and Control, 2013.
[10] César Lincoln C. Mattos, Zhenwen Dai, Andreas Damianou, Jeremy Forth, Guilherme A. Barreto, and Neil D. Lawrence. Recurrent Gaussian processes. arXiv preprinl 1Xii::5111.06444: 2016.
[11] Pashupati Hegde, Çağatay Yıldz, Harri Lähdesmäki, Samuel Kaski, and Markus Heinonen. Bayesian inference of ODEs with Gaussian processes. arXiv preprint arXiv:2106. 10905, 2021.
[12] Katharina Ensinger, Friedrich Solowjow, Sebastian Ziesche, Michael Tiemann, and Sebastian Trimpe. ftructur--peseevvii Gaussian process dynamics. arXiv preprint 102.v:2102.0166 2021.
[13] Arno Solin, Ella Tamir, and Prakhar Verma. Scalable inference in SDEs by direct matching of the Fokker-PlackkKoologgon equation. In Advances in Neural Information Processing Systems, 2021.
[14] James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisen- roth. Efficiently sampling functions from Gaussian process posteriors. International Confer- ence Machine Learning, 2020.
[15] Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi. and Harri Lähdesmäki. Learning unknown ODE models with Gaussian processes. In International Conference Machine Learning, 2018.
[16] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, 2018.
[17] Yulia Rubanova, Ricky T. Q Chen, and David K Duvenaud. Latent ordinary differential equations for srreuulrrl-samppldd time series. In Advances in Neural Information Processing Systems, 2019.

---

[18] Patrick Kidger, James Morrill, James Foster, and Terry J. Lyons. Neural controlled differential equations for irregular time series. In Advances in Neural Information Processing Systems, 2020.
[19] Cagatay Yildiz, Markus Heinonen, and Harri Lähdesmäki. Continuous-timg model-based reinforcement learning. In International Conference on Machine Learning, 2021.
[20] Yuan Yin, Vincent Le Guen, Jérémie Dona, Emmanuel de Bézenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting, In International Conference on Learning Representations, 2021.
[21] Manuel Haußmann, Sebastian Gerwinn, Andreas Look, Barbara Rakitsch, and Melih Kandemir. Learning partially known stochastic dynamics with empirical PAC Bayes. In International Conference Artificial Intelligence and Statistics, 2021.
[22] Zhaozhi Qian, William R Zame, Mihaela van der Schaar, Lucas M Fleuren, and Paul Elbers. Integrating expert ODEs into neural ODEs: Pharmacology and disease progression. In Advances in Neural Information Processing Systems, 2021.
[23] Carl E. Rasmussen and Christopher K.I. Williams. Gaussian processes for machine learning. MIT Press, 2006.
[24] Joaquin Settfonnrr-Candell and Carl Edward Rasmussen. A unifying view of sparse approxi- mate Gaussian process regression. Journal of Machine Learning Research, 2005.
[25] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo innutt.. In Advances in Neural Information Processing Systems, 2006.
[26] Sainbayar Sukhbaatar, Rob Fergus, et al. Learning multiagent communication with backpropa- gation. In Advances in Neural Information Processing Systems, 2016.
[27] Roger Frigola. Bayesian time series learning with Gaussian processes. PhD thesis, University of Cambridge, 2015.
[28] Cagatay Yildiz, Markus Heinonen, and Harri Lahdesmaki. ODE2VAE: Deep generative second order ODEs with Bayesian neural networks. In Advances Neural Information Processing Systems, 2019.
[29] Moritz Dichl and Sébastien Gros. Numerical Optimal Control. University of Freiburg, 2017.
[30] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 2017.
[31] James Hensman, Nicolo Fusi, and Neil Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, 2013.
[32] Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Marc Deisenroth. Deep Gaussian processes with importance-eeighted variational inference. In International Conference Machine Learning, 2019.
[33] Michalis Titsias and Miguel Lázaro-Gredille conjugate inference. In International Conference on Machine Learning, 2014.
Doubly stochastic variational Bayes for non-
[34] Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. In AAAI Conference Artificial Intelligence, 2020.
[35] Zijie Huang, Yizhou Sun, and Wei Wang. Learning continuous system dynamics from irregularl--smpple partial observations. Advances in Neural Information Processing Systems, 2020.
[36] Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-time models of latent stochastic dynamical systems, International Conference on Machine Learning, pages 1726-1734 PMLR, 2019.
12

---

[37] Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. GRU-ODE-Bai DaDeessys: Con- tinuous modeling of sporadicall--obseve time series. In Advances Neural Information Processing Systems, 2019.
[38] Raj Dandekar, Karen Chung, Vaibhav Dixit, Mohamed Tarek, Aslan Garcia- Valadez, Kr- ishna Vishal Vemula, and Chris Rackauckas. Bayesian neural ordinary differential equations. arXiv preprint arXiv:2012.07244, 2020.
[39] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In International Conference on Machine Learning, 2017.
[40] Michael Poli, Stefano Massaroli, Clayton M Rabideau, Junyoung Park, Atsushi Yamashita, Ha- jime Asama, and Jinkyoo Park. Continuous-depth neural models for dynamic graph prediction. arXiv preprint arXiv:2 arXiv:2106.1/58/ 2021.
[41] Kai Xu, Akash Srivastava, Dan Gutfreund, Felix Sosa, Tomer Ullman, Josh Tenenbaum, and Charles Sutton. A bayesian-symbolie approach to reasoning and learning in intuitive physics. Advances in Neural Information Processing Systems, 34:2478-2490, 2021.
[42] Juntang Zhuang, Nicha Dvomek, Xiaoxiao Li, Sckhar Tatikonda, Xenophon Papademetris, and James Duncan. Adaptive checkpoint adjoint method for gradient estimation in neural ODE, In International Conference on Machine Learning, 2020.
[43] 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 Vetito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems, 2019.
[44] Ilya Sutskever, Geoffrey Hinton, and Graham Taylor. The recurrent temporal restricted Boltz- mann machine. Advances in Neural Information Processing Systems, 2008,
[45] A.H. Karimi, J.v. Kügelgen, B. Schölkopf, and Valera. Algorithmic recourse under imperfect causal knowledge. Advances in Neural Information Processing Systems, 2020.
13


Reviews:
[{"review_id": "LVz3X0IjQfk", "reviewer": "Reviewer_RykX", "summary": "This manuscript proposes a method to infer the dynamics of interacting objects while incorporating uncertainty estimates. They use a gaussian process based ODE method, while modeling the dynamics for a set of interacting variables as a combination of independent and interacting terms. The authors also consider the case where they have partial observations of the system. They show the utility of their method on the bouncing balls and charges datasets.", "questions": "Please include figures such as in Figure 2 for the comparison between GPODE and I-GPODE. It is unclear why the ELL is much worse using GPODEs.\n\nWhile the test in Section 5.2.3 is well performed, it would be good to list how I-GPODE does on three interacting objects as well as on the individual objects.\n\nIt would have been good to show the results of this method in a real-world scenario.", "limitations": "The authors have a limited discussion of the limitations and potential negative societal impact of their work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is quite clearly written. The authors detail out the method and the relationship to other methods succinctly. The method they develop is an extension to the GP-ODE method, with additional constraints to take into account the interactions and independent variables. As such, it has limited originality. This also limits the significance of the paper.\n\nWhile the authors show the advantages of using their method as opposed to the GP-ODE (lower MSE and better ELL) on both datasets, it is unclear why the I-GPODE would have a better ELL than the GPODE, since the main change that the authors are proposing is an additional structure. Moreover, their method does not seem to be as accurate as the GDE (I-NODE) method in the MSE (Tables 3 and 5), which we also see in Figure 2. This dampens the enthusiasm for the paper.\n\nThe test performed in Section 5.2.3 is a good one - they first learn the dynamics of three interacting objects, then simulate just one object, and see that competing methods do not learn the individual dynamics as well. This may be because the authors are explicitly learning the individual as opposed to the interactive dynamics.", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "67oOdUEyjcj", "reviewer": "Reviewer_gdCd", "summary": "The paper takes a Bayesian approach to model dynamical systems with Gaussian processes and combines it with a message passing modeling of the dynamics.", "questions": "Can you show that the functions you learn $f_s,f_b$ learn the write functions in the sense that if you disable the interactions then $f_b$ alone would return the correct dynamics?", "limitations": "NA", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The presentation is pretty clear.\n- The method makes sense as it combines two well established approaches.\n\nWeaknesses:\n- Limited novelty. It is a relatively small modification of GPODE.\n- Limited significance. I-NODE achieves better results so main gain is calibration which is a small gain.\n- Results are presented in a wasteful, table 2 and 3 should be joint together for better clarity.\n- The claim in the first line in the abstract that \"for the first time uncertainty-aware modeling of continuous-time dynamics of interacting objects\" is false. E.g. Xu et al \"A Bayesian-Symbolic Approach to Reasoning and Learning in Intuitive Physics\".\n\nSmall remark: you have $p(Y_{1:N}|\\cdot)$ in eq. 6 which isn't clear", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "63ZNtkH4YA7", "reviewer": "Reviewer_F5hb", "summary": "The authors propose a novel uncertainty-aware modeling of continuous-time dynamics of interacting objects. They model the trajectory of interacting systems using latent GP-ODEs by incorporating the interaction term explicitly into the differential function, and use inducing points to scale their variational inference. They show on several datasets that their method outperform other state-of-the-art approaches, that either don't model interactions between the systems explicitly or aren't uncertainty-aware.", "questions": "I really enjoyed the paper and the proposed model, and have only a few minor points.\n\nIn tables 4 and 5, why is I-GPODE-L performing so poorly (I-GPODE-L-S > I-GPODE >> I-GPODE-L)? \nThe improvement of I-GPODE-L-S over I-GPODE is discussed but I would like to have some explanation on why I-GPODE-L is performing worse than I-GPODE. \n\nIn Figure 2 and Appendix Figure 5, why do the shaded lines not follow the inferred trajectory? \nThe links to citation and figures are not clickable, I am not sure if this was intentional or not. On line 123, you also forget to include the section number in \"Finally, please see Section for an investigation\"\n\nEquation 6 with generalization to multiple sequences in the appendix\nIn the Appendix table, I would recommend highlighting the best result in bold, as is done in the main paper. The tables in the appendix are quite hard to read.\n\nI would recommend citing [Learning interpretable continuous-time models of latent stochastic dynamical systems, Duncker et al.., ICML 2019], that is the first paper proposing GP-SDE. This work was prior to [Scalable Inference in SDEs by Direct Matching of the Fokker–Planck–Kolmogorov Equation, Solin et al., NeurIPS 2021] and differed in their Inference algorithm. ", "limitations": "I don't see any potential negative societal impact of their work. ", "rating": 9, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This is a great paper, as it is the first proposed method for uncertainty-aware modeling of continuous-time dynamics of interacting objects. This is an important contribution to the field of latent variable modeling. The paper is very clear, well-written, and results support the claims made in the paper.", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "2f9RWgV0zM", "reviewer": "Reviewer_QEA6", "summary": "The authors introduce a GP interaction model (I-GPODE) which is capable of modelling the continuous time dynamics of a set of interacting objects, whilst providing a robust quantification of uncertainty using latent GP-ODEs. The kinematics of each individual object and the interactions between all objects present are modelled separately; this decomposition allows the complex continuous-time dynamics to be handled by the model. The work leverages recent advances in neural ODEs and efficient sampling from GPs using pathwise updates. Results presented are better or comparable than competing methods in terms of MSE, but considerably better than all competing methods in terms of ELL.", "questions": "As mentioned in the Weaknesses section above, how are the performance metrics and computational performance affected by an increasing number of objects? What is the upper limit in terms of what is feasible to model? I think a relatively brief set of results considering a scenario where $N_a \\gg 5$ would make a nice addition to the results section. ", "limitations": "There is some discussion of the limitations of the work (e.g. overfitting occurring if the dimensionality of the ODE state is larger than necessary), however I think this could have been expanded a little. There is no discussion of negative societal impact, although I don’t feel that this is an issue for the subject matter in question.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": " ### Strengths\nThe I-GPODE introduced in this work is a novel contribution, addressing the problem of *uncertainty-aware* modelling of interacting object dynamics. The authors do benchmark against competing models which can also be applied to this problem setting, and some of these approaches do outperform the I-GPODE in terms of MSE in some cases, however the ELLs presented throughout the results section make a very compelling case that the I-GPODE significantly improves upon existing work in terms of robust uncertainty quantification in this setting.\n\nThe introduction, background and related work sections give a very good overview of where this work sits in terms of the wider research area; Table 1 also helped lend a lot of context. The generative model itself and the variational inference scheme employed are both laid out in a clear and concise manner; Figure 1 does a good job of quickly outlining the key ideas before you even get to the main contributions in Section 3.\n\nThe submission is well written and presented, easy to follow and technically sound. Overall, this is a great piece of work which addresses an open problem, therefore I recommend acceptance.\n\n### Weaknesses\nWhilst there is discussion of how the model complexity scales with the number of evaluation time points, I think some discussion of how the number of objects considered affects both the performance metrics and the computational feasibility of this approach (e.g. computational complexity and/or runtimes) would have been interesting to see. \n", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept