Year: 2022

Abstract:
Dealing with missing values and incomplete time series is a labor-intensive, tedious, inevitable task when handling data coming from real-world applications. Effective spatio-temporal representations would allow imputation methods to reconstruct missing temporal data by exploiting information coming from sensors at different locations. However, standard methods fall short in capturing the nonlinear time and space dependencies existing within networks of interconnected sensors and do not take full advantage of the available – and often strong – relational information. Notably, most state-of-the-art imputation methods based on deep learning do not explicitly model relational aspects and, in any case, do not exploit processing frameworks able to adequately represent structured spatio-temporal data. Conversely, graph neural networks have recently surged in popularity as both expressive and scalable tools for processing sequential data with relational inductive biases. In this work, we present the first assessment of graph neural networks in the context of multivariate time series imputation. In particular, we introduce a novel graph neural network architecture, named GRIN, which aims at reconstructing missing data in the different channels of a multivariate time series by learning spatio-temporal representations through message passing. Empirical results show that our model outperforms state-of-the-art methods in the imputation task on relevant real-world benchmarks with mean absolute error improvements often higher than 20%.

heading: 1 INTRODUCTION
text: Imputation of missing values is a prominent problem in multivariate time-series analysis (TSA) from both theoretical and practical perspectives (Little & Rubin, 2019). In fact, in a world of complex interconnected systems such as those characterizing sensor networks or the Internet of Things, faulty sensors and network failures are widespread phenomena that cause disruptions in the data acquisition process. Luckily, failures of these types are often sparse and localized at the single sensor level, i.e., they do not compromise the entire sensor network at once. In other terms, it is often the case that, at a certain time step, missing data appear only at some of the channels of the resulting multivariate time series. In this context, spatio-temporal imputation methods (Yi et al., 2016; Yoon et al., 2018b) aim at reconstructing the missing parts of the signals by possibly exploiting both temporal and spatial dependencies. In particular, effective spatio-temporal approaches would reconstruct missing values by taking into account past and future values, and the concurrent measurements of spatially close neighboring sensors too. Here, spatial similarity does not necessarily mean physical (e.g., geographic) proximity, but rather indicates that considered sensors are related w.r.t. a generic (quantifiable) functional dependency (e.g., Pearson correlation or Granger causality – Granger, 1969) and/or that are close in a certain latent space. Relational information, then, can be interpreted as a set of constraints – linking the different time series – that allows replacing the malfunctioning sensors with virtual ones.
Among different imputation methods, approaches based on deep learning (LeCun et al., 2015; Schmidhuber, 2015; Goodfellow et al., 2016) have become increasingly popular (Yoon et al., 2018a; Cao et al., 2018; Liu et al., 2019). However, these methods often completely disregard available relational information or rely on rather simplistic modifications of standard neural architectures tailored ∗Equal contribution. Correspondence to andrea.cini@usi.ch.
for sequence processing (Hochreiter & Schmidhuber, 1997; Chung et al., 2014; Bai et al., 2018; Vaswani et al., 2017). We argue that stronger, structural, inductive biases are needed to advance the state of the art in time series imputation and allow to build effective inference engines in the context of large and complex sensor networks as those found in real-world applications.
In this work, we model input multivariate time series as sequences of graphs where edges represent relationships among different channels. We propose graph neural networks (GNNs) (Scarselli et al., 2008; Bronstein et al., 2017; Battaglia et al., 2018) as the building block of a novel, bidirectional, recurrent neural network for multivariate time series imputation (MTSI). Our method, named Graph Recurrent Imputation Network (GRIN), has at its core a recurrent neural network cell where gates are implemented by message-passing neural networks (MPNNs; Gilmer et al., 2017). Two of these networks process the input multivariate time series in both forward and backward time directions at each node, while hidden states are processed by a message-passing imputation layer which is constrained to learn how to perform imputation by looking at neighboring nodes. In fact, by considering each edge as a soft functional dependency that constraints the value observed at corresponding nodes, we argue that operating in the context of graphs introduces a positive inductive bias for MTSI. Our contributions are manifold: 1) we introduce a methodological framework to exploit graph neural networks in the context of MTSI, 2) we propose a novel, practical and effective implementation of a GNN-based architecture for MTSI, and 3) we achieve state-of-the-art results on several and varied MTSI benchmarks. Our method does not rely on any assumption on the distribution of the missing values (e.g., presence and duration of transient dynamics and/or length of missing sequences) other than stationarity of the underlying process. The rest of the paper is organized as follows. In Section 2 we discuss the related works. Then, in Section 3, we formally introduce the problem settings and the task of MTSI. We present our approach to MTSI in Section 4, by describing the novel framework to implement imputation architectures based on GNNs. We proceed with an empirical evaluation of the presented method against state-of-the-art baselines in Section 5 and, finally, we draw our conclusions in Section 6.

