CGCNImp: a causal graph convolutional network for multivariate time series imputation

Attributes causality modeling

Determining complex correlation dependencies is a key problem in the process of imputing with multivariate time series data. Here we use the neural Granger causality (Tank et al., 2021; Tank et al., 2018) to model the correlation dependency between the attributes of the multivariate time series. Let x_t ∈ R^d be a d-dimensional stationary time series and assume we have observed the process at n timestamps T = (t₀, t₁, …, t_n−1). The basic idea of neural Granger causality is to gauge the extent to which the past activity of one time series is predictive of another time series. Thus, let h_t ∈ R^d represents the d-dimensional hidden state at time t, which represents the historical context of the time series for predicting a component x_ti. The hidden state at time t + 1 is updated recursively (1)${\displaystyle h_{t+1}=f\left(x_t,h_t\right)}$ where f is some nonlinear function that depends on the particular recurrent architecture. We opted for an LSTM to model the recurrent function f due to its effectiveness at modeling complex time dependencies. The standard LSTM model takes the form: (2)${\displaystyle f_t=\sigma \left(W_f{x}_t+U_f{h}_{t-1}\right)i_t=\sigma \left(W_i{x}_t+U_i{h}_{t-1}\right)o_t=\sigma \left(W_o{x}_t+U_o{h}_{\left(t-1\right)}\right)c_t=f_t\odot c_{t-1}+i_t\odot tanh\left(W_c{x}_t+U_c{h}_{t-1}\right)h_t=o_t\odot tanh\left(c_t\right)}$ where ⊙ denotes component-wise multiplication and i_t, f_t, and o_t represent input, forget and output gates, respectively. These control how each component of the state cell c_t, is updated and then transferred to the hidden state used for prediction h_t.W_f, W_i, W_o, W_c, U_f, U_i, U_o, U_c are the parameters that need to learn by LSTM. The output for series i is given by a linear decoding of the hidden state at time t: (3)${\displaystyle x_{ti}=g_i\left(x_{<t}\right)+e_{ti}}$ where the dependency of g_i on the full past sequence x_<t is due to recursive updates of the hidden state. The LSTM model introduces a second hidden state variable c_t, referred to as the cell state, giving the full set of hidden parameter as (c_t, h_t).

In Eq. (2) the set of input maxtrices, (4)${\displaystyle W={\left({\left(W_f\right)}^T,{\left(W_i\right)}^T,{\left(W_o\right)}^T,{\left(W_c\right)}^T\right)}^T}$ controls how the past time series x_t, influences the forget gates, input gates, output gates, and cell updates, and, consequently, the update of the hidden representation. A group lasso penalty across the columns of W can be selected to indicate which Granger series causes series i during estimation. The loss function for modeling the attribute correlation dependencies is as follows: (5)${\displaystyle L_{NG}={min}_{W,U,W_o}\sum _{t=2}^T{\left(x_{it}-g_i\left(x_{<t}\right)\right)}^2+\lambda \sum _{j=1}^d\left|\right|W|{|}_1}$ where U = (((U_f)^T, (U_i)^T, (U_o)^T, (U_c)^T)^T) . The adjacent matrix A, which is produced by neural Granger causality, is stated as: (6)${\displaystyle A_{ij}=\left|\right|W_{g_j}^i|{|}_F^2}$

This adjacency matrix is used in the graph convolution network component discussed in the next subsection.

Attributes correlation dependency modeling

