Graph Neural Ordinary Differential Equations

Introduction

Introducing appropriate inductive biases on deep learning models is a well–known approach to improving sample efficiency and generalization performance (Battaglia et al., 2018). Graph neural networks (GNNs) represent a general computational framework for imposing such inductive biases when the problem structure can be encoded as a graph or in settings where prior knowledge about entities composing a target system can itself be described as a graph (Li et al., 2018c; Gasse et al., 2019; Sanchez-Gonzalez et al., 2018).

GNNs have shown remarkable results in various application areas such as node classification (Zhuang and Ma, 2018; Gallicchio and Micheli, 2019), graph classification (Yan et al., 2018) and forecasting (Li et al., 2017; Wu et al., 2019) as well as generative tasks (Li et al., 2018b; You et al., 2018). A different but equally important class of inductive biases is concerned with the type of temporal behavior of the systems from which the data is collected i.e., discrete or continuous dynamics. Although deep learning has traditionally been a field dominated by discrete models, recent advances propose a treatment of neural networks equipped with a continuum of layers (Haber and Ruthotto, 2017; Chen et al., 2018). This view allows a reformulation of the forward and backward pass as the solution of the initial value problem of an ordinary differential equation (ODE). Such approaches allow direct modeling of ODEs and can guide discovery of novel general purpose deep learning models.

In this work we propose the system–theoretic framework of graph neural ordinary differential equations (GDEs) by defining ODEs parametrized by GNNs. GDEs are designed to inherit the ability to impose relational inductive biases of GNNs while retaining the dynamical system perspective of continuous–depth models. We validate GDEs experimentally on a static semi–supervised node classification task as well as spatio–temporal forecasting tasks. GDEs are shown to outperform their discrete GNN analogues: the different sources of performance improvement are identified and analyzed separately for static and dynamic settings.

We extend the GDE framework to the spatio–temporal setting and formalize a general autoregressive GDE model as a hybrid dynamical system. The structure–dependent vector field learned by GDEs offers a data–driven approach to the modeling of dynamical networked systems (Lu and Chen, 2005; Andreasson et al., 2014), particularly when the governing equations are highly nonlinear and therefore challenging to approach with analytical methods. Autoregressive GDEs can adapt the prediction horizon by adjusting the integration interval of the ODE, allowing the model to track the evolution of the underlying system from irregular observations.

In general, no assumptions on the continuous nature of the data generating process are necessary in order for GDEs to be effective. Indeed, following recent work connecting different discretization schemes of ODEs (Lu et al., 2017) to previously known architectures such as FractalNets (Larsson et al., 2016), we show that GDEs can equivalently be utilized as high–performance general purpose models. In this setting, GDEs offer a grounded approach to the embedding of classic numerical schemes inside the forward pass of GNNs.

Graph Neural Ordinary Differential Equations

We begin by introducing the general formulation of GDEs.

Without any loss of generality, the inter–layer dynamics of a GNN node feature matrix can be represented in the form:

To reduce stiffness of learned vector fields, alleviating the computational burden of adaptive ODE solvers, the node features can be augmented in several ways (Dupont et al., 2019; Massaroli et al., 2020) by concatenating additional dimensions or prepending input layers to the GDE.

Symbolically, the output of the GDE is obtained by the following

Note that applying an output layer or network to $\mathbf{Y}$ before passing it to downstream applications is generally beneficial.

We restrict the integration interval to $\operatorname{\mathcal{S}}=$, given that any other integration time can be considered a rescaled version of $\operatorname{\mathcal{S}}$. Following (Chen et al., 2018) we use the number of function evaluations (NFE) of the numerical solver utilized to solve (1) as a proxy for model depth. In applications where $\operatorname{\mathcal{S}}$ acquires a specific meaning (i.e forecasting with irregular timestamps) the integration domain can be appropriately tuned to evolve GDE dynamics between arrival times (Rubanova et al., 2019) without assumptions on the functional form of the underlying vector field, as is the case for example with exponential decay in GRU–D (Che et al., 2018).

