# Joint Forecasting and Interpolation of Graph Signals Using Deep Learning

Gabriela Lewenfus, Wallace Alves Martins, Symeon Chatzinotas, and Björn
Ottersten Ms. Lewenfus is with the Federal University of Rio de Janeiro
(UFRJ), (e-mail: gabriela.lewenfus@gmail.com); Prof. Martins is with UFRJ (on
leave) and with the University of Luxembourg (as Research Associate), (e-mail:
wallace.martins@smt.ufrj.br); Prof. Chatzinotas is with the University of
Luxembourg (SnT), (e-mail: symeon.chatzinotas@uni.lu); and Prof. Ottersten is
with the University of Luxembourg (SnT), (e-mail: bjorn.ottersten@uni.lu)This
work was financed in part by the ERC project AGNOSTIC and by the Brazilian
research agencies CNPq and Faperj.© 2020 IEEE. Personal use of this material
is permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for
advertising or promotional purposes, creating new collective works, for resale
or redistribution to servers or lists, or reuse of any copyrighted component
of this work in other works.

###### Abstract

We tackle the problem of forecasting network-signal snapshots using past
signal measurements acquired by a subset of network nodes. This task can be
seen as a combination of multivariate time-series prediction and graph-signal
interpolation. This is a fundamental problem for many applications wherein
deploying a high granularity network is impractical. Our solution combines
recurrent neural networks with frequency-analysis tools from graph signal
processing, and assumes that data is sufficiently smooth with respect to the
underlying graph. The proposed approach outperforms state-of-the-art deep
learning techniques, especially when only a small fraction of the graph
signals is accessible, considering two distinct real world datasets:
temperatures in the US and speed flow in Seattle. The results also indicate
that our method can handle noisy signals and missing data, making it suitable
to many practical applications.

###### Index Terms:

Multivariate time series, forecasting and interpolation, deep learning,
recurrent neural networks (RNNs), graph signal processing (GSP)

##  I Introduction

Spatiotemporal (ST) prediction is a fundamental abstract problem featuring in
many practical applications, including climate analyses [1], transportation
management [2], neuroscience [3], electricity markets[4], and several
geographical phenomenon analyses [5]. The temperature in a city, for instance,
is influenced by its location, by the season, and even by the hour of the day.
Another example of data with ST dependencies is the traffic state of a road,
since it is influenced by adjacent roads and also by the hour of the day. ST
prediction boils down to _forecasting_ (temporal prediction) and
_interpolation_ (spatial prediction). The former refers to predicting some
physical phenomenon using historical data acquired by a network of spatially-
distributed sensors. The latter refers to predicting the phenomenon with a
higher spatial resolution. In this context, ST data can be seen as a network
signal in which a time series is associated with each network element; the
dynamics (time-domain evolution) of the time series depends on the network
structure (spatial domain), rather than on the isolated network elements only.
The interpolation is useful to generate a denser (virtual) network.

Classical predictive models generally assume independence of data samples and
disregard relevant spatial information [6], [7]. Vector autoregressive (VAR)
[8], a statistical multivariate model, and machine learning (ML) approaches,
such as support vector regression (SVR) [9] and random forest regression [10],
can achieve higher accuracy than classical predictive models; yet, they fail
to fully capture spatial relations. More recently, some progress has been made
by applying neural networks111In this paper the word “network” can refer to a
neural network in the context of deep learning or a physical network that is
represented by a graph. (NNs) to predict ST data [11, 2, 12, 1, 13, 14]. NNs
have the capacity of not only mapping an input data to an output, but also of
learning a useful representation to improve the mapping accuracy [15].
Nonetheless, the fully-connected architecture of these NNs fails to extract
simultaneous spatial and temporal features from data.

In order to learn spatial information from these multivariate time series,
some works have combined convolutional NNs (CNNs) with recurrent NNs (RNNs),
such as long short-term memory (LSTM) [16, 17, 18, 19, 20]. However, CNNs are
restricted to grid-like uniformly structured data, such as images and videos.
To overcome this issue, and inspired by graph222Graphs are mathematical
structures able to represent rather general datasets, including ST data with
irregular domains, as in sensor networks. signal processing (GSP), some works
have developed convolutions on graph-structured data (graph signals) [21, 22,
23, 24, 25], which have been used in combination with either RNN, time
convolution, and/or attention mechanisms to make predictions in a variety of
applications. These works are summarized in TABLE I.

TABLE I: Summary of recent works that use deep learning to predict ST data. First column defines the application. Second and third columns refer to the spatial and temporal techniques employed, respectively. “Conv.” means temporal convolution; AE means auto encoder; RBM means restricted Boltzman machine; “other” encompasses other predictive strategies, such as attention mechanisms Application | GSP | Temporal |  Paper
---|---|---|---
traffic | Yes | RNN |  [25], [26, 27, 28, 29, 30, 31]
Conv. |  [32, 33, 34, 35, 36, 37]
other |  [38], [39]
No | RNN |  [17], [20], [40, 41, 42, 43, 44],
other |  [2] ,[19], [45, 46, 47, 48, 49]
wind | Yes | RNN |  [50]
No | RNN |  [14]
Other |  [51, 52, 53]
meteorological | Yes | AE |  [54]
No | RNN |  [16]
Other |  [55, 5]
body-motion | Yes | RNN |  [56, 57, 58]
neuroscience | Yes | Conv. |  [59]
No | RBM |  [60]
semantic | Yes | RNN |  [61]

GSP theory has been applied to analyze/process many irregularly structured
datasets in several applications [62, 63]. An import task addressed by GSP is
interpolation on graphs, i.e., (spatially) predicting the signals on a subset
of graph nodes based on known signal values from other nodes [64]. In general,
graph interpolation is based on local or global approaches. Local methods,
such as $k$-nearest neighbors ($k$-NNs) [65], compute the unknown signal
values in a set of network nodes using values from their closest neighbors,
being computationally efficient. Global methods, on the other hand,
interpolate the unknown signal values at once and can provide better results
by taking the entire network into account at the expensive of a higher
computational burden [64, 66]. Many GSP-based interpolation techniques have
been proposed [67, 68, 64, 69, 70, 71, 72, 71, 73, 74]. Due to the irregular
structure of graph signals, the interpolation problem may become ill-
conditioned, calling for efficient selection strategies for obtaining optimal
sampling sets [75, 76]. In fact, the problem of interpolating a graph signal
(GS) can also be addressed as semi-supervised classification [77, 78, 79, 80,
81] or regression [82, 83, 84] tasks. More recently, deep learning (DL)
solutions have also been developed [24, 85].

In many applications, affording to work with large networks may be
impractical; for example, placing many electrodes at once in the human cortex
may be unfeasible. Installation and maintenance costs of devices can also
limit the number of sensors deployed in a network [84]. Thus, developing a
predictive model capable of forecasting (temporal prediction) and
interpolating (spatial prediction) time-varying signals defined on graph nodes
can be of great applicability. This problem can be regarded as a semi-
supervised task, since only part of the nodes are available for training.
Other works have addressed this problem: in [86] the graph is extended to
incorporate the time dimension and a kernel-based algorithm is used for
prediction; this approach, therefore, relies on the assumption of smoothness
in the time domain, which is not reasonable for many applications, such as
traffic flow prediction. In [50], the ST wind speed model is evaluated in a
semi-supervised framework in which only part of the nodes are used for
training the model, while interpolation is performed only in test phase.
Therefore, the parameters learned during the training phase do not take into
account the interpolation aspect.

Two straightforward solutions to deal with the problem of forecasting and
interpolating sampled GSs are:333We can regard a GS that will be interpolated
as a sampled version of a GS defined over a denser set of nodes belonging to a
virtual graph. (i) applying a forecasting model to the input GS and then
interpolating the output; or (ii) interpolating the sampled GS and then
feeding it to the forecasting model. These solutions tackle the ST prediction
task separately and may fail to capture the inherent coupling between time and
space domains. In this paper, a graph-based NN architecture is proposed to
handle ST correlations by employing GSP in conjunction with a gated-recurrent
unit (GRU). Thus, we address the inherent nature of ST data by jointly
forecasting and interpolating the underlying network signals. A global
interpolation approach is adopted as it provides accurate results when the
signal is smooth in the GSP sense, whereas an RNN forecasting model is adopted
given its prior success in network prediction. Herein, not only the sampled GS
is input to a predictive model but also its spectral components, which carry
spatial information on the underlying graph. The major contribution of our
proposed model is, therefore, the ability to learn ST features by observing a
few nodes of the entire graph.

Considering the proposed learning framework, we introduce four possible
classes of problems:

  * •

supervised applications, where the labels of all nodes are available for
training but only a fixed subset of graph nodes can be used as input to the
model in the test phase;

  * •

semi-supervised application, wherein only data associated with a subset of
nodes are available for training and computing gradients;

  * •

noise-corrupted application, in which all nodes are available during the
entire process, but additive noise corrupts the network signals;

  * •

missing-value application, where a time-varying fraction of nodes are
available for testing, but all nodes can be used for training.

The proposed approach achieves the best results in most tested scenarios
related to the aforementioned applications, as compared to DL-based
benchmarks.

The paper is organized as follows: Section II presents some fundamental
aspects of GSP, focusing on the sampling theory that will be used to build the
interpolation module of the proposed learning model. Section III describes the
new learning framework. Section IV describes four classes of applications that
can benefit from the proposal. Section V presents the numerical results and
related discussions. Section VI contains the concluding remarks of the paper.

##  II GSP Background

Let $\mathcal{G}\triangleq(\mathcal{V},\mathcal{E},\mathbf{A})$ be a weighted
undirected and connected graph, where
$\mathcal{V}\triangleq\\{v_{1},\ldots,v_{N}\\}$ is the set of $N$ nodes,
$\mathcal{E}$ is the set of edges, and $\mathbf{A}$ is the $N\times N$
adjacency matrix containing edge weights $A_{mn}$. The adjacency matrix can be
a similarity matrix or be built based on prior information, such as nodes’
locations in the physical network. A GS is a real-valued scalar function
$x:\mathcal{V}\to\mathbb{R}$ taking values on the graph nodes, and it will be
represented by the $N$-dimensional vector $\mathbf{x}$ with entries
$\left[\mathbf{x}\right]_{n}=x_{n}=x(v_{n})$.

The diagonal matrix $\mathbf{D}\in\mathbb{R}^{N\times N}$ is the degree matrix
in which $D_{nn}=\sum_{m}A_{nm}$ measures the connectivity degree of each
node.

