Autoregressive GNN-ODE GRU Model for Network Dynamics

Complex networks, as the modeling of complex systems, are ubiquitous in real life, involving different disciplines and domains, such as physics [1], biology [2], sociology [3], computer science [4] and control science [5] . Complex networks are essential carriers of information and have attracted much attention from researchers in recent years [6]. Researchers attempt to utilize various tools to exploit and make good use of the information hidden behind the network. Besides, what interests researchers more is the dynamic process on the complex networks [7], that is, the evolution of the complex systems. Each unit in the network interacts with other units along the edges according to its specific dynamic rules, and the state of each unit evolves over time. Usually, the dynamic rules of each unit can be represented by the ordinary differential equation(ODE) or partial differential equation(PDE). For instance, in the gene regulatory network, the expression level of one gene is affected or restricted by the expression level of other genes, and the expression levels of genes change dynamically [8]. In the mutualistic interaction network, the number of a species is not merely affected by the ecological environment and external factors but also by other partners with mutualistic interactions [7]. The fundamental and essential part of the dynamic process on the network is the governing equation which controls and models the dynamic laws of all the units. Discovering governing equations helps to understand the dynamics of the physical process, which also provides underlying guidance to the advancement of technology in return [9].

However, in many current dynamic systems on networks, governing equations are relatively unknown or partially known, which causes great difficulties in understanding and quantitative descriptions of the system. With the development of technology, states of the dynamic systems are measurable and traceable with instruments and sensors, and abundant time series data are observed, which gives rise to a new mode of data-driven discovery in governing equations or system identification [10]. Earlier researches include classic linear approaches [11], symbolic regression [12], dynamic modeling decomposition [13], nonlinear regression [14], and nonlinear Laplacian spectral analysis [15]. Sparse regression also makes good advancement in this field [16]. Recent efforts are devoted to leveraging artificial neural network (ANN) and deep neural network(DNN) to tackle the time series data and identify the system model [16]. Champion et al. [17] present a deep auto-encoder neural network to model the coordinate transformation and utilize the sparse identification of nonlinear dynamics algorithm for parsimonious modeling. Nevertheless, these methods are susceptible to the system size and are not applicable to the complex dynamic process and large networks composed of thousands of units. In practice, node states of the network can not be observed ideally at regular intervals, and the time series data is usually irregular and sparse in most cases, which also causes problems for these methods. Nowadays, instead of identifying the dynamic model in mathematical functions, researchers try to simulate, fit and learn the network dynamics in the form of a black box without need of identifying or knowing the specific model formulas. NDCN is proposed to use the Ordinary Differential Equation system combined with the Graph Neural Network block to learn the continuous network dynamics process [18]. The GNN block is used to model the complicated and nonlinear interactions between node states, and the real-time dynamic changes of node states are caught or tracked by the ODE system. The real-time states of nodes can be achieved similarly by solving an ODE initial-value problem. The depth of GNN blocks is analogous to the length in the time domain. NDCN uses a single ODE system to model the network dynamics, and the predicted states of nodes are obtained by solving the same ODE initial-value problem. However, using a single ODESolver, NDCN can not handle the long-term network dynamics well. The deviation between the predicted values and real states of nodes may get larger and larger over time if the predicted states of nodes in the early period do not fit well with the true observations. In addition, NDCN does not make full use of all the observations effectively. The observations of node states other than the initial values are not engaged in the training process as prior knowledge. This is probably one of the main reasons for this bad performance in the predictions of the long-term network dynamics. And this situation will become even worse when the observations are very sparse. It is difficult for a single ODE system to learn or fit the real dynamics of node states with sparse observations. To mitigate or overcome the effect of the sparsity of observations on NDCN and improve the ability of the long-term prediction, instead of directly using ODE systems to model the entire network dynamics process like NDCN, we model the network dynamics in segments using the ODE equations in an autoregressive way, which can reduce the difficulty of learning the dynamic process. In this way, the true observations can be fully involved in the training process as prior knowledge and effectively correct the prediction bias at each timestamp to ensure that the initial value for the next prediction is as accurate as possible. By this means, we impose a strong constraint on both ends of the ODESolver to reduce the prediction error and force the ODE system to easily capture or fit the real dynamics.

