Multivariate Time Series Forecasting with Latent Graph Inference
Authors: Victor Garcia Satorras
Syama Sundar Rangapuram
Tim Januschowski
Venue: ICLR 2022 Submitted (2022)
Keywords: Time Series Forecasting
Graph Neural Networks
Graph Inference
Multivariate Time Series

Abstract: This paper introduces a new architecture for multivariate time series forecasting that simultaneously infers and leverages relations among time series. We cast our method as a modular extension to univariate architectures where relations among individual time series are dynamically inferred in the latent space obtained after encoding the whole input signal. Our approach is flexible enough to scale gracefully according to the needs of the forecasting task under consideration. In its most straight-forward and general version, we infer a potentially fully connected graph to model the interactions between time series, which allows us to obtain competitive forecast accuracy compared with the state-of-the-art in graph neural networks for forecasting. In addition, whereas previous latent graph inference methods scale O(N^2) w.r.t. the number of nodes N (representing the time series), we show how to configure our approach to cater for the scale of modern time series panels. By assuming the inferred graph to be bipartite where one partition consists of the original N nodes and we introduce K nodes (taking inspiration from low-rank-decompositions) we reduce the time complexity of our procedure to O(NK). This allows us to leverage the dependency structure with a small trade-off in forecasting accuracy. We demonstrate the effectiveness of our method for a variety of datasets where it performs better or very competitively to previous methods under both the fully connected and bipartite assumptions.

Under review as a conference paper at ICLR 2022
MULTIVARIATE TIME SERIES FORECASTING WITH LA- TENT GRAPH INFERENCE
Anonymous authors Paper under double-blind review
ABSTRACT
This paper introduces a new architecture for multivariate time series forecasting that simultaneously infers and leverages relations among time series. We cast method as modular extension to univariate architectures where relations among individual time series are dynamically inferred in the latent space obtained after encoding the whole input signal. Our approach flexible enough to scale gracefully according to the needs of the forecasting task under consideration. In its most straight-forwart and general version, we infer potentially fully connected graph to model the interactions between time series, which allows us to obtain competitive forecast accuracy compared with the state-of-the-art in graph neural networks for forecasting. In addition, whereas previous latent graph inference methods scale ( w.r.. the number of nodes N (representing the time series), we show how to configure our approach to cater for the scale of modern time series panels. By assuming the inferred graph to be bipartite where one partition consists of the original N nodes and introduce K nodes (taking inspiration from low-rank-decompen decompositions) we reduce the time complexity our procedure ONN)). This allows to leverage the dependency structure with small trade-off in forecasting accuracy. We demonstrate the effectiveness of our method for variety of datasets where performs better or very competitively to previous methods under both the fully connected and bipartite assumptions.
INTRODUCTION
Time Series Forecasting (TSF) has been widely studied due to its practical significance in wide variety of applications such as climate modelling (Mudelsee 2019). supply chain management in retail (Larson] 2001, Böse et al. 2017). market analysis in finance (Andersen et al., 22055), traffic control (Li et al. 2017) and medicine (Kaushik et al., 2020). Petropoulos er al. (2020) provide non-systematic overview of further applications. In TSF, given sequence of data points indexed over time, we aim to estimate its future values based on previously observed data. Data often multivariate, meaning that multiple variables vary time, each variable may not only depend its own historical values, but also on other variables' past. Efficiently modelling the dependencies between these variables is still an open problem.
Multivariate Time Series (MTS) methods aim to leverage the dependencies between variables in order improve the forecasting accuracy. Some classical MTS forecasting algorithms such as Vector Autoregression (VAR) (Litkepohh] 20055 or Gaussian Processes (Roberts et al. 2013) only consider linear dependencies among variables. A natural way model non linelr dependencies in deep learning is via Graph Neural Networks (GNNs) (Bruna et al., 2013) Defferrard et al.) 2016 Kipf and Welling 2016). In fact, GNNs have been successfully applied for the Multivariate Time Series forecasting (Li et al. 2017, Yu et al. 2017, Chen et al. (020a), leveraging the relations among different series. But these methods require a pre-defined adjacency matrix which may only be available in some specific datasets, for example, traffic datasets, where it can be constructed from the spatial structure of a city. More recently, a family of methods that do not require a pre-defined adjacency matrix have been proposed (Wu et al. 2019 Franceschi al.) 2019 Wu et al., 2020 Shang et al. 202T). In this case, a latent graph representation is inferred while forecasting, allowing operate larger variety of datasets. Our work belongs to this category.

---

Under review as a conference paper at ICLR 2022
However, inferring all pairwise relations may at higher computational cost, the number of relations scales quadratically O(N2) w.r.. the number of nodes. This limits the scalability to large datasets. Additionally, the above mentioned latent graph inference methods perform message passing update at every time step iteration which can also be expensive.
To overcome these limitations, we propose Time Series forecasting that is more efficient than previous algorithms while achieving better or competitive performance. We cast the latent graph inference as a modular and easy-to-implement extension to current univariate models. The graph is dynamically inferred per inputted data stream allowing more flexible representation than a static graph for the whole dataset. Additionally, we optionally reduce the complexity from ((V2) (Fully Connected Assumption) to O(N (Bipartite Assumption) where K < N for small trade off in performance.
latent graph inference algorithm for Multivariate
BACKGROUND
2.1 TIME SERIES FORECASTING
In time series forecasting we want to estimate a future time series X++1:T given its past Xtaa where to < s indexes over time, and (optionally) some context information c. For the multivariate case assume the time series is composed of N variates at time such that Xto:T {xi,to:T..... XN.to:T) E RNxT +1+ this section we distinguish two main categories of time series forecasting methods, Global Univariate and Multivariate.
Global Univariate methods: In this case we only use the past of each univariate to predict its future. However, the model weights 0n are shared across all univariate time series. More formally:
1: = fx( C;; Ow)
where & denotes the estimated values, i E [1....,N) N) indexes over multivariates and f.)-) is the estimator function with learnable parameters 0, shared across time series. Conditioning on the past of each univariate may limit the performance of the forecasting algorithm compared to multivariate ones. Despite that, it simplifies the design of fo and already provides reasonable results. A popular example of univariate models would be (Salinas et al., 20200.
Multivariate methods: Multivariate methods condition on all past data (all N variates) and directly predict the multivariate target. More formally:
81+1.7 fm (Xto:,
Different variables may be correlated and/or depend on the same con-foure con-founders. For example, in retail forecasting, PPE masks and antibacterial soaps jointly increased in demand during the early days of the COVID-19 pandemic. In traffic forecasting, an increase of the outcome traffic flow in given neighborhood may result in an increase of the income traffic flow on another one. Modelling these dependencies may improve the forecasting accuracy, but it may come at a cost of higher complexity and hence expensive algorithms, specially when trying model all pairwise interactions between variates.
2.2 GRAPH NEURAL NETWORKS
Graph Neural Networks (GNNs) Bruna et al. 2013 Defferrard et al. 2016 Kipf and Welling 2016 operate directly graph structured data. They have gained lot of attention in the last years due to their success in a large variety of domains which benefit from modelling interactions between different nodes/entities. In the context of multivariate time series, GNNs can be used to model the interactions between time series. In this work we consider the type of GNN introduced by (Gilmer al. 2017). Given graph g (V, with nodes U; V and edges eij E we define graph convolutional layer as:
mij de (h) h!)
Σ mij JEA(()
h on (h/, m;)
m;
Where d, and dn are the edge and node functions, usually approximated as Multi Layer Perceptrons (MLPs), hy E af is the of-dimmnsional embedding of node i layer / and m,y is the edge
2

---

Under review as a conference paper at ICLR 2022
embedding that propagates information from node 2/ to U.. A GNN is constructed by stacking multiple of these Graph Convolutional Layers h GCLjh' EJ. Additionally, in Equation S we include any E (0, 1) which is a scalar value that performs the edge inference or attention over the neighbors similarly to Veličković al. (2017). As done in (Satorras et al. 2021), we choose this value to be computed as the output of a function aij = 0a (m:y) where Φa is composed of just a linear layer followed by a sigmoid activation function.
3 RELATED WORK
Time series forecasting has been extensively studied in the past due to its practical significance
with a number of recent overview articles available (et al. 2020, [Bnidiis et al.] 2020; Lim and Zohren 202T. Traditionally, most classical methods are univariate nature (see e.g., Hyndman and Athanasopoulos (2017) for an overview). While some of these have multi-variate extensions (e.g.. ARMA and VARMA models), they are limited by the amount of related time series information they can incorporate. Dynamic factor models (Geweke) 1977, Wang et al., 2019a) are fore-runners of a family of models that has recently received more attention, the so-called global models (Januschowski and Hyndman 2022). These global models estimate their parameters
et al. 2019; Montero-Manne over an entire panel of time series, so thereby taking advantage of cross time series learning, but still produce univariate forecast. Many such global models have been proposed building on the
main neural network architectures like RNNs (Salinas et al.) 2020 Liberty et al. 2020; Bandara
ct al. 2019). CNNs (Wen et al. 2017 Chen et al. 2020b), Transformers (Li et al. 2019] Lim et al.
2021: Eisenach et al. (020)) and also combining classical probabilistic models with deep learning
(Rangapuram al.J 2018 Kurle al. 2020 de Bézenac al.) 2020). However, these global models
do not explicitly model the relationship between the time series in the panel.
Most recently, global multi variate forecasting models have received attention. in particular models that attempt to capture the relationship of the time series via a multi-variat likelihood (Rasul et al.] 2020 2021 de Bézenac al. 2020) Salmas al. 2019). Here, attempt to capture the multi- variate nature of many modern forecasting problems primarily by using a multi-variate time series as input. For this, a natural way model and exploit the relationship between time series is via Graph Neural Networks (Bruna et al. 2013 Defferrard et al.] (016)] Kipf and Welling, 2016) which have been successfully applied between entities/nodes can benefit the prediction task. Even in those cases where edges are not explicitly provided in the dataset, attention or latent graph can be inferred from the node embeddings such that GNNs can still leverage the structure of the data. Some examples of latent graph inference or attention are (Wang et al.] 2019b) in point clouds, (Franceschi et al. 2019) in semi-supervised graph classification, (Ying et al. 2018) in hierarchical graph representation learning, (Kipf et al. 2018) in modelling dynamical systems, segmentation, (Garcia and Bruna 2017 Kossen al. 202T) in image classification, (Cranmer et al.
wide variety of deep learning domains where exploiting relations
Kazi et al. 2020 in zero-shot learning and 3D point cloud
2020) in inferring symbolic representations and (Fuchs et al. 2020 Satorras et al.) 2021) in molecular property prediction.
In Multivariate Time Series (MTS) forecasting we can leverage dependencies between time series by exchanging information among them. Lai et al. (2018): (hhih al. (2019) are some of the first deep learning approaches designed to exploit those pair- wise dependencies. More recent methods, (Li al. 2017 Yu all 2017 Seo et al. 2018 Zhao et al. 2019) are built in the intersection of Graph Neural Networks and time series forecasting but they require a pre-defined adjacency matrix. Lately, new methods that infer a latent graph from the node embeddings have been introduced in MTS
forecasting (Wu et al. 2020; Shang et al. 2021] Cao et al. 202T). and in MTS anomaly detection (Zhang et al.] 2020; Deng and Hood 2021). These methods can be applied to any dataset even when explicitly defined adjacency matrix. But this comes with a limitation, inferring the edges can be expensive since those scale quadratically O(N2) w.r.t the number of nodes/variables N
there is not
or ((N3) in (Coo et al. 1)
Our approach is related to (Wu et al. 2020 Shang et al. 2021). but in contrast a) Our latent graph is dynamically inferred for each inputted data stream instead of a static graph for the whole dataset which allows a more flexible graph representation. b) It is modular since it can be added as an extension to standard univariate methods after encoding the input signal Xtoit this also makes the
3

