Discrete Graph Structure Learning for Forecasting Multiple Time Series
Authors: Chao Shang
Jie Chen
Jinbo Bi
Venue: ICLR 2021 Poster (2021)
Keywords: Time series forecasting
graph neural network
graph structure learning

Abstract: Time series forecasting is an extensively studied subject in statistics, economics, and computer science. Exploration of the correlation and causation among the variables in a multivariate time series shows promise in enhancing the performance of a time series model. When using deep neural networks as forecasting models, we hypothesize that exploiting the pairwise information among multiple (multivariate) time series also improves their forecast. If an explicit graph structure is known, graph neural networks (GNNs) have been demonstrated as powerful tools to exploit the structure. In this work, we propose learning the structure simultaneously with the GNN if the graph is unknown. We cast the problem as learning a probabilistic graph model through optimizing the mean performance over the graph distribution. The distribution is parameterized by a neural network so that discrete graphs can be sampled differentiably through reparameterization. Empirical evaluations show that our method is simpler, more efficient, and better performing than a recently proposed bilevel learning approach for graph structure learning, as well as a broad array of forecasting models, either deep or non-deep learning based, and graph or non-graph based.

Published as a conference paper at ICLR 2021
DISCRETE GRAPH STRUCTURE LEARNING FOR FORE- CASTING MULTIPLE TIME SERIES
Chao Shang"
Jie Chen
University of Connecticut chao. wwwngguconnnnnnn..oo edu
MIT-IBM Watson AT Lab, IBM Research chenjie bm. com
Jinbo Bi
University of Connecticut jinbo. www.ucin.eeeec...... edu
ABSTRACT
Time series forecasting is extensively studied subject in statisties, economics, and computer science. Exploration of the correlation and causation among the variables in multivariate time series shows promise in enhancing the perfor- mance of a time series model. When using deep neural networks as forecasting models, we hypothesize that exploiting the pairwise information among multiple (multivariate) time series also improves their forecast. If an explicit graph struc- ture is known, graph neural networks (GNNs) have been demonstrated powerful tools to exploit the structure. In this work, we propose learning the structure si- multaneously with the GNN if the graph is unknown. We cast the problem learning probabilistic graph model through optimizing the mean performance over the graph distribution. The distribution is parameterized by a neural network so that discrete graphs can be sampled differentiably through reparameterization. Empirical evaluations show that our method is simpler, more efficient, and better performing than recently proposed bilevel learning approach for graph structure learning, as well as a broad array forecasting models, either deep or non-deep learning based, and graph or non-graph based.
INTRODUCTION
Time series data are widely studied in science and engineering that involve temporal measurements. Time series forecasting is concerned with the prediction of future values based on observed ones in the past. It has played important roles in climate studies, market analysis, traffic control, and energy grid management (Makridakis et al.) 1997) and has inspired the development of various predictive models that capture the temporal dynamics of the underlying system. These models range from early autoregressive approaches (Hamilton, 1994 Asteriou & Hall 2011 to the recent deep learning methods (Seo et al. 2016 Li et al.] 2018 et al. 2018 Zhao al. 2019).
Analysis of univariate time series (a single longitudinal variable) has been extended to multivariate time series and multiple (univariate or multivariate) time series. Multivariate forecasting models find strong predictive power in stressing the interdependency (and even causal relationship) among the variables. The vector autoregressive model (Hamilton) [994] is an example of multivariate analysis, wherein the coefficient magnitudes offer hints into the Granger causality (Granger 1969) of one variable to another.
For multiple time series, pairwise similarities connections among them have also been explored to improve the forecasting accuracy (Yu et al. 2018). An example is the traffic network where each node denotes time series recording captured by particular sensor. The spatial connections of the roads offer insights into how traffic dynamics propagates along the network. Several graph neural
"This work was done while C. Shang was an intern at MIT-IBM Watson AI Lab, IBM Research. To whom correspondence should be addressed.

---