heading: 2 RELATED WORKS
text: Time series imputation There exists a large literature addressing missing value imputation in time series. Besides the simple and standard interpolation methods based on polynomial curve fitting, popular approaches aim at filling up missing values by taking advantage of standard forecasting methods and similarities among time series. For example, several approaches rely on k-nearest neighbors (Troyanskaya et al., 2001; Beretta & Santaniello, 2016), the expectation-maximization algorithm (Ghahramani & Jordan, 1994; Nelwamondo et al., 2007) or linear predictors and statespace models (Durbin & Koopman, 2012; Kihoro et al., 2013). Low-rank approximation methods, such as matrix factorization (Cichocki & Phan, 2009), are also popular alternatives which can also account for spatial (Cai et al., 2010; Rao et al., 2015) and temporal (Yu et al., 2016; Mei et al., 2017) information. Among linear methods, STMVL (Yi et al., 2016) combines temporal and spatial interpolation to fill missing values in geographically tagged time series.
More recently, several deep learning approaches have been proposed for MTSI. Among the others, deep autoregressive methods based on recurrent neural networks (RNNs) found widespread success (Lipton et al., 2016; Che et al., 2018; Luo et al., 2018; Yoon et al., 2018b; Cao et al., 2018). GRU-D (Che et al., 2018) learns how to process sequences with missing data by controlling the decay of the hidden states of a gated RNN. Cao et al. (2018) propose BRITS, a bidirectional GRU-D-like RNN for multivariate time series imputation that takes into account correlation among different channels to perform spatial imputation. Other successful strategies in the literature have been proposed that exploit the adversarial training framework to generate realistic reconstructed sequences (Yoon et al., 2018a; Fedus et al., 2018; Luo et al., 2018; 2019). Notably, GAIN (Yoon et al., 2018a) uses GANs (Goodfellow et al., 2014) to learn models that perform imputation in the i.i.d. settings. Luo et al. (2018; 2019) aim, instead, at learning models that generate realistic synthetic sequences and exploit them to fill missing values. Miao et al. (2021) use an approach similar to GAIN, but condition the generator on the predicted label for the target incomplete time series. Concurrently to our work, Kuppannagari et al. (2021) developed a graph-based spatio-temporal denoising autoencoder for spatio-temporal data coming from smart grids with known topology. Liu et al. (2019), instead, uses adversarial learning to train a multiscale model that imputes highly sparse time series in a hierarchical fashion. However, we argue that none of the above-cited methods can take full
advantage of relational information and nonlinear spatio-temporal dependencies. Most importantly, the above methods do not fully exploit the flexibility and expressiveness enabled by operating in the context of graph processing.
Graph neural networks for TSA Graph neural networks have been exploited in TSA mostly in spatio-temporal forecasting methods. The idea behind most of the methods present in the literature is to modify standard neural network architectures for sequential data by relying on operators that work in the graph domain. For example, Seo et al. (2018) propose a GRU cell where gates are implemented by spectral GNNs (Defferrard et al., 2016); Li et al. (2018) propose an analogous architecture replacing spectral GNNs with a diffusion-convolutional network (Atwood & Towsley, 2016). Note that these models are different w.r.t. approaches that use recurrent networks to propagate information graph-wise (Scarselli et al., 2008; Li et al., 2016). Yu et al. (2017) and Wu et al. (2019; 2020b) propose, instead, spatio-temporal convolutional neural networks that alternate convolutions on temporal and spatial dimensions. Similar approaches have also been studied in the context of attention-based models (Vaswani et al., 2017) with spatio-temporal Transformer-like architectures (Zhang et al., 2018; Cai et al., 2020). Another particularly interesting line of research is related to the problem of learning the graph structure underlying an input multivariate time series (Kipf et al., 2018; Wu et al., 2020b; Shang et al., 2020). While previously mentioned approaches focus on multivariate time series prediction, other methods aim at predicting changes in graph topology (Zambon et al., 2019; Paassen et al., 2020). Conversely, methods such as Temporal Graph Networks (Rossi et al., 2020) are tailored to learn node embeddings in dynamical graphs. Finally, recent works have proposed GNNs for imputing missing features in the context of i.i.d. data. Among the others, Spinelli et al. (2020) propose an adversarial framework to train GNNs on the data reconstruction task, while You et al. (2020) propose a bipartite graph representation for feature imputation. Lately, GNNs have also been exploited for spatial interpolation (Appleby et al., 2020; Wu et al., 2020a) – sometimes referred to as kriging (Stein, 1999). To the best of our knowledge, no previous GNN-based method targeted missing value imputation for generic multivariate time series.

