GSA-Forecaster: Forecasting
Graph-Based Time-Dependent Data with
Graph Sequence Attention

The first challenge is to model spatial dependency. Graph-based time-dependent data is a set of dependent time series. Spatial dependency refers to the dependency between the different time series. Two time series show spatial dependency when they often change accordingly to each other. To illustrate this, consider the hourly taxi demand around New York Penn Station and Grand Central Terminal, two major train stations in Manhattan, New York City. As Fig. 6 shows, when the taxi demand near New York Penn Station increases (decreases), the demand near Grand Central Terminal will generally increase (decrease) with a 0.82 Pearson correlation coefficient. Therefore, the taxi demand at these two locations exhibits spatial dependency. Spatial dependency can be short-range (between adjacent locations) and long-range (between distant locations), and it is not necessarily related to the physical distance between locations. For example, though the physical distance between New York Penn Station and Grand Central Terminal is larger than the distance between either of them to the Empire State Building, their taxi demands correlate more (i.e., exhibit long-range spatial dependency) than either of their taxi demands correlates with the demand near the Empire State Building (0.33 and 0.48 correlation coefficients, respectively).

Leveraging spatial dependency can lead to better forecasting graph-based time-dependent data, as demonstrated in prior work [9]. For example, if we know there is spatial dependency between time series A and B, we can encode the data in time series A by considering the data in time series B, alleviating noise or fluctuations caused by random factors; we can also leverage their spatial dependency to refine the predictions of future data in time series A by considering the predictions of future data in time series B, and vice versa.

Following practice in [9], we use a graph to represent spatial dependency, where a node of the graph indexes a time series, and an edge corresponds to spatial dependency between two time series. Such graph is not provided as prior knowledge in many applications (e.g., in predicting the taxi demand in Manhattan). Models need to account for this graph and leverage it for forecasting.

The second challenge relates to capturing temporal dependency. Temporal dependency refers to relations between data at different time instants. Data at two time instants show temporal dependency if they are similar to each other. Capturing temporal dependency allows a model to leverage similar historical data/trend to better forecast future data.

We find that graph-based time-dependent data can show both short-range and long-range temporal dependency. For example, Fig. 7 shows that hourly taxi demand near New York Penn Station manifests both daily dependency (short-range; 24 hours between similar data) and weekly dependency (long-range; 168 hours between similar data). Models must be able to capture both short-range and long-range temporal dependency for accurate forecasting.

The third challenge is to account for auxiliary information that can augment existing data used for forecasting. Some auxiliary information may impact the evolution of graph-based time-dependent data. Accounting for such information helps improve the accuracy of forecasting. We illustrate this with the taxi demand around the Javits Center, a major convention center in New York City. Fig. 8 shows the impact of weather (temperature) on the average hourly taxi demand between January 1, 2009 and June 30, 2016, for different hours of the day. We can see that when there is notable taxi demand (e.g., $>$ 20), a decrease in temperature leads to an increase in the taxi demand, which may be because colder weather encourages people to choose more convenient ways of travel (e.g., taxi). Major exogenous factors like weather can cause taxi demand changes (but these exogenous factors are not included in the taxi demand dataset). Models need to consider them as helpful auxiliary information for better predicting taxi demand.

To encode spatial dependency, there are two key problems to be solved: (1) determining the graph that represents spatial dependency when the graph is implicit or not provided; (2) leveraging the graph for prediction.

Determining graph. The statistical model of Gaussian Markov random fields [31] is adopted to learn the spatial dependency graph from data. This way, we can infer spatial dependency based on the variation of the data at each node. This is usually better than estimating spatial dependency based on ad-hoc metrics or heuristics. For instance, in the example shown in Fig. 6, using a physical distance metric to estimate spatial dependency (i.e., assuming that nodes that are physically adjacent to each other are dependent) may mistakenly consider the taxi demand at New York Penn Station and the Empire State Building as spatially dependent because they are physically adjacent, and may mistakenly neglect the spatial dependency between the taxi demand at New York Penn Station and Grand Central Terminal just because of the longer distance between them. For this example, by exploiting the statistical correlation between time series, the theory of Gaussian Markov random fields successfully discovers the spatial dependency between the taxi demand at New York Penn Station and Grand Central Terminal [9].

Gaussian Markov random fields model a graph signal $\mathbf{x}_{t}$ that consists of the data $x_{t}^{i}$ at each node $v_{i}$ at a certain time $t$, i.e., $\mathbf{x}_{t}=\left[{x_{t}^{1}},\cdots,{x_{t}^{N}}\right]^{\top}$, as a random vector drawn from a multivariate Gaussian distribution. The probability density function of $\mathbf{x}_{t}$ is:

where $\bm{\mu}$ and $Q$ are the mean and precision matrix of the distribution.