Published as a conference paper at ICLR 2021
network (GNN) approaches (Seo et al.
2016 Li et al. 2018 Yu et al. 2018 Zhao al. 2019) have been proposed recently to leverage the graph structure for forecasting all time series simultaneously.
The graph structure however is not always available it may be incomplete. There could be sev- cral reasons, including the difficulty in obtaining such information or a deliberate shielding for the protection of sensitive information. For example, data set comprising sensory readings of the nation-wide energy grid is granted access to specific users without disclosure of the grid structure. Such practical situations incentivize the automatic learning of the hidden graph structure jointly with the forecasting model.
Because GNN approaches show promise in forecasting multiple interrelated time series, in this pa- per we are concerned with structure learning methods applicable to the downstream use of GNNs. A prominent example is the recent work of |Franceschi et al.) (2019) (named LDS), which a meta-learning approach that treats the graph as a hyperparameter in a bilevel optimization frame- work (Franceschi et al.] (2017). Specifically, let Xumin and N,u denote the training and the validation sets of time series respectively, A E {0, 1 _ denote the graph adjacency matrix of the time se- ries, w denote the parameters used in the GNN, and and F denote the the loss functions used dur- ing training and validation respectively (which may not be identical). LDS formulates the problem as learning the probability matrix 0 E [0, 1|mxn, which parameterizes the element wise Bernoulli distribution from which the adjacency matrix A is sampled:
min (x~Ber(0) (A,A, (f). Xva)),
s.t. w(0) argmin Eew)Be(0)) [LAA, w, Xmmn)]
Formulation gives bilevel optimization problem. The constraint (which by itself is an optimiza- tion problem) defines the GNN weights as a function of the given graph, so that the objective is to optimize such graph only. Note that for differentiabilityy one does not directly operate on the discrete graph adjacency matrix A, but on the continuous probabilities instead.
LDS has two drawbacks. First, its computation is expensive. The derivative of w with respect 0 is computed by applying the chain rule on a recursive-dynamics surrogate of the inner optimization argmin. Applying the chain rule on this surrogate equivalent to differentiating an RNN, which is either memory intensive if done in the reverse mode or time consuming if done in the forward mode, when unrolling deep dynamics. Second, it is challenging to scale. The matrix 0 has (n²) entries to optimize and thus the method is hard to scale to increasingly more time series.
In light of the challenges of LDS, we instead advocate a unilevel optimization:
min W
(AA-Ber((((w)) (F(A, w, Xumin)).
Formulation 2 trains the GNN model as usual, except that the probabilities & (which parameterizes the distribution from which A is sampled), is by itself parameterized. We absorb these parameters, together with the GNN parameters, into the notation w. We still use a validation set X for usual hyperparameter tuning, but these hyperparameters are not 0 as treated by (T. In fact, formulation (T) may need second validation set to tune other hyperparameters.
The major distinction of our approach from LDS is the parameterization ((w). as opposed to an inner optimization w().. In our approach, a modeler owns the freedom to design the parameterization and better control the number of parameters as n² increases. To this end, time series representation learning and link prediction techniques offer ample inspiration for modeling. In contrast, LDS is more agnostic as no modeling needed. The effort, instead, lies in the nontrivial treatment of the inner optimization (in particular, its differentiation).
As such, our approach is advantageous in two regards. First, its computation less expensive, because the gradient computation of a unilevel optimization is straightforward and efficient and implementations are mature. Second, it better seales, because the number of parameters does not grow quadratically with the number of time series.
We coin our approach GTS (short for "graph for time series"), signaling the usefulness of graph structure learning for enhancing time series forecasting. It is important to note that the end purpose of the graph is to improve forecasting quality, rather than identifying causal relationship of the series or recovering the ground-truth graph, if any. While causal discovery of multiple scalar variables is an
2

---

Published as a conference paper at ICLR 2021
established field, identifying causality among multiple multivariate time series requires nontrivial extension that spans beyond the current study. On the other hand, the graph, either learned or pre- existing, serves additional information that helps the model better capture global signals and apply on each series. There does not exist golden measure for the quality of the learned graph except forecasting accuracy. For example, the traffic network does not necessarily offer the best pairwise relationship GNN can exploit for forecasting traffic series. Nevertheless, to robustify GTS we incorporate regularization that penalizes significant departure from one's prior belief. If a certain "ground-truth graph is believed, the learned graph will be healthy variation of it for more accurate forecast.
2 RELATED WORK
Time series forecasting has been studied for decades by statisticians. Is out of the scope of this pa- per to comprehensively survey the literature, but we will focus more on late developments under the deep learning context. Early textbook methods include (vector) autoregressive models (Hamilton) 1994), autoregressive integrated moving average (ARIMA) (Asteriou & models (HMM) (Baum & Petrie 1966), and Kalman filters (Zarchan Musofi 20000. Generally speaking, these are linear models that use window of the past information to predict the next time step, although nonlinear versions with parameterization are subsequently developed.
hidden Markov
A notable nonlinear extension was the RNN (Williams et al., 1986), which later evolved into LSTM (Hochreiter Schmidhuber 1997). BiLSTM (Schuster & Paliwal 1997). and GRU (Cho et al. 2014). which addressed several limitations of the vanilla RNN, such as the vanishing gradient problem. These architectures are hard to parallelize because of the recurrent nature of the forward and backward computation. More recently, Transformer (Vaswani et al. (017) and BERT (Devlin) jet al. 2019) were developed to address parallelization, by introducing attention mechanisms that si- multancously digested past (and future) information. Although these models are more heavily used for sequence data under the context of natural language processing, they are readily applicable for time series well (Shih et al. 2019 Li
al. 2019).
Graph neural networks (Zhang et al., 2018 Zhou et al.) 2018 Wu et al. 2019) emerged quickly in deep learning handle graph-structured data. Typically, graph nodes are represented by feature vec- tors, but for the case of time series, a number of specialized architectures were recently developed; see, e.g., GCRN (Seo al., 2016), DCRNN (Li et al. 2018). STGCN (Yu et al., 2018), and T- GCN (Zhao et al. 2019). These architectures essentially combine the temporal recurrent processing with graph convolution augment the representation learning of the individual time series.
Graph structure learning (not necessarily for time series) appears in various contexts and thus meth- ods span broad spectrum. One field of study is probabilistic graphical models and casual infer- ence, whereby the directed acyclic structure is enforced. Gradient-based approaches in this context include NOTEARS (Zheng et al., 2018). DAG-GNN (Yu et al., 2019), and GraN-DAG (Lachapelle)
et al. 2020). On the other hand, general graph may still be useful without resorting to causality. LDS Franceschi et al. 2019) is a meta-learning approach that demonstrates to improve the perfor- mance on node classification tasks. MTGNN (Wu et al. (2200) parameterizes the graph as degree-k graph, which is learned end-to-end with a GNN for forecasting time series. We, on the other hand, vatent-variable approach and learns a latent graph for forecasting system dynamics. Our approach is closely related to NRI and will compare with it in the following section after introducing the technical details.
allow more general structural prior for the graph. NRI (Kipf et al. 2018) adopts
3 METHOD
In this section, we present the proposed GTS method, elaborate the model parameterization, and describe the training technique. We also highlight the distinctions from NRI (Kipf et al.) 20018).
Let us first settle the notations. Denote by X the training data, which is three dimensional tensor, with the three dimensions being feature, time, and the n series. Superscript refers to the series and subscript refers to time; that is, X denotes the i-th series for all features and time and X; denotes the f-th time step for all features and series. There are in total S time steps for training. The model will window of T steps forecast the next 7 steps. For each valid denote by
3

---

Published as a conference paper at ICLR 2021
entire data X
learned structure 0
sampled graph
feature extractor
link predictor
sampling
recurrent graph convolution
recurrent graph convolution
recurrent graph convolution
recurrent graph convolution
t+
T+
windowed data
forecast
Figure 1: GTS architecture.
X44T+144T- = [(A, w, X the model, which forecasts X44T+14+4++ from observations X++1++T, through exploiting the graph structure and being parameterized by w. Using to denote the loss function between the prediction and the ground truth, a typical training objective reads
E. w, X(41+4)), X++T+1+44+7).
Three remaining details are the parameterization of A, the model f, and the loss
3.1 GRAPH STRUCTURE PARAMETERIZATION
The binary matrix A 10, by itself is challenging parameterize, because requires differentiable function that outputs discrete values 0/1. A natural idea is to let A be random variable of the matrix Bernoulli distribution parameterized by 0 E [o, _]"x", that Aij independent for all the (i, j) pairs with A. ~ Ber(0:5). Here, fij is the success probability of Bernoulli distribution. Then, the training objective needs be modified to
E EAwBer(0) Σ> (((A, w, X(11+T), X [
(4)
As hinted in Section we further parameterize 0 as 0(w). because otherwise the n² degrees of freedom in 8 render the optimization hard to scale. Such parameterization, however, im- poses challenge on differentiability. if the expectation (4) is evaluated through sample average: the gradient of (4) does not flow through A in usual Bernoulli sampling. Hence, we apply the Gumbel reparameterization trick proposed by Jang et al. (2017) and Maddison et al. (2017): Ay sigmoid (log (0ty/ 04)) - (g); L )ss), where Iip Bi Gumbel 1) for all i, When the temperature + 0, A. = 1 with probability 0; and 0 with remaining probability. In practice, we anneal progressively training such that tends to zero.
For the parameterization of 8, we use a feature extractor to yield feature vector for each series and link predictor that takes in pair of feature vectors and outputs link probability. The feature extractor maps matrix X to vector 2 for each i. Many sequence architectures can be applied; we opt for simple one. Specifically, we perform convolution along the temporal dimension, vectorize along this dimension, and apply fully connected layer to reduce the dimension: that is, 2 (C(vec(Conv(("))
Note that the feature extractor is conducted on the entire sequence rather than window of T time steps. Weights shared among all series.
The link predictor maps pair of vectors (2f,2)) 20) to a scalar 0.y E [0, We concatenate the two vectors and apply two fully connected layers achieve so; that is, 0g FC(FC(=*||=/)) The last activation needs be sigmoid. See the top part of Figure
3.2 GRAPH NEURAL NETWORK FORECASTING
The bottom part of Figure I is the forecasting model f. We use a souuunccet---eeueneee (seq2seq) model (Sutskever al. 2014) map X (4414 X X for each series Seq2seq is
4

---

Published as a conference paper at ICLR 2021
typically recurrent model, but with graph structure available among the series, leverage recurrent graph convolution to handle all series simultaneously, as opposed to the usual recurrent mechanism that treats each series separately.
Specifically, for each time step t', the seq2seq model takes X² for all series as input and updates the internal hidden state from - to Hy. The encoder part of the seq2seq performs recurrent updates from = t + to = t + T, producing H²+T as a summary of the input. The decoder part uses Hz+T to continue the recurrence and evolves the hidden state for another steps. Each hidden state
Hy, 1' = t T + 1 1 + T + т, simultaneously serves as the output X and the input to the next time step.
The recurrence that accepts input and updates hidden states collectively for all series uses a graph convolution to replace the usual multiplication with a weight matrix. Several existing architectures serve this purpose (e.g.. GCRN (Seo et al. 2016), STGCN (Yu al. 20018), and T-GCN (Zhao al.) 2019h, but is designed for directed graphs:
the diffusion convolutional GRU defined in DCRNN (Li al. 2018) because
Re = sigmoid (WR *A [Xv || Hy-1] | + bR),
C = tanh A [X (Ry @ He-il+ bc), Hy Up Hr-1 + (1 Uv) Cv.
Up sigmoid (Wo *A [Xy || Hy-1]+ bu).
where the graph convolution * defined as
We x Y - (w,, A)* + 2 A2))
with Do and Dy being the out-degree and in-degree matrix and || being concatenation along the feature dimension. Here, 2 22 be for Q = R, U, C are model parameters and the diffusion degree K a hyperparameter.
We remark that subsequent experiment corroborates, this GNN model can be replaced by other similar ones (e.g.. T-GCN), such that the forecast performance remains similar while still being superior over all baselines. In comparison, the more crucial part of our proposal is the structure learning component (presented the preceding subsection), without which falls back to a model either using no graphs or needing a supplied one, both performing less well.
3.3 TRAINING, OPTIONALLY WITH A PRIORT KNOWLEDGE OF THE GRAPH
The base training loss (per window) is the mean absolute error between the forecast and the ground truth
(+T base ( (R4T+IITHH, X4T+L14T4+) L [& - -el-
Additionally, we propose regularization that improves graph quality, through injecting priori knowledge of the pairwise interaction into the model. Sometimes an actual graph among the time series is known, such as the case of traffic network mentioned in Section [ Generally, even an explicit structure is unknown, neighborhood graph (such a ANN graph) may still serve as reason- able knowledge. The use ANN encourages sparsity k: is small, which circumvents the drawback of & constraints that cannot be easily imposed because the graph is not As such, we use the cross-entropp between and the a priori graph A" as the regularization:
variable to optimize.
freg -Eg AAi log Osj - (1 --yy) log(1 0zy).
(5)
The overall training loss is then E, Pamce + Alegg with X > 0 being the regularization magnitude.
3.4 COMPARISON WITH NRI
GTS appears similar to NRI (Kipf et al., 2018) on the surface, because both compute a pairwise structure from multiple time series and use the structure to improve forecasting, In these two meth- ods, the architecture to compute the structure, as well as the one to forecast, bare many differences; but these differences are only secondary. The most essential distinction is the number of structures. To avoid confusion, here we say "structure" (0) rather than "graph" (A) because there are combina- torially many graph samples from the same structure. Our approach produces given one set of n series. On the contrary, the autoencoder approach adopted by NRI produces dif- ferent structures given different encoding inputs. Hence, feasible use of NRI can only oceur in the
single structure
5

---

Published as a conference paper at ICLR 2021
following two manners. (a) A single set of n series is given and training is done windowed data, where each window will produce a separate structure. (b) Many sets are given and training is done through iterating each set, which corresponds to separate structure. Both cases different from our scenario, where a single set of time series is given and a single structure is produced.
Fundamentally, NRI is variational autoencoder and thus the inference of the structure is amor- tized inference: under setting (b) above, the inferred structure a posterior given a set of series. The amortization uses an encoder parameterization to free off the tedious posterior inference whenever new set of series arrives. Moreover, under the evidence lower bound (ELBO) training objective, the prior graph, each edge of which takes value uniformly in [0. 1]. our case, on the contrary, single structure is desired. Thus, amortized inference is neither necessary nor relevant. Furthermore, one may interpret the priori information A for regularization as a "structural prior;" however, for each node pair it offers a stronger preference on the existence/absence of an edge than uniform probability.
4 EXPERIMENTS
In this section, conduct extensive experiments to show that the proposed method GTS outper- forms a comprehensive set of forecasting methods, including one that learns a hidden graph structure (LDS, adapted for time series). We also demonstrate that GTS is computationally efficient and is able to learn a graph close to the a priori knowledge through regularization, with little compromise on the forecasting quality.
4.1 SETUP
Data sets. We experiment with two benchmark data sets METR-LA and PEMS-BAY from Li et al. (018)) and a proprietary data set PMU. The first two are traffic data sets with given graphs serving as ground truths; we perform no processing and follow the same configuration as in the referenced work for experimentation. The last one is a sensor network of the U.S. power grid without a given grid topology. For details, see Appendix Section A For all data sets, we use a temporal 70/10/20 split for training, validation, and testing, respectively.
Baselines. We compare with a number of forecasting methods:
1. Non-deep learning methods: historical average (HA), ARIMA with Kalman filter (ARIMA), tor auto-regression (VAR), and support vector regression (SVR). The historical average accounts for weekly seasonality and predicts for day by using the weighted average of the same day in the past few weeks.
Deep learning methods that treat each series separately (i.e., no graph): feed-forward neural network (FNN) and LSTM.
3. GNN method applied the given graph (or ANN graph for PMU): DCRNN (Li et al. 2018
4. GNN methods that simultaneously learn graph structure. We use LDS (Franceschi et al. 2019) learn the graph, wherein the forecast model is DCRNN. We name the method "LDS" for short. Additionally, we compare with NRI (Kipf al. 2018).
5. Variant of GTS: We maintain the graph structure parameterization but replace the DCRNN fore- cast model by T-GCN (Zhao et al. 2019). We name the variant "GTSv."
Except LDS and NRI, all baselines follow the configurations presented in L Li 2018). For LDS, we follow Franceschi et al. (2019). For NRI, we follow Kipf et al. (0018.
Evaluation metrics. All methods evaluated with three metries: root mean square emror (RMSE), and mean absolute percentage error (MAPB).
absolute error (MAE),
For details hyperparameter setting and training platform, Appendix Section B

---

Published as a conference paper at ICLR 2021
4.2 RESULTS
Forecasting quality. We first evaluate the performance GTS through comparing with all the aforementioned baselines. Because of the significant memory consumption of NRI, this method is exccuted on only the smaller data set PMU. The tasks are to forecast 15, 30, and 60 minutes.
Table 1: Forecasting error (METR-LA).
Metric
HA
ARIMA
VAR
SVR
FNN
LSTM
DCRNN
LDS
GTSv
GTS
min mmh
MAE
4.16
3,99
4.42
3.99
3.44
2.77
2.61
2.39
RMSE
7.80
8.21
7.89
8.45
7.94
6,30
5.38
5.35
4.83
4.41
МАPE
13.0%
9.6%
10.2%
9.3%
9.9%
9.6%
7.3%
7.1%
6.8%
6.0%
MAE
4.16
5.15
5.41
5.05
4.23
3.77
3.15
3.14
2.94
2.65
RMSE
7.80
10.45
9.13
8.17
7.23
6.45
5.06
МАРE
13.0%
12.7%
12.7%
12.1%
12.9%
10.9%
8.8%
8.6%
7.9%
7.0%
min
MAE
4.16
6.90
6.52
6.72
4.49
4.37
3.60
3,63
3,46
2.99
RMSE
7.80
13.23
10.11
13.76
8.69
8.69
7.59
7.67
6.67
5.85
MAPE
13.0%
17.4%
15.8%
16.7%
14.0%
13.2%
10.5%
10.34%
9.9%
8.3%
Table Jummarries the results for METR-LA. A few observations follow. (1) Deep learning meth- ods generally outperform non-deep learning methods, except historical average that performs on par with deep learning in some metrics. Scasonality is strong indicator of the repeating traffic patterss and not surprisingly HA performs reasonably well despite simplicity. (2) Among the deep learning methods, graph- based models outperform non-graph models. This result corroborates the premise of this work: graph structure is helpful. (3) Among the graph-based methods, LDS performs slightly better than DCRNN. The difference between these two methods is that the latter employs the given graph, which may may not imply direct interactions, whereas the former learns graph in the data-driven manner. Their performances however are quite similar. (4) The most encouraging result is that the proposed method GTS significantly outperforms LDS and hence DCRNN, GTS learns graph structure through parameterization, rather than treating as phrameerr which the case in LDS. (5) The performance of the variant GTSv stays between GTS and LDS. This observa- tion corroborates that the proposed structure learning component contributes more crucially to the overall performance than does a careful choice of the GNN forecasting component.
To dive into the behavioral difference between GTS and DCRNN, we plot Figure forecasting examples. One
that both methods produce smooth series. In the top example, overall the GTS curve is closer to the moving average of the ground truth than is the DCRNN curve (see c.g., the left part and the U shape). In the bottom example, the GTS curve better captures the sharp dip toward the end of the series. both examples, there exist several short but deep downward spikes. Such anomalous data captured by neither methods.
(Seconds)
1600
1562
LOS-DCRNN
milehh
Our Model
1230
1000
Secane
00:06
03:00
06:00
09:00
12:00
15:00
18:00
21:00
00:00
800
03-05-2012
600
Training
200
163
mile/n
0
METR-LA
PEMS-BAY
Grous Tun
Figure 3: Training time per epoch. Not learn- ing graph (DCRNN) is the fastest to train and learning graph by using LDS needs orders of magnitude more time. Our model learns graph with favorable overhead.
00:00
03:00
00
00
12:00
15:00
18:00
21.00
00:00
03-05-2012
Figure 2: One-day forecast (METR-LA)
7

---

Published as conference paper at ICLR 2021
Additionally, we summarize the results for (see Appendix Section C for the former). The observations are rather similar to those of METR- LA. Our model produces the best prediction in all scenarios and under all metrics. Additionally, for the PMU data set, NRI performs competitively, second to GTS/GTSv and better than LDS in most of the cases.
and PMU in Table and 2 respectively.
Table 2: Forecasting error (PMU).
Metric
FNN
LSTM
DCRNN
LDS
NRI
GTSv
GTS
MAE MAE(x10
1.23
1.02
0.71
0.49
0.66
0.26
0.24
RMSE (x10-3
1.28
1.63
1.42
1.26
0.27
0.20
0.19
15
MAPE
0.20%
0.21%
0.09%
0.07%
0.14%
0.05%
0.04%
MAE (x10
1.42
1.11
1.08
0.81
0.71
0.31
0.30
RMSE (x10
2.06
1.91
1.79
0.30
0.23
60 30
MAPE
0.23%
0.20%
0.15%
0.12%
0.15%
0.05%
0.05 %
MAE (x
1,88
1.79
1.78
1.45
0.83
0.39
0.41
RMSE (x10-3
2.58
2.75
2.65
2.54
0.46
0.32
0,30
MAPE
0.29%
0.27%
0.24%
0.22%
0.17%
0.07 %
0.07 %
Computational efficiency. We compare the training costs of the graph-based methods: DCRNN, LDS, and GTS. See Figure DCRNN is the most efficient to train, since no graph structure learn- ing involved. To learn the graph, LDS needs orders of magnitude more time than does DCRNN. Recall that LDS employs a bilevel optimization ( which is computationally highly challenging. In contrast, the proposed method GTS learns the graph structure as byproduct of the model train- ing 2 Its training time is approximately three times of that of DCRNN, favorable overhead compared with the forbidding cost of L.DS.
Effeet of regularization. We propose in Section 3.3 using regularization to incorporate priori knowledge of the graph. One salient example knowledge is sparsity, which postulates that many node pairs barely interact. We show the effect of regularization on the data set PMU with the use of a &NN graph as knowledge. The task is 15-minute forecasting and results (expected degree and MAE) are plotted in Figure 4 The decreasing curves in both plots indicate that using smaller k or increasing the regularization magnitude produces sparser graphs. The bars give the MAEs, all around 2.4e-4, indicating equally good forecasting quality. (Note that the MAE for LDS is 4.9e-4.)
20
Degree
Degree
40
MAE
MAE
Degree
MAE (10-
Degree
MAE (10-
10
20
3
10
1
k=0
k=30
1=10
k=5
A=
A=2
A=10
A=20
Fix regularization magnitude and vary k.
(b) Fix and vary regularization
The case means no &NN regularization.
Figure Effect of regularization (PMU), "Degree" means expected degree of the graph.
Learned structures. To examine the learned structure 0, we further show its difference from the given graph adjacency matrix A" (binary) and visualize one particular example in Figure 5 The difference is defined as freg/ freg/2 (average cross entropy; see (5). One reads that when X = 20, the difference 0.34. indicates that the learned probabilities in & are on average 0.3 away from the cutoff threshold for 0 such
entries of A", because log(1 0.3) 0.34. When using 0.5 difference possibly results in false-positive edges (existing in 9 but not in A"; orange dashed) and false-negative edges (existing in A but not in 0: none in the example).
8

---

Published as a conference paper at ICLR 2021
Note that the regularization strength À weighs the forecasting error (MAE) and the cross entropy in the loss function. When À = O, the training loss is not regularized, yielding optimal forecast results reported in Table 2 When À oc, effectively enforces 0 to be identical AA and hence the model reduces to DCRNN, whose forecasting performance is worse than our model. The interesting question is when À interpolates between the two extremes, whether it is possible to find sweet spot such that forecasting performance is close to our model but meanwhile 0 is close to A". Figure 5 suggests positively. We stress that our model does intend to learn a "ground-truti" graph (e.g., the traffic network or the power grid); but rather, learn a structure that GNN can exploit to improve forecast.
A =
À x=
À
10
À
20
CE
1.93
1.87
0.53
0).34
15 min
MAE
2,47e-4
2.74c-4
2.83e-4
2.87c-4
CE
1.93
1.87
0.53
0..4
30 min
MAE
3.02e-4
3.26e-4
3,44e-4
3,59e-4
CE 1.93
1.87
0.53
0.34
60 min
MAE
4.14e-4
4.33e-4
4.78e-4
5.12e-4
Figure 5: Learned structures (PMU). Left: difference between A and 0 (in terms of average cross entropy) as À varies. Right: A overlaid with 0 (orange edges exist in & but not in A") when À = 20.
Other structural priors. In the PMU data set, we use synthetic ANN structure prior A due to the lack of known graph. For METR-LA and PEMS-BAY, however, such as graph can be con- structed based on spatial proximity (Li al. 2018). We show in Figure effect of regularization for these data sets. Similar to the findings of PMU, moving between 0 and oc interpolates two extremes: the best forecast quality and recovery of A". With reasonable choice of X (e.g., 0.3), the forecast quality degrades only slightly but the learned structure rather close to the given A3, judged from the average cross entropy.
10
3.50
10a
Cross Entropy
CE
Cross Entropy
CE
1.7
MAE
3.25
MAE
100
100
MAE
1.6
3.00
MA
10-1
10-1
2.75
1.5
10-2
2.50
10-2
1.4
0.3
A
A
Figure 6: Effect of regularization (left: METR-LA; right: PEMS- BAY; 60 min forecast). Average cross entropy measures departure from A".
5 CONCLUSIONS
We have presented a time series forecasting model that learns graph structure among multiple time series and forecasts them simultancously with a GNN. Both the graph and the GNN are learned end-to-end, maximally exploiting the pairwise interactions among data streams. The graph struc- ture is parameterized by neural networks rather than being treated as (hyper)parameter, hence significantly reducing the training cost compared with a recently proposed bilevel optimization ap- proach LDS. We conduct comprehensive comparisons with number of baselines, including non- deep learning methods and deep learning methods (which either ignore the pairwise interaction, use given graph, or learn graph by using LDS), and show that our approach attains the best forecast- ing quality. We also demonstrate that regularization helps incorporate a priori knowledge, rendering the learned graph healthy variation of the given for accurate forecast.

---

Published as a conference paper at ICLR 2021
ACKNOWLEDGMENT AND DISCLAIMER
This material is based upon work supported by the Department of Energy under Award Number(s) DE-0E0000910. C. Shang was also supported by National Science Foundation grant IIS-1718738 (to J. Bi) during this work. J. Bi was additionally supported by National Institutes of Health grants K02-DA043063 and R01-DA051922. This report was prepared as an account of work sponsored by agencies of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal li- ability or responsibility for the accuracy, completeness, usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommenda- tion, favoring by the United States Government any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.
REFERENCES
Dimitros Asteriou and Stephen G. Hall. ARIMA models and the Box-Jenkins methodology. In Applied Econometrics, pp. 265-286. Palgrave MacMillan, 2011.
L. E. Baum and T. Petrie. Statistical inference for probabilistic functions of finite state markov chains. The Annals of Mathematical Statistics, 37(6):5544-1565 1966.
Kyunghyun Cho, Bart van Merrienboera, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Hol- ger Schwenk, and Yoshua Bengio. Leamning phrase representations using RNN encoder-decoder for statistical machine translation. In EMNLP. 2014.
Jacob Devlin, Ming- Wei Chang, Kenton Lee, and Kristina Toutanova. BERT: Pre-training of deep bidirectional transformers for language understanding. In NAACL-HLT, 2019.
Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and reverse gradient-based hyperparameter optimization. In ICML, 2017.
Luca Franceschi, Mathias Niepert, Massimiliano Pontil, and Xiao He. Learning discrete structures for graph neural networks. In JCML, 2019.
C. W. J. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3):4244333 1969.
James D. Hamilton. Time Series Analysis. Princeton University Press, 1994.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, 9(8): 1735-1780, 1997.
H. V. Jagadish, Johannes Gehrke, Alexandros Labrinidis, Yannis Papakonstantinou, Jignesh M. Pa- tel, Raghu Ramakrishnan, and Cyrus Shahabi. Big data and its technical challenges. Commun. ACM, 57(7):86-99 94, 2014.
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with Gumbel-Softmax. In ICLR, 2017.
Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In ICML, 2018.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural DAG learning. In ICLR, 2020.
Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan. Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting. In NewrlPS, 2019.
10

---

Published as a conference paper at ICLR 2021
Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural net- work: Data-driven traffic forecasting. In ICLR, 2018.
Chris J. Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. In ICLR, 2017.
Spyros G. Makridakis, Steven C. Wheelwright, and Rob Hyndman. Forecasting: Methods and Applications. Wiley, 3 edition edition, 1997.
Mike Schuster and Kuldip K. Paliwal. Bidirectional recurrent neural networks. IEEE Transactions on Signal Processing, 45(11):2673-266 1997.
Youngjoo Seo, Michaël Defferrard, Pierre Vandergheynst, and Xavier Bresson. Structured sequence modeling with graph convolutional recurrent networks. arXiv:161 2016.
Shun- Yao Shih, Fan- Keng Sun, and Hung yi Lee. Temporal pattern attention for multivariate time series forecasting. Machine Learning volume, 108:1421-1441, 2019.
Ilya Sutskever, Oriol Vinyals, and Quoc Viet Le. Sequence to sequence learning with neural net- works. In NIPS, 2014.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In NIPS, 2017.
Ronald J. Williams, Geoffrey E. Hinton, and David E. Rumelhart. Learning representations by Back-propagating errors. Nature. 323(6088):533-536 1986.
Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive survey on graph neural networks. Preprint arXiv: 1ee::191..000 900566,,, 2019.
Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang. Con- necting the dots: Multivariate time series forecasting with graph neural networks. In KDD, 2020.
Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. In IJCAI, 2018.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In ICMI., 2019.
Paul Zarchan and Howard Musoff. Fundamentals of Kalman Filtering. American Institute of Aero- nautics and Astronautics, Incorporated, 2000.
Ziwei Zhang, Peng Cui, and Wenwu Zhu. Deep learning on graphs: A survey. Preprint 1nXiv:181..04200 2018.
Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li. T-GCN: temporal graph convolutional network for traffic prediction. IEEE Transactions on Intelligent Transportation Systems, 2019.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Con- tinuous optimization for structure learning. NeurlPS, 2018.
Jie Zhou, Ganqu Cui, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. Preprint 1tXiv:8812008434. 2018.
11

---

Published as a conference paper at ICLR 2021
Appendices
A ADDITIONAL DETAILS OF DATA SETS
METR-LA is a traffic data set collected from loop detectors in the highway of Los Angles, CA (Ja- gadish et al., 2014). contains 207 sensors, each of which records four months of data at the frequency of five minutes. A graph of sensors is given; it was constructed by imposing radial basis function on the pairwise distance of sensors at a certain cutoff. For more information see et al. 2018). We perform no processing and follow the same configuration as in Lt al. (2018)
PEMS-BAY is also a traffic data set, collected by the California Transportation Agencies Perfor- mance Measurement System. It includes 325 sensors in the Bay Area for a period of six months, at same five-minute frequency. Construction of the graph the same as that of METR-LA. No processing is performed.
PMU contains time series data recorded by the phasor measurement units (PMUs) deployed across the U.S. power grid. We extract one month of data (February 2017) from one interconnect of the grid, which includes 42 PMU sources. Each PMU records number of state variables and we use only the voltage magnitude and the current magnitude. The PMUs sample the system states at high rates (either 30 or 60 Hertz). We aggregate every five minutes, yielding a data frequency the same the above two data sets. Different from them, this data set offers neither the grid topology, the sensor identities, nor the sensor locations. Hence, a "ground truth" graph is unknown.
However, it is highly plausible that the PMUs interact in nontrivial manner, since some series highly correlated whereas others not much. Figure 7 shows three example series. Visually, the first series appears more correlated to the second one than to the third one. For example, in the first two series, the blue curves (the variable ip_m) are visually seasonal and synchronous. Moreover, inside the purple window, the red curves (the variable vp.m) the first two series show three down- ward spikes, which are missing in the third series. Indeed, the correlation matrix between the first -0.04) and that between the first and the third series (0.18 -0.10) ). Such an 0.22 0.22
two series ( (-0.31 0.96
observation justifies graph structure learning among the PMUs.
important note a few processing steps of the data set because of its noisy and incomplete nature. The data set contains a fair amount of unreliable readings (e.g., outliers). Hence, consult domain experts and set lower and upper bounds to filter out extremely small and large values. Ac- counting for missing data, within every five minutes take the mean of the available readings if any, or impute with the mean of the entire series.
B ADDITIONAL DETAILS OF EXPERIMENT SETTING
Hyperparameters. Several hyperparameters are tuned through grid search: initial learning rate 10.1, 0.01, 0.001), dropout rate 10. 0.2, 0.3), embedding size of LSTM (32, 64, 128, 256). the k value in KNN {5, 10, 20, 30). and the weight of regularization {0. 1, 2, 5, 10, 20). For other hyperparameters, the convolution kernel size in the feature extractor is 10 and the decay ratio of learning rate is 0.1. After tuning, the best initial learning rate for METR-LA and PEMS-BAN is 0.01 and for PMU is 0.001. The optimizer is Adam.
Because the loss function is expectation (see (1) and (2), the expectation is computed as an average of 10 random samples. Such an averaging is needed only for model evaluation. In training, one random sample suffices because the optimizer is stochastic optimizer.
Platform. We implement the models in Porch.. All experiments are run on one compute node of an IBM Power9 server. The compute node contains 80 CPU cores and four NVidia V100 GPUs, but because of scheduling limitation we use only one GPU.
Code is available htttts////gtttbb...c///aaaaaaaas
12

---

Published as a conference paper at ICLR 2021
C ADDITIONAL RESULTS FOR FORECASTING QUALITY
See Table for the forecast results of PEMS-BAY. The observations are rather similar to those of Tables and in Section 4.2 In particular, GTS produces the best prediction in all scenarios and under all metrics.
0.880
0.010
0.009
0.875
0.008
0.870
0.007
0.005
0.860
0.855
FeB 02
2017
Feb
2017
Feb 04 2017
Fes 05 2017 Feb 06 2017
Fab 07 2017
Feb 08 2017
Feb 09 2017
0.003
time
0.013
0.870
0.012
0.011
0.010
0.009
0.008
0.860
0.007
0.006
0.855
Fab
0.005
time
0.876
0.008
0.874
0.007
0.006
0.872
0.870
0.868
0.886
0.864
Feb 2017 Fes 2017
2017 Fet 05 2017
Feb 06 2017 time
Feb 07 2017
Feb 2017 Feb 2017
0.000
Figure 7: Example time series from the PMU data set. Data have been standardized and the vertical axes not show the raw value.
Table 3: Forecasting enror (PEMS-BAY)
Metrie
ARIMA
VAR
SVR
FNN
LSTM
DCRNN
LDS
GTSv
GTS
min 15 mi
MAE
2.88
1.62
1.74
1.85
2.20
2.05
1,38
1.33
1.16
1.12
RMSE
5.59
3.30
3.16
3.59
4.42
4.19
2.95
2.81
2.27
2.16
MAPE
6.8%
3.5%
3.6%
3.8%
5.2%
4.8%
2.9%
2.8%
2.3%
2.3%
MAE
2,88
2.33
2.32
2.48
2.30
2.20
1.74
1.67
1.44
1.34
RMSE
5.59
4.76
4.25
5.18
4.63
4.55
3,97
3.80
2.97
2.74
MAPE
6.8%
5.4%
5.0%
5.5%
5.4%
5.2%
3.9%
3.8%
3.1%
2.9%
min
MAE
2.88
3,38
2.93
3.28
2.46
2.37
2.07
1.99
1.81
1.58
RMSE
5.59
6.50
5.44
7.08
4.98
4,96
4.74
4.59
3.78
3.30
MAPE
6.8%
8.3%
6.5%
8.0%
5.9%
5.7%
4.9%
4.8%
4.1%
3.6%
13

---

Published as a conference paper at ICLR 2021
D UPDATES OF TABLES 1 2 AND 3
Our implementation had been developed based on the PyTorch version of DCRNN httpps: 7/github. com/ chnsh/DCRN? http:/DDDCNN_PIIoccc_c PyTorch). was brought our attention recently that this version calculated the evaluation metrics MAE/RMSE/MAPE in a manner slightly different from that used report the results of DCRNN the official publication (https:/
com/ Chnsh/DCRNV PyTorch/ Sarrchisssuee isssusssss We updated Tables I and 3 by correcting the cal- culations be consistent with the official DCRNN results. See Tables 5 and 6 Despite the correction, observations and conclusions regarding the comparison of different methods remain un- changed.
Table 4: (Correction of Table
Forecasting error (METR-LA).
Metric
HA
ARIMA
VAR
SVR
FNN
LSTM
DCRNN
LDS
GTSv
GTS
min
MAE
4.16
3.99
4.42
3.99
3.99
3.44
2.77
2.75
2.74
2.64
RMSE
7.80
8.21
7.89
8,45
7.94
6.30
5.38
5.35
5.09
4.95
MAPE
13.0%
9.6%
10.2%
9.3%
9.9%
9.6%
7.3%
7.1%
7.3%
6.8%
60 min 30 min
MAE
5.15
5.41
5.05
4.23
3.77
3.14
3.01
RMSE
7.80
10.45
9.13
10.87
8.17
7.23
6.45
6.45
6.02
5,85
МАРЕ
13.0%
12.7%
12.7%
12.1%
12.9%
10.9%
8.8%
8.6%
8.7%
8.2%
MAE
4.16
6.90
6.52
6.72
4.49
4.37
3.60
3,63
3.53
3.41
RMSE
13.23
13.76
8.69
7.67
MAPE
13.0%
17.4%
15.8%
16.7%
14.0%
13.2%
10.5%
10.3%
10.3%
9.9%
Table 5: (Correction of Table
Forecasting emror (PMU).
Metric
FNN
LSTM
DCRNN
LDS
NRI
GTSv
GTS
60 min min 15 min
MAE (x10
1.23
1.02
0.71
0.49
0.66
0.35
0.26
RMSE (x10
1.28
1.63
1.42
1.26
0.27
0.22
0.20
MAPE
0.20%
0.21%
0.09%
0.07%
0.14%
0.06%
0.04%
(x10
1.42
0.81
0.71
0.45
0.37
RMSE (x10
1.81
2.06
1.91
1.79
0.30
0.29
0.26
MАPE
0.23%
0.20%
0.15%
0.12%
0.15%
0.07%
0.06% %
MAE (x10
1.88
1.79
1.78
1.45
0.83
0.63
0.59
RMSE (x100
2.58
2.75
2.65
2.54
0.46
0.43
0.41
MAPE
0.29%
0.27%
0.24%
0.22%
0.17%
0.11%
0.10%
Table 6: (Correction of Table
Forecasting error SIEMS-BAY)
Metric
ARIMA
VAR
SVR
FNN
LSTM
DCRNN
LDS
GTSv
GTS
30 min 15 min
MAE
2.88
1.62
1.74
1.85
2.20
2.05
1.38
1.33
1.35
1.32
RMSE
5.59
3,30
3.16
3.59
4.42
4.19
2.95
2.81
2.64
2.62
MAPE
6.8%
3.5%
3.6%
3.8%
5.2%
4.8%
2.9%
2.8%
2.9%
2.8%
MAE
2.33
2.32
2.48
2.30
2.20
1.74
1.67
1.64
RMSE
5.59
4.76
4.25
5.18
4.63
4.55
3.97
3.80
3.45
3.41
MAPE
6.8%
5.4%
5.0%
5.5%
5,4%
5.2%
3.9%
3.8%
3,9%
3.6%
min
MAE
2.88
3.38
2.93
3.28
2.46
2.37
2,07
1.99
1.99
1.91
RMSE
5.59
6.50
5.44
4.98
4.96
4.74
4.59
4.05
3.97
MAPE
6.8%
8.3%
6.5%
8.0%
5.9%
5.7%
4.9%
4.8%
4.7%
4.4%
14


Reviews:
[{"review_id": "o422tncF9cr", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Paper summary:\n\nThis paper proposes an approach for time series forecasting that learns the graph structure among multiple (multivariate) time series simultaneously with the parameters of a Graph Neural Network (GNN). The problem is formulated as learning a probabilistic graphical model by optimizing the expectation over the graph distribution, which is parameterized by a neural network and encapsulated in a single differentiable objective. Empirical evidence suggests that the proposed GTS obtains superior forecasting performance to both deep and non-deep learning based, as well as graph and non-graph based, competitor forecasting models. In addition, GTS appears to be more computationally efficient compared to LDS, a recently proposed meta-learning graph-based approach.\n\n##########################################################################\n\nStrong points:\n1. A time series forecasting model is proposed to automatically learn a graph structure among multiple time series and forecast them simultaneously using a GNN.\n\n2. The graph structure and the parameters of the GNN are learned simultaneously in a joint end-to-end framework.\n\n3. The graph structure is parameterized by neural networks rather than being treated as a (hyper)parameter, thus significantly reducing the training cost compared with the recently proposed bilevel optimization approach LDS. \n\n4. A structural prior-based regularization is incorporated in GTS. In case a “ground-truth” graph is provided upfront, this may serve as a healthy variation of such a graph for the purpose of more accurate forecast.\n\n5. Extensive experiments are conducted in which the proposed GTS is compared to a number of baselines, including a recently proposed graph structure learning approach, and deep or non-deep learning based (as well as graph or non-graph based) forecasting models.\n\n6. The experimental results demonstrate that GTS outperforms its competitor approaches in terms of forecasting accuracy and is more efficient that the recently proposed LDS.\n\n7. Generally, the paper is well written, while the notation is clear and easy to follow.\n\n##########################################################################\n\nImprovement points:\n1. In section 3.4. (Comparison with NRI), the authors state that the “structural prior” $A^{a}$ offers a stronger preference on the existence/absence of each edge than a uniform distribution over all edges. This seems a bit unclear, thus I would encourage the authors to elaborate a bit more on this difference between GTS and NRI w.r.t. the structural prior.\n\n2. In the case of the PMU dataset, despite the fact that the grid topology is not provided, the authors still consider a certain structural prior by constructing a kNN graph among the PMUs. I am wondering whether the correlation between the series (mentioned briefly in Appendix A) is used for the graph construction or another distance/similarity metric is considered?\n\n3. Two variables are recorded by the 42 PMUs, however each node in the constructed graph (shown in Fig. 5) corresponds to one PMU. In case a single node corresponds to a single PMU, then I wonder how the similarity between two PMUs’ recordings is calculated across the two variables (voltage magnitude and current magnitude)?\n\n4. The authors construct the PMU dataset by extracting only one month of data. However, a single month of PMU data would not allow for capturing certain long-term seasonalities (for instance, the PMU recordings are typically impacted by outages that occur more frequently in certain seasons or periods in the year). Is this perhaps due to data unavailability? If that is not the case, I would ask the authors to clarify the reasoning behind the decision to extract the data for February 2017?\n\n5. In Tables 2 & 3 (Appendix C), some of the MAPE values obtained by GST are bolded even though the same percentages are reported for GTSv. In such cases, I would suggest the authors to either bold the MAPEs obtained by both GTS and GTSv, or present the MAPE values using more decimal places.\n\n6. There are several minor textual errors throughout the paper that can be easily addressed. Some of them are summarized as follows:\n- The term “LDS” is initially used at the beginning of page 2, but is not defined earlier in the text.\n- In the third paragraph on page 2, “computation is expensive” should be replaced by “its computation is expensive”.\n- I am wondering whether the training loss $L$ should be used instead of the validation loss $F$ in Eq. (2)? If so, correct accordingly, otherwise disregard this comment.\n- In the next-to-last paragraph on page 2, consider replacing “it is better scaled” with “it scales better”.\n- In the last paragraph of the Related Work section, “of node classification tasks” should be replaced by “on node classification tasks”.\n- At the beginning of section 3, the term “NRI” is used, but is not defined earlier in the text.\n- On page 5, consider replacing the abbreviation “ELBO” with “evidence lower bound (ELBO)”. \n- In the first paragraph on page 7, “treating it a (hyper)parameter as in LDS” could be replaced with “treating it *as* a (hyper)parameter *which is the case* in LDS”.\n- In the second paragraph on page 8, replace “regularization $\\lambda$” with “regularization strength $\\lambda$”. In the same sentence, consider adding “the” before both “forecasting error” and “cross entropy”.\n\n##########################################################################\n\nQuestions during rebuttal period:\nPlease address the aforementioned remarks/questions.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "HTeDz0UymjB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper considers learning both graph structures and NNs for time series data, similar to the idea of LDS (Franceschi et al., 2019). Observing the computation and scalability issues with LDS, authors propose a unilevel optimization form wrt. the mean performance over the graph distribution. This is done via NNs, with input being the observed sequence, to output a real matrix whose elements are then treated as weights for the Gumbel trick. NN structures, training procedure, etc. mostly follow existing works.\n\nOverall, the paper is well presented, easy to understand, with a simple and somewhat effective modification over LDS. I generally like simple ideas with sufficient insights and explanations (though there is not much in this work), but I'm not sure if the empirical improvement is sufficient. I recommend a weak acceptance for now and may change my score after reading other reviews.\n\nI only have one question: the proposed idea is not restricted to time-series case. So how does it perform for non-time-series data? It would be a big benefit if the proposed idea also helps in a more general case than the present scope.\n\n** after reading response ** I thank authors for replying to my question. I maintain the previous rating.\n\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "07INhWEkW-Q", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 4, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes an approach for multivariate time series forecasting by trying to estimate dependence across dimensions via a learned graph structure.  The dimensions are considered as nodes in a graph, and the problem is mapped to learning a discrete graph structure that can help with downstream forecasting task.\nThe paper shows that a graph neural network (GNN) can be leveraged even though an explicit structure is unknown, to improve forecasting performance. This is achieved while learning the graph structure and forecasting architectures in an end-to-end fashion. The proposed approach is computationally efficient compared to a bilevel optimization approach where a discrete graph structure is learnt in a meta-learning framework. The approach is further claimed to be able to incorporate apriori knowledge of the graph structure by proposing a regularizer that ensures that the learned graph structure stays close to the known graph structure. The proposed approach improves forecasting performance in comparison to several strong baseline methods on three real-world datasets.\nIn general, the paper is well-written and easy to follow. \n\nAttempts have been made in the past for learning such a discrete graph structure from data. The authors mention LDS [2] and NRI [3] as closest to their work. The authors attempt to explicitly compare the proposed approach to NRI. The authors claim that \"The most essential distinction is the number of structures\": one structure is learned in the proposed approach while many structures are learned in NRI. From what I could follow, this single structure is achieved by using the entire multivariate time series data to obtain a feature vector for each dimension (series) via a neural network instead of using window-wise data. In this sense, this appears to be a simplification of NRI, rather than being something novel and different.\nThe proposed setup and the approach are different and novel compared to NRI as the \"Amortized inference is not desired nor relevant\": I am not sure how this makes the proposed approach non-trivial given NRI? \nFurthermore, in contrast to LDS, the key contribution of the proposed approach is to get rid of the bilevel optimization. But then, that also seems to rely mainly on the Gumbel reparameterization trick which has been used in NRI for forecasting albeit for a slightly different setting.\n\nAnother very closely related approach to the proposed one is that in [1], which the authors seem to be unaware of. \nOne of the main claims of the proposed approach is an attempt to learn the graph structure and forecasting model in an end-to-end learning fashion. However, this problem has already been attempted in [1].\nTherefore, it is difficult to comment on the novelty and contribution of the paper without a comparison with [1], especially since most of the benchmark tasks and datasets used in this paper are present in [1] as well. \n\nI have some concerns regarding the empirical evaluation:\n1. Can the observations in Fig. 2 be attributed to the graph learning part? Despite the fact that the only difference between DCRNN and the proposed method seems to be the graph structure learning part, it is still not obvious qualitatively as to why the observations of Fig. 2 can be attributed to the graph learning part, e.g. why is \"GTS curve better captures the sharp dip toward the end of the series\" attributable to the graph learning part qualitatively or as per domain knowledge? I think an empirical analysis on a synthetic dataset to support such claims related to ablation could be useful.\n2. In Fig. 4a, the regularization seems to induce sparsity and has been observed by the authors as well: \"increasing the regularization magnitude produces sparser graphs\". But the efficacy of such regularization on forecasting performance is not clear as k=0 (no kNN regularization) seems to have the best forecasting performance in Fig. 4a. This seems to imply that using kNN graph knowledge is not adding any value.\nSimilar observations can be made in Fig. 4b, where increasing $\\lambda$ leads to increasing MAE. As such, the effect or usefulness of regularization is not clear.\n3. The analysis on regularization and learned structures is done using a kNN graph as apriori knowledge for the PMU dataset. Rather than relying on another data-driven graph structure (kNN graph) as ground truth, I wonder if it would be useful to do such analysis on the public datasets (METR-LA and PEMS-BAY) for which the ground truth structures are actually known. As such, the evaluations on \"effect of regularization\" and \"learned structures\" do not seem conclusive.\n4. In Tables 2 and 3, some of the bolds also apply to GTSv but are missing.\n\nGiven the above points, the originality of the work and the contributions are not clear.\n\nOther minor points:\n1. Related Work can also benefit from more precise references to papers that use the mentioned architectures for time series forecasting. The current references are too generic.\n2. typo: hyperparemeter\n\n\nReferences:\n[1] Connecting the Dots: Multivariate Time Series Forecasting with Graph Neural Networks, Wu et. al, KDD2020. https://dl.acm.org/doi/abs/10.1145/3394486.3403118\n[2] LDS: Learning discrete structures for graph neural networks, Franceschi et. al, ICML, 2019\n[3] NRI: Neural relational inference for interacting systems. Kipf et. al, ICML, 2018.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a graph neural network-based approach to forecasting multiple time series. It incorporates structure learning similar in some ways to NRI (Kipf et al.) and a recurrent graph convolution forecaster given the inferred graph based on DCRNN (Li et al.). The paper shows consistently improved performance over several different kinds of baselines (classical, deep learning, and graph-based deep learning) on three datasets, two of which are public and one proprietary.

Reviewer 1 thought the paper was “well presented and easy to understand”. I agree. The reviewer liked the simplicity of the approach but wondered whether the empirical improvement was sufficient. The reviewer asked about applicability outside the time series domain and the authors provided a satisfying response.

Reviewer 3 thought that simultaneously learning graph structure and forecasting was an “understudied topic”. The reviewer pointed out several strengths including: the end-to-end nature of the approach, the reduction in training cost due to the direct parameterization of the adjacency matrix structure, the optimal structural regularization scheme, and the extensive experimentation. Like R1, they thought the paper was well written. They suggested several points of improvement, which were mainly seeking clarification. They made a good point regarding the PMU dataset in that only one month was considered, which was not enough to capture long-term seasonalities. This seems like a limitation to me. The authors responded to each of the points in turn. Regarding the point about the PMU dataset, the authors stated that the data was extremely noisy, so they settled on a month that was comparatively clean.

The review from R4 was not as positive as the other reviews. R4 proposed that the present work may be a “simplification of NRI” rather than being novel/different because NRI was a windowed-based approach, while GTS was based on the entire time series. The reviewer also pointed out what appears to be a highly relevant paper (Wu et al.); though this appeared in KDD 2020 and in my opinion it’s understandable if the authors missed it. Finally, R4 raised several issues with the empirical evaluation. They make a very good point that the analysis on regularization used the kNN graph where a “ground truth” graph was not available. Why not evaluate the regularizer on the other datasets where ground truth is available? The authors responded with an updated paper addressing this point. The authors responded to other points of criticism and a fairly extensive debate ensued.
The key points of the debate were:
A difference in opinion between the significance in departure between the structure learning mechanism in this work (GTS) and NRI (Kipf et al.)
Overlap between a recently proposed paper (Wu et al., KDD 2020). I am fairly sympathetic with the authors here. That work is fairly late-breaking and they do point out a key difference: "the advantage of our method over MTGNN is clear: we can get structures we desire, not restricted to a degree-k graph. This advantage is an important contribution, because MTGNN hard-wires the graph parameterization and enforces that all nodes have the same number of (out) neighbors. What if one desires instead a graph that approximately obeys spatial proximity, like that in the case of METR-LA? Our method yields such a graph."
Some minor concerns with Fig. 2 and the structural prior/regularization analysis

I have read the paper, and while R4’s concerns are legitimate, I think this paper is clearly over the bar. I support this paper’s acceptance and ask the authors to take the reviews into consideration when revising their paper. As R4 suggests, they could add a controlled set of experiments on a perhaps synthetic dataset to show that the proposed GTS is better than NRI. If the scalability of NRI is a concern, then this can be highlighted to the same extent as LDS.
Decision: Accept (Poster)