heading: 3 PRELIMINARIES
text: Sequences of graphs We consider sequences of weighted directed graphs, where we observe a graph Gt withNt nodes at each time step t. A graph is a couple Gt = 〈Xt,Wt〉, where Xt ∈ RNt×d is the node-attribute matrix whose i-th row contains the d-dimensional node-attribute vector xit ∈ Rd associated with the i-th node; entry wi,jt of the adjacency matrix Wt ∈ RNt×Nt denotes the scalar weight of the edge (if any) connecting the i-th and j-th node. Fig. 1 exemplifies this modelling framework. We assume nodes to be identified, i.e., to have a unique ID that enables time-wise consistent processing. This problem setting can be easily extended to more general classes of graphs with attributed edges and global attributes. In this work, we mainly focus on problems where the topology of the graph is fixed and does not change over time, i.e., at each time step Wt = W and Nt = N .
Any generic multivariate time series fits the above framework by letting each channel of the sequence (i.e., each sensor) correspond to a node and using the available relation information to build an adjacency matrix. If no relational information is available, one could use the identity matrix, but this would defeat the purpose of the formulation. A more proper choice of Wt can be made using any standard similarity score (e.g., Pearson correlation) or a (thresholded) kernel. A more advanced approach instead could aim at learning an adjacency directly from data by using, for instance, spatial attention scores or resorting to graph learning techniques, e.g., Kipf et al. (2018). From now on, we assume that input multivariate time series have homogeneous channels, i.e., sensors are of the same
type. Note that this assumption does not imply a loss in generality: it is always possible to standardize node features by adding sensor type attributes and additional dimensions to accommodate the different types of sensor readings. Alternatively, one might directly model the problem by exploiting heterogeneous graphs (Schlichtkrull et al., 2018).
Multivariate time series imputation To model the presence of missing values, we consider, at each step, a binary mask Mt ∈ {0, 1}Nt×d where each row mit indicates which of the corresponding node attributes of xit are available in Xt. It follows that, m i,j t = 0 implies x i,j t to be missing; conversely, if mi,jt = 1, then x i,j t stores the actual sensor reading. We denote by X̃t the unknown ground truth node-attribute matrix, i.e., the complete node-attribute matrix without any missing data. We assume stationarity of missing data distribution and, in experiments, we mostly focus on the missing at random (MAR) scenario (Rubin, 1976). We neither make assumptions on the number of concurrent sensor failures, nor on the length of missing data blocks, i.e., multiple failures extended over time are accounted for. Clearly, one should expect imputation performance to scale with the number of concurrent faults and the time length of missing data bursts.
The objective of MTSI is to impute missing values in a sequence of input data. More formally, given a graph sequence G[t,t+T ] of length T , we can define the missing data reconstruction error as
L ( X̂[t,t+T ], X̃[t,t+T ],M [t,t+T ] ) = ∑t+T
h=t ∑Nt i=1
〈mih, ` ( x̂ih, x̃ i h ) 〉
〈mih,mih〉 , (1)
where x̂ih is the reconstructed x̃ i h; M [t,t+T ] and m i h are respectively the logical binary complement of M[t,t+T ] and mih, `( · , · ) is an element-wise error function (e.g., absolute or squared error) and 〈 · , · 〉 indicates the standard dot product. Note that, in practice, it is impossible to have access to X̃[t,t+T ] and, as a consequence, it is necessary to define a surrogate optimization objective by, for example, using a forecasting loss or generating synthetic missing values. In the context of trainable, parametric, imputation methods, we consider two different operational settings. In the first one, named in-sample imputation, the model is trained to reconstruct missing values in a given fixed input sequence X[t,t+T ], i.e., the model is trained on all the available data except those that are missing and those that have been removed from the sequence to emulate additional failures for evaluation. Differently, in the second one (referred to as out-of-sample imputation), the model is trained and evaluated on disjoint sequences. Note that in both cases the model does not have access to the ground-truth data used for the final evaluation. The first operational setting simulates the case where a practitioner fits the model directly on the sequence to fill up its gaps. The second, instead, simulates the case where one wishes to use a model fitted on a set of historical data to impute missing values in an unseen target sequence.

