{'fraction_non_alphanumeric': 0.034834765557413924, 'fraction_numerical': 0.014995097756502682, 'mean_word_length': 4.7510779436152575, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Graph Neural Ordinary Differential Equations

We introduce the framework of continuous--depth graph neural networks (GNNs). Graph neural ordinary differential equations (GDEs) are formalized as the counterpart to GNNs where the input-output relationship is determined by a continuum of GNN layers, blending discrete topological structures and differential equations. The proposed framework is shown to be compatible with various static and autoregressive GNN models. Results prove general effectiveness of GDEs: in static settings they offer computational advantages by incorporating numerical methods in their forward pass; in dynamic settings, on the other hand, they are shown to improve performance by exploiting the geometry of the underlying dynamics.


Introduction
Introducing appropriate inductive biases on deep learning models is a well-known approach to improving sample efficiency and generalization performance . 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;. Figure 1: Graph neural ordinary differential equations (GDEs) model vector fields defined on graphs, both in cases when the structure is fixed or changes in time, via a continuum of graph neural network (GNN) layers.
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 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;. 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.
Blending graphs and differential equations 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.
Sequences of graphs 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.
GDEs as general-purpose models 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.

General Framework
Definition of GDE Without any loss of generality, the inter-layer dynamics of a GNN node feature matrix can be represented in the form: where X e ∈ R n×h is an embedding of X 1 , F G is a matrix-valued nonlinear function conditioned on graph G and Θ(s) ∈ R p is the tensor of trainable parameters of the s-th layer. Note that the explicit dependence on s of the dynamics is justified in some graph architectures, such as diffusion graph convolutions (Atwood and Towsley, 2016). A graph neural differential ordinary equation (GDE) is defined as the following Cauchy problem: where F G : S × R n×h × R p → R n×h is a depth-varying vector field defined on graph G.
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.
Well-posedness Let S := [0, 1]. Under mild conditions on F, namely Lipsichitz continuity with respect to H and uniform continuity with respect to s, for each initial condition (GDE embedded input) X e , the ODE in (1) admits a unique solution H(s) defined in the whole S. Thus there is a mapping Ψ from R n×h to the space of absolutely continuous functions S → R n×h such that H := Ψ(X e ) satisfies the ODE in (1). This implies the the output Y of the GDE satisfies Y = Ψ(X e )(1).
Symbolically, the output of the GDE is obtained by the following Note that applying an output layer or network to Y before passing it to downstream applications is generally beneficial.
Integration domain We restrict the integration interval to S = [0, 1], given that any other integration time can be considered a rescaled version of S. Following  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 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).
GDE training GDEs can be trained with a variety of methods. Standard backpropagation through the computational graph, adjoint sensitivity method (Pontryagin et al., 1962) for O(1) memory efficiency , 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.

Static Models
Graph convolution differential equations 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: where L G ∈ R n×n is the graph Laplacian and σ is as nonlinear activation function. We denote with C G the graph convolution operator, i.e. C G H(s) = L G H(s)Θ(s). A general formulation of the continuous counterpart of GCNs, graph convolution differential equation (GCDE), is therefore obtained by defining F G as a multilayer convolution, i.e.
Note that the Laplacian L 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  are also compatible with the continuous-depth formulation.
Additional models and considerations 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.

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 ≡ t and can be adapted depending on the requirements. For example, given a time window ∆t, the prediction performed by a GDE assumes the form: regardless of the specific GDE architecture employed. Here, GDEs represent a natural model class for autoregressive modeling of sequences of graphs {G t } and seamlessly link to dynamical network theory. This line of reasoning naturally leads to an extension of classical spatiotemporal architectures in the form of hybrid dynamical systems (Van Der Schaft and Schumacher, 2000;Goebel et al., 2009), i.e., systems characterized by interacting continuous and discretetime dynamics. Let (K, >), (T , >) be linearly ordered sets; namely, K ⊂ N \{0} and T is a set of time instants, T := {t k } k∈K . We suppose to be given a state-graph data stream which is a sequence in the form {(X t , G t )} t∈T . Let us also define a hybrid time domain as the set I := k∈K ([t k , t k+1 ], k) and a hybrid arc on I as a function Φ such that for each k ∈ K, t → Φ(t, k) is absolutely continuous in {t : (t, j) ∈ dom Φ}. Our aim is to build a continuous model predicting, at each t k ∈ T , the value of X t k+1 , given (X t , G t ). The core idea is to have a GDE smoothly steering the latent node features between two time instants and then apply some discrete operator, resulting in a "jump" of H which is then processed by an output layer. Solutions of the proposed continuous spatio-temporal model are therefore hybrid arcs. Autoregressive GDEs The solution of a general autoregressive GDE model can be symbolically represented by: where F, G, K are GNN-like operators or general neural network layers 2 and H + represents the value of 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 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 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 G with other common recurrent modules, such as vanilla RNNs or LSTMs (Hochreiter and Schmidhuber, 1997). It should be noted that the operators F, G, 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.
2 More 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. Figure 3: Node embedding trajectories defined by a forward pass of GCDE-dpr5 on Cora, Citeseer and Pubmed. Color differentiates between node classes. 4.1 Transductive Node Classification Experimental setup 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 · 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.
Models and baselines 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  PyTorch package to solve and backpropagate through the ODE solver.
Continuous-depth models in static tasks 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  compared to equivalent discrete models. While there is no reason to expect performance improvements solely from the transition away from discrete architectures, continuousdepth 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.
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. Resilience to integration time For each integration time S ∈ [1, 5, 10], 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).