In this paper, we propose an Autoregressive GNN-ODE GRU Model(AGOG) to learn and capture the continuous network dynamics in a data-driven manner. Considering the sparsity of the observations of node states, we use the Autoregressive model to make a one-step ahead prediction based on the observation history to guarantee the continuity and accuracy of node states. The GNN module is utilized to model the complicated and nonlinear network dynamics, and this module includes two parts that imply the interactions between nodes and the dependence of dynamic evolution, respectively. The dynamics of hidden states of node states is formulated by the ODE system. The augmented ODE system [19] is used to map the hidden state of node states into the continuous time domain, and the hidden state is solved in the expanded space, which could improve generalization and stability and achieve lower losses. And the hidden states are updated through observations by GRUCell. The objective function is defined to minimize the reconstruction error and continuity error which describe the accuracy of prediction and the continuity of dynamic processes, respectively. To evaluate the performance of our model, in consideration of these three tasks: interpolation reconstruction, extrapolation prediction, and regular sequences prediction, we visualize the learned network dynamics and test our model and other baseline methods in different network dynamics with various underlying networks. The results demonstrate that our model can learn and capture the continuous network dynamics precisely and stably. Our model consistently outperforms other baseline methods in most cases and can still achieve comparable performance in other cases.

The rest of this paper is organized as follows. The following Section is the Related Work. Our model is described in Section $3$. The experiments are executed in Sections $4$ and $5$. This paper is summarized in the last Section.

The proposal of Neural ODE [20] causes a sensation and ignites the enthusiasm of researchers. Different from discrete hidden layers in the neural networks, Neural ODE uses the ODE system as a fundamental component, and the output is calculated through the ODESolver. The continuous neural dynamics of Neural ODE determines that it can handle or incorporate time-series data, especially with non-uniform intervals, which inspires further research greatly [21, 22, 23]. Rubanova et al. [24] define two models: autoregressive ODE-RNN and Latent ODE based on variational auto-encode, to handle the time series with irregular intervals. Latent ODE uses ODE-RNN as the recognition network to identify the approximate posterior of the initial state, and the ODESolver is then utilized to achieve the hidden state at any time point. Lechner et al. [25] theoretically prove that ODE-RNN is plagued by gradient vanishing and gradient explosion and introduce ODE-LSTM to address this problem and capture the long-term dependency in irregularly sampled time-series data. Inspired by the Neural ODE, Brouwer et al. [26] modify the Gated Recurrent Unit in the continuous time domain and combine it with the Bayesian update network to model the sporadic observations. Herrera et al. [27] consider the lack of theoretical guarantee for the predictive capabilities of methods mentioned above and then propose Neural Jump ODE to model the conditional expectation for stochastic processes in continuous time.

The powerful capabilities of Graph Neural Networks(GNN) have been proved in various tasks, such as node classification, link prediction, and graph classification [28, 29]. In recent years, multiple variants of GNN spring up such as GCN [30], GAE [31], GAT [32] and GraphSage [33]. GCN [30] can be divided into two classes: spectral-based and spatial-based approaches. Spectral-based approaches define the convolutional operator in the spectral domain, and the convolutional operator is deemed as the noise filter in the graph signal processing [34]. Researchers consider designing different convolutional operators to capture different structural properties of graphs. Spatial-based approaches define message passing rules or convolutional operators based on graph topology, and the node representation or node message is updated by its neighbors’ and its own [35, 36]. To improve the performance of GNN, some tricks or modules are proposed. Skip connection is tried to alleviate the over-smoothing and make GNN deeper [37]. Sampling methods consider reducing the size of node neighbors by conducting sampling for each node or each convolution layer [38]. Pooling methods focus on designing pooling layers to learn high order or graph level hierarchy representation of graphs [39]. Besides, GNN combined with gate mechanisms is also proposed to use LSTM or GRU in the message passing to improve the long-term passing of node information [40].

