Temporal Graph Neural Networks for Irregular Data
This paper proposes a temporal graph neural network model for forecasting of
graph-structured irregularly observed time series. Our TGNN4I model is designed
to handle both irregular time steps and partial observations of the graph. This
is achieved by introducing a time-continuous latent state in each node,
following a linear Ordinary Differential Equation (ODE) defined by the output
of a Gated Recurrent Unit (GRU). The ODE has an explicit solution as a
combination of exponential decay and periodic dynamics. Observations in the
graph neighborhood are taken into account by integrating graph neural network
layers in both the GRU state update and predictive model. The time-continuous
dynamics additionally enable the model to make predictions at arbitrary time
steps. We propose a loss function that leverages this and allows for training
the model for forecasting over different time horizons. Experiments on
simulated data and real-world data from traffic and climate modeling validate
the usefulness of both the graph structure and time-continuous dynamics in
settings with irregular observations.
Joel Oskarsson, Per Sidén, Fredrik Lindsten
# Temporal Graph Neural Networks for Irregular Data

###### Abstract

This paper proposes a temporal graph neural network model for forecasting of graph-structured irregularly observed time series. Our TGNN4I model is designed to handle both irregular time steps and partial observations of the graph. This is achieved by introducing a time-continuous latent state in each node, following a linear Ordinary Differential Equation (ODE) defined by the output of a Gated Recurrent Unit (GRU). The ODE has an explicit solution as a combination of exponential decay and periodic dynamics. Observations in the graph neighborhood are taken into account by integrating graph neural network layers in both the GRU state update and predictive model. The time-continuous dynamics additionally enable the model to make predictions at arbitrary time steps. We propose a loss function that leverages this and allows for training the model for forecasting over different time horizons. Experiments on simulated data and real-world data from traffic and climate modeling validate the usefulness of both the graph structure and time-continuous dynamics in settings with irregular observations.

## 1 Introduction

Many real-world systems can be modeled as graphs. When data about such systems is collected over time, the resulting time series has additional structure induced by the graph topology. Examples of such temporal graph data is the traffic speed in the road network (Li et al., 2018) and the spread of disease in different regions (Rozemberczki et al., 2021). Building accurate machine learning models in this setting requires taking both the temporal and graph structure into account.

While many works have studied the problem of modeling temporal graph data (Wu et al., 2020), these approaches generally assume a constant sampling rate and no missing observations. In real data it is not uncommon to have irregular or missing observations due to non-synchronous measurements or errors in the data collection process. Dealing with such irregularities is especially challenging in the graph setting, as node observations are heavily interdependent. While observations in one node could be modeled using existing approaches for irregular time series (Rubanova et al., 2019; Schirmer et al., 2022), the situation becomes complicated when irregular observations in different nodes occur at different times.

In this paper we tackle two kinds of irregular observations in graph-structured time series: (1) irregularly spaced observation times, and (2) only a subset of nodes being observed at each time point. We propose the TGNN4I model for time series forecasting. The model uses a time-continuous latent state in each node, which allows for predictions to be made at any time point. The latent dynamics are motivated by a linear Ordinary Differential Equation (ODE) formulation, which has a closed form solution. This ODE solution corresponds to an exponential decay (Che et al., 2018) together with an optional periodic component. New observations are incorporated into the state by a Gated Recurrent Unit (GRU) (Cho et al., 2014). Interactions between the nodes are captured by integrating Graph Neural Network (GNN) layers both in the latent state updates and predictive model. To train our model we introduce a loss function that takes into account the time-continuous model formulation and irregularity in the data. We evaluate the model on forecasting problems using traffic and climate data of varying degrees of irregularity and a simulated dataset of periodic signals.

## 2 Related Work

GNNs are deep learning models for graph-structured data (Gilmer et al., 2017; Wu et al., 2020). By learning representations of nodes, edges or entire graphs GNNs can be used for many different machine learning tasks. These models have been successfully applied to diverse areas such as weather forecasting (Lam et al., 2022) molecule
generation (Zang and Wang, 2020) and video classification (Kosman and Di Castro, 2022). Temporal GNNs model also time-varying signals in the graph. This extension to graph-structured time series is achieved by combining GNN layers with recurrent (Li et al., 2018), convolutional (Wu et al., 2019; Yu et al., 2018) or attention (Guo et al., 2019) architectures. While the graph is commonly assumed to be known a priori, some approaches also explore learning the graph structure jointly with the temporal GNN model (Zhang et al., 2022; Wu et al., 2020).

Time series forecasting is a well-studied problem and a vast amount of methods exist in the literature. Traditional methods in the area include ARIMA models, vector auto-regression and Gaussian Processes (Box et al., 2015; Roberts et al., 2013). Many deep learning approaches have also been applied to time series forecasting. This includes Recurrent Neural Networks (RNNs) (Lazzeri, 2020), temporal convolutional neural networks (Chen et al., 2020) and Transformers (Giuliari et al., 2021).

The latent state of an RNN can be extended to continuous time by letting the state decay exponentially in between observations (Che et al., 2018). Such decay mechanisms have been used for modeling data with missing observations (Che et al., 2018), doing imputation (Cao et al., 2018) and parametrizing point processes (Mei and Eisner, 2017). Another way to define time-continuous states is by learning a more general ODE. In neural ODE models (Chen et al., 2018; Kidger, 2021) the latent state is the solution to an ODE defined by a neural network. Neural ODEs have been successfully applied to irregular time series (Rubanova et al., 2019) and can be used to define the dynamics of temporal GNNs (Fang et al., 2021; Poli et al., 2021). Poli et al. (2021) use such a model for graph-structured time series with irregular time steps, but consider the full graph to be observed at each observation time. Also the GraphCON framework of Rusch et al. (2022) combines GNNs with a second order system of ODEs. In GraphCON the time-axis of the ODE is however aligned with the layers of the GNN. This makes the framework suitable for node- and graph-level predictive tasks, rather than time-series modeling.

Another related body of work is concerned with using GNNs for data imputation in time series (Cini et al., 2022; Gordon et al., 2021; Omidshafei et al., 2022). These methods generally do not assume that the time series come with some known graph structure. Instead, the GNN is defined on some graph specifically constructed for the purpose of performing imputation.

