Year: 2021

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(N2) 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.

text: 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(N2) 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.

heading: 1 INTRODUCTION
text: Time Series Forecasting (TSF) has been widely studied due to its practical significance in a 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., 2005), traffic control (Li et al., 2017) and medicine (Kaushik et al., 2020). Petropoulos et al. (2020) provide a non-systematic overview of further applications. In TSF, given a sequence of data points indexed over time, we aim to estimate its future values based on previously observed data. Data is often multivariate, meaning that multiple variables vary over time, each variable may not only depend on 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 to improve the forecasting accuracy. Some classical MTS forecasting algorithms such as Vector Autoregression (VAR) (Lütkepohl, 2005) or Gaussian Processes (Roberts et al., 2013) only consider linear dependencies among variables. A natural way to model non-linear 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., 2020a), 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 et al., 2019; Wu et al., 2020; Shang et al., 2021). In this case, a latent graph representation is inferred while forecasting, allowing to operate on a larger variety of datasets. Our work belongs to this category.
However, inferring all pairwise relations may come at a higher computational cost, as the number of relations scales quadratically O(N2) w.r.t. 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 a new latent graph inference algorithm for Multivariate 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 a more flexible representation than a static graph for the whole dataset. Additionally, we optionally reduce the complexity from O(N2) (Fully Connected Assumption) to O(NK) (Bipartite Assumption) where K ⌧ N for a small trade off in performance.

heading: 2 BACKGROUND

heading: 2.1 TIME SERIES FORECASTING
text: In time series forecasting we want to estimate a future time series xt+1:T given its past xt0:t where t0  t  T indexes over time, and (optionally) some context information c. For the multivariate case we assume the time series is composed of N variates at a time such that xt0:T = {x1,t0:T , . . . ,xN,t0:T } 2 RN⇥T t0+1. In 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 ✓u are shared across all univariate time series. More formally:
x̂i,t+1:T = fu(xi,t0:t, ci; ✓u) (1)
where x̂ denotes the estimated values, i 2 {1, . . . , N} indexes over multivariates and fu(·) is the estimator function with learnable parameters ✓u 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 f✓ and already provides reasonable results. A popular example of univariate models would be (Salinas et al., 2020).
Multivariate methods: Multivariate methods condition on all past data (all N variates) and directly predict the multivariate target. More formally:
x̂t+1:T = fm(xt0:t, c). (2)
Different variables may be correlated and/or depend on the same 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 a 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 more expensive algorithms, specially when trying to model all pairwise interactions between variates.

heading: 2.2 GRAPH NEURAL NETWORKS
text: Graph Neural Networks (GNNs) (Bruna et al., 2013; Defferrard et al., 2016; Kipf and Welling, 2016) operate directly on graph structured data. They have gained a 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 et al., 2017). Given a graph G = (V, E) with nodes vi 2 V and edges eij 2 E , we define a graph convolutional layer as:
mij = e(h l i,h l j) mi =
X
j2N (i)
↵ijmij h l+1 i = h(h l i,mi) (3)
Where e and h are the edge and node functions, usually approximated as Multi Layer Perceptrons (MLPs), hli 2 Rnf is the nf-dimensional embedding of a node vi at layer l and mij is the edge
embedding that propagates information from node vj to vi. A GNN is constructed by stacking multiple of these Graph Convolutional Layers hl+1 = GCL[hl, E ]. Additionally, in Equation 3 we include ↵ij 2 (0, 1) which is a scalar value that performs the edge inference or attention over the neighbors similarly to Veličković et al. (2017). As done in (Satorras et al., 2021), we choose this value to be computed as the output of a function ↵ij = ↵(mij) where ↵ is composed of just a linear layer followed by a sigmoid activation function.