The precision matrix $Q$ characterizes the relations between data at different nodes. It can be estimated by solving graphical lasso [33], a sparse-penalized maximum likelihood estimator. After estimating the precision matrix $Q$, Forecaster uses it to calculate the conditional correlation between the data of every two nodes [34, 9]. A conditional correlation with large absolute value reveals the spatial dependency between two nodes and hence an edge is introduced to connect them.

Leveraging graph. Forecaster embeds graph structure into the architecture of Transformer, enabling the modified Transformer to encode spatial dependency. To do this, Forecaster proposes sparse linear layers (neural layers that are sparse and represent graph structure) and uses them as substitutes for the fully-connected layers within Transformer (while keeping the nonlinear activation functions and other components). Fig. 9 shows an example of a graph with three nodes $v_{1}$, $v_{2}$, and $v_{3}$ and two sparse linear layers (${l}^{\textrm{th}}$ and ${l+1}^{\textrm{th}}$ layers) designed based on this graph.⁸⁸8Different sparse linear layers can be designed based on this graph, e.g., by changing the number of neurons assigned to different nodes. This is a user choice. To design a sparse linear layer, Forecaster first assigns neurons to the nodes of the graph for encoding their data. In Fig. 9b, as an example, one neuron of the $l^{\textrm{th}}$ layer and two neurons of the ${l+1}^{\textrm{th}}$ layer are assigned to each node in Fig. 9a. Then, Forecaster connects the neurons assigned to the same node or graph-dependent nodes, enabling the learnable weight on each connection to encode the impact of spatial dependency. In Fig. 9b, Forecaster connects the neurons assigned to nodes $v_{1}$ and $v_{2}$ (and connects the neurons assigned to nodes $v_{2}$ and $v_{3}$) as these nodes are strongly dependent. However, as nodes $v_{1}$ and $v_{3}$ are not dependent, their neurons are not connected. In this way, the encoding for the data at a node is only impacted by its own encoding and the encoding of its dependent nodes, which captures the spatial dependency between nodes.

Sparse linear layers are a type of graph convolutional network similar to GCN [35] and TAGCN [36, 37]. All of these networks enrich the encodings of a node by the encodings of its adjacent nodes. The only major difference is, instead of treating adjacent nodes equally, sparse linear layers learn different weights for different adjacent nodes, considering that a node may depend with different strengths on its different adjacent nodes.

Forecaster leverages the standard attention mechanism, which has been employed in Transformer and other alike architectures for NLP tasks [25, 26, 27, 28, 29, 30], to capture temporal dependency within graph-based time-dependent data.⁹⁹9We give a simplified review of how Forecaster leverages the standard attention mechanism. More details are in [9].

Forecaster first encodes the graph signal $\mathbf{x}_{t}$ at each time instant $t$ with encoding $\mathbf{e}_{t}\in\mathbb{R}^{d_{\mathrm{model}}}$, where ${d_{\mathrm{model}}}$ represents the dimension of the encoding. To predict the encoding $\mathbf{e}_{t+1}$ of the next graph signal $\mathbf{x}_{t+1}$, Forecaster uses the encoding of the current graph signal as its initial estimate $\widetilde{\mathbf{e}}_{t+1}$:

Then, Forecaster employs the standard attention mechanism to find historical graph signals that are similar to the initial estimate and uses them to develop a refined estimate. Intuitively, developing an estimate based on these multiple similar graph signals alleviates the impact of noise of each individual graph signal on the estimate,¹⁰¹⁰10The noise refers to the impact of random factors on graph signals, e.g., a random taxi pickup. as suggested by the central limit theorem [38]. This improves the initial estimate that is based solely on the current graph signal.

Specifically, Forecaster compares the initial estimate $\widetilde{\mathbf{e}}_{t+1}$ with the encoding of each historical graph signal $\mathbf{e}_{t+i}$ ($i=-T+1,\ldots,0$; $T$ is the length of the history), and obtains their similarity scores $s_{t+1,\,t+i}^{(h)}$ from multiple perspectives (a perspective is also termed as a head [25, 9]; $h=1,\ldots,H$, $H$ is the number of heads):

where $s_{t+1,\,t+i}^{(h)}\in\mathbb{R}$ is the similarity score between $\widetilde{\mathbf{e}}_{t+1}$ and $\mathbf{e}_{i}$ under head $h$; $W_{Q}^{(h)},W_{K}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$ are learnable parameter matrices for head $h$ that are implemented as a sparse linear layer, respectively; $\left\langle\cdot,\>\cdot\right\rangle$ is the inner product between two vectors, measuring their similarity.