The closest work to ours found in the literature is the LG-ODE model of Huang et al. (2020). They consider the same types of irregularities, but are motivated more by a multi-agent systems perspective. LG-ODE is based on an encoder-decoder architecture and trained by maximizing the Evidence Lower Bound (ELBO). The encoder builds a spatio-temporal graph of observations and aggregates information using an attention mechanism (Vaswani et al., 2017). The decoder then extrapolates to future times by solving a neural ODE. The encoder-decoder setup differs from our TGNN4I model that sequentially incorporates observations and auto-regressively makes predictions at every time point. Because of the multi-agent motivation Huang et al. (2020) are also more focused on smaller graphs with few interacting entities, but longer forecasting horizons. An extended version of LG-ODE, called CG-ODE (Huang et al., 2021), also aims to learn the graph structure in the form of a weighted adjacency matrix.

## 3 A Temporal GNN for Irregular Observations

### Setting

Consider a directed or undirected graph \(\mathcal{G}=(V,E)\) with node set \(V\) and edge set \(E\). Let \(\{t_{i}\}_{i=1}^{N_{t}}\) be a set of (possibly irregular) time points s.t. \(0<t_{1}<\cdots<t_{N_{t}}\). We will here present our model for a single time series, but in general we have a dataset containing multiple time series. Let \(\mathcal{O}_{i}\subseteq V\) be the set of nodes observed at time \(t_{i}\). If \(n\in\mathcal{O}_{i}\), we denote the observed value as \(\mathbf{y}_{i}^{n}\in\mathbb{R}^{d_{y}}\) and any accompanying input features as \(\mathbf{x}_{i}^{n}\in\mathbb{R}^{d_{x}}\). We let \(\mathbf{y}_{i}^{n}=\mathbf{x}_{i}^{n}=\mathbf{0}\) if \(n\notin\mathcal{O}_{i}\). Note that this general setting encompasses a spectrum of irregularity, from single node observations (\(|\mathcal{O}_{i}|=1\)\(\forall i\)) to fully observed graphs (\(\mathcal{O}_{i}=V\)\(\forall i\)). A table of notation is given in appendix A.

The problem we consider is that of forecasting. At future time points we want to predict the value at each node, given all earlier observations. Since observations are irregular and we want to make predictions at arbitrary times, we need to consider models that can make predictions for any time in the future.

We consider a model where at each node \(n\) a latent state \(\mathbf{h}^{n}(t)\in\mathbb{R}^{d_{h}}\) evolves over continuous time. We define the dynamics of \(\mathbf{h}^{n}(t)\) by: (1) how \(\mathbf{h}^{n}(t)\) evolves in between observations, and (2) how \(\mathbf{h}^{n}(t)\) is updated when node \(n\) is observed. If node \(n\) is observed at time \(t_{i}\) we incorporate this observation into the latent state using a GRU cell (Cho et al., 2014). This information can then be used for making predictions at future time points. An overview of our model is given in Figure 1.

### Time-continuous Latent Dynamics

Consider a time interval \(|t_{i},t_{j}|\) where node \(n\) is not observed. During this interval we define the latent state of node \(n\) by the sum \(\mathbf{h}^{n}(t)=\mathbf{\bar{h}}_{i}^{n}+\mathbf{\bar{h}}^{n}(t)\). The first part \(\mathbf{\bar{h}}_{i}^{n}\) is constant over the time interval, constituting a base level around which the state evolves. The dynamics of \(\mathbf{\bar{h}}^{n}(t)\) are
dictated by a linear ODE of the form

\[d\tilde{\mathbf{h}}^{n}(t)=A\tilde{\mathbf{h}}^{n}(t)\:dt \tag{1}\]

with \(A\in\mathbb{R}^{d_{h}\times d_{h}}\) and initial condition \(\tilde{\mathbf{h}}^{n}(t_{i})=\tilde{\mathbf{h}}_{i}^{n}\). Over this interval the ODE has a closed form solution (Arrowsmith and Place, 1992) given by

\[\tilde{\mathbf{h}}^{n}(t)=\exp(\delta_{t}A)\tilde{\mathbf{h}}_{i}^{n} \tag{2}\]

where \(\exp\) is the matrix exponential function and we define \(\delta_{t}=t-t_{i}\). Assuming that all eigenvalues of \(A\) are unique, we can use its eigen-decomposition1\(A=Q\Lambda Q^{-1}\) to write

Footnote 1: Recall that the eigen-decomposition of a diagonalizable matrix \(A\) is given by \(A=Q\Lambda Q^{-1}\), where \(\Lambda\) is a diagonal matrix containing the eigenvalues of \(A\) and the columns of \(Q\) are the corresponding eigenvectors (Searle and Khuri, 1982).

\[\tilde{\mathbf{h}}^{n}(t)=Q\exp(\delta_{t}\Lambda)Q^{-1}\tilde{\mathbf{h}}_{i}^{n}. \tag{3}\]

Since \(Q\) contains an eigen-basis of \(\mathbb{R}^{d_{h}}\) it can be viewed as a transition matrix, changing the basis of the latent space. While we could in principle learn \(Q\), we note that the basis of the latent space has no physical interpretation and we can without loss of generality choose it such that \(Q=I\).

Next, if we make the assumption that \(A\) has real eigenvalues the resulting dynamics are given by

\[\tilde{\mathbf{h}}^{n}(t)=\exp(-\delta_{t}\operatorname{diag}(\mathbf{\omega}_{i}^{n} ))\tilde{\mathbf{h}}_{i}^{n} \tag{4}\]

where \(\mathbf{\omega}_{i}^{n}>0\) is a parameter vector representing the negation of the eigenvalues. We arrive at an exponential decay in-between observations, a type of dynamics used with GRU-updates in existing works (Che et al., 2018; Cao et al., 2018). The positive restriction on \(\mathbf{\omega}_{i}^{n}\) ensures the stability of the dynamical system in the limit, which for the complete state means that \(\mathbf{h}^{n}(t_{j})\to\tilde{\mathbf{h}}_{i}^{n}\) as \(t_{j}\to\infty\).