The work of combining GNN and Neural ODE or PDE has also been proposed. NDCN [18] uses the GNN-based ODE system to learn the nonlinear and complex continuous dynamics on the network and achieve predictions of node states at an arbitrary time. The GNN module is utilized to model the network dynamics equation, and the ODE system is used to map the discrete node states into the continuous time domain. Besides, NDCN is also suitable for the node classification task. The continuous dynamics of ODE corresponds to the discrete layers of GNN, and the continuous time can be interpreted as the depth of discrete layers. Similarly, CGNNs [41] generalize GNN with discrete layers into a continuous physical process. Inspired by the idea of Pagerank and graph diffusion methods, CGNNs use the ODE system to define the continuous dynamics of node embedding, and the derivative of node embedding is derived from the embedding of its neighbors and itself and its initial values. CGNNs can effectively improve over-smoothing, help make GNN deeper, and learn long-term dependency. Grand [42] is also motivated by the continuous diffusion process and model graph learning as a discretization of a basic PDE. And Grand uses different spatial and temporal discretization formalism and diffusivity functions to model the diffusion on the graph. However, CGNNs and Grand are both derived from the graph diffusion mechanisms, designed for the node classification task, and not particularly targeted for modeling the dynamic processes between nodes, and the terminal states or the final representation of nodes are achieved in a short time. Moreover, STGODE [43] applies CGNN to the Traffic Flow Forecasting. STGODE uses the continuous GNN-ODE to capture the spatial-temporal dependency and alleviate over-smoothing.

Consider that a continuous dynamic process on a network is governed by:

where $x(t)$ are the states of nodes and $G=(V,E)$ represents the underlying network structure that describes the interaction relationships between nodes. And $f$ is the governing function which dominates the dynamic law of the evolution of node states. The number of nodes of $G$ is $n$, and $x(t)\in\mathbb{R}^{n\times k}$ where $k$ denotes the dimension of node states. And $k$ can be any integer greater than zero.

For this continuous dynamic process on the network, we can have a series of observations of node states, $\mathbf{x}=\{(t_{0},x_{0}),(t_{1},x_{1}),...,(t_{i},x_{i}),...,(t_{T},x_{T})\}$, where each $t_{i}\in\mathbb{R}$ denotes the timestamp of the observation $x_{i}\in\mathbb{R}^{n\times k}$ and $t_{0}<...<t_{i}<...<t_{T}$. We want to realize predictions of the node states on $G$ at any given time which may be larger or smaller than $t_{T}$. To tackle this problem, we propose an Autoregressive GNN-ODE GRU model to provide a novel solution. Autoregressive models make further predictions based on the observation history and specify the conditional probability $p(x_{i}|x_{i}-1,...,x_{0})$. One-step ahead predictions are achieved by solving an ODE initial-value problem. And the new observations as prior knowledge are combined with the hidden states for the predictions at the next timestamp. The diagram of our model is shown in Fig. 1, and the framework of our model is specified as follows.

Suppose that the initial observation of node states is $x_{0}\in\mathbb{R}^{n\times k}$, and the corresponding observation timestamp is $t_{0}$. We first map $x_{0}$ into a hidden state $h(t_{0})\in\mathbb{R}^{n\times d}$:

in which $W_{e}\in\mathbb{R}^{k\times d}$ and $b_{e}\in\mathbb{R}^{n\times d}$.

The dynamic of hidden states is specified by the following ordinary differential equation system:

in which the function $g_{\theta}(h(t),t)$ formulates the latent dynamic of hidden states with the parameter $\theta$. Consider that the dynamic process is modeled or relies on the network structure. The function $g_{\theta}$ is designated expressly by a GNN module:

where $\Phi$ is the normalized Laplacian matrix, $\Phi=D^{-\frac{1}{2}}(D-A)D^{-\frac{1}{2}}$ in which $A$ denotes the adjacency matrix of $G$ and $D$ denotes the degree matrix of the network. On the right side of Eq. 4, the left part models the interactions of hidden states between nodes, while the right part specifies the dependence of the dynamic evolution of the hidden state on the current hidden state of nodes. The left part simulates the influence of dynamic interactions between nodes on the state of nodes, and the right part represents the self-dynamic of nodes which implies the evolution rule of the node states influenced by their own states.