After this step, Forecaster passes the similarity scores through a softmax layer [39] and obtains the attention scores ${\alpha}_{t+1,\,t+i}^{(h)}$ between the initial estimate $\widetilde{\mathbf{e}}_{t+1}$ and the encoding of each historical graph signal $\mathbf{e}_{t+i}$:

where ${\alpha}_{t+1,\,t+i}^{(h)}\in\left(0,1\right)$, $\sum_{k=-T+1}^{0}{\alpha_{t+1,\,t+k}^{(h)}}=1$.

The attention scores characterize the amount of attention that Forecaster pays to each historical graph signal. Forecaster uses the attention scores as weights for historical graph signals to predict the next graph signal:

The resulting $\widehat{\mathbf{e}}_{t+1}$ is a refined estimate for the encoding of the next graph signal; $W_{O}\in\mathbb{{R}}^{d_{\mathrm{model}}\times d_{\mathrm{model}}}$ is a learnable parameter matrix that combines the update $\Delta\widetilde{\mathbf{e}}_{t+1}^{(h)}$ from each head $h$; $W_{V}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$ is a learnable parameter matrix for calculating the updates; $W_{O}$ and $W_{V}$ are implemented as a sparse linear layer, respectively; $\cdot\|\cdot$ represents a concatenation of two vectors.

After getting $\widehat{\mathbf{e}}_{t+1}$, Forecaster decodes it and obtains the predicted next graph signal $\widehat{\mathbf{x}}_{t+1}$. This way, Forecaster completes a single-step forecasting (i.e., forecasting the next graph signal). Forecaster repeats the above process (by advancing the current timestamp $t$ to future timestamps) for multi-step forecasting.

The standard attention mechanism directly relates historical graph signals with predicted graph signals, capturing temporal dependency better than RNN or CNN-based approaches. However, it is still not sufficient to fully capture temporal dependency. To illustrate this, we revisit the example shown in Fig. 4 and perform the standard attention mechanism (described in Section 3.2) using this example. For simplicity, we assume that standard attention has a single head (i.e., $H=1$) and parameter matrices $W_{Q}$ and $W_{K}$ are equal to the identity matrix (i.e., $W_{Q}=W_{K}=I$); and Forecaster encodes the graph signal $x_{t}$ at each time instant with encoding $\mathbf{e}_{t}=\left[e_{t}^{1},\cdots,e_{t}^{n}\right]^{\top}$, such that

where $n=8$, $m=0$, $e_{t}^{i}=1\,\textrm{or}\,-1$.

Fig. 10 shows the resulting attention scores when Forecaster predicts the next graph signal (data point $r^{\prime}$). We can see that Forecaster pays attention to the historical graph signals (data points $a$, $b$, $c$, …, $r$) that are similar to the current graph signal (data point $r$), rather than those (data points $c^{\prime}$ and $g^{\prime}$) that are similar to the next one. The reason behind this is: (1) Forecaster uses the current graph signal as an initial estimate for the next one, introducing error in the initial estimate; (2) the standard attention mechanism compares only the inaccurate initial estimate with the historical graph signals, causing the initial erroneous estimate to induce Forecaster to pay attention to the wrong data. As a result, Forecaster fails to capture the temporal dependency between the data, leading to erroneous prediction.

GSA-Forecaster performs multi-step forecasting: predicting $T^{\prime}$ future graph signals based on $T$ historical graph signals and $T+T^{\prime}$ historical/future auxiliary information. The task can be formalized as learning a function $f\left(\cdot\right)$ that performs the following mapping:

where $\mathbf{x}_{t}$ is the graph signal at time $t$, $\mathbf{x}_{t}=\left[{x_{t}^{1}},\cdots,{x_{t}^{N}}\right]^{\top}\in\mathbb{{R}}^{N}$, with $x_{t}^{i}$ the data at node $i$ at time $t$; $\mathbf{{a}}_{t}\in\mathbb{{R}}^{A}$ is the auxiliary information at time $t$, $A$ the dimension of the auxiliary information.¹¹¹¹11For simplicity, we assume that different nodes share the same auxiliary information. It is easy to generalize our model to other cases. More specifically, GSA-Forecaster performs the multi-step forecasting in an autoregressive manner: after predicting the graph signals $\{\mathbf{x}_{t+i}\mid i=1,\ldots,k-1\}$, GSA-Forecaster uses them together with other known information to predict the next graph signal $\mathbf{x}_{t+k}$. GSA-Forecaster iterates this process (i.e., $k=1,\ldots,T^{\prime}$) until predicting all the $T^{\prime}$ future graph signals.¹²¹²12The same computations that appear in multiple iterations are performed only once. This way, GSA-Forecaster transforms the multi-step forecasting task into a series of single-step tasks of forecasting $\mathbf{x}_{t+k}$. In the rest of Section 4, we focus on how GSA-Forecaster solves this single-step forecasting task. We use the following conventions: $t$ as current timestamp, $\widehat{\mathbf{x}}_{t+k}$ as our prediction result for $\mathbf{x}_{t+k}$, and, by abuse of notation, $\{\mathbf{x}_{t+i}\mid i=1,\ldots,k-1\}$ as our previously predicted graph signals.¹³¹³13In the rest of this section, the actual values of the graph signals $\{\mathbf{x}_{t+i}\mid i=1,\;\ldots,\;k-1\}$ are not involved. Given this, and to simplify our equations, we use the same notations to refer to the predicted values of these graph signals.