heading: 3 RELATED WORK
text: 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; Benidis et al., 2020; Lim and Zohren, 2021). Traditionally, most classical methods are univariate in 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 et al., 2019; Montero-Manso and Hyndman, 2022). These global models estimate their parameters over an entire panel of time series, so thereby taking advantage of cross time series learning, but still produce a 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 et al., 2019), CNNs (Wen et al., 2017; Chen et al., 2020b), Transformers (Li et al., 2019; Lim et al., 2021; Eisenach et al., 2020) and also combining classical probabilistic models with deep learning (Rangapuram et al., 2018; Kurle et al., 2020; de Bézenac et 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-variate likelihood (Rasul et al., 2020; 2021; de Bézenac et al., 2020; Salinas et al., 2019). Here, we attempt to capture the multivariate nature of many modern forecasting problems primarily by using a multi-variate time series as input. For this, a natural way to model and exploit the relationship between time series is via Graph Neural Networks (Bruna et al., 2013; Defferrard et al., 2016; Kipf and Welling, 2016) which have been successfully applied to a wide variety of deep learning domains where exploiting relations between entities/nodes can benefit the prediction task. Even in those cases where edges are not explicitly provided in the dataset, attention or a 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, (Kazi et al., 2020) in zero-shot learning and 3D point cloud segmentation, (Garcia and Bruna, 2017; Kossen et al., 2021) in image classification, (Cranmer et al., 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); Shih et al. (2019) are some of the first deep learning approaches designed to exploit those pair-wise dependencies. More recent methods, (Li et al., 2017; Yu et al., 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., 2021), and in MTS anomaly detection (Zhang et al., 2020; Deng and Hooi, 2021). These methods can be applied to any dataset even when there is not an 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 or O(N3) in (Cao et al., 2021).
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 xt0:t, this also makes the
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 O(NK) where K ⌧ N with the bipartite assumption.

heading: 4 METHOD
text: In Section 2.1 we 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 1. We can break down the Univariate Model in two main steps fu = fenc fdec such that xi,t0:t fenc ! zi fdec ! x̂i,t+1:T , where fenc encodes the input signal xi,t0,t (and optionally some context information ci) into an embedding zi and fdec estimates the future signal from this embedding. In our method, we include a multivariate aggregation module AGG in between fenc and fdec that propagates information among nodes in the latent space z = {z1, . . . , zN}. This aggregation module takes as input the embedding zi = fenc(xi,t0:t, ci), and outputs a modified embedding ẑ = AGG(z) where information has been propagated among nodes. These new embeddings are then passed as input to the decoder fdec. The resulting algorithm is:
Univariate Encoder zi = fenc(xi,t0:t, ci) (4) Multivariate extension ẑ = AGG(z) (5)
Univariate Decoder x̂i,t+1:T = fdec(ẑi) (6)
Notice the overall model is multivariate but fenc and fdec remain univariate, allowing a modular extension from current univariate methods to multivariate. Additionally, in contrast to the most recent methods (Wu et al., 2020; Shang et al., 2021), our model does not propagate information among nodes at every time step [t0, . . . , t] 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 zi as context information ci = 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 a standard GNN. It results in good performance, but its complexity scales O(N2) 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.2 as the aggregation module AGG. We assume a fully connected graph G = {V, E} where each time series embedding zi is associated to a node of the graph i 2 V , specifically, the embedding zi is provided as the input h0i to the GNN from Equation 3, then the GNN runs for L layers and the output node embedding hLi is provided as the input ẑi of the decoder. Moreover, despite the fully connected assumption, the GNN infers attention weights ↵ij 2 (0, 1) (eq. 3) 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, we name this model a Fully Connected Graph Neural Network (FC-GNN). In summary, the
aggregation module AGG under this FC-GNN setting, is able to 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" or multiplied by an inferred scalar value ↵ij .
Bipartite graph assumption | BP-GNN In the previous section we 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 a bipartite graph assumption which reduces the complexity to O(NK), K being a parameter of choice K ⌧ N . To accomplish this, we define a bipartite graph G = (Y,U , E) with two sets of nodes Y and U . Here Y is a set of N nodes corresponding to the N dimensions of the multivariate time series and U is a set of K auxiliary nodes. Each time series embedding zi is associated to one of the nodes i 2 Y . Additionally, we define embeddings u = {u1, . . .uK} associated to the auxiliary set U . Edges E interconnect all nodes between the two subsets {Y , U}, but there are no connections among nodes of the same subset. This results in 2NK edges, with a computation complexity of O(NK).
The algorithm works in the following way. We input into the GNN the union of the two node subsets V = Y[U . Specifically, the input embedding h0 defined in Equation 3 is the concatenation of the time series embeddings z with the auxiliary node embeddings u (i.e. h0 = z||u) 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 ! U , next the other way around U ! Y . Notice the auxiliary nodes U do not contain any information about the data until they have received messages from Y , this is why we chose this schedule. A visual illustration of the whole algorithm is presented in Figure 1, where the bipartite graph is depicted at the top right of the image (BP-GNN).
We have conceptually defined the Bipartite Graph Neural Network. In Table 1 we introduce the equations that formally define it as an extension of the standard GNN Equation 3. Notice that it can be simply formulated as a two steps process where the indexes i, j belong to each one of the subsets U or V depending on the direction of the messages (V ! U or U ! V). Additionally, we used different learn-
able parameters between Step 1 and Step 2 in the modules e, ↵ and h since it resulted in better accuracy than sharing parameters. In Appendix A.1 we provide the BP-GNN equations from Table 1 in its full form. Following, we define the adjacency matrices corresponding to the two message passing steps (assuming all inference parameters ↵ij = 1):
A1 = 0N⇥N 0N⇥K 1K⇥N 0K⇥K , A2 = 0N⇥N 1N⇥K 0K⇥N 0K⇥K , Ã = A2A1 = KN⇥N 0N⇥K 0K⇥N 0K⇥K (7)
The two adjacency matrices A1 and A2 are non-symmetric and they define a directed message passing scheme. A1 refers to Y ! U and A2 refers to U ! Y . After running these two propagation steps, all time series nodes Y will have received messages from all other time series nodes Y as in the fully connected graph case, but the number of messages have been reduced (iff 0  K < N/2) due to the factorization of the adjacency matrix. The first N rows and columns in matrix Ã from Equation 7 represent the sum of all paths that communicate the time series nodes Y among them after the two updates A2A1.
Architecture details As explained in this Section 4, our method is composed of three main modules, the encoder fenc, the decoder fdec and the aggregation function AGG. We choose to use relatively simple networks as encoder and decoder. The decoder fdec is defined as a Multi Layer Perceptron (MLP) with a single hidden layer and a residual connection in all experiments. The encoder fenc 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 a larger encoding length and the translation equivariance of CNNs showed to be more beneficial. The encoder fenc, first encodes the input signal xi,t0:t to an embedding vector by using the mentioned MLP or CNN, and then concatenates a unique
identifier (ci = id) to the obtained embedding vector resulting in zi. ci could (optionally) include additional context information if it was provided in the dataset. The combination of these simple networks with our proposed aggregation module AGG fully defines our model. All architecture details are explained in detail in Appendix A.2.
The a aggregation module AGG was 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 e, h and ↵. e consists of a two layers MLP, h is a one layer MLP with a skip connection from he input to the output and ↵ 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 Mean Squared Error between the estimated values x̂t+1:T from Equation 6 and the ground truth as L = l(x̂t+1:T ,xt+1:T ).

heading: 5 EXPERIMENTS
text: 5.1 DATASETS AND BASELINES
We first evaluate our method in METR-LA and PEMS-BAY datasets from (Li et al., 2017) which record traffic speed statistics on the highways of Los Angeles county and the Bay Area respectively. We also evaluate in the publicly available Solar-Energy, Traffic, Electricity and Exchange-Rate. Specifications for each dataset are presented in Table 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 a single model is trained for all time series but they are treated independently without message exchange among them. These baselines include a simple linear Auto Regressive model (AR) and a 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 a multilayer perceptron defined by h in Equation 3.
• Multivariate with a known graph: These methods require a 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 et al., 2017), STGCN (Yu et al., 2017) and MRA-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), TPALSTM (Shih et al., 2019), MTGNN (Wu et al., 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 predefined adjacency matrix. Comparisons to NRI (Kipf et al., 2018) can be found in previous literature (Shang et 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://github. com/chaoshangcs/GTS.

heading: 5.2 MAIN RESULTS
text: In this section we evaluate our method in METR-LA and PEMS-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
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 a batch size 16 in a Tesla V100-SXM GPU. Further implementation details are provided in Appendix B.1.
Results are reported in tables 11 and 4, 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 information 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 FCGNN 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 to the dynamical inference. Time results are presented in Table 4. BP-GNN is the most efficient algorithm in both METR-LA and PEMSBAY by a large margin. FC-GNN is also more efficient than previous methods in both datasets but it is still limited by the O(N2) scalability. In Section 5.3 we will evaluate in a 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 O(N). In this timing benchmark (Table 4) we included all methods from the previous Table 11 that have a publicly available implementation in METR-LA and PEMS-BAY datasets.

heading: 5.3 SINGLE STEP FORECASTING
text: 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.,
2019; Wu et al., 2020) firstly proposed by TLST-skip (Lai et al., 2018). In this experiment the network 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 (MAE). All datasets have been split in 60%/20%/20% for training/val/test respectively. As proposed in (Lai 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 2 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.1. Results with standard deviations are provided in Appendix C.3
Results in Table 5 are consistent with the previous experiment. FC-GNN outperforms all previous works in most metrics. There is a 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 metrics and very close to FC-GNN. Regarding the timing results (Table 6), BP-GNN is the most efficient graph inference method by a 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-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 FCGNN 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 2 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).

heading: 5.4 INFERRED GRAPH ANALYSIS
text: 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.
(a) "Cycle Graph Gaussians" samples each series value xi,t 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 xi,t ⇠ N ( xi 1 mod N,t 5; 2), where = 0.9 and = 0.52. (b) We call the second dataset "Noisy Sinusoids", inspired by the Discrete Sine Transformation we generate arbitrary signals as the sum of different sinusoids plus gaussian noise. We share the same sinusoids among the first half 1  i  5 of the panel and among the second half 6  i  10 while sampling different gaussian noise for each signal i. 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.1. Both BP-GNN and FC-GNN consist of a 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 t. For BP-GNN we construct the adjacencies as the first N rows and columns of Ã = A2A1 (Eq. 7, 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.
Results are reported in Figure 2. FC-GNN perfectly matches the ground truth adjacencies and it 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 to obtain very competitive performance (even for small K) by inferring suboptimal matrices with denser connections. In the "Noisy Sinusoids" dataset, the benefit of increasing K is smaller, this is coherent with the fact the ground truth matrix is 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.

heading: 5.5 CHOOSING THE NUMBER OF AUXILIARY NODES K
text: In the previous section we evaluated BP-GNN for different K values, and we found small values were 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 K values. Therefore, we chose a relatively small number of auxiliary nodes K = 4 that resulted in a good
trade-off accuracy vs computation. In Figure 3, we plot the MAE for different K values in METR-LA for a 15 min time horizon, and we see the best performances are already achieved at a small K values. In the middle plot of Figure 3, we report the forward running time in (s) for a batch size 16 for METR-LA as we sweep different K values. Finally, the right plot of Figure 3 shows the scalability of the number of edges w.r.t. to the number of nodes for O(N2) (e.g. FC-GNN) and O(2N4) (i.e. BP-GNN (K=4). This illustrates the vast difference in scalability as the number of nodes increases.

heading: 6 CONCLUSIONS
text: We presented a 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 a latent graph in the space of the embeddings of the univariate time series and can be inserted in many neural-network 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.

Review summaries:
The paper proposes using graph neural net (GNN) operations to combine per-series embeddings, to enable multivariate forecasting.  Specifically, 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.  The 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.  The 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.

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.