Missing Data Estimation in Temporal Multilayer Position-Aware Graph Neural Network (TMP-GNN)

Abstract

GNNs have been proven to perform highly effectively in various node-level, edge-level, and graph-level prediction tasks in several domains. Existing approaches mainly focus on static graphs. However, many graphs change over time and their edge may disappear, or the node/edge attribute may alter from one time to the other. It is essential to consider such evolution in the representation learning of nodes in time-varying graphs. In this paper, we propose a Temporal Multilayer Position-Aware Graph Neural Network (TMP-GNN), a node embedding approach for dynamic graphs that incorporates the interdependence of temporal relations into embedding computation. We evaluate the performance of TMP-GNN on two different representations of temporal multilayered graphs. The performance is assessed against the most popular GNNs on a node-level prediction task. Then, we incorporate TMP-GNN into a deep learning framework to estimate missing data and compare the performance with their corresponding competent GNNs from our former experiment, and a baseline method. Experimental results on four real-world datasets yield up to $58\%$ lower $\mathrm{ROC}\mathrm{AUC}$ for the pair-wise node classification task, and $96\%$ lower MAE in missing feature estimation, particularly for graphs with a relatively high number of nodes and lower mean degree of connectivity.

1. Introduction

Graph neural networks (GNN) have recently been used as a standard in developing machine learning methods for graphs. These graphs have been formed from sources such as transportation networks [1,2,3,4], brain networks [5], social media community networks, etc. The GNN architecture has effectively combined the node/edge features and graph topology to build a distributed representation. The resulting representation can be used to solve node-level, edge-level [6] and graph-level prediction tasks [5]. Several recent works have been applied to represent nodes, edges and graphs in many domains, such as bioinformatics, natural language processing and transportation. In [7], the authors used a hybrid method of a decomposed LSTM and a GCN to predict next period prescription. The authors in [8] proposed a GCN-based Relational Fusion Network (RFN) for speed limit classification and speed estimation in road networks. The study introduced the notion of volatile homophily in the road segments of adjacent regions and enumerated the substantial differences between the road networks and other networks previously learned by GCN. The underlying differences include a smaller number of node attributes, potential abrupt behavioral changes among the neighboring segments and a lower mean node degree. Another study [9] learned the graph representation by applying LSTM autoencoders on graph sequences generated from random walks, and utilized it for the classification task.

The goal of node embedding methods is to identify a vector representation that captures the node location within a broader topological structure of the graph. Most node embeddings learned from GNN architectures focus on single static graphs. These methods assume that the number/position of nodes as well as their interactions do not change over time. However, in many applications, we tackle spatio-temporal measurements collected over time, wherein the necessity of time-varying graphs emerges. In such graphs, the number of nodes, their connecting edges and the edge weights vary from time to time; thus, dynamic node embedding should be learned accordingly from the graph. Node embedding approaches in such graphs can be mainly categorized as temporal-first [10,11] and structural-first [12,13] according to [14], in which the former and the latter prioritize learning temporal and structural information in their framework first and then feed static topological and dynamic temporal features into their model, respectively. Time-varying graphs can be defined as continuous- or discrete-time dynamic graphs, in which the latter can be represented as a sequence of interdependent time layers. Each layer is a graph building from existing nodes and weighted edges corresponding to a given time. For instance, the study in [5] proposed a graphSAGE-based graph representation learning for time-variant EEG signals to apply brain state classification, where the dynamic graph is modeled as a sequence of individual graphs per one time stamp. Another work [8] combines temporal, spatial and semantic views to predict taxi demands at a future time stamp.

One straightforward way of computing node embedding in discrete time-varying graphs is to use a static GNN-based node embedding for each individual time layer and aggregate the results of corresponding layers through recurrent neural network variants. However, this method would compute an embedding for each layer independently, and consequently ignore the inter-layer correlation. In this paper, we present TMP-GNN, a Temporal Multilayer Position-Aware GNN-based node embedding, which is an extension of its static version, Position-Aware GNN (P-GNN) [15]. The goal of P-GNN is to learn position-aware node embeddings that utilize the local network structure and the global network position of a given node with respect to randomly selected nodes called anchor sets, which enables us to distinguish among isomorphic nodes—nodes that are based in very different parts of the graph but have topologically the same structure. The distinguishing power will improve if the node/edge features are available. The resulting embeddings can be later used to approximately calculate the shortest path distance among the embedded nodes in the graph. The contributions of the paper are summarized as follows:

(1) We learn the short-term temporal dependencies, global position and feature information of the graph jointly through our TMP-GNN embedding component and utilize the derived representation in a missing data estimation framework.

(2) Instead of modeling a dynamic graph via multi-graph, we exploit a supra-adjacency matrix to encode a temporal graph with its intra-layer and inter-layer coupling in a single graph, which facilitates faster learning and a higher area under the receiver operating characteristic curve $(\mathrm{ROC}\mathrm{AUC})$ for pair-wise node classification.

(3) We deploy the concept of conditional centrality derived from eigenvector-based centrality to distinguish nodes of higher influence and integrate it in message aggregation across the graph.

(4) We identify a new experimental setting to enhance the training of the multi-graph deployment of dynamic graphs so that it can better trace the behavioral changes of nodes in the graph.

(5) We use hidden states learned from bi-directional GRU (bi-GRU) to learn the long-term temporal dependencies in both forward and backward directions to estimate missing values.

(6) We conduct several experiments using four real-world datasets, with a wide range of node numbers, degrees of connectivity, edge dynamics and areas of study. The results illustrate that TMP-GNN improves the area under the ROC curve significantly for classification tasks as compared to the best baseline, and it reduces the MAE of missing data estimation.

2. Notations and Preliminaries

Before presenting the details of our proposed architecture, we provide some background on temporal multilayer graphs.

2.1. Notation

Figure 1 illustrates a temporal multilayer graph. We can represent the graph as ${\mathcal{G}}^{\left(t\right)}=(\mathcal{V},{\rho }^{(\mathcal{E}\times t)},{\mathcal{E}}^{\left(t\right)},\tilde{\mathcal{E}})$, where $\mathcal{V}$ are the set of nodes $\{v,\forall v\in \mathcal{V}\}$. Here, ${\rho }^{(\mathcal{E}\times t)}$ indicate whether a given edge is present at a given time t. ${\mathcal{E}}^{\left(t\right)}$ shows intra-layer weighted edges at time t, ${\mathcal{E}}^{\left(t\right)}=\{e_{uv}^{\left(t\right)},\forall u,v\in \mathcal{V}\}$. $\tilde{\mathcal{E}}$ indicates inter-layer edges $\{(t,t^{\prime },{\tilde{x}}_{t{t}^{\prime }})\in \tilde{\epsilon }\}$, where ${\tilde{x}}_{t{t}^{\prime }}$ is the edge weight between time layers $(t,t^{\prime }).$ The multilayer graph is also denoted by a sequence of adjacency matrices ${\mathit{A}}^{\left(\mathbf{t}\right)}\in {\mathbb{R}}^{N\times N}$, where $A_{uv}^{\left(t\right)}$ denotes the directed edge from node u to v in time layer t, as follows:

$$A_{uv}^{\left(t\right)}=\left\{\begin{array}{cc}{x}_{uv}^{\left(t\right)},\hfill & \mathrm{if}u,v \mathrm{are} \mathrm{connected} \mathrm{at} \mathrm{time} t\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(1)

where $x_{uv}^{\left(t\right)}$ is the feature vector associated with $e_{uv}^{\left(t\right)}$. In scenarios where the node attributes are available, another set is defined as $W^{\left(t\right)}=\{w_v^{\left(t\right)},\forall t\in {\mathbb{N}}_T\}$, where $w_v^{\left(t\right)}$ is the feature vector associated with node v in time layer t.

The set $\tilde{\mathcal{E}}$ encodes the weight and topology of the coupling of individual instances of the same nodes between pairs of time layers $(t,t^{\prime })$.

The inter-layer coupling described above is diagonal and uniform [16], which means that the existence of inter-layer edges is restricted between separate instances of the same nodes from one layer to another. For instance, there is no edge from node v in t to node u in $t^{\prime }$. The corresponding layers are uniformly coupled, meaning that the inter-layer edge weight between two layers is identical for all nodes in those layers. Please note that this assumption can be generalized so that different inter-layer edge weights are assigned to different subsets of nodes; that is, ${\tilde{x}}_{t{t}^{\prime }}^u\ne {\tilde{x}}_{t{t}^{\prime }}^v$.

There are two possible ways to form the coupling between the layers: directed and undirected chain. The former couples the adjacent time layers by neglecting the directionality of time as in the following equation:

$${\tilde{A}}_{uv}^{\left(t{t}^{\prime }\right)}=\left\{\begin{array}{cc}\hfill 1& \hfill \forall u,v\mathrm{if}u=v,|t^{\prime }-t|=\delta \\ \hfill 0& \mathrm{otherwise}\hfill \end{array}\right.$$

(2)

where $\delta$ is the sampling rate. The latter not only respects the time direction, but also includes weighted undirected connections between the same nodes of all pairs of time layers as in:

$${\tilde{A}}_{uv}^{\left(t{t}^{\prime }\right)}=\left\{\begin{array}{cc}\hfill 1+\gamma ,\forall u,v& \mathrm{if}u=v,|t^{\prime }-t|=\delta \hfill \\ \gamma \hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(3)

where $\tilde{\mathit{A}}$ is a $T\times T$ inter-layer matrix, and $\gamma$ is the node teleportation probability [17]. In this paper, we have used the former approach for inter-layer coupling. The inter-layer coupling of the same nodes decrease as $|t^{\prime }-t|$ increases, and we adjust $\delta$ according to the data characteristics and inter-sequence correlation in our experiments.

A node embedding function, e.g., SAGE, can be represented as $f:{\mathcal{V}}^{\left(t\right)}\to {\mathcal{Z}}^{\left(t\right)}$, which maps a given node v to d-dimensional vector $z$ in layer t, where ${\mathcal{Z}}^{\left(t\right)}=\{{\mathit{z}}_v^{\left(t\right)},\forall v,t\in {\mathbb{N}}_{\left|\mathcal{V}\right|},{\mathbb{N}}_T\},{\mathit{z}}_v^{\left(t\right)}\in {\mathbb{R}}^d$. Once the node embedding is computed, the corresponding edge embedding can be calculated as ${\mathit{z}}_{uv}^{\left(t\right)}=g({\mathit{z}}_u^{\left(t\right)},{\mathit{z}}_v^{\left(t\right)})$, where the g function is equal to mean in the case of our study.

2.2. Supracentrality Matrix for Temporal Multilayer Position-Aware Graph

As explained, a temporal graph is represented through a sequence of adjacency matrices, each of which refers to one layer of a dynamic network at a specific point of time. Then, we construct a supracentrality matrix $\mathbb{C}\left(\omega \right)$ by linking centrality matrices across time layers through a weighted inter-layer parameter called $\omega$, which is used to adjust the extent of coupling strength among pairs of time layers. The entries of the dominant eigenvector of $\mathbb{C}\left(\omega \right)$ show joint centrality, the importance of every node–layer pair $(v,t)$. Additionally, marginal and conditional centralities are defined to represent the relative importance of a node compared to other nodes at time layer t.

The supracentrality framework is mainly focused on eigenvector-based centrality, which is obtained by calculating the centralities as the elements of the dominant eigenvector corresponding to the largest-magnitude eigenvalue of a centrality matrix $\mathcal{C}\left(A\right)$, which is defined as a function of network adjacency matrix $A$. We have selected a centrality measure from one of the most popular choices to be equal to adjacency matrix ($\mathcal{C}\left(A\right)=A$).

$\mathbb{C}\left(\omega \right)$, a supracentrality matrix for a dynamic graph, is a group of matrices formed as below:

$$\mathbb{C}\left(\omega \right)=\widehat{\mathbb{C}}+\omega \widehat{\mathbb{A}},$$

(4)

where

$$\begin{array}{c}\hfill \widehat{\mathbb{C}}=\left[\begin{array}{ccccc}{\mathbf{C}}^{\left(1\right)}& 0& 0& \cdots & 0\\ 0& {\mathbf{C}}^{\left(2\right)}& 0& \cdots & 0\\ 0& 0& {\mathbf{C}}^{\left(3\right)}& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& 0& 0& \cdots & {\mathbf{C}}^{\left(T\right)}\end{array}\right],\end{array}$$

and

$$\begin{array}{c}\hfill \widehat{\mathbb{A}}=\left[\begin{array}{ccccc}{\tilde{A}}^{\left(11\right)}\mathit{I}& {\tilde{A}}^{\left(12\right)}\mathit{I}& {\tilde{A}}^{\left(13\right)}\mathit{I}& \cdots & {\tilde{A}}^{\left(1T\right)}\mathit{I}\\ {\tilde{A}}^{\left(21\right)}\mathit{I}& {\tilde{A}}^{\left(22\right)}\mathit{I}& {\tilde{A}}^{\left(23\right)}\mathit{I}& \cdots & {\tilde{A}}^{\left(2T\right)}\mathit{I}\\ {\tilde{A}}^{\left(31\right)}\mathit{I}& {\tilde{A}}^{\left(32\right)}\mathit{I}& {\tilde{A}}^{\left(33\right)}\mathit{I}& \cdots & {\tilde{A}}^{\left(3T\right)}\mathit{I}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ {\tilde{A}}^{\left(T1\right)}\mathit{I}& {\tilde{A}}^{\left(T2\right)}\mathit{I}& {\tilde{A}}^{\left(T3\right)}\mathit{I}& \cdots & {\tilde{A}}^{\left(TT\right)}\mathit{I}\end{array}\right].\end{array}$$

Here, $\mathbb{C}\left(\omega \right)$ consists of two components in (4): $\widehat{\mathbb{C}}$ and $\omega \widehat{\mathbb{A}}$. The former (diagonal) component, $Diag[{\mathbf{C}}^{\left(1\right)},\cdots ,{\mathbf{C}}^{\left(T\right)}]$, represents a set of T weighted centrality matrices of individual layers, and the latter (off-diagonal block), $\widehat{\mathbb{A}}=\tilde{\mathit{A}}\otimes \mathit{I}$, encodes the uniform and diagonal coupling with strength parameter $\omega$ between the time layers, where $\tilde{\mathit{A}}$ is defined in (2); ⊗ is the Kronecker product, and $\mathit{I}$ is a $N\times N$ identity matrix, since we only consider the coupling among the same nodes between consecutive pairs of layers. We study the dominant eigenvector $\mathbb{V}\left(\omega \right)$ of $\mathbb{C}\left(\omega \right)$ corresponding to the largest eigenvalue ${\lambda }_{max}$ as in:

$$\mathbb{C}\left(\omega \right)\mathbb{V}\left(\omega \right)={\lambda }_{max}\left(\omega \right)\mathbb{V}\left(\omega \right).$$

(5)

The elements in the dominant right eigenvector of $\mathbb{V}\left(\omega \right)$ are interpreted as scores that measure the importance of node–layer pairs $(v,t)$. As such, the eigenvalue entity ${\mathrm{W}}_{v,t}\left(\omega \right)$ in the following Equation (6), called the joint centrality, is a d-dimensional vector, reflecting the centrality of node v at layer t, where ${\mathbb{V}}_{N(t-1)+n}\left(\omega \right)$ is the $N(t-1)+n$-th entity of the largest eigenvalue $\mathbb{V}\left(\omega \right)$ and n refers to node order $v_n$, which is omitted here for simplicity.

$${\mathrm{W}}_{v,t}\left(\omega \right)={\mathbb{V}}_{\mathrm{N}(t-1)+n}\left(\omega \right).$$

(6)

Inspired by probability theory, the authors in [16] defined marginal layer centrality $\left(\mathrm{MLC}\right)$ and conditional centrality in the following Equations (7) and (8). $\mathrm{MLC}$ indicates the average joint centralities over all nodes at time layer t. The conditional centrality of node v shows the importance of the node relative to other nodes at layer t.

$${\mathrm{MLC}}_v\left(\omega \right)=\sum _{v=1}^{\left|\mathcal{V}\right|}{\mathrm{W}}_{v,t}\left(\omega \right),$$

(7)

$$C{C}_v\left(\omega \right)=\frac{{\mathrm{W}}_{v,t}\left(\omega \right)}{{\mathrm{MLC}}_v\left(\omega \right)}.$$

(8)

Choice of $\omega$ for Steady-State Node Ranking

A number of factors should be taken into account while choosing the value of $\omega$. The choice of $\omega$ should be appropriate for the targeted application or the research question.

When $\omega$$\to \infty$, the centrality measures change slowly over time. On the other hand, in the presence of small $\omega$, the $C{C}_v\left(\omega \right)$ fluctuates from one time layer to the other.

We considered $\omega$ at a strong coupling regime. If $\omega$ is too large (e.g., $\omega$ > 100), the $C{C}_v\left(\omega \right)$ becomes stationary, which is equal to the average $C{C}_v\left(\omega \right)$ over all time layers [18]. Moreover, the dominant eigenvalue of $\mathbb{C}\left(\omega \right)$ converges to the same values of $\widehat{\mathit{A}}$. That is, ${\lambda }_{max}\left(\omega \right)\to {\tilde{\mu }}_1$, where ${\lambda }_{max}$ and ${\tilde{\mu }}_1$ are the dominant eigenvalues of $\mathbb{C}\left(\omega \right)$ and $\widehat{\mathit{A}}$, respectively.

Through the power iteration method of eigenvector measurements [19], we found the minimum value of $\omega$ that satisfies the above requirements (Algorithm 1). In the next section, we utilize stationary ${CC}_v\left(\omega \right)$ to distinguish highly important nodes belonging to P-GNN anchor sets, and then exploit it to aggregate the information across the anchor sets.

Algorithm 1 Find the appropriate value of $\omega$

Require: Initialize the value of $\omega$ for $\omega \ge 0$ (we start from $\omega =10$), construct $\mathit{B}=\omega \widehat{\mathbb{A}}$,

$\mathit{C}=\mathbb{C}\left(\omega \right)$

choose a random number for $b_k$, $c_{k^{\prime }}$ as the initial value for the largest eigenvector of $\mathit{B}$

and $\mathit{C}$. We set $b_k$, $c_{k^{\prime }}$ = $[11\cdots 1]$ with its length equal to the number of columns of $\mathit{B}$, $\mathit{C}$,

which is $N\times T$ in our case.

$b_k\leftarrow \frac{b_k}{\left|\right|b_k\left|\right|}$, $c_{k^{\prime }}\leftarrow \frac{c_{k^{\prime }}}{\left|\right|c_{k^{\prime }}\left|\right|}$

${\lambda }_k\leftarrow b_k^{*}\mathit{B}{b}_k$, ${\mu }_{k^{\prime }}\leftarrow c_k^{*}\mathit{C}{c}_k$

k = 0

while$|{\lambda }_k-{\mu }_{k^{\prime }}|<0.01$ do

while $|{\lambda }_{k+1}-{\lambda }_k|<0.01$ do

$k\leftarrow k+1$

$b_k=\frac{\mathit{B}{b}_k}{\left|\right|\mathit{B}{b}_k\left|\right|}$

${\lambda }_k=b_k^{*}\mathit{B}{b}_k$

end while

return ${\lambda }_k$, $b_k$

while $|{\mu }_{k^{\prime }+1}-{\mu }_{k^{\prime }}|<0.01$ do

$k^{\prime }\leftarrow k^{\prime }+1$

$c_{k^{\prime }}=\frac{\mathit{C}{c}_{k^{\prime }}}{\left|\right|\mathit{C}{c}_{k^{\prime }}\left|\right|}$

${\mu }_{k^{\prime }}=c_{k^{\prime }}^{*}\mathit{C}{c}_{k^{\prime }}$

end while

return ${\mu }_{k^{\prime }}$, $c_{k^{\prime }}$

end while

return $\omega$

2.3. TMP-GNN: Multilayer Position-Aware-Based Graph Neural Network

Several GNN methods are distinguished by the way in which they aggregate nodes’ information from their neighborhood to form a node representation at each layer k. The approach can be summarized as two functions, AGGREGATE and COMBINE [20], in which the former decides on how the information is aggregated from adjacent nodes, and the latter updates the node representation from the previous layer $k-1$ to layer k. The underlying functions can be combined and jointly represented for some GNNs, as follows:

$$\begin{array}{cc}\hfill {\mathit{a}}_{v_n}^{\left(k\right)}& ={\mathrm{AGGREGATE}}^{\left(k\right)}(\{{\mathit{h}}_{u_n}^{(k-1)}:u_n\in \mathrm{neighbors}\left(v_n\right)\}),\hfill \\ \\ \hfill {\mathit{h}}_{v_n}^{\left(k\right)}& ={\mathrm{COMBINE}}^{\left(k\right)}({\mathit{h}}_{v_n}^{(k-1)},{\mathit{A}}_{v_n}^{\left(k\right)}),\hfill \end{array}$$

(9)

where $h_v^{\left(k\right)}$ is the feature representation vector of node v at iteration k/k-th layer of GNN, and $\mathrm{Neighbors}\left(v\right)$ is the set of nodes adjacent to v in $\mathcal{G}$.

We use a position-aware GNN, which, instead of aggregating the information from the nearest neighbors, aggregates the positional and feature information of each node with a randomly selected number of nodes called anchor sets. Then, the computed message is aggregated across the anchor sets for each node to incorporate the global topological information of the graph in node embedding. Figure 2 illustrates in detail how the node representation is learned through one P-GNN layer. The goal is to find the best position-aware embedding ${\mathit{z}}_v^t$ with minimum distortion for a given node v at time layer t.

We have made the following modifications to P-GNN:

• Generalization of P-GNN to Time-Varying Graphs:
  We adopt the input of P-GNN as supracentrality matrix $\mathbb{C}\left(\omega \right)$, which represents a temporal multilayer graph with $N\times T$ number of nodes. We compute an embedding for all nodes in all time layers, since $e_{uv}/x_{uv}/w_{uv}$ defined in (1) can change from t to $t^{\prime }$. The embedding ${\mathit{z}}_v^t$ will then be aggregated from an RNN-based representation to estimate missing data.
• Modification of Message Computation Function (F):
  In an Intelligent Transportation System (ITS), the average speed of neighboring nodes might correlate more or less at different time layers due to a variety of factors, i.e., different types of residential zone, special events, accidents, etc. Using an attention mechanism while computing an anchor set’s message with respect to a given node v can alleviate misinformative messages from the anchor set to influence node v’s embedding. Therefore, we use the attention mechanism to learn the relative weights between the feature vector of v and its nearest neighbor from the anchor set. The fact at different degrees could apply to other application domains as well. As such, we modify our message computation function to incorporate the attention mechanism.
  In P-GNN, we have multi-level aggregations, as demonstrated in Figure 2. First, the message of each node v is computed with regard to each anchor set $R_j$ through function F. From each anchor set, only the nodes with an up to two-hop distance to node v are considered for message computation, as shown in the following equation:

  $$\begin{array}{cccc}\hfill {\mathcal{M}}_j& =\left\{F(u,v,{\mathit{h}}_u^{k-1},{\mathit{h}}_v^{k-1})\right|u\in R_j,R_j\subset \mathcal{R},\mathcal{R}\sim \mathcal{V},\hfill & & \hfill \\ \hfill d(v,R_j)& ={min}_{u\in R_j}d(u,v)\le 2\},\hfill & & \hfill \\ \hfill {\mathit{M}}_v^k\left[j\right]& =\mathrm{MEAN}\left({\mathcal{M}}_j\right),\hfill \end{array}$$

  (10)

  where ${\mathit{M}}_v^k\left[j\right]$ indicates the aggregated message of v corresponding to $R_j$, which is the average of individual outputs of $F(u,v,{\mathit{h}}_u^{k-1},{\mathit{h}}_v^{k-1})$, where F is defined as:

  $$\begin{array}{c}\hfill F(u,v,{\mathit{h}}_u^{k-1},{\mathit{h}}_v^{k-1})=\frac{1}{d_{sp}^q(u,v)+1}{\mathit{a}}_v,\end{array}$$

  (11)

  where $d_{sp}^q(u,v)$ is defined as below:

  $$\begin{array}{c}\hfill d_{sp}^q(u,v)=\left\{\begin{array}{cc}{d}_{sp}(u,v)\hfill & \mathrm{if}{d}_{sp}(u,v)\le q,\hfill \\ \infty ,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

  (12)

  where $d_{sp}(u,v)$ is the shortest path between a pair of nodes $(u,v)$, and we write ${\mathit{a}}_v=w_c\left(\mathrm{CONCAT}({\mathbf{c}}_v,{\mathit{h}}_v)\right)$, where ${\mathbf{c}}_v$ is denoted as ${\mathbf{c}}_v={\sum }_{u^{\prime }}{\alpha }_{v{u}^{\prime }}{\mathit{h}}_{u^{\prime }}$ and ${\alpha }_{v{u}^{\prime }}$ is computed as follows:

  $$\begin{array}{ccc}\hfill {\alpha }_{v{u}^{\prime }}& =& \frac{exp(\mathrm{score}({\mathit{h}}_v^{k-1},{\mathit{h}}_{u^{\prime }}^{k-1}))}{{\sum }_{u^{\prime }}exp(\mathrm{score}({\mathit{h}}_v^{k-1},{\mathit{h}}_{u^{\prime }}^{k-1}))},\hfill \end{array}$$

  (13)

  where the score is calculated as:

  $$\begin{array}{cc}\hfill \mathrm{score}({\mathit{h}}_v^{k-1},{\mathit{h}}_{u^{\prime }}^{k-1})& =V^T\left(tanh(W_1{\mathit{h}}_v^{k-1}+W_2{\mathit{h}}_{u^{\prime }}^{k-1})\right),\hfill \end{array}$$

  where V, $W_1$, $W_2$ and $w_c$ are trainable weights. ${\mathbf{c}}_v$ and ${\mathit{a}}_v$ are the context vector and attention coefficient, respectively. This is inspired by [22].
• Modification of $AG{G}_R$:
  In P-GNN, ${\mathit{M}}_v\left[j\right]$, associated with anchor set j, is averaged across the anchor sets to generate ${\mathit{h}}_v$. We choose to differentiate nodes based on their conditional centrality in stationary status $\left(C{C}_v\left(\omega \right){|}_{\omega \to \infty }\right)$, as higher conditional eigenvector centrality indicates a higher influence of a given node v and its surrounding nodes compared to the ones with lower eigenvector centrality. Additionally, corresponding informative anchor sets contain at least 2-hop neighbor(s) of node v, wherein the ones with higher $C{C}_v\left(\omega \right)$ deserve to have higher weights for aggregation. As such, we substitute $AG{G}_R$ by the weighted mean of ${\mathit{M}}_v\left[j\right]$, where the weights are proportionate to stationary $C{C}_v\left(\omega \right)$, as follows:

  $$\begin{array}{c}\hfill {\mathit{h}}_v^{\left(k\right)}=\frac{1}{J}\sum _j{r}_j{\mathit{M}}_v\left[j\right],\forall j\in [1,J]\end{array}$$

  (14)

  where $r_j$ is calculated as:

  $$\begin{array}{ccc}\hfill r_j& =& \frac{{\sum }_{u^{\prime }\in R_j}C{C}_{{u^{\prime }}_j}\left(\omega \right)}{{\sum }_{j^{\prime }=1}^J{\sum }_{u\in R_{j^{\prime }}}C{C}_{u_{j^{\prime }}}\left(\omega \right)}{|}_{\omega \to \infty },\hfill \\ & & d(v,u^{\prime })\le 2,\hfill \end{array}$$

  (15)

  where $\mathrm{j}$ is the number of anchor sets. Calculation of the conditional node centrality in stationary condition is implemented separately from the main algorithm, and the result has been used in $AG{G}_R$ to aggregate the computed messages across the anchor sets.
  Anchor sets are selected randomly based on Bourgain’s theorem; they come in different sizes that distribute exponentially. $\mathit{M}$ is computed for node v from anchor set $R_j$, if $R_j$ hits v, which means $d(v,R_j)\le 1$. Large anchor sets have a higher probability of hitting v, but are less informative of the positional information of the node, as v hits at least one of many nodes in the anchor sets. On the other hand, small anchor sets have a lower chance of hitting v; however, they provide positional information with high certainty [15]. In terms of message communication, each node in a P-GNN layer shares a message with $\mathcal{O}(nlog{}^{2}n)$ number of anchor sets for a graph with n nodes. Suppose that we have a temporal graph with $n^{\prime }=\left(\right|\mathcal{V}|\times T)$ and average anchor set size of m; the total communicated messages will be $\mathcal{O}(m{n}^{\prime }log{}^2{n}^{\prime })$. This is because each node in TMP-GNN only communicates with nodes that are up to two-hop distant within each anchor set $R_j$, so that m is limited ($1\le m<<n^{\prime }$).

At the end of the last P-GNN layer, the node embedding ${\mathit{z}}_v$ is calculated after applying a nonlinear $\sigma$ to an inner product of a weight vector $\mathit{w}$ and ${\mathit{M}}_v$ as in ${\mathit{z}}_v=\sigma \left(\mathit{w}{\mathit{M}}_v\right)$. The output embedding of for all nodes at all layers t can be reformed into a supraembedding matrix $\mathbb{Z}$ as follows:

$$\mathbb{Z}=\left[\begin{array}{ccccc}{z}_{v_1}^1& z_{v_1}^2& z_{v_1}^3& \cdots & z_{v_1}^T\\ z_{v_2}^1& z_{v_2}^2& z_{v_2}^3& \cdots & z_{v_2}^T\\ z_{v_3}^1& z_{v_3}^2& z_{v_3}^3& \cdots & z_{v_3}^T\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ z_{v_{\left|\mathcal{V}\right|}}^1& z_{v_{\left|\mathcal{V}\right|}}^2& z_{v_{\left|\mathcal{V}\right|}}^3& \cdots & z_{v_{\left|\mathcal{V}\right|}}^T\end{array}\right],$$

(16)

where ${\mathit{z}}_v^t$ is the embedding vector of node v at time t. In some scenarios, edge embedding is needed, rather than node embedding. We estimate $e_{(}uv)$ embedding by averaging the embedding of the ending nodes as:

$${\mathit{z}}_{uv}=\frac{1}{2}({\mathit{z}}_u+{\mathit{z}}_v)$$

(17)

We utilize the result of TMP-GNN embedding for missing data estimation. Before proceeding to our proposed architecture, we briefly review the other components of our framework for missing data.

2.4. Bi-Directional Recurrent Neural Network (Bi-GRU)

We use bi-directional GRU as a component in our proposed architecture to estimate missing points in the temporal multilayer graph $\mathcal{G}$. The Bi-GRU [23] (Figure 3) captures the temporal correlation within a time layer in both forward and backward directions. The input to Bi-GRU is a triplet array ${\mathit{Y}}_e$, ${\mathit{M}}_e$ and ${\Delta }_e$ that is produced from the graph feature set, $\mathbb{X}=\{x_{11}^1,x_{11}^2,\dots ,x_{uv}^t,\dots ,x_{uv}^T,\dots ,x_{u^{\prime }{v}^{\prime }}^T\}$, where $x_{uv}^t$ is the feature associated with $e_{uv}$ at time layer t. We sometimes refer to $x_{uv}$ as $x_e$ interchangeably, where e could be any edge in the graph. We also assume that the graph nodes are constant. However, the number of edges and their corresponding features may change from one layer to another. We set $x_{uv}^t$ to 0, if $e_{uv}$ does not exist at time t (${\rho }^{(e_{uv}\times t)}=0$). We also choose to randomly remove missing points and set it to 0 in ${\mathit{Y}}_e$. ${\mathit{M}}_e$ is called the mask array containing 0 s and 1 s, which indicates the coordinates of missing and observed points, respectively, corresponding to ${\mathit{Y}}_e$. Additionally, each element in ${\Delta }_e$ illustrates the time difference between the current and the last layer at which the measurement is recorded. ${\Delta }_e$ is defined to handle the different sampling rates associated with data heterogeneity from different sources [24,25]. Each edge contains $\mathrm{D}$ streams of features. We use $x_{e_d}$ to represent stream d of the feature associated with e, where e is an edge connecting a pair of nodes in the graph $\mathcal{G}$.

Our goal is to find the best estimate ${\widehat{x}}_{e_d}^t$ with minimum RMSE for a particular missing point through solving the following optimization problem and finding the function f as defined below.

$$\begin{array}{ccc}& & \underset{\mathrm{f}}{min}{\mathbf{E}}_{\mathcal{F}}\left[\sum _{t=1}^T\sum _{d=1}^{\mathrm{D}}(1-m_{e_d}^t)\mathcal{L}({\mathcal{X}}_{e_d}^t,{\widehat{\mathcal{X}}}_{e_d}^t)\right],\hfill \end{array}$$

(18)

where the loss function is defined as $\mathcal{L}({\mathcal{X}}_{e_d}^t,{\widehat{\mathcal{X}}}_{e_d}^t)=({\{x_{e_d}^t-f_{e_d}^t(\mathcal{X},\mathcal{T}))}^2,\mathcal{X}=\{x_{e_d}^t,\forall e,\forall d,\forall t\},$$\mathcal{T}=\{1,2,\dots ,T\}$. $\mathbf{E}$; f and $\mathcal{F}$ indicate expected values for the desired estimator function and unknown probability distribution that the records are sampled from, respectively. ${\mathcal{X}}_{e_d}^t$ shows stream d of all edges’ features at time t, and $m_{e_d}$ represents one element of the mask array ${\mathit{M}}_e$ defined above.

2.5. E-TMP-GNN: Extended TMP-GNN for Missing Data Estimation

Figure 4 and Figure 5 illustrate our proposed hybrid architecture of TMP-GNN and Bi-GRU. We propose two architectures: E-TMP-GNN I and E-TMP-GNN II. In E-TMP-GNN I, the input is supracentrality matrix $\mathbb{C}\left(\omega \right)$, which encodes the intra-layer and inter-layer edge weights of multilayer graph $\mathcal{G}$. The TMP-GNN’s output $\mathbb{Z}$ encodes the topological structure and positional and feature information across the nodes within a layer, and among the same nodes between consecutive layers. Thereafter, the resulting d-dimensional embedding $\mathbb{Z}$ is used to calculate edge embedding as in (16), which is then passed to the input of Bi-GRU to merge with existing edge features and form the below array.

$$\begin{array}{c}\hfill {\mathbb{X}}_{\mathrm{TMP}-\mathrm{PGNN}}=\{x_{e_{d_i}}^{\left(t\right)},{\mathit{z}}_{e_{d_i}}^{\left(t\right)}\forall i\in {\mathbb{N}}_{\mathrm{D}},\forall e\in {\mathcal{E}}^{\left(t\right)},{\rho }^{(e\times t)}=1\}\end{array}$$

(19)

The Bi-GRU is an RNN variant wherein the output of the previous layer is a part of the current layer’s input. This characteristic allows the information to propagate step by step [23] and enables the Bi-GRU to capture long-term temporal correlations of edges in two directions, which is advantageous for missing data imputation. We use one layer of Bi-GRU and form a triplet array of ${\mathit{M}}_e$, ${\mathit{Y}}_e$ and ${\Delta }_e$, as explained earlier in this section. An initial estimation, ${\tilde{x}}_{e_d}^t$, is made by applying a nonlinear activation function to the last hidden states of Bi-GRU, $h_{hybri{d}_{e_d}}^{\left(t\right)}$, as below:

$${\tilde{x}}_{e_d}^t=\sigma (U_d[{\overrightarrow{h}}_{hybri{d}_{e_d}}^{\left(t\right)};{\overleftarrow{h}}_{hybri{d}_{e_d}}^{\left(t\right)}]+g_d),$$

(20)

where the two arrows indicate the forward and backward direction of Bi-GRU, respectively. $U_d$ and $g_d$ are weight and bias vectors. In order to capture the correlation among sequences of $x_{e_{d_i}}{|}_{i=1}^{\mathrm{D}}$, a fully connected layer with a hidden dimension equal to the number of feature streams $\mathrm{D}$ is placed afterward. The output of this layer ${\widehat{x}}_{e_d}^t$ is the final estimate of $x_{e_d}^t$.

$${\widehat{\mathit{x}}}^{\left(t\right)}=\sigma (\mathit{w}{\mathit{h}}^{\left(t\right)}+{\beta }_0),$$

(21)

where ${\mathit{h}}^{\left(t\right)}$ is defined as follows:

$${\mathit{h}}^{\left(t\right)}=\varphi (U{\mathit{x}}^{\left(t\right)}+V{\mathit{Y}}^{\left(t\right)}+{\gamma }_0),$$

(22)

where ${\mathit{Y}}^{\left(t\right)}$ is denoted as $[{\mathit{x}}^{\left(t\right)},{\mathit{M}}^{\left(t\right)}]$.

The Bi-GRU and FC components apply interpolation and imputation, respectively, as they capture the dynamics within a time layer $\left(\{{\mathit{x}}_{e_d}^t,\forall t\in {\mathbb{N}}_T\}\right)$, and across the streams of features $(\{{\{x,{\mathit{z}}^{T\ast }\}}_{e_d}^{\left(t\right)},d=\{1,2,\cdots ,D\}\})$ simultaneously. The combination of Bi-GRU and FC, named Multi-directional RNN, is inspired by [24] (see Figure 6).

Our second proposed architecture is illustrated in Figure 4. As can be seen, the two pipelines have distinct inputs. The input to TMP-GNN is the temporal multilayer network represented as suppracentrality matrix $\mathbb{C}\left(\omega \right)$, and we feed Bi-GRU through $\mathbb{X}$. The last $h_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}$ output from TMP-GNN is then passed through a fully connected layer to reduce the representation dimension as in the following equation:

$${\mathit{h}\prime }_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}=\sigma ({\mathit{w}}_p{\mathit{h}}_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}+{\mathit{b}}_p),$$

(23)

where ${\mathit{w}}_p$ and ${\mathit{b}}_p$ are weight and bias vectors, respectively. The output is then merged by ${\mathit{h}}_{\mathrm{Bi}-{\mathrm{GRU}}_e}^{\left(t\right)}=[{\overrightarrow{\mathit{h}}}_{\mathrm{Bi}-{\mathrm{GRU}}_e}^{\left(t\right)};{\overleftarrow{\mathit{h}}}_{\mathrm{Bi}-{\mathrm{GRU}}_e}^{\left(t\right)}]$ through softmax attention, ${\mathit{h}}_{{\mathrm{pgb}}_e}$ is generated as follows:

$$\begin{array}{ccc}\hfill {\mathit{w}}_p^{\prime },{\mathit{w}}_b^{\prime }& =& \mathrm{Softmax}({\mathit{w}}_p,{\mathit{w}}_b),\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\mathit{h}}_{{\mathrm{pgb}}_e}^{\left(t\right)}& =& {\mathit{w}}_p^{\prime }{\mathit{h}\prime }_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}+{\mathit{w}}_b^{\prime }{\mathit{h}}_{\mathrm{Bi}-{\mathrm{GRU}}_e}^{\left(t\right)}\hfill \end{array}$$