As Fig. 11(a) shows, the architecture of GSA-Forecaster consists of embeddings, an encoder, and an decoder. The encoder and the decoder employ two types of graph sequence attention with similar mathematical formalism but different goals (i.e., GSA (Filtering) for noise filtering and GSA (predicting) for data prediction, respectively) as shown in Fig. 11(b). To capture spatial dependency, GSA-Forecaster uses sparse linear layers as shown in Section 3.1 that reflect the graph structure throughout its architecture.

The goal of GSA (predicting) in each decoder layer is to predict the encoding of the next graph signal. To do this, it looks for the previous graph signals that show temporal dependency with the next graph signal, and it uses these previous graph signals to predict the next graph signal. It takes the following entities as input as shown in Fig. 11:

• •

  entry C: previous data sequence
  $\left\{\mathbf{e}_{t+i}^{\prime}\mid i=-T+1,\ldots,k-1\right\}$

• •

  entry D: estimated encoding of the next graph signal ¹⁴¹⁴14For the first decoder layer, the entry D input $\mathbf{e}_{t+k}^{\prime}$ is equal to the initial estimate $\mathbf{e}_{t+k-1}$. For each of the other decoder layers, the entry D input $\mathbf{e}_{t+k}^{\prime}$ is equal to the refined estimate made and output by the previous decoder layer and may be different from the initial estimate.
  $\mathbf{e}_{t+k}^{\prime}$

• •

  entry E: auxiliary information
  $\left\{\mathbf{e}_{t+i}^{a}\mid i=-T+1,\ldots,k\right\}$

It outputs a refined estimate for the encoding of the next graph signal $\widehat{\mathbf{e}}_{t+k}$. The procedures for calculating $\widehat{\mathbf{e}}_{t+k}$ are described below.

First, GSA (predicting) estimates the similarity between the next graph signal $\mathbf{x}_{t+k}$ and each previous graph signal $\mathbf{x}_{t+i}$ by calculating their similarity scores. With such scores, we can find the previous graph signals that the next graph signal truly depends on. The similarity scores are calculated based on the following equation:

where $s_{t+k,\,t+i}^{(h)}$ is the similarity score between the next graph signal $\mathbf{x}_{t+k}$ and the previous graph signal $\mathbf{x}_{t+i}$ under head $h$,¹⁵¹⁵15Same as prior work [9, 25], we also compare graph signals under multiple perspectives, where a perspective is termed as a head. $i=-T+M,\ldots,k$, $h=1,\ldots,H$, $M$ is the temporal neighborhood size, $H$ is the number of heads; $w,{w_{A}},{w_{P}}\in\mathbb{R}^{+}$ are learnable parameters; $\left|\cdot\right|$ calculates the norm of a vector; $\left\langle\cdot,\>\cdot\right\rangle$ is the inner product between two vectors, measuring their similarity.

The similarity score $s_{t+k,\,t+i}^{(h)}$ consists of three terms. The first term approximates the similarity between the next graph signal $\mathbf{x}_{t+k}$ and the previous graph signal $\mathbf{x}_{t+i}$ by comparing their temporal neighborhoods. Here, the temporal neighborhood (TN) at time $t+j$ ($j=-T+M,\ldots,k$) contains the encodings of graph signals $\mathbf{e}^{\prime}_{t+j-M+1,}$ $\ldots_{,}$ $\mathbf{e}^{\prime}_{t+j}$ (see Fig. 12 for illustration). As the first term is based on comparing multiple graph signals within the two temporal neighborhoods, the inaccuracy or noise of individual graph signal (e.g., inaccuracy in the estimate $\mathbf{e}_{t+k}^{\prime}$) has less impact on it. Therefore, it more accurately and robustly characterizes the temporal dependency between graph signals. In the first term, $W_{Q}^{(h)},W_{K}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$ are learnable parameter matrices for comparing graph signals. They are implemented as sparse linear layers.

The second term of $s_{t+k,\,t+i}^{(h)}$ approximates the similarity by comparing auxiliary information $\mathbf{e}_{t+k}^{a}$ and $\mathbf{e}_{t+i}^{a}$. Intuitively, similar auxiliary information has similar impact on the graph signals and thus should be accounted for. In the second term, $W_{QA}^{(h)},W_{KA}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{aux}}}{H}\times d_{\mathrm{aux}}}$ are learnable parameter matrices implemented as fully-connected layers.

