Temporal Graph ODEs for Irregularly-Sampled Time Series

1 Introduction

Representation learning for graphs has been gaining increasing attention over recent years. Such popularity often builds on the fact that complex phenomena are frequently understood as systems of interacting entities described as a graph. For such a reason, learning on graph-structured data through Deep Graph Networks (DGNs) [Bacciu et al., 2020; Wu et al., 2021] has been adopted for solving problems in a variety of fields, such as biology, social science, and sensor networks [Gilmer et al., 2017; Zitnik et al., 2018; Monti et al., 2019; Derrow-Pinion et al., 2021; Gravina et al., 2022; Bacciu et al., 2024].

Graph-based processing methods turned out to be extremely effective in processing spatio-temporal data too [Li et al., 2018; Wu et al., 2019b; Bai et al., 2021; Chen et al., 2022; Zambon and Alippi, 2022; Marisca et al., 2022; Jiang and Luo, 2022; Errica et al., 2023; Cini et al., 2023a]. Such a scenario is a setting where the temporal modeling accounts for functional dependencies - assimilated, in a broad sense, as spatial relationships - existing among the interacting entities. Real-world complex problems described as temporal graphs, e.g., those associated with social interactions, call for novel methods that can move beyond the common assumptions found in most of the methods proposed until now. Indeed, such problems require dealing with mutable relational information, irregularly and severely under-sampled data.

∗Corresponding author.

Some recent works propose to model input-output data relations as a continuous dynamic described by a learnable ordinary differential equation (ODE), instead of discrete sequences of layers commonly used in deep learning. Neural ODE-based approaches have been exploited to model non-temporal data, including message-passing functions for learning node-level embeddings [Poli et al., 2019; Chamberlain et al., 2021; Eliasof et al., 2021; Rusch et al., 2022; Gravina et al., 2023]. Notably, relying on ODEs has shown promising for modeling complex temporal patterns from irregularly and sparsely sampled data [Chen et al., 2018; Rubanova et al., 2019; Kidger et al., 2020].

In this paper, we formulate Temporal Graph Ordinary Differential Equation (TG-ODE), a general continuous-time modeling framework for temporal graphs that encompasses some methods from the literature as its specific instances, and demonstrate that TG-ODE can be an effective design choice to operate with irregularly and sparsely sampled observations. TG-ODE is designed through the lens of ODEs for effective learning of irregularly sampled temporal graphs. With TGODE, a differential equation is learned directly from data to solve a downstream task and predictions are trajectories obtained by numerical integration of the learned ODE.

The key contributions of this work can be summarized as follows:

• (i) we introduce TG-ODE, a general modeling framework suited for handling irregularly sampled temporal graphs;
• (ii) we introduce new benchmarks of synthetic and realworld scenarios for evaluating forecasting models on irregularly sampled temporal graphs; and
• (iii) we conduct extensive experiments to demonstrate the benefits of our method and show that TG-ODE outperforms state-of-the-art DGNs on all benchmarks.

Finally, we stress that, other than the outstanding empirical performance achieved by even simple TG-ODE instances, the framework allows us to reinterpret many state-of-the-art DGNs as a discretized solution of an ODE, thus facilitating their extension to handle graph streams with irregular sampling.

2 Problem Statement

We consider a dynamical system of interacting entities u ∈ V, which we refer to as nodes, that is described by a Cauchy problem defined on an ODE of the form

with initial condition X(0) = X0. System state X(t) = {xu(t) : u ∈ V(t)} is a function of time t and collects the node-level states xu(t) ∈ Rdx associated with each node u ∈ V(t). The node set V(t) is allowed to vary over time. We denote the dynamic set of edges encoding the node relations as E(t) ⊆ V(t)×V(t). E(t) = {euv(t) : (u, v) ∈ E} is a set of edge-level attributes evu(t) ∈ Rde defining the type or strength of the interaction between nodes u, v ∈ V(t). The system can also be driven by vector zu(t) ∈ Rc accounting for exogenous variables relevant to the problem at hand, such as weather conditions, hour of the day, or day of the week. Accordingly, for all u ∈ V(t), we write

to emphasize the local dependencies of node state xu(t) at a time t from its neighboring nodes v ∈ Nu(t) = {v ∈ V(t) : (v, u) ∈ E(t)} at the corresponding time.