So, the next hidden state at $t_{1}$, $h^{\prime}(t_{1})$, can be achieved by the ODESolver:

This initial-value problem is solved by numerical methods, e.g., the fourth-order Runge-Kutta method, the high-order Dormand-Prince DOPRI5 method, and the Euler method. We use the ODE solvers from the torchdiffeq package (https://github.com/rtqichen/torchdiffeq), and details can be found in  [20]. To make the ODESolver more powerful and expressive, we adopt the Augmented Neural ODE [19] to offer a novel solution to the problems above. The space of the hidden state is expanded, and the ODE is solved from $\mathbb{R}^{n\times d}$ to $\mathbb{R}^{n\times(d+p)}$ where $p$ represents the augmented dimension. The augmented space is defined as $a(t)\in\mathbb{R}^{n\times p}$, and Eq. 3 is updated to:

For simplicity, the hidden state with the augmented space is denoted as $h_{a}(t)$, and $h_{a}(t_{0})=[h(t_{0}),\mathbf{0}]$. Intuitively, the hidden state $h(t)$ is merged with the augmented zero space $a(t)$, and the ODE is solved in the expanded space. With the additional space, the ODE flow can easily flow to the augmented space, and the learned hidden state is represented in a higher dimension and has a better expressivity and precision. Simultaneously, the trained function $G_{\theta}$ becomes more smooth and achieves better generalization and lower time cost. Then the predicted hidden state $h^{\prime}(t_{1})\in\mathbb{R}^{n\times(d+p)}$ at $t_{1}$ is achieved by:

And the prediction of node states at $t_{1}$, $x^{\prime}_{1}$, is achieved through a linear layer:

where $W_{o}\in\mathbb{R}^{(d+p)\times k}$ and $b_{o}\in\mathbb{R}^{n\times k}$.

The real node states at $t_{1}$ are denoted as $x_{1}$. We use real observations as the prior knowledge to revise and update the obtained hidden state, $h^{\prime}(t_{1})\in\mathbb{R}^{n\times(d+p)}$. We also map $x_{1}$ into the hidden state as $h(t_{1})\in\mathbb{R}^{n\times d}$ and utilize the GRUCell to update hidden state $h^{\prime}(t_{1})$ at $t_{1}$ given current true observation $x_{1}$:

In detail, Eq. 11 is formulated as:

in which $W_{h}\in\mathbb{R}^{(d+p)\times d}$, $U_{h}\in\mathbb{R}^{d\times d}$, $b_{wh},b_{uh}\in\mathbb{R}^{n\times d}$. The updated hidden state at $t_{1}$ is denoted as $\hat{h}(t_{1})\in\mathbb{R}^{n\times d}$. The augmented hidden state $h_{a}(t_{1})$ at $t_{1}$ is expressed as $[\hat{h}(t_{1}),\mathbf{0}]$, and the updated prediction of node states at $t_{1}$ from the augmented updated hidden state $h_{a}(t_{1})$, can be calculated from:

Then $h_{a}(t_{1})$ can be the initial value of the ODESolver to get the hidden state at the next moment $t_{2}$, $h^{\prime}(t_{2})$:

And the prediction $x^{\prime}_{2}$ at $t_{2}$ can be achieved through $h^{\prime}(t_{2})$. Repeating these procedures above, we can obtain the prediction of node states at the corresponding time. There are two types of outputs, $x^{\prime}_{i}$ and $\hat{x}_{i}$, defined as the prediction and the updated prediction to show the distinction, respectively. $x^{\prime}_{i}$ is obtained from the solution of the ODESolver $h^{\prime}(t_{i})$, through a linear layer based on Eq. 9, while $\hat{x}_{i}$ is achieved from the augmented updated hidden state $h_{a}(t_{i})$, through a linear layer based on Eq. 13. The augmented updated hidden state $h_{a}(t_{i})$ is obtained in two steps. The predicted hidden state of the ODESolver solution $h^{\prime}(t_{i})$ is updated by the GRUCell through the real observations $x_{i}$ at the same time to get the updated hidden states $\hat{h}(t_{i})$, then $h_{a}(t_{i})$ is achieved after the augmentation of the $\hat{h}(t_{i})$. And $x^{\prime}_{i}$ is deemed as the final result that we need.

The loss function is defined as below:

where $|\cdot|$ denotes the $l_{1}$-normal loss, that is the mean element-wise absolute value difference. Our aim is to minimize this objective function, Eq. 15. The above function has two terms: the left part is the reconstruction error, while the right part is the continuity error. It is easy to understand the reconstruction error which requires the predicted value to be closer to the true states. As for the continuity error, in the dynamic of the hidden states, we expect that the updated hidden states should be consistent with the predicted hidden states from the ODESolver at the same timestamp. The updated hidden states are then used as the initial value for the predictions at the next timestamp. We use the true observations to correct the hidden state and hope that the updated hidden state will not deviate too much from the original hidden state. By introducing the true observations, we impose a strong constraint on both ends of the ODESolver and force the ODE system to capture the real dynamics easily. Moreover, the continuity error implies the constraints on the smoothness and continuity of the trajectory of the predicted states. $W_{e},W_{o},\theta_{0},\theta_{1},b_{e},b_{o},b_{0}$, and GRUCell are all trainable parameters, and they are updated through backpropagation of the gradient of the loss. The loss function is optimized by the Adam method, and the pseudocode of our Autoregressive GNN-ODE GRU model is presented in Algorithm 1. For clarity, our model is abbreviated as AGOG(Autoregressive GNN-ODE GRU Model).

We evaluate our model on different structures of networks with various network dynamics. Complex network models include: 1) Community network [45], 2) Grid network, 3) Erdós and Rényi random network [46], 4) Power-law network [47], and 5) WS small-world network [48]. For each type of network models, we investigate three kinds of network dynamics processes. And the real benchmark dynamic of node states is generated according to its underlying network structure and rules of dynamic evolution by the Dormand-Prince method. Regarding the experimental setup, we generally follow the settings in NDCN [18], except for setting a larger time interval for observations, because we consider the sparsity of observations in real life. In the simulation, for each network dynamics on each network structure, we randomly sample 120 snapshots of node states with irregular time interval, $\mathbf{x}=\{(t_{0},x_{0}),(t_{1},x_{1}),...,(t_{119},x_{119}),t_{0}<t_{1}...<t_{119}\}$, where sampling intervals between adjacent datapoints are not equal. The last $20$ datapoints $\{(t_{i},x_{i})\}_{i=100,...,119}$ are used for the extrapolation prediction task. As for the first 100 datapoints, we randomly choose $P\%$ datapoints for training, and the residual datapoints are targeted for the interpolation reconstruction task. A smaller $P$ means that the observed data is more sparsity. To illustrate the superiority of our model, we use the MAE(Mean Absolute Error), and Normalized L1 loss metrics to measure the performance of our model on these two tasks. The observed node states may have zero, so the MAPE(Mean Absolute Percentage Error) is not applicable here. Instead, the Normalized L1 loss is selected, which is normalized by the mean value of all observations.