On the other hand, if we instead allow \(A\) to have complex eigenvalues we can decompose \(A\) using the real Jordan form (Horn and Johnson, 2012). This makes \(\Lambda\) block-diagonal with \(2\times 2\) blocks

\[C_{j}=\left[\begin{array}{cc}a_{j}&-b_{j}\\ b_{j}&a_{j}\end{array}\right] \tag{5}\]

corresponding to complex conjugate eigenvalues \(a_{j}\pm b_{j}i\). If we compute the matrix exponential for this \(\Lambda\) we end up with a combination of exponential decay and periodic dynamics (Arrowsmith and Place, 1992). The solution gives dynamics that couple each pair of dimensions as

\[\begin{split}\left[\tilde{\mathbf{h}}^{n}(t)\right]_{j:j+1}& =\exp(\delta_{t}a_{j})\times\\ &\left[\begin{array}{cc}\cos(b_{j}\delta_{t})&-\sin(b_{j}\delta_{t})\\ \sin(b_{j}\delta_{t})&\cos(b_{j}\delta_{t})\end{array}\right]\left[\tilde{ \mathbf{h}}_{i}^{n}\right]_{j:j+1}.\end{split} \tag{6}\]

We parametrize also these dynamics with \(\mathbf{\omega}_{i}^{n}=[-a_{1},-a_{3},\ldots,-a_{d_{h}-1},b_{1},b_{3},\ldots,b_{ d_{h}-1}]^{\intercal}>0\). Note that when \(b_{j}\to 0\) these periodic dynamics reduce to the exponential decay in Eq. 4, but with the parameter for \(-a_{j}\) shared across pairs of dimensions. We consider both the exponential decay and the more general periodic dynamics as options for our model. The periodic dynamics can naturally be seen as advantageous for modeling periodic data, as we will explore empirically in section 4.4.

### Incorporating Observations from the Graph

When node \(n\) is observed at time \(t_{i}\) the observation is incorporated into the latent state by a GRU cell, incurring an

Figure 1: **Left:** Example graph with four nodes and their latent states, here one-dimensional for illustration purposes. Node observations are indicated with \(\bigcirc\). **Right:** Schematic diagram of the GRU update and predictive model \(g\). Everything inside the shaded area happens instantaneously at time \(t_{i}\). The GRU cell outputs the new initial state \(\tilde{\mathbf{h}}_{i}^{n}+\tilde{\mathbf{h}}_{i}^{n}\), the static component \(\tilde{\mathbf{h}}_{i}^{n}\) and the ODE parameters \(\mathbf{\omega}_{i}^{n}\). The parameters \(\mathbf{\omega}_{i}^{n}\) define the dynamics until the next observation. The prediction \(\tilde{\mathbf{g}}_{i}^{n}\) here is based on information from all earlier time points.
instantaneous jump in the latent dynamics to a new value \(\mathbf{\bar{R}}_{i}^{n}+\mathbf{\bar{h}}_{i}^{n}\). Inspired by the continuous-time LSTM of Mei and Eisner (2017), we extend the GRU cell to output also the parameters \(\mathbf{\omega}_{i}^{n}\), which define the dynamics of \(\mathbf{\bar{h}}^{n}(n)\) for the next time interval.

So far we have considered each node of the graph as separate entities, but in a graph-based system future observations of a node can depend on the history of the entire graph. To capture this we let the state of each node depend on observations and states in its graph neighborhood. This is achieved by introducing GNNs Gilmer et al. (2017); Wu et al. (2020) in the GRU update Zhao et al. (2018). We replace matrix multiplications with GNNs, taking inputs both from the node \(n\) itself and from its neighborhood \(\mathcal{N}(n)=\{m|(m,n)\in E\}\). The type of GNN we use is a simple version of a message passing neural network Gilmer et al. (2017), defined as

\[\begin{split}\mathrm{GNN}&\Big{(}\mathbf{h}^{n},\{\mathbf{h }^{m}\}_{m\in\mathcal{N}(n)}\Big{)}\\ &=W_{1}\mathbf{h}^{n}+\frac{1}{|\mathcal{N}(n)|}\sum_{m\in\mathcal{N} (n)}\!\!\!e_{m,n}W_{2}\mathbf{h}^{m}\end{split} \tag{7}\]

where the matrices \(W_{1},W_{2}\) are learnable parameters shared among all nodes and \(e_{m,n}\) is an edge weight associated with the edge \((m,n)\). The use of edge weights allows for incorporating prior information about the strength of connections in the graph. We can additionally stack multiple such GNN layers, append fully-connected layers and include non-linear activation functions in between. The inclusion of multiple GNN layers makes the GRU update dependent on a larger graph neighborhood.

In a standard GRU cell the input and previous state are first mapped to three new representations using matrices \(U\) and \(W\)Cho et al. (2014). These representations are then used to compute the state update. In our GNN-based GRU update the matrices are replaced by GNNs and we require seven such intermediate representations to update \(\mathbf{\hat{h}}_{i}^{n}\), \(\mathbf{\bar{h}}_{i}^{n}\) and \(\mathbf{\omega}_{i}^{n}\), as shown below. The node states are combined as

\[[\mathbf{u}_{i,1}^{n},\dots,\mathbf{u}_{i,7}^{n}]^{\intercal}=\mathrm{GNN}^{U}\Big{(} \mathbf{h}^{n}(t_{i}),\{\mathbf{h}^{m}(t_{i})\}_{m\in\mathcal{N}(n)}\Big{)} \tag{8}\]

where the resulting vector is split into seven equally sized chunks. Another GNN is then used for the combined observations and input features \(\mathbf{\tilde{x}}_{i}^{n}=[\mathbf{y}_{i}^{n},\mathbf{x}_{i}^{n}]^{\intercal}\),

\[[\mathbf{v}_{i,1}^{n},\dots,\mathbf{v}_{i,7}^{n}]^{\intercal}=\mathrm{GNN}^{W}\Big{(} \mathbf{\tilde{x}}_{i}^{n},\{\mathbf{\tilde{x}}_{i}^{m}\}_{m\in\mathcal{N}(n)}\Big{)}. \tag{9}\]