(25)

and an initial estimate of the missing value is made as:

$${\tilde{x}}_{e_d}^{\left(t\right)}=\sigma (Q^d{\mathit{h}}_{{\mathrm{pgb}}_{e_d}}+{\alpha }_0^d)$$

(26)

A fully connected neural network is used at the end to finalize the prediction of missing values as follows:

$${\widehat{x}}_{e_d}^{\left(t\right)}=\sigma ({\mathit{w}}_1{\mathit{h}}_{{\mathrm{pgb}}_e}^d+{\gamma }_1)$$

(27)

The E-TMP-GNN I aims at extracting additional features out of the embedding that is yielded from TMP-GNN, and use it to further enrich the edge feature sets. We demonstrate in Section 3 that the added streams provide a significant improvement in MAE compared to the baseline. TMP-GNN II, however, aims at reducing the number of feature streams by implementing a softmax layer that obtains ${\mathit{h}}_{\mathrm{Bi}-{\mathrm{GRU}}_e}^{\left(t\right)}$ and ${\mathit{h}}_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}$ as inputs, since ${\mathit{h}}_{\mathrm{TMP}-{\mathrm{GNN}}_e}^{\left(t\right)}$ is expected to contain valuable information about the temporal correlation among the same nodes of consecutive time layers.

3. Performance Evaluation