We express any solution of ODE (1) as the temporal graph

defined for t ≥ 0. However, we assume to observe system (1) only as a (discrete) sequence of snapshot graphs

that arrive at irregular timestamps, i.e., the sampling is not uniform and, in general, ti - ti-1 ̸= ti+1 - ti. Each snapshot Gt = (Vt, Et,Xt,Et) corresponds to an observation of the system state at a specific timestamp t ∈ R. Edge set Et contains the functional relations between the nodes u, v ∈ Vt at time t (thus, Et can vary over time), while Xt = {xt,v : v ∈ Vt} gathers the observed node states, with xt,v the state of node v associated to the observed graph snapshot at time t. Edge attributes, if present, are denoted as eu,v ∈ Et. Figure 1 visually summarizes this concept.

In this paper, we address the problem of learning a model of the differential equation underlying the observed data, which is subsequently exploited to provide estimates of unobserved system's node states and make forecasts. To ease readability, in the following we drop the time variable t.

3 Temporal Graph Ordinary Differential Equation

To learn the function F in (2), we consider a family of models

parameterized by vector θ, and optimized so that the solution ˆx of the differential equation

minimizes the discrepancy with the observed sequence of graphs in (4). We follow the message-passing paradigm [Gilmer et al., 2017] and instantiate (6) as

with message function ϕM, aggregation operator ρ (e.g., the mean), and update function ϕU . Note that functions ϕU , ϕM, and ρ have learnable parameters θ, which are shared across all nodes and time steps. An example of a message passing operator is when ϕM(xu, xv, evu) = xv, ρ is the mean, and ϕU is a dense feed-forward network of its inputs. We refer to the above framework in (7) as Temporal Graph Ordinary Differential Equation (TG-ODE).

We observe that, since our framework relies on ODEs, it can naturally deal with snapshots that arrive at an arbitrary time. Indeed, the original Cauchy problem can be divided into multiple sub-problems, one per snapshot in the temporal graph. Here, the i-th sub-problem is defined for all u ∈ Vt as

in the time span between the two consecutive timestamps, i.e., t ∈ [ti-1, ti], where ψ is a function that combines the i-th observed state of the node u related to the snapshot graph Gti-1 in (4) (i.e., xti-1,u) and the prediction ˆxu(ti-1) obtained by solving (8) at the previous step. When given, we consider the true - potentially variable - topology E(t) to define the neighborhoods for t ∈ [ti-1, ti], otherwise, we set E(t) ≡ Eti-1 for every t, i.e., equal to the last observed topology associated with Gti-1 . Accordingly, we optimize θ in order to minimize the mean of some loss L,

where prediction ˆxu(ti) at time ti is obtained by solving (8).