Most of the graph convolutions in the literature are based on the Laplacian
matrix $\mathbf{L}\triangleq\mathbf{D}-\mathbf{A}$, which is a symmetric semi-
definite matrix, for a symmetric adjacency ${\bf A}$. Let $\mathbf{U}$ be the
matrix of orthonormal eigenvectors of $\mathbf{L}$. The graph Fourier
transform (GFT) of the GS $\mathbf{x}$ is
$\hat{\mathbf{x}}\triangleq\mathbf{U}^{\rm T}\mathbf{x}$ and the graph
frequencies are considered as the eigenvalues of $\mathbf{L}$,
$\lambda_{1},...,\lambda_{N}\geq 0$. The eigenvectors of the adjacency matrix
can also be used as the Fourier basis [87]. For the reader’s convenience,
TABLE II contains the main notations that will be used in this work.

###  II-A GS Sampling

Let $\mathcal{S}\triangleq\\{s_{1},...,s_{M}\\}\subset{\mathcal{V}}$ be a
subset of nodes with $M\leq N$ nodes; the vector of measurements
$\mathbf{x}_{\mathcal{S}}\in\mathbb{R}^{M}$ is given by
$\mathbf{x}_{\mathcal{S}}=\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{x}$, where
the sampling operator

$\displaystyle\left[\boldsymbol{\Psi}_{\mathcal{S}}\right]_{mn}=\left\\{\begin{array}[]{cc}1,&\textrm{if\;}v_{n}=s_{m}\\\
0,&{\rm otherwise}\end{array}\right.$ (3)

selects from $\mathcal{V}$ the nodes in $\mathcal{S}$. The interpolation
operator $\boldsymbol{\Phi}_{\mathcal{S}}$ is an $N\times M$ matrix such that
the recovered signal is
$\tilde{\mathbf{x}}=\boldsymbol{\Phi}_{\mathcal{S}}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{x}$.
If $\tilde{\mathbf{x}}=\mathbf{x}$, the pair of sampling and interpolation
operators $(\boldsymbol{\Phi}_{\mathcal{S}},\boldsymbol{\Psi}_{\mathcal{S}})$
can perfectly recover the signal $\mathbf{x}$ from its sampled version. As the
rank of $\boldsymbol{\Phi}_{\mathcal{S}}\boldsymbol{\Psi}_{\mathcal{S}}$ is
smaller or equal to $M$, this is not possible for all
$\mathbf{x}\in\mathbb{R}^{N}$ when $M<N$. However, perfect reconstruction can
be achieved for a class of bandlimited GS.

The GS $\mathbf{x}_{\rm b}$ is said $\mathcal{F}$-bandlimited if
$\left[\hat{\bf x}_{\rm b}\right]_{n}=0\;\forall n$ such that
$\lambda_{n}\not\in\mathcal{F}\subset\\{\lambda_{1},\ldots,\lambda_{N}\\}$,
that is, the frequency content of $\mathbf{x}_{\rm b}$ is restricted to the
set of frequencies $\mathcal{F}$. Some works also restrict the support of the
frequency content and consider that a GS $\mathbf{x}_{\rm b}$ is
$\omega$-bandlimited if $\left[\hat{\bf x}_{\rm b}\right]_{n}=0\;\forall n$
such that $\lambda_{n}>\omega$[88]. In this paper, a bandlimited signal is a
sparse vector in the GFT domain. The following theorem guarantees the perfect
reconstruction of an $\mathcal{F}$-bandlimited GS for some sampling sets.

###### Theorem 1

([71],[68]) If the sampling operator $\boldsymbol{\Psi}_{\mathcal{S}}$
satisfies

$\displaystyle{\rm
rank}(\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}})=|\mathcal{F}|=K,$
(4)

then $\mathbf{x}_{\rm
b}=\boldsymbol{\Phi}_{\mathcal{S}}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{x}_{\rm
b}$ as long as
$\boldsymbol{\Phi}_{\mathcal{S}}=\mathbf{U}_{:,\mathcal{F}}\boldsymbol{\Sigma}$,
with $\boldsymbol{\Sigma}$ satisfying
$\boldsymbol{\Sigma}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}}=\mathbf{I}_{K}$
and $\mathbf{U}_{:,\mathcal{F}}$ a submatrix of $\mathbf{U}$ with columns
restricted to the indices associated with the frequencies in $\mathcal{F}$.

The condition in (4) is also equivalent to

$\displaystyle{\rm
SV}_{\max}(\mathbf{U}_{\overline{\mathcal{S}},\mathcal{F}})\leq 1,$ (5)

where ${\rm SV}_{\max}(.)$ stands for the largest singular value [89] and
$\overline{\mathcal{S}}=\mathcal{V}/\mathcal{S}$. This means that no
$\mathcal{F}$-bandlimited signal over the graph $\mathcal{G}$ is supported on
$\overline{\mathcal{S}}$.

In order to have
$\boldsymbol{\Sigma}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{V}_{:,\mathcal{F}}=\mathbf{I}_{K}$
we must have $M\geq K$, since ${\rm rank}(\mathbf{U}_{:,\mathcal{F}})=K$. If
$M\geq K$, $\boldsymbol{\Sigma}$ is the pseudo-inverse of
$\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}}$ and the
interpolation operator is

$\displaystyle\boldsymbol{\Phi}=\mathbf{U}_{:,\mathcal{F}}(\mathbf{U}_{:,\mathcal{F}}^{\rm
T}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}})^{-1}\mathbf{U}_{\mathcal{S},\mathcal{F}}^{\rm
T}\,.$ (6)

Since $\mathbf{U}$ is non-singular, there is always at least a subset
$\mathcal{S}$ such that the condition in (4) is satisfied. Nonetheless, for
many choices of $\mathcal{S}$,
$\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}}$ can be full rank
but ill-conditioned, leading to large reconstruction errors, especially in the
presence of noisy measurements or in the case of approximately bandlimited GS.
To overcome this issue, optimal sampling strategies, in the sense of
minimizing reconstruction error, can be employed [71]. Note that
$\boldsymbol{\Phi}$ depends on both $\mathcal{S}$ and $\mathcal{F}$, but this
dependence is omitted in order to simplify the notation.

###  II-B Approximately Bandlimited GS

In practice, most GSs are only approximately bandlimited[79]. A GS is
approximately $(\mathcal{F},\epsilon)$-bandlimited if [89]

$\displaystyle\mathbf{x}=\mathbf{x}_{\rm b}+\boldsymbol{\eta},$ (7)

where $\mathbf{x}_{\rm b}$ is an $\mathcal{F}$-bandlimited GS and
$\boldsymbol{\eta}$ is an $\overline{\mathcal{F}}$-bandlimited GS such that
$\|\boldsymbol{\eta}\|_{2}<\epsilon$. If signal $\mathbf{x}$ is sampled on the
subset $\mathcal{S}$ and recovered by the interpolator in (6), the error
energy of the reconstructed signal is upper bounded by

$\displaystyle\|\tilde{\mathbf{x}}-\mathbf{x}\|_{2}\leq\frac{\|\boldsymbol{\eta}\|_{2}}{{\rm
cos}(\theta_{\mathcal{S},\mathcal{F}})},$ (8)

where $\theta_{\mathcal{S},\mathcal{F}}$ is the maximum angle between the
subspace of signals supported on $\mathcal{S}$ and the subspace of
$\mathcal{F}$-bandlimited GS. It can be shown that ${\rm
cos}(\theta_{\mathcal{S},\mathcal{F}})={\rm
SV}_{\min}(\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}})$;
therefore, in order to minimize the upper bound of the reconstruction error in
(8), the set $\mathcal{S}$ should maximize ${\rm
SV}_{\min}(\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}})$.
Finding the optimal set $\mathcal{S}$ is a combinatorial optimization problem
that can require an exhaustive search in all possible subsets of $\mathcal{V}$
with size $M$. A suboptimal solution can be obtained by the greedy search in
[71, Algorithm 1].

TABLE II: Notations Notation |  Definition
---|---
$\mathcal{G}$ |  graph
$\mathcal{V}$ |  entire set of graph nodes
$\mathcal{S}$ |  subset of graph nodes
$\mathcal{F}$ |  subset of graph spectrum
$\overline{\mathcal{S}}$ |  the complement of the set $\mathcal{S}$
$\mathbf{L}$ |  Laplacian matrix
$\mathbf{U}$ |  matrix of Laplacian eigenvectors
$\mathbf{U}_{:,\mathcal{F}}$ |  submatrix of $\mathbf{U}$ with columns in the set $\mathcal{F}$
$\mathbf{U}_{\mathcal{S},\mathcal{F}}$ |  submatrix of $\mathbf{U}$ with columns in the set $\mathcal{F}$ and rows in $\mathcal{S}$
$\boldsymbol{\Psi}_{\mathcal{S}}$ |  sampling operator $\mathcal{V}\to\mathcal{S}$
$\boldsymbol{\Phi}_{\mathcal{S}}$ |  interpolation operator $\mathcal{S}\to\mathcal{V}$
$\mathbf{x}_{\mathcal{S}}$ |  signal $\mathbf{x}$ restricted to the set $\mathcal{S}$
$\hat{\mathbf{x}}_{\mathcal{F}}$ |  the frequency content of the GS $\mathbf{x}$ restricted to the set $\mathcal{F}$

##  III Joint Forecasting and Interpolation of GSs

In this section, we propose an ST neural network to jointly interpolate graph
nodes and forecast future signal values. More specifically, the task is to
predict the future state $\mathbf{x}^{t+p}$ of a network given the history
$\mathbf{X}_{\mathcal{S}}^{t}=\\{\mathbf{x}_{\mathcal{S}}^{t},...,\mathbf{x}_{\mathcal{S}}^{t-\tau+1}\\}$.444The
GS at timestamp $t$ is denoted by bold lowercase letter, $\mathbf{x}^{t}$,
whereas the history set containing the sampled GSs in previous timestamps is
denoted by bold capital letter, $\mathbf{X}_{\mathcal{S}}^{t}$. Thus, the
input signal is a GS composed by $M$ nodes and the output GS is a network-
signal snapshot composed by $N\geq M$ nodes. For now, to describe the learning
model’s architecture, we shall assume $p=1$.

###  III-A Gated-Recurrent Unit

The proposed learning architecture employs a GRU cell [90] as the basic
building block. The GRU cell is composed by a hidden state $\mathbf{h}^{t}$,
which allows the weights of the GRU to be shared across time, as well as by
two gates $\mathbf{q}^{t}$ and $\mathbf{r}^{t}$, which modulate the flow of
information inside the cell unit. Fig. 1 depicts the architecture of a GRU.
The gates are given by:

$\displaystyle\mathbf{q}^{t}$
$\displaystyle=\sigma(\mathbf{W}_{q}\mathbf{x}^{t}+\mathbf{V}_{q}\mathbf{h}^{t-1}+\mathbf{b}_{q}),$
(9) $\displaystyle\mathbf{r}^{t}$
$\displaystyle=\sigma(\mathbf{W}_{r}\mathbf{x}^{t}+\mathbf{V}_{r}\mathbf{h}^{t-1}+\mathbf{b}_{r}),$
(10)