Convolutional neural networks (CNNs) can derive local correlation features but can only be used in Euclidean space. GCNs, however, are semi-supervised models that can handle arbitrary graph-structured data. As such, they have received widespread attention. GCNs can include spectrum and/or spatial domain convolutions. In this study, we use spectrum domain convolutions. In the Fourier domain, spectral convolutions on graphs are defined as the multiplication of a signal x with a filter g_θ: g_θ ∗ x = U g_θ(U^T x). Here U is the matrix of eigenvectors of the normalized graph Laplacian L = I_N − $D^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}=U$ λ U^T, U^Tx is the graph Fourier transform of x, A ∈ R^d∗d is an adjacency matrix and λ is diagonal matrix of its eigenvalues. In multivariate time series, x can also be a X ∈ R^n∗d, where d refers to the number of features and n refers to the time internals. Given the adjacent matrix A which is produced by neural Granger causality, GCNs can perform the spectrum convolutional operation to capture the correlation characteristics. The GCN model can be expressed as: (7)${\displaystyle H_G=\sigma \left({\tilde{W}}^{-\frac{1}{2}}\tilde{A}{\tilde{W}}^{-\frac{1}{2}}X\theta \right)}$ where $\tilde{A}=A+I_N$ is an adjacent matrix with self-connection structures, I_N is an identity matrix, $\tilde{W}$ is a degree matrix, H_G ∈ R^n∗d is the output of GCN which is the input of the temporal auto-correlation dependency modeling, θ is the parameter of GCN, and σ(⋅) is an activation function used for nonlinear modeling.

Temporal auto-correlation dependency modeling

Obtaining complex temporal auto-correlation dependencies is another key problem with imputation of multivariate time series data. In particular, sometimes the input decay may not fully capture the missing patterns since not all missingness information can be represented in decayed input values. Due to its effectiveness at modeling complex time dependencies, we choose to model the temporal dependencies using an LSTM (Graves, 2012). However, to properly learn the characteristics of the original incomplete time series dataset, we find that the time lag between two consecutive valid observations is always changing due to the nil values. Further, the time lags between observations are very important since they follow an unknown non-uniform distribution. These changeable time lags remind us that the influence of the past observations should decay with time if a variable has been missing for a while.

Thus, a time decay vector α is introduced to control the influence of the past observations. Each value of α should be greater than 0 and smaller than 1 with the larger the δ, the smaller the decay vector. Hence, the time decay vector α is modeled as a combination of δ: (8)${\displaystyle {\alpha }_t=1/e^{max\left(0,W_{\alpha }{\delta }_t+b_{\alpha }\right)}}$ where W_α and b_α are parameters that need to be learned. Once the decay vector has been derived, the hidden state in the LSTM h_t−1 is updated in an element-wise manner by multiplying the decay vector α_t to fit the decayed influence of the past observations. Thus, the update functions of the LSTM are as follows: (9)${\displaystyle h_{t-1}^{{}^{\prime }}={\alpha }_t\odot h_{t-1}{i}_t=\sigma \left(W_i\left[h_{t-1}^{{}^{\prime }};H_{G_t}\right]+b_i\right)f_t=\sigma \left(W_f\left[h_{t-1}^{{}^{\prime }};H_{G_t}\right]+b_f\right)s_t=f_t\odot s_{t-1}+i_t\odot tanh\left(W_s\left[h_{t-1}^{{}^{\prime }};H_{G_t}\right]+b_s\right)o_t=\sigma \left(W_o\left[h_{t-1}^{{}^{\prime }};H_{G_t}\right]+b_o\right)h_t=o_t\odot tanh\left(s_t\right)}$ where W_f, W_i, W_o, W_c, b_f, b_i, b_o, b_s are the parameters that need to be learned by the by LSTM and H_G is the output of the attribute correlation dependency modeling.

Attentive neural networks have recently demonstrated success in a wide range of tasks and, for this reason, we use one here. Let H_L be a matrix consisting of output vectors H_L = [h₁, h₂, …, h_n] ∈ℝ^T×d that the LSTM layer produced, where n is the time series length. The representation β_ij of the attention score is formed by a weighted sum of these output vectors: (10)${\displaystyle {\beta }_{ij}=\frac{exp\left(tanh\left(W\left[h_i|h_j\right]\right)\right)}{\sum _{k=1}^{T}exp\left(tanh\left(W\left[h_k|h_j\right]\right)\right)}}$ (11)${\displaystyle H^{{}^{\prime }}=Atte{n}_L\times H_L}$ where the $Atte{n}_L=\left[\begin{array}{c}\hfill {\beta }_{11},{\beta }_{12},\dots ,{\beta }_{1T}\hfill \\ \hfill {\beta }_{21},{\beta }_{22},\dots ,{\beta }_{2T}\hfill \\ \hfill \cdots \hfill \\ \hfill {\beta }_{n1},{\beta }_{n2},\dots ,{\beta }_{nT}\hfill \\ \hfill \hfill \end{array}\right]$ is the attention score.