GDEs can be trained with a variety of methods. Standard backpropagation through the computational graph, adjoint sensitivity method (Pontryagin et al., 1962) for $\mathcal{O}(1)$ memory efficiency (Chen et al., 2018), or backpropagation through a relaxed spectral elements discretization (Quaglino et al., 2019). Numerical instability in the form of accumulating errors on the adjoint ODE during the backward pass of Neural ODEs has been observed in (Gholami et al., 2019). A proposed solution is a hybrid checkpointing–adjoint scheme commonly employed in scientific computing (Wang et al., 2009), where the adjoint trajectory is reset at predetermined points in order control the error dynamics.

Taxonomy of GDEs

In the following, we taxonomize GDEs models distinguishing them into static and spatio–temporal (autoregressive) variants.

Based on graph spectral theory (Shuman et al., 2013; Sandryhaila and Moura, 2013), the residual version of graph convolution network (GCN) (Kipf and Welling, 2016) layers are in the form:

Note that the Laplacian $\mathbf{L}_{\operatorname{\mathcal{G}}}$ can be computed in different ways, see e.g. (Bruna et al., 2013; Defferrard et al., 2016; Levie et al., 2018; Zhuang and Ma, 2018). Diffusion–type convolution layers (Li et al., 2017) are also compatible with the continuous–depth formulation.

We include additional derivation of continuous counterparts of common static GNN models such as graph attention networks (GAT) (Veličković et al., 2017) and general message passing GNNs as supplementary material.

2 Spatio–Temporal Models

For settings involving a temporal component (i.e., modeling dynamical systems), the depth domain of GDEs coincides with the time domain $s\equiv t$ and can be adapted depending on the requirements. For example, given a time window $\Delta t$, the prediction performed by a GDE assumes the form:

The solution of a general autoregressive GDE model can be symbolically represented by:

where $\mathbf{F},\mathbf{G},\mathbf{K}$ are GNN–like operators or general neural network layersMore formal definitions of the hybrid model in the form of hybrid inclusions can indeed be easily given. However, the technicalities involved are beyond the scope of this paper. and $\mathbf{H}^{+}$ represents the value of $\mathbf{H}$ after the discrete transition. The evolution of system (4) is indeed a sequence of hybrid arcs defined on a hybrid time domain. A graphical representation of the overall system is given by the hybrid automata as shown in Fig. 2. Compared to standard recurrent models which are only equipped with discrete jumps, system (4) incorporates a continuous flow of latent node features $\mathbf{H}$ between jumps. This feature of autoregressive GDEs allows them to track the evolution of dynamical systems from observations with irregular time steps. In the experiments we consider $\mathbf{G}$ to be a GRU cell (Cho et al., 2014), obtaining graph convolutional differential equation–GRU (GCDE–GRU). Alternatives such as GCDE–RNNs or GCDE–LSTMs can be similarly obtained by replacing $\mathbf{G}$ with other common recurrent modules, such as vanilla RNNs or LSTMs (Hochreiter and Schmidhuber, 1997). It should be noted that the operators $\mathbf{F},\mathbf{G},\mathbf{K}$ can themselves have multi–layer structure.

Experiments

We evaluate GDEs on a suite of different tasks. The experiments and their primary objectives are summarized below:

Semi–supervised node classification on static, standard benchmark datasets Cora, Citeseer, Pubmed (Sen et al., 2008). We investigate the usefulness of the proposed method in a static setting via an ablation analysis that directly compares GCNs and analogue GCDEs solved with fixed–step and adaptive solvers.

Trajectory extrapolation task on a synthetic multi–agent dynamical system. We compare Neural ODEs and GDEs, providing a motivating example for the introduction of additional biases in the form of second–order models (Yıldız et al., 2019; Massaroli et al., 2020).

Traffic forecasting on an undersampled version of PeMS (Yu et al., 2018) dataset. We measure the performance improvement obtained by a correct inductive bias on continuous dynamics and robustness to irregular timestamps.