where
$\\{\mathbf{W}_{q},\mathbf{V}_{q},\mathbf{W}_{r},\mathbf{V}_{r}\\}\subset\mathbb{R}^{M\times
M}$ are matrices whose entries are learnable weights,
$\\{\mathbf{b}_{q},\mathbf{b}_{r}\\}\subset\mathbb{R}^{M}$ are the bias
parameters, and $\sigma(\cdot)$ is the sigmoid function.

The update of the hidden state $\mathbf{h}^{t}$ is a linear combination of the
previous hidden state and the candidate state $\mathbf{c}^{t}$:

$\displaystyle\mathbf{c}^{t}$
$\displaystyle=\sigma(\mathbf{W}_{c}\mathbf{x}^{t}+\mathbf{V}_{c}(\mathbf{h}^{t-1}\odot\mathbf{r}^{t})+\mathbf{b}_{c}),$
(11) $\displaystyle\mathbf{h}^{t}$
$\displaystyle=\mathbf{q}^{t}\odot\mathbf{c}+(1-\mathbf{q}^{t})\odot\mathbf{h}^{t-1},$
(12)

with $\odot$ being the element-wise multiplication. Similar to LSTM [91], the
additive update of the hidden state can handle long-term dependencies by
avoiding a quick vanishing of the errors in back-propagation, and by not
overwriting important features. The GRU structure was chosen to compose the
forecasting module of the proposed method because its performance is usually
on pair with LSTM, but with a lower computational burden [92]; nonetheless, it
could be replaced by an LSTM or any other type of RNN.

$\times$element-wisemultiplication$+$element-
wiseaddition$\sigma$$\sigma$Tanh$\times$+$\times$$\times$$-1$$\mathbf{h}^{t-1}$$\mathbf{x}^{t}$$\mathbf{c}^{t}$$\mathbf{z}^{t}$$\mathbf{r}^{t}$$\mathbf{x}^{t+1}$$\mathbf{h}^{t}$
Figure 1: GRU cell.

###  III-B Forecasting Module

The proposed learning model, named spectral graph GRU (SG-GRU), combines a
standard GRU cell applied to the vertex-domain GSs comprising
$\mathbf{X}_{\mathcal{S}}^{t}$ with a GRU cell applied to the frequency-domain
versions of the latter GSs comprising $\hat{\mathbf{X}}_{\mathcal{F}}^{t}$.
The GRU acting on frequency-domain signals is named here spectral GRU (SGRU),
and has the same structure as the standard GRU, except for the dimension of
weight matrices and bias vectors, which are $K\times K$ and $K$, respectively.
The dimension of the hidden state is therefore $K$.

Assuming that the entire GS $\mathbf{x}$ is
$(\mathcal{F},\epsilon)$-bandlimited, most of the information about it is
expected to be stored in $\hat{\mathbf{x}}_{\mathcal{F}}$. Then, given an
admissible operator $\boldsymbol{\Psi}_{\mathcal{S}}$, one has

$\displaystyle\|\mathbf{x}-\boldsymbol{\Phi}_{\mathcal{S}}\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{x}\|_{2}\leq\frac{\epsilon}{{\rm
SV}_{\min}(\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}})}.$ (13)

The choice of $\mathcal{F}$ will be further discussed in the experiments
described in Section V.

The SGRU module in the proposed learning framework is able to predict the
(possibly time-varying) graph-frequency content of the network signals. This
is key to fully grasping the underlying spatial information embedded in the
graph frequency content. Besides, it is worth pointing out that the proposed
SGRU is slightly different from simply combining the spectral graph
convolution (SGC) in [93] with a GRU: in both cases, the input signal is
previously transformed to the Fourier domain, but in the SGRU a standard GRU,
composed by matrix-vector multiplications, is applied to the transformed
signal, whereas in the latter case, the SGC computes the graph convolution,
which is an element-wise vector multiplication. Thus, the SGRU is able to
better capture the temporal relations among different spectral components.

###  III-C GS Interpolation

The outputs from the GRU, $\mathbf{y}_{\mathcal{S}}^{t+1}$, and from the SGRU
,$\hat{\mathbf{z}}_{\mathcal{F}}^{t+1}$, are interpolated by
$\boldsymbol{\Phi}_{\mathcal{S}}$ and
$\hat{\boldsymbol{\Phi}}_{\mathcal{S}}=\boldsymbol{\Phi}_{\mathcal{S}}\mathbf{U}_{:,\mathcal{F}}$,
respectively. The resulting $N$-dimensional vectors $\mathbf{y}^{t+1}$ and
$\mathbf{z}^{t+1}$ are stacked in a single vector of size $2N$ which is
processed by a fully connected (FC) layer to yield

$\displaystyle\tilde{\mathbf{x}}^{t+1}=\Theta(\mathbf{y}^{t+1},\mathbf{z}^{t+1}),$
(14)

as illustrated in Fig. 2.

$\mathbf{X}_{\mathcal{S}}^{t}$GRUGFTSGRU$\boldsymbol{\Phi}$$\boldsymbol{\hat{\Phi}}$Vec$\boldsymbol{\Theta}$$\mathbf{x}^{t+1}$$\hat{\mathbf{X}}_{\mathcal{S}}^{t}$$\mathbf{y}_{\mathcal{S}}^{t+1}$$\hat{\mathbf{z}}_{\mathcal{F}}^{t+1}$$\mathbf{y}^{t+1}$$\mathbf{z}^{t+1}$
Figure 2: Proposed SG-GRU model. The input GS follows two routes in parallel:
in the upper route, the GRU followed by interpolation is applied to the GS; in
the bottom route, the GS is transformed to frequency-domain before being
processed by the SGRU module and thereafter being interpolated. The outputs of
these two parallel processes are stacked into a single vector, represented by
operation “Vec”, and fed to an FC layer.

###  III-D Loss Function

The loss function employed is the (empirical) mean square error (MSE). In a
supervised scenario, the signal from all nodes are available for training,
thus enabling the use of the entire GS $\mathbf{x}^{t+1}$ as label to compute
the loss function. Given the batch size $T_{\rm b}$, the loss function for the
supervised training $\mathcal{L}_{\rm s}$ is

$\displaystyle\mathcal{L}_{\rm s}=\frac{1}{T_{\rm
b}N}\sum\limits_{t=1}^{T_{\rm
b}}\|\tilde{\mathbf{x}}^{t+1}-\mathbf{x}^{t+1}\|_{2}^{2}.$ (15)

In the semi-supervised task, on the other hand, only the sampled ground-truth
signal $\mathbf{x}^{t+1}_{\mathcal{S}}$ can be accessed. In order to achieve
better predictions on the unknown nodes, we propose to interpolate the sampled
ground-truth signal by $\boldsymbol{\Phi}_{\mathcal{S}}$ before computing the
MSE, yielding

$\displaystyle\mathcal{L}_{\rm ss}=$ $\displaystyle\frac{1}{T_{\rm
b}N}\sum\limits_{t=1}^{T_{\rm
b}}\|\tilde{\mathbf{x}}^{t+1}-\mathcal{I}(\mathbf{x}^{t+1}_{\mathcal{S}})\|_{2}^{2},$
(16)

where