Experimental configuration. The number of nodes of the generated networks is set as $400$. The initial node state $x_{0}$ is set as the same for each case. The training percentage of observed snapshots $P\%$ is set as $10\%$, and $30\%$, respectively. With regard to hyperparameters in our model, the hidden dimension $d$ of hidden state $h(t_{i})$ is set as $20$, $h(t_{i})\in\mathbb{R}^{n\times 20}$, and the augmented dimension $p$ is set as $5$. And our model is trained for $800$ epochs with the learning rate $0.01$ in all cases. The Euler method is used for the ODESolver.

We use NDCN method [18] as the baseline to compare with our model. NDCN uses a simple GNN module to model the interactions between node states, and the discrete node states are mapped into the continuous time domain by the ODE system. The main goal of NDCN is state prediction on the network. There are also three variants of NDCN: No-Encode, No-Graph, and No-Control. For succinctness, we only consider NDCN in this paper because of its best performance between it and its variants. Moreover, there are also some GNN-ODE methods that are not designed for predictions of node states, but we consider them as baselines here as well. These methods define the evolution function of node embedding as the ODE function, and the final embedding is achieved by the ODESolver in a finite time. Here we regard the node embedding as the hidden state of nodes states. We use the ODE functions of these GNN-ODE methods to model the network dynamics, respectively. CGNNs [41] expand the discrete GNN layers into a continuous physical process inspired by the diffusion process on graphs and define an evolution equation of node embedding that the derivative of node embedding is the combination of itself, the embedding of its neighbors, and its initial value. CGNNs have two different variants, which consider whether each dimension of node embedding evolves independently or not, and we denote them as CGNN and WCGNN, respectively, which are all considered here. Similarly, Grand [42] is also motivated by the continuous diffusion process and model GNN as a discretization of a basic PDE. And Grand uses different spatial or temporal discretization formalisms and diffusivity functions to model the diffusion on the graph. The diffusivity is modeled by the scaled dot product attention mechanism. The more general non-linear form of GRAND [42] is tested here, denoted as GRAND. We use the default parameters for NDCN as stated in  [18]. The hidden size is set as $20$ for CGNN, WCGNN, and GRAND. And the dimension of the query and key matrices is both set as $20$ for GRAND. The loss function of CGNN, WCGNN, and GRAND is the same as that of NDCN.