The code will be open–sourced after the review phase and is included in the submission.

The first task involves performing semi–supervised node classification on static graphs collected from baseline datasets Cora, Pubmed and Citeseer (Sen et al., 2008). Main goal of these experiments is to perform an ablation study on the source of possible performance advantages of the GDE framework in settings that do not involve continuous dynamical systems. The $L_{2}$ weight penalty is set to $5\cdot 10^{-4}$ on Cora, Citeseer and $10^{-3}$ on Pubmed as a strong regularizer due to the small size of the training set (Monti et al., 2017). We report mean and standard deviation across $100$ training runs. Since our experimental setup follows (Kipf and Welling, 2016) to allow for a fair comparison, other baselines present in recent GNN literature can be directly compared with Table 1.

All convolution–based models are equipped with a latent dimension of $64$. We include results for best performing vanilla GCN baseline presented in (Veličković et al., 2017). To avoid flawed comparisons, we further evaluate an optimized version of GCN, GCN∗, sharing exact architecture, as well as training and validation hyperparameters with the GCDE models. We experimented with different number of layers for GCN∗: (2, 3, 4, 5, 6) and select $2$, since it achieves the best results. The performance of graph convolution differential equation (GCDE) is assessed with both a fixed-step solver Runge–Kutta (Runge, 1895; Kutta, 1901) as well as an adaptive–step solver, Dormand–Prince (Dormand and Prince, 1980). The resulting models are denoted as GCDE–rk4 and GCDE–dpr5, respectively. We utilize the torchdiffeq (Chen et al., 2018) PyTorch package to solve and backpropagate through the ODE solver.

Ensuring a low error solution to the ODE parametrized by the model with adaptive–step solvers does not offer particular advantages in image classification tasks (Chen et al., 2018) compared to equivalent discrete models. While there is no reason to expect performance improvements solely from the transition away from discrete architectures, continuous–depth allows for the embedding of numerical ODE solvers in the forward pass. Multi–step architectures have previously been linked to ODE solver schemes (Lu et al., 2017) and routinely outperform their single–step counterparts (Larsson et al., 2016; Lu et al., 2017). We investigate the performance gains by employing the GDE framework in static settings as a straightforward approach to the embedding of numerical ODE solvers in GNNs.

The variants of GCDEs solved with fixed–step schemes are shown to outperform or match GCN∗ across all datasets, with the margin of improvement being highest on Cora and Citeseer. Introducing GCDE–rk2 and GCDE–rk4 is observed to provide the most significant accuracy increases in more densely connected graphs or with larger training sets. In particular, GCDE–rk4 outperforming GCDE–rk2 indicates that, given equivalent network architectures, higher order ODE solvers are generally more effective, provided the graph is dense enough to benefit from the additional computation. Additionally, training GCDEs with adaptive step solvers naturally leads to deeper models than possible with vanilla GCNs, whose layer depth greatly reduces performance. However, the high number of function evaluation (NFEs) of GCDE–dpr5 necessary to stay within the ODE solver tolerances causes the model to overfit and therefore generalize poorly. We visualize the first two components of GCDE–dpr5 node embedding trajectories in Figure 3. The trajectories are divergent, suggesting a non–decreasing classification performance for GCDE models trained with longer integration times. We provide complete visualization of accuracy curves in Appendix C.

For each integration time $S\in$, we train $100$ GCDE-dpr5 models on Cora and report average metrics, along with $1$ standard deviation confidence intervals in Figure 4. GCDEs are shown to be resilient to changes in $S$; however, GCDEs with longer integration times require more training epochs to achieve comparable accuracy. This result suggests that, indeed, GDEs are immune to node oversmoothing (Oono and Suzuki, 2019).

2 Multi–Agent Trajectory Extrapolation

We evaluate GDEs and a collection of deep learning baselines on the task of extrapolating the dynamical behavior of a synthetic mechanical multi–particle system. Particles interact within a certain radius with a viscoelastic force. Outside the mutual interactions, captured by a time–varying adjacency matrix $\mathbf{A}_{t}$, the particles would follow a periodic motion. The adjaciency matrix $\mathbf{A}_{t}$ is computed along the trajectory as:

where $\mathbf{x}_{i}(t)$ is the position of node $i$ at time $t$. Therefore, $\mathbf{A}_{t}$ results to be symmetric, $\mathbf{A}_{t}=\mathbf{A}_{t}^{\top}$ and yields an undirected graph. The dataset is collected by integrating the system for $T=5s$ with a fixed step–size of $dt=1.95\cdot 10^{-3}$ and is split evenly into a training and test set. We consider $10$ particle systems. An example trajectory is shown in Figure 5. All models are optimized to minimize mean–squared–error (MSE) of 1–step predictions using Adam (Kingma and Ba, 2014) with constant learning rate $0.01$. We measure test mean average percentage error (MAPE) of model predictions in different extrapolation settings. Extrapolation steps denotes the number of predictions each model $\Phi$ has to perform without access to the nominal trajectory. This is achieved by recursively letting inputs at time $t$ be model predictions at time $t-\Delta t$ i.e $\hat{\mathbf{Y}}_{t+\Delta t}=\phi(\hat{\mathbf{Y}}_{t})$ for a certain number of extrapolation steps, after which the model is fed the actual nominal state $\mathbf{X}$ and the cycle is repeated until the end of the test trajectory. For a robust comparison, we report mean and standard deviation across 10 seeded training and evaluation runs. Additional experimental details, including the analytical formulation of the dynamical system, are provided as supplementary material.

As the vector field depends only on the state of the system, available in full during training, the baselines do not include recurrent modules. We consider the following models:

A 3–layer fully-connected neural network, referred to as Static. No assumption on the dynamics

A vanilla Neural ODE with the vector field parametrized by the same architecture as Static. ODE assumption on the dynamics.

A 3–layer convolution GDE, GCDE. Dynamics assumed to be determined by a blend of graphs and ODEs

A 3–layer, second–order (Yıldız et al., 2019; Massaroli et al., 2020) GCDE and referred to as GCDE-II. GCDE assumptions in addition to second–order ODE dynamics.

A grid hyperparameter search on number of layers, ODE solver tolerances and learning rate is performed to optimize Static and Neural ODEs. We use the same hyperparameters for GDEs.

Figure 6 shows the growth rate of test MAPE error as the number of extrapolation steps is increased. Static fails to extrapolate beyond the 1–step setting seen during training. Neural ODEs overfit spurious particle interaction terms and their error rapidly grows as the number of extrapolation steps is increased. GCDEs, on the other hand, are able to effectively leverage relational information to track the system: we provide complete visualization of extrapolation trajectory comparisons in the supplementary material. Lastly, GCDE-IIs outperform first–order GCDEs as their structure inherently possesses crucial information about the relative relationship of positions and velocities that is accurate with respect to the observed dynamical system.

3 Traffic Forecasting

We evaluate the effectiveness of autoregressive GDE models on forecasting tasks by performing a series of experiments on the established PeMS traffic dataset. We follow the setup of (Yu et al., 2018) in which a subsampled version of PeMS, PeMS7(M), is obtained via selection of 228 sensor stations and aggregation of their historical speed data into regular 5 minute frequency time series. We construct the adjacency matrix $\mathbf{A}$ by thresholding of the Euclidean distance between observation stations i.e. when two stations are closer than the threshold distance, an edge between them is included. The threshold is set to the 40${}^{\text{th}}$ percentile of the station distances distribution. To simulate a challenging environment with missing data and irregular timestamps, we undersample the time series by performing independent Bernoulli trials on each data point. Results for 3 increasingly challenging experimental setups are provided: undersampling with $30\%$, $50\%$ and $70\%$ of removal. In order to provide a robust evaluation of model performance in regimes with irregular data, the testing is repeated $20$ times per model, each with a different undersampled version of the test dataset. We collect root mean square error (RMSE) and MAPE. More details about the chosen metrics and data are included as supplementary material.