heading: 4 GRAPH RECURRENT IMPUTATION NETWORK
text: In this section, we present our approach, the Graph Recurrent Imputation Network (GRIN), a graphbased, recurrent neural architecture for MTSI. Given a multivariate time series X[t,t+T ] with mask M[t,t+T ], our objective is to reconstruct missing values in the input sequence by combining the information coming from both the temporal and spatial dimensions. To do so, we design a novel bidirectional graph recurrent neural network which progressively processes the input sequence both forward and backward in time by performing two stages of imputation for each direction. Then, a feed-forward network takes as input the representation learned by the forward and backward models and performs a final – refined – imputation for each node of the graph and step of the sequence. More precisely, the final imputation depends on the output of two GRIN modules whose learned representations are finally processed (space and time wise) by a last decoding multilayer perceptron (MLP). An overview of the complete architecture is given in Fig. 2. As shown in the figure, the two modules impute missing values iteratively, using at each time step previously imputed values as input. We proceed by first describing in detail the unidirectional model, and then we provide the bidirectional extension.
Unidirectional model Each GRIN module is composed of two blocks, a spatio-temporal encoder and a spatial decoder, which process the input sequence of graphs in two stages. The spatio-temporal encoder maps the input sequence X[t,t+T ] to a spatio-temporal representation H[t,t+T ] ∈ RNt×l by
exploiting an ad-hoc designed recurrent GNN. The spatial decoder, instead, takes advantage of the learned representations to perform two consecutive rounds of imputation. A first-stage imputation is obtained from the representation by using a linear readout; the second one exploits available relational, spatial, information at time step t. In particular, the decoder is implemented by an MPNN which learns to infer the observed values at each i-th node – xit – by refining first-stage imputations considering – locally – Ht−1 and values observed at neighboring nodes.
Spatio-temporal Encoder In the encoder, the input sequence X[t,t+T ] and mask M[t,t+T ] are processed sequentially one step at a time, by means of a recurrent neural network with gates implemented by message-passing layers. Any message-passing operator could be used in principle. In particular, given zit,k−1, i.e., the node features vector at layer k− 1, we consider the general class of MPNNs described as
MPNNk(zit,k−1,Wt) = γk ( zit,k−1, ∑ j∈N (i) ρk ( zit,k−1,z j t,k−1 )) = zit,k, (2)
whereN (i) is the set of neighbors of the i-th node in Gt, γk and ρk are generic, differentiable, update and message functions (e.g., MLPs), and Σ is a permutation invariant, differentiable aggregation function (e.g., sum or mean). Note that several definitions of neighborhood are possible, e.g., one might consider nodes connected by paths up to a certain length l. For the sake of simplicity, from now on, we indicate with MPNN(zit,Wt) the forward pass of a generic K-layered message-passing neural network. In the following, we use MPNNs as the building blocks for our spatio-temporal feature extractors. To learn the dynamics of the system, we leverage on gated recurrent units (GRUs; Cho et al., 2014). As previously mentioned, similarly to Seo et al. (2018) and Li et al. (2018), we implement the GRU gates by relying on the message-passing layers defined above. At the node level, the elements of the message-passing GRU (MPGRU) can be described as:
rit = σ ( MPNN ([ x̂ i(2) t ||mit||hit−1 ] ,Wt )) (3)
uit = σ ( MPNN ([ x̂ i(2) t ||mit||hit−1 ] ,Wt )) (4)
cit = tanh ( MPNN ([ x̂ i(2) t ||mit||rit hit−1 ] ,Wt )) (5)
hit = u i t hit−1 + (1− uit) cit (6)
where rit,u i t are the reset and update gates, respectively, h i t is the hidden representation of the i-th node at time t, and x̂i(2)t is the output of the decoding block at the previous time-step (see next paragraph). The symbols and || denote the Hadamard product and the concatenation operator, respectively. The initial representation Ht−1 can either be initialized as a constant or with a learnable embedding. Note that for the steps where input data are missing, the encoder is fed with predictions from the decoder block, as explained in the next subsection. By carrying out the above computation time and node wise, we get the encoded sequence H[t,t+T ].
Spatial Decoder As a first decoding step, we generate one-step-ahead predictions from the hidden representations of the MPGRU by means of a linear readout
Ŷ (1) t = Ht−1Vh + bh, (7)
where Vh ∈ Rl×d is a learnable weight matrix and bh ∈ Rd is a learnable bias vector. We then define the filler operator as
Φ(Yt) = Mt Xt + M t Yt; (8) intuitively, the filler operator replaces the missing values in the input Xt with the values at the same positions in Yt. By feeding Ŷ (1) t to the filler operator, we get the first-stage imputation X̂ (1) t such that the output is Xt with missing values replaced by the one-step-ahead predictions Ŷ (1) t . The resulting node-level predictions are then concatenated to the mask Mt and the hidden representation Ht−1, and processed by a final one-layer MPNN which computes for each node an imputation representation sit as
sit = γ ( hit−1, ∑ j∈N (i)/i ρ ([ Φ(x̂ j(1) t )||h j t−1||m j t ])) . (9)
Notice that, as previously highlighted, the imputation representations only depend on messages received from neighboring nodes and the representation at the previous step. In fact, by aggregating only messages from the one-hop neighborhood, the representations sit are independent of the input features xit of the i-th node itself. This constraint forces the model to learn how to reconstruct a target input by taking into account spatial dependencies: this has a regularizing effect since the model is constrained to focus on local information. Afterward, we concatenate imputation representation St with hidden representation Ht−1, and generate second-stage imputations by using a second linear readout and applying the filler operator:
Ŷ (2) t = [St||Ht−1]Vs + bs; X̂(2)t = Φ(Ŷ (2)t ) (10)
Finally, we feed X̂(2)t as input to the MPGRU (Eq. 3 – 6) to update the hidden representation and proceed to process the next input graph Gt+1.
Bidirectional Model Extending GRIN to account for both forward and backward dynamics is straightforward and can be achieved by duplicating the architecture described in the two previous paragraphs. The first module will process the sequence in the forward direction (from the beginning of the sequence towards its end), while the second one in the other way around. The final imputation is then obtained with an MLP aggregating representations extracted by the two modules:
ŷit = MLP ([ si,fwdt ||hi,fwdt−1 ||si,bwdt ||hi,bwdt+1 ]) , (11)
where fwd and bwd denote the forward and backward modules, respectively. The final output can then be easily obtained as X̂[t,t+T ] = Φ(Ŷ[t,t+T ]). Note that, by construction, our model can exploit all the available relevant spatio-temporal information, since the only value explicitly masked out for each node is xit. At the same time, it is important to realize that our model does not merely reconstruct the input as an autoencoder, but it is specifically tailored for the imputation task due to its inductive biases. The model is trained by minimizing the reconstruction error of all imputation stages in both directions (see Appendix A).

heading: 5 EMPIRICAL EVALUATION
text: In this section, we empirically evaluate our approach against state-of-the-art baselines on four datasets coming from three relevant application domains. Our approach, remarkably, achieves stateof-the-art performance on all of them.
• Air Quality (AQI): dataset of recordings of several air quality indices from 437 monitoring stations spread across 43 Chinese cities. We consider only the PM2.5 pollutant. Prior works on imputation (Yi et al., 2016; Cao et al., 2018) consider a reduced version of this dataset, including only 36 sensors (AQI-36 in the following). We evaluate our model on both datasets. We use as adjacency matrix a thresholded Gaussian kernel (Shuman et al., 2013) computed from pairwise geographic distances.

heading: D M MAE MSE MRE (%) MAE MSE MRE (%)
text: For the air quality datasets, we adopt the same evaluation protocol of previous works (Yi et al., 2016; Cao et al., 2018) and we show results for both the in-sample and out-of-sample settings. For the traffic and energy consumption datasets, we consider only the out-of-sample scenario (except for matrix factorization which only works in-sample). We simulate the presence of missing data by considering 2 different settings: 1) Block missing, i.e, at each step, for each sensor, we randomly drop 5% of the available data and, in addition, we simulate a failure with probability pfailure = 0.15% and sample its duration uniformly in the interval [min steps,max steps], where min steps and max steps are the number of time steps corresponding respectively to 1 and 4 hours in the traffic case and 2 hours and 2 days for CER-E; 2) Point missing, i.e., we simply randomly mask out 25% of the available data. We split all the datasets into training/validation/test sets. We use as performance metrics the mean absolute error (MAE), mean squared error (MSE) and mean relative error (MRE; Cao et al., 2018) computed over the imputation window. For all the experiments, we use as messagepassing operator the diffusion convolution introduced by Atwood & Towsley (2016). We consider BRITS (Cao et al., 2018) as the principal competing alternative among non-adversarial deep autoregressive approaches, as it shares architectural similarities with our methods. As additional baselines we consider: 1) MEAN, i.e., imputation using the node-level average; 2) KNN, i.e., imputation by averaging values of the k = 10 neighboring nodes with the highest weight in the adjacency matrix Wt; 3) MICE (White et al., 2011), limiting the maximum number of iterations to 100 and the number of nearest features to 10; 4) Matrix Factorization (MF) with rank = 10; 5) VAR, i.e., a Vector Autoregressive one-step-ahead predictor; 6) rGAIN, i.e., an unsupervised version of SSGAN (Miao et al., 2021) which can be seen as GAIN (Yoon et al., 2018a) with bidirectional recurrent encoder and decoder; 7) MPGRU, a one-step-ahead GNN-based predictor similar to DCRNN (Li et al., 2018).
We provide further comment and in depth details on baselines and datasets, together with additional experiments on synthetic data in the appendix.