We observe that for most ODEs it is not possible to compute analytical solutions. For such a reason, it is common practice to resort to numerical approximations that leverage discretization strategies (e.g., forward Euler's method). In this case, the solution is computed through iterative applications of the method over a discrete set of points in the time interval and, as observed in [Chen et al., 2018; Haber and Ruthotto, 2017], the process can be assimilated to that of a Recurrent Neural Network (RNN). For simplicity, here we employ the forward Euler's method according to which a solution to (8) is obtained by the recursion

starting from initial condition x0u = xu(0) = ψ(xti-1,u, ˆxu(ti-1)) and is reiterated until ϵℓ ≥ ti - ti-1. In (10), ϵ ≪ ti - ti-i is the step size, while ℓ indicates the generic iteration step, and tℓ = ti-1 + ϵℓ. Finally, a solution ˆx(t) to (8) in the interval [ti-1, ti] is provided by the discretization ˆxu(ti-1 + ϵℓ) = xℓu , for all ℓ, and interpolated elsewhere. The process is visually summarized at the bottom of Figure 2.

We acknowledge that not all resulting ODEs allow unique solutions and yield numerical stable problems. Generally, numerical stability is associated with the ODE to solve rather than the input data. Thus, a proper design of the considered family of ODEs in (7) can prevent stability issues. Indeed, the solution of a Cauchy problem exists and is unique if the differential equation is uniformly Lipschitz continuous in its input and continuous in t, as states in the Picard-Lindel¨of theorem [Coddington and Levinson, 1955]. Thus, different implementations of (7) should address the continuous behavior of the differential equation. We note that this theorem holds for our model if the underlying neural network has finite weights and uses Lipshitz non-linearities, e.g., the tanh.

By (10) and the generality of the message passing in (7), we observe that TG-ODE allows us to cast basically any standard DGN through the lens of an ODE for temporal graphs with irregular timestamps. Secondly, we stress that, even though TG-ODE is solved here by means of the forward Euler's method, other discretization methods can still be utilized. To conclude, our framework can be implemented using the aggregation function that is most suitable for the given task and the discretization method that best fits the computational resources and problem at hand, such as [Veliˇckovi´c et al., 2018; Hu et al., 2020; Choi et al., 2023; Eliasof et al., 2024a]. As a demonstration of this, in Section 5

we explore the neighborhood aggregation scheme proposed in [Du et al., 2017]. Thus, (7) can be reformulated as

where σ is an activation function, K is the number of hops in the neighborhood, θk is the k-th weight matrix, N k u is the k-hop neighborhood of u, and α(k) u,v is a normalization term. For instance, α(k) u,v = ( ˆ d(k) v ˆ d(k) u )-1/2 weighs according to the degrees ˆ d(k) v and ˆ d(k) u of nodes v and u in the k-hop graph; other choices can include edge attributes as well. We note that θk is a parameter specific to the k-hop neighborhood of node u. Thus, it allows the model to learn different transformation patterns at different distances from the considered node u.

4 Related works

Deep Graph Networks for static graphs. In the static graph domain, most DGNs can be generalized by the concepts introduced by MPNN [Gilmer et al., 2017], which is designed to capture and propagate information between nodes in a graph through a message-passing mechanism. Thus, each node updates its state by exchanging information with its neighboring nodes. MPNN can be formulated as

where ϕU and ϕM are respectively the update and message functions, which are responsible for computing neighbor states and updating the node's state accordingly. Most of the state-of-the-art DGNs can be derived by (12). Indeed, the definition of the message and update function allows implementing DGNs with different properties [Kipf and Welling, 2017; Veliˇckovi´c et al., 2018; Hamilton et al., 2017; Xu et al., 2019; Defferrard et al., 2016; Hu et al., 2020; Du et al., 2017].

We can extend the abovementioned DGNs to the domain of temporal graphs by selecting appropriate operators in (7) and, in turn, it allows us to tailor the TG-ODE to exploit relational inductive biases and fulfill given application requirements. For instance, by considering a graph attentional operator [Veliˇckovi´c et al., 2018] in (7), we can implement an anisotropic message passing within the temporal graph during the update of node states.

Deep Graph Networks for temporal graphs. Given the sequential structure of temporal graphs, a natural choice for many methods has been to extend Recurrent Neural Networks to graph data. Indeed, most of the models presented in the literature can be summarized as a combination of DGNs and RNNs.

Some approaches adopt a stacked architecture, where DGNs and RNNs are used sequentially [Seo et al., 2018; Pareja et al., 2020], enabling to separately model spatial and temporal dynamics. Other approaches integrate the DGN inside the RNN [Li et al., 2019; Chen et al., 2022; Seo et al., 2018; Li et al., 2018; Zhao et al., 2020; Bai et al., 2021], allowing to jointly capture the temporal evolution and the spatial dependencies in the graph. We refer to [Cini et al., 2023b] and [Gravina and Bacciu, 2024] for a deeper discussion.

Differently from these approaches, which are intrinsically designed to deal with regular time series, TG-ODE can naturally handle arbitrary time gaps between observations. This makes our framework more suitable for realistic scenarios, in which data are irregular over time.

Continuous Dynamic Models. NeuralODE [Chen et al., 2018] (NODE) has emerged as an effective class of neural network models suitable for learning systems' continuous dynamics, drawing a connection between RNNs and ODEs. Despite the similarity with RNNs, such architectures can deal with irregular time series since the continuously-defined dynamics can naturally incorporate data that arrive at arbitrary times [Chen et al., 2018; Rubanova et al., 2019].

Inspired by the NODE approach, GDE [Poli et al., 2019] links DGNs for static graphs with ODEs. In the static graph domain, ODE-based architectures have been proposed with different aims, such as preserving long-range dependencies [Gravina et al., 2023], reducing the computational complexity of message passing [Wang et al., 2021; Wu et al., 2019a], and mitigating the over-smoothing phenomena [Eliasof et al., 2021; Rusch et al., 2022; Kang et al., 2024].

In the temporal domain, TDE-GNN [Eliasof et al., 2024b] employs higher-order ODEs to capture the temporal graph dynamic, while NDCN [Zang and Wang, 2020] extends GDE to learn continuous-time dynamics on both static and temporal graphs by diffusing input features until the termination time. MTGODE [Jin et al., 2022] adopts an ODEbased approach to deduce missing graph topologies from the time-evolving node features in regularly sampled temporal graphs. Differently, [Huang et al., 2020] and [Huang et al., 2021] propose an ODE-based model in the form of a variational auto-encoder for learning latent dynamics from sampled initial states. To infer missing observations, the methods consider both past and future neighbors' information. This prevents them from being used in an online setting, where data becomes available in a sequential order. Lastly, STGNCDE [Choi et al., 2022] employs a stacked architecture of two neural controlled differential equations to model temporal and spatial information, respectively. In the STG-NCDE's paper, irregular data are considered, yet they are handled by making them regular via interpolation.

In contrast to these approaches, in this paper we explicitly address irregularly-sampled temporal graphs and we propose a simple model to showcase the effectiveness and efficiency of the TG-ODE framework in working with such data, eliminating the need for additional strategies, such as interpolation. It should be noted that while many of the ODE-based approaches mentioned earlier can be viewed as instances of the introduced TG-ODE framework, our model is specifically designed to demonstrate the benefits of this approach.

5 Experiments

We provide an empirical assessment of our method against related temporal DGN models from the literature 1 . First, we test the efficacy in handling dynamic graphs with irregularly sampled time series by evaluating the models on several heat diffusion scenarios (see Section 5.1). Afterward, we assess and discuss the performance on real graph benchmarks on traffic forecasting problems (see Section 5.2). We report in Table 1 the grid of hyper-parameters employed in our experiments by each method. We carried out the experiments on 7 nodes of a cluster with 96 CPUs per node.

5.1 Heat Diffusion

In this section, we focus on simulating the heat diffusion over time on a graph. The data is composed of irregularly sampled graph snapshots providing the temperature of the graph's nodes at the given timestamp. We address the task of predicting the nodes' temperature at future (irregular) timestamps.

1 We release the code implementing our methodology and reproducing our empirical analysis at https://github.com/gravins/ TG-ODE.

Datasets

In our experiment, we consider a grid graph consisting of 70 nodes, each of which is characterized by an initial temperature xu(0) randomly sampled in the range between 0 and 0.2. We randomly alter the initial temperature profile by generating hot and cold spikes located at some nodes. A hot spike is characterized by a temperature between 10 and 15, while a cold spike is between -15 and -10. Each altered node has a 40% chance of being associated with a cold spike and 60% with a hot spike. We considered two different experimental scenarios depending on the number of altered nodes. In the first scenario, we alter the temperature of a single node. In the second one, we alter the temperature of one third of the graph's nodes. We will refer to these settings as single-spike and multi-spikes, respectively.

We collected the ground truth by simulating the heat diffusion equation through the forward Euler's method with step size ϵ = 10-3 . Figure 3 illustrates two snapshot graphs from the simulated heat diffusion. The training set consists of 100 randomly selected timestamps over the 1000 steps used to simulate the diffusion process. The validation and test sets are generated from two different simulations similar to the one used for building the training set. However, validation and test sets are obtained through 500-step simulations, and only 50 of them are kept as validation/test sets. We simulated seven different diffusion functions, i.e., -LX(t), -L2X(t), -L5X(t), -tanh(L)X(t), -5LX(t), -0.05LX(t), and -(L+N0,1)X(t). Here, N0,1 stands for a noise sampled from a standard normal distribution, and L is the normalized graph Laplacian.

TG-ODE and Baseline Models

We explored the performance of TG-ODE leveraging the aggregation scheme in [Du et al., 2017] and the forward Euler's method as discretization procedure, for simplicity. Thus, the nodes' states for the entire snapshot are updated as

where K corresponds to the number of neighborhood hops and θk is the k-th weight matrix. We recall that other choices of aggregation and discretization schemes are possible. We compared our method with six common DGNs for dynamic graphs: A3TGCN [Bai et al., 2021], DCRNN [Li et al., 2018], TGCN [Zhao et al., 2020], GCRN-GRU [Seo et al., 2018], GCRN-LSTM [Seo et al., 2018], and NDCN [Zang and Wang, 2020]. We note that whenever we used the NDCN model with embedding dimension set to none (see Table 1), the resulting model corresponds to DNND [Liu et al., 2023].

Moreover, we considered two additional baselines: NODE [Chen et al., 2018] and LB-baseline. NODE represents an instance of our approach that does not take into account node interactions. Instead, LB-baseline returns the same node states received as input (i.e., , the prediction of Xti+1 is ˆX(ti+1) = Xti ) and provides a lower bound on the performance we should expect from the learned models.

We designed each model as a combination of three main components. The first is the encoder which maps the node input features into a latent hidden space; the second is the temporal graph convolution (i.e., TG-ODE or the DGN baselines) or the NODE baseline; and the third is a readout that maps the output of the convolution into the output space. The encoder and the readout are Multi-Layer Perceptrons that share the same architecture among all models in the experiments.

To allow all considered baseline models to handle irregular timestamps, we used a similar strategy employed for TGODE. Specifically, we selected the unit of time, τ, and then we iteratively applied the temporal graph convolution for a number of steps equal to the ratio between the time difference between two consecutive timestamps and the time unit, i.e., #steps = (ti+1 - ti)/τ .

We performed hyper-parameter tuning via grid search, optimizing the Mean Absolute Error (MAE). We trained the models using the Adam optimizer for a maximum of 3000 epochs and early stopping with patience of 100 epochs on the validation error.

Results

We present the results on the heat diffusion tasks in Table 2 and Table 3, using the log10(MAE) as performance metric in both the single-spike and multi-spikes scenarios. The first observation is that TG-ODE has outstanding performance compared to literature models and the baseline. Despite its simpler architecture, our method produces an error that is significantly lower than the runner-up in each task. In the singlespike setting, TG-ODE achieves a log10(MAE) that is on average 308% to 628% better than the competing models in each task.

Interestingly, not all DGN-based models are capable of improving the results of the LB-baseline. This situation suggests that such approaches attempt to merely learn the mapping function between inputs and outputs rather than learning the actual latent dynamics of the system. Such behavior becomes more evident in the more complex multi-spike scenario. Here, our method achieves up to almost 2080% better log10(MAE) score and more literature models fail in improving the performance with respect to the LB-baseline. These results indicate that capturing the latent dynamics is fundamental, in particular, when the time intervals between observations are not regular over time. We conclude that such methods from the literature might not be suitable for more realistic settings characterized by continuous dynamics and sporadic observations.

Finally, we observe that GCRN-GRU and GCRN-LSTM generate the highest error levels, while DCRNN, NDCN, and NODE are the best among the baselines. Since literature models use RNN architectures to learn temporal patterns, it is reasonable to assume that the poor performance might be due to the limited capacity of RNNs to handle non-uniform time gaps between observations. In contrast, ODE-based models (NODE, NCDN and ours) demonstrate enhanced learning capabilities in this scenario. The performance gap between our model and the considered baselines is an indication that the diverse spatial patterns learned by different DGN architectures can heavily impact the performance of the performed tasks.

5.2 Graph Benchmarks

This section introduces a set of graph benchmarks whose objective is to assess traffic forecasting performance from irregular time series; similar to the heat diffusion tasks, we predict the future node values given only the past history.

Datasets

We considered six real-world graph benchmarks for traffic forecasting: MetrLA [Li et al., 2018], Montevideo [Rozemberczki et al., 2021], PeMS03 [Guo et al., 2022], PeMS04 [Guo et al., 2022], PeMS07 [Guo et al., 2022], and PeMS08 [Guo et al., 2022]; we report additional details about the datasets in Table 4. We used a modified version of the original datasets where we employed irregularly sampled observations. We will refer to the datasets by using the subscript

"i" - e.g., MetrLAi - to make apparent the difference from the original versions.

We generated irregular time series by randomly selecting a third of the original graph snapshots for most of the experiments; ratios from 3% to 94% are studied in Figure 5. We considered a temporal data splitting in which 80% of the previously selected snapshots are used as training set, 10% as validation set, and the remaining as test set.

TG-ODE and Baseline Models

For these experiments, we considered the same models, baseline and architectural choices of the heat diffusion experiments. Since NODE does not take into account interactions between nodes for its predictions, we choose not to include it as a baseline in this scenario. Hyper-parameter tuning has been performed by grid search, optimizing the MAE. Optimizer settings are the same as for the previous experiments.

Results

Table 5 reports the traffic forecasting results in terms of MAE. Similarly to the heat diffusion scenario, TG-ODE shows a remarkable performance improvement compared to literature models, achieving an MAE that is up to 202% better than the runner-up model. Moreover, as reported in Figure 4, we observe that TG-ODE is 2× to 13× faster than the other approaches under test. NDCN is the sole method matching the speed of our approach. However, it's noteworthy that NDCN utilizes only one neighbor hop, thereby simplifying the final computation.

We observe that the baseline performs poorly in these benchmarks, suggesting that such tasks are more complex than the ones based on heat diffusion. Despite all DGN models outperforming the LB-baseline, they still produce an error that is on average double than that of TG-ODE, highlighting the added value of our approach when dealing with temporal graphs characterized by irregular sampling.

Finally, we comment that the DCRNN and NDCN suffered from gradient issues in most of the tasks. We believe this is due to their inability to learn the latent dynamics of the system when the models' outputs are not computed over a regular time series.

Impact of the Sample Sparsity

To demonstrate the effectiveness of our approach, we study the prediction performance under different sparsity levels. We consider here the PeMS04 dataset. In this analysis, we systematically decreased the number of considered graph snapshots in the time series. This reduction makes the resulting task more challenging than the original one, as the snapshots become more sparse over time - the expected difference ti+1 - ti gets larger. Additionally, the model has fewer data to learn the task, thereby amplifying the task complexity. We generated the irregular time series by randomly selecting 500, 1000, 2000, 4000, 8000, or 16000 graph snapshots from the original dataset (from 3% to 94% of the original data), resulting in varying degrees of sparsity. For each dataset size, we used an 80/10/10 temporal data split and performed hyperparameter tuning, as previously done in this section.

Figure 5 illustrates the performance of our method, TGODE, at various degrees of data sparsity. As expected, we observe that as the number of samples increases, the test MAE decreases. Notably, TG-ODE maintains robust performance even with higher degrees of sparsity, with a decrease in performance by only ∼7 points when reducing the size from 16000 to 4000 samples. While the prediction error is indeed relatively large when considering only 3-7% of the original data, we comment that it is still substantially better than that of the LB-baseline and comparable to that of the other DGN's which used 33% of the data. Overall, the model exhibits excellent performance even in situations of high to extreme sparsity (i.e., less than 8000 samples). The observed outcome supports the effectiveness of TG-ODE, emphasizing its potential for real-world applications with irregularly sampled temporal graphs.

6 Conclusions

We have presented Temporal Graph Ordinary Differential Equation (TG-ODE), a new general framework for effectively learning from irregularly sampled temporal graphs. Thanks to the connection between ODEs and neural architectures, TGODE can naturally handle arbitrary time gaps between observations, allowing to address a common limitation of DGNs for temporal graphs, i.e., the restriction to work solely on regularly sampled data.

To demonstrate the benefits of our approach, we conducted extensive experiments on ad-hoc benchmarks that include several synthetic and real-world scenarios. The results of our experimental analysis show that our method outperforms state-of-the-art models for temporal graphs by a large margin. Furthermore, our method benefits from a faster training, thus suggesting scalability to large networks.

Despite being appealing for many realistic application setups, we acknowledge that not all resulting ODEs allow unique solutions and yield numerical stability problems, thus requiring some additional care from the user.

Looking ahead to future developments, we intend to broaden the investigation of more sophisticated numerical methods to solve the learned temporal graph ODE, e.g., using adaptive multistep schemes [Ascher and Petzold, 1998]. Extending the proposed framework to the problem of reconstructing missing data is another interesting research direction to consider.

Acknowledgments

This work has been supported by EU-EIC EMERGE (Grant No. 101070918), by the EU NextGenerationEU programme under the funding schemes PNRR-PE-AI (PE00000013) FAIR - Future Artificial Intelligence Research, and by the Swiss National Science Foundation project HORD-GNN (FNS 204061).