In order to measure performance gains obtained by GDEs in settings with data generated by continuous time systems, we employ a GCDE–GRU–dopri5 as well as its discrete counterpart GCGRU (Zhao et al., 2018). To contextualize the effectiveness of introducing graph representations, we include the performance of GRUs since they do not directly utilize structural information of the system in predicting outputs. Apart from GCDE–GRU, both baseline models have no innate mechanism for handling timestamp information. For a fair comparison, we include timestamp differences between consecutive samples and sine–encoded (Petneházi, 2019) absolute time information as additional features. All models receive an input sequence of $5$ graphs to perform the prediction.

Non–constant differences between timestamps result in a challenging forecasting task for a single model since the average prediction horizon changes drastically over the course of training and testing. Traffic systems are intrinsically dynamic and continuous in nature and, therefore, a model able to track continuous underlying dynamics is expected to offer improved performance. Since GCDE-GRUs and GCGRUs are designed to match in structure we can measure this performance increase from the results shown in Table 2. GCDE–GRUs outperform GCGRUs and GRUs in all undersampling regimes. Additional details and prediction visualizations are included in Appendix C.

Related work

There exists a concurrent line of work (Xhonneux et al., 2019) introducing a GNN variant evaluated on static node classification tasks where the output is the analytical solution of a linear ODE. Sanchez-Gonzalez et al. (2019) proposes using graph networks (GNs) (Battaglia et al., 2018) and ODEs to track Hamiltonian functions, whereas (Deng et al., 2019) introduces a GNN version of continuous normalizing flows (Chen et al., 2018; Grathwohl et al., 2018) for generative modeling, extending (Liu et al., 2019). Our goal is developing a unified system–theoretic framework for continuous–depth GNNs covering the main variants of static and spatio–temporal GNN models. Our work provides extensive experimental evaluations on both static as well as dynamic tasks with the primary aim of uncovering the sources of performance improvement of GDEs in each setting.

Discussion

Several lines of work concerned with learning the graph structure directly from data exist, either by inferring the adjacency matrix within a probabilistic framework (Kipf et al., 2018) or using a soft–attention (Vaswani et al., 2017) mechanism (Choi et al., 2017; Li et al., 2018a; Wu et al., 2019). In particular, the latter represents a commonly employed approach to the estimation of a dynamic adjacency matrix in spatio–temporal settings. Due to the algebraic nature of the relation between the attention operator and the node features, GDEs are compatible with its use inside the GNN layers parametrizing the vector field. Thus, if an optimal adaptive graph representation $\mathbf{S}(s,\mathbf{H})$ is computed through some attentive mechanism, standard convolution GDEs can be replaced by $\dot{\mathbf{H}}=\sigma\left(\mathbf{S}\mathbf{H}\bm{\Theta}\right).$

GDE variants operating on sequences of dynamically changing graphs can, without changes to the formulation, directly accommodate addition or removal of nodes as long as its number remains constant during the flows. In fact, the size of parameter matrix $\bm{\Theta}$ exclusively depends on the node feature dimension, resulting in resilience to a varying number of nodes.

Conclusion

In this work we introduce graph neural ordinary differential equations (GDE), the continuous–depth counterpart to graph neural networks (GNN) where the inputs are propagated through a continuum of GNN layers. The GDE formulation is general, as it can be adapted to include many static and autoregressive GNN models. GDEs are designed to offer a data–driven modeling approach for dynamical networks, whose dynamics are defined by a blend of discrete topological structures and differential equations. In sequential forecasting problems, GDEs can accommodate irregular timestamps and track the underlying continuous dynamics, whereas in static settings they offer computational advantages by allowing for the embedding of black–box numerical solvers in their forward pass. GDEs have been evaluated on both static and dynamic tasks and have been shown to outperform their discrete counterparts.

References

A.1 General Static Formulation

For clarity and as an easily accessible reference, we include below a general formulation table for the static case

A.2 Computational Overhead