We divide this section to three parts. First, we review our real-world datasets and their associated characteristics in depth. Then, we discuss two potential inputs to our proposed TMP-GNN pipeline and their impact on node embedding performance. Thereafter, we evaluate the performance of E-TMP-GNN I and E-TMP-GNN II on missing data estimation. We implemented E-TMP-GNN I and E-TMP-GNN II in python using Networkx (https://github.com/networkx, accessed on 9 April 2022), PyTorch (https://github.com/pytorch, accessed on 9 April 2022) and Tensorflow (https://github.com/tensorflow/tensorflo, accessed on 9 April 2022). We built in and added various components to the P-GNN implementation [15].

3.1. Datasets

We run our experiment on four datasets: TomTom, COVID-19 Mobility (Mobility) [26], PhD Exchange (PhD), and the Seattle Inductive Loop Detector (Loop Detector).

$\mathbf{TomTom}\mathbf{Dataset}$: This dataset [27] contains space mean speed and average travel time during peak congestion for hundreds of approximately one-km-long road segments across the Greater Toronto Area (GTA). The measurements are collected on a one-minute basis during congestion, where the average space mean speed is less than $70\%$ of the vehicular speed during free-flow conditions. We select a three-hour interval from 4:30 p.m. to 7:30 p.m. on Thursday, 8 September 2016, which is expected to have a high proportion of collected records. During this interval, we select road segments with a significant number of mutual timestamps, where measurements for a sufficient number of the segments are conducted. From the point of view of the multilayer network, we select 10 layers of the GTA road network of 1867 nodes and 985 edges, where each node and edge is the start/ending point of a road segment and a road segment, respectively. We define each edge feature as the ratio of $\frac{\mathrm{space}\mathrm{mean}\mathrm{speed}}{\mathrm{free}-\mathrm{flow}\mathrm{speed}}$ and $\frac{\mathrm{free}-\mathrm{flow}\mathrm{travel}\mathrm{time}}{\mathrm{travel}\mathrm{time}}$ that is in range of [0, 1]. To reduce the correlation between the consecutive records and better assess the performance of pair-wise node classification in TMP-GNN and missing data estimation, the records are downsampled by a rate of $\frac{8}{1}$ ($\delta =\frac{1}{8}$). The edge weights are also normalized so that the features are comparable during the learning process. Moreover, the weighted adjacency matrix is created using the highly accurate geolocations of the road segments in the GTA. Once the node embedding is calculated through TMP-GNN, the edge embedding is estimated by averaging the embedding of the corresponding nodes that the edge is connected to as in (17). The resulting embedding will then be used in missing data estimation.

$\mathbf{PhD}\mathbf{Dataset}$: This case of study utilizes a network that represents the exchange of PhD graduates between universities in the fields of mathematical science. The features we use in our experiments include the graduation year of the students who obtained a doctoral degree, his/her official academic advisor(s) and the degree-granting university. The aim of former studies was to estimate the flow of doctorate exchanges between universities, individuals who graduate from university a and then are hired at the university b. We consider the years 1950–2010, which result in $T=61$ time layers, and there exists a set of $N=231$ connected universities during this time span in the united states. In order to build the graph, directed intra-layer edges are created to represent a doctoral degree obtained by a graduate from university a at year t, who is later hired as a professor at university b, and at least advised one student there. In our case, the edge weight indicates the number of graduate doctorates from university a in year t who was later hired at university b as a faculty member. The edge direction opposes the flow of PhD graduates moving between universities $(b\to a)$ [16,28].

$\mathbf{Loop}\mathbf{Detector}\mathbf{Dataset}$: Our third case of study is the data collected by inductive loop detectors on freeways in the Seattle area [3,29]. The dataset contains temporal sequences of the time mean speed of the freeway system. The speed information at a desired reference point (milepost) is averaged from multiple loop detectors on the main lanes of the same direction at the specific reference point. Measurements are recorded every 5 minutes. The adjacency matrix calculated in [29] is used to convert the data into an undirected graph where an edge indicates that a pair of reference points are connected, without specifying the direction of connectivity. We have downsampled the dataset by a rate of $\frac{2}{1}$ ($\delta =\frac{1}{2}$) and selected 11 time layers of the graph.

$\mathbf{Mobility}\mathbf{Dataset}$: Our last case of study is the human mobility flow dataset in the U.S. during the COVID-19 epidemic [26]. To build the graph, a directed intra-layer edge is constructed to represent the flow of people between pairs of states. Each edge indicates the ratio of visits to the population count among the states. We focused on January 1st–22nd 2020, which led to 22 time layers.

Table 1. Characteristics of the study datasets.

Characteristics	TomTom	Loop Detector	PhD	Mobility
$\left|\mathcal{V}\right|$	1867	323	231	72
$\mathrm{static}:\left|\mathcal{E}\right|,\mathrm{dynamic}$ *: ${\sum }_{t=1}^T{\mathcal{E}}^{\left(t\right)}$	985	1001	10,365 *	2692
$\left|\mathcal{V}\right|\times T$	18,670	3553	14,091	1584
Largest Connected Component (LCC)	50	323	13,847	1144
No. of Isolated Nodes	18,620	3230	244	440

* indicates the graph with a dynamic number of edges over time.

demonstrates the characteristics associated with the graphs built from the 4 selected datasets.

Among all, TomTom and Mobility have the highest and the lowest number of nodes, respectively. From an edge density perspective, we have ordered the graph from the most connected to the least connected as follows: Mobility, Loop Detector, PhD and TomTom. From a topological point of view, TomTom and Loop Detector are static graphs. That is, ${\mathcal{E}}^{\left(t\right)}={\mathcal{E}}^{\left(t^{\prime }\right)},t\ne t^{\prime }$ and the change is in the value of $x_{uv}^{\left(t\right)}/w_v^{\left(t\right)}$ only. However, the PhD and Mobility graphs change throughout the time layers in terms of $|{\mathcal{V}}^{\left(t\right)}|,x_{uv}^{\left(t\right)}$ and $w_v^{\left(t\right)}$. Among the two, PhD has the highest dynamic rate, and turns into the sparsest graph in multiple time layers. Figure 7 demonstrates $|{\mathcal{E}}^{\left(t\right)}|$ per layer t for these graphs.

The third row in the table shows the number of nodes in the temporal multilayered graph, explained in Section 2.2, which is equal to $\left(\right|\mathcal{V}|\times T)$. TomTom and PhD have the largest numbers of nodes in the temporal graph given the number of available time layers. The fourth row indicates the size of the largest connected component (LCC) in https://www.overleaf.com/project/6251d321650dd619f849d51c (accessed on 9 April 2022) for the dynamic graph. The parameter is considered as an important topological invariant of the graph and can be computed using a variant of depth-first search for directed and undirected graphs [30]. Here, graphs can be ordered according to ratio of $\frac{\mathrm{LCC}}{\left|\mathcal{V}\right|}$ in descending order as follows: Loop Detector, PhD, Mobility, and TomTom. The entire Loop detector undirected graph is considered as one connected component; TomTom as a directed graph has multiple weakly connected components, with each component as small as $2\%$ of the total $\left|\mathcal{V}\right|$. We calculate the number of isolated nodes as the nodes that are not connected to the largest component by this equation: $\mathrm{no}.\mathrm{of}\mathrm{isolated}\mathrm{nodes}=N\times T-\mathrm{LCC}$, which could be another indicator of nodes with lower node degrees and, thus, sparser graphs. TomTom has the largest number of isolated nodes, as expected.

3.2. Potential Inputs to TMP-GNN: Single Graph vs. Multiple Graphs

There are several ways in which we can feed P-GNN:

• Supracentrality Graph $\mathbb{C}\left(\omega \right)$: In this case, we represent the centrality of a temporal multilayer graph via a supracentrality matrix that is analogous to a single graph with $N\times T$ number of nodes. The main advantage of this approach is to illustrate the graph with all weighted intra-layer and inter-layer edges. However, the matrix size might be relatively large in the presence of a large number of T. In the next section, we demonstrate that this approach outperforms the one with a multiple-graph input.
• Multiple Graphs: In this approach, a sequence of centrality/adjacency matrices associated with each time layer, $\{{\mathit{A}}^{\left(1\right)},{\mathit{A}}^{\left(2\right)},\cdots ,{\mathit{A}}^{\left(T\right)}\}$, is fed into P-GNN. That is, the multilayer graph is treated as individual instances of a single graph without consideration of inter-layer coupling.
• The input to P-GNN can be a single graph with the entry of its corresponding adjacency matrix as follows:

  $$a_{uv}=\left\{\begin{array}{cc}\hfill [x_{uv}^{\left(1\right)},x_{uv}^{\left(2\right)},\cdots ,x_{uv}^{\left(T\right)}]& :\mathrm{if}u,v \mathrm{are} \mathrm{connected},\hfill \\ 0\hfill & :\mathrm{otherwise}.\hfill \end{array}\right.$$

  (28)
  In this approach, we treat the multilayer graph as a single graph where each element of the adjacency matrix is a time series of features associated with an $e_{uv}$. In order to capture the temporal characteristics of a v/$e_{uv}$, a recurrent neural network is required before applying function F in Figure 2, which leads to high computational complexity. We have not implemented this method.

$\mathbf{Training}$: There are two different ways to train the time-varying graph to compute node embedding. In the case of a single graph, graph nodes are split into training, evaluation, and testing, with $80\%$, $10\%$, and $10\%$ respectively. The classification accuracy on the test set is recorded once the best performance on the evaluation is reached. As for the multiple graphs, there are two techniques to train the graph: first, $80\%$ of the graphs can be used for training and the remaining for evaluation and the test set. Second, $80\%$ of the nodes of all available graphs are used for training, and $10\%$ of the remaining nodes in all time layers of graphs are used for testing, with an equal number for the validation set. Although both techniques can be utilized for time-invariant multiple graphs, e.g., Protein [15,31], we exploit the second choice for our multilayer graph task (Figure 8). It enables us to capture change trajectories throughout the time layers, whereas the first approach splits the graphs without considering the time ordering, and also the behavioral changes of any of the nodes in $20\%$ of the graphs are not taken into account. To the best of our knowledge, this is the first time that such an experimental setting is used for node embedding in multiple graphs.

Inductive vs. Transductive Learning: There are two different settings to learn the node embedding: inductive learning and transductive learning. In the latter setting, the graph is trained with fixed node ordering. In this technique, the model needs to be retrained when some additional nodes are added, the ordering has changed or a new graph is given. Moreover, the node attributes are augmented as one-hot identifiers that could restrict the generalization ability of the model [32]. The authors in [15] assess the transductive learning performance of P-GNN on a link prediction task, in terms of whether a pair of nodes are connected. On the other hand, in inductive learning, the learned positional information can be transferred to an unseen graph. In particular, instead of augmenting the one-hot identifiers of node attributes, an order-invariant scalar could be assigned to nodes whenever node attributes are not available. We choose to train TMP-GNN in an inductive learning setting and utilize the pair-wise node classification task, as in [15], to demonstrate the performance of TMP-GNN embedding.

As discussed in (2), our goal is to compute the embedding for time-varying graphs with potential changes in edge existence ${\rho }^{(\mathcal{E}\times t)}$ or intra-layer edge weights represented in ${\mathit{A}}^{\left(t\right)}$. Therefore, pair-wise node classification would be a more challenging task than link prediction, since the node class is determined according to the status of its surrounding edges and/or the status of other isomorphic nodes that are dynamic in such graphs. When the learning task requires node positional information too, only using structure-aware embedding is not sufficient.

Figure 9 illustrates the $\mathrm{ROC}\mathrm{AUC}$ of TMP-GNN for pair-wise node classification given the second and first type of input, the sequences of individual graphs with each graph associated with one layer and a single supragraph, respectively (Section 3.2). All datasets demonstrate significantly better performance in the presence of a single-graph input versus a multi-graph input. It is highlighted that a single graph has better represented the chronological property of a multilayer temporal graph. The difference in performance improvement achieved by TMP-GNN is significantly larger in the case of a single supragraph for the same reason. In the case of TomTom, Loop Detector and PhD, the other four GNNs perform relatively the same. As for Mobility, all five GNN-based approaches perform well, with relatively negligible differences in classification accuracy. This could be due to the lower number of nodes plus lower dynamics in the number of edges per time layer $|{\mathcal{E}}^{\left(t\right)}|$, as shown in Figure 2, and stronger connectivity within time layers in comparison to other case studies that are confined to the expressive power of TMP-GNN, being less distinct from others. TMP-GNN demonstrates the largest difference in AUC improvement for TomTom from multi-graph to a single supragraph. As noted in Section 3.1, it has the highest number of nodes, with over 5 times more than PhD, the second largest graph, in conjunction with the lowest connectivity degree, and the largest number of isolated nodes, which further underlines the discriminating power of TMP-GNN. It is worth mentioning that the TomTom temporal graph, with a relatively larger number of nodes, requires higher communication message passing in the TMP-GNN layer, as described in Section 2.3, and thus higher computational effort compared to the other three datasets.

For three out of four case studies in multi-graph implementation (Figure 9a), the best performer is TMP-GNN; the second most competent is SAGE. The reason could be the similarity in message computation function F, where the feature information of two nodes $({\mathit{h}}_v,{\mathit{h}}_u)$ or the attended message from u towards node v are concatenated as in (11).

Figure 10 reveals the impact of the availability of the edge feature on $\mathrm{ROC}\mathrm{AUC}$ for all five GNN-based models for the Loop Detector dataset. From a topological point of view, Loop Detector is an undirected static graph. When edge features are not available, all GNNs perform relatively the same. When taking the time-varying edge attributes $x_{uv}^{\left(t\right)}$ into account, it turns into a dynamic graph wherein the node class is also dependent on time-varying attributes, and TMP-GNN demonstrates its superiority over the other widely used structure-aware embedding. This is where the power of TMP-GNN is highlighted for time-varying graphs.

Table 2. Performance Comparison of E-TMP-GNN I, E-TMP-GNN II vs. counterpart methods.

Estimate	TomTom	Loop Detector	PhD	Mobility
${\widehat{\eta }}_{\mathrm{Multi}-\mathrm{graph}}^{\mathrm{ROC}\mathrm{AUC}}(\%)$	7–42	0–22	18–35	6–20
${\widehat{\eta }}_{\mathrm{Sin}\mathrm{gle}\mathrm{supragraph}}^{\mathrm{ROC}\mathrm{AUC}}(\%)$	28–41	20–25	30–37	6–58
${\widehat{\eta }}_{\mathrm{E}-\mathrm{TMP}-\mathrm{GNN}\mathrm{I}}^{\mathrm{MAE}}(\%)$	59–69	0–12	17–27	0
${\widehat{\eta }}_{\mathrm{E}-\mathrm{TMP}-\mathrm{GNN}\mathrm{II}}^{\mathrm{MAE}}(\%)$	94–96	29–96	0	0–15

shows the percentage of $\mathrm{ROC}\mathrm{AUC}$ improvement achieved by TMP-PGNN for multi-graph and single-supragraph implementations, as well as the decrease in MAE achieved by E-TMP-GNN I and II compared to the baselines for all four datasets. We define ${\widehat{\eta }}_{\mathrm{y}}^{\mathrm{metric}}$ in the following equation as a measure to calculate the improvement:

$$\begin{array}{ccc}\hfill & & {\widehat{\eta }}_{\mathrm{y}}^{\mathrm{metric}}=\frac{\left|{\eta }_{\mathrm{y}}^{\mathrm{metric}}-{\eta }_{\mathrm{y}’}^{\mathrm{metric}}\right|}{{\eta }_{\mathrm{y}’}^{\mathrm{metric}}}\times 100,\hfill \end{array}$$

(29)

where $\eta$ notes the performance, metric indicates $\mathrm{ROC}\mathrm{AUC}$ or MAE according to the underlying experiment, and $|\ast |$ is the absolute value. y represents TMP-GNN, E-TMP-GNN I, or E-TMP-GNN II, and y’ is the compared baseline method. It is noted that 0 indicates zero or no improvement.

In Figure 11, we evaluate the performance of our proposed architectures, TMP-GNN I and TMP-GNN II, on the estimation of missing data, and compare it to an RNN variant architecture named M-RNN (Figure 6) and the second best classifier demonstrated in Figure 9b, which is CGN, GAT, GIN, and SAGE for the TomTom, PhD, Loop, and Mobility datasets, respectively. We use Mean Absolute Error (MAE) as a measure of estimation accuracy and choose to remove missing values on a temporal dimension. The points are chosen according to uniform distribution in conjunction with a determined missing threshold $\tau$.

All case studies experience a higher MAE as $\tau$ increases. Moreover, either TMP-GNN I or TMP-GNN II outperform both M-RNN and the counterpart GNN for each case. This is due to the high expressive power of additional features, learning through TMP-GNN representation. The underlying time-varying features proved to be significantly more efficient for missing data estimation than the ones derived from the second-best node classifier in the previous experiment, since it takes position information and inter-layer coupling into account simultaneously; the latter makes the architecture more powerful than the ones utilizing the RNN variant for capturing long-term temporal dependencies. TomTom has the highest MAE compared to others, followed by PhD, Loop and Mobility, respectively. This follows the same order as their corresponding multi-layer graphs when ranked according to their degree of connectivity; the sparser the graph is, the larger the estimation error would be. In addition, the percentage of improvement is larger for TomTom and PhD compared to the other two, which means that the proposed method is more effective for graphs with weaker connectivity.

For three out of four datasets (TomTom, Loop, Mobility), E-TMP-GNN II outperforms E-TMP-GNN I. In the case of the PhD dataset, E-TMP-GNN I is the best performer. The nature of being dynamic is more highlighted in the case of PhD and Mobility, since there is a significant change in the number of edges per time point (Figure 9), in addition to the variation in edge weights existing in the two other datasets. In TMP-GNN I, the time-aware node embedding generated by TMP-GNN will play a more effective role once added to the existing features of the dataset so that it is fed into the interpolation component implemented by Bi-GRU, whereas E-TMP-GNN II reduces the dimension of TMP-GNN’s time-aware embedding, and then combines it with Bi-GRU’s generated representation in an attentive way. In summary, dimension reduction and the softmax function might reduce the effectiveness of the time-aware embedding, later affecting the level of MAE decrease in missing data estimation. The more time-varying the dataset is, the more efficient time-aware representation is needed to better estimate missing data points.

4. Conclusions

We propose TMP-GNN, a generalized version of P-GNN node embedding, for temporal multilayer graphs. It takes short-term temporal correlations, as well as feature and positional information, into account. We also incorporate the notion of eigenvector-based centrality to distinguish nodes with a higher influence on their neighbors. Then, we extend the design to E-TMP-GNN I and II to tackle the missing data challenge in dynamic networks with multiple streams of node/edge features. Conducted experiments on four real-world datasets with diverse characteristics demonstrate significant improvements in node-level classification and missing data estimation tasks. It would be interesting to consider formulating the missing data problem directly into an edge representation task for future work.