The third term of $s_{t+k,\,t+i}^{(h)}$ approximates the similarity by comparing temporal positional encodings $\mathbf{e}_{k}^{p}$ and $\mathbf{e}_{i}^{p}$. In some applications, graph signals at some relative temporal positions may often show temporal dependency (e.g., the current taxi demand usually depends on the taxi demand a week ago). We use the learnable temporal positional encodings $\{\mathbf{e}_{i}^{p}\mid i=-T+M,\ldots,k\}$ to capture this kind of temporal dependency that often hold at some relative temporal positions. In the third term, $W_{QP}^{(h)},W_{KP}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{pos}}}{H}\times d_{\mathrm{pos}}}$ are learnable parameter matrices implemented as fully-connected layers.

After getting the similarity scores, we pass them through a softmax layer to calculate the attention scores, which characterize the amount of attention that is paid to each previous graph signal for predicting the next graph signal:

where $\alpha_{t+k,\,t+i}^{(h)}\in\left(0,1\right)$ is the attention score for the previous graph signal $\mathbf{x}_{t+i}$, $i=-T+M,\ldots,k-1$.

To illustrate the attention scores, we perform GSA (predicting) on the example shown in Fig. 4. For simplicity, we let $M=10,H=1,w=10,w_{A}=w_{P}=0$, matrices $W_{Q}^{(h)}=W_{K}^{(h)}=I$ (i.e., identity matrix) and reuse the encoding strategy described in Eq. 6. Fig. 12 shows the attention scores. We can see that when predicting the next graph signal (data point $r^{\prime}$), GSA (predicting) pays attention to only the previous graph signals (data points $c^{\prime}$ and $g^{\prime}$) similar to the next graph signal, successfully capturing temporal dependency.

In some cases, especially when strong non-stationarity happens, the next graph signal is not similar to any of the previous graph signals. In these circumstances, we cannot predict the next graph signal based on the previous ones that show temporal dependency with it. Instead, we should focus on the recent graph signals and try to find the pattern of how they evolve over time. We term this pattern as the recent trend. We use RNN to learn the recent trend and predict based on it. To this end, we extend the attention scores to reflect the relative importance between capturing temporal dependency and capturing the recent trend:

where ${\alpha}_{t+k,\,t+i}^{\prime(h)}\in\left(0,1\right)$ is the extended attention score for the previous graph signal $\mathbf{x}_{t+i}$ ($i=-T+M,\ldots,k-1$), which measures the amount of attention paid to $\mathbf{x}_{t+i}$ for capturing temporal dependency; while ${\alpha}_{t+k,\,t+k}^{\prime(h)}\in\left(0,1\right)$ is the extended attention score that measures the amount of attention paid to the recent trend. $\sum_{j=-T+M}^{k}{\alpha}_{t+k,\,t+j}^{\prime(h)}=1$.

The difference between Eqs. 10 and 9 is that the denominator on the right hand side of Eq. 10 includes the term associated with ${s}_{t+k,\,t+k}^{(h)}$ (the similarity score between the next graph signal and itself). When the next graph signal is not similar to any of the previous graph signals, ${s}_{t+k,\,t+k}^{(h)}$ is much larger than other similarity scores, and, as a result, ${\alpha}_{t+k,\,t+k}^{\prime(h)}\approx 1$, ${\alpha}_{t+k,\,t+i}^{\prime(h)}\approx 0$, $i=-T+M,\ldots,k-1$, which means we should focus on the recent trend to predict.

With the extended attention scores, we calculate the refined estimate for the next graph signal $\widehat{\mathbf{e}}_{t+k}$ as:

where $\Delta{\mathbf{e}}_{t+k}^{\prime(h)}$ is the update to the original estimate $\mathbf{e}_{t+k}^{\prime}$ under head $h$.

$\Delta{\mathbf{e}}_{t+k}^{\prime(h)}$ consists of two terms. The first term uses the extended attention score ${\alpha}_{t+k,\,t+i}^{\prime(h)}$ as the weight for the encoding of each previous graph signal $\mathbf{e}_{t+i}^{\prime}$, capturing the impact of temporal dependency ($W_{V}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$ is a learnable parameter matrix implemented as a sparse linear layer).