As is the case for other models sharing the continuous–depth formulation (Chen et al., 2018), the computational overhead required by GDEs depends mainly by the numerical methods utilized to solve the differential equations. We can define two general cases for fixed–step and adaptive–step solvers.

In the case of fixed–step solvers of k–th order e.g Runge–Kutta–k (Runge, 1895), the time complexity is $O(nk)$ where $n:=S/\epsilon$ defines the number of steps necessary to cover $[0,S]$ in fixed–steps of $\epsilon$.

For general adaptive–step solvers, computational overhead ultimately depends on the error tolerances. While worst–case computation is not bounded (Dormand and Prince, 1980), a maximum number of steps can usually be set algorithmically.

A.3 Additional GDEs

Let us consider a single node $v\in\operatorname{\mathcal{V}}$ and define the set of neighbors of $v$ as $\operatorname{\mathcal{N}}(v):=\{u\in\operatorname{\mathcal{V}}~{}:~{}(v,u)\in\operatorname{\mathcal{E}}\lor(u,v)\in\operatorname{\mathcal{E}}\}$. Message passing neural networks (MPNNs) perform a spatial–based convolution on the node $v$ as

where, in general, $\mathbf{h}^{v}(0)=\mathbf{x}_{v}$ while $\mathbf{u}$ and $\mathbf{m}$ are functions with trainable parameters. For clarity of exposition, let $\mathbf{u}(\mathbf{x},\mathbf{y}):=\mathbf{x}+\mathbf{g}(\mathbf{y})$ where $\mathbf{g}$ is the actual parametrized function. The (5) becomes

and its continuous–depth counterpart, graph message passing differential equation (GMDE) is:

Graph attention networks (GATs) (Veličković et al., 2017) perform convolution on the node $v$ as

Similarly, to GCNs, a virtual skip connection can be introduced allowing us to define the graph attention differential equation (GADE):

where $\alpha_{vu}$ are attention coefficient which can be computed following (Veličković et al., 2017).

Appendix B Spatio–Temporal GDEs

We include a complete description of GCGRUs to clarify the model used in our experiments.

Following GCGRUs (Zhao et al., 2018), we perform an instantaneous jump of $\mathbf{H}$ at each time $t_{k}$ using the next input features $\mathbf{X}_{t_{k}}$. Let $\mathbf{L}_{\operatorname{\mathcal{G}}_{t_{k}}}$ be the graph Laplacian of graph $\operatorname{\mathcal{G}}_{t_{k}}$, which can computed in several ways (Bruna et al., 2013; Defferrard et al., 2016; Levie et al., 2018; Zhuang and Ma, 2018). Then, let

Finally, the post–jump node features are obtained as

Appendix C Additional experimental details

We carried out all experiments on a cluster of 4x12GB NVIDIA® Titan Xp GPUs and CUDA 10.1. The models were trained on GPU.

C.1 Node Classification

All models are trained for $2000$ epochs using Adam (Kingma and Ba, 2014) with learning rate $lr=10^{-3}$ on Cora, Citeseer and $lr=10^{-2}$ on Pubmed due to its training set size. The reported results are obtained by selecting the lowest validation loss model after convergence (i.e. in the epoch range $1000$ – $2000$). Test metrics are not utilized in any way during the experimental setup. For the experiments to test resilience to integration time changes, we set a higher learning rate for all models i.e. $lr=10^{-2}$ to reduce the number of epochs necessary to converge.

SoftPlus is used as activation for GDEs. Smooth activations have been observed to reduce stiffness (Chen et al., 2018) of the ODE and therefore the number of function evaluations (NFE) required for a solution that is within acceptable tolerances. All the other activation functions are rectified linear units (ReLU). The exact input and output dimensions for the GCDE architectures are reported in Table 3. The vector field $\mathbf{F}$ of GCDEs–rk2 and GCDEs–rk4 is parameterized by two GCN layers. GCDEs–dopri5 shares the same structure without GDE–2 (GCN). Input GCN layers are set to dropout $0.6$ whereas GCN layers parametrizing $\mathbf{F}$ are set to $0.9$.