Multi-Agent Trajectory Extrapolation
Experimental setup 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 A t , the particles would follow a periodic motion. The adjaciency matrix A t is computed along the trajectory as: where x i (t) is the position of node i at time t. Therefore, A t results to be symmetric, A t = A t 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 · 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 meansquared-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 Φ 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 − ∆t i.eŶ t+∆t = φ(Ŷ t ) for a certain number of extrapolation steps, after which the model is fed the actual nominal state 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.
Models and baselines 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. Results 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.

Traffic Forecasting
Experimental setup 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 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 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.

Models and baselines
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.  Results 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)  and ODEs to track Hamiltonian functions, whereas (Deng et al., 2019) introduces a GNN version of continuous normalizing flows 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
Unknown or dynamic topology 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 S(s, H) is computed through some attentive mechanism, standard convolution GDEs can be replaced byḢ = σ (SHΘ) .
Addition or removal of nodes 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 Θ 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.

A Graph Neural Ordinary Differential Equations
Notation Let N be the set of natural numbers and R the set of reals. Scalars are indicated as lowercase letters, vectors as bold lowercase, matrices and tensors as bold uppercase and sets with calligraphic letters. Indices of arrays and matrices are reported as superscripts in round brackets.
Let V be a finite set with | V | = n whose element are called nodes and let E be a finite set of tuples of V elements. Its elements are called edges and are such that ∀e ij ∈ E, e ij = (v i , v j ) and v i , v j ∈ V. A graph G is defined as the collection of nodes and edges, i.e. G := (V, E). The adjeciency matrix A ∈ R n×n of a graph is defined as If G is an attributed graph, the feature vector of each v ∈ V is x v ∈ R nx . All the feature vectors are collected in a matrix X ∈ R n×nx . Note that often, the features of graphs exhibit temporal dependency, i.e. X := X t .

A.1 General Static Formulation
For clarity and as an easily accessible reference, we include below a general formulation table for the static case

Graph Neural Ordinary Differential Equations
Note that the general formulation provided in (4) can similarly serve as a reference for the spatiotemporal case.

A.2 Computational Overhead
As is the case for other models sharing the continuous-depth formulation , 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.
Fixed-step 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/ defines the number of steps necessary to cover [0, S] in fixed-steps of .
Adaptive-step 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
Message passing GDEs Let us consider a single node v ∈ V and define the set of neighbors of v as N (v) := {u ∈ V : (v, u) ∈ E ∨(u, v) ∈ E}. Message passing neural networks (MPNNs) perform a spatial-based convolution on the node v as where, in general, h v (0) = x v while u and m are functions with trainable parameters. For clarity of exposition, let u(x, y) := x + g(y) where g is the actual parametrized function. The (5) becomes and its continuous-depth counterpart, graph message passing differential equation (GMDE) is: Attention GDEs Graph attention networks (GATs) (Veličković et al., 2017) perform convolution on the node v as where α vu are attention coefficient which can be computed following (Veličković et al., 2017).

B Spatio-Temporal GDEs
We include a complete description of GCGRUs to clarify the model used in our experiments.

B.1 GCGRU Cell
Following GCGRUs (Zhao et al., 2018), we perform an instantaneous jump of H at each time t k using the next input features X t k . Let L G t k be the graph Laplacian of graph 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 where Θ xz , Θ hz , Θ xr , Θ hr , Θ xh , Θ hh are matrices of trainable parameters, σ is the standard sigmoid activation and 1 is all-ones matrix of suitable dimensions.  Table 3: General architecture for GCDEs on node classification tasks. GCDEs applied to different datasets share the same architecture. The vector field F is parameterized by two GCN layers. GCDEs-dopri5 shares the same structure without GDE-2 (GCN). Evolution of the adjacency matrix, 5 dt increments Figure 8: Snapshots of the evolution of adjacency matrix A t throughout the dynamics of the multiparticle system. Yellow indicates the presence of an edge and therefore a reciprocal force acting on the two bodies C.1 Node Classification

C Additional experimental details
Training hyperparameters 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.
Architectural details SoftPlus is used as activation for GDEs. Smooth activations have been observed to reduce stiffness  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 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 F are set to 0.9.

C.2 Multi-Agent System Dynamics
Dataset Let us consider a planar multi agent system with states x i (i = 1, . . . , n) and secondorder dynamics:ẍ and  The force f ij resembles the one of a spatial spring with drag interconnecting the two agents. The term −x i , is used instead to stabilize the trajectory and avoid the "explosion" of the phase-space. Note that f ij = −f ji . The adjaciency matrix A t is computed along a trajectory which indeed results to be symmetric, A t = A t and thus yields an undirected graph. Figure 8 visualizes an example trajectory of A t . The extrapolation is terminated after 5 steps and the nominal state is fed to the model.
We collect a single rollout with T = 5, dt = 1.95 · 10 −3 and n = 10. The particle radius is set to r = 1. Additional results 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
Dataset and metrics 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 y, andŷ ∈ R p is the set of vectorized target and prediction of models respectively. and · 1 denotes Hadamard division and the 1-norm of vector. where (·) 2 and √ · denotes the element-wise square and square root of the input vector, respectively. y t andŷ t denote the target and prediction vector.
Architectural details 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 F tasked with evolving the hidden features between arrival times. F is parametrized by 2 GCN layers, one with tanh activation and the second without activation. ReLU is used as the general activation function.
Training hyperparameters 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.
Additional results 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.

Station 200
Observation GRU GCGRU GCDE-GRU Figure 12: Traffic data prediction results of 50% undersampling. GCDE-GRUs are able to evolve the latents between timestamps and provide a more accurate fit.