The second term reflects the impact of the recent trend, which is captured by GRU [41], a popular type of RNN. The GRU takes the encodings of graph signals in the recent temporal neighborhood $\{\mathbf{e}_{t+i}^{\prime}\mid i=k-M+1,\ldots,k-1\}$ as inputs, and it outputs its update to the original estimate of the next graph signal. Specifically, we first reset the GRU to its default state and then input $\mathbf{e}^{\prime}_{t+k-M+1}$, $\mathbf{e}^{\prime}_{t+k-M+2}$, $\cdots$, and $\mathbf{e}^{\prime}_{t+k-1}$ to the GRU one by one. After that, the GRU outputs its update. Note that, for different $k$ ($k=1,\cdots,T^{\prime}$), the inputs to the GRU are different. Therefore, we need to repeat the above procedures (i.e., reset state and input corresponding graph signals) for each $k$.

We concatenate the update $\Delta{\mathbf{e}}_{t+k}^{\prime(h)}$ of all the heads ($\cdot\|\cdot$ represents a concatenation of two vectors) and use a learnable parameter matrix $W_{O}\in\mathbb{{R}}^{d_{\mathrm{model}}\times d_{\mathrm{model}}}$ (implemented as a sparse linear layer) to aggregate them. Finally, the aggregated update is added to the original estimate $\mathbf{e}^{\prime}_{t+k}$, and we obtain the final estimate $\widehat{\mathbf{e}}_{t+k}$.

The goal of GSA (filtering) in each encoder layer is to filter out noise in the encodings of historical graph signals. GSA (filtering) inputs the following entities as shown in Fig. 11:

• •

  entry A: the encodings of historical graph signals
  $\{\mathbf{e}_{t+i}\mid i=-T+1,\ldots,0\}$

• •

  entry B: auxiliary information
  $\left\{\mathbf{e}_{t+i}^{a}\mid i=-T+1,\ldots,0\right\}$

It outputs the filtered encodings of historical graph signals $\{\mathbf{e}_{t+i}^{\prime}\mid i=-T+1,\ldots,0\}$. The procedures of calculating the filtered encodings are described below.

First, GSA (filtering) estimates the similarity between the true values of any two graph signals $\mathbf{x}_{t+i}$ and $\mathbf{x}_{t+j}$ by calculating their similarity scores. Here, true values refer to the values of graph signals as if noise were removed, which are of course unknown. With the similarity scores, we can add graph signals with similar true values together, letting their noise cancel each other to filter out noise, as the central limit theorem suggests [38]. The similarity scores are calculated as:

where $s_{t+i,\,t+j}^{(h)}$ is the similarity score between the true values of graph signals $\mathbf{x}_{t+i}$ and $\mathbf{x}_{t+j,}\,$ $i,\,j=-T+1,\ldots,0$, $h=1,\ldots,H$, $H$ is the number of heads; $l_{1}=-\mathrm{min}(M_{1},\,T-1+\mathrm{min}(i,\,j))$, $l_{2}=\mathrm{min}(M_{2},\,-\mathrm{max}(i,\,j))$, $M_{1}$ and $M_{2}$ are temporal neighborhood parameters. Simiar to the case before, $W_{Q}^{(h)},\,W_{K}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$, $W_{QA}^{(h)},\,W_{KA}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{aux}}}{H}\times d_{\mathrm{aux}}}$, $W_{QP}^{(h)},\,W_{KP}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{pos}}}{H}\times d_{\mathrm{pos}}}$, and $w,w_{A},w_{P}\in\mathbb{R}^{+}$ are learnable parameters, and $\{\mathbf{e}_{i}^{p}\mid i=-T+1,\ldots,0\}$ are learnable temporal positional encodings.

The similarity score $s_{t+i,\,t+j}^{(h)}$ consists of three terms. The first term approximates the similarity between the true values of the graph signals $\mathbf{x}_{t+i}$ and $\mathbf{x}_{t+j}$ by comparing their temporal neighborhoods. Here, the temporal neighborhood at time $t+n$ ($n=-T+1,\ldots,0$) contains the encodings of graph signals $\mathbf{e}_{t+n-\mathrm{min}(n+T-1,\,{M_{1}}),}$ ${\ldots}_{,}$ $\mathbf{e}_{t+n+\mathrm{min}(-n,\,{M_{2}}),}$ where parameters $M_{1}$ and $M_{2}$ limit the size of the temporal neighborhood. As comparing temporal neighborhoods involves comparing multiple graph signals, the first term is less impacted by the noise from individual graph signals, and it thus more accurately and robustly reflects the similarity. The second and third terms capture the impact of auxiliary information and temporal positional information, respectively.

After obtaining the similarity scores, GSA (filtering) passes them through a softmax layer to calculate the attention scores:

where the attention score $\alpha_{t+i,\,t+j}^{(h)}\in\left(0,1\right)$ characterizes the amount of attention paid to the graph signal $\mathbf{x}_{t+j}$ for filtering the graph signal $\mathbf{x}_{t+i,}\,$ $i,\,j=-T+1,\ldots,0$.