---

Under review as a conference paper at ICLR 2022
or
FC-GNN
BP-GNN
Xt,to
Zi
Zi
X1,t+1:77
fenc
fdec
AGG
Figure 1: Illustration of the presented method under fully connected and bipartite graph assumptions.
graph operation cheaper since the message exchange is only done at time step t. c) Optionally reduces the number of edges from O(N2) to ((NK) where K < N with the bipartite assumption.
METHOD
In Section 2.1 discussed two different sorts of forecasting methods: Global Univariate Models and Multivariate Models. We cast our multivariate algorithm as a modular extension of the univariate case from Equation We can break down the Univariate Model in main steps fu Sene d such fune
fuc, Xi,++::7, where fene encodes the input signal Xi.to.6 (and optionally some context information c;) into an embedding z; and faee estimates the future signal from this embedding. In method, we include multivariate aggregation module AGG in between Sene and faee that propagates information among nodes in the latent space {zı,- ... }. This aggregation module takes as input the embedding z, (xne(Xi,tott C;), and outputs modified embedding à AGG(z) (z) where information has been propagated among nodes. These new embeddings are then passed as input to the decoder fdee- The resulting algorithm is:
that Xi.foo Z₁
Univariate Encoder 2; = feno Ci)
(4)
Multivariate extension
AGG(z))
(5)
Univariate Decoder & = (2,)
(6)
Notice the overall model is multivariate but fene and faee remain univariate, allowing modular extension from current univariate methods to multivariate. Additionally, in contrast to the most recent methods (Wu al. 2020 Shang et al. 202T) our model does not propagate information among nodes at every time step [o...... but only in the AGG module. This makes the algorithm cheaper since the message propagation step is usually an expensive operation. Additionally, we experimented including a unique identifier of the time series in the encoded signal as context information c, = id, which resulted in a significant improvement in accuracy as we will show in the experiments section. In the following subsections we propose two different Graph Neural Network solutions for the AGG module, each has its benefits and weaknesses regarding performance and scalability.
Fully connected graph assumption FC-GNN
From the two assumptions, this is the most straightforward to implement given standard GNN, results in good performance, but its complexity scales O(N²) with respect to the number of nodes N. Despite this, it still resulted in faster computation times than previous O(N2) methods. Here, we directly use the Graph Neural Network defined in Section 2 the aggregation module AGG. We assume fully connected graph G = {V, E} where each time series embedding z; is associated to a node of the graph
V. specifically, the embedding z; is provided as the input hę to the GNN from Equation 3 then the GNN runs for L layers and the output node embedding h is provided as the input 2, of the decoder. Moreover, despite the fully connected assumption, the GNN infers attention weights a,j (0, (eq. B for each edge and input sample. This is equivalent to dynamically inferring a latent graph where edges are associated to soft values in the range (0, 1). In the experiments section, name this model Fully Connected Graph Neural Network (FC-GNN). In summary, the
4

---

Under review as a conference paper at ICLR 2022
aggregation module AGG under this FC-GNN setting, is able model non-linear relations among time series by exchanging pairwise messages mij (eq. 3 among them, these messages are non-linear embeddings of the edges and they are "gated" multuplied by an inferred scalar value aij.
Bipartite graph assumption BP-GNN
In the previous section presented a GNN method that exchanges information among time series under a fully connected graph assumption which computationally scales O(N2). In this section we introduce bipartite graph assumption which reduces the complexity to ON) K being parameter of choice K < N. To accomplish this, we define a bipartite graph g = () (V,,,E) with two sets of nodes Y and U. Here y is set of N nodes corresponding to the N dimensions of the multivariate time series and U is set of K auxiliary nodes. Each time series embedding z; y. Additionally, we define embeddings = {uI..... uK
associated one of the nodes associated to the auxiliary set U. Edges & interconnect all nodes between the two subsets {V,u), but there are no connections among nodes of the same subset. This results 2NK edges, with a computation complexity of O)N
The algorithm works in the following way. We input into the GNN the union of the two node subsets == Ju.. Specifically, the input embedding h defined in Equation 3 the concatenation of the time series embeddings z with the auxiliary node embeddings u (i.e. he z 2Ilu) where u are free learnable parameters initialized as Gaussian noise. Then, messages follow an asynchronous schedule, first they are sent from the time series nodes to the auxiliary nodes y - UI next the other way around + V. Notice the auxiliary nodes U/ do not contain any information about the data until they have received messages from V this is why we chose this schedule. A visual illustration of the whole algorithm is presented in Figure M where the bipartite graph is depicted the top right of the image (BP-GNN).
We have conceptually defined the Bipartite Graph Neural Network. In Table | we introduce the equa- tions that formally define it extension of the standard GNN Equation 3 Notice that it can be simply formulated as two steps process where the indexes i, j belong to each one of the subsets 21 or V depending the direction of the messages (V 4 or 21 > V). Additionally, we used different learn- able parameters between Step I and Step in the modules de, &a and on since it resulted in better accuracy than sharing parameters. In Appendix [A.] we provide the BP-GNN equations from Table in its full form. Following, we define the adjacency matrices corresponding the two message passing steps (assuming all inference parameters aij- 1):
Step y + 2
Step U + y
Equation
Equation where
where
Eeu & jEY
iey & j jeu
Table 1: BP-GNN formulation.
ONNN A
ONKK Ox*K
A2
ONNN OKXN
1NNK Ox*к
Ā= A2A1
KNXN 0KxN
ONKK Oσ*
TRX
The two adjacency matrices A₁ and Ag are non-symmetric and they define a directed message passing scheme. A₁ refers y/ 21 and A₂ refers to + y. After running these two propagation steps, all time series nodes y' will have received messages from all other time series nodes ) as in the fully connected graph case, but the number of messages have been reduced (iff S K < N/2) due the factorization of the adjacency matrix. The first N rows and columns in matrix Á from Equation 7 represent the sum of all paths that communicate the time series nodes y among them after the two updates AzAr
Architecture details
As explained in this Section 4 our method composed of three main modules, the encoder fenes the decoder fdee and the aggregation function AGG. We choose to use relatively simple networks as encoder and decoder. The decoder fiee is defined as a Multi Layer Perceptron (MLP) with a single hidden layer and residual connection in all experiments. The encoder fene is also defined as an MLP for METR-LA, PEMS-BAY and our synthetic datasets and as a Convolutional Neural Network (CNN) for the other datasets since these require larger encoding length and the translation equivariance of CNNs showed to be more beneficial. The encoder fenes first encodes the input signal Xioto to embedding vector by using the mentioned MIP CNN, and then concatenates unique
5

---

Under review as a conference paper at ICLR 2022
identifier (c, id) to the obtained embedding vector resulting in 2;- c, could (optionally) include additional context information it was provided in the dataset. The combination of these simple networks with proposed aggregation module AGG fully defines our model. All architecture details are explained in detail in Appendix AA.2
The aggregation module AGG defined under two different assumptions in, the Fully Connected Graph assumption (FC-GNN) and the Bipartite Graph assumption (BG-GNN). In both cases the GNN is fully parametrized by the networks de. Φn and da. &e consists of a two layers MLP, ch is a one layer MLP with a skip connection from he input to the output and do is just a linear layer followed by a Sigmoid activation function. All these architecture choices are explained in more detail in Appendix A.2, In all experiments, the loss was computed as the Mcan Squared Error between the estimated values X++1:7 from Equation and the ground truth as f = 1(x++1:T; x(-ITT).
EXPERIMENTS
5.1 DATASETS AND BASELINES
We first evaluate our method in METR-LA and PEMS-BAY datasets from (i et al. 2017) which record traffic speed statistics on
ANodes
# Samples
Context length
Pred. length
METR-LA PEMS-BA BAY
34,272
12
325
52,116
12
12
the highways of Los Angeles county and the Bay Area respectively. We also evaluate in the publicly available Solar-Energe Traffic, Electricity and Exchange-Raue Rate. Specifica- tions for each dataset presented in Table
Solar-Fnergy
137
$2,560
168
I
Traffic
862
17,544
168
I
Electricity
26,304
Exchange-Rate
7,588
168 I
Table 2: Dataset specifications.
2 where #Nodes is the number of time series in the panel, #Samples is the number of time steps in each time series, Context length is the length of the input window and Pred. Length is the length of the predicted window. We compare to a variety of baselines including previous works and variations of our proposed method. We distinguish the following three main types of baselines (Univariate, Multivariate with known graph and Multivariate with graph inference):
Univariate: In this case single model is trained for all time series but they treated independently without message exchange among them. These baselines include a simple linear Auto Regressive model (AR) and variation of our FC-GNN method that we called NE-GNN where all edges have been removed such that the aggregation function AGG becomes equivalent to multilayer perceptron defined by On in Equation
Multivariate with known graph: These methods require previously defined graph, therefore they are restricted to those datasets where an adjacency matrix can be pre-defined (e.g. METR-LA, PEMS-BAY). From this group we compare to DCRNN (Li al. 7017p. STGCN (Yu al. 2017
and HIRA-BGCN (Chen et al. 2020a).
Multivariate with graph inference or attention: These methods exchange information among different time series by simultaneously inferring relations among them or by attention mechanisms, From this group we compare to LDS (Franceschi et al., 2019). LST-Skip (Lai et al. 2018), TPA- LSTM (Shih et al., 2019). MTGNN (Wu et al.J 2020) and GTS (Shang et al. 2021). Graph WaveNet (Wu et al., 2019) also belongs to this group but unlike the others it jointly uses a pre- defined adjacency matrix. Comparisons to NRI (Kipf et al., 2018) can be found in previous
literature (Shang al. 2021: Zügner et al. 2021). We additionally include a variation of our FC-GNN without unique node ids, we denote it by (w/o id) beside the model name. GTS numbers have been obtained by averaging over 3 runs its official implementation.
https://qitts
Sonn/hhaoshangcs//CC
5.2
MAIN RESULTS
In this section we evaluate our method in METR-LA and PEMS FEMS-BAY datasets. For this experiment we used the training setup from GTS (Shang et al. 2021) which uses the dataset partitions and evaluation metrics originally proposed in (Li et al. 2017). The model has been trained by minimizing the Mean Absolute Error (MAE) between the predicted and the ground truth samples. The reported

---

Under review as a conference paper at ICLR 2022
METR-LA
PEMS-BAY
15 min
30 min
60 min
15 min
30 min
60 min
MAE RMSE MAPE
|MAE
RMSE
MAPE
MAE RMSE
MAPE
MMA
RMSE
MAPE
MAE
RMSE
MAPE
|AA RMSE MAPE
Lincar / AR
3.81
8.80
9.13%
4.94
11.14
12.17%
6.30
12.91
16.72%
1.59
3.41
3.27%
22.15
4.87
4.77
%
22.97
6.65
7.03%
DCRNN
2.77
5.38
7.30%
3.15
6.45
8.80%
3.60
7.60
10.50%
1.38
2.95
2.90%
1.74
3.97
3.90%
2.07
4.74
4,90%
STGCN
5.74
3.47
7.24
9.40
12.7%
1.36
2.96
2.9%
4.27
4.2%
2.49
5.69
MRA-BGCN
2.67
5.12
6.8%
3.06
6.17
8.3%
3.49
7.30
10.0%
1.29
2.72
2.9%
1.61
3.67
3.8%
1.91
4.46
4.6%
Graph WaveNet*
2.69
5.15
6.90%
3.07
6.22
8.37%
3.53
7.37
10.01%
1.30
2.74
2.73%
1.63
3.70
3.67%
1.95
4.52
4.63%
LDDS
2.75
5.35
7.1%
3.14
6.45
8.6%
3.63
7.67
10.34%
1.33
2.81
2.8%
1.67
3.80
3.8%
1.99
4.59
4.8 %
MIGGN
5.18
6.86%
3,05
8.19%
7.23
9.87%
1.32
2.77%
3.74
3.69%
1.94
4.49
4.53%
GTS
5.19 2.64
6.79%
3.06
6,30
8.24%
3.56
7.55
9.95%
1,35
2.84
2.85%
1.67
3.82
3.80%
1.96
4,53
4.62%
Ablation study
NE- GNN (uoo id)
2.80
5.73
7.50%
3.40
7.15
9.74%
4.22
8.79
13.06%
1.40
3.03
2.92%
1.85
4.25
4.21%
2.39
5.50
5.93%
FC-GNN (woo
5.65
7.39%
3.36
9.59%
8.64
12.70%
1.39
2.88%
4.18
4.11%
2.32
5.35
5.71%
NE. NE.GNN
2.69
5.57
7.21%
3.14
6.74
9.01%
3.62
7.88
10.94%
1.36
2.88
2.86%
1.72
3.91
3.93 %:
2.07
4.79
5.04%
Our models
FC GNN
2.60
5.19
6.78%
2.95
6.15
8.14%
3.35
7.14
9.73%
1.33
2,82
2.79%
1.65
3.75
3.72%
1,93
4.46
4.53%
BP-GNN (K=1)
2.64
5.37
7.07%
3.02
6.42
8.46%
7.32
9.91%
1.33
2,80%
3.78
3.75%
1.94
4.46
4.57%
Table 3: Benchmark on METR-LA and PEMS-BAY datasets. Mean Absolute Error (MAE), Root Mean Squared Error (RMSH) and Mean Absolute Percentage Error (MAPE) are reported for different time horizons (15, 30, 60) minutes. Results have been averaged over 5 runs.
metrics are MAE, Root Mean Squared Error (RMSE) and Mean Absolute Percentage Error (MAPE) from (Li et al. 2017). All metrics have been averaged over 5 runs. All our models (FC-GNN, BP-GNN and NE-GNN) contain 2 graph convolutional layers, 64 features in the hidden layers. Swish activation functions (Ramachandran et al., 2017) and have been trained with batch size 16. The number of auxiliary nodes for BP-GNN was set to K 4. Time experiments report the average forward pass in seconds for batch size 16 in a Tesla V100-SXM GPU. Further implementation details are provided in Appendix B.
Results are reported in tables 11 and results with standard deviations are reported in Appendix B.3
FC-GNN outperforms other methods in most metrics while being computationally cheaper than previous works. On the other hand, BP-GNN performs very competitively w.r.t. previous works (even outperforming all previous methods in some metrics) but with a vast decrease in computation. Furthermore, these performances are achieved without providing the structure in- formation of the city (i.e. pre-defined adjacency) unlike in those methods that require it or in Graph Wavenet that optionally uses it. Additionally, notice the performance gap between FC- GNN and NE-GNN is larger when including a unique identifier of the nodes. This means the network can better leverage the information exchange among nodes when they are uniquely identified, but it still benefits from the message passing scheme when they are not uniquely identified (w/o id) thanks the dynamical inference. Time results are presented in Table 4 BP-GNN the most efficient algorithm in both METR- BAY by large margin. FC-GNN is also more efficient than previous methods both datasets but it is still limited by the O(N2) scalability. In Section 5.3 we will evaluate in larger dataset showing that BP-GNN IS even more efficient w.r.t. to other methods as the graph becomes larger thanks to its linear scalability ONN. In this timing benchmark (Table we included all methods from the previous Table TT that have a publicly available implementation in METR-LA and PEMS-BAY datasets.
Forward Time (s)
METR-LA PEMS-BAY
and PEMS-
Linear
.0002
.0022
DCRNN
.2559
.754
Graph WaveNet MIGNN GTS
.0500
.0673
.0160
.0371
.0869
.087
NE-GNN GNN BP-GNN (K=4)
.0033
.0047
.0108
.0533
.0044
.0066
Table 4: Forward time in seconds for different methods.
5.3 SINGLE STEP FORECASTING
In this section, we evaluate our method in the publicly avilable Solar-Energy, Traffic, Electricity and Exchange- Rate datasets. In contrast to METR-LA and PEMS-BAY, these datasets do not contain spatial information from which a graph can be pre-defined, therefore, methods that rely on a known graph are not directly applicable. We use the same training settings as (Lai et al. 2018 Shih et al.
7

---

Under review as a conference paper at ICLR 2022
Dataset
Solar Energy
Trailic
Electricity
Exchange
Horizon
Horizon 12
Horizon
Horizon
Mothods
Metries #Top?
3
12
24
24
3
6
24
3
12
24
RSE
(0) .2435 .3790
.5911
.6699
.5991
.6118
.6252
.6293
0995
.1035
.1050
.1054
.0228
0229
.0353
0445
AR
CORR
(0)
9710 9263
.8107
.5314
.7752
.7568
.7544
.7519
8845
8632
.8591
8595
.934
.9566
.956
.9557
RNN-GRU
RSE
(0)
.143 .259
.3254
.4643
4777
.4893
.950
.4973
0864
.0931 .1007
.1007
0226
0280
.0356
.0449
CORR
(1)
9843
9690
.9677
.8870
8721
.690
8614
.5888
9283
9135
-9077
9119
9735
9658
.9111
.9554
LST-skip
RSE
(0)
1843 2559
.3254
4663
4777
4893
.950
4973
0864
0031 1007
1007
0226
0280
.0356
0449
CORR
(U)
.9433
.9667
88870
8690
_6614
9283
9935
.907
9119
9911
.9544
TRA-LSTM
RSE
(3)
.1803
.2347 .3234 4389
4487
4658
.644
.766
0823
0916 0964
1006
.0174
.0241
.0341 .0444
CORR
(4)
.9550
.942
.9877
.901
8812
.817
.8717
8629
9439
9937
.2500
.9333
.9900
.9799
.9564
.9811
MTGNN
RSE
(5)
.1778
.2488
.3109
.4270
4162
.4754
4461
.4535
0744
.0878
.9166
.0953
.094
.0259
.0399
.0456
CORR
(6)
99852
9726
.9509
.033
8963
.8667
.8794
.8810
9474
.9316
.2778
.2344
.9786
.9708
.9551
.9372
NE-GNN
RSE
(2)
.1898
2580
.3472
.441
4212
4586
.4679
.4433
00722
.0117
.966
.994
.0175
.0244
.0338
.0447
CORR
(U)
99899
.9633
.9367
.8915
.8748
.8700
9494
.9308
.262
.9535
.9352
FC-GNN
RSE
(12)
.1651
.202
.2981
.397
4057
.4395
.4624
.4620
.0732
.0907
.0915
.1979
.0174
00455
.0344
.0450
CORR
(12)
.9876
.765
9951
9148
.850
88764
.8751
.9521
.9404
9381
9294
9772
9685
95388
.9499
BP-GNN
RSE
(11)
.1704
.257
.3072
4050
4095
.4470
4640
.664
.0740
.8998
0940
0980
0175
.0244
.03.9
.0442
CORR
(10)
.9865
9751
9522
9138
8999
8880
8744
.823
9519
9396
9945
9288
.9699
9984
.9530
.9660
Table 5: Benchmark on Solar-Energy, Traffic, Electricity and Exchange-Rate Root Relative Squared
Error (RSE) and Empirical Correlation Coefficient (CORR) are reported for different horizons {3, 6, 12, ). All results have been averaged over 5 runs. #Top2 column counts how many metrics in each row are in the top 2 (i.e. bold).
2018). Sj. In this experiment the network 2019 Wu et al. 2020) firstly proposed by TLST-skip Lai al,
is trained to predict only one time step into the future (single step forecasting) with a given horizon (3, 6, 12, or 24) by minimizing the Mean Absolute Error (MAB). All datasets have been split in 60%/20%6/20 22% for training/val/test respectively. As proposed in (La et al. 2018) we use the Root Relative Squared Error (RSE) and Empirical Correlation Coefficient (CORR) as evaluation metrics (both defined in Appendix C.2 All our models (FC-GNN, BP-GNN and NE-GNN) contain graph convolutional layers and 128 features in the hidden layers The timing results have been run with batch size 16. Further implementation details are provided in Appendix C.T Results with standard deviations provided in Appendix
Results in Table 5 are consistent with the previous experiment. FC-GNN outperforms all previous works in most metrics. There is significant improvement w.r.t. its cousin baseline NE-GNN, this demonstrates that sharing information among nodes is beneficial (except in Exchange-Rate dataset which only contains 8 nodes). BP-GNN performs better than previous methods in most metries and very close to FC-GNN. Regarding the timing results (Table C. BP-GNN is the most efficient graph inference method by large margin in most datasets. The larger the number of nodes in the dataset, the larger the computational improvement of BP-GNN w.r.t. other methods due to its linear scalability. For example, in Solar- Solar-Energy (137 nodes), BP-GNN is 1.43 times faster than FC-GNN and 1.92 times faster than MTGNN. On the other hand, in a larger dataset as Electricity (321 nodes) BP-GNN is 6.38 times faster than FC-GNN and 9.18 times faster than MTGNN. In the largest dataset, Traffic (862 nodes), BP-GNN becomes 34.58 times faster than FC-GNN and 97.59 times faster than MTGNN, although in traffic, MTGNN and FC- GNN did not fit in the GPU for a batch of 16 (also because to the O(N2) complexity). ), and we had to pass the samples in batches of which eliminates part of the GPU parallelization. Running BP-GNN in batches of 2, would result in it being 7.71 and 18.92 times faster than FC-GNN and MTGNN respectively. Finally, in such small graphs as Exchange Rate, there is no computational benefit in using the bipartite assumption since the number of edges for both the bipartite and the fully connected graphs becomes the same (for K=4 and N=8).
Exchange
Solar
Electricity
Traffic
Nodes
8
137
321
862
MIGNN
.0062
.0146
.0771
1.1808
NE-GNN
.0034
.0599
.0067
.0117
FC-GNN
.0533
.0109
.0536
.4184
BP-GNN
.0076
.0766
.0084
.0121
Table 6: Forward time in seconds for dif- ferent methods in Solar-Energe Traffic, Electricity and Exchange- datasets.
5.4 INFERRED GRAPH ANALYSIS
In this section we study the adjacency matrices inferred by our FC-GNN and BP-GNN methods. For this purpose, we created two synthetic datasets of N = 10 time series and T 10.000 timesteps each.
8

---

Under review as a conference paper at ICLR 2022
(a) "Cycle Graph Gaussians" samples each series value x from the past (t-5) of another series (i-1 mod N) from the panel. The resulting adjacency matrix is a directed cycle graph. More formally, the dataset is generated from the following gaussian distribution X,, N (5x;x-11 mod N,t-5; o?). where B 0.9 and 0.5². (b) We call the second dataset "Noisy Sinusoids", inspired by the Discrete Sine Transformation we generate arbitrary signals as the sum different sinusoids plus gaussian noise. We share the same sinusoids among the first half 1 S i
of the panel and among the second half 6 A i S 10 while sampling different gaussian noise for each signal This creates strong dependencies among the signals within each half. The expected adjacencies are plotted in Figure 2 Further details and visualizations of these two datasets are provided in Appendix D. Both BP-GNN and FC-GNN consist of single graph convolution from which adjacencies are obtained. We predict the next time step into the future and optimize the MAE during training. Since our graph inference mechanism is dynamic we average them over 10 timesteps For BP-GNN we construct the adjacencies the first N rows and columns of Á Ay (Eq, Sec. 4 Notice that small K values in BP-GNN may enforce sharing auxiliary nodes and in consequence a more dense adjacency A. Further implementation details are described in Appendix D.2
Expected Adjacency
NE-GNN
FC-GNN
BP-GNN K=10
BP-GNN K=5
SP.GNN
BP-GNN
K=1
Cycle Graph Gaussian
Adjacency
0,4023
0,4098
0,4149
0,4995
0.4955
Adjacency
Nolay Sinusoids
Figure 2: Inferred adjacency matrices and MAE losses in the proposed synthetic datasets.
Results reported in Figure FC-GNN perfectly matches the ground truth adjacencies and has the lowest MAE test loss. In the "Cycle Graph Gaussians" dataset, BP-GNN infers a similar matrix to the ground truth when provided with enough auxiliary nodes K 10, but the matrix becomes denser when reducing K, which is expected since different messages are forced to share the same auxiliary nodes. Despite that, in both datasets, BP-GNN manages obtain very competitive performance (even for small K) by inferring suboptimal matrices with denser connections. In the "Noisy Sinusoids" dataset, the benefit of increasing is smaller, this is coherent with the fact the ground truth matrix denser. In all cases the proposed models significantly outperform their non-edge variant NE-GNN. Further visualizations for different random seeds are provided in Appendix D.3
5.5
CHOOSING THE NUMBER OF AUXILIARY NODES K
0.00
GIN-2)
2.68
0.007
Forqee Smme 14
Teaff:
2.07
0.066
Sde 02-
rew--Br
Nummerr
265
0000
2.t
263
Murbbee aurtliary rasctes
Number
Nummbr radies
Figure 3: (Left) Mean Absolute Error METR-LA for different K values (10 runs average). (Middle) Running times in METR-LA. (Right) Number of edges w.r.t to the number of nodes.
In the previous section we evaluated BP-GNN for different K values, and we found small values enough to achieve good performance in the denser adjacency dataset (Noisy Sinusoids). In real world scenarios, we found a similar behavior, where results are competitive with small values. Therefore, we chose relatively small number of auxiliary nodes K 4 that resulted in good

---

Under review as a conference paper at ICLR 2022
trade-off accuracy computation. In Figure plot the MAE for different K values in METR-LA for 15 min time horizon, and we see the best performances are already achieved at small K values. In the middle plot of Figure E] we report the forward running time in (s) for batch size 16 for METR-LA as we sweep different K values. Finally, the right plot of Figure S shows the scalability of the number of edges w.r.t. the number of nodes for O(N2) (e.g. FC-GNN) and ( (i.e. BP-GNN (K=4). This illustrates the vast difference in scalability as the number of nodes increases.
6 CONCLUSIONS
We presented novel approach for multi-variate time series forecasting which lends itself to easy integration into existing univariate approaches by ways of a graph neural network (GNN) component. This GNN infers latent graph in the space of the embeddings of the univariate time series and can be inserted in many neural-networt based forecasting models. We show that this additional graph/dependency structure improves forecasting accuracy both in general and in particular when compared to the state of the art. We alleviate typical scalability concerns for GNNs via allowing the introduction of auxiliary nodes in the latent graph which we construct as a bipartite graph that bounds the computational complexity. We show that this leads to computationally favorable properties as expected while the trade-off with forecasting accuracy is manageable.
REFERENCES
Andersen, G., Bollerslev, T., Christoffersen, P., and Diebold, X. (2005). Volatility forecasting.
Bandara, K., Shi, P., Bergmeir, C., Hewamalage, H., Tran, Q. and Seaman, B. (2019). Sales demand forecast e-commerce using a long short-term memory neural network methodology. In International Conference Neural Information Processing, pages 462-474. Springer.
Benidis, K., Rangapuram, S. S., Flunkert, V., Wang, B., Maddix, D. C., Türkmen, A. C.,, Gasthaus, J., Bohlke-Schneider M., Salinas, D., Stella, L., Callot, L., and Januschowski, T. (2020). Neural forecasting: Introduction and literature overview. CoRR, abs/2004. 10240.
Böse, J..H., Flunkert, V., Gasthaus, J., Januschowski, T., Lange, D., Salinas, D., Schelter, S., Seeger, M., and Wang. Y. (2017). Probabilistic demand forecasting at scale. Proceedings of the VLDB Endowment, 10(12): (0(12):1694-1700 199--1005
Bruna, J., Zaremba, W., Szlam, A., and LeCun, (2013). Spectral networks and locally connected networks on graphs. arXiv preprint arXiv: (XXiv:3312.620:
Cao, D., Wang, Y., Duan, J., Zhang, C., Zhu, X., Huang, C., Tong, Y., Xu, B., Bai, Tong, J., et al. (2021). Spectral temporal graph neural network for multivariate time-series forecasting, arXiv preprint 1131:2203..0711.
Chen, W., Chen, L., Xie, Y., Cao, W., Gao, Y., and Feng, X. (2020a). Multi-range attentive bicomponent graph convolutional network for traffic forecasting. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3529-3536.
Chen, Y., Kang, Y., Chen, Y., and Wang, Z. (2020b). Probabilistic forecasting with temporal convolutional neural network. Neurocomputing. 399:491-501.
Cranmer, M., Sanchez-Gonzales A., Battaglia, P., Xu, R., Cranmer, K., Spergel, D., and Ho. S. (2020). Discovering symbolic models from deep learning with inductive biases. arXiv preprini arXiv:2006. / 1287.
de Bézenac, E., Rangapuram. S. S., Benidis, K.,
M., Kurle, R., Stella, L., Hasson, H., Gallinari, P., and Januschowski, T. (2020). Normalizing kalman filters for multivariate time series analysis. Advances in Neural Information Processing Systems, 33.
Defferrard, M., Bresson, X., and Vandergheynst, P. (2016). Convolutional neural networks graphs with fast localized spectral filtering. In Advances in neural information processing systems, pages 3844-3852.
10

---

Under review as a conference paper at ICLR 2022
Deng, A. and Hooi, B. (2021). Graph neural network-based anomaly detection in multivariate time series. In Proceedings of the 35th AAAI Conference on Artificial Intelligence, Vancouver, BC, Canada, pages 2-9.
Eisenach, C., Patel, Y., and Madeka, D. (2020). Mqtransformer: Multi-horizon forecasts with context dependent and feedback-aware attention. arXiv preprint arXiv:2009.1477
et al., P. (2020). Forecasting: theory and practice. arXiv preprint arXiv:2012.03854.
Franceschi, L.,, Niepert, M., Pontil, M., and He, X. (2019). Learning discrete structures for graph neural networks. In International conference on machine learning, pages 1972-1982. PMLR.
Fuchs, F. B., Worrall, D. E., Fischer, V., and Welling, M. (2020). Se (3)-transformers: 3d roto- translation equivariant attention networks. arXiv preprint arXiv:2006. 10503.
Garcia, V. and Bruna, J. (2017). Few-shot learning with graph neural networks. arXiv preprint arXiv:1711.04043.
Geweke, J. (1977). The dynamic factor analysis of economic time series. Latent variables in socio-economic models.
Gilmer, J., Schoenholz, S. S.. Riley, F., Vinyals, O., and Dahl, G. E. (2017). Neural message passing for quantum chemistry. In International conference on machine learning, pages 1263-1272. PMLR.
Hyndman, R. J. and Athanasopoulos, G. (2017). Forecasting: Principles and practice. www. otexts. org/fpp. 987507109.
Januschowski, T., Gasthaus, J., Wang. Y., Salinas, D., Flunkert, V., Bohlke-Schneiderr M., and Callot, L. (2019). Criteria for classifying forecasting methods. International Journal of Forecasting.
Kaushik, S., Choudhury, A., Sheron, P. K., Dasgupta, N., Natarajan, S., Pickett, L. A., and Dutt, V.
(2020), Ai healtheare: time-series forecasting using statistical, neural, and ensemble architectures. Frontiers in big data, 3:4.
Kazi, A., Cosmo, L., Navab, N., and Bronstein, M. (2020). Differentiable graph module (dgm) for graph convolutional networks. arXiv preprint arXiv:2002.0499)
Kipf, T., Fetaya, E.,, Wang, K.-C., Welling, M., and Zemel, R. (2018). Neural relational inference for interacting systems. In International Conference on Machine Learning, pages 2688-26977 PMLR.
classiiiitiiacccssiiictttiiio Kipf, T. N. and Welling, M. (2016). Semi-supervised classification with graph convolutional networks. arXiv preprint 06.900070.
Kossen, J., Band, N., Lyle, C., Gomez, A. N., Rainforth, T., and Gal, Y. (2021). Self-attention between datapoints: Going beyond individual input-output pairs in deep learning. arXiv preprint 1.Xiv:2106.02584,
Kurle, R., Rangapuram, S. S., de Bézenac, E., Günnemann, S., and Gasthaus, J. (2020). Deep boobblackwellieed particle filters for time series forecasting. Advances in Neural Information Processing Systems, 33.
Lai, G., Chang, W.-C., Yang, Y., and Liu, H. (2018). Modeling long-and short-term temporal patterns with deep neural networks. The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 95-104.
Larson, D. (2001). Designing and managing the supply chain: concepts, strategies, and case studies. Journal of Business Logistics, 22(1):259.
Li. S., Jin, X., Xuan, Y., Zhou, X., Chen, W., Wang, Y-X.,, and Yan, X. (2019). Enhancing the locality and breaking the memory bottleneck of transformer time series forecasting. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
Li, Y., Yu. R., Shahabi, C., and Liu, Y. (2017). Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. arXiv preprint 170.0197...

---

Under review as a conference paper at ICLR 2022
Liberty, E.,, Karnin, Z., Xiang, B., Rouesnel, L., Coskun, B.., Nallapati, R., Delgado, J., Sadoughi, A., Astashonok, A., Das, P., Balioglu, C., Charkravarty, S., Jha, M., Gaultier, P., Januschowski, T., Flunkert, V., Wang, B., Gasthaus, J., Rangapuram, S., Salinas, D., Schelter, S., Arpin, D., and Smola, A. (2020). Elastic machine learning algorithms in Amazon SageMaker. In Proceedings of the 2020 International Conference Management of Data, SIGMOD '20, New York, NY, USA. ACM.
Lim, B., Ank, S. Ó., Loeff, N., and Pfister, T. (2021). Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting.
Lim, B. and Zohren, S. (2021). Time-series forecasting with deep learning: survey. Philosoph- ical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2194):0220209
Lütkepohl, H. (2005). New introduction to multiple time series analysis. Springer Science & Business Media.
Montero-Manso, P. and Hyndman, R. (2022). Principles and algorithms for forecasting groups of time series: Locality and globality. International Journal of Forecasting.
Mudelsee, M. (2019). Trend analysis of climate time series: A review of methods. Earth-science reviews, 190:310-324
Petropoulos, F., Apiletri, D., Assimakopoulos, V., Babai, M. Z.., Barrow, D. K., Bergmeir, C.,, Bessa, R. J., Boylan, J. E., Browell, J., Carnevale, C., Castle, J. L., Cirillo, P., Clements, M. P., Cordeiro, C., Oliveira, F. C., Bacts. S. D., Dokumentov, A., Fiszeder, P.. Franses, H. Gilliland. M.. Gönül, M. S., Goodwin, P., Grossi, L.,, Frushka-Cockayner Y., Guidolin, M., Guidolin, M., Gunter, U., Guo, X., Guseo, R., Harvey, N., Hendry, D. F., Hollyman, R., Januschowski, T., Jeon, J., Jose, V. R. R., Kang, Y., Koehler, A. B., Kolassa, S., Kourentzes, N., Leva, S., Li, F., Litsiou, K., Makridakis, S., Martinez, A. B., Meeran, S., Modis, T., Nikolopoulos, K., Önkal, D., Paccagnini, A., Panapakidis, I., Pavía, J. M., Pedio, M., Pedregal, J., Pinson, P., Ramos, P., Rapach, D. E., Reade, J. J., Rostami-Tabar B., Rubaszek, M., Sermpinis, G., Shang. H. L., Spiliotis, E., Syntetos, A. A., Talagala, P. D., Talagala, T. S., Tashman, L., Thomakos, D., Thorarinsdottir. T., Todini, E., Arenas, J. R. T., Wang, X., Winkler, R. L.., Yusupova, A., and Ziel, F. (2020). Forecasting: theory and practice. arXiv preprint 20120220/..03854.
Ramachandran, P., Zoph, B., and Le, Q. V. (2017). Searching for activation functions. arXiv preprint arXiv:
Rangapuram, S. S., Seeger, M. W., Gasthaus, J., Stella, L.., Wang. Y., and Januschowski, T. (2018).
Deep state space models for time series forecasting. In Advances in Neural Information Processing Systems, pages 7785-7794.
Rasul, K., Seward, C., Schuster, I., and Vollgraf, R. (2021). Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. Proceedings of ICLR 2021.
Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U., and Vollgraf, R. (2020). Multi-variate probabilis- tic time series forecasting via conditioned normalizing flows. arXiv preprint 10Xiv:2002.06103.
Roberts, S., Osborne, M., Ebden, M., Reece, S., Gibson, N., and Aigrain, S. (2013), Gaussian pro- cesses for time-series modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371((184):20)) 10550.
Salinas, D., Bohike-Schneiderr, M., Callor, L., Medico, R., and Gasthaus, J. (2019). High-dimensional multivariate forecasting with low-rank gaussian copula processes. In Advances in Neural Informa- tion Processing Systems 32.
Salinas, D., Flunkert, V., Gasthaus, J., and Januschowski, T. (2020). Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 1113))118111911
Satorras, V. G., Hoogeboom, E., and Welling, M. (2021). E (n) equivariant graph neural networks. arXiv preprint ar 102.0:2102.09844
12

---

Under review as a conference paper at ICLR 2022
Seo, Y., Defferrard, M., Vandergheynst, P., and Bresson, X. (2018). Structured sequence modeling with graph convolutional recurrent networks. In International Conference on Neural Information Processing, pages 362-373, Springer.
Shang, C., Chen, J., and Bi, (2021). Discrete graph structure learning for forecasting multiple time series. arXiv preprint arXiv:2101.06861
Shih, S.-Y., Sun, F.-K., and Lee, H.-y. (2019). Temporal pattern attention for multivariate time series forecasting. Machine Learning, 108(8): 1421-1444
Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., and Bengio, Y. (2017). Graph attention networks. arXiv preprint arXiv:1710.1000 10903.
Wang, Y., Smola, A., Maddix, D., Gasthaus, J., Foster, D., and Januschowski, T. (2019a). Deep factors for forecasting. In International Conference on Machine Learning, pages 6607-6617.
Wang. Y., Sun, Y., Liu, Z., Sarma, S. E., Bronstein, M. M., and Solomon, J. M. (2019b). Dynamic graph enn for learning on point clouds. Acm Transactions On Graphics (tog). 38(5):1-12.
Wen, R., Torkkola, K., Narayanaswamy, B., and Madeka, D. (2017). multi-horizon quantile recurrent forecaster. arXiv preprint
Wu, Z., Pan, S., Long, G., Jiang, J., Chang, X., and Zhang, C. (2020). Connecting the dots: Multivariate time series forecasting with graph neural networks. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 753-763.
Wu, Z., Pan, S., Long, G., Jiang, J., and Zhang. C. (2019). Graph wavenet for deep spatial-temppral graph modeling. arXiv preprint arXiv: 19Xi::99060002//1
Ying. R., You, J., Morris, C.,, Ren, X., Hamilton, W. L., and Leskovec, J. (2018). Hierarchical graph representation learning with differentiable pooling. arXiv preprint arXiv:1806.0888 08804.
Yu, B., Yin, H., and Zhu, Z. (2017). Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. arXiv preprint arXiv: 1709. 748754
Zhang, R., Hao, Y. Yu, D., Chang, W.-C., Lai, G., and Yang, (2020). Correlation-anl awweee change- point detection via graph neural networks. In International Conference on Neural Information Processing, pages 555-567. Springer.
Zhao, L., Song. Y., Zhang, C., Liu, Y., Wang, P., Lin, T., Deng, M., and Li, H. (2019). T-gen: temporal graph convolutional network for traffic prediction. IEEE Transactions Intelligent Transportation Systems, 21(9):3848-388 -3588.
Zügner, D., Aubet, F.-.., Satorras, V. G., Januschowski, T., Günnemann, S., and Gasthaus, J. (2021). A study of joint graph inference and forecasting. arXiv preprint arXiv:2109 04979.
13


Reviews:
[{"review_id": "V2QWia_3B4s", "reviewer": "Reviewer_kzU4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Among the strengths we can think of the method, and the representation. Namely: \nMethods that account for relations due to confounders or other relations are important to consider and directly incorporating them in time series analysis  could avoid using regularization during model fitting. \nThe latent graph inference idea could be useful to automatically identify the relations existent in a multi-variate time series model more broadly; for instance, in partially observed settings. \nWith respect to the weaknesses, while the paper is relatively well organized, there are many details that are not clear. In particular, the architecture is presented and its components listed, but the connection of these components and the overall goal of the paper is not discussed in enough detail as to justify the choices. For instance how are colinearities identified by the architecture? How are non-linear dependencies discovered and summarized by the architecture? What other confounding can the architecture detect?", "summary_of_the_paper": "This paper presents a method to combine information of multivariate time series by extending univariate architectures. The technique uses a graph representation to represent interactions by assuming a bipartite structure, which allows the technique to scale the representation to reduce the complexity from $O(N^2)$ to $O(N K)$ using K additional nodes. The architecture represents each time series as one of N nodes and associate K embeddings (nodes) and thus there are connections between the N nodes and the K nodes but not within each group of nodes. The architecture encodes the confounders, co-linearities, and other information among the time series to improve forecasting. The experiments show the performance (and time efficiency) of the technique against several baselines on METR-LA and PEMS-BAY datasets. The experiments also show the performance on single-step forecasting on four publicly available datasets. Finally, synthetic datasets were used to evaluate the adjacency matrices in control scenarios. \n", "summary_of_the_review": "The paper has some strengths that could be useful in applications of multivariate time series forecasting. The graph representation could be useful not only for prediction but also for analysis and interpretability of results. However, I feel that more can be said about the architecture and how it affects the type of relations that it is capable of represent.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "4xd6tTfKFto", "reviewer": "Reviewer_CU8M", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is fairly well-written and its motivation is clear. It is also technically correct. The authors tried their proposed method on a few small datasets.\n\nThe FC-GNN version proposed in the submitted paper is neither novel nor significant. The main argument of the paper is reducing the computational complexity of the relational inference, which I think is a valid motivation for this paper, but I would have rather expected to try a more challenging dataset with a large number of nodes for achieving the same goal, i.e. see if the reduction in computational complexity in this way is helpful in a practical manner or not. Inferring relations between a small number of features, i.e. 8 to 800, is not a challenging problem, and other papers such as FNP (the functional neural processes, NeurIPS 2019) and BayReL (Bayesian Relational Learning, NeurIPS 2020) inferred a much larger relational graph (structured latent space) with a simple method $O(N^2)$, however, it is really challenging to learn a graph with 20,000 nodes.\n\nEven the BP-GNN model proposed in this paper is learning a Gaussian mixture in my opinion, as it will learn the connection weight of each node with Gaussian auxiliary nodes, leading to having a Gaussian mixture latent space. I believe the authors at least need to show that adding such a computational learning $O(NK)$ is improving the performance compared to a mixture of Gaussian latent space as a baseline. \n\nApart from that, the authors are learning one graph for a period of time. In my opinion, the inferred graph only captures high-level information, it might miss the short-term dependencies. An ablation study on this might be very helpful similar to those that have been done on NLP.\n\nI also believe the authors still are able to compare their model with NRI. In addition, for those methods \"Multivariate with a known graph\", the authors can construct ad-hoc graphs such as identity and/or fully connected graphs to compare the performance.  \n\nI don’t find the plots in Figure 2 very insightful and it is quite hard to see any informative insights regarding BP-GNN in my opinion. I encourage the authors to create an ablation study based on the BP-GNN and see if they can reconstruct the bipartite graph. Even for the fully connected graph, I prefer to see a real example that you have the graph and see if the method can learn the graph. The authors should be able to find some data in traffic applications for example. \n\nThe computational complexity needs to be discussed in the main paper as it is the main argument of the paper. \n\nComparisons with GRNN and RNN methods are needed in which you have the true graph or you do not have it to see if learning the graph dependency is helpful or not. \n\nThe standard deviation of the performance should be reported, considering the results are very close to each other. ", "summary_of_the_paper": "The paper proposes to infer a relational latent space in order to do multi-variate time series forecasting. The main contribution is reducing the computational complexity of the inferred graph by converting the fully connected graph learning problem to a bipartite graph learning method in which the nodes are connected to K auxiliary nodes. ", "summary_of_the_review": "Overall, I think that this paper would need to be extremely solid on the experimental side with potentially further experiments, e.g. using a larger number of features, robustness to missing data, etc, in order to make up for the somewhat limited contributions in terms of novelty.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "0AMfyQkdcVi", "reviewer": "Reviewer_xjBF", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Strengths:**\n1. The approach is interesting and has some novel elements, especially using a latent bipartite graph as a step towards better scalability. \n\nEnabling multivariate forecasting with unknown relationships among series while also addressing the scalability challenge is an important problem in forecasting, and this is a promising direction.\n\nAdditionally I really liked how the authors posed the idea as a simple add-on to univariate methods, to be able to tack on multivariate modeling to a given univariate method - though it may have been interesting to actually test this idea out more in the experiments (as in, take some other or recent state of the art global univariate models and see if using this approach adding the graph neural net to it would help - for the given datasets and especially if there are some that really depend on multivariate models).\n\n\n2. A large collection of datasets and experiments were performed, and an extensive set of graph neural net based forecasting methods compared with.\n\n\n3. Additional useful ablation study, hyper parameter sensitivity study, and adjacency matrix result analyses was also provided.  This provides more interesting information and understanding of the proposed approaches.\n\n\n**Weaknesses:**\n1.  The complete method details and procedure is not clear from the description, and code is also not provided, so as it is the work does not seem fully reproducible.\n\na.) For the auxilliary nodes, U, it says they are initialized with random noise, but nothing else is said about them - i.e., how they are used in training and inference.  In particular, is this initialization done once at the beginning of training / initializing the entire model, and the values of U are fixed from then on to those values (for the rest of training and subsequent test prediction)?\nOr are the embedding values updated as part of training as well?  Or are different random values used in each time window?  This is not clear.\n\nb.) The description is a little confusing in a few places.  In particular, in Equation 3 and later references back to it - it would be helpful to provide the actual definitions of the various phi functions (at least in the appendix) and dimensions / example dimensions of  various inputs and outputs such as the embeddings.\n\nc.) It would be helpful to understand the actual training process, even if in the appendix or said briefly - i.e., presumably sliding windows were used across the full input time series to generate each sample, presumably for every step size of 1?  Presumably 1 batch is just different time windows, containing all time series, so inputs are BxTxN where B is batch size, T is time window size, and N is the number of time series?\n\nd.) Earlier when describing past work, context c is mentioned.  How is context c incorporated into the proposed model, as it seems absent from the method description?  What about exogenous series (i.e., context or features that vary with time but are not part of the series we are predicting)?  \n\n\n2. Experiment details are missing (e.g., based on the description I assume proper hyper parameter selection was not done) and results seem not too convincing, aside from speedup using the bipartite graph formulation.\n\na.) Various missing experiment details:\ni. When describing the data, what is \"Enc. Length\" and \"Dec. Length\" in Table 2?  Is it \"Enc. Length\" the history window size used in each case, and \"Dec. Length\" the prediction window size?   If so that means for the first two datasets the next 12 values are being predicted, how then are different horizons evaluated / what is the meaning of the different horizons in the results table?  For the second set of datasets, why not evaluate predicting multiple time points (such as next 48 hours) as is commonly done in other past work, as opposed to single time point prediction?\n\nii. It would also be best to report the periodicity of the time series - i.e., is it daily, hourly, etc...\n\niii. Also details about the evaluation procedure are missing - in particular, was a sliding test window used to evaluate the performance, how many test windows were used, was re-training done before each, etc.?  Some of these may be able to be found in the pointed to prior work (saying it is the same setup) but really these details should be included in the appendix, and what was specifically done here for all particulars cannot be known from just the setup (such as if retraining is done or not).\n\niv. What context (c) or exogenous series were used in the experiments?\n\nb.) As far as I can tell, since specific and different hyper parameters are reported in the appendix for each data set, it seems like the hyper parameters were used that gave the best results on the test set, which if so can be misleading and not reflect the actual performance.  Hyper parameters should be selected from a grid (or using HPO) based on validation data only.  I saw no mention of hyper parameter selection being performed in this manner, and I did not see reported the set of hyper parameters searched over for each method, so I have to assume it was not done.  If it was done, please report the set of hyper parameters searched over for each method, data set, and horizon.  Code would have also been useful to verify this here.\n\nc.) Although 5 different runs were performed and averages across those were reported, no standard deviations are reported. The std. dev. should be reported as well, especially since in many cases there is only a tiny difference in scores compared to past methods.\n\nd.) It is hard to read through the details of the results since there are so many, but from what I can tell it does not seem like the proposed method significantly improves over past methods - results seem practically the same or only tiny improvements.  It might help to present overall average results (for example mean and std. scores across all datasets as a final column, and plotting the different average metric scores with error bars vs. dataset size / number of series or horizon), percent improvements, and especially significance analyses (e.g., a critical difference diagram I think could be most useful here given the number of different methods and dataset+horizons - e.g., see https://mirkobunse.github.io/CriticalDifferenceDiagrams.jl/dev/ or https://github.com/hfawaz/cd-diagram).\n\ne.) Comparisons to other classes of models would be nice to include as well - in particular good, recent global univariate and non-graph multivariate methods.\nIt would be nice to have representatives from the state of the art of global univariate methods (such as the recent transformer approaches such as Li et al. 2019 or more recent works that improved on this) in the comparison to see how they stack up, and validate the hypothesis that the multivariate modeling is even necessary in all these datasets.\nSimilarly it would be nice to see the comparison with recent non-graph neural net approaches, as they are just as relevant here, and there are recent approaches with state of the art results, some mentioned here, like, Salinas et al. 2019 and other more recent ones and ones focusing on scalable multivariate forecasting (mentioned below).  This would validate the hypothesis that using graphs is necessary / could help over these other approaches.\n\n\n3. Novelty is somewhat limited.\n\nAs pointed out there are several works learning graph structure as part of forecasting (and some works on hierarchical structure, though not applied to forecasting - e.g., Ying, Rex, et al. \"Hierarchical graph representation learning with differentiable pooling.\" NeurIPS 2018).  One additional missed work for learning graph structure as part of multivariate forecasting is the following:\n* Cao, Defu, et al. \"Spectral Temporal Graph Neural Network for Multivariate Time-series Forecasting.\" NeurIPS 2020.\n\nArguably the main novelty comes from introducing the latent bipartite graph in order to improve the scalability, however, the are other past works that also focus on improving scalability for multivariate forecasting, and that can learn the relationships scalably.\nWhile I feel the approach proposed here is novel given the usage of graphs, I also feel this idea is in line with other prior work that also consider somewhat similar ways of enabling scalable multivariate forecasting (e.g., using a smaller set of latent node = series), and it may be worth mentioning this connection and approaches targeted at scalable multivariate forecasting.  \n\nOne is mentioned - Salinas et al. 2019 - that similarly uses per-series model components to embed each series and combine them with a low-rank Gaussian copula, which enables more scalable modeling and training (even enabling sub-sampling among the series).  \nAnother closely related line of work targeted at scalability is to use a smaller set of latent global series (similar in concept to the latent nodes introduced here since a node corresponds to a time series) - and modeling learns these latent global series (or how to derive them from the inputs) and the input series relationship to these latent global series to get predictions for individual time series.  These similarly enable improving the complexity of the forecast modeling part from O(N^2) to O(NK^2) using the K latent series.  E.g.:\n* Sen et al. \"Think globally, act locally: a deep neural network approach to high-dimensional time series forecasting.\" NeurIPS 2019.\n* Nguyen et al. \"Temporal Latent Auto-Encoder: A Method for Probabilistic Multivariate Time Series Forecasting.\" AAAI 2021.\n\nIt would be nice to compare and contrast against these as well, at least in text to call out the related line of work, and pros and cons.  Ideally it would also be nice to see a comparison in the experiments to some of these or else another recent state of the art multivariate forecasting method that is not considering graph structure / graph neural nets.\n", "summary_of_the_paper": "The paper proposes using graph neural net (GNN) operations to combine per-series embeddings, to enable multivariate forecasting.\n\nSpecifically, N individual series are separately encoded for a given time window to get representations per series.  These representations are then updated with a GNN - either assuming a fully connected graph (with edge weights computed as part of the model), or using a bipartite graph (using a smaller set of K << N auxilliary nodes).  I.e., after multiple layers of the GNN, the representations are updated with information from the other series' representations.  Finally, the final representations are passed through per-series decoders.\n\nThe latent bipartite graph formulation enables more efficient operation of the GNN component as instead of computing O(N^2) messages only O(NK) need to be computed in a given pass.\n\nThe authors compare the proposed approach with other Graph based forecasting approaches (including ablated versions of the proposed approach) on 6 datasets - 2 where there is a given graph structure (so past work requiring the graph structure to be known can be used) and 4 where there is not.  They demonstrate competitive performance of the proposed approach (including the bipartite graph approach) and significant speed up of forward passes using the bipartite formulation especially for larger N.  They also examine the inferred adjacency matrices for different methods on some synthetic data.  Additionally, they show hyper parameter sensitivity results for varying K (number of latent nodes in the bipartite graph) for one dataset.", "summary_of_the_review": "Overall I think the idea is interesting and in particular the direction of a latent bipartite graph as a graph neural net way of making multivariate forecasting more scalable is novel, however it could be useful to point out and contrast to related alternate approaches for scalable multivariate forecasting, which are similar but not using graph neural nets.\n\nHowever there are several questions about the method and questions and concerns about the experiments that I feel need to be addressed before I can recommend the paper for acceptance - i.e., specifics of the method so it could be reproduced and understood by the readers, and details of the experiments to determine how hyper parameters were chosen and if the results show significant improvements.\n\nWithout additional clarifications and details around the experiments and the method I am leaning on the side of rejection.  However, if these can be addressed I could change my recommendation.\n\n \n\n***Update after responses and discussions:***\n\nI appreciate the responses, and read the author responses, updates, and other reviews and comments.  I tend to agree with the other reviewers' points and don't feel all concerns are addressed with the updates. \n\nAs pointed out the key novelty is really using the bipartite graph to increase computational efficiency, and as other reviewers have argued little analyses around this is done, larger data sets should be included, and training time scaling reported.  Despite the authors claim that larger data sets do not exist, and aside from simulated data, some larger ones exist and have been used in several other papers. E.g., the large wiki dataset has > 100,000 time series (and used in one of the scalable time series forecasting papers I cited and other past work it cites), and the recent M5 dataset and competition also has 30,490 inter-related time series at the lowest level, and there are several forecasting papers published about it as well.  Additionally more can be found with some searching / reading. The authors would have at least seen the wiki dataset if they looked at other prior papers including what I mentioned in my review. \n\nIn particular, as the main benefit they are claiming from their method is enabling scaling for multivariate forecasting to large amounts of time series, they ideally should consider related work on scaling forecasting for large collections of time series and at least discuss this area of work and contrast to their work, which they did not, despite the recommendation. They also should compare to the state-of-the-art multivariate forecasting methods (including without considering graphs, to demonstrate the graph approach is helpful), and especially need to compare to those designed to scale to many time series, which also was not done.\n\nI also agree with the other reviewer - to only compare and report forward evaluation run times leaves uncertainty about practical performance, as training is the biggest issue in most cases.  Further, there seems little reason not to report it and explain what it reveals unless the benefits are not as much in this case - which would also be interesting and fine just needs to be reported and some analyses provided. Ideally the analyses around training and forward time should be performed for varying number of time series across multiple datasets, e.g., plotting on x-axis the number of time series, and on y-axis the training and forward time (and nice to also see the error metrics as well per number of time series).  It could even be combined in the same tables or plots for forward time by showing average training time per step and average forward time per step, for example, on dual y-axis if space is a concern) vs. number of time series.   This could potentially also be generated from larger datasets if needed by sub-sampling series, or just using the natural number of time series in different datasets.\n\nThe fact training time isn't reported makes one wonder if there could be some reasons it's not providing as much scaling benefits for training, e.g., maybe the train time is actually not much if any better, due to the extra propagation steps having more of an impact at training time as perhaps slowing down the back propagation updates (which we don't see from forward time alone), and also from the fact that all time series still must be present at once so there could be some dominating factor on the data side that doesn't change.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies multivariate time series forecasting by making relational inference in a latent space. It attempts to address the important issue of reducing the computational complexity of the inferred graph. This motivation is well articulated.

Despite its merits, concerns have been raised regarding the relatively weak evaluation without using datasets involving more many nodes to demonstrate the scalability of the proposed method, which is a major selling point of the paper. As such, while the motivation of the work is clear, its experimental evaluation is not thorough enough to demonstrate the scalability of the proposed method.

The authors made the remark in their response that they are not aware of any public time series dataset of this size (which is not agreed by another reviewer who pointed out that some much larger datasets were used in other papers). Note that it is not uncommon in other work to use synthetic datasets to evaluate the scalability as well as other properties of the proposed methods.

Moreover, clarity of the presentation also has room for improvement.

The paper has potential for publication in a top venue if the comments and suggestions are incorporated to revise the paper.
Decision: Reject