heading: 5.1 RESULTS
text: Empirical results show that GRIN can achieve large improvements in imputation performance on several scenarios, as well as increased flexibility. In fact, differently from the other state-of-theart baselines, GRIN can handle input with a variable number of dimensions. Tab. 1 shows the experimental results on the air quality datasets. In the in-sample settings, we compute metrics using as imputation the value obtained by averaging predictions over all the overlapping windows; in the out-of-sample settings, instead, we simply report results by averaging the error over windows. GRIN largely outperforms other baselines on both settings. In particular, in the latter case, GRIN decreases MAE w.r.t. the closest baseline by more than 20% in AQI. Interestingly, GRIN consistently outperforms BRITS in imputing missing values also for sensors corresponding to isolated (disconnected) nodes, i.e., nodes corresponding to stations more than 40 km away from any other station (see B.1): this is empirical evidence of the positive regularizations encoded into GRIN. Our method achieves more accurate imputation also in the 36-dimensional dataset, where we could expect the graph representation to have a lower impact. Results for the traffic and smart grids datasets are shown in Tab. 2. In the traffic dataset, our method outperforms both BRITS and rGAIN by a wide margin in all the considered settings while using a much lower number of parameters (see A). In the traffic datasets, on average, GRIN reduces MAE by ≈ 29% w.r.t. BRITS and, in particular, in the Point missing setting of the PEMS-BAY dataset, the error is halved. In CER-E, GRIN consistently outperforms other baselines. Besides show-
ing the effectiveness of our approach in a relevant application field, this experiment also goes to show that GRIN can be exploited in settings where relational information is not readily available.
Finally, Tab. 3 show results – in terms of MAE – of an ablation study on the out-of-sample scenario in AQI, METR-LA (in the Block Missing settings), and CER-E (Point Missing setting). In particular, we compare GRIN against 3 baselines to assess the impact of the spatial decoder and of the bidirectional architecture. The first baseline is essentially a bidirectional MPGRU
where values are imputed by a final MLP taking as inputs hfwdt−1 and h bwd t+1 , while the second one has an analogous architecture, but uses hidden representation and time step t (for both directions) and, thus, behaves similarly to a denoising autoencoder. As reference, we report the results of the unidirectional MPGRU. Results show that the components we introduce do contribute to significantly reduce imputation error. It is clear that spatial decoding and the bidirectional architecture are important to obtain accurate missing data reconstruction, especially in realistic settings with blocks of missing data. Interestingly, the denoising model suffers in the Block missing scenario, while, as one might expect, works well in the Point Missing setting. For additional results and discussion about scalability issues, we refer to the appendix of the paper.
5.2 VIRTUAL SENSING
As a final experiment, we provide a quantitative and qualitative assessment of the proposed method in virtual sensing. The idea (often studied in the context of kriging – see Section 2) is to simulate the presence of a sensor by adding a node with no available data and, then, let the model reconstruct the corresponding time series. Note that for the approach to work several assumptions are needed: 1) we have to assume that the physical quantity being monitored can be reconstructed from observations at neighboring sensors; 2) we should assume a high-degree of homogeneity of sensors (e.g., in the case of air quality stations we should assume that sensors are placed at the same height) or that the features characterizing each neighboring sen-
sor (e.g., placement) are available to the model. In this context, it is worth noting that, due to the inductive biases embedded in the model, GRIN performs reconstruction not only by minimizing reconstruction error at the single node, but by regularizing the reconstructed value for imputation at neighboring sensors. We masked out observed values of the two nodes of AQI-36 with highest (station no. 1014) and lowest (no. 1031) connectivity, and trained GRIN on the remaining part of the data as usual. Results, in Fig. 3, qualitatively show that GRIN can infer the trend and scale for unseen sensors. In terms of MAE, GRIN scored 11.74 for sensor 1014 and 20.00 for sensor 1031 (averages over 5 independent runs).

heading: 6 CONCLUSIONS
text: We presented GRIN, a novel approach for MTSI exploiting modern graph neural networks. Our method imputes missing data by leveraging the relational information characterizing the underlying network of sensors and the functional dependencies among them. Compared against state-of-the-art baselines, our framework offers higher flexibility and achieves better reconstruction accuracy on all the considered scenarios. There are several possible directions for future works. From a theoretical perspective, it would be interesting to study the properties that would guarantee an accurate reconstruction. Furthermore, future work should study extensions able to deal with a non-stationary setting and further assess applications of GRIN in virtual and active sensing.

heading: REPRODUCIBILITY STATEMENT
text: Code to reproduce experiments presented in the paper is provided as supplementary material together with configuration files to replicate reported results. All datasets, except CER-E, are open and downloading links are provided in the supplementary material. The CER-E dataset can be obtained free of charge for research purposes (see appendix). For experiments where failures are simulated, we use random number generators with fixed seed for missing data generation to ensure reproducibility and consistency among experiments and baselines.

