Cheap and Deterministic Inference for Deep State-Space Models of Interacting Dynamical Systems
Authors: Andreas Look
Barbara Rakitsch
Melih Kandemir
Jan Peters
Venue: Accepted by TMLR ()

Abstract: Graph neural networks are often used to model interacting dynamical systems since they
gracefully scale to systems with a varying and high number of agents. While there has been
much progress made for deterministic interacting systems, modeling is much more challenging
for stochastic systems in which one is interested in obtaining a predictive distribution
over future trajectories. Existing methods are either computationally slow since they rely
on Monte Carlo sampling or make simplifying assumptions such that the predictive distribution
is unimodal. In this work, we present a deep state-space model which employs graph
neural networks in order to model the underlying interacting dynamical system. The predictive
distribution is multimodal and has the form of a Gaussian mixture model, where the
moments of the Gaussian components can be computed via deterministic moment matching
rules. Our moment matching scheme can be exploited for sample-free inference leading to
more efficient and stable training compared to Monte Carlo alternatives. Furthermore, we
propose structured approximations to the covariance matrices of the Gaussian components
in order to scale up to systems with many agents. We benchmark our novel framework
on two challenging autonomous driving datasets. Both confirm the benefits of our method
compared to state-of-the-art methods. We further demonstrate the usefulness of our individual
contributions in a carefully designed ablation study and provide a detailed empirical
runtime analysis of our proposed covariance approximations.

Published in Transactions on Machine Learning Research (03/2023)
Cheap and Deterministic Inference for Deep State-Space Models of Interacting Dynamical Systems
Andreas Look Bosch Center for Artificial Intelligence
ondreas.look@bosch.com
Melih Kandemir
kandemarr@imada..du.. dk
University of Southern Denmark
Barbara Rakitsch
barbara. hwwphamwwwkiitscchbooccc....
Bosch Center for Artificial Intelligence
Jan Peters
peters@ias
wttp:/iiiiiinffmmmmmammm informatik.....mmmmm tu tu-darmstadt..em de
Technical University Darmstadt
Reviewed on OpenReview: // net/ forum? id- UggdBy40v5
Abstract
Graph neural networks are often used to model interacting dynamical systems since they gracefully scale to systems with a varying and high number of agents. While there has been much progress made for deterministic interacting systems, modeling is much more challeng- ing for stochastic systems which one interested in obtaining a predictive distribution over future trajectories. Existing methods are either computationally slow since they rely on Monte Carlo sampling or make simplifying assumptions such that the predictive distri- bution is unimodal. In this work, present deep state-space model which employs graph neural networks in order to model the underlying interacting dynamical system. The pre- dictive distribution is multimodal and has the form of a Gaussian mixture model, where the moments of the Gaussian components can be computed via deterministic moment matching rules. Our moment matching scheme can exploited for sample-free free inference, leading to more efficient and stable training compared to Monte Carlo alternatives. Furthermore, we propose structured approximations to the covariance matrices of the Gaussian components in order to scale up to systems with many agents. We benchmark our novel framework on two challenging autonomous driving datasets. Both confirm the benefits of our method compared to sate-of-thhearr methods. We further demonstrate the usefulness of our indi- vidual contributions in carefully designed ablation study and provide a detailed runtime analysis of our proposed covariance approximations. Finally, we empirically demonstrate the generalization ability of our method by evaluating its performance on unseen scenarios.
1 Introduction
Many dynamical systems, such traffic flow Li al. 2018 Yu et al. 2018), fluid dynamics (Ummenhofer et al. 2019) human motion (Jain et al., 2016), involve interactions between agents. Graph Neural Networks powerful tool in these settings, since they allow learning the dynamics of interacting systems from data only. Recent research has made great advances for modeling deterministic complex systems by being able to extrapolate from systems with a small number of agents and short time horizons to systerns with a high number of agents and long time horizons. These methods have been successfully applied a number of physical systems covering fluids, rigid solids and deformable materials (Sanchez- Gonzalzz al. 2020). However, for many real-world applications, predicting single future trajectory for each agent is not. enough, since the stochasticity in the dynamical system
(GNN) (attaglia et al. 2018) have recently emerged

---