C.2 Multi–Agent System Dynamics

Let us consider a planar multi agent system with states $\mathbf{x}_{i}$ ($i=1,\dots,n$) and second–order dynamics:

The force $\mathbf{f}_{ij}$ resembles the one of a spatial spring with drag interconnecting the two agents. The term $-\mathbf{x}_{i}$, is used instead to stabilize the trajectory and avoid the ”explosion” of the phase–space. Note that $\mathbf{f}_{ij}=-\mathbf{f}_{ji}$. The adjaciency matrix $\mathbf{A}_{t}$ is computed along a trajectory

which indeed results to be symmetric, $\mathbf{A}_{t}=\mathbf{A}_{t}^{\top}$ and thus yields an undirected graph. Figure 8 visualizes an example trajectory of $\mathbf{A}_{t}$.

We collect a single rollout with $T=5$, $dt=1.95\cdot 10^{-3}$ and $n=10$. The particle radius is set to $r=1$.

Node feature vectors are $4$ dimensional, corresponding to the dimension of the state, i.e. position and velocity. Neural ODEs and Static share an architecture made up of 3 fully–connected layers: $4n$, $8n$, $8n$, $4n$ where $n=10$ is the number of nodes. The last layer is linear. We evaluated different hidden layer dimensions: $8n$, $16n$, $32n$ and found $8n$ to be the most effective. Similarly, the architecture of first order GCDEs is composed of 3 GCN layers: $4$, $16$, $16$, $4$. Second–order GCDEs, on the other hand, are augmented by $4$ dimensions: $8$, $32$, $32$, $8$. We experimented with different ways of encoding the adjacency matrix $\mathbf{A}$ information into Neural ODEs and $Static$ but found that in all cases it lead to worse performance.

We report in Figure 9 test extrapolation predictions of $5$ steps for GDEs and the various baselines. Neural ODEs fail to track the system, particularly in regions of the state space where interaction forces strongly affect the dynamics. GDEs, on the other hand, closely track both positions and velocities of the particles.

C.3 Traffic Forecasting

The timestamp differences between consecutive graphs in the sequence varies due to undersampling. The distribution of timestamp deltas (5 minute units) for the three different experiment setups (30%, 50%, 70% undersampling) is shown in Figure 11.

As a result, GRU takes 230 dimensional vector inputs (228 sensor observations + 2 additional features) at each sequence step. Both GCGRU and GCDE–GRU graph inputs with and 3 dimensional node features (observation + 2 additional feature). The additional time features are excluded for the loss computations. We include MAPE and RMSE test measurements, defined as follows:

where $(\cdot)^{2}$ and $\sqrt{\cdot}$ denotes the element-wise square and square root of the input vector, respectively. $\mathbf{y}_{t}\ \text{and}\ \hat{\mathbf{y}}_{t}$ denote the target and prediction vector.

We employed two baseline models for contextualizing the importance of key components of GCDE–GRU. GRUs architectures are equipped with 1 GRU layer with hidden dimension 50 and a 2 layer fully–connected head to map latents to predictions. GCGRUs employ a GCGRU layer with 46 hidden dimension and a 2 layer fully–connected head. Lastly, GCDE–GRU shares the same architecture GCGRU with the addition of the flow $\mathbf{F}$ tasked with evolving the hidden features between arrival times. $\mathbf{F}$ is parametrized by 2 GCN layers, one with tanh activation and the second without activation. ReLU is used as the general activation function.

All models are trained for 40 epochs using Adam(Kingma and Ba, 2014) with $lr=10^{-2}$. We schedule $lr$ by using cosine annealing method (Loshchilov and Hutter, 2016) with $T_{0}=10$. The optimization is carried out by minimizing the mean square error (MSE) loss between predictions and corresponding targets.

Training curves of the models are presented in the Fig 10. All of models achieved nearly 13 in RMSE during training and fit the dataset. However, due to the lack of dedicated spatial modeling modules, GRUs were unable to generalize to the test set and resulted in a mean value prediction.