heading: ACKNOWLEDGMENTS
text: This research is funded by the Swiss National Science Foundation project 200021 172671: “ALPSFORT: A Learning graPh-baSed framework FOr cybeR-physical sysTems”. The authors wish to thank the Institute of Computational Science at USI for granting access to computational resources.

heading: APPENDIX

heading: A DETAILED EXPERIMENTAL SETTINGS
text: In this appendix, we give more details on the experimental settings used to evaluate our approach. We train all the models by sampling at random 160 batches of 32 elements for each epoch, we fix the maximum number of epochs to 300 and we use early stopping on the validation set with a patience of 40 epochs. All methods are trained using a cosine learning rate scheduler with initial value of 0.001, decayed over the 300 training epochs. During training, we randomly mask out an additional 5% of the input data for each batch to foster robustness to noise and missing data.
For GRIN we minimize the following loss function is
L = L ( Y[t,t+T ],X[t,t+T ],M[t,t+T ] ) + L ( Y
(1),fwd [t,t+T ] ,X[t,t+T ],M[t,t+T ]
) + L ( Y
(2),fwd [t,t+T ] ,X[t,t+T ],M[t,t+T ] ) + L ( Y
(1),bwd [t,t+T ] ,X[t,t+T ],M[t,t+T ]
) + L ( Y
(2),bwd [t,t+T ] ,X[t,t+T ],M[t,t+T ]
) ,
where each L ( · , · , · ) is of the form of Eq. 1 and the element-wise error function is MAE. Note that here we are using X[t,t+T ] and M[t,t+T ] instead of X̂[t,t+T ] and M [t,t+T ].
For BRITS, we use the same network hyperparameters of Cao et al. (2018) for the AQI-36 dataset. To account for the larger input dimension, for the other datasets we increase the number of hidden neurons in the RNNs cells to 128 for AQI/METR-LA and 256 for PEMS-BAY/CER-E. The number of neurons was tuned on the validation sets. For rGAIN we use the same number of units in the cells of the bidirectional RNN used by BRITS, but we concatenate a random vector (sampled from a uniform distribution) of dimension z = 4 to the input vector in order to model the sampling of the data generating process. To obtain predictions, we average out the outputs of k = 5 forward passes. For VAR we used an order of 5 and trained the model with SGD. Since the VAR model needs past 5 observations to predict the next step, we pad each sequence using the mean for each channel. Here we used a batch size to 64 and a learning rate of 0.0005. The order was selected with a small search in the range [2, 12]: we found out a window size of 5 to be ideal for all the considered datasets. For GRIN we use the same hyperparameters in all the datasets: a hidden dimension of 64 neurons for both the spatio-temporal encoder and the spatial decoder and of 64 neurons for the MLP. We use diffusion convolution as message-passing operation, with a diffusion step k = 2 in the spatio-temporal encoder and k = 1 in the temporal decoder. Note that, due to the architectural differences, the other neural network baselines have a number of parameters that is far higher than GRIN (depending on the considered dataset, up to ≈ 4M against ≈ 200K). For MPGRU we use the same hyperparameters of GRIN (64 units for both the spatio-temporal encoder and the decoder).
For data processing we use the same steps of Li et al. (2018), data are normalized across the feature dimension (which means graph-wise for GRIN and node-wise for BRITS/rGAIN/VAR). Data masked out for evaluation are never used to train any model.
All the models were developed in Python (Van Rossum & Drake, 2009) using the following opensource libraries:
• PyTorch (Paszke et al., 2019);
• numpy (Harris et al., 2020);
• Neptune1 (neptune.ai, 2021);
• scikit-learn (Pedregosa et al., 2011);
• fancyimpute (Rubinsteyn & Feldman).
1https://neptune.ai/
The implementation of the diffusion convolutional operator was adapted from the Graph-WaveNet codebase 2. For the implementation of BRITS, we used the code provided by the authors3. The code to reproduce the experiments of the paper is available online4.

heading: B DATASETS
text: In this appendix, we provide more details on the datasets that we used to run experiments. Tab. 4 shows detailed statistics for graph structure associated with each dataset, while Fig. 4 shows the corresponding adjacency matrices. Tab. 5 shows missing data statistics. In the following subsections, we go deeper into details for each dataset.