In Eq. (11), We reconstruct the missing value by some linear transformation of the hidden state H′ at time t. Hence the reconstruction loss is formulated as: (12)${\displaystyle L_{reg}=\sum _{x\in D}{∥x\otimes m-\hat{x}\otimes m∥}_2}$

x represents the input multivariate time series data, $\hat{x}$ represents the imputed multivariate time series data and m means the masking matrix. The expression in Eq. (12) is the masked reconstruction loss that calculates the squared errors between the original observed data x and the imputed sample. Here, it should be emphasized that when calculating the loss, we only calculate the observed data as previously described in Cao et al. (2018), Luo et al. (2018), Luo et al. (2019) and Liu et al. (2022).

Noise reduction and smoothness imputation

In the past, reconstructions were performed directly, which ignores the noise in the actual sampling process. However, in real-world multivariate time series data, when time series are collected the observations may be contaminated by various types of error or noise. Hence, these imputation values may be unreliable. To ensure the reliability of the imputation results, a total variation reconstruction regularization term is applied to the reconstruction results. The method is based on a smoothing function where neighbors with similar values are used to smooth the time series. When applied to time series data, abrupt changes in trend, spikes, dips and the like can all be fully preserved. Compared to a two-norm smoothing constraint, the total variation reconstruction term can ensure smoothness without losing the dynamic performance of the time series (Boyd & Vandenberghe, 2004). Equation (13) applies the total variation reconstruction term to the reconstruction results. As a result, noise in the original data is reduced and completion accuracy is improved. The reconstruction loss is formulated as: (13)${\displaystyle L_{SL}=\sum _{j=1}^M\sum _{i=1}^N\left|\right|{\hat{x}}_{i+1}^j-{\hat{x}}_i^j|{|}_1}$

where M is the number of time series, that is, the number of variables, and N is the length of each time series.

The total object function of our model is: (14)${\displaystyle L_{loss}=\alpha \ast L_{NG}+\beta \ast L_{reg}+\theta \ast L_{SL}}$ where α, β, θ indicate the weights among different parts of the total loss. We optimize Eq. (14) by Adam optimizer.

KDD CUP 2018 dataset

The KDD dataset comes from the KDD CUP Challenge 2018. The dataset, which is is a public meteorologic dataset, has about 15% missing values. It was collected hourly between 2017/1/20 to 2018/1/30 in Beijing, collecting air quality and weather data. Each record contains 12 attributes, for example CO, weather, temperature etc. In our experiment, we select 11 common features as the comparison methods described in the next section. We split this dataset for every hour. For every 48 h, we randomly drop p percent of the dataset, and then we impute these time series with different models and calculate the imputation accuracy by RMSE and MAE where p ∈ {10, 20, 30, 40, 50, 60, 70, 80, 90}.

Bird migration dataset in China

The Birds Migration Dataset collects migration trace data which comes from the project Strategic Priority Research Program of the Chinese Academy of Sciences. The dataset was collected hourly between 2017/12/30 to 2018/5/10 based on the species Anser fabalis and Anser albifrons. Each record contains 13 attributes which are longitude, latitude, speed height, speed velocity, heading, temperature etc. The dataset is about 10% missing values. We select 10 common features comprising longitude, latitude, speed height, speed velocity, heading, temperature etc. for our experiments. We split this dataset into 5 min time series, and for every 5 min, we randomly drop p percent of the dataset, and then we impute these time series with different models and calculate the imputation accuracy by RMSE and MAE between original values and imputed values where p ∈ {10, 20, 30, 40, 50, 60, 70, 80, 90}.