Note that while all nodes might not be observed at time \(t_{i}\), \(\mathbf{\tilde{x}}_{i}^{m}=\mathbf{0}\) for any unobserved \(m\) and thus do not contribute to the sum over neighbors in the GNN. With the combined information from the graph neighborhood the full GRU update is computed as

\[\mathbf{r}_{i}^{n} =\sigma(\mathbf{v}_{i,1}^{n}+\mathbf{u}_{i,1}^{n}+\mathbf{b}_{1}) \tag{10a}\] \[\mathbf{z}_{i}^{n} =\sigma(\mathbf{v}_{i,2}^{n}+\mathbf{u}_{i,2}^{n}+\mathbf{b}_{2})\] (10b) \[\mathbf{q}_{i}^{n} =\tanh(\mathbf{v}_{i,3}^{n}+(\mathbf{r}_{i}^{n}\odot\mathbf{u}_{i,3}^{n})+\bm {b}_{3})\] (10c) \[\mathbf{\bar{h}}_{i}^{n}+\mathbf{\hat{h}}_{i}^{n} =(\mathbf{1}-\mathbf{z}_{i}^{n})\odot\mathbf{h}^{n}(t_{i})+\mathbf{z}_{i}^{n} \odot\mathbf{q}_{i}^{n}\] (10d) \[\mathbf{\bar{r}}_{i}^{n} =\sigma(\mathbf{v}_{i,4}^{n}+\mathbf{u}_{i,4}^{n}+\mathbf{b}_{4})\] (11a) \[\mathbf{\bar{z}}_{i}^{n} =\sigma(\mathbf{v}_{i,5}^{n}+\mathbf{u}_{i,5}^{n}+\mathbf{b}_{5})\] (11b) \[\mathbf{\bar{q}}_{i}^{n} =\tanh(\mathbf{v}_{i,6}^{n}+(\mathbf{\bar{r}}_{i}^{n}\odot\mathbf{u}_{i,6}^{n} )+\mathbf{b}_{6})\] (11c) \[\mathbf{\bar{h}}_{i}^{n} =(\mathbf{1}-\mathbf{\bar{z}}_{i}^{n})\odot\mathbf{\bar{h}}_{i-1}^{n}+\mathbf{\bar {z}}_{i}^{n}\odot\mathbf{\bar{q}}_{i}^{n}\] (11d) \[\mathbf{\omega}_{i}^{n} =\log(\mathbf{1}+\exp(\mathbf{v}_{i,7}^{n}+\mathbf{u}_{i,7}^{n}+\mathbf{b}_{7})) \tag{12}\]

where \(\sigma\) is the sigmoid function and \(\mathbf{b}_{1}\)-\(\mathbf{b}_{7}\) learnable bias parameters. In Eq. 11d we let \(\mathbf{\bar{h}}_{i}^{n}=\mathbf{\bar{h}}_{k-1}^{n}\) if \(n\notin\mathcal{O}_{k}\). Eq. 10a-10d correspond to one GRU update, using the decayed state \(\mathbf{h}^{n}(t_{i})\). Eq. 11a-11d define a separate GRU update, but for the decay target \(\mathbf{\bar{h}}_{i}^{n}\). Note that we get \(\mathbf{\hat{h}}_{i}^{n}\), the initial value of \(\mathbf{\tilde{h}}^{n}(t)\), implicitly from the difference between Eq. 10d and Eq. 11d. Finally Eq. 12 computes the parameters \(\mathbf{\omega}_{i}^{n}\) defining the dynamics of \(\mathbf{\bar{h}}^{n}(t)\) up until the next observation of node \(n\). The parameters of the GRU cell are shared for all nodes in the graph. To capture any node-specific properties we parametrize initial states \(\mathbf{h}^{n}(0)\) separately for all nodes and learn these jointly with the rest of the model.

### Predictions

The time-continuous dynamics ensure that there is a well-defined latent state \(\mathbf{h}^{n}(t)\) in each node at each time point \(t\). The value of the time series can then be predicted at any time by applying a mapping \(g\colon\mathbb{R}^{d_{h}}\to\mathbb{R}^{d_{y}}\) from this latent state to the prediction \(\mathbf{\tilde{y}}_{j}^{n}\).

The addition of GNNs into the GRU update makes the latent state dynamics of each node dependent on historical observations in its neighborhood. However, since the GRU updates happen only when a node is observed, information from observed neighbors might not be incorporated immediately in the latent state. Consider three consecutive time points \(t_{i}<t_{i+1}<t_{i+2}\) s.t. \(n\in\mathcal{O}_{i}\), \(n\notin\mathcal{O}_{i+1}\). Then any observation \(\mathbf{y}_{i+1}^{m}\) for \(m\in\mathcal{O}_{i+1}\cap\mathcal{N}(n)\) will not be taken into account by the model for the prediction \(\mathbf{\hat{y}}_{i+2}^{n}\), as that prediction is based on a latent state with only information from time \(t_{i}\). To remedy this we choose also the predictive model \(g\) to contain one or more GNN layers,

\[\mathbf{\tilde{y}}_{j}^{n}=\mathrm{GNN}^{g}\Big{(}[\mathbf{h}^{n}(t_{j}),\mathbf{x}_{j}^{n}]^ {\intercal},\big{\{}[\mathbf{h}^{n}(t_{j}),\mathbf{x}_{j}^{m}]^{\intercal}\big{\}}_{m\in \mathcal{N}(n)}\Big{)}. \tag{13}\]

This way \(g\) takes the latent states and input features of the whole neighborhood into account for prediction. We name the full proposed model Temporal Graph Neural Network for Irregular data (**TGNN4I**).
### Loss Function