Finally, GSA (filtering) uses the attention scores as weights to add graph signals with similar true values together, filtering out their noise. The filtered encoding $\mathbf{e}_{t+i}^{\prime}$ ($i=-T+1,\ldots,0$) is computed as:

where $\Delta{\mathbf{e}}_{t+i}^{(h)}$ is the update to the original encoding $\mathbf{e}_{t+i}$ under head $h$; $W_{O}\in\mathbb{{R}}^{d_{\mathrm{model}}\times d_{\mathrm{model}}}$ and $W_{V}^{(h)}\in\mathbb{{R}}^{\frac{d_{\mathrm{model}}}{H}\times d_{\mathrm{model}}}$ are learnable parameter matrices implemented as sparse linear layers.

The GSA-Forecaster described earlier uses a single scalar to characterize the similarity between two graph signals (see Eqs. 8 and 12). In other words, different nodes share the same attention score. We call this type of attention mechanism graph-level attention, as the whole graph shares the same attention score. The graph-level attention may not work well with datasets with strong non-stationarity. For example, unexpected events may temporarily impact the data associated with a node such that the node does not have the same attention score as other nodes some times. For such datasets, it is ideal that we incorporate a node-level attention mechanism on top of the graph-level attention to deal with non-stationarity by allowing different nodes to have different attention scores.

Fig. 13 illustrates how we extend GSA-Forecaster. Our design consists of two types of GSA-Forecaster — GSA-Forecaster (Graph-Level Attention) and GSA-Forecaster (Node-Level Attention).

GSA-Forecaster (Graph-Level Attention) is described in Sections 4.1– 4.4. It takes a series of graph signals $\{\mathbf{x}_{t+i}\mid i=-T+1,\ldots,k-1\}$ and a series of auxiliary information $\{\mathbf{a}_{t+i}\mid i=-T+1,\ldots,k\}$ as inputs, and predicts the next graph signal $\widehat{\mathbf{x}}_{t+k,\textrm{graph}}$. Note that a graph signal is a vector that stacks the data associated with each node. For example, the predicted next graph signal $\widehat{\mathbf{x}}_{t+k,\textrm{graph}}=\left[{\widehat{x}_{t+k,\textrm{graph}}^{1}},\cdots,{\widehat{x}_{t+k,\textrm{graph}}^{N}}\right]$, where $\widehat{x}_{t+k,\textrm{graph}}^{j}$ is the predicted next data at node $j$.

GSA-Forecaster (Node-Level Attention) is similar to the design described in Sections 4.1– 4.4 with the only major modification that it predicts for each node separately. For predicting the next data associated with node $j$, it takes only the previous data at node $j$ ($\{x_{t+i}^{j}\mid i=-T+1,\ldots,k-1\}$) as inputs, and it outputs its prediction $\widehat{x}_{t+k,\textrm{node}}^{j}$. Its architecture stays the same except that all the sparse linear layers are replaced with dense layers, because only the data associated with a single node passes through its architecture each time. This way, different nodes can have different attention scores when non-stationarity occurs. We reuse the same GSA-Forecaster (Node-Level Attention) for predicting the data of different nodes. As GSA-Forecaster (Node-Level Attention) is mainly designed for coping with non-stationarity, it does not leverage temporal positional encodings that help capture temporal dependency during stationarity (i.e., set $w_{P}=0$ in Eqs. 8 and 12).

After obtaining the predictions ($\widehat{x}_{t+k,\textrm{graph}}^{j}$ and $\widehat{x}_{t+k,\textrm{node}}^{j}$) made by the two different GSA-Forecasters, the Joiner module combines the two predictions and outputs the final prediction $\widehat{x}_{t+k}^{j}$. Though Joiner can use sophisticated strategies to combine the two predictions, here, for simplicity, it uses the following strategy:

where $\gamma$ is a hyperparameter that is set based on the validation set. We repeat this process to predict the next data associated with each node.

The node-level attention is an optional extension that we only apply for datasets with strong non-stationarity. In our evaluation, we apply this extension for the METR-LA and PEMS-BAY datasets.

Similar to prior work [4, 19, 9], our evaluation uses the root-mean-square error (RMSE) [43] and the mean absolute percentage error (MAPE) [43] to measure the accuracy of the forecasting results. RMSE characterizes the absolute error, while MAPE characterizes the relative error. Lower RMSE and MAPE mean better accuracy.

For the NYC Taxi dataset, following the practice in prior work [19, 9], we set a threshold on taxi demand when calculating MAPE: if $x_{t}^{i}<20$ ($x_{t}^{i}$ refers to the hourly taxi demand at location $i$ at time $t$), exclude $x_{t}^{i}$ from the calculation of MAPE.