heading: B.1 AIR QUALITY
text: Air pollution is nowadays a ubiquitous problem. The Urban Computing project (Zheng et al., 2014; 2015) published several datasets containing real measurements of different indices affecting human life in urban spaces. We consider as benchmark the dataset regarding the air quality index (AQI). The complete dataset contains hourly measurements of six pollutants from 437 air quality monitoring stations, spread over 43 cities in China, over a period of one year (from May 2014 to April 2015). Prior works on imputation (Yi et al., 2016; Cao et al., 2018) considered a reduced version of this dataset, including only 36 sensors (AQI-36). This dataset is particularly interesting as a benchmark for imputation due to the high rate of missing values (25.7% in AQI and 13.2% in AQI-36). Along with Yi et al. (2016), we consider as the test set the months of March, June, September and December. We consider both the in-sample and out-of-sample scenarios. In latter case, we do not consider windows overlapping with any of the test months. We use the same procedure of Yi et al. (2016) to simulate the presence of missing data for evaluation.
2https://github.com/nnzhan/Graph-WaveNet 3https://github.com/caow13/BRITS 4https://github.com/Graph-Machine-Learning-Group/grin
AQI CER-E PEMS-BAY METR-LA
We select windows of data of length T = 24 for AQI and T = 36 for AQI-36 (in line with Cao et al. (2018)). To evaluate the imputation performances, we mask out from the test set and use as groundtruth the value xi,jt if: (1) the value is not missing (m i,j t = 1) and (2) the value is missing at the same hour and day in the following month. Besides air quality readings, the dataset provides geographic coordinates of each monitoring station. To obtain an adjacency matrix from the geographic distances between nodes, we use a thresholded Gaussian kernel (Shuman et al., 2013): the weight wi,jt = w i,j of the edge connecting i-th and j-th node is
wi,j =
{ exp ( −dist(i,j) 2
γ
) dist (i, j) ≤ δ
0 otherwise , (12)
where dist ( · , · ) is the geographical distance operator, γ controls the width of the kernel and δ is the threshold. We set γ to the standard deviation of geographical distances in AQI-36 in both datasets. We set δ so that it corresponds to a distance of ≈ 40 km.

heading: B.2 TRAFFIC
text: The study of traffic networks is key for the development of intelligent transportation systems and a relevant application field for network science. While previous works (Yu et al., 2017; Li et al., 2018; Wu et al., 2019; Shang et al., 2020) have assessed spatio-temporal deep learning methods for the traffic forecasting task, we focus on reconstruction. We use as benchmark the PEMS-BAY and METR-LA datasets from Li et al. (2018). PEMS-BAY contains 6 months of data from 325 traffic sensors in San Francisco Bay Area, while METR-LA contains 4 months of sensor readings from 207 detectors in the Los Angeles County Highway (Jagadish et al., 2014); for both datasets, the sampling rate corresponds to 5 minutes.
We use input sequences of 24 steps, which correspond to 2 hours of data. For adjacency, we use a thresholded Gaussian kernel applied to geographic distances following previous works Wu et al. (2019). We split the data into three folds, using 70% of them for training and the remaining 10% and 20% for validation and testing, respectively.

heading: B.3 SMART GRIDS
text: We consider data from the Irish Commission for Energy Regulation (CER) Smart Metering Project5. We select only the subset of the available smart meters monitoring energy consumption of small and medium-sized enterprises (SMEs), i.e., 485 time series with samples acquired every 30 minutes. Note that access to dataset can be obtained free of charge for research purposes.
We build an adjacency matrix by extracting a k-nearest neighbor graph (with k = 10) from the similarity matrix built by computing the week-wise correntropy (Liu et al., 2007) among time series. As in the traffic case, we use a 70%/10%/20% split for training, validation and testing and use a window size of 24 steps. Data were normalized using standard scaling as in the previous settings, and we did not perform additional preprocessing steps.
5https://www.ucd.ie/issda/data/commissionforenergyregulationcer

heading: C ADDITIONAL RESULTS
text: In this appendix, we show an additional experiment in a controlled environment, comparison against additional baselines, additional ablation studies, and sensitivity analyses.

heading: C.1 SYNTHETIC DATA
text: In this experiment, we test our method on the simulated particle system introduced by Kipf et al. (2018)6. We simulate the trajectories of N = 10 particles in a (10× 10) box with elastic collision. Each particle carries a either positive or negative charge q ∈ {1,−1}. Two particles attract each other if they have opposite sign, otherwise they repel. Interaction forces between two particles are ruled by Coulomb’s law. We collect two datasets, each containing 5000 independent simulations of T = 36 steps each. In the first dataset, particles have the same charge in every simulation. In the second one, we sample the charges uniformly at random at the beginning of every simulation. In both scenarios, the initial location and velocity of the particles are drawn randomly. At each step, we randomly remove blocks of consecutive readings with a probability pfailure = 2.5% and a length sampled uniformly from the interval [4, 9]. Here, a reading consists of the (x, y) coordinates of the particle’s position. We further mask out 2.5% of positions at random. The percentage of values not masked out is≈ 74%. For evaluation purposes, we generate another mask using the same missing data distribution and use the masked values as ground-truth for evaluation. We split dataset in training/validation/test folds using 70%/10%/20% splits, respectively.
We test our method (GRIN) and BRITS in both synthetic datasets. We use 32 units for the hidden layer of BRITS (≈ 25K parameters) and 16 units for both the encoder and decoder of GRIN (≈ 10K parameters). Results are reported in Tab. 6. Both the methods take as input only the particles’ positions, with no information about the charges. As can be seen, consistently with what observed by Kipf et al. (2018), relational representations are impressively effective in this scenario. Our method outperforms the baseline by more than an order of magnitude in terms of MSE. Surprisingly, BRITS is more accurate in the setting with varying charge. Our hypothesis is that the added stochasticity acts as a regularization and forces BRITS to learn a more general model.

heading: C.2 EMPIRICAL COMPARISON AGAINST MATRIX FACTORIZATION WITH SIDE INFORMATION
text: As mentioned in Section 2, several matrix factorization approaches – often studied in the context of recommender systems – can be regularized by considering priors on the spatio-temporal struc-
6https://github.com/ethanfetaya/NRI
ture of the data. Intuitively, spatial regularization is achieved by imposing soft constraints on the smoothness of the interpolated function w.r.t. nodes of an underlying graph (Cai et al., 2010; Rao et al., 2015). Temporal regularization can be obtained by imposing analogous constraints modelling temporal dependencies as – eventually weighted – edges of a graph. In temporal regularized matrix factorization (TRMF; Yu et al., 2016), similarly, coefficients of an autoregressive model are used as temporal regularizer.
Tab. 7 shows a comparison of different matrix factorization approaches on imputation in the air quality datasets (where we considered the in-sample setting in Section 5). For TRMF we used an implementation adapted from the Transdim repository7, while for graph regularized matrix factorization (GMRF) we use a custom implementation of the method proposed by Cai et al. (2010). We fixed the rank to be equal to 10 (as the one used in all the experiments for standard MF) and tuned the regularization coefficients on a validation set. Results do show that introducing spatial and temporal regularization improve w.r.t. vanilla MF; however, deep learning methods – and even linear VAR predictors – achieve far superior reconstruction accuracy here. Arguably, low-rank approximation methods might instead have an edge in a low-data regime: this type of analysis is, however, out of the scope of this work.

heading: C.3 SCALABILITY
text: With reference to a standard bidirectional GRU, using MPNNs to implement the cell’s gates increases the computational complexity by a factor that scales with the number of edges O(E) – if using an efficient sparse implementation – or with the number of nodes squared O(N2). Luckily, this overhead can be amortized as most of the computation can be parallelized. Research on scalable and memory efficient GNNs is a very active field (e.g., Hamilton et al., 2017; Frasca et al., 2020): depending on the task, the designer can opt for massage passing operators that meet the application requirements in terms of performance, time and space constraints.

heading: C.4 ABLATION STUDY
text: Here we provide two different ablation studies, the first one on the architecture of GRIN and the second one on the graph structure.

heading: C.4.1 ARCHITECTURAL ABLATIONS
text: Tab. 8 shows additional results for the ablation study presented in Section 5. Consistently with what we already observed, the spatial decoder and bidirectional architecture improve performance and appear particularly relevant in settings with blocks of missing data.

heading: C.4.2 GRAPH STRUCTURE ABLATIONS
text: Here we study how exploiting the relational structure of the problem affects the accuracy of the reconstruction. In particular, we run two additional experiments on the METR-LA dataset (Block missing settings), where instead of using as adjacency matrix the thresholded kernel in Eq. 12, we use (1) a fully connected graph (W = 1) and (2) a graph with no edges (W = I). To provide node – i.e., sensor – identification, we use learnable embeddings as additional node features. Results are
7https://github.com/xinychen/transdim
shown in Tab. 9, performance for BRITS are reported as reference. It is clear that the constraints posed by the graph structure do have an impact on the accuracy of missing data imputation and, at the same time, that spatial information is relevant for the task.

heading: C.5 SENSITIVITY ANALYSIS
text: Finally, in this subsection we carry out an assessment of performance degradation w.r.t. the amount of missing data. Before discussing results, there are a few remarks that are worth bringing up regarding imputation in highly sparse settings. In the first place, GRIN, as well as a large portion of the state-of-the-art baselines, is an autoregressive model, which means that it might be subject to error accumulation over long time horizons. Furthermore, here, consistently with Section 5, we consider the out-of-sample setting which is particularly challenging in the sparse data regime. That being said, GRIN achieves remarkable performance also in this benchmark.
We train one model each for GRIN and BRITS by randomly masking out 60% of input data for each batch during training, then, we run the models on the test set by using evaluation masks with increasing sparsity (note that this causes a distribution shift in evaluation). For each level of sparsity, evaluation is repeated 5 times by sampling different evaluation masks. Results are reported in Tab. 10 and Fig. 5 shows that GRIN outperforms BRITS in all the considered scenarios.
Table 10: Performance with different amounts of missing data. Results averaged over 5 different evaluation masks in the out-sample setting. (P) indicates the Point Missing setting.

heading: METR-LA (P)
text: % Missing 10 20 30 40 50 60 70 80 90
GRIN 1.87 ± 0.01 1.90 ± 0.00 1.94 ± 0.00 1.98 ± 0.00 2.04 ± 0.00 2.11 ± 0.00 2.22 ± 0.00 2.40 ± 0.00 2.84 ± 0.00 BRITS 2.32 ± 0.01 2.34 ± 0.00 2.36 ± 0.00 2.40 ± 0.00 2.47 ± 0.00 2.57 ± 0.01 2.76 ± 0.00 3.08 ± 0.00 4.02 ± 0.01
10 20 30 40 50 60 70 80 90 Missing data (%)
2.0
2.5
3.0
3.5
4.0
M A
E
GRIN BRITS
Figure 5: The plot shows graphically the results in Tab. 10.

Review summaries:
This paper proposes to leverage GNN that takes available relational information for multivariate time series imputation. Each feature dimension of the time series is regarded as a node in GNN and the message-passing is used to implement the update function of bidirectional RNN. On three benchmarks, this new method achieves better performance than the previous STOA method BRITS.

A novel recurrant graph neural network, for spatio-temporal encoding, coupled with a message passing neural network, for spatial decoding, is proposed for imputing missing values in a multivariate time series; The architecture uses temporal and spatial properties of the data. It seems that the approach of gated recurrent neural networks (and related architectures using gated recurrent units, GRUs) is extended to include the graph topology through neighbourhood message passing. Experiments show impressive results on three benchmark datasets, and comparison with simple and state-of-the-art baselines. A short ablation study excludes the spatial decoder and omits the bi-directional architecture to empricially show that they have a significant impact on the imputation performance. Another short experiment tests the ability of this architecture to replace entire missing sensors, using an air pollution dataset with multiple sensors, results are convincing.

In this paper, the authors propose a GNN based imputation method for Time Series.  They argue that past works for time series imputation only take time into account, and not the spacial relationships between time series. The model (called GRIN) works by encoding the data via a GNN, and performing 2 rounds of imputation. The GNN itself acts as gates in a GRU that accounts for the sequential information. This is repeated for a "forward" and "backward" pass, and a final MLP infers the missing value. Experiments on multiple datasets show that the proposed method outperforms baselines.