In order to make predictions for arbitrary future time points we introduce a suitable loss function based on the time-continuous nature of the model. Let \(\hat{\mathbf{y}}_{i\to j}^{m}\) be the prediction for node \(m\) at time \(t_{j}\), based on observations of all nodes at times \(t\leq t_{i}\). Define also a time-continuous weighting function \(w\colon\mathbb{R}^{+}\to\mathbb{R}^{+}\) and the set \(\tau_{m,i}=\{j:m\in\mathcal{O}_{j}\wedge i<j\}\) containing the indices of all times after \(t_{i}\) where node \(m\) is observed. We do not include predictions from the first \(N_{\text{init}}\) time steps in the loss, treating this as a short warm-up phase. The loss function for one graph-structured time series is then

\[\mathcal{L}_{\ell}=\frac{1}{N_{\text{obs}}}\sum_{m\in V}\sum_{i=N_{\text{init }}+1}^{N_{t}}\sum_{j\in\tau_{m,i}}\frac{\ell\big{(}\hat{\mathbf{y}}_{i\to j}^{m}, \mathbf{y}_{j}^{m}\big{)}w(t_{j}-t_{i})}{j-N_{\text{init}}-1} \tag{14}\]

where \(\ell\) is any loss function for a single observation and \(N_{\text{obs}}=\sum_{i=N_{\text{init}}+2}^{N_{t}}\lvert\mathcal{O}_{i}\rvert\) the total number of node observations. We use \(\mathcal{L}_{\text{MSE}}\) with Mean Squared Error (MSE) as \(\ell\), but the framework is fully compatible with other loss functions as well. This includes general probabilistic predictions with a negative log-likelihood loss. Dividing by \(j-N_{\text{init}}-1\), the number of times observation \(j\) has been predicted, guarantees that later observations are not given a higher total weight.

The weighting function \(w\) allows for specifying which time-horizons that should be prioritized by the model. This choice is highly application-dependent and should capture which predictions that are of interest when later deploying the model in some real-world setting. If we care about all time horizons, but want to prioritize predictions close in time, a suitable choice could be \(w(\Delta_{t})=\exp\bigl{(}-\frac{\Delta_{t}}{\Omega}\bigr{)}\). It might also be desirable to focus predictive capabilities on a specific \(\Delta_{t}\). If we want predictions around \(\Delta_{t}=\mu\) to be prioritized we can for instance use a Gaussian kernel as

\[w(\Delta_{t})=\exp\biggl{(}-\Bigl{(}\tfrac{\Delta_{t}-\mu}{\Omega}\Bigr{)}^{2 }\biggr{)}. \tag{15}\]

A limitation of the proposed loss function is the quadratic scaling in the number of time steps, as predictions are made from all times to all future observations. This especially requires large amounts of memory for nodes that are observed at many time steps. However, for many sensible choices of \(w\) predictions far into the future have a close to \(0\) impact on the loss. In practice we can utilize this to approximate \(\mathcal{L}_{\ell}\) by only making predictions \(N_{\text{max}}\) time steps into the future. This approximation explicitly corresponds to setting \(\tau_{m,i}=\{j:m\in\mathcal{O}_{j}\wedge i<j\leq i+N_{\text{max}}\}\) and changing the denominator in Eq. 14 to \(\min(N_{\text{max}},j-N_{\text{init}}-1)\). Alternatively, we can select a weight function with finite support which implies that many terms in Eq. 14 will be exactly zero. This does however require more bookkeeping than the aforementioned truncation method.

## 4 Experiments

The TGNN4I model was implemented2 using PyTorch and PyG (Fey and Lenssen, 2019). We evaluate the model on a number of different datasets. See appendix D and E for details on the pre-processing and experimental setups used. As the loss function \(\mathcal{L}_{\text{MSE}}\) captures errors throughout an entire time series we adopt this also as our evaluation metric. Given that the loss weighting \(w\) used for training accurately represents how we value predictions at different time horizons, it is natural to use the same choice for evaluation. In our experiments we rescale each time series so that \(t\in[0,1]\) and use \(w(\Delta_{t})=\exp\bigl{(}-\tfrac{\Delta_{t}}{0.04}\bigr{)}\). By inspecting \(w\) and the time steps in the data we also choose a suitable \(N_{\text{max}}=10\).