The key hyperparameters of GSA-Forecaster are in Table I. GSA-Forecaster uses the following loss function:

where a constant weight $\eta$ is used such that RMSE and MAPE have almost equal contributions to the loss function.

We compare GSA-Forecaster against the following state-of-the-art deep learning models:

DCRNN [4] is an RNN- and GCN-based model. It uses RNN to capture temporal dependency and uses graph convolutional network (GCN) to capture spatial dependency.

Graph WaveNet [22] is a CNN- and GCN-based model. It uses dilated CNN to model temporal dependency and uses GCN to model spatial dependency.

AGCRN [23] is an RNN- and GCN-based model. It fits the graph structure of GCN during training.

GMAN [44] is an attention mechanism-based model. It applies attention over spatial and temporal dimensions.

Transformer [25] is an attention mechanism-based model. It uses the standard attention mechanism to capture temporal dependency.

Forecaster [9] is also an attention mechanism-based model. It is built upon Transformer, leveraging its capability of capturing temporal dependency. It embeds the graph structure into its architecture to handle spatial dependency.

Table II compares GSA-Forecaster with the baseline models on the accuracy of forecasting hourly taxi demand in New York City — overall accuracy and accuracy of the predictions in each of the future three hours. We run each model five times and report the mean and standard deviation of their forecasting accuracy. As the table shows, the baseline model Forecaster outperforms the other baselines DCRNN, Graph WaveNet, GMAN, AGCRN, and Transformer. This is because (1) the standard attention mechanism in Forecaster captures temporal dependency better than the RNN in DCRNN and AGCRN as well as the CNN in Graph WaveNet; (2) Forecaster captures spatial dependency (by embedding the graph structure into its architecture), while Transformer fails to do so; (3) Forecaster leverages the theory of Gaussian Markov random fields to learn the dependency graph from data rather than using a pre-defined graph like GMAN does. Our proposed GSA-Forecaster achieves the best forecasting accuracy as it further enhances Forecaster’s capability to capture temporal dependency by using graph sequence attention (vs. the standard attention used by Forecaster). Compared with the best baseline Forecaster, GSA-Forecaster achieves more accurate predictions for each of the future three hours and reduces the overall RMSE and MAPE by 6.7% and 5.8%, respectively. Additionally, we ablate auxiliary information and test its impact on GSA-Forecaster. We find that auxiliary information improves the forecasting accuracy of GSA-Forecaster. However, even if auxiliary information is ablated, GSA-Forecaster still significantly outperforms the baseline models.

To provide a visual example for the performance of the models, Fig. 14 (a) – (e) show their three-hour-ahead forecasting results for the taxi demand around New York Penn Station during November 8 – 14, 2015. We compare the forecasting error of DCRNN, Graph WaveNet, Transformer and Forecaster with that of GSA-Forecaster and show their difference in Fig. 14 (f) – (i) (a positive difference means the baseline model has larger error). We can see that, overall, GSA-Forecaster has better accuracy. In this example, compared with the best baseline Forecaster, GSA-Forecaster has smaller error 63% of the time; on average, GSA-Forecaster has 13% less error.

In Table III, we compare GSA-Forecaster with the baseline models on the accuracy of forecasting traffic speed in Los Angeles (the METR-LA dataset) and the Bay Area (the PEMS-BAY dataset)—accuracy of the predictions in the future 15 minutes (i.e., short term), 30 minutes (i.e., mid-term), and 60 minutes (i.e., long-term). Here, we report their RMSE for forecasting non-zero traffic speeds (i.e., RMSE (non-zero speed)) and all traffic speeds (i.e., RMSE (all speed)). The latter takes into account the accuracy of forecasting zero traffic speeds, which has significant societal implications, e.g., accurately predicting zero traffic speeds helps government agencies and transportation departments to take measures to reduce traffic congestion. We also report MAPE for forecasting non-zero traffic speeds (i.e., MAPE (non-zero speed); it is invalid to calculate MAPE for zero traffic speeds). As the table shows, GSA-Forecaster outperforms DCRNN, Transformer, Forecaster, GMAN, and AGCRN in all the categories. Compared with the best baseline Graph WaveNet, GSA-Forecaster achieves better performance in 13 categories out of the total 18 categories. GSA-Forecaster significantly outperforms Graph WaveNet on the RMSE of forecasting all speeds (i.e., RMSE (all speed)) by up to 25%, while slightly losing in the short-term prediction for non-zero speeds by at most 2%. The advantage of GSA-Forecaster comes from its better capability of handling non-stationarity, for example, when unexpected events cause traffic congestion and zero traffic speed occurs. To provide an overall comparison across different metrics and different forecasting categories, we calculate the average normalized error of a model as follows and average it across different categories (15-, 30-, and 60-minute ahead predictions):