Under ideal conditions, the node states can be observed at equal time intervals. In this section, we try to use our model to capture the network dynamics based on regular sequences of node states. In the experiment, we also consider these continuous network dynamics with different underlying networks above. Regarding the experimental setup, we generally follow the settings in NDCN [18], except for setting a larger time interval for observations, because we consider the sparsity of observations in real life. For each case, we sample 80 snapshots of node states with equal intervals, $\mathbf{x}=\{(t_{0},x_{0}),(t_{1},x_{1}),...,(t_{79},x_{79}),$ $\{t_{0}<...<t_{79}\}$, where time intervals between adjacent observations are equal. The first $80\%$ observations are used for training, while the last $20\%$ are considered as the test set for the regular sequences prediction task. For evaluation, we also employ the MAE(Mean Absolute Error) and Normalized L1 loss metrics to validate the performance of methods.

Experimental configuration. The basic experimental settings are set as the same as that in the previous Section. Consider the time interval between observations is equal. Beyond the baseline methods mentioned above, we also choose temporal-GNN algorithms as the baseline, which are also considered in  [18]. Temporal-GNN algorithms are the combinations of GCN and RNN blocks, and they can not be applied to the interpolation reconstruction task. Only the extrapolation prediction task is considered here. The hidden size of the GCN and RNN blocks are set as $10$ and $5$, respectively. The default parameters are used for NDCN, and the parameters of our model and other GNN-ODE based methods remain the same as in Section $4$.

Revealing the continuous dynamics on complex networks is essential for understanding complex systems. But this task is very hard due to the relatively unknown and partially known equations and high dimensions of complex systems. In addition, the observations of node states are generally non-uniform and very sparse in real life, which also poses a great challenge. In this paper, we propose an Autoregressive GNN-ODE GRU Model(AGOG) to capture the continuous network dynamics based on the observed node states in a data-driven manner. The further predictions are based on the observation history. With the predictions achieved, the true observations at the same time as the prior information are used to update the hidden states for the predictions at the next step. We impose a strong constraint to make the observations into the training process and make the ODE system learn the real continuous network dynamics more easily and accurately. We use a GNN module to simulate and fit the network dynamics, and this module has two parts that describe the interactions between nodes and the dependence of the dynamic evolution, respectively. The hidden state of nodes is specified by the ODE system, and we utilize the augmented ODE to map the hidden state into continuous time. The augmented ODE solves the hidden states in the augmented space, which could achieve lower losses and improve generalization and stability. Then the hidden states are updated through observations by GRUCell. To evaluate the validity of our model, we visualize the learned network dynamics by our model and test it in three tasks: interpolation reconstruction, extrapolation prediction, and regular sequences prediction. The results show that our model can consistently outperform other baseline methods or achieve comparable performance in these dynamics with different underlying networks. Our model can accurately learn the continuous network dynamics with a small error. In future work, we consider introducing the attention mechanism into the learning of network dynamics.