Footnote 2: Our code and datasets are available at [https://github.com/joeloskarsson/tgnn4i](https://github.com/joeloskarsson/tgnn4i).

We consider three versions of our TGNN4I model: (**static**) with a constant latent state in-between observations (\(\tilde{\mathbf{h}}^{n}(t)=\tilde{\mathbf{h}}_{i}^{n}\ \forall t\!\in\!\lvert t_{i},t_{j}\rvert\)), (**exponential**) with the exponential decay dynamics from Eq. 4 and (**periodic**) with the combined decay and periodic dynamics from Eq. 6. In all our experiments the training time of a single model on an NVIDIA A100 GPU is less than an hour and for the smallest dataset (METR-LA) not more than 20 minutes.

### Baselines

We compare TGNN4I to multiple baseline models. As a simple starting point we consider a model that always predicts the last observed value in each node for all future time points (**Predict Previous**). Che et al. (2018) propose the GRU-D model for irregular time series, which we extend with our parametrization of the exponential decay and include as a baseline. GRU-D does not use the graph structure explicitly, so there are two ways to adapt this model to our setting. We can view the entire graph-structured time series as one series with \((\lvert V\rvert d_{y})\)-dimensional vectors at each time step (**GRU-D (joint)**). Alternatively, we can view the time series in each node as independent (**GRU-D (node)**), which is essentially the same as TGNN4I with all edges in the graph removed. Two **Transformer** baselines are also included, used in the same (joint) and (node) configurations. In these models the irregular observations are handled through attention masks and the use of timestamps in the sinusoidal positional encodings. We also compare against the **LG-ODE** model of Huang et al. (2020) using the code provided by the authors. We follow their proposed training procedure, where the model encodes the first half of each time series and has to predict the second half. When computing \(\mathcal{L}_{\text{MSE}}\) using LG-ODE we encode all observations up to \(t_{i}\) and decode from that time point in order to get each \(\hat{\mathbf{y}}_{i\to j}^{m}\). More details on the baseline models are given in appendix C. An attempt was also made to adapt the RAINDROP model of Zhang et al. (2022) to our fore
casting setting. We were however unable to get useful predictions without making major changes to the model and it is therefore not included here.

### Traffic Data

We experiment on the PEMS-BAY and METR-LA datasets, containing traffic speed sensor data from the California highway system (Li et al., 2018; Chen et al., 2001). To be able to control the degree of irregularity, we start from regularly sampled data and choose subsets of observations. We use the versions of the datasets pre-processed by Li et al. (2018). Each dataset is split up into time series of 288 observations (1 day). PEMS-BAY contains 180 such time series with 325 nodes and METR-LA 119 time series with 207 nodes. We include the time of day and the time since the node was last observed as input features \(\mathbf{x}_{i}^{n}\). In order to introduce irregularity in the time steps we next subsample each time series by keeping only 72 of the 288 observations. These \(N_{t}=72\) observation times are the same for all nodes. However, from these subsampled time series we furthermore sample subsets containing 25%-100% of all \(N_{t}\times|V|\) individual node observations. This results in irregular observation time points and a fraction of nodes observed at each time. Our additional pre-processing prevents us from a direct comparison with Li et al. (2018), as their method does not handle irregular observations.

We report results for both datasets in Table 1 and highlight the best performing models on METR-LA in Figure 2. GRU-D (joint) has a hard time modeling all nodes jointly, often not performing better than the simple Predict Previous baseline. The Transformer models achieve somewhat better results, but still not competitive with TGNN4I. We additionally note that the Transformers can be highly sensitive to the random seed used for initialization, something that we have not observed for other models. Comparing TGNN4I and GRU-D (node) in Figure 2 it can be noted that the importance of using the graph structure increases when there are fewer observations. Out of the different versions of TGNN4I the exponential and periodic dynamics show a clear advantage over the static one, with the largest difference for the most sparsely observed data. We have observed that the periodic models output only low frequencies, resulting in dynamics and results similar to the model with exponential decay. While the periodic dynamics in Eq. 6 have fewer degrees of freedom when reduced to pure exponential decay, this does not seem to hurt the performance in this example.

We found that the training time of LG-ODE scales poorly to large graphs, limiting us to only training a single model for each dataset. The predictions are however quite poor, especially on the PEMS-BAY data. While the model seems to learn something more than just predicting the mean, it is not competitive with our TGNN4I model. We believe that the poor performance of LG-ODE can be explained by a combination of multiple things: (1) The LG-ODE model is primarily designed for data with clear continuous underlying dynamics, which might not match this type of traffic data. (2) When the model is trained as proposed by Huang et al. (2020), it can require large amounts of data. For some of the experiments in the original paper 20 000 sequences are used for training, while we use less than 150. (3) The slow training has limited possibilities for exhaustive hyperparameter tuning on our datasets. Training one LG-ODE model on the PEMS-BAY data takes us over 50 hours.

### USHCN Climate Data

Irregular and missing observations are common problems in climate data (Schneider, 2001). The United State Historical Climatology Network (USHCN) daily dataset contains over 50 years of measurements of multiple climate variables from sensor stations in the United States (Menne et al., 2015). We use the pre-processing of De Brouwer et al. (2019) to clean and subsample the data. The target variables chosen are minimum and maximum daily temperature (\(T_{\rm min}\) and \(T_{\rm max}\)), which we model as separate datasets. While existing works (De Brouwer et al., 2019; Schirmer et al., 2022) have treated time series from different sensor stations as independent, we model also the spatial correlation by constructing a 10-nearest-neighbor graph using the sensor positions. Each full dataset contains 186 time series of length \(N_{t}=100\) on a graph of 1123 nodes. The pre-processed USHCN data is sparsely observed with only around 5% of potential node observations present.

We report results on both datasets in Table 2. Due to the large size of the graph it was not feasible to apply the LG-ODE model here. We note that for these datasets the (joint) baselines clearly outperform the (node) versions. For GRU-D this is the opposite of what we saw in the traffic data. This can be explained by the fact that climate data

Figure 2: Test \(\mathcal{L}_{\rm MSE}\) on METR-LA traffic data for the best performing models from Table 1. Shaded areas correspond to 95% confidence intervals based on re-training models with 5 random seeds.
has strong spatial dependencies. The (joint) models can to some extent learn to pick up on these, while for (node) no information can flow between nodes. The best results are however achieved by TGNN4I, showing the added benefit of utilizing the spatial graph.

### Synthetic Periodic Data

In the previous experiments, using periodic dynamics with TGNN4I has not added any value. Instead, the learned dynamics have been largely similar to just using exponential decay. This should to some extent be expected, as none of the previous datasets show any clear periodic patterns at the considered time scales. To investigate the possible benefits of the periodic dynamics we instead create a synthetic dataset with periodic signals propagating over a graph.

The synthetic dataset is based on a randomly sampled directed acyclic graph with 20 nodes. We define a periodic base signal

\[\rho^{n}(t)=\sin(\phi^{n}t+\eta^{n}) \tag{16}\]

with random parameters \(\phi^{n}\) and \(\eta^{n}\) for each node \(n\). The target signal \(y^{n}(t)\) in each node is then defined through

\[\kappa^{n}(t) =\rho^{n}(t)+\frac{0.5}{|\mathcal{N}(n)|}\sum_{m\in\mathcal{N}(n)} \kappa^{m}(t-0.05) \tag{17a}\] \[y^{n}(t) =\kappa^{n}(t)+\epsilon^{n}(t) \tag{17b}\]

where \(\epsilon^{n}(t)\) is Gaussian white noise with standard deviation 0.01. The target signal in each node depends on the base signal in the node itself and the signals in neighboring nodes at a time lag of 0.05. To construct one time series we sample \(y^{n}(t)\) at 70 irregular time points on \([0,1]\). In total we sample 200 such time series and keep 50% of the node observations in each.

We train versions of the GRU-D (node) and TGNN4I models with different latent dynamics on the synthetic data.

\begin{table}
\begin{tabular}{l c c c c} \hline \hline  & \multicolumn{4}{c}{**PEMS-BAY**} \\ \cline{2-5}
**Model** & 25\% & 50\% & 75\% & 100\% \\ \hline Predict Previous & 26.32 & 18.60 & 15.25 & 13.50 \\ GRU-D (joint) & 18.79\(\pm\)0.07 & 18.27\(\pm\)0.10 & 17.93\(\pm\)0.08 & 17.75\(\pm\)0.12 \\ GRU-D (node) & 8.79\(\pm\)0.06 & 6.62\(\pm\)0.02 & 5.82\(\pm\)0.06 & 5.49\(\pm\)0.06 \\ Transformer (joint) & 12.05\(\pm\)1.19 & 13.13\(\pm\)2.59 & 12.21\(\pm\)1.95 & 11.09\(\pm\)1.38 \\ Transformer (node) & 16.49\(\pm\)0.17 & 14.44\(\pm\)0.48 & 13.20\(\pm\)0.56 & 13.16\(\pm\)1.23 \\ LG-ODE & 27.00 & 24.93 & 24.71 & 23.52 \\ TGNN4I (static) & 7.41\(\pm\)0.09 & 5.98\(\pm\)0.07 & 5.29\(\pm\)0.08 & 4.89\(\pm\)0.05 \\ TGNN4I (exponential) & **7.10\(\pm\)0.07** & **5.78\(\pm\)0.05** & 5.23\(\pm\)0.03 & **4.87\(\pm\)0.09** \\ TGNN4I (periodic) & **7.10\(\pm\)0.09** & 5.80\(\pm\)0.08 & **5.22\(\pm\)0.09** & **4.87\(\pm\)0.02** \\ \hline \hline \multicolumn{5}{c}{**METR-LA**} \\ \cline{2-5}
**Model** & 25\% & 50\% & 75\% & 100\% \\ \hline Predict Previous & 9.86 & 7.54 & 6.52 & 6.04 \\ GRU-D (joint) & 8.38\(\pm\)0.05 & 8.03\(\pm\)0.04 & 7.89\(\pm\)0.03 & 7.80\(\pm\)0.02 \\ GRU-D (node) & 4.36\(\pm\)0.08 & 3.62\(\pm\)0.07 & 3.28\(\pm\)0.08 & 3.16\(\pm\)0.04 \\ Transformer (joint) & 5.70\(\pm\)1.41 & 7.17\(\pm\)1.66 & 5.95\(\pm\)1.90 & 6.11\(\pm\)1.80 \\ Transformer (node) & 7.01\(\pm\)0.31 & 6.34\(\pm\)0.24 & 5.84\(\pm\)0.23 & 5.96\(\pm\)0.50 \\ LG-ODE & 8.51 & 7.35 & 6.71 & 6.24 \\ TGNN4I (static) & 3.86\(\pm\)0.02 & 3.31\(\pm\)0.02 & 3.03\(\pm\)0.02 & 2.88\(\pm\)0.02 \\ TGNN4I (exponential) & **3.68\(\pm\)0.05** & **3.18\(\pm\)0.03** & **2.97\(\pm\)0.03** & **2.86\(\pm\)0.04** \\ TGNN4I (periodic) & 3.69\(\pm\)0.02 & 3.19\(\pm\)0.04 & 3.01\(\pm\)0.05 & 2.88\(\pm\)0.03 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Test \(\mathcal{L}_{\text{MSE}}\) (multiplied by \(10^{2}\)) for the traffic datasets with different fractions of node observations. Where applicable we report mean \(\pm\) one standard deviation across 5 runs with different random seeds. The lowest mean \(\mathcal{L}_{\text{MSE}}\) for each dataset and observation percentage is marked in bold.

\begin{table}
\begin{tabular}{l c c} \hline \hline  & \(T_{\min}\) & \(T_{\max}\) \\ \hline Predict Previous & 16.88 & 17.18 \\ GRU-D (joint) & 8.03\(\pm\)0.23 & 7.97\(\pm\)0.19 \\ GRU-D (node) & 13.12\(\pm\)0.03 & 13.67\(\pm\)0.04 \\ Transformer (joint) & 7.36\(\pm\)0.41 & 7.37\(\pm\)0.28 \\ Transformer (node) & 15.68\(\pm\)0.32 & 15.74\(\pm\)0.34 \\ TGNN4I (static) & 6.97\(\pm\)0.05 & 6.86\(\pm\)0.04 \\ TGNN4I (exponential) & **6.72\(\pm\)0.04** & **6.60\(\pm\)0.04** \\ TGNN4I (periodic) & **6.72\(\pm\)0.05** & 6.63\(\pm\)0.03 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Test \(\mathcal{L}_{\text{MSE}}\) (multiplied by \(10^{2}\)) for the two USHCN climate datasets.
Also our other baselines are included for comparison. Results are reported in Table 3. For both GRU-D (node) and TGNN4I we see a large difference between the different types of latent dynamics. The periodic dynamics seem to help the model to keep track of the base signal in the node and its neighborhood in order to achieve accurate future predictions. While this is a synthetic example, periodic behavior is prevalent in much time series data and being able to explicitly model this in the latent state can be highly advantageous. Attempts were made to also train the GRU-D (joint) model on this dataset, but it failed to pick up on any patterns and ended up only predicting a constant value for all nodes and times.

### Loss Weighting

To investigate the impact of the loss weighting function \(w\) we trained four TGNN4I models on the subsampled PEMS-BAY dataset with 25% observations. We used exponential dynamics and considered the weighting functions

\[w_{1}(\Delta_{t}) =1 \tag{18a}\] \[w_{2}(\Delta_{t}) =\exp\!\left(-\frac{\Delta_{t}}{0.04}\right)\] (18b) \[w_{3}(\Delta_{t}) =\exp\!\left(-\!\left(\frac{\Delta_{t}-0.1}{0.02}\right)^{2}\right)\] (18c) \[w_{4}(\Delta_{t}) =\mathbb{I}_{\{\Delta_{t}\in[0.18,0.22]\}}. \tag{18d}\]

Figure 3 shows the test MSE for the trained models at different \(\Delta_{t}\) in the future, as well as plots of the weighting functions. As an example, the prediction \(\hat{\mathbf{y}}_{i\to j}^{m}\) has \(\Delta_{t}=t_{j}-t_{i}\).

We note that the choice of \(w\) can have substantial impact on the error of the model at different time horizons. The exponential weighting in \(w_{2}\) makes the model focus heavily on short-term predictions. This results in better predictions for low \(\Delta_{t}\). At the shortest time horizon the exponential weighting yields an 11% improvement over the model trained with constant \(w_{1}\), but this comes at the cost of far higher errors for long-term predictions. Interestingly the \(w_{1}\) model gives better predictions at all time horizons than the models with \(w_{3}\) and \(w_{4}\), which focus on predictions at some specific time ahead. We believe that there can be a feedback effect benefiting the constant weighting, where learning to make good short-term predictions also aid the learning of long-term prediction, for example by finding useful intermediate representations. A drawback of weighting with \(w_{1}\) is however that since the loss never approaches 0 there is no properly motivated choice of \(N_{\text{max}}\) for our loss approximation. As the model trained with \(w_{4}\) still gives good predictions in the interval \([0.18,0.22]\), there can still be practical reasons to choose such a weighting. With this choice we could reduce the innermost sum in Eq. 14 to only those \(j\):s that lie in the interval of interest.

## 5 Discussion

We have proposed a temporal GNN model that can handle both irregular time steps and partially observed graphs. By defining latent states in continuous time our model can make predictions for arbitrary time points in the future. In this section we discuss some details and limitations of the approach, and also give some pointers to interesting directions for future work.

### Efficient Implementation

In order to efficiently implement the training and inference of TGNN4I there is a key design choice between (1) storing everything in dense matrices and utilizing massively parallel GPU-computations, and (2) utilizing sparse representations in order to avoid computing values that will never be

Figure 3: MSE (Top) for predictions at time \(\Delta_{t}\) in the future for models trained with different loss weighting \(w\) (Bottom). The weighting functions are described in Eq. 18a–18d. To compute the MSE all predictions in the test set were binned based on their \(\Delta_{t}\), with bin width 0.02. The bottom subplot also shows a histogram of the number of predictions in each bin.

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & **Static** & **Exponential** & **Periodic** \\ \hline GRU-D (node) & 8.88\(\pm\)0.36 & 3.13\(\pm\)0.06 & 2.81\(\pm\)0.04 \\ TGNN4I & 15.12\(\pm\)0.05 & 2.91\(\pm\)0.17 & **1.95\(\pm\)0.11** \\ \hline Predict Prev. & 27.52 & \\ Transformer (joint) & 23.19\(\pm\)0.38 & \\ Transformer (node) & 15.39\(\pm\)0.05 & \\ LG-ODE & 16.61\(\pm\)0.23 & \\ \hline \hline \end{tabular}
\end{table}
Table 3: Test \(\mathcal{L}_{\text{MSE}}\) (multiplied by \(10^{2}\)) for synthetic data.
used. Which of these is to be preferred depends on the sparsity of observations in the data. Our implementation follows the massively parallel approach, using binary masks for keeping track of \(\mathcal{O}_{i}\) at each time point. In order to scale TGNN4I to massive graphs in real world scenarios it could be interesting to consider a version of the model distributed over multiple machines, perhaps directly connected to sensors producing the input data. A sparse implementation would be strongly preferred for such an extension.

### Linear and Neural ODEs

The dynamics of TGNN4I are defined by a linear ODE and additionally restricted by the assumption of unique eigenvalues in \(A\). This has the benefit of a closed form solution that is efficient to compute, but also limits the types of dynamics that can be learned. Our approach can be contrasted with Neural ODEs (Chen et al., 2018), that allow for learning more expressive dynamics. Neural ODEs do however lack closed form solutions and require using numerical solvers (Kidger, 2021). This incurs a trade-off between speed and numerical accuracy. In experiments we have compared TGNN4I with the LG-ODE model (Huang et al., 2020), which uses a Neural ODE decoder. The slow training of the LG-ODE model can to a large extent be attributed to the numerical ODE solver. While more complex latent dynamics can be useful for some datasets, it can also be argued that simpler dynamics can be compensated with a high enough latent dimension \(d_{h}\) and a flexible enough predictive model \(g\)(Schirmer et al., 2022).

### Societal and Sustainability Impact

While our contributions are purely methodological, many applications of graph-based and spatio-temporal data analysis have a clear societal impact. Our example applications of traffic and climate modeling both have potential to aid efforts of transforming society in more sustainable directions, such as those described in the United Nations sustainable development goals 11 and 13 (Rolnick et al., 2022; United Nations, 2015). Traffic modeling allows us to both understand travel behavior and predict future demand. This can enable optimizations of transport systems, both improving the experience of travelers and reducing the environmental impact. Integrating machine learning methods with climate modeling has potential to speed up simulations and increase our understanding of the climate around us. However, it is surely also possible to find applications of our method with a damaging impact on society, for example through undesired mass-surveillance.

### Future Work

We consider a setting where the graph structure is both known and constant over time. In some practical applications it is not obvious how to construct the graph describing the system. To tackle this problem our model could be combined with approaches for also learning the graph structure (Stankovic et al., 2020; Zhang et al., 2022). Extending our method to dynamic graphs, that evolve over time, would not require any major changes and could be an interesting direction for future experiments.

Our focus has been on forecasting, but the model could also be trained for other tasks. The time-continuous latent state in each node could be used for imputing missing observations or performing sequence segmentation. Also classification tasks are possible, either classifying each node separately or the entire graph-structured time series.

While our model can produce predictions at arbitrary time points, an extension would be to also predict the time until the next observation occurs. One way to achieve this would be to let the latent state parametrize the intensity of a point process (Mei and Eisner, 2017; Jia and Benson, 2019). Building such point process models on graphs could be an interesting future application of our model.

## Acknowledgments

This research is financially supported by the Swedish Research Council via the project _Handling Uncertainty in Machine Learning Systems_ (contract number: 2020-04122), the Excellence Center at Linkoping-Lund in Information Technology (ELLIIT), and the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The computations were enabled by the Berzelius resource at the National Supercomputer Centre, provided by the Knut and Alice Wallenberg Foundation. All datasets available online were accessed from the Linkoping University network.