$\displaystyle\displaystyle\left[\mathcal{I}(\mathbf{x_{\mathcal{S}}})\right]_{n}=\left\\{\begin{array}[]{cc}x_{n},&{\rm
if\;}v_{n}\in\mathcal{S}\\\
\left[\boldsymbol{\Phi}_{\mathcal{S}}\mathbf{x}^{t+1}_{\mathcal{S}}\right]_{n},&{\rm
otherwise.}\end{array}\right.$ (19)

###  III-E Computational Complexity

The SG-GRU consists of two GRU cells, refereed as GRU and SGRU, which compute
$6$ matrix-vector multiplications each. The dimensions of the weight matrices
in these recurrent modules applied on the vertex and frequency domains are
$M^{2}$ and $K^{2}$, respectively, where $K$ was set to $\frac{M}{3}$ in this
paper (this choice will be further discussed in Section V). The input of the
SGRU is the sampled GS in the frequency domain, obtained by applying the
truncated GFT, which is a $K\times M$ matrix. This transform can be pre-
computed, avoiding the matrix vector multiplication during the loop
recurrence. In this case the input of the network becomes a signal with
dimension $M+K$. The output of the GRU and the SGRU are, thereafter,
interpolated by $N\times M$ and $N\times K$ matrices, respectively, which are
pre-computed before running the model. Finally, an FC layer is applied to the
interpolated signals, costing $2N^{2}$ flops. Note that the truncated GFT, the
interpolations, and the FC layers are out of the recurrence loop and do not
increase the computational cost if a larger sequence length $\tau$ is used.
Thus, the computational cost per iteration of the SG-GRU is

$\displaystyle KM+6\tau(M^{2}+K^{2})+N(K+M)+2N^{2}\;\text{[flops]}.$ (20)

##  IV Applications

The proposed learning architecture in Fig. 2 can handle both supervised and
semi-supervised scenarios. In the supervised case, measurements from the $N$
network nodes are available in the training step but not necessarily for
testing. This supervised scenario covers many different applications; a case
in point is a weather station network wherein the temperature sensors are
working during a period of time, but then, suddenly, some of them are shut
down due to malfunctioning or maintenance cost reduction. In the semi-
supervised case, on the other hand, only part of the nodes appear in the
training set and can, therefore, be used to compute gradients. Again, the
semi-supervised scenario also covers many practical applications; for
instance, when a sensor network is deployed with a limited number of nodes to
reduce the related costs, but a finer spatial resolution is desirable, which
can be obtained by a virtual denser sensor network.

Considering these two basic scenarios, we can conceive four specific types of
applications:

###  IV-A Supervised Application

Input GS is composed by $M\leq N$ nodes but labels of all $N$ nodes are used
to compute the loss function in (15). As mentioned before, this learning model
can be applied to situations in which all the $N$ sensors are temporarily
activated and, afterwards, $N-M$ sensors are turned off.

###  IV-B Semi-supervised Application

Both input GS and labels are composed by $M<N$. Thus, only the $M$ in-sample
are available to train the model using the loss function in (16). In this
application, it is desired to predict the state of a static network with $N$
nodes, considering that only $M<N$ sensors are deployed.

###  IV-C Noise-corrupted Application

Input GS is composed by all the $N$ nodes with signals corrupted by
uncorrelated additive noise, and the labels are the entire ground-truth GS.
This application allows working with the proposed learning model when the
sensors’ measurements are not accurate. In this case, only the denoising
capacity of the proposed model is evaluated, hence no sampling is performed
over the input data.

###  IV-D Missing-value Application

Input GS is composed by all the $N$ nodes but, at each time instant, a
fraction of the $N$ values measured by the sensor network are randomly chosen
to be replaced by NaN (not a number). It is worth highlighting that this
application is different from the (pure) supervised application in Section
IV-A. In the supervised scenario, the set of known nodes, $\mathcal{S}$, is
fixed across time, whereas the application of missing values considers
different sets of known signal values at each time instant $t$. In other
words, we have a supervised scenario with a time-dependent sampling set
$\mathcal{S}^{t}$. The labels are the entire ground-truth GS. This setup
evaluates the performance of the proposed SG-GRU when some of the sensors’
measurements are missing, which could be due to transmission failures in a
wireless network.

##  V Numerical Experiments

In this section, we assess the performance of the proposed SG-GRU scheme in
two real datasets. The simulation scenarios are instances of the four
applications described in Section IV.

###  V-A Dataset Description

The proposed learning model was evaluated on two distinct multivariate time-
series datasets: temperatures provided by the Global Surface Summary of the
Day Dataset (GSOD), which can be accessed at [94], and the Seattle Inductive
Loop Detector Dataset (SeattleLoop) [40].

####  V-A1 Global Surface Summary of the Day Dataset

The GSOD dataset consists in daily temperature measurements in oC from $2007$
to $2013$, totalling $2,557$ snapshots, in $430$ weather stations distributed
in the United States.555Weather stations in the Alaska and in Hawai were not
considered. The source provides more weather stations but only $430$ worked
fully from $2007$ until $2013$. These stations are spatially represented by a
$10$-nearest-neighbor graph with nonzero edge weights given by [95]:

$\displaystyle
A_{nm}\\!=\\!\frac{e^{-(d_{nm}^{2}+h_{nm}^{2})}}{\sqrt{\sum_{j\in\mathcal{N}_{n}}e^{-(d_{nj}^{2}+h_{nj}^{2})}}\sqrt{\sum_{j\in\mathcal{N}_{m}}e^{-(d_{mj}^{2}+h_{mj}^{2})}}},$
(21)

in which $\mathcal{N}_{n}$ is the set of neighboring nodes connected to the
node indexed by $n$, whereas $d_{nm}$ and $h_{nm}$ are, respectively, the
geodesic distance and the altitude difference between weather stations indexed
by $n$ and $m$. The adjacency matrix is symmetric and the diagonal elements
are set to zero.

####  V-A2 Seattle Inductive Loop Detector Dataset

The SeattleLoop dataset contains traffic-state data collected from inductive
loop detectors deployed on four connected freeways in the Greater Seattle
area. The $323$ sensor stations measure the average speed, in miles$/$hour,
during the entire year of $2015$ in a $5$-minute interval, providing $105,120$
timesteps. This dataset is thus much larger than GSOD. The graph adjacency
matrix provided by the source [40] is binary and the GS snapshots are barely
bandlimited with respect to the graph built on this adjacency matrix. To build
a network model in which the SeattleLoop time series is
$(\mathcal{F},\epsilon)$-bandlimited with a reasonably small $\epsilon$, the
nonzero entries of the binary adjacency matrix were replaced by the radial-
basis function

$\displaystyle A_{nm}={\rm
e}^{-\frac{\|\mathbf{x}_{n}-\mathbf{x}_{m}\|^{2}}{10}},$ (22)

where $\mathbf{x}_{n}$ and $\mathbf{x}_{m}$ are time series, containing $1000$
time-steps, corresponding to nodes $v_{n}$ and $v_{m}$, respectively.

###  V-B Choice of Frequency Set $\mathcal{F}$

The larger the set $\mathcal{F}$ the more information about the input signal
is considered in the model. However, the interpolation using (6) is admissible
only if $|\mathcal{F}|=K\leq M$ [71]. Moreover, if $K$ increases, the smaller
singular value of $\mathbf{U}_{\mathcal{S},\mathcal{F}}$ tends to decrease,
leading to an unstable interpolation. Since the GSs considered in this paper
are approximately bandlimited, using $K$ close to $M$ accumulates error during
the training of the network. Based on validation loss, $K$ was set to
$\frac{M}{3}$.

When all nodes are available for training, that is, in the applications
described in Sections IV-A, IV-C, and IV-D, $\mathcal{F}$ is chosen as the $K$
Laplacian eigenvalues corresponding to the dominant frequency components (the
ones with highest energy) of signals measured at the first $100$ days. In the
semi-supervised approach, on the other hand, the spectral content of the
entire GS is unknown. Since the GSs considered in this paper are usually
smooth, in the sense that most of their frequency content is supported on the
indices associated with the smaller Laplacian eigenvalues, the set
$\mathcal{F}$ was chosen as the $K$ smallest eigenvalues $\lambda_{n}$ in this
scenario. The set $\mathcal{F}$ used in the application described in Section
IV-B is, therefore, slightly different from the set $\mathcal{F}$ used in the
scenarios related to the applications of Sections IV-A, IV-C, and IV-D.

###  V-C Competing Learning Techniques

Recently many DL-based models were shown to outperform classical methods in
the task of predicting ST data. Nonetheless, to the best of our knowledge,
only [50] addresses the problem of predicting ST data by training a learning
model with $M<N$ nodes, with the aim of reducing the training time duration.
Therefore, the performance of our proposed method is here compared with DL-
based models from the literature that do not actually handle sampled input
GSs. Thus, we adapted the DL-based models from the literature by combining
them with an interpolation strategy, such as $k$-NN and the GSP-based
interpolator $\boldsymbol{\Phi}_{\mathcal{S}}$. In this context, the
interpolation can be performed either: (i) before running the forecasting
technique, so that the input of the competing DL-based model will be the
entire GS; or (ii) after running the forecasting technique, so that the input
of the competing DL-based model will be a sampled GS, thus requiring fewer
learnable parameters.

We use as benchmark some LSTM-based NNs, which were shown to perform well in
the strict forecasting task (i.e., time-domain prediction) on the SeattleLoop
dataset in comparison with other baseline methods, such as ARIMA and SVR [31].
In addition, we also consider the ST graph convolution network (STGCN)
proposed in [32] as benchmark. In summary, the competing techniques (adapted
to deal with sampled GSs) are:

  1. (i)

LSTM: simple LSTM cell;

  2. (ii)

C1D-LSTM: a $1$D convolutional layer followed by an LSTM cell;

  3. (iii)

SGC-LSTM: the SGC from [93] followed by an LSTM;

  4. (iv)

TGC-LSTM: a traffic graph convolution based on the adjacency matrix combined
with LSTM [31];666Code from
https:$//$github.com$/$zhiyongc$/$Graph$\\_$Convolutional$\\_$LSTM .

  5. (v)

STGCN: a combination of the graph convolution from [22] with a gated-temporal
convolution [32]; Hyperparameters were set as in [32] since they lead to
smaller MSE in the validation set (filter sizes were evaluated from the set
$\\{16,32,64\\}$).777Code from
https:$//$github.com$/$VeritasYin$/$Project$\\_$Orion.

As mentioned before, the above competing techniques do not tackle joint
forecasting and interpolation tasks. Thus, they were combined with an
interpolation technique. The output of methods (i)-(iv) were interpolated by
$\boldsymbol{\Phi}_{\mathcal{S}}$, whereas a $1$-hop neighborhood
interpolation was applied before the method (v), that is, each unknown value
$x^{t}_{n}$ was set as

$\displaystyle[x^{t}_{n}]_{\rm
unknown}=\frac{1}{|\mathcal{N}_{n}|}\sum_{m\in\mathcal{N}_{n}}x_{m}^{t}\,.$
(23)

Unlike LSTM-based methods, the interpolation in (23) provided better results
when combined with STGCN to handle the sampled input GS. The TGC-LSTM was only
applied to the SeattleLoop dataset since it uses a free-flow reachability
matrix, being specifically designed for traffic networks.

TABLE III summarizes the competing learning techniques along with the
corresponding interpolation methods. “GSP interpolation” stands for
interpolation by $\boldsymbol{\Phi}_{\mathcal{S}}$ and “$1$-hop interpolation”
stands for interpolation by averaging records on the $1$-hop neighborhood.
Also the order followed by each procedure is indicated: “interpolation first”
stands for interpolating the data before running the learning model, whereas
“model first” refers to running the model and then interpolating its output.

TABLE III: Summary of competing techniques Model | Interpolation | Order of procedures
---|---|---
LSTM | GSP interpolation | model first
C1D-LSTM | GSP interpolation | model first
SGC-LSTM | GSP interpolation | model first
TGC-LSTM | GSP interpolation | model first
STGCN |  $1$-hop interpolation | interpolation first

###  V-D Figures of Merit

The prediction performance was evaluated by the root mean square error (RMSE)
and the mean absolute error (MAE):

$\displaystyle{\rm RMSE}=\frac{1}{T_{\rm t}}\sum_{t=1}^{T_{\rm
t}}\left(\sqrt{\frac{1}{N}\sum_{n=1}^{N}(e_{n}^{t})^{2}}\right),$ (24)

$\displaystyle{\rm MAE}=\frac{1}{T_{\rm t}}\sum_{t=1}^{T_{\rm
t}}\left(\frac{1}{N}\sum_{n=1}^{N}|e_{n}^{t}|\right),$ (25)

where $T_{\rm t}$ is the number of test samples and $e_{n}^{t}$ is the
prediction error of the $t^{\rm th}$ test sample and $n^{\rm th}$ node. In the
noisy setup, the mean absolute percentage error (MAPE) was also evaluated:

$\displaystyle{\rm MAPE}=\frac{1}{T_{\rm t}}\sum_{t=1}^{T_{\rm
t}}\left(\frac{1}{N}\sum_{n=1}^{N}\frac{|e_{n}^{t}|}{|x_{n}^{t}|}\right)\times
100\%\,.$ (26)

The error metrics MAE and RMSE have the same units as the data of interest,
but RMSE is more sensitive to large errors, whereas MAE tends to treat more
uniformly the prediction errors.

###  V-E Experimental Setup

In the applications described in Sections IV-A and IV-B, $75\%$, $50\%$, and
$25\%$ from the $N$ nodes in $\mathcal{V}$ were selected to compose the set
$\mathcal{S}$ using a greedy method of [96], called E-optimal design, with the
set $\mathcal{F}$ corresponding to the first $M$ smallest Laplacian
eigenvalues. This choice of $\mathcal{F}$ relies on the smoothness of the
underlying GS, that is, nodes near to each other are assigned with similar
values. The same sampling sets were used for both supervised and semi-
supervised training. All the experiments were conducted with a time window of
length $\tau=10$. The prediction length was $p=1$ and $p=3$ samples ahead for
the GSOD dataset, that is, $1$ day and $3$ days, respectively, and $p=1$ and
$p=6$ samples to SeattleLoop, that is $5$ and $30$ minutes, respectively.

The datasets were split into: $70\%$ for training, $20\%$ for validation, and
$10\%$ for test. Batch size was set to $T_{\rm b}=40$ and the learning rate
was $10^{-4}$, with step decay rate of $0.5$ after every $10$ epochs. Training
was stopped after $100$ epochs or $5$ non-improving validation loss epochs.
The input of the model was normalized by the maximum value in the training
set. The model was trained by the RMSprop [97] with PyTorch default parameters
[98]. The network was implemented in PyTorch 1.4.0 and experiments were
conducted on a single NVIDIA GeForce GTX 1080.

###  V-F Results: Supervised Application

TABLE IV and TABLE V show the MAE and RMSE in the supervised application. The
proposed method outperformed all competitors in virtually all scenarios. When
the sample size decreases, the performance gap increases compared to the
benchmarks. On the GSOD dataset, the SG-GRU performed much better than the
other strategies. We can see that, as the temperature GS is approximately
($\mathcal{F}$,$\epsilon$)-bandlimited with small $\epsilon$, the SG-GRU
successfully captures spatial correlations by predicting the GSs’ frequency
content.

TABLE IV: MAE and RMSE of supervised prediction applied to the GSOD dataset |  | $M=0.75N$ | $M=0.50N$ | $M=0.25N$
---|---|---|---|---
| Methods | MAE | RMSE | MAE | RMSE | MAE | RMSE
$p$=1 | SG-GRU | $\boldsymbol{1.66}$ | $\boldsymbol{2.21}$ | $\boldsymbol{1.73}$ | $\boldsymbol{2.28}$ | $\boldsymbol{1.74}$ | $\boldsymbol{2.31}$
LSTM | 2.37 | 3.11 | 2.35 | 3.09 | 2.52 | 3.31
C1D-LSTM | 2.32 | 3.02 | 2.40 | 3.15 | 2.66 | 3.49
SGC-LSTM | 3.15 | 4.15 | 3.20 | 4.25 | 3.23 | 4.28
STGCN | 2.20 | 2.98 | 2.44 | 3.29 | 2.40 | 3.22
TABLE V: MAE and RMSE of supervised prediction applied to the SeattleLoop dataset |  | $M=0.75N$ | $M=0.50N$ | $M=0.25N$
---|---|---|---|---
| Methods | MAE | RMSE | MAE | RMSE | MAE | RMSE
$p$=1 | SGGRU | $\boldsymbol{2.79}$ | $\boldsymbol{4.16}$ | $\boldsymbol{3.02}$ | $\boldsymbol{4.59}$ | $\boldsymbol{3.38}$ | $\boldsymbol{5.40}$
LSTM | 3.15 | 4.79 | 3.64 | 5.59 | 4.45 | 7.03
C1D-LSTM | 3.25 | 4.95 | 3.70 | 5.70 | 4.49 | 7.08
SGC-LSTM | 3.59 | 5.57 | 3.97 | 6.14 | 4.60 | 7.26
TGC-LSTM | 3.03 | 4.59 | 3.54 | 5.45 | 4.40 | 6.98
STGCN | $\boldsymbol{2.79}$ | 4.32 | 3.11 | 4.82 | 3.65 | 6.10

###  V-G Results: Semi-supervised Application

The loss function in (16) was used for training the SG-GRU and the LSTM-based
methods. For the STGCN, the interpolation of the target GS in (19) was
replaced by the $1$-hop interpolation. TABLE VI and TABLE VII show the result
of the SG-GRU and the competing approaches on the SeattleLoop and GSOD
datasets, respectively. Fig. 3b plots the outputs of the SG-GRU and LSTM
methods, in the second semester of 2013 over the ground-truth signal, for a
weather station out of the sampling set, highlighted in Fig. 3a, considering a
situation with $50\%$ of known nodes. The SG-GRU outperformed the competing
methods in the GSOD dataset. Since temperature GSs are highly smooth in the
graph domain, the GSP interpolation, which is based on the assumption that the
GS is bandlimited, provides good reconstruction. The energy of SeattleLoop
dataset, on the other hand, is not as concentrated as the GSOD dataset,
leading to a larger reconstruction error. Even with this limitation on the
prior smoothness assumption, the SG-GRU outperformed the STGCN combined with
$1$-hop interpolation and the TGC-LSTM combined with GSP interpolation when
the sampling set size is $25\%$ or $50\%$ of the total number of nodes. It is
worth mentioning that the STGCN and the TGC-LSTM are learning models designed
specifically for traffic forecasting. When the horizon of prediction is $30$
minutes, then the SG-GRU achieved the smallest errors among all other methods.
This could be due to simultaneous ST features extraction by the SGRU module.
Fig. 4 depicts the predicted speed by SG-GRU, STGCN, and TGC-LSTM for an
unknown sensor with $p=1$, $M=0.50N$, and during the day $11/24/2015$. As can
be seen, the SG-GRU was able to better fit many points in the speed curve. It
is worth mentioning that, despite the STGCN having poorly fitted the curve in
Fig. 4b, it actually achieved higher accuracy on the known samples.

TABLE VI: MAE and RMSE of semi-supervised prediction applied to the GSOD dataset |  | $M=0.75N$ | $M=0.50N$ | $M=0.25N$
---|---|---|---|---
| Methods | MAE | RMSE | MAE | RMSE | MAE | RMSE
$p$=1 | SG-GRU | $\boldsymbol{1.77}$ | $\boldsymbol{2.38}$ | $\boldsymbol{1.88}$ | $\boldsymbol{2.53}$ | $\boldsymbol{2.06}$ | $\boldsymbol{2.76}$
LSTM | 2.35 | 3.03 | 2.41 | 3,16 | 2.72 | 3.54
C1D-LSTM | 1.83 | 2.44 | 2.00 | 2.65 | 2.24 | 2.97
SGC-LSTM | 2.75 | 3.66 | 2.84 | 3.76 | 3.01 | 3.97
STGCN | 2.34 | 3.2 | 3.75 | 5.02 | 6.92 | 8.65
$p$=3 | SG-GRU | $\boldsymbol{2.84}$ | $\boldsymbol{3.76}$ | $\boldsymbol{2.90}$ | $\boldsymbol{3.85}$ | $\boldsymbol{2.99}$ | $\boldsymbol{3.94}$
LSTM | 2.88 | 3.83 | 2.95 | 3.92 | 3.04 | 4.03
C1D-LSTM | 2.88 | 3.84 | 2.96 | 3.92 | 3.05 | 4.03
SGC-LSTM | 3.12 | 4.15 | 3.16 | 4.20 | 3.28 | 4.36
STGCN | 3.33 | 4.40 | 4.28 | 5.53 | 6.95 | 8.48

(a) US weather stations from GSOD dataset.

(b) Predicted temperature on a single sensor.

Figure 3: (a) Graph of sensors in the GSOD dataset. The known ($50\%$) and
unknown ($50\%$) nodes are colored by blue and gray, respectively. The red
node, which does belong to $\mathcal{S}$, indicates the weather station whose
temperature predictions are shown in (b); and (b) output of the SG-GRU and the
LSTM over the ground-truth temperature in the $2^{\rm nd}$ semester of 2013
measured by the node highlighted in red in (a).

(a) SG-GRU and TGC-LSTM predictions.

(b) SG-GRU and STGCN predictions.

Figure 4: Predicted signal of the sensor i005es16920 using a subset with $50\%$ of the nodes for the SG-GRU, TGC-LSTM, and STGCN. The evaluated sensor was absent in the sampling set $\mathcal{S}$. TABLE VII: MAE and RMSE of semi-supervised approaches applied to the SeattleLoop dataset |  | $M=0.75N$ | $M=0.50N$ | $M=0.25N$
---|---|---|---|---
| Methods | MAE | RMSE | MAE | RMSE | MAE | RMSE
$p$=1 | SG-GRU | 2.98 | 4.60 | $\boldsymbol{3.53}$ | $\boldsymbol{5.55}$ | $\boldsymbol{4.50}$ | $\boldsymbol{7.28}$
LSTM | 3.06 | 4.73 | 3.61 | 5.66 | 4.56 | 7.34
C1D-LSTM | 3.09 | 4.77 | 3.67 | 5.74 | 4.61 | 7.4
SGC-LSTM | 3.46 | 5.38 | 3.86 | 5.99 | 4.65 | 7.44
TGC-LSTM | 3.01 | 4.61 | 3.64 | $\boldsymbol{5.55}$ | 4.82 | 7.75
STGCN | $\boldsymbol{2.88}$ | $\boldsymbol{4.65}$ | 3.72 | 6.46 | 5.67 | 10.3
$p$=6 | SG-GRU | $\boldsymbol{3.87}$ | $\boldsymbol{6.18}$ | $\boldsymbol{4.18}$ | $\boldsymbol{6.61}$ | $\boldsymbol{4.88}$ | $\boldsymbol{7.77}$
LSTM | 3.96 | 6.34 | 4.31 | 6.81 | 4.98 | 7.94
C1D-LSTM | 3.96 | 6.37 | 4.3 | 6.83 | 5.02 | 7.97
SGC-LSTM | 4.12 | 6.63 | 4.44 | 6.98 | 5.03 | 7.98
TGC-LSTM | 4.91 | 7.89 | 5.17 | 8.23 | 8.29 | 12.5
STGCN | 4.54 | 6.82 | 4.56 | 7.90 | 6.11 | 10.8

###  V-H Results: Noise-corrupted Application

In many real situations, sensors’ measurements can be contaminated with noise,
which may worsen forecasting accuracy. Therefore, to deal with these
situations, it is important to develop robust algorithms. Consider a GS
$\mathbf{x}$ with standard deviation $\sigma_{x}$ and a measurement Gaussian
noise, uncorrelated across both time and graph-domain, $\boldsymbol{\eta}$
with standard deviation $\sigma_{\eta}$. The noisy GS is
$\tilde{\mathbf{x}}=\mathbf{x}+\boldsymbol{\eta}$ if the whole network is
measured or
$\tilde{\mathbf{x}}_{\mathcal{S}}=\boldsymbol{\Psi}_{\mathcal{S}}\mathbf{x}+\boldsymbol{\eta}$,
if only the subset $\mathcal{S}$ is measured.

To evaluate the robustness of the proposed learning scheme, both SeattleLoop
and GSOD datasets were corrupted by additive Gaussian noise with zero mean and
standard deviation (std) $\sigma_{\eta}=0.5\sigma_{x}$ and
$\sigma_{\eta}=0.1\sigma_{x}$, where $\sigma_{x}$ is the std of the entire
dataset: $10^{\rm o}{\rm C}$ for GSOD dataset and $12.74$ miles$/$h for
SeattleLoop dataset. In this experiment, nodes were not sampled and only the
capability of handling noisy input was evaluated. TABLE VIII and TABLE IX show
MAE, RMSE, and MAPE888Temperatures in the GSOD dataset were converted to
Fahrenheit before computing MAPE to avoid division by zero. of the forecasting
models respectively evaluated on $100$ and $30$ simulations of each of these
noisy scenarios.

In the GSOD dataset, the proposed model achieved reasonable error levels in
the presence of noisy measurements: for instance, MAE and RMSE increased $9\%$
and $7\%$ in comparison with the supervised situation with $M=0.75N$ when the
additive noise has std $\sigma_{\eta}=0.1\sigma_{x}$. Many GS denoising
approaches are based on attenuating high frequencies of the GS [99, 100]. The
SGRU module of the proposed model promotes the smoothness of the predicted GS
similarly: it runs a predictive algorithm over a restricted subset of the
graph frequency content, $\mathcal{F}$, and thereafter computes the inverse
GFT considering only this restricted subset,

In the SeattleLoop dataset, the MAE and RMSE evaluated on the proposed model
increased $4\%$ and $2\%$, respectively, in comparison with the supervised
situation with $M=0.75N$ when the additive noise has std
$\sigma_{\eta}=0.1\sigma_{x}$, This is a highly acceptable result, even though
the STGCN achieved lower errors.

TABLE VIII: MAE, RMSE and MAPE ($\%$) of forecasting applied to the GSOD with noise corruption |  | $\sigma_{\eta}=0.1\sigma_{x}$ | $\sigma_{\eta}=0.5\sigma_{x}$
---|---|---|---
| Methods |  MAE |  RMSE |  MAPE |  MAE |  RMSE |  MAPE
$p$=1 | SG-GRU |  $\boldsymbol{1.81}$ |  $\boldsymbol{2.36}$ |  $\boldsymbol{7.70}$ |  $\boldsymbol{2.01}$ |  $\boldsymbol{2.61}$ |  $\boldsymbol{8.52}$
LSTM |  1.98 |  2.59 |  8.55 |  2.11 |  2.75 |  9.70
C1D-LSTM |  1.90 |  2.49 |  8.07 |  2.03 |  2.65 |  8.65
SGC-LSTM |  2.94 |  3.89 |  13.9 |  2.95 |  3.91 |  13.9
STGCN |  2.19 |  2.94 |  10.3 |  2.66 |  3.48 |  12.3
$p$=3 | SG-GRU |  $\boldsymbol{2.85}$ |  $\boldsymbol{3.79}$ |  $\boldsymbol{13.41}$ |  $\boldsymbol{2.89}$ |  $\boldsymbol{3.83}$ |  $\boldsymbol{13.5}$
LSTM |  2.88 |  3.83 |  13.6 |  2.93 |  3.88 |  13.8
C1D-LSTM |  2.86 |  3.8 |  13.4 |  2.91 |  3.85 |  13.6
SGC-LSTM |  3.18 |  4.21 |  15.2 |  3.16 |  4.2 |  15.2
STGCN |  3.17 |  4.23 |  15.6 |  3.21 |  4.25 |  15.7
TABLE IX: MAE, RMSE and MAPE ($\%$) of forecasting applied to the SeattleLoop with noise corruption  |  | $\sigma_{\eta}=0.1\sigma_{x}$ | $\sigma_{\eta}=0.5\sigma_{x}$
---|---|---|---
| Methods |  MAE |  RMSE |  MAPE |  MAE |  RMSE |  MAPE
$p$=1 | SG-GRU |  2.91 |  4.27 |  13.7 |  3.13 |  4.66 |  16.1
LSTM |  3.21 |  4.85 |  18.2 |  3.45 |  5.20 |  19.4
C1D-LSTM |  3.30 |  5.05 |  19.4 |  3.48 |  5.30 |  20.5
SGC-LSTM |  3.96 |  6.28 |  28.3 |  4.07 |  6.44 |  29.9
TGC-LSTM |  2.88 |  4.24 |  13.0 |  3.19 |  4.74 |  14.2
STGCN |  $\boldsymbol{2.63}$ |  $\boldsymbol{3.85}$ |  $\boldsymbol{11.3}$ |  $\boldsymbol{3.08}$ |  $\boldsymbol{4.39}$ |  $\boldsymbol{12.9}$
$p$=6 | SG-GRU |  $\boldsymbol{3.74}$ |  $\boldsymbol{6.02}$ |  26.3 |  $\boldsymbol{3.84}$ |  $\boldsymbol{6.19}$ |  26.1
LSTM |  3.99 |  6.35 |  28.6 |  4.07 |  6.47 |  27.1
C1D-LSTM |  3.99 |  6.39 |  28.6 |  4.07 |  6.49 |  28.1
SGC-LSTM |  4.56 |  7.27 |  34.6 |  4.58 |  7.29 |  34.7
TGC-LSTM |  3.79 |  6.09 |  $\boldsymbol{26.1}$ |  3.92 |  6.28 |  $\boldsymbol{25.3}$
STGCN |  3.77 |  6.18 |  26.4 |  4.00 |  6.33 |  26.0

###  V-I Results: Missing-value Application

Another common problem in real time-series datasets are missing values, which
could occur due to sensor’s malfunctioning or failure in transmission. To
evaluate the performance of the SG-GRU in this situation, $10\%$ of both
SeattleLoop and GSOD datasets were randomly set to NaN. Before applying the
forecasting methods, each NaN value, $x^{t}_{n}$, was interpolated by the
$1$-hop interpolation in (23). TABLE X and TABLE XI show the numerical results
of this scenario considering two forecasting horizons on the GSOD and
SeattleLoop datasets, respectively. The forecasting accuracy decreases when
there are missing values, as expected. For instance, in the GSOD dataset, MAE
and RMSE increased $6\%$ and $4\%$ in comparison with the supervised situation
with $M=0.75N$. In the SeattleLoop dataset, MAE increases about $10\%$ whereas
the RMSE decreases about $8\%$. The GFT in the proposed model (and also in
combination with the LSTM-based models) tends to smooth the output signal,
reducing large deviations and consequently the RMSE. Nonetheless it can
slightly increase the forecasting error across many nodes, leading to the
increase in MAE.

TABLE X: MAE, RMSE and MAPE ($\%$) of forecasting applied to the GSOD dataset with $10\%$ of missing values |  | $p=1$ | $p=3$
---|---|---|---
| Methods |  MAE |  RMSE |  MAPE |  MAE |  RMSE |  MAPE
| SG-GRU |  $\boldsymbol{1.75}$ |  $\boldsymbol{2.3}$ |  $\boldsymbol{7.53}$ |  $\boldsymbol{2.87}$ |  $\boldsymbol{3.77}$ |  $\boldsymbol{13.1}$
| LSTM |  2.56 |  3.41 |  12.3 |  2.94 |  3.93 |  14.1
| C1D-LSTM |  2.53 |  3.36 |  11.9 |  2.92 |  3.9 |  13.9
| SGC-LSTM |  3.51 |  4.81 |  20.2 |  3.22 |  4.34 |  16.8
| STGCN |  2.10 |  2.87 |  9.97 |  3.22 |  4.31 |  15.2
TABLE XI: MAE, RMSE and MAPE ($\%$) of forecasting applied to the SeattleLoop dataset with $10\%$ of missing values |  | $p=1$ | $p=6$
---|---|---|---
| Methods |  MAE |  RMSE |  MAPE |  MAE |  RMSE |  MAPE
| SG-GRU |  3.10 |  $\boldsymbol{3.85}$ |  4.60 |  $\boldsymbol{6.17}$ |  14.9 |  $\boldsymbol{23.7}$
| LSTM |  3.37 |  4.07 |  5.08 |  6.49 |  19.6 |  27.3
| C1D-LSTM |  3.44 |  4.06 |  5.23 |  6.49 |  19.6 |  27.3
| SGC-LSTM |  4.01 |  4.58 |  6.34 |  7.31 |  16.3 |  35.1
| TGC-LSTM |  3.15 |  3.91 |  4.70 |  6.26 |  13.5 |  24.8
| STGCN |  $\boldsymbol{2.60}$ |  3.91 |  $\boldsymbol{3.95}$ |  6.26 |  $\boldsymbol{11.8}$ |  24.9

###  V-J Computational Cost and Efficiency

In the SeattleLoop Dataset, the epoch duration of SG-GRU was, on average,
$8.5$ s, whereas the more complex approaches, TGC-LSTM and STGCN, took around
$40$ s and $84$ s per epoch, respectively. In the GSOD dataset, which is much
shorter than the SeattleLoop, the average epoch duration of SG-GRU, LSTM, and
STGCN were $0.20$ s, $0.25$ s, and $2.5$ s, respectively. TABLE XII shows the
average training time, including pre-processing and data preparation, as well
as test phases for the 3 semi-supervised scenarios applied on the SeattleLoop
and GSOD datasets, with $p=1$. The SG-GRU required more epochs to converge
than STGCN, but it still trains faster than the STGCN and also than the other
competing approaches.

TABLE XII: Average computational time in seconds | SeattleLoop | GSOD
---|---|---
Methods | Training | Test | Training | Test
SG-GRU | 414.68 | 4.89 | 11.18 | 0.01
LSTM | 1134.0 | 5.63 | 30.74 | 0.01
C1D-LSTM | 1319.0 | 5.90 | 36.90 | 0.03
SGC-LSTM | 2770.3 | 6.10 | 39.89 | 0.04
TGC-LSTM | 1027.1 | 5.48 | - | -
STGCN | 725.58 | 12.6 | 83.92 | 0.12

###  V-K Final Remarks on the Results

The consistently better results obtained by the SG-GRU for the GSOD dataset
steem from the smoothness of the temperature GS with respect to the graph
domain; SG-GRU relies on the assumption of bandlimited GSs. Therefore, SG-GRU
is a promising approach to predict spatially smooth GSs. It is worth
mentioning that the choice of the adjacency matrix is fundamental for a good
performance, since it eventually defines the smoothness of the GSs. In the
SeattleLoop dataset, which is not really smooth, the SG-GRU outperformed both
the STGCN and the LSTM-based approaches when the sample size was small and the
prediction time horizon was $30$ minutes, thus indicating that the SG-GRU can
capture ST dependencies by taking the network frequency content into account.
Moreover, SG-GRU has low computational cost and can be boosted with more
recurrent or fully connected layers, when sufficient computational resources
are available.

##  VI Conclusion

This work presented a new deep learning technique for jointly forecasting and
interpolating network signals represented by graph signals. The proposed
scheme embeds GSP tools in its basic learning-from-data unit (SG-GRU cell),
thus merging model-based and deep learning approaches in a successful manner.
Indeed, the proposal is able to capture spatiotemporal correlations when the
input signal comprises just a small sample of the entire network.
Additionally, the technique allows reliable predictions when input data is
noisy or some values are missing by enforcing smoothness on the output
signals. As future works, we envisage the use of the proposed SG-GRU as part
of an anomaly detector in network signals, in which the anomalous sensors’
measurements are characterized by large deviations from the neighboring
sensors.

## References

  * [1] A. G. Salman, B. Kanigoro, and Y. Heryadi, “Weather forecasting using deep learning techniques,” in _International Conference on Advanced Computer Science and Information Systems (ICACSIS)_ , Oct 2015, pp. 281–285.
  * [2] Y. Lv, Y. Duan, W. Kang, Z. Li, and F. Wang, “Traffic flow prediction with big data: A deep learning approach,” _IEEE Transactions on Intelligence Transportation Systems_ , vol. 16, no. 2, pp. 865–873, Apr 2015\.
  * [3] S. M. Smith, K. L. Miller, G. Salimi-Khorshidi, M. Webster, C. F. Beckmann, T. E. Nichols, J. D. Ramsey, and M. W. Woolrich, “Network modelling methods for FMRI,” _NeuroImage_ , vol. 54, no. 2, pp. 875–891, Jan 2011.
  * [4] L. Li, K. Ota, and M. Dong, “When weather matters: IoT-based eletric. load forecasting for smart grid,” _IEEE Communications Magazine_ , vol. 55, no. 10, pp. 46–51, Oct 2017.
  * [5] E. Racah, C. Beckham, T. Maharaj, S. Ebrahimi Kahou, M. Prabhat, and C. Pal, “ExtremeWeather: A large-scale climate dataset for semi-supervised detection, localization, and understanding of extreme weather events,” in _Advances in Neural Information Processing Systems 30 (NIPS)_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds. Curran Associates, Inc., Dec 2017, pp. 3402–3413.
  * [6] B. M. Williams and L. A. Hoel, “Modeling and forecasting vehicular traffic flow as a seasonal ARIMA process: Theoretical basis and empirical results,” _Journal of Transportation Engineering_ , vol. 129, no. 6, pp. 664–672, Oct 2003.
  * [7] L. Cai, Z. Zhang, J. Yang, Y. Yu, T. Zhou, and J. Qin, “A noise-immune Kalman filter for short-term traffic flow forecasting,” _Physica A: Statistical Mechanics and its Applications_ , vol. 536, p. 122601, Dec 2019.
  * [8] S. R. Chandra and H. Al-Deek, “Predictions of freeway traffic speeds and volumes using vector autoregressive models,” _Journal of Intelligence Transportation Systems_ , vol. 13, no. 2, pp. 53–72, Apr 2009.
  * [9] O. Kramer and F. Gieseke, “Short-term wind energy forecasting using support vector regression,” in _Soft Computer Models in Industrial and Environmental Applications, $6^{\rm th}$ International Conference (SOCO)_, vol. 87. Springer, 2011, pp. 271–280.
  * [10] G. Leshem and Y. Ritov, “Traffic flow prediction using AdaBoost algorithm with random forests as a weak learner,” _Journal International Journal of Intelligence Technology_ , vol. 2, pp. 1305–6417, Jan 2007.
  * [11] W. Huang, G. Song, H. Hong, and K. Xie, “Deep architecture for traffic flow prediction: deep belief networks with multitask learning,” _IEEE Transactions on Intelligence Transportation Systems_ , vol. 15, no. 5, pp. 2191–2201, Apr 2014.
  * [12] Yuhan Jia, Jianping Wu, and Yiman Du, “Traffic speed prediction using deep learning method,” in _IEEE International Conference on Intelligence Transportation Systems (ITSC)_ , Dec 2016, pp. 1217–1222.
  * [13] A. Salman, Y. Heryadi, E. Abdurahman, and W. Suparta, “Weather forecasting using merged long short-term memory model (LSTM) and autoregressive integrated moving average (arima) model,” _Journal of Comput. Science_ , vol. 14, pp. 930–938, July 2018.
  * [14] S. Liang, L. Nguyen, and F. Jin, “A multi-variable stacked long-short term memory network for wind speed forecasting,” in _IEEE International Conference on Big Data_ , Dec 2018, pp. 4561–4564.
  * [15] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MA, USA,USA: MIT press, 2016.
  * [16] X. SHI, Z. Chen, H. Wang, D.-Y. Yeung, W.-k. Wong, and W.-c. WOO, “Convolutional LSTM network: A machine learning approach for precipitation nowcasting,” in _Advances in Neural Information Processing Systems 28_ , C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, Eds. Curran Associates, Inc., 2015, pp. 802–810.
  * [17] H. Yu, Z. Wu, S. Wang, Y. Wang, and X. Ma, “Spatiotemporal recurrent convolutional networks for traffic prediction in transportation networks,” _Sensors_ , vol. 17, no. 7, pp. 1–16, June 2017.
  * [18] C. J. Huang and P. H. Kuo, “A deep CNN-LSTM model for particulate matter (${\rm Pm}_{2.5}$) forecasting in smart cities,” _Sensors_ , vol. 18, no. 2220, July 2018.
  * [19] Y. Wu, H. Tan, L. Qin, B. Ran, and Z. Jiang, “A hybrid deep learning based traffic flow prediction method and its understanding,” _Transportation Research Part C: Emerging Technologies_ , vol. 90, pp. 166–180, May 2018.
  * [20] W. Li, W. Tao, J. Qiu, X. Liu, X. Zhou, and Z. Pan, “Densely connected convolutional networks with attention LSTM for crowd flows prediction,” _IEEE Access_ , vol. 7, pp. 140 488–140 498, Sep 2019.
  * [21] M. Henaff, J. Bruna, and Y. LeCun, “Deep convolutional networks on graph-structured data,” _eprint arXiv:1506.05163_ , June 2015.
  * [22] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _Advances in Neural Inf. Process. Systems 29_ , D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, Eds. Curran Associates, Inc., 2016, pp. 3844–3852.
  * [23] R. Levie, F. Monti, X. Bresson, and M. M. Bronstein, “CayleyNets: Graph convolutional neural networks with complex rational spectral filters,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 1, pp. 97–109, Nov 2019.
  * [24] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _eprint arXiv:1609.02907_ , Sep 2016.
  * [25] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in _6 th International Conference on Learning Representations (ICLR)_, Fev 2018, pp. 1–16.
  * [26] L. Zhao, Y. Song, M. Deng, and H. Li, “Temporal graph convolutional network for urban traffic flow prediction method,” _eprint arXiv:1811.05320_ , Dec 2018.
  * [27] J. Zhang, X. Shi, J. Xie, H. Ma, I. King, and D.-Y. Yeung, “GaAN: gated attention networks for learning on large and spatiotemporal graphs,” _eprint arXiv:1803.07294_ , Mar 2018.
  * [28] B. Wang, X. Luo, F. Zhang, B. Yuan, A. L. Bertozzi, and P. J. Brantingham, “Graph-based deep modeling and real time forecasting of sparse spatio-temporal data,” _eprint arXiv:1804.00684_ , Apr 2018.
  * [29] X. Geng, Y. Li, L. Wang, L. Zhang, Q. Yang, J. Ye, and Y. Liu, “Spatiotemporal multi-graph convolution network for ride-hailing demand forecasting,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, pp. 3656–3663, July 2019.
  * [30] W. Chen, L. Chen, Y. Xie, W. Cao, Y. Gao, and X. Feng, “Multi-range attentive bicomponent graph convolutional network for traffic forecasting,” _eprint arXiv:1911.12093_ , Nov 2019.
  * [31] Z. Cui, K. Henrickson, R. Ke, and Y. Wang, “Traffic graph convolutional recurrent neural network: A deep learning framework for network-scale traffic learning and forecasting,” _IEEE Transactions on Intelligence Transportation Systems_ , Nov 2019.
  * [32] B. Yu, H. Yin, and Z. Zhu, “Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting,” in _Proceedings of the 27 th International Joint Conference on Artificial Intelligence_, July 2018, pp. 3634–3640.
  * [33] M. Wang, B. Lai, Z. Jin, Y. Lin, X. Gong, J. Huang, and X. Hua, “Dynamic spatio-temporal graph-based cnns for traffic prediction,” _eprint arXiv:1812.02019_ , Mar 2020.
  * [34] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang, “Graph WaveNet for deep spatial-temporal graph modeling,” _International Joint Conference on Artificial Intelligence (IJCAI)_ , pp. 1907–1913, Aug 2019.
  * [35] K. Lee and W. Rhee, “DDP-GCN: Multi-graph convolutional network for spatiotemporal traffic forecasting,” _eprint arXiv:1905.12256_ , May 2019\.
  * [36] Z. Diao, X. Wang, D. Zhang, Y. Liu, K. Xie, and S. He, “Dynamic spatial-temporal graph convolutional neural networks for traffic forecasting,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, July 2019, pp. 890–897.
  * [37] S. Fang, Q. Zhang, G. Meng, S. Xiang, and C. Pan, “GstNet: Global spatial-temporal network for traffic flow prediction,” in _Proceedings of the 28 th International Joint Conference on Artificial Intelligence (IJCAI)_, Aug 2019, pp. 10–16.
  * [38] C. Song, Y. Lin, S. Guo, and H. Wan, “Spatial-temporal synchronous graph convolutional networks: A new framework for spatial-temporal network data forecasting,” https://github.com/wanhuaiyu/STSGCN/blob/master/paper/AAAI2020-STSGCN.pdf, 2020\.
  * [39] S. Guo, Y. Lin, N. Feng, C. Song, and H. Wan, “Attention based spatial-temporal graph convolutional networks for traffic flow forecasting,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, July 2019, pp. 922–929.
  * [40] Z. Cui, R. Ke, Z. Pu, and Y. Wang, “Deep bidirectional and unidirectional LSTM recurrent neural network for network-wide traffic speed prediction,” _eprint arXiv:1801.02143_ , Jan 2018.
  * [41] X. Wang, C. Chen, Y. Min, J. He, B. Yang, and Y. Zhang, “Efficient metropolitan traffic prediction based on graph recurrent neural network,” _arXiv:1811.00740_ , Nov 2018.
  * [42] B. Liao, D. McIlwraith, J. Zhang, T. Chen, C. Wu, S. Yang, Y. Guo, and F. Wu, “Deep sequence learning with auxiliary information for traffic prediction,” _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 537–546, July 2018.
  * [43] B. Liao, J. Zhang, M. Cai, S. Tang, Y. Gao, C. Wu, S. Yang, W. Zhu, Y. Guo, and F. Wu, “Dest-ResNet: A deep spatiotemporal residual network for hotspot traffic speed prediction,” in _Proceedings of the $26^{\rm th}$ ACM International Conference on Multimedia_, Oct 2018, pp. 1883–1891.
  * [44] D. Han, J. Chen, and J. Sun, “A parallel spatiotemporal deep learning network for highway traffic flow forecasting,” _International Journal of Distributed Sensor Networks_ , vol. 15, no. 2, Feb 2019.
  * [45] Z. Lv, J. Xu, K. Zheng, H. Yin, P. Zhao, and X. Zhou, “LC-RNN: A deep learning model for traffic speed prediction,” in _International Joint Conference on Artificial Intelligence (IJCAI)_ , July 2018, pp. 3470–3476.
  * [46] Y. Lin, X. Dai, L. Li, and F. Y. Wang, “Pattern sensitive prediction of traffic flow based on generative adversarial framework,” _IEEE Transactions on Intelligence Transportation Systems_ , vol. 20, no. 6, pp. 2395–2400, Aug 2019.
  * [47] H. F. Yang, T. S. DIllon, and Y. P. P. Chen, “Optimized structure of the traffic flow forecasting model with a deep learning approach,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 28, no. 10, pp. 2371–2381, Oct 2017.
  * [48] K. Zhang, L. Zheng, Z. Liu, and N. Jia, “A deep learning based multitask model for network-wide traffic speed prediction,” _Neurocomputing_ , vol. 396, pp. 438 – 450, Apr 2020.
  * [49] X. Ouyang, C. Zhang, P. Zhou, H. Jiang, and S. Gong, “Deepspace: An online deep learning framework for mobile big data to understand human mobility patterns,” _arXiv:1610.07009_ , Mar 2016.
  * [50] M. Khodayar and J. Wang, “Spatio-temporal graph deep neural network for short-term wind speed forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 10, no. 2, pp. 670–681, Apr 2019.
  * [51] M. Khodayar, O. Kaynak, and M. E. Khodayar, “Rough deep neural architecture for short-term wind speed forecasting,” _IEEE Transactions on Ind. Inform._ , vol. 13, no. 6, pp. 2770–2779, July 2017.
  * [52] Q. Zhu, J. Chen, L. Zhu, X. Duan, and Y. Liu, “Wind speed prediction with spatio-temporal correlation: A deep learning approach,” _Energies_ , vol. 11, p. 705, Mar 2018.
  * [53] C. Y. Zhang, C. L. Chen, M. Gan, and L. Chen, “Predictive deep boltzmann machine for multiperiod wind speed forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 6, no. 4, pp. 1416–1425, July 2015.
  * [54] M. Khodayar, S. Mohammadi, M. E. Khodayar, J. Wang, and G. Liu, “Convolutional graph autoencoder: A generative deep neural network for probabilistic spatio-temporal solar irradiance forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 2, pp. 571–583, Apr 2020.
  * [55] S. Rasp and S. Lerch, “Neural networks for postprocessing ensemble weather forecasts,” _Monthly Weather Review_ , vol. 146, no. 11, pp. 3885–3900, Oct 2018.
  * [56] C. Li, Z. Cui, W. Zheng, C. Xu, and J. Yang, “Spatio-temporal graph convolution for skeleton based action recognition,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , Apr 2018.
  * [57] C. Wu, X.-J. Wu, and J. Kittler, “Spatial residual layer and dense connection block enhanced spatial temporal graph convolutional network for skeleton-based action recognition,” in _The IEEE International Conference on Comput. Vision (ICCV) Workshops_ , Oct 2019.
  * [58] C. Si, W. Chen, W. Wang, L. Wang, and T. Tan, “An attention enhanced graph convolutional LSTM network for skeleton-based action recognition,” in _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * [59] S. Gadgil, Q. Zhao, E. Adeli, A. Pfefferbaum, E. V. Sullivan, and K. M. Pohl, “Spatio-temporal graph convolution for functional MRI analysis,” _eprint arXiv:2003.10613_ , Mar 2020.
  * [60] H. Huang, X. Hu, J. Han, J. Lv, N. Liu, L. Guo, and T. Liu, “Latent source mining in FMRI data via deep neural network,” _IEEE International Symposium on Biomedical Imaging (ISBI)_ , pp. 638–641, July 2016.
  * [61] Y. Li, D. Tarlow, M. Brockschmidt, and R. Zemel, “Gated graph sequence neural networks,” _arXiv:1511.05493_ , Sep 2015.
  * [62] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Processing Magazine_ , vol. 30, no. 3, pp. 83–98, May 2013.
  * [63] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 808–828, May 2018.
  * [64] S. K. Narang, A. Gadde, E. Sanou, and A. Ortega, “Localized iterative methods for interpolation in graph structured data,” in _IEEE Global Conference on Signal and Information Processing (GlobalSIP)_ , Dec 2013, pp. 491–494.
  * [65] J. Chen, H.-r. Fang, and Y. Saad, “Fast approximate knn graph construction for high dimensional data via recursive lanczos bisection.” _Journal of Machine Learning Research_ , vol. 10, no. 9, Sep 2009.
  * [66] _Interpolation of graph signals using shift-invariant graph filters_ , Aug 2015\.
  * [67] _Multiscale Wavelets on Trees, Graphs and High Dimensional Data: Theory and Applications to Semi Supervised Learning._ , 2010.
  * [68] I. Pesenson, “Sampling in Paley-Wiener spaces on combinatorial graphs,” _Transactions of the American Mathematical Society_ , vol. 360, no. 10, pp. 5603–5627, May 2008.
  * [69] A. Anis, A. Gadde, and A. Ortega, “Towards a sampling theorem for signals on arbitrary graphs,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2014, pp. 3864–3868.
  * [70] S. Chen, F. Cerda, P. Rizzo, J. Bielak, J. H. Garrett, and J. Kovačević, “Semi-supervised multiresolution classification using adaptive graph filtering with application to indirect bridge structural health monitoring,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 11, pp. 2879–2893, Mar 2014.
  * [71] S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete signal processing on graphs: Sampling theory,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 24, pp. 6510–6523, Aug 2015.
  * [72] S. K. Narang and A. Ortega, “Local two-channel critically sampled filter-banks on graphs,” _IEEE International Conference on Image Processing_ , pp. 333–336, Sep 2010.
  * [73] M. J. Spelta and W. A. Martins, “Online temperature estimation using graph signals,” _XXXVI Simpósio Brasileiro de Telecomunicaçoes e Processamento de Sinais (SBrT)_ , pp. 154–158, Sep 2018.
  * [74] M. J. M. Spelta and W. A. Martins, “Normalized LMS algorithm and data-selective strategies for adaptive graph signal estimation,” _Signal Processing_ , vol. 167, p. 107326, Feb 2020.
  * [75] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for bandlimited graph signals using graph spectral proxies,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 14, pp. 3775–3789, Mar 2016\.
  * [76] L. F. O. Chamon and A. Ribeiro, “Greedy sampling of graph signals,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 1, pp. 34–47, Jan 2018.
  * [77] D. I. Shuman, M. Faraji, and P. Vandergheynst, “Semi-supervised learning with spectral graph wavelets,” _Proceedings of the International Conference on Sampling Theory and Applications (SampTA)_ , May 2011.
  * [78] _Active Semi-Supervised Learning Using Sampling Theory for Graph Signals_ , ser. KDD ’14. Association for Computing Machinery, 2014.
  * [79] S. Chen, R. Varma, A. Singh, and J. Kovačević, “Signal recovery on graphs: Fundamental limits of sampling strategies,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 2, no. 4, pp. 539–554, Oct 2016.
  * [80] S. Parisot, S. I. Ktena, E. Ferrante, M. Lee, R. Guerrero, B. Glocker, and D. Rueckert, “Disease prediction using graph convolutional networks: Application to autism spectrum disorder and Alzheimer’s disease,” _Medical Image Analysis_ , vol. 48, pp. 117 – 130, Aug 2018.
  * [81] G. Cheung, W. Su, Y. Mao, and C. Lin, “Robust semi-supervised graph classifier learning with negative edge weights,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 4, no. 4, pp. 712–726, Mar 2018.
  * [82] W. Huang, A. G. Marques, and A. Ribeiro, “Collaborative filtering via graph signal processing,” in _$25^{\rm th}$ European Signal Processing Conference (EUSIPCO)_, Sep 2017, pp. 1094–1098.
  * [83] W. Huang, A. G. Marques, and A. R. Ribeiro, “Rating prediction via graph signal processing,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 19, pp. 5066–5081, Aug 2018.
  * [84] K. Manohar, B. W. Brunton, J. N. Kutz, and S. L. Brunton, “Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns,” _IEEE Control Systems Magazine_ , vol. 38, no. 3, pp. 63–86, June 2018.
  * [85] B. Xu, H. Shen, Q. Cao, Y. Qiu, and X. Cheng, “Graph wavelet neural network,” _7 th International Conference on Learning Representations (ICLR)_, pp. 1–13, Apr 2019.
  * [86] D. Romero, V. N. Ioannidis, and G. B. Giannakis, “Kernel-based reconstruction of space-time functions on dynamic graphs,” _IEEE Journal on Selected Topics in Signal Processing_ , vol. 11, no. 6, pp. 856–869, Sep 2017.
  * [87] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 7, pp. 1644–1656, Jan 2013.
  * [88] S. K. Narang and A. Ortega, “Downsampling graphs using spectral theory,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Mar 2011, pp. 4208–4211.
  * [89] P. Lorenzo, S. Barbarossa, and P. Banelli, _Chapter 9 - Sampling and Recovery of Graph Signals_ , P. M. Djurić and C. Richard, Eds. Academic Press, 2018.
  * [90] K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, “On the properties of neural machine translation: Encoder-decoder approaches,” _eprint arXiv:1409.1259_ , Sep 2014.
  * [91] A. Graves, “Generating sequences with recurrent neural networks,” _eprint arXiv:1308.0850_ , Aug 2013.
  * [92] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” _arXiv:1412.3555_ , Dec 2014\.
  * [93] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” _eprint arXiv:1312.6203_ , Dec 2013.
  * [94] National climactic data center. (202001, Feb. 1). [Online]. Available: ftp://ftp.ncdc.noaa.gov/pub/data/gsod
  * [95] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs: Frequency analysis,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 12, pp. 3042–3054, Apr 2014.
  * [96] B. J. Winer, “Statistical principles in experimental design,” _McGraw-Hill Book Company_ , 1962.
  * [97] G. Hinton, N. Srivastava, and K. Swersky, “Lecture 6e rmsprop: Divide the gradient by a running average of its recent magnitude,” 2012.
  * [98] Pytorch RMSprop. (2020, Feb. 1). [Online]. Available: https://pytorch.org/docs/stable/optim.html?highlight=rms#torch.optim.RMSprop
  * [99] D. I. Shuman, P. Vandergheynst, and P. Frossard, “Chebyshev polynomial approximation for distributed signal processing,” in _2011 International Conference on Distributed Computing in Sensor Systems and Workshops (DCOSS)_ , June 2011, pp. 1–8.
  * [100] N. Tremblay and P. Borgnat, “Subgraph-based filterbanks for graph signals,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 15, pp. 3827–3840, Mar 2016.