Published in Transactions on Machine Learning Research (03/2023)
has significant consequences. For instance, in autonomous driving, the driver's intention (e.g. overtaking, turning, lane changing) is hidden factor that may induce different modes of driving trajectories.
Prediction: Mode 1 Weight: 0.12
Prediction: Mode Weight: 0.74
Prediction: Mode 3 Weight: 0.14
Ground Truth
Beve
Figure We approximate the predictive distribution of a latent Graph Deep Mate-Space Model (GDSSM) as a Gaussian mixture distribution via deterministic moment matching rules. Given historical information in the form of observed Irajectories (dashed lines), our proposed GDSSM architecture predicts the future dynamics while taking interactions between traffic participants into account. We present the true future trajectorics (left-most plot) solid lines. For cach mode and traffic participant, we show the predicted 95% confidence interval of our model (three right-most plots). Our model accounts for interactions, for example, the brown vehicle is only entering the roundabout if the blue vehiele is staying in the roundabout (Mode 1). If the blue vehicle is leaving the roundabout, the entering lane for the brown vehicle is blocked by the blue vehicle (Mode 2, Mode 3) and the brown car has to wait. The ground truth data and the map come from the rounD dataset (Krajewski et al. 2020
In the deep learning literature, multiple architectures have been successfully applied time-series data with two prominent classes being a. recurrent methods that apply fixed transition model repeatedly each time step (e.g. Hochreiter (1991) Chung et al. (2014) ), b) history-based methods that aggregate information from the past using either convolutional filters (e.g. Bai et al. (2018). Oord et al. (2016)) or attention modules (e.g- Vaswani et al. (2017). Li et al. (2019)). Recurrent methods capture the internal state of the system at each time point a latent state Zr. Predicting the output this manner respects the causal order of the dynamical system, i.c. the latent state I, of time point t is needed order to compute the latent state I141 at the next time point + State-space models (e.g. [aärkää (2013) belong to the class of recurrent methods. They are defined by two probability distributions: the transition model, p(#++1 [2:a,, that propagates the latent state forward in time and the emission model, p(Hr(zx)) that maps the latent state into the observational space. Obtaining multistep-ahead predictions for state-space models requires propagating distributions in the latent space. For general transition models, this operation cannot be done in closed-form and existing methods cither use numerical integration schemes (c.g. Sooln et al. [2021), Look et al. (2022)) or Monte Carlo hods (e.g. Krishnan al. (2017), Bayer et al. (2021)). For interacting systems, the latent space grows linearly with the number of agents, requiring a high number of MC samples, which can make these methods prohibitively slow. To the best of our knowledge, mumerical integration schemes have not been explored for multi-agent dynamical systems. In contrast, probabilistic history- based methods directly predict the distribution over future trajectories, mitigating the sampling overhead. However, this approach make the learning problem hard as the model needs to learn the future distribution for multiple time steps ahead. To account for this, complex models (c.g. large neural networks) are often necessary. This can prevent their usage in embedded systems with limited memory capacity.
In this work, we present a novel approach for modeling stochastic dynamical systems with interacting agents that is able to generate expressive multi-modal predictive distributions over future trajectories in runtime and memory efficient manner. We model the unknown stochastic dynamical system Model (DSSM) in which the shared dynamics of all agents are modcled in a joint latent space using GNNs. Our model belongs to the family of recurrent neural networks and we replace the expensive MC operalions during training and testing by introducing a novel deterministic moment matching scheme. Prior work Look al.) (002)) on moment matching for dynamical systems does not consider interacting systems and is restricted to unimodal processes. We overcome the first limitation by applying GNNs in the transition model and the second limitation by placing Gaussian Mixture Model (GMM) over the initial latent state.
Deep State-Space
2

---

Published in Transactions on Machine Learning Research (03/2023)
For each mixture component, we independently apply our moment matching rules in order to arrive multimodal predictive distributions over future trajectories. In autonomous driving, the initial latent state is often estimated from historical information and can provide information about the drivers' intentions (Tang & Salakhutdinov 2019). Conditioned on the initial latent state, the predictive distribution can often be accurately modeled with an unimodal distribution (Cui et al., 2019 Chai et al. 2019). Finally, as there exists a wide varicty of dense traffic scenarios, the high number of agents can result prohibitively large GMM covariance matrices as their size grows quadratically with the number of traffic participants. We address this problem by proposing structured covariance approximations.
We summarize our contribution as follows:
We derive output moments for GNN layers, which makes GNNs applicable moment matching algorithus. This leads to the first deterministie inference scheme for deep state-space models for interacting systems.
We introduce a GMM distribution over the initial latent states that results in multimodal predictive distributions over future trajectories.
We propose structured approximations to the GMM covariance matrices that can reduce the com- putational complexity of our approach from ((M)) O(M)), where M is the number of agents.
In our experiments, we benchmark our proposed model on two challenging autonomous driving datasets. Our results demonstrate that our deterministic model has strong empirical performance compared to state- of-the-art alternatives. We visualize the predictive output distribution for a real-world traffic scenario on roundabout with multiple agents in Fig. [ The futurr distribution is highly multi-modal, as traffic participants can leave the roundabout from several exits. Our model is capable of predicting multiple modes, which we efficiently approximate as GMM, and takes interactions into account using GNNs.
To gain further insights into our model and inference scheme, we carefully examine the impact of the individual contributions of our work in an ablation study. Next, we provide an empirical runtime study of our covariance approximations. Our findings indicate that sparse covariance approximations reduce the computational complexity by a factor of up to 100, which makes them favourable for applications with limited computational resources. We conclude our experiments by studying the generalization capabilities of our model on out-of-distribution data, e.g. traffic environments that have not been observed during training.
2 Background
In this chapter, we first provide background on (deep) state-space models for single-agent systems. We proceed by giving a small recap on moment matching rules that allow us propagate the latent dynamics forward in time in a deterministic manner. Next, we review graph neural networks for interaction modeling. State-space models and graph neural networks form the basis for new model for stochastic dynamical interacting systems, which we introduce in Sec. 3 The deterministic moment matching rules build the foundation of our training and testing routines that we describe in Sec.
2.1 Deep State-Space Models
State-Space Models (SSM) are a model class for dynamical systems (e.g, Schön al. (0011), Särkkä (2013)
that assume that each Dy Dy.dimnnsion observed variable 3e E R. is emitted by latent D, dimeennsion latent variable E D=, The latents are coupled via first-order Markovian dynamics, e.g. the state at time point only depends on the state of the previous time point æ;-1. Typically the observed state yt-1 does not contain all necessary information in order lo reliably predict the next observed state yt- Consider the case of traffic forecasting in which the observed state Ye contains the position of a vehicle but not its velocity or acceleration data. We can then supply the latent state v, with the missing information in order to allow for accurate forecasts about the next time point. Consequently, SSMs are flexible model class that allows to make reliable forecasts about complex systems.
3

---

Published in Transactions on Machine Learning Research (03/2023)
A Deep State-Space Model (DSSM) is non-linear SSM in which the transition model, that maps the latent state from the current to the next time point, and the emission model, that maps the latent state to the outputs, are realized by neural networks. They come in handy for applications in which the true underlying dynamics have an unknown non-lincar functional form and must be estimated from data. Assuming additive Gaussian noise (e.g. Krishnan al. (2017)) their generative model can be written down as follows
zo ~ p()o())
(1)
z N (1ol ( diag (I)),
1 1, ..... T t 1 T
(2)
1/2 N ((ulg(rs), (I)))))),
(3)
where I RD is the context variable that encodes auxiliary information, such as historical or relational information. The mean update [xr, RD x R + R2 R= governing the deterministic component of the transition model, is parameterized by neural net with arbitrary architecture. Similarly, the variance update I) D= x R2 - T is parameterized by another neural net, which models the stochasticity of the system. Both f and L are neural networks, which are parameterized by Oy and Oz. The emission model follows a Gaussian distribution with mean g(zx) R. + RD, and variance [(xx)) RD + R where g and I are both neural networks with arbitrary architecture and paramcters fg and 0r.
Assuming additive Gaussian noise allows us to interpret the transition model [Eq. (2) as discretized neural stochastic differential equation (Tzen & Raginsky 2019 Look et al. 2022). While we do not pursue this line of work any further, we note that this connection allows for saraight-forwan extensions to irregularly sampled time series. Finally, there exists work that couples state-space models with recurrent neural networks (Chung
et al. 2015 Fraccaro et al. 2016) whose gating mechanism can help learning long-term effects.
2.2
Transition Kernel
In this section, we provide background on how to compute the t-step transition kernel, p( I) with t > 1, which allows us to propagate the latent state forward in time for general state-space models. It is defined by the following recurrence
(10e = plare/te D)p(s I)dat-1;
(4)
where p(zo|Z) follows Eq. (L.
It is worth noting that Eq. cannot be computed in closed-form since the
distribution p(# No, I) has to be propagated through the non-lincar transition model p(æ; [ [Eq.
2.2.1 Overview
Various approximations to the transition kernel Eq. (41 have been proposed that can be roughly split into
two groups: (a) MC sampling based approaches (Brandt & Santa-Clara] 2002 Pedersen 1995 Elerian et all
2001) and (b) deterministic approximations based on assumed densities (Särkkä et al. 2015 While MC
based approaches can, in the limit of infinitely many samples, approximate arbitrarily complex distribu- tions, they are often slow in practice and their convergence is difficult to assess. In contrast, deterministic approaches often build on the assumption that the t-step transition kernel can be approximated by Gaussian distribution. This assumption can be justified if the transition model can be locally linearly approximated and the observations are sufficiently densely sampled. The transition kernel then computed in an iterative manner by applying a Gaussian approximation at each time step and propagating the moments along the time direction using a numerical integration scheme (Särkkä & Sarmavuori 2013 Särkkä et al., 2015).
Prior work in the context of neural SDEs (Look et al. 2022 proposes to first discretize the differential equations and afterwards apply moment matching through the neural network layers as numerical integration scheme. Since the resulting algorithm propagales moments through time and neural network layers, it is termed Bidimensional Moment Matching (BMM). This approach was shown to be superior over standard numerical intogration schemes in terms of compute and accuracy. Since the transition model in SSMs can be interpreted as discretized SDEs (Särkkä al.] 2215),, we build our work on their approach. Moment propagation through neural network layers has also been applied previously in the context of expectation
4

---

Published in Transactions on Machine Learning Research (03/2023)
propagation Hernandez-Lobato & Adams 2015 Ghosh et al.
2016
deterministic variational inference (Wu
et al. 2019 and evidential deep learning
(Haussmann et al.
2020
2.2.2
Bidimensional Moment Matching
In this section, we recapitulate the original Bidimensional Moment Matching (BMM) algorithm of Look
et al. (2022) and its use for propagating the latent state forward in time [Eq. (2). BMM approximates the transition kernel p( I) by combining horizontal moment matching along the time axis with vertical moment matching across the neural network layers.
Horizontal Moment Matching In order to facilitate the computation of the transition kernel, we replace p(zg/z0, I) for all time steps t = l..... with a Gaussian distribution
p0 = / 2 ((az)l
Iddx--1
(5)
with mean (I) and covariance (I). This approximation simplifies the problem to calculating the first two moments of the transition kernel and is assumed to work well if the dynamics can be locally approximated by a linear model, which is the case for many applications. Assuming that the one-step transition kernel p(a+/10-11 I) follows Eq. (2), the mean (T) and covariance Ex( can be computed as a function of prior moments (Z) and 1 (I) (Look al.) 2022)
(e( Mr- (I) +
(6)
[I = (Z) + Co + Co Cov[ + diag (E|[(
where Cov(ør-1, ((X-1, denotes the cross-covariance between the random vectors in the arguments.
Vertical Moment Matching The mean ( and covariance (x) of the transition function, as well as the expected variance update E))( can be computed result of moment propagation through neural network layers. For many common layers, including affine transformations and ReLU activation functions, the corresponding output moments can be either computed in closed-form good approximations are available in the literature (Wu et al.) 2019). In contrast, approximating the cross- covariance Cov [r,, fo(x,, Z)] cannot be achieved using moment matching rules since we cannol. decompose the cross-cov term into layerwise operations. Instead, we resort to Stein's Lemma using
Covirr, f(x;, I)] = Cov Ef (xr,,
(7)
where the expected Jacobian can be approximated as
(Look et al. 2022
E
B
(8)
Above, Jf denotes the Jacobian of layer 1 at time step The expectation of the Jacobian is analytically available or can be closely approximated for common layer types.
2.3 Graph Neural Networks
Graph neural networks (GNNs) have emerged as powerful method for interaction modeling (Battagiia ct al. 2018 Hamilton et al. 2017 Gilmer et al. 2017). Given set of agents and relational information in form of a graph, each agent corresponds to one node in the graph that is equipped with a set of fealures. The relation between the agents is encoded via the edges and information exchange between the agents takes place by sending messages along the edges. By performing multiple rounds of paaaaggg,..sing. information can flow along the graph. This allows for interactions between non-adjacent agents, provided that path between the agents exists.

---

Published in Transactions on Machine Learning Research (03/2023)
More formally, define the structure of our GNN follows. For M agents, GNN receives inputs set of node features (*m) Jm=1: where I E RDDD, and xm E IRR:, and set of edges & = {e JMM
mmm''= which is part of the context variable I E RP=, The edge attribute em.m' has a binary encoding, where emmmn if agent m and agent m' are related. The GNN output is an update of the node features, i.c. GN, with 2 RMD,; and consists of the following two steps that may be repeated multiple times:
For each agent m, receive message Nm E IR D= by aggregrating information from neighboring agents:
N AGG
13) R.D:, =
where (e"""e""."" 1} denotes the set of all neighbours of node m. The aggregation operation AGG is permutation invariant, i.c. it does not change when the ordering of the inputs is swapped and generalizes to a varying number of inpuls. A commonly used aggregation operation that we also apply in our work is the mean function.
2. For each agent m, update the node information:
T UPDATE(2" Nmi I],
(10)
where UPDATE(I" ,NN, RD X R. R + DD is typically implemented by neural network.
A simple form of interacting dynamical system takes the features of each agent at its current position and connects agents that are within pre-defined radius with edges. The GNN operation updates the position and velocity information of each agent by taking information of the adjacent traffic participants into account.
3 Graph Deep State-Space Models
We aim to model stochastic dynamical interactions between agents following complex behavioural patterns, such road trallic interactions. We exlend deep state-space models to interacting systems by proposing Graph Deep State-Space Models (GDSSM), which use graph neural networks in the transition model to capture interactions between agents. We define our probabilistic model in this section and then introduce a novel scheme for efficient and deterministic training and predictions in the subsequent section.
We are interested in modeling the dynamics of M interacting agents with deep state-space models by using coupled latent space. In other words, instcad of using D.-dimensional latent space for cach agent, we assume that the agents share latent space of size MD,. Since (i) the number of agents can vary between scenes, (ii) the transition model should be agnostic to the order of the agents and (iii) it is challenging to parameterize high-dimensional latent spaces, we use GNNs in the transition model.
More formally, we denote the state of agent m at time step t as xt" E Rx and the set of all state variables as x, = {a(") m==1- The dynamics of x, follow Eq. 2 where the mean update f(x,, RRD., RDz MD. and the variance update L(x;, I) RDD, R MDD, implemented with the help of graph neural networks
fla), MM,
a(.l,aMi,2)
f(2nI)) [iaM No II]
L(x
(11)
NNx II)
The agent-specific mean update is denoted by f(a?" NN..nn II P. x R. x P2 RD. and the variance by Mar Nm I) D= x RD= x RDr - R.f Both implement the update function in general graph neural networks [Eq. [10]]. whereas No contains aggregated information of the states from all neighboring agents Eq. A deterministic variant. of our model, e.g. setting T(x,, I) 0, has been successfully used for

---

Published in Transactions on Machine Learning Research (03/2023)
learning surrogate models for complex physical systems (Sanalez-Gonzalez et al., 2020 We further assumm that it is sufficient. to couple the latent dynamics across the agents and keep the emission model [Eq. (3) independent across agents. Note that our transition model consists of a single aggregation and update step which is sufficient if the data is densely sampled such that the information flow between agents is fast compared to the evolution of the state dynamics. However, is also straight- forward to extend the model to multiple message- passsgggpasinggggg steps per time point by stacking multiple GNN layers in the mean and variance update function.
Furthermore, we note that although the transition noise factorizes across agents, correlations between agents emerge since the mean and the variance depend not only on the state of the m-th agent, but also on the states of all neighboring agents. After a aggregation steps, our GNN model accounts for correlations wet agent m and agent m' provided that they are connected by a path that is at most a steps long. In contrast, methods, that only take the state of the m-th agent into account, do not lead to any correlations, while methods, that take the state of all other agents into account, lead a fully correlated covariance matrix after one time step.
In order to complete the probabilistic description of our model, we further specify the distribution of the initial latent state To E RMD, with Gaussian Mirture Model (GMM)
(2) TV (Z)), No ~ NNvoluo (Z), diag( Lo, (Z))),
(12)
where Cat([Ti (Z).... TV (I)) is a categorical distribution with V mixture components. Each component is specified by its weight (I) RRr + R,, mean (I, R + RMD=, and diagonal covariance (Z, RDr + R.D. The weights T1:V form standard V-simplex. We use a GNN, which we refer to as the embedding function h(Z): _D + RV +V+2VMD in order to model the initial state distribution
(I)
h(I)=
Po.1.V (I)
(13)
(T.
We assume that the context variable contains relational information as set of edges as well as historical information for each agent in the form of observed trajectory. In sense, the embedding function acts hereby filter, which learns distribution over the initial latent state from past observations.
In autonomous driving, the initial latent state can be connected to the drivers' intention (Tang & Salakhut- dinov 2019). The context information I is often not sufficient to rule out different hypotheses about the future, e.g. does the car behind us want to overtake in the next five seconds or not. Using a mixture model allows us to incorporate different hypotheses into the model in a principled manner, which will ultimately lead to highly multinodal predictive distributions.
State-space models and graph neural networks have been previously combined for multi-agent trajectory forecasting by Yang et al. (2020). In contrast to our work, the authors (i) use different model definition by applying recurrett neural networks and a non- Gan-saassian density in the transition model and (ii) perform Monte Carlo sampling during inference which can lead to slow convergence. In the next chapter, we show that our model definition allows for more efficient training by performing deterministic moment matching rules. We compare against Monte Carlo alternatives in our experiments.
4 Deterministic Approximations for GDSSMs
In this chapter, we present our novel inference scheme for GDSSMs that enables efficient and determinis- tie training and predictions. In Sec. 4.1 we first give short overview over existing inference techniques and compare it to our scheme which aims at directly maximizing the predictive log-likelihood of future trajectorics. The prodictive log-likelihood cannot be computed in closod-form. We propose an efficient and deterministic approximation to it in Sec. 4.2 that relies on bidimensional moment matching. Our moment matching rules necessitate the computation of output moments and expected Jacobians of graph neural
7

---

Published in Transactions on Machine Learning Research (03/2023)
network layers. We present ouput moments for commonly used layers in Sec. 4.3 As our algorithm approx- imates the output distribution at each time step with a Gaussian mixture distribution over all agents, the resulting covariance matrix for each mixture component can become computationally intractable for a large number of traffic participants. We address this pain point in Sec. 4.4 by proposing sparse approximations to the covariance.
4.1 Parameter Inference
Classical inference methods for state-space models aim at directly maximizing the log-likelihood of the data
T
log. -T. log / p( (I)
II
p(a (2 Typ(y are dd drr,
(14)
where ((1) defined Eq. (2) and ())) Eq. 1 This quantity can only be computed closed-form if the emission and transition model are linear Gaussians. In our case, the transition model is parameterized by graph neural network and the emission model by standard neural network. Both functions are highly non-linear and render an analytical solution to Eq. (14 infeasible.
Therefore, most existing methods apply either a particle filter (Schön et al., 2015), variational inference Krishnan et al.
Bayer et al. 2021) or a combination of both (Naesseth et al. 2018 Maddison et al. 2018 order to approximate the log- liclliiclhood. All of these approaches have common
2017
2017 Le et al. that they require learning a proposal distribution q(Fo,..., FT .B......r7) In particle filtering, the proposal distribution is recursively defined by the importance function q( (T0, -......--,--. - 91, y1...... yt) and is optimal in terms of variance by setting q=za/=o..... ...,, ..., 21 ,l.....,3.. 21) p( ye) (e.g. Doucet et al. (2000). Variational inference aims to find the best approximation to the true posterior within the chosen variational
family by minimizing the Kullback-Leiblee (KL) divergence between the true and approximate posterior (Blei et al. 2017), while variational sequential Monte Carlo seeks to minimize the KL. divergence on an extended sampling space (Le et al. 2018).
However, the proposal distribution is only used as an auxiliary tool during inference. For many prediction tasks (Djuric et al. 2020 Jain et al. 2019 Chai et al. 2019), the quantity of interest is the predictive
likeeiodddood (PLL)
T
.......... PLL(11 .0..........)))) I)
Σ logp(v((Z)
(15)
1=1
T
Σ log / p(ao|T)p(aaaag T)p(ye drddda,
where the transition kernel p(ue|=0, Z) is defined in Eq.
In contrast to the standard training objective, the predictive ligeliiiihooo propagates the latent state forward in time without receiving any feedback from the observations; mimicking the behavior during test time. Since we want to use the same objective during training and test time, we opt for directly maximizing the predictive logelikel during training. The observation that using one-step ahead predictions in training is not sufficient. order obtain reliable multi-step ahead predictions during testing has also been made Bengio et al. (2015) where the authors proposc a schoduled sampling strategy in order gradually switch from single to multi-step ahead predictions during training. Multi-step ahead training has also been successfully applied for maato-temporal forecasting Pal et al. 2021 and model- planning (Hafner et al. 2019).
The PLL aims at optimizing the average marginal log- logelikelihood logr(s]]Z), while the standard objective aims at optimizing the joint likelihood log tog II). We deem the choice between the two objectives application dependent: The PLL most useful for tasks that can be solved by assessing the marginal distributions only. An important and large application class that falls into this category in the context of autonomous driving in which the marginal distributions of neighboring traffic participants at given time
8

---

Published in Transactions on Machine Learning Research (03/2023)
horizon often sufficient in order to control the car (e.g. Herman et al. (2022) As consequence, many papers in the autonomous driving literature report as evaluation metrics the performance of their met at fixed time intervals which matches our PLL objective (e.g. Djuric et al.
(2020); Jain et al. (2019); Chai ct al. (0191).. ). In contrast, optimizing the joint log-likelihodd is preferred for tasks that require sampling realistic looking trajectories, as it is done for instance sentence generation (Vaswani et al. 2017
4.2 Approximating the Predictive Log-Likelihood using Bidimensional Moment Matching
We are interested in the predictive log- likelihoddod PL Vr|Z), which describes the predictive log- likelihood of all traffic participants up to time step In order to calculate it, we need to solve the nested set of integrals given in Eq. (15). The BMM algorithm allows us approximate the transition kernel p(ax|Z) J p(ara (oo [)p(2 Tydo in case that the initial state zo has Gaussian distribution. However, in our model formulation the initial latent state zo follows GMM to allow for multimodality. In order to account for that, we approximate the marginal latent distribution p(a+lZ) as
[ = Σ To To(D)((a,o)
(16)
where each mixture component [T) is approximated with the BMM algorithm p(
N (Z), (Z)). Assuming GMM at the initial state allows to efficiently obtain multimodal predic- tions while being computationally efficient: We obtain multimodal distributions by explicitly marginalizing
over the mixture components and we compute each component efficiently by applying the BMM algorithm. Finally, we approximate the marginal distribution p(ve (Z) as a GMM by another round moment matching
p(ul [I) = / (ingg
(17)
nu (Z), Br,(()).
where (I) and (I) are the mean and covariance of the v-th mixture component at the t-th time step. These two moments are available as
ag.0 (Z) = Elg(
Br,, = Cov[g(Zr,v)) + diag (x:)()
(18)
which is direct outcome of the law of the unconscious statistician. We present the pseudocode for computing p(y.(Z) using our method in Algorithm
Algorithm 2 further shows how we can optimize our model
parameters employing the PLL [Eq. 15
training objective.
In our method, we approximate p(=,, (Z) and p(u (Z) by applying moment matching to each Gaussian mixture component independently. We note that this approximation becomes exact for locally linear transition and emission functions as stated Thm. [
Theorem 1. The marginal distribution P( is analytically computed as
P( = Σ m. (II
for a GDSSM with the below generative model
(It([=](I),.... Ty (Z)),
Zo N(Ho, diag(Eo Co,, (I))),
2. N (xalate11 + 1(1, v, I) ding (I(I,, Z))), 3/ ~ N (velg(t, T),,, diog (T v1,,,22) ())))
1 = 1 T T
D(t, ltt, I) are time t, component u, and contezt depending matrices with where f(t, v. I), Ltt, v g(, v, v) appropriate dimensionalily.

---

Published in Transactions on Machine Learning Research (03/2023)
Proof. The proof is straightforward as the output moments of the transition and emission function are analytically tractable. We provide details in App. A
Algorithm Ridimensional Moment Matching in Latent Space (BMMLS)
Inputs:
Prediction horizon
Context variable
D. Mean update
L(x; Iz (z)
Covariance update
((m):0a)
Mean emission
((
b Covariance emission
h(Z
Embedding function
Outputs: Approximate marginal distribution P(ar|Z)
Ho, (Z), Si (T)) 0,)
GMM at initial step, Eq.
13
mixture component E (1..... V)
for time step E (0........................ ... .
do
Horizontal Moment Matching
Ms+1. + (x(2) + H[
Bq
E++1,> Ex(2) Cov( Covlzs, ((s,, Covla,, f 2:0,11 diag (BI,(H,, end for
Eq.
a + Blg(zr,.:0p))
Eq.
18
Br. (2) Covlg(ar.wi 0 + diag (B()());0r)))
Eq.
end for return Σ (Z),
Algorithm 2 Deterministic Training of a GDSSM
Inputa:
Dataset
1:0,)
Mean update
Lx0,)
Covariance update
g(m;;
Mran emission
[(#::0r)
b. Covariance emission
M(I; 9x) Outputs: Optimized weights 9 (0y, B2, 0g. Ar. 9n)
Emhedding function
while not converged do
I. #1.T T for MI. do
b Sample and Bur from dataset.
P(m|IZ) = (s I,
b Approximate marginal distribution p(M(IZ) vla algorithm
PLL(n)....
D
Caleulate predictive log-likellbood 13
8 - 9+ apt.yyr13)
Update weights 8
08
end while return 8
For general non-linear transition and emission functions, the quality of this scheme depends on two factors: (a) how well the transition and emission function can be approximated in locally linear fashion and (b) if the covariance matrix Er, is small enough over the time horizon t. We observe in our experiments that the approximation works well and further illustrate its behavior that Gaussian mixture models have also been successfully used for filtering problems (Alspach & Sorenson, 1972 where, similar as in our work, moment matching is performed for each component individually in the predict step.
small toy dataset in Fig. 2 Finally, note
4.3
Output Moments of Graph Neural Network Layers
In order to use the the BMM framework, we need to be able to calculate the first two output moments as well as the expected Jacobian of graph neural network layers. In the following, we derive the analytic expression of the output moments and expected Jacobian for the common graph neural net layers: (i) nodc-wise affine transformation, and (ii) mean aggregation. Output moments for the ReLU activation provided in Wu et al. 2019). For completeness, we also provide the output moments for the ReLU activation in App. B
Let 2, E Pa. be the node features at layer of node m at time step and xt {ab Jm=1 the set of
all node features with of R.D... For the sake of brevity, we have omitted here the index of the mixture
10

---

Published in Transactions on Machine Learning Research (03/2023)
1.00
0.34
0.33
0.34
0.33
0.00
Vallue
Valuu
Valum
Time secords
Time seconds
Time seconds
Time seconds
one-dimensional, multi-modal toy problem. Solid black lines represent the ground truth dynamics with three different modes of equal probability. For our model, we report for each predicted mode the mean, 95% confidence interval and mixture weight. We set the dimension the latent space to D2 = and vary (from left. to right) the number of modes V in the model from to 4. For V < 3, the number of modes is set too small, and we observe mode covering behavior. For V 3, our model can recover the ground truth dynamics. If the number of components is too high (V=4), the mixture weights of redundant components are set to
Figure 2: Predictions of our model after training
component
We denote mean and covariance of a graph with M nodes layer and time step as
B
Corlet a
Cov(ae 1.1 M
B[ay]
Cov[ur]=
; J.M xz
(19)
Coolll
Covlat M L.M.
J.m. J.m' Dr. A typical GNN architecture (sce Sec. 2.3) consists where B[a?"""] E PRa and Cov(a, Sht of alternating aggregation steps, in which information of all neighbors is collected, and update steps, in which the features of the node are updated. For the aggregation step, we derive the output moments for the commonly used mean aggregation operation in Sec. 4.3.3 For the update step, we assume that the neural network is built as a sequence of affine transformations and nonlincaritics. The output moments of nonlinear activations are applied independently across agents. As consequence, their rules do not change when used in the GNN setting and we can use the derivations from Wu et al. (2019). Hence it remains open to derive the output moments for affine transformations, as they used in GNN context, which we tackle in Sec. 4.3.2 The mean aggregation operation and the node-wise affine operation used in the update step are special cases of a standard alfine layer, which we review in Sec. 4.3.
4.3.1 Standard Affine Transformation
Suppose apply affine transformation node at layer / with state 2/""
1+1. = A + + Tt
(20)
with weight matrix wi and bias B'. The output moments are analytically tractable as
Bla 1+1 1+ W'EL
Cov [rt 1+1, [.]] W Cov [of
The expected Jacobian of the affine transformation reads as J W!.
4.3.2 Node-Wirs Affine Transformation
The node-wise affine transformation applics each node layer with state z I. the same transformation simultaneously with weight matrix W and bias b/ The node-wise affine transformation be interpreted
11

---

Published in Transactions on Machine Learning Research (03/2023)
standard affine transformation acting on the set of all nodes z/ as
W + 8' W 42 + b'
W
0
0
b
W
0
+
= (LM & 1 x + (1M & b',
(21)
W W
J.M
+ M/
W
L.M
W
where TM is the identity matrix with shape M x M, 1Ms is vector of ones with shape M x 1, and & is the Kronecker product. The output moments of node- wise affine transformation are
Elat
Covlat ( W W"Coolz'(
Similarly as for the standard affine transformation, the expected Jacobian of node-wise affine transformation is analytically available as J W
4.3.3 Mean Aggregation
A commonly used aggregation operation is the mean aggregator, which calculates the message L.N. m at time step at layer as
to node
J.Nm
2t W
(22)
m' m'
Let L.W be the set of all messages, i.c. written a linear transformation
2t
INN
{zh.Vm
I. Im-
The mcan aggregation can be equivalently
I.N (A & 1D..)
(23)
Above, A E RM denotes the row normalized adjacency matrix, which summarizes the edge information & in matrix format and ID.,, denotes the identity matrix with dimensionality Da.1 Dz,1- The Kronecker product expands the adjacency matrix accordingly to the Dg.-dimensional node features. Hence, the mean aggregation corresponds to linear transformation with a weight matrix consisting of M x M blocks, where each block is diagonal matrix of shape Dæ x Da- Its moments are analytically available as
Ela IN AEE=']
Covlat I.V
(24)
The expected Jacobian is available as J -/=AA = Å.
4.4 Sparse Covariance Approximation
For settings with large number of agents M or with high-dimensional state 2.77 the application of the BMM algorithm can become computationally expensive. In the following, we review the computational complexity of the BMM algorithm. We assume that the GNN model consists of a mean aggregation step followed by multiple node-wise affine transformations and nonlinearities, which is the same architecture that we employ later on in our experiments. The mean aggregation is done for each of the DN latent states independently, and we denote the maximum hidden layer width with H.
Computing the nonlinearities is cheap as the operation acts elementwise and their effect on the runtime can be neglected during this analysis. The other two operations (see Sec. 4.3.2 and Sec. 4.3.3) can be described by affine operations for which the weight matrices are heavily structured; the mean aggregation step corresponds to weight matrix consisting of M x M diagonal blocks, the node-wise affine transformation to block- diagonal weight matrix with M blocks of shape H x H. Propagation of the full covariance mairix through a neural network (forward cost) has the computational complexity of 2 H2 D² +M² HDD++M' D2), where the first term is due to the cost of the H H-dimensional node-wise affine transformations in the hidden layers, the second and third term are due to the cost of the H D-dimensional node-wise affine transformation after the aggregation operation, and the fourth term is due to the aggregation operation.
12

---

Published in Transactions on Machine Learning Research (03/2023)
Diagonal Covariance
A
A
B
Off- Off-Diagonal Covariance
C
Diagonal
D
Cross-Covariance
Off- Off-Diagonal
Cross- coos--Cvvreennce
A B D
Figure 3: Groupings of the covariance matrix for graph structured data. Consider the example graph with M agents, which is depicted in the left panel. Each agent has dimensionality Dz. In the middle panel, we show the covariance matrix between all agents, which consists of 5 x 5 blocks, where each block has dimensionality Dæ x Dr.
The computational cost of the expected Jacobian is in the following. Let the expected Jacobian of the aggregation operation be B[J and the product of the expected Jacobians of the subsequent neural net layers [[J]"*] 11 E[J]]. The expected Jacobian of the neural net layers E[Jret] is block-diagonal with M blocks of shape D3 D3 and its computation takes O(MH + M 2H D2) time. Multiplying B[JØ] with B[Jret] results in a fully populated matrix where each entry can be computed by a single dot product, due to the structure of its factors, and its computation contributes with O(M?D)) to the runtime.
DZ) and give its derivation
Consoquently, the total cost is dominated by the forward cost, using the full covariance matrix Since the computational cost quickly becomes intractable due the cubie dependence with respect to the number of agents in the forward pass, we next propose different sparse approximations to the covariance matrix. Independent of the chosen approximation, the cost of the expected Jacobian remains unchanged, as it does not depend on the covariance matrix.
when M3 whenn
Full: Model the full covariance matrix.
Forward cost: O(M + + M
Total cost: H + M2 + MD + M))
Main Diagonal: Keep the diagonal entries in the covariance blocks, which corresponds to the blue line in Fig. 3
Forward cost: 0(MH² + MHD + M²).
Total cost: OMMH3 MH"D, + D²)
Main Blocks: Keep the block-diagonal blocks the covariance matrix, which corresponds to the orange blocks and blue lines Fig.
Forward cost: 3 2 Dr MHD? + M²D)). Total cost: OMM M D. MHD M2D)).
All Diagonals:
Structure the covariance matrix in blocks of shape M M and keep the diagonal entries in each block, which corresponds the blue and red lines in Fig. 3 Forward cost: H + M + M²).).
Total cost: H + M + M + MH MH + M²D2).
Note that setting the off-diagonal blocks to zero, as done in Main Blocks and Main Diagonal, corresponds to independence assumption between the agents and leads to runtime reduct ion from ((M3) to O(M)). In contrast, assuming diagonal structure within each block, as performed in Main Diagonal and All Diagonals, corresponds to an independence assumption between the features and leads to runtime reduction from ((H)) to O(H²) the forward pass.
For reference, taking Monte Carlo approach has the computational complexity O SMHD. +MPD)) where 8 the mumber of Monte Carlo samples.
13

---

Published in Transactions on Machine Learning Research (03/2023)
Velicle 1
Log Covariance
pesmaain 601
Vehicee Vehicle
Vahicke
5
6
Vehicle
Vehicle
Vehicle
(a) Example Scene.
(b) Covariance at second. (c) Covariance at 3 seconds. (d) Covariance at 5 seconds.
We visualize an example scene from the rounD dataset (Krajewski et al. 2020 and the covariance matrix in the latent space at different time steps for the case of a unimodal initial distribution. In the left plot, dashed lines represent the observed history and solid lines represent the true future trajectories.
Figure 4:
Finally, it is important to note that the covariance matrix only has the same structure the graph after the first time step. Agents that are not connected via an edge can still have non-zero cross-covariance at time step I, provided that they are connected by a path that is at most steps long. In our applications, this leads to non-sparse covariances after a few time steps, since the number of agents is small compared to the time horizon. We present the covariance matrix at three different time steps for an exemplary scene in Fig. For a short prediction horizon of one second, the covariance matrix has an approximately diagonal shape. As the prediction horizon increases, the covariance matrix becomes dominated by its diagonal entries. We remark that the could further increase the information spread across agents by using architecture with multiple aggregation steps within the GNN module if required.
complex and is longer
5 Experiments
We provide experiments two challenging autonomous driving datasets. The first experiment (Sec. 5.1 conducts ablation study, while the second experiment (Sec. 5.2) benchmarks model against state- of-the-art methods. We provide runtime analysis and a benchmark of our proposed sparse covariance approximations in Sec. 5.3 In Sec. 5.4, we analyse the generalization capabilities of our model by bench- marking on novel unseen traffic environments.
We provide details about the training procedure in App. and the architocture for the embedding, tran- sition, and emission functions in App. H Accompanying code is available under https://githubb.com
As evaluation metrics, we use the Root-Mcan-Squaal Error (RMSE) and the predictive negative log- log-likelihood (NLL) at fixed time intervals (1s, 2s, 3s, 4s, 5s). For more details about the evaluation, we refer to App. D
5.1 rounD
The IrounD dataset (Kajeewski et al. 2020 consists vehicle trajectories recorded at different roundabouts in Germany. As the roundabouts involve many interactions among vehicles, we expect the predictive distri- butions multimodal and highly complex. We use the recordings from the roundabout in Neuweiler near Aachen for training and testing purposes. The dataset consists of 13,129 tracked objects recorded at 25 Hz in 22 sessions, amounting to total recording time of 6.6 hours. We remove pedestrians, bicycles, well parked vehicles from the dataset, as their influence negligible. After datasct curation, we are left with 12,715 tracked objects. We downsample the recordings by a factor of five and construct a dataset consisting of eight-seeond--long segments with 50 % overlap, resulting in 5405 snippets. We use the first three seconds as the track history, which part the context variable I, and the following five seconds as the prediction horizon. The first recording scssions, corresponding 4,314 snippets, are used for training and validation. The final four recording sessions, corresponding 1, 1,091 snippets, are used for testing.
the vehicle behaviour patterns in the roundabouts is
14

---

Published in Transactions on Machine Learning Research (03/2023)
We build the connection graph by connecting vehicles with an Euclidean distance of less than 30 meters. The dataset has been recorded with drones and we keep the original global coordinate system.
5.1.1 Baselines
We compare our method with multiple baselines models. For each baseline, we remove one key assumption of our model. We cite papers that employ similar ideas appropriate. We did not reimplement these works but performed an ablation study in which we replaced specific components of our model in the interest of fair comparison. Furthermore, we compare against different training strategics as an alternative to maximizing the PLL and analyze the multimodality of our model.
(1) Model:
(L.i) GDSSM: Our model as proposed in Sec.
(1.61) Non Recurrent GNN (e.g. Casas et al. 2020 Herman al. (2022) This architecture receives the context variable I, performs one round of message passing, and subsoquently outputs one normal distribution for each of the next five seconds without using a recurrent architecture.
(1. No Latent Noise (e.g. Sanchez-Gonzaler et al. (2020)): ): We remove the noise from the latent dynamics [Eq. 22 while keeping the emission model unchanged. The uncertainly can no longer be propagated forward in time as the emission model acts independently for each time point.
(i.i0) Linearily (e.g- et al. (2020)): ) We remove all non-linearities from the latent dynamics. Note that keep the non-linee non-linear emission model in order map the dynamics into latent space in which the system can be linearly approximated.
(i.v) No Interactions (e.g, Krishnan al. (2017)): Our model with a diagonal adjacency matrix, i.e. we remove all edges from the graphs. This model neglects interactions between traffic participants.
(it) Modes:
(ú.i) V Modes: We vary the number of mixture components V in the GMM prior [Eq. 12 in order to lest the ellect of multimodality.
(ii.i) oo Modes: This alternative does not follow an assumed density approach, i.c. the marginal latent distribution p[zx[Z) is not approximaled as a GMM with a bounded number of modes. Instead p(az[Z) is approximated by simulating a large number of trajectories, where each trajectory can follow different mode (Brandt Santa-Clara 2002). We set the number of particles to 100.
(ii Objectives:
(ii..) PLL/Det.: We train the model on the predictive logelikel [Eq.
15
and use our deterministic
approximations for GDSSM as proposed in Sec. 4.2 and Sec. 4.3
(M..)) PLL/MO: We train the model on the predictive log-likelihood Egq 15
and take assumed density GMM. The intractable integrals
approach by approximating the marginal distribution [Eq. 117 are solved via Monte Carlo (MC) integration. One forward pass through our model amounts ap- proximately to the computational cost of 12 Monte Carlo simulations. For a fair comparison, we use during training 16 particles, which more costly than training with our proposed moment propaga- tion algorithm, and test with 100 particles. We use the same amount of particles for all MC based methods.
(il.it) ELBO/MC (e.g. Krishnan et al. 2017) ): The model is trained by maximizing the Evidence Lower Bound (ELBO). We give a description of this loss function App. E. The approximate postcrior filtering distribution that models the latent state as normal distribution. We propagate the
latent state forward in time direction via Monte Carlo sampling.
15

---

Published in Transactions on Machine Learning Research (03/2023)
(t MCO/MC (e.g. Maddison et al. 2017 ): The model is trained by maximizing the Monte Carlo Objective (MCO), which combines particle filters with variational inference. We give description of this loss function in App. E The proposal distribution is a filtering distribution and we propagate the particles in time direction via Monte Carlo sampling.
5.1.2 Results
We provide benchmark results of all methods in Tab. T First, we compare our deterministic training and testing scheme (GDSSM PLL/Det.) against its Monte Carlo alternative (GDSSM PLL/MC). Though Monte Carlo based training is more costly than training with BMM, the Monte Carlo results are significantly outperformed by our method. Our results indicate that deterministic approach leads to more effective approximations compared to Monte Carlo sampling despite the approximation error we obtain by our de- terministic moment malching scheme. One potential explanation for our finding is that the Monte Carlo approaches suffer under a high variance since the latent space for multi-agent space grows linearly with the number of agents.
When changing the training objective from the PLL to ELBO/MC or MCO/MC, we observe that the performance is comparable PLL/Me for the prediction horizon of 1 second. For longer prediction horizons, MCO/MC degrade quicker compared to
the performance of the models that are trained on ELBO/MC the model trained on PLL. We believe that this behavior can be explained by the mismatch wetween training and testing objectives. When training on MCO/MC or ELBO/MC, the model obtains feedback from the observations via the proposal distribution, while during testing it has to produce multi-step ahead predictions without receiving any feedback from the observations.
Next, we study if our method can capture multimodality by increasing the number of components in the GMM prior. It worth noting that GMMs can approximate arbitrarily complex distributions when the number of components is chosen high enough. In our experiments, we observe that increase of components significantly decreases the NLL., making the GMM prior a vital ingredient of our method and suggesting that the true predictive distribution is highly multi-modal. If the mumber of components is chosen too small, our model adjusts its uncertainty predictions accordingly. For example, we observe in Fig, Tthe uncertainly of the orange agent to increase as the vehicle is close to the exit of the roundabout. The high predictive uncertainty can be explained by two potential future outcomes: the orange vehicle can leave the roundabout or stay inside. Consequently, our model learns to compensate if the number of components is picked too low, which in turn allows us to trade accuracy for computational runtime. When using tailored implementations, one can easily scale up to a larger number of components, since their computations can be parallelized without any hurdles. We further note that the RMSE is less affected by the choice of the number of modes and stays within two standard errors. Since the RMSE value measures the error between the true value and the expected value of the predictive distribution, this result suggests that the expected value can already be modeled well by predictive distributions with very few modes (actually, one mode seems to suffice for this). We further provide the minRMSE values in App F that decrease with increasing number of modes, indicating that the different modes correspond to diverse set of plausible trajectories.
Next, we compare our model to simpler alternative in which the dynamics are assumed to be linear, which allows calculating the moments of the transition model exactly and in closed form (Särkkä & Solin 2019 In contrast, our approach, GDSSM, approximates the non-linear dynamics in a local linear way by using deterministic moment matching results. We find that our approach achieves lower RMSE and NLL, which can most likely be altributed to the higher modeling flexibility of our proposed model class.
Our ablation study further shows that discarding latent noise and removing interactions between traffic participants results in higher RMSE and NLL. Compared to modeling the dynamics in non-recurrent manner, our model achieves a similar RMSE and outperforms its competitor in terms of NLL.
5.2
NGSIM
The Next Generation Simulation (NGSIM) dataset (Halkias & Colyar 2007) consists of vehicle trajectories recorded at 10 Hz. at two different. highways, US-101 and 1-80, in the United States. The dataset is commonly
16

---

Published in Transactions on Machine Learning Research (03/2023)
Table 1: RounD results. We provide RMSE and NLL (mean + standard error over 10 runs).
Now Recur. GNN
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
GDSSM
No Lat. Noise
Lineaz Modlo
No Interaction.
Mode PLL,/Det.
Mades ELBO/MC
Modls MCO/MC
1 Modle PLL//Det
Modle PLL/Det.
Mode PLL/MC
Made PLL,/Det
Mades
Moxles
Modes
Pit.DDet.
PLL/Det.
PLL/Det.
PLI/Des
1s
0.72
+ 0.02
1,53 + 0.04
0.95 + 0.34
1.22 = 0.08
0:58 = 0.02
11311 0.00
0.90 = 0.02
0.79 + 0.002
0.76 = 0.02
0.74 0.03
0.73 0.03
RAISE
1.00
3.42-
0.00
2.45 0.00
0.08
2.12 0.03
2.00 0.02
1.87 0.03
1.84
0.04
35
3.54
0.03
5.99
1 0.14
4.55 1 0.18
3.88 0.07
3.88 0.05
3.89
0.06
3.60 0.04
3.36
1 0.03
3.35 0.03
3.35 0.03
3.33 0.04
46
5.26
+ 0.34
1:24 = 0:26
7,82 + 0.27
0.58
01.08
1.13 = 0.09
0.39
0.13
5.37 + 0.00
5.08 d: 0.04
5.15 =: 0.00
5.14 2: 0.08
5,06 0.04
7:20
11:24
0.10
8.83 0:13
7.65 0.06
7:34
7:34
1x
1.30
0.01
2.62 0.04
2.70 + 0.06
2.96 0.08
1.85 0.08
1.67 0.07
2.82 0.03
1.48 + 0.05
1.34 + 0.08
1.36 +
0.05
NL.
2.
3.25
+ 0.02
4.45 + 0.07
4.21 + 0.03
3.80
0.10
3.93 = 0.13
3.37 0.08
4.11 = 0.02
2.91 + 0.00
2.93 = 0.07
2.33 # 0.06
2.79 0.00
4:40
0.00
0.03
0.13
0.19
4.07
0.09
3.67
4.01
5.13
0.02
6.59
1 0.13
6.01
1 0.03
8.17 0.16
6.01 0.22
4.03 0.09
5.01
0.07
4.46 1 0.03
4.52 0.11
4.21 0.15
4.34 0.18
5.
5.71
0.02
6,98
0.14
6.73 + 0.04
7.26
0.20
0,80
0.22
0.01 0.10
5.41 = 0.00
5.05 4: 0.04
4.82 0.24
4,39 0.15
4.27 0.23
used for benchmarking traffic forecasting methods and allows us to compare our method against prior art, We adopt the experimental setup of Deo & Trivedi] (2018) and use both highway scenarios. As provided in the original publication, we employ a local coordinate system that is centered around an ego-vehicle, We split the scenarios into three 15 minute long time spans resembling mild, moderate, and congested traffic conditions and downsample each trajectory by factor of two. The test set consists of a fourth of all trajectories randomly sampled from both locations. As in the rounD experiment, we split each trajectory into eight-second-looig segments, where the first three seconds are used as the track history and the following five seconds as the prediction horizon.
Similarly as in Diehl et al. (2019): Lenz et al.] (2017); Wheeler & Kochenderfer (2016). we introduce a connection graph based on the lane position of each vehicle. Each vehicle has at most six connections to other vehicles, which are the nearest vehicles in /ront/behind the cone/left/righh lane. The vehicles
the outermost lanes have a maximum of four neighbours, as there are no neighbours to the left or right.
5.2.1 Baselines
We compare our approach with the following methods:
(1) Constant Velocity [Mereat et al., 2019): This method uses a linear state-space model together with a Kalman filter in order to make predictions. It does not take interactions between agents into account.
(ii) Convolational Social (CS)-LSTM (Deo & Trivedi 2018): Interactions between vehicles are modeled by introducing a grid and applying a convolutional layer on top. Dynamics are modeled by a deterministic LSTM, which predicts the mean and the variance of a normal distribution at each time step.
(wv) Tentooll (ST)-LSTM (Chen et al. 2020): Interactions are modeled with a apatio-temporal graph structure, i.c. graph with a position and time depending component. Dynamics are modeled by a deterministic LSTM that predicts the mean and the variance of a normal distribution at each time step. To the best of our knowledge, this is the only GNN based method that also reports NLL.
(in) Multiple Futures Prediction (MFP) (Tang& & Salakhutdinov 2019): A recurrent model with deterministic transition dynamics. Stochasticity is introduced via the initial state and noisy observations that are fed back into the dynamical model. Interactions are modeled by an attention module.
5.2.2 Results
We provide benchmark results of our proposed model in Tab. 2 Similar to the experiments on the rounD dataset in Sec. 5.1 we observe that (i) deterministic training and testing is more efficient than its Monte Carlo based alternative and (ii) increasing the number of modes improves the performance.
Next, we compare our method, using one component in the GMM prior only, with all other unimodal prediction methods. We can observe that our approach outperforms its competitors in terms of NLL. Furthermore, our method gives similar RMSE as other methods and only outperformed by MFP. However, and in contrast to the other methods in the benchmark, MFP reports the best RMSE over 5 Monte Carlo samples, which makes it difficult to draw fair conclusion.
17

---

Published in Transactions on Machine Learning Research (03/2023)
We subsequently increase the number of modes for model and for MFP. In terms of NLL, our method achieves superior results when it comes to long-term predictions (3s, 4s, 58), while MFP performs better for short-term predictions (Is, 28). Accurate long-term prediction of the agents in a driving scene is crucial for high-lovel autonomous driving, as an accurate environment model is a prerequisite for precise planning of driving controls. For instance, advanced driver-assistance systems usually use time horizons between 38 and 5s for driver warnings and emergency brakes, while autonomous cars aim for a time horizon for 5s or longer in order ensure safe and comfortable rides (Philipp & Goehring 2019 We provide minRMSE values for different number of modes measure of predictive diversity in App.
Table 2: NGSIM results. We provide RMSE and NLL averages and standard errors over 10 runs. For MFP, we report the best RMSE values over 5 Monte Carlo samples.
Const, Vel. I Mode
CS-LSTM 1 Mode
ST-LSTM 1 Mode
MFP 1 Mode
MFP Modes
GDSSM 1 Mode PLL/MC
GDSSM 1 Mode PLL/Det.
GDSSM
GDSSM 3 Modes PLL/Det.
GDSSM 4 Modes PLL/Det.
2 Modes PLL/Det.
18
0.75
0.61
0.51
0.54
0.54
0.56 1: 0.00
0.53 1 0.01
0.55 1 0.00
0.55
0.01
0.57 1: 0.02
RMSE
28
1.81
1.27
1.21
1.16
1.17
1.27 + 0.02
1.18 + 0.01
1.22 + 0.01
1.23
0.02
1.24 + 0.04
38
3.16
2.09
2.01
1.90
1.91
0.03
1.98 + 0.02
2.02 + 0.02
2.05
0.03
2.06 +
-180
3.10
3.01
2.78
2.75
3.16
0.01
2.99
0.03
3.03
0.04
3.06
0.05
3.05 0.06
58
6.70
4.37
4.31
3.83
3.78
4.47 - 0.05
1.29 = 0.04
1.33 = 0.07
4.33
0.07
4.35 = 0.08
18
0.80
0.58
0.90
0.73
-0.65
0.64 1: 0.04
0.19 1 0.02
-0.12
0.02
.0.15 1. 0.03
-0.16 1: 0.03
28
2.30
2.14
2.41
2.33
1.19
2.10 + 0.10
1.61 + 0.02
1.37 + 0.02
1.35 + 0.02
1.32 + 0.02
NLL
36
3.21
3.03
3.25
3.17
2.28
2.92 + 0.11
2.42 + 0.02
2.23 + 0.02
2.23 + 0.02
2.19 + 0.02
48
3.89
3.68
3.61
3.77
3.06
3.54 0.10
3.02 = 0.02
2,88 - 0.02
2.88 0.02
2.82 + 0.02
58
4.44
4.22
4.36
3.69
0.10
3.50
0.02
0.02
3.36
5.3
Covariance Approximations
We provide a detailed runtime analysis for different covariance approximations Sec. 5.3.1 and study their impact on the predictive performance in Sec. 5.3.2
5.3.1 Runtime
We visualize the runtime different covariance approximations, as well as the runtimo of the Monte Carlo alternative in Fig. 5 as a function of input dimensionality and number of agents. We use the same NSDE architecture as in our experiments on the rounD and NGSTM dataset. For the Monte Carlo alternative, we visualize the runtime for 16 particles, as we use the same number of particles for training in Sec. 5.1 and
5.2
We first confirm that propagation of the full covariance matrix is more costly than any of the proposed approximations. In fact, our sparse covariance approximations can reduce the runtime up to factor of 100 for systems with a large number of agents and a high input dimensionality. As we derived in Sec. 4.4 when using the BMM algorithm with a full covariance matrix or the all diagonals approximation, the computational cost shows cubic dependence on the number of agents. In contrast, the main diagonal approximation, the main blocks approximation, as well MC based predictions have a quadratic dependence on the number of agents. This makes these two approximations an attractive alternative to MC based predictions, when systems with a high number of agents need to be modeled with a limited computational budget.
For systems with a moderate input dimensionality (+s 8) and number of agents (>> 16), all of our proposed approximations require only up to 5 ms in order to compute the distribution at the next time point, which corresponds to the cost of 5-10 Monte Carlo simulations.
5.3.2 Benchmark
Next, we study the effect of different covariance approximations on the performance. We report the results for the case unimodal GDSSM, The results are depicted Tab. We report here our main findings.
First, modeling the full covariance matrix results in the best performance terms of lowest RMSE and NLL. Second, the all diagonals covariance approximation performs the best among the covariance approximations.
18

---

Published in Transactions on Machine Learning Research (03/2023)
Full
Main Diagonal
Main Blocks
All Diagonals
Monte Carlo
102
2 Agents
2 Agents
101
2 Agents
2 Agents
10a
2 Agents
Time in
Time in
Time in
10-1
Time
10
10-3
10
21
2
2
23
2
27
2
2
2
2
2
27
2
2
Input Dimensionality
Input Dimensionality
Input Dimensionality
Input Dimensionality
Input Dimensionality
Figure 5: Wallclock time for output moment calculation for the latent dynamics using GNNs with three hidden layers of size 24 for different covariance approximations (from left to right). For each approximation, we plot the runtime as function of the input dimensionality D2 and number of agents M. The GNNs are initialized at random and we report the average runtime over 100 repetitions.
It achieves comparable RMSE as the full solution, but falls slightly behind in NLL when the system is highly interactive (see rounD dataset). In this setting, it outperforms the other two sparse approximations with respect to NLL. This behavior might be explained by the assumptions made in the covariance approxima- tion: it is the only approximation that allows for correlations between agents as its structure only neglects dependencies between latent features.
Third, the differences performances between the full solution and the different approximations with respect to RMSE lie between one and two standard errors. Modeling the main diagonal only can thus be sufficient for applications with low computational resources and high demand in accuracy by accepting slight loss in calibration. For these applications, the runtime can be reduced from O(M3) to O(M²).
Table 3: Test performance for different covariance approximations the rounD and NGSIM datasct for the uninodal case. We provide average and standard error over 10 runs.
NGSIM
Full
Main Diagonal
Main Blockss
All Diagonals
Full
Main Diagonal
Main Blocks
All Diagonals
18
0.79 + 0.02
0.82 + 0.02
0.83 = 0.01
0.78 = 0.02
0.53 # 0.01
0.54 # 0.01
0.55 = 0.01
0.53 + 0.01
RMSE
28
1.87
0.02
1.88 + 0.02
1.89
0.05
1.87 + 0.02
1.18 + 0.01
1.18 + 0.02
1.19 + 0.02
1.19 0.02
3.36
0.03
+ 0.02
0.06
3.34 +
+ 0.02
1.99
+ 0.03
1.99
5.08
0.04
5.07
0.04
5.09
0.07
5.05 1 0.03
2.99
0.03
2.98
3.07
0.04
2.97
0.03
58
7.21
+ 0.05
7.25
+ 0.06
7.30
+ 0.08
7.25 + 0.05
1.29 + 0.04
4.31 : 0.05
4.54 4 0.06
4.32 1: 0.04
18
1.48
0.05
1.77 +: 0.06
1.79
0.05
1.69 +: 0.03
0.19
+: 0.02
0.24 4: 0.04
0.22 4 0.08
0.18
0.03
0.03
0.06
0.03
1.59
3.87
0.02
4.27
0.03
4.28
0.05
4.18 1 0.02
2.42
0.02
2.44
2.46
1 0.03
2.43
0.02
48
4.46
- 0.08
5.00
+ 0.02
5.01
= 0.05
1.92 + 0.02
3.02
- 0.02
3.01 $ 0.02
3.06
- 0.01
3.01 :: 0.02
5.05
1 0.04
5.69 1 0.02
5.68 1 0.06
5.64 1: 0.05
3.50 1 0.02
3.59 1 0.02
3.62 1 0.03
3.57 1 0.02
5.4
Testing
We analyze the generalization capabilities of our model by testing it on out-of-distribution data, e.g. traffic environments that have not. been observed during training.
For this experiment, we reuse the rounD datasct KKajewskii et al.] 2020) since consists of record- ings at three different roundabouts. Here, each roundabout corresponds to separate traffic environment. We select one recording from each roundabout. (KackerstraBe in Aachen, Thiergarten in Alsdorf and
19

---

Published in Transactions on Machine Learning Research (03/2023)
Neuweiler near Aachen) and apply the same data curation steps in Sec. 5.1 In order to generalize between different traffic environements, we change the experimental setup follows:
Local coordinate system: We transform the data into local coordinate system, which is centered at the cgo-vchicle and is oriented to the heading direction of the ego-vehicle.
Include map information to the context variable Z. We extract a local map around the ego- vehicle, which spans a rectangle with a length of 74 meters and a width 44 meters. Afterwards, we apply a binary masking to the map which divides the image into drivable and non-drivable areas. Our task is to jointly model all vehicles, which lie within this rectangle.
More details on the preprocessing can be found in App. G
5.4.1 Results
Next, we analyze the gencralization capabilitics of our model by comparing the following strategies: (1)
training on two traffic environments and testing on third distinct traffic environment, (2) directly training the model the test traffic environment and (3) training on all traffic environments. Tn order to enable fair comparison, the data for each traffic environment is split into non-overlapping training and test sets.
We present the results for different traffic environments in Tab. The first column in Tab. 4 de- scribes the RMSE and NLL when we train our model on the traffic environments Thiergarten (T) and Neuweiler (N) and test it on the traffic environment Kackertstraße (K). The second column describes the predictive performance by using scenes from the same traffic environment (K) for training and testing. The third column describes the predictive performance by training all three traffic environments (KNT) and testing on the traffic environment K. The remaining columns benchmark the generalization capabilities of our model on the traffic environments T and N and are set up in an analogous fashion.
In all experiments, we observe that the performance increases from (1) using different training envi- ronments for training and testing, to (2) using the same enviroument for training and testing and (3) using all environments for training. The difference in performance between (1) and (2) is moderate for the locations K and T demonstrating that our model is capable of generalizing unseen traffic environments during testing. For the location N, the performance difference is increased which can be explained by studying the locations in more detail: N is multi-lane roundabout with congested traffic, the roundabouts K and T are single-lane with moderate traffic. In consequence, the out-of-domain test on the traffic environment is more challenging. We visualize exemplary predictions of our model GDSSM Fig.
Table 4: Test performance on different traffic environments on the rounD dataset using our method GDSSM (1 mode). We vary for each test traffic environment the training traffic environments. We provide averages and standard errors over 10 runs. Kackertstraße- K, Thiergarten--T Neuweiler-N.
Test Train
K
N
TN
K
RTN
KN
KTN
KT
N
KTN
1s
2.12 + 0.09
1.88 + 0.05
1.55
+ 0.04
1.42 + 0.04
1.49 + 0.05
1.21 + 0.04
2.98+
0.12
2.02 + 0.07
1.55 + 0.04
RAISE
4.56
0.13
3.73 +
3.34
0.09
2.82
0.05
2.53
0.09
2.27
+
5.82=
0.19
3.38
0.06
3.02 +
7.51
0.24
6.05
4: 0.17
5.62
0.20
4.42
0.12
3.67
0.12
3.35
4 0.07
8.96
0.24
4.50
0.13
4.91
0.18
48
10.57
= 0.34
8.65
+ 0.26
0.17 7.99
6.75 4: 0.18
5.33
+ 0.15
5.00 + 0.16
12.44
0:33
681
+ 0.14
6.91 + 0.18
13.00
0.43
10.73
7.59
8.77
0.19
4.11 à. 0.07
4.02 1 0.04
3.13
0.06
3.35
0.06
3.16
0.06
2,88 1 0.06
4.75 0.07
3.84
0.04
3.21 1: 0.08
2x
5.38
+ 0.10
4,84 + 0.03
4.86
0.00
4.22
0.04
3.99 + 0.07
3.65 + 0.04
5.82+ 0.09
4,60 + 0.05
4.42 + 0.09
6.55
0.17
6.42
0.17
5.49
0.10
4.98
7.32
0.17
5.65
0.04
48
7.71
à 0.20
7.07 1 0.14
7.57
0.21
7.20
à: 0.08
6.38 4 0.19
5.80 4: 0.13
9.65 0.18
6.70
4: 0.11
6.71 :: 0.08
0.18 8.63
8.14 + 0.18
0.15 8.08
0.12 8.54
7.46
0.21
+ 0.11
9.11+ 0.18
7.54
0.22
7.04 +
20

---

Published in Transactions on Machine Learning Research (03/2023)
Prediction: Mode Weight: 0.40
Prediction: Mode Weight: 0.60
Predictian: Made
Prediction: Mode: 2
Weight: 0.57
Weight: 0.43
nel
nvel
(a) Kackertstraße (K)
(b) Thiergarten (T)
Figure 6: Predictions of our model GDSSM (2 modes) on out-of-domain traffic enviromments, i.e. the training and testing traffic environments are distinct. Given the history of each traffic participant (dashed lines), we visualize the predicted 95% confidence interval. Solid lines represent the true future trajectory.
6 Conclusion
In this work, we have proposed GDSSMs in which the latent dynamics of the agents are coupled via. GNNs in order to capture interactions among multiple agents. We derived moment matching rules for GNN layers that allow for deterministic inference and introduced a GMM prior over the initial latent states in order to allow for multimodal predictions. Both together lead to an efficient and stable algorithm that is able to produce complex and nonlinear predictive distributions. We confirmed that our novel method shows strong empirical performance on two challenging autonomous driving datasets. Finally, we proposed sparse approximations to the covariance matrix considering the computational limits of real-world vehicle control units. Depending the required calibration, approximations can lead to significant reduction of the runtime without impeding accuracy.
In future work, we seek to increase the robustness of our proposed model lowards novel and unseen traflic scenarios. One way could be to incorporate epistemic uncertainty into our model formulation by placing a prior over the weights of the GNN. To achieve this, we could combine our model with recent advances variational inference in order find an approximation the intractable weight posterior. Another research direction of interest. modeling of irregular and partially observed dynamical systems. Prior work uses continuous time encoder networks well as continuous time transition model in latent space (Rubanoval ct al. 2019 Brouwer et al. 2019b. Following this vein of work, an extension of our model towards continuous time networks seemss a promising direction for modeling interactive systems.
References
Daniel Alspach and Harold Sorenson. Nonlinear Bayesian estimation using Gaussian sum approximations. IEEE transactions on automatic control, 17(4), 1972.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling. arXiv, abs/1803.01277 2018.
Mayank Bansal, Alex Krizhevsky, and Abhijit S. Ogale. ChauffeurNet: Learning Drive by Imitating the Best and Synthesizing the Worst. arXiu, abs/1812.03077) 2018.
Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanchez-Gonzaler Vinícius Flores Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Çaglar Gülçehre, Н. Francis Song, Andrew Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey R. Allen, Charles Nash, Victoria Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmoct Kohli, Matthew Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu. Relational inductive biases, deep learning, and graph networks. arXiv, 0bs/1806.011 2018.
Justin Bayer, Maximilian Soelch, Atanas Mirchev, Baris Kayalibay, and Patrick van der Smagt. Mind the Gap when Conditioning Amortised Inference in Sequential Latent-Variable Models. In ICLR, 2021.
21

---

Published in Transactions on Machine Learning Research (03/2023)
Samy Bengio, Oriol Vinyals, Navdeep Jaitly, and Noam Shazeer. Scheduled Sampling for Sequence Prediction with Recurrent Neural Networks. In NewrIPS, 2015.
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational Inference: A Review for Statisticians.
Journal of the Americon Statistical Association, 112(518), 2017.
Michael W. Brandt and Pedro Santa-Clara, Simulated Likelihood Estimation of Diffusions with an Appli- cation to Exchange Rate Dynamics in Incomplete Markets. Journal of Financial Economics, 63(274), 2002.
Edward De Brouwer, Jaak Simm, Adamn Arany, and Yves Moreau. GRU-ODE-Bayes: Continuous Modeling of Sporadically--Oservell Time Series. In Neur/PS, 2019.
Yuri Burda, Roger Baker Grosse, and Ruslan Salakhutdinov. Importance Weighted Autoencoders. In ICLR, 2016.
Sergio Casas, Cole Gulino, Simon Suo, Katie Luo, Renjie Liao, and Raquel Urtasun. Implicit Latent Variable Model for Ccene-Consistent Motion Forecasting. In ECCV, 2020.
Yuning Chai, Benjamin Sapp, Mayank Bansal, and Dragomir Anguelov. MultiPath: Multiple Probabilistic Anchor Trajectory Hypotheses for Behavior Prediction. arXiv, abs/1910.054419, 2019.
Guangxi Chen, Ling Hu, Qieshi Zhang, Ziliang Ren, Xiangvang Gao, and Jun Cheng. ST-LSTM: Spatio-
Temporal Graph Based Long Short-Term Memory Network For Vehicle Trajectory Prediction. In IEEE ICIP, 2020.
Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio. Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling. arxiv, abs/1411.3355. 2014.
Junyoung Chung, Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron C Courville, and Yoshua Bengio. A Recurrent Latent Variable Model for Sequential Data. In NewrlPS, 2015.
Henggang Cui, Vladan Radosavljevic, Fang-Chieh Chou, Tsung-Han Lin, Thi Nguyen, Tuu-Kuo Huang, Jeff Schneider, and Nemanja Djuric. Multimodal Trajectory Predictions for Autonomous Driving using Deep Convolutional Networks. In ICRAA, 2019.
Nachiket Deo and M. Trivedi. Convolutional Social Pooling for Vehicle Trajectory Prediction. In CVPR Workshop, 2018.
Frederik Diehl, Thomas Brunner, Michael Le, and Alois Knoll. Graph Neural Networks for Modelling Traffic Participant Interaction. In IEEE Intelligent Vehicles Symposium, 2019.
Nemanja Djurie, Vladan Radosavljevic, Henggang Cui, Thi Nguyen, Fang-Chiet Chou, TIsung-Han Lin, Nitin Singh, and Jeff Schneider. Uncortainty-aware Short-term Motion Prediction of Traffic Actors for Autonomous Driving. In IEEE WACV, 2020.
Arnaud Doucet, Simon Godsill, and Christophe Andrieu. On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and computing, 10(3), 2000.
Ola Flerian, Siddhartha Chib, and Neil Shephard. Likelihood Inference for Discretely Observed Nonlinear Diffusions. Econometrica, (9(,), 2001.
Marco Fraccaro, Søren Kaae Sønderby, Ulrich Paquet, and Ole Winther. Sequential Neural Models with Stochastic Layers. In NewrIPS, 2016.
Yarin Gal and Zoubin Ghahramani. Dropout as Bayesian Approximation: Representing Model Uncertainty in Deep Learning. In ICML, 2016.
Pascal Germain, Francis Bach, Alexandre Lacoste, and Simon Lacoste-Julien. PAC-Bayesian Theory Meets Bayesian Inference. In NewrIPS, 2016.
22

---

Published in Transactions on Machine Learning Research (03/2023)
Soumya Ghosh, Francesco Delle Fave, and Jonathan Yedidia. Assumed Density Filtering Methods for Learn- ing Bayesian Neural Networks. In AAAT, 2016.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural Message Passing for Quantum Chemistry. In JCML, 2017.
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-31) 2007.
Danijar Hafner, Timothy P. Lillicrap, Ian Fischer, Ruben Villegas, David Ha, Honglak Lee, and James Davidson. Learning Latent Dynamics for Planning from Pixels. In ICML, 2019.
John Halkias and James Colyar. Us highway 101 dataset. Technical report, Federal Highway Administration, 2007. URL httpp://roosp..tl. wts. htt.b//..o////iee/////////////// http://rrrasp.mm...ooooorrrr///r
William L. Hamilton, Rex Ying, and Jure Leskovec. Inductive Representation Learning on Large Graphs. In NewrIPS, 2017.
Trevor Hastie, Robert Tibshirani, Jerome Friedman, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction, volume Springer, 2009.
Manuel Haussmann, Sebastian Gerwinn, and Melih Kandemir. Bayesian Evidential Deep Learning with PAC Regularization. arXiv, abs/1906.00816, 2020.
Michael Herman, Jörg Wagner, Vishnu Prabhakaran, Nicolas Möser, Hanna Ziesche, Waleed Ahmed, Lutz Bürkle, Ernst Kloppenburg, and Claudius Gläser. Pedestrian Behavior Prediction for Automated Driving: Requirements, Metrics, and Relevant Features. IEEE Transactions on Intelligent Transportation Systems, 23(9), 2022.
Jose Miguel Hornandez-Lobato and Ryan Adams. Probabilistic Backpropagation for Scalable Learning of Bayesian Neural Networks. In ICML, 2015.
Sepp Hochreiter. Untersuchungen zu dynamischen neuronalen Netzen. Diploma, Technische Universitäl München, 91(1), 1991.
Ajay Jain, Sergio Casas, Renjie Liao, Yuwen Xiong, Song Feng, Sean Segal, and Raquel Urtasun. Discrete Residual Flow for Probabilistic Pedestrian Behavior Prediction. In CoRL, 2019.
Ashesh Jain, Amir Roshan Zamir, Silvio Savareso, and Ashutosh Saxena. Structural-RNN: Deep Learning on Spatio-Temppraa Graphs. In CVPR, 2016.
Robert Krajewski, Tobias Moers, Julian Bock, Lennart Vater, and Lutz Eckstein. The rounD Dataset: A Drone Dataset of Road User Trajectories at Roundabouts in Germany. In ITSC, 2020.
Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured Inference Networks for Nonlinear State Space
Models. In AAAI, 2017.
Tuan Anh Le, Maximilian Igl, Tom Rainforth, Tom Jin, and Frank Wood. Auto-Encoding Sequential Monte Carlo. ICLR, 2018.
David Lenz, Frederik Diehl, Michael Truong-Le, and Alois C. Knoll. Deep Neural Networks for Markovian Interactive Scene Prediction in Highway Scenarios. In IEFE Intelligent Vehicles Symposium, 2017.
Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan. Enhancing the Locality and Breaking the Memory Bottleneck of Transformer on Time Series Forecasting. In NewrIPS, 2019.
Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting. In ICLR, 2018.
23

---

Published in Transactions on Machine Learning Research (03/2023)
Yunzhu Li, Hao He, Jiajun Wu, Dina Katabi, and Antonio Torralba. Learning Compositional Koopman Operators for Model-Based Control. In ICLR, 2020.
Andreas Look, Melih Kandemir, Barbara Rakitsch, and Jan Peters. Neural SDEs. IEEE TPAMI, 2022.
Deterministic Approximation
Chris J Maddison, John Lawson, George Tucker, Nicolas Heess, Mohammad Norouzi, Andriy Mnih, Arnand Doucet, and Yee Teh. Filtering Variational Objectives. In NeurIPS, 2017.
Jcan Mcreat, Nicole Zoghby, Guillaume Sandou, Dominique Beauvois, and Guillermo Gil. Inertial Single Vehicle Trajectory Prediction Baselines and Applications with the NGSIM Dataset. arXio, abs/1908..14777 2019.
Christian Naesseth, Scott Linderman, Rajesh Ranganath, and David Blei. Variational Sequential Monte Carlo. In AISTATS, 2018.
Aaron den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu. WaveNet: A Generative Model for Raw Audio. arXiv, abs/1609.03499, 2016.
Soumyasundar Pal, Liheng Ma, Yingxue Zhang, and Mark Coales. RNN with Particle Flow for Probabilistic Spatio-temporal Forecasting. In ICML, 2021.
Asger Roer Pedersen. A New Approach to Maximum Likelihood Estimation for Stochastic Differential Eqqatioon Based on Discrete Observations. Scandinavion Journal of Statistics, 22(1), 1995.
Andreas Philipp and Daniel Goehring, Analytic Collision Risk Calculation for Autonomous Vehicle Naviga- tion. In ICRA, 2019.
Yulia Rubanova, Ricky T. Q. Chen, and David K Duvenaud. Latent Ordinary Differential Equations for Irregelarly-Sampled Time Series. In NewrIPS, 2019.
Alvaro Sanchez-Gonzale Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to Simulate Complex Physics with Graph Networks. In ICML, 2020.
Simo Särkkä. Bayesian Fillering and Smoothing. Cambridge University Press, 2013.
Simo Särkkä and Juha Sarmavuori. Gaussian Filtering and Smoothing for Continuous-Diserete Dynamic Systems. Signal Processing, 93(2), 2013.
Simo Särkkä and Arnold Solin. Applied Stochastic Differential Equations. Cambridge University Press., 2019.
Simo Särkkä, Jouni Hartikainen, Isambi Sailon Mbalawata, and Heikki Haario. Posterior Inference Parameters of Stochastic Differential Equations via Non-Linear Gaussian Filtering and Adaptive MCMC. Statistics and Computing, 25(2), 2015.
Thomas B Schön, Adrian Wills, and Brett Ninness. System Identification of Nonlinear State-Space Models.
Automatica, 47(1), 2011.
Thomas B Schön, Fredrik Lindsten, Johan Dahlin, Johan Wägberg, Christian A Nacsseth, Andreas Svensson, and Liang Dai. Sequential Monte Carlo Methods for System Identification. IPAC, 48(28), 2015.
Arno Solin, Ella Tamir, and Prakhar Verma. Scalable Inference in SDEs by Direct Matching of the Fokker- Planck-Kolmogorov Equation. In NewrIPS, 2021.
Charlie Tang and Russ R Salakhutdinov. Multiple Futures Prediction. In NewrIPS, 2019.
Belinda Tzen and Maxim Raginsky. Neural Stochastic Differential Equations: Deep Latent Gaussian Models in the Diffusion Limit. arXiv, abs/1905.09683, 2019.
24

---

Published in Transactions on Machine Learning Research (03/2023)
Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian Fluid Simulation with Continuous Convolutions In ICLR, 2019.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Milia Polosukhin. Attention Is All You Need. In NewrIPS, 2017.
Tim A. Wheeler and Mykel J. Kochenderfer. Factor Graph Scene Distributions for Automotive Safety Analysis. In ITSC, 2016.
Anqi Wu, Sebastian Nowozin, Edward Meeds, Richard E. Turner, Jose Mignel Hernandez-Lobato, and Alexander L. Gaunt. Deterministic Variational Inference for Robust Bayesian Neural Networks. In ICLR, 2019.
Fan Yang. Ling Chen, Fan Zhou, Yusong Gao, and Wei Cao. Relational State-Space Model for Stochastic Multi-Object Systems. In TCLR, 2020.
Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-Tempores Graph Convolutional Networks: Deep Learn- ing Framework for Traffic Forecasting. In IJCAL, 2018
25

---

Published in Transactions on Machine Learning Research (03/2023)
A Proof of Theorem i
Theorem 1. The marginal distribution p(y: (I)
is analytically computed as
P(y Σ To (I)
for a GDSSM with the below generative model
Cat((IT(II)) (I)...., TV (I)),
zu ~ ((H,, ding(Eo,
2, (x + f(1, v, -1- diog (L(t, u, Z))), y1 ~ N (yrlg(t, v, I),,, diag vZ)),
t = I. T T
where f(1,, v; T(t, v, g(, v, I), uI, are time i, component u, and contect I depending matrices with appropriate dimensionality.
Proof. The proof is straightforward as the output moment of the transition and emission function are ana- lytically available
/ oo
=
Σ To( /
V
di,
Σ f
Σ (I)N a (I), (I))
The moments at time step t of a linear time depending dynamical system are available as
Särkkä 2013
(T) II +
t==1
(I) =
II (I + f(l,
f(f,,,I)) v.z1)
diaggEo,,(2
f(t ((F,
Σ
II (I+ (I+ v0.21) Liag(L(t', #'t'41
II (I+ f(t", o,)) Lt" =0"+41
+ diag(L(t, v)
We obtain the same expression via the BMM algorithm, which is easy to prove by inserting the locally linear system into Eq. G Finally, mean az (Z) and covariance (Z) of the output at time step t are available as
Särkkä 2013
(Z) Egg(t, v, D)x] = g(t, v, I),, (I),
(I) = Covig(t, v, [)x + diag (E)) = g(t, v, I)E,, (I)g(t, v, I)T - diag (T v. 2)).
26

---

Published in Transactions on Machine Learning Research (03/2023)
B
Output Moments for the ReLU Activation Function
Suppose, we apply the ReLU activation function layer and time step to the input of
zt ++ = max(0,
Output Moments We can approximate the output moments of xz+1 as function of the input moments Wu et al. 2019
E mmmax(0,22)) (0, R diag(Cov((z)) SR
B[r]]/ diag(Cov(a?))
Cov Cov|max(( (0, ,r)] V diag((oowss))diiggggoo
F (B[r)), Cov(z!])
where SR(x) = ((x) + (x) with and representing the PDF and CDF of standard Gaussian variable. The function [x]] Cov( (2()) is a shorthand notation for
F(E[//] Cov[e:)) = A (B[a)), Cov(z!)) exp [-Q (Bja)), [ov[x!
The function A(B[x]], [ov(x!]) can be estimated
A(B[a/), Cov[z = SR())SR(()) + (()())))))
where p Cov[zi]/
is dimensionless matrix and dh
E[x!// / diag(Cov((z)) is dimensionless vector.
The i, j-th element of ((E[a)), Cov(z!) can be estimated as:
Plas
aresin cresin(pt..) (Ph.,) plaus de. ch
Ikus 2x
Q [Eve[), 20ken (1 + Prass)
(1d2)) +
loog
Phi., Pkas
where of. arcsin(p)) + ph 1+ 1++11-ph 2 0.kk
The operator @ denotes the elementwise disivision and
② denotes the elementwise multiplication.
Expected Jacobian Since activation functions are applied clement-wise off-diagonal entrics of the expected gradient zero. The diagonal of the Jacobian of the ReLU function is the expected Heaviside step function Wu et al. (2019)
diag E V max (0, x
For more detailed derivation, refer to the work of Wu et al. 2019
C
Training Details
We train all models with the ADAM optimizer and stochastic mini-batches. We use batch size of 4 and a learning rate of 0.0001. In order to accelerate training of multi-modal GDSSMs we initialize the transition and observation neural nets with the pretrained versions of the uni-modal GDSSMs.
C.1 DounD
We train deterministic models (GDSSM Det. and Non Recur. GNN) for 50k weight updates. For the MC based models (GDSSM MC) we use 100k weight updates. The dataset contains 4,314 training and 1,091 testing snippets. We do not separate validation dataset. we observed no overfitting.
27

---

Published in Transactions on Machine Learning Research (03/2023)
C.2 NGSIM
We train all models for 1000k updates the NGSIM dataset. The dataset contains 5,922k /60k// 860k/1,5054 aaaiin/vliaatio///n snippets. We validate the models on 1K random minibatches from the validation dataset after every 10k weight updates.
C.3 Out- DfDtiiibibbtti Testing
We train all models for 50k weight updates. The dataset consists of three sub-datasets with 254/ 251/ 137 snippets, where 80% of each sub-dataset are used for training and the remaining 20% for testing. Due to the small size of the sub-datasets observed overfitting, We address this by using the last 20% of each training sub-dataset for validation. We validale the model after every 1k weight updates.
D Evaluation Metrics
In the following, we give a brief overview on different evaluation metrics. We provide the negative log- likelihood (NLL) and the Root-Mean-Square Error (RMSE) in the main of the paper, and minRMSE in App.
Root-MeannSquare Error
The (root) mean-square error at time point is defined as
M(n) (yin
Din) 7 (ui,, Vin) RMSE (t) MSSe (t),
MSF Bass (1)
Σ
(25)
N
M(n)
where N is the number of snippets in the datasct, (() is the number of agents the n-th snippet, yf,, is the true location of the m-th agent at time point of the n-th snippet, and 7in the corresponding predicted location.
Unfortunately, it is not straight-forward to extend RMSE to probabilistic forecasts. To the best of our knowledge, none of the existing variants is a strictly proper scoring rule (Gneiting & Raftery 2007 where the latter is optimal if and only if the predictive distribution matches the true distribution. An example of a strictly proper scoring rule that we also provide in the paper is the predictive NLI:
N M(x)
logp(Bi..1 (Z)
NLL(t)= Σ
(26)
M(n)
For deterministic GDSSMs, we approximate PT directly with our method that we introduced in Sec. 4.2 For GDSSMs that do not follow an assumed density approach, we approximate p(u?" [Z) via. Monte
Carlo integration.
In the following, discuss different RMSE variants in more detail.
Bayes Predictor: common evaluation metric used in the ML community (e.g.
Gal & Ghahramani
2016)) is to compute
1 N
M(n) (y
ElOr) - E M(n)
RMSE(t)
(27)
where we assume that the probabilistic predictor will perform model averaging before making its final prodiction. This score only evaluates the goodness of the first moment and docs not take any information of higher moments into account. We provide its values in the main part of the paper.
Gibbs Predictor: An alternative way for extending deterministic loss functions to the probabilistic sce- nario be found in the PAC-Bayes setting (e.g. Germain et al. (2016)) ) by considering the average loss in
28

---

Published in Transactions on Machine Learning Research (03/2023)
the risk function. In our case, this would lead to
N M(n)
(v D(,)) ( Hin)
RMSEGibba(t)
E
(28)
M()
where we take expectation with respect the prodictions Vi However, we discourage take this for- mulation as evaluation metric since it penalizes any form of variance in the predictor and favors Dirac distribution around the expected value of the true distribution. Note that this result can be tightly linked to the bias-variance trade-off in statistical learning theory (e.g. Hastie et al. (2009)).
Min Predictor: A popular alternative in the traffic forecasting literature (e.g. Tang & Salakhutdinor 2019h) is to compute the minimal error over a set of potential predictions
M(a) Σ
(y Vi,,,s) )T (vi, Di,,,)
minRMSE(1)
min
(29)
M(n)
where (0)h,)))= the set of S potential predictions. This score favors models that produce set of diverse predictions where at least one candidate is close to the true outcome but ignores the goodness of the remaining S predictors. While gives some useful information about the model capacity, does not give complete representation of the test-time performance when the best mode is unknown. We provide its values in App.
E
Alternative Parameter Inference Methods
One commonly used inference method in the context of machine learning circumvents maximizing the log- likelihood and instead maximizes the Evidence Lower Bound (ELBO)
ELBO VT (Z) = =B (( log p(uI..... WT, 3oo IT [Z)
(30)
(000 -.00-..,37. FT)
that involves learning an approximation (x IT) to the intractable smoothing distribution ... VT, I) Krishnan et al. 2017
A tighter bound to the log-likelihood log bog.(u...... yr I) can be obtained by calculating the impor- tance weighted log-likelihood, which we refer to as the Monte Carlo Objective (MCO) (Maddison et al. 2017 Burda et al. 2016
...... MCO .r.......I]]] IZ) Ey(xa.....zr)) log plvn..... Vr Fo, FT. II]
(31)
(F ..... TT,k)
where K is the number of Monte Carlo samples and Fe, is the k-th sample time step t. For state-space models, recent work combined the MCO with particle filters (Naesseth et al. 2018 Maddison et al. 2017
Le et al. 2018
F Extended Results
We provide minRMSE values for the rounD and NGSIM datasct. We set the proposal predictions to the mean values of each mode and the min operator selects for each agent the proposal that produces the lowest error per snippet over the complete prediction horizon. This score favors models in which the proposal predictions are diverse as long as at least one candidate is close to the true outcome. As shown in Tab. the minRMSE decreases as we increase the number of modes, which indicates that the different modes correspond to diverse set of plausible trajectories.
29

---

Published in Transactions on Machine Learning Research (03/2023)
Table 5: minRMSE values as a function of the number of modes on the rounD and NGSIM dataset. We provide average and standard error over 10 runs. For RMSE and NLL results, please refer Tab.
and
Tab. 2
rounD
NGSIM
Mode
2 Modes
3 Modes
4 Modes
1 Mode
2 Modes
Modes
4 Modes
0.79 + 0.02
0.75
0.04
0.76 +: 0.03
0.69
0.02
0.53 0.01
0.46 + 0.01
0.35 0.01
0.32 + 0.01
minRMSB
1.87 = 0.02
185 0.06
1.83 + 0.07
1.68 0.05
1.18 0.01
1.05 0.01
0.50 0.01
0.71 + 0.01
3.36 0.03
3.46
1 0.11
3.05 = 0.16
2.75
0.12
1.98 0.02
1.73 0.02
1.33
1.18 1 0.03
5.08
0.04
5.07
0.13
4.55
- 0.28
4.04
0.45
2.99 0.03
2.60 0.04
2.02
1.75 +\  0.04
58
7.24 + 0.05
6.29
0.23
5.95
+: 0.47
4.43 0.53
4.29 0.04
3.69 0.06
2.88 0.05
2.46 + 0.05
G
Experimental Setup for
Testing
G.1 Dataset Construction
We construct dataset, which consists of three different traffic environments. We select one recording from the roundabout in Kackortstraße (K) in Aachen, one recording from the roundabout in Thiergarten (T) in Alsdorf, and one recording from the roundabout in Neuweiler (N) near Aachen. We use the same preprocessing procedure as in Sec. [5.1] We remove pedestrians, bicycles, and parked vehicles from each traffic environment. There remain 319/ 264/ 389 tracked objects over a time span of 0.3/ 0.3/ 0.15 hours the traffic environments K/ T/ N. We downsample the recordings by factor 5 and then construct for each traffic environment dataset which consists of 8s long snippets with 50% overlap. The first three seconds are used as the track history and the following five seconds the prediction horizon. We obtain 254/ 251/ 137 snippets for the traffic environments K/ T// N. For each traffic environment we use the first 80% snippets for training and the remaining 20% snippets for testing.
G.2 Map Processing
The map processing follows existing work in the domain of traffic forecasting Bansal et al. 2018 Herman ct al. (2022). Given an ego-vehicle and its history, we first calculate the heading of the vehicle. The heading is calculated as the average heading direction over the last 0.2 observed seconds. The position and the heading direction jointly define the Region-of-Interea (ROI) on the map. The ROI a rectangle with a length of 74 meters and width of 44 meters, which is centered at the ego-vehicle and oriented according the heading of the eeo-vehicle. We further convert the RGB-image into a binary road image. We visualize
the processed map information in Fig. 7b
Eqo Vehidee Hiadding
eeve
(a) Full map with ego- vehicle.
(b) Rotated, cropped, and masked map information.
Figure 7: Processing of map information.
30

---

Published in Transactions on Machine Learning Research (03/2023)
H Network Architectures
H.1 Network Architectures without World Models
Below we present the neural network architectures for the experiments Sec.
5.2 5.3
State
Adjacency
State M x 4
Adjacency
State M x
M x
M x M
M x M
Mean Aggregator
Mean Aggregator
FC-24 + ReLU
FC-24 + ReLU
FC-24 + ReLU
FC
FC-24 + ReLU
FC + ReLU
FC
State M 4
State
State
M x 4
M x 2
(a) Mean update function
(b) Variance update function
(c) Observation mean function
History
Adjacency M x M
State
Observation
Adjacency
M x 30
M x
M x 2
M
FC-30 Tanh
FC-30 + Tanh
Mcan Aggregator
Mean Aggregator
FC-64 + Tanh
FC-64 + Tanh
FC
FC+ Exp
FC
FC+ Exp
Mean
Covariance
Mean
Covariance
M x V x
X
x 4
M x 4
(d) Embedding function
(e) Approximate posterior
Figure 8: Architectures without map information. For fully connected layers we give the number of output neurons.
The embedding function receives the history of all M agents. This history seconds long with time step of 0.2 seconds and consists of two dimensional coordinates. After flattening, the input is vector of size 30. The output of the embedding function is GMM with V mixture components. We use the mean aggregator in the embedding, mean and variance update function. The mean aggregator calculates the message lo agent m accordingly to Eq. 22 After the message of is calculated, is concatenated with the state x!". Mean and variance update functions are neural networks, which conduct at each prediction step one round of message passing and then caleulate the output. Our emissions model uses a neural network for the mean function g(x)) and a constant vector for Q(x:). The emission model maps the state of each agent back to the observed space and does not depend interactions. The approximate posterior is used
31

---

Published in Transactions on Machine Learning Research (03/2023)
in Sec. 5.1 for the baselines that trained the ELBO MCO. In contrast maximizing the PLL does not necessitate approximate posterior. The approximate posterior network takes the previous state and the current observation in order to approximate the latent distribution at the current time step.
H.2 Network Architectures with Map Information
These neural net architectures have been used in our experiments in Sec. 5.4 It closely follows the archi- tectures in Sec. H.1 We add an additional neural net, which encodes the masked map of size 500x300 into a flattened 256-dimensional vector. This map embeediing is used as an additional input to the embedding function.
State
Adjacency M M
State
Adjacency M x M
State
M x 16
M x 16
M 16
Mean Aggregator
Mean Aggregator
FC-50 ReLU
FC-50 + ReLU
FC-50 + ReLU
FC
FC-50 + ReLU
FC + ReLU
FC
State
State
State
M 16
M x 16
M 2
(a) Mean update function History+ Map
(b) Variance update function
(c) Observation mean function
Adjacency M M
Masked Map
M 30+ 256
500 x 300
FC-512 + Tanh
Con-
FC-512 + Tanh
Conv-6-12-ReelU+uMMaaaP
Mean Aggregator
Mon--12-44-RLLU-AAA
FC-512 + Tanh
FC
FC
FC+ Exp
Mean
Covariance
Encoded Map
M X x 16
x V x 16
256
(d) Embedding function
(e) Map encoder
Figure 9: Architectures with map information. For fully connected layers we give the number of output neurons. For convolutional layers we give the number input channels and output channels, We use a kernel size of 3 and stride 2. For MaxPool layers we use kernel size of and stride 2.
32


Reviews:
[{"review_id": "ae0lvBEs50", "reviewer": "Reviewer_CaF5", "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\nI love the simplicity and motivation behind the method! Moreover, I think the writing is top-notch and does a fantastic job describing the many working pieces of the proposed method. Moreover, the ablation method is **very** extensive and demonstrates not only the importance of individual pieces (non-linear vs linear dynamics, GNN vs no agent interaction, deterministic vs MC) but also provides an extensive evaluation of several different approximations for the covariance. \n\n# Weaknesses\n\nIn this section, I will discuss weaknesses, questions, and comments. Firstly, it seems at each time instant, only one round of message passing is done, while to my understanding, multiple rounds of message passing can be done at each time point. Is there a reason this was done? It would be nice if this were discussed. \n\nNext, on page 5 the authors state: *\"Furthermore, we note that although the transition noise factorizes across agents, our model is capable of modeling correlations between agents since the mean and the variance depend not only on the state of the m-th agent, but also on the states of all neighboring agents. In consequence, after $a$ aggregation steps, our model accounts for correlations between agent m and agent m′ provided that they are connected by a path that is at most a steps long.\"*. This isn't unique to the proposed approach but rather is a property of using transition functions for the mean and variance: for instance, if I used \n$$ p(x_t \\mid x_{t-1}) = \\mathcal{N}(x_t \\mid Ax_{t-1}, \\textrm{diag}\\left(\\exp( Bx_{t-1})) \\right) $$\nthen this would allow correlations between agents to propagate. \n\nLastly, there are no training details whatsoever for any of the experiments, making the paper nearly impossible to reproduce. Neither in the main paper or the appendix, do I see any details on batch size, learning rate, optimizer, etc. ", "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": "", "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 authors introduce Graph Deep State Space Models (GDSSM) to learn state space models over $M$ interacting agents. The contributions of this work are twofold: 1) the combination of classical deep-state space models with graph neural networks and 2) a deterministic training scheme based on moment matching, combined with a multi-modal prior, allows for efficient training of the model while allowing for expressivity.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "1. Eq. 14 is only true if you model the transition dynamics is $x_t = x_{t-1} + f(x_{t-1})$ e.g, $f$ is learning the velocity. As a simple counter-example, if $x_t = Ax_{t-1}$ and $x_t \\sim \\mathcal{N}(\\mu, \\Sigma)$ then $\\mathbb{E}[x_t] = A\\mu$ while if $x_t = x_t + Ax_{t-1}$ then $\\mathbb{E}[x_t] = \\mu + A\\mu$. Please update the previous state space model equations to reflect this. \n\n2. For completeness, it would be nice if the authors reviewed how to propagate means and covariances through the common nonlinearities, or at least for the non-linearities used in the paper.\n\n3. Some key references are missing. Mainly, there have been previous works on training dynamics models on multi-step ahead predictive log likelihood: https://arxiv.org/abs/1811.04551 and https://arxiv.org/abs/2106.06064.\n\n4. Training details need to be added to make this work reproducible. \n\n5. In section 5.2.1, the authors state: *\"There exists a large body of prior work [...] which applies GNNs for traffic forecasting on the NGSIM\ndataset and provide quantification of predictive uncertainty\" but none of the methods compared against use GNNs! For completeness, it would be nice if the authors also compared against a previously proposed method or, at the very at least, include GNN variants of the methods they did compare against, e.g., the output of the GNN can be fed into an LSTM as input.", "broader_impact_concerns": "There are no concerns on the ethical implications.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9ponOJLt3t", "reviewer": "Reviewer_ZrYM", "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\n* The problem is important and well-founded.\n* The proposed method explores an interesting direction.\n\n## Major Weaknesses\n\nUnfortunately, overall the paper is unconvincing because of a lack of clarity along with questionable methods and evaluation. The major issues I see are as follows:\n\n### “For many prediction tasks, the quantity of interest is the predictive log-likelihood”\n\nPlease provide a citation (other than Look et al.) and further explanation for this. Why would it be desirable to only match the marginals of the data? The structure in the full joint distribution is extremely important -- without it, trajectories generated by the model may lack coherent temporal structure and instead only match the marginals at each timestep.\n\nWhy is the likelihood of the initial observation $y_0$ dropped from the ‘predictive log likelihood’?\n\nThe scheduled sampling paper only supports the idea that it is important to have matching train loss and test evaluation. The text of your paper seems to suggest it supports more -- please make it clear which claim the scheduled sampling paper is used to support.\n\nFinally, you claim that the predictive log likelihood is the quantity of interest, but then report negative log likelihood in evaluation. As defined in your paper, the negative log likelihood and predictive log likelihoods are quite different. After the emphasis on the PLL in the main text, why only report the NLL? It makes the claims in the methods section seem suspect.\n\n### Monte Carlo baseline methods\n\nYou repeatedly claim MC-based methods are slow to converge and perform worse. You provide scant details about your MC-based baseline, however. More information about the training method used is need to determine if the baseline is correct and strong. The current state of the art methods for training nonlinear dynamical system models use sequential Monte Carlo, such as NASMC and VSMC/FIVO/AESMC.\n\nNASMC: Gu et al. 2015, https://arxiv.org/abs/1506.03338\n\nVSMC: Naesseth et al. 2017, https://arxiv.org/abs/1705.11140\n\nFIVO: Maddison et al. 2017, https://arxiv.org/abs/1705.09279\n\nAESMC: Le et al. 2017, https://arxiv.org/abs/1705.10306\n\n### Your moment matching scheme is possibly unfounded for mixtures of Gaussians\n\nYour model uses a mixture of Gaussians as a prior over its initial latent state, and claims to approximate nonlinear functions of this distribution (e.g. the application of neural network layers) using moment matching. This scheme is inadequately described and is possibly incorrect. It seems to be based broadly on the assumption that moment matching can be applied to mixtures of Gaussians by applying it individually to each component. Fitting a mixture of Gaussians via moment matching requires a large number of moments, not just two as used in this work. With a maximum entropy assumption (very common) the distribution that matches two moments is the Gaussian, not a mixture of Gaussians. This casts doubt on this work's approximation of each latent state using a mixture of Gaussians.\n\nPlease clarify how your moment matching scheme fits a mixture of Gaussians at each step.\n\n### General issues with evaluation\n\nYour evaluation metrics are unclear, and possibly incorrect. You claim to evaluate the negative log likelihood of your models, however the negative log likelihood is intractable for nonlinear latent variable models as stated in the main text of your paper. Details on your approximation scheme for the negative log likelihood should be included. A good method for this setting is SMC with a large number of particles.\n\nSimilarly, you should integrate out the latents in your model to compute the RMSE, but you are unclear on this point. You state in tables 1 and 2, “For the case of predictors with multiple modes we calculate the RMSE for each mode and report the lowest”. This methodology is incorrect, and could largely explain the gains claimed by incorporating multimodal state priors.\n\n## Minor Weaknesses\n\n### Incorrect and overly-broad claims.\n\nMany statements made in the introduction are overly broad or incorrect. Please narrow them or provide citations. For example:\n\n**“Deterministic systems, such as complex physical simulators”**\n\nMany simulators are not deterministic and include a random component.\n\n**“There are two established model families that can account for model uncertainty in dynamical systems\"**\n\nThis is overly broad. Stochastic dynamical systems research is incredibly diverse, with over a century of history. Furthermore, the two classes a) and b) that follow are unclear and seem to overlap. It would be good to clarify these categories and narrow your claims.\n\n**“However, propagating a distribution through a non-linear recurrent system, as it is required in this case, cannot be performed in closed-form, and existing methods build on Monte Carlo (MC) simulations.”**\n\nAgain, this is not true in general. Please narrow your claims, clarify what you mean, and provide citations.\n\n**“​​and well-calibrated predictions compared to state-of-the-art alternatives”**\n\nYou did not check the calibration of your model, so please do not make this claim.\n\n**\"In autonomous driving, the initial latent state can be estimated from historical information and is closely linked to the drivers’ intentions.”**\n\nThis is overly broad and must be untrue in general. The content of a latent state is model-specific.\n\n**“They come in handy for applications in which the true underlying dynamics are not known and must be estimated from data.”**\n\nThis is also true for linear state-space models.\n\n**“All of these approaches have in common that a good approximation to the smoothing distribution, i.e. p(x0, . . . , xT |y0, . . . , yT , I), is central for the algorithm to succeed.”**\n\nThis is not true.\n\n### Typos and grammatical errors\n\nThere are numerous typos and grammatical errors. For example, in the introduction:\n\n**“Allow to learn”**\n\n**“In order to account for this increasing model complexity, large network architectures are necessary that can restrict its usage on embedded systems with limited memory capacity”**\n\n**“Our model is capable of predicting multiple modes, which we efficiently approximate as a GMM and takes interaction into account by using GNNs.”**\n\nPlease fix those and others.", "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": "", "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 work proposes a model and inference scheme for modeling multiple interacting agents over time. The model is based on graph neural networks, which are used to explicitly encode constraints on which agents can interact. The inference scheme is based on bidirectional moment matching (BMM). The authors specialize BMM to the graph neural network setting and evaluate the model on traffic datasets, finding that their method improves over baselines.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "The main requested changes are:\n\n1. Clarify how your moment-matching scheme fits a mixture of Gaussians.\n2. Clarify your evaluation methodology, specifically how you estimated the negative log likelihood and why you picked the best mode for computing the RMSE.\n3. Clarify your baseline implementations.\n\nThese clarifications are critical in securing my recommendation for acceptance. However, clarification alone is insufficient. It is possible that clarifying these items will convince me that the author's method or evaluation are incorrect.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "HAyU5o5YhV", "reviewer": "Reviewer_vpQ1", "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": "Strong aspects:\n\n* There is a good motivation for the model, and for using graph neural networks to define transition dynamics for interacting agents\n* The model itself is a sensible extension of existing work, with an interesting deterministic inference approach\n* The experimental results are clearly presented, and demonstrate the utility of the method\n\nWeak aspects:\n\n* The contributions of the paper are not as well delineated from the existing work as they could be; in particular the organization of sections 3 and 4 is not ideal. Section 3 is introducing particular modeling choices (e.g., parameterizations of mean functions and variance functions; choice of GMM as prior over initial state), and then Section 4 covers a mix of background material (e.g. in section 4.2) as well as contributions (e.g. the moment computations for graph networks in section 4.4). This could be much easier to follow if there were a clear separation between (a) existing work, (b) modeling assumptions in this paper, (c) contributions required to apply existing work / inference algorithms to this model.\n* I'm afraid I do not understand the way in which the parameters of the GNNs itself are learned in this model. While much of section 4 describes the inference algorithm, there is not much information on the learning algorithm for the network parameters. Due to the approximations made in the moment propagation steps, it doesn't seem to me that the final log likelihood approximation (e.g. as computed in Alg. 1) would be appropriate for directly running gradient descent steps for updating the GNN parameters. (Please correct me if I am wrong!) I also don't see this explicitly stated. Presumably updating the GNN parameters either requires marginalizing out x, or else sampling / imputing values of x, and then taking gradient descent steps?\n* I am not familiar with this dataset and would like reassurance that this is how the data has been used in previous work. While the experiments do show that the model is successful at the chosen task, it seems like has been far simplified in terms of its generalization goals: essentially, restricting to a single roundabout (so no notion of generalization to different locations), and to 8 second segments at 5hz (i.e. with 15 frames as context, predicting an additional 25 frames). Does this match the sort of predictive task undertaken in other work looking at this data?", "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": "", "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 introduces a state space model where the transition dynamics are parameterized by a graph neural network. This is particularly relevant in multi-agent systems, where such a model is capable of learning complex interactions between agents. An inference algorithm is proposed, which adapts an assumed density approach known as bidimensional moment matching to the particular structure of the model, which requires propagating moments through a series of nonlinearities followed by an aggregation step.\n\nThe proposed model is evaluated on a dataset of vehicular traffic in a roundabout, with ablation studies (some of which reduce to previous work) successfully demonstrating the utility of the approach.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "The proposed changes, which would make me happy to raise my recommendation, essentially track the three points in the \"weaknesses\" section above.\n\n1. Please re-organize sections 3 and 4 to more clearly separate the model, and the paper's methodological contributions, from the background material on inference and BMM\n2. Please make explicit the objective that is used for learning the GMM parameters, and how that is handled algorithmically\n3. Please clarify the experimental setting, in the context of other work on this dataset\n\n", "broader_impact_concerns": "No concerns regarding broader impact", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This article addresses an interesting topic, how to model interacting stochastic systems and how to make inference for such systems. The paper has now been reviewed by three reviewers: two recommend accept (ZrYM  and CaF5) and one leans towards acceptance (vpQ1). The authors have updated substantially the paper at the revision stage and all the reviewers are now happy with the re-organized manuscript and the writing. They agree that the problem addressed is important and that the proposed methodology is principled and explores an interesting direction. Initial questions about some of the technical aspects of the methodology (moment matching) have been addressed satisfactorily by the authors and all reviewers do believe the work is technically sound and well evaluated.

I thus believe the material should be published as it is. The authors should make their code accessible for reproducibilty as promised in their response to reviewers.
Decision: Accept (poster)