Year: 2019

Abstract:
Many real-world applications involve multivariate, geo-tagged time series data: at each location, multiple sensors record corresponding measurements. For example, air quality monitoring system records PM2.5, CO, etc. The resulting time-series data often has missing values due to device outages or communication errors. In order to impute the missing values, state-of-the-art methods are built on Recurrent Neural Networks (RNN), which process each time stamp sequentially, prohibiting the direct modeling of the relationship between distant time stamps. Recently, the self-attention mechanism has been proposed for sequence modeling tasks such as machine translation, significantly outperforming RNN because the relationship between each two time stamps can be modeled explicitly. In this paper, we are the first to adapt the self-attention mechanism for multivariate, geo-tagged time series data. In order to jointly capture the self-attention across different dimensions (i.e. time, location and sensor measurements) while keep the size of attention maps reasonable, we propose a novel approach called Cross-Dimensional Self-Attention (CDSA) to process each dimension sequentially, yet in an order-independent manner. On three real-world datasets, including one our newly collected NYCtraffic dataset, extensive experiments demonstrate the superiority of our approach compared to state-of-the-art methods for both imputation and forecasting tasks.

heading: 1 INTRODUCTION
text: Various monitoring applications, such as those for air quality (Zheng et al. (2015)), health-care (Silva et al. (2012)) and traffic (Jagadish et al. (2014)), widely use networked observation stations to record multivariate, geo-tagged time series data. For example, air quality monitoring systems employ a collection of observation stations at different locations; at each location, multiple sensors concurrently record different measurements such as PM2.5 and CO over time. Such time series are important for advanced investigation and also are useful for future forecasting. However, due to unexpected sensor damages or communication errors, missing data is unavoidable. It is very challenging to impute the missing data because of the diversity of the missing patterns: sometimes almost random while sometimes following various characteristics.
Traditional data imputation methods usually suffer from imposing strong statistical assumptions. For example, Scharf & Demeure (1991) and Friedman et al. (2001) fit a smooth curve on observations in either time series (Ansley & Kohn (1984); Shumway & Stoffer (1982)) or spatial distribution (Friedman et al. (2001); Stein (2012)). Deep learning methods (Li et al. (2018); Che et al. (2018); Cao et al. (2018); Luo et al. (2018a)) have been proposed to capture temporal relationship based on RNN (Cho et al. (2014b); Hochreiter & Schmidhuber (1997); Cho et al. (2014a)). However, due to the constraint of sequential computation over time, the training of RNN cannot be parallelized and thus is usually time-consuming. Moreover, the relationship between each two distant time stamps cannot be directly modeled. Recently, the self-attention mechanism as shown in Fig. 1(b) has been proposed by the seminal work of Transformer (Vaswani et al. (2017)) to get rid of the limitation of sequential processing, accelerating the training time substantially and improving the performance significantly on seq-to-seq tasks in Natural Language Processing (NLP) because the relevance between each two time stamps is captured explicitly.
In this paper, we are the first to adapt the self-attention mechanism to impute missing data in multivariate time series, which cover multiple geo-locations and contain multiple measurements as
Under review as a conference paper at ICLR 2020
Completed Value Available ObservationMissing Value
Data Imputation Model (a) Location (L) Time (T) Measurement (M) Location (L) Time (T) Measurement (M) (c)
Temporal Only  (Traditional)
Cross‐Dimensional  (Ours)
(L) (T) (M)
(L)
(T)
(M)
Query
Key
Value
Weighted Sum
(b)
Attention (Q,K,V)Attention Map
shown in Fig. 1(a). In order to impute a missing value in such unique multi-dimensional data, it is very useful to look into available data in different dimensions (i.e. time, location and measurement), as shown in Fig. 1(c), to capture the intra-correlation individually. To this end, we investigate several choices of modeling self-attention across different dimensions. In particular, we propose a novel Cross-Dimensional Self-Attention (CDSA) mechanism to capture the attention crossing all dimension jointly yet in a decomposed manner. In summary, we make the following contributions:
(i) We are the first to apply the self-attention mechanism to the multivariate, geo-tagged time series data imputation task, replacing the conventional RNN-based models to speed up training and directly model the relationship between each two data values in the input data.
(ii) For such unique time series data of multiple dimensions (i.e. time, location, measurement), we comprehensively study several choices of modeling self-attention crossing different dimensions. Our proposed CDSA mechanism models self-attention crossing all dimensions jointly yet in a dimension-wise decomposed way, preventing the size of attention maps from being too large to be tractable. We show that CDSA is independent with the order of processing each dimension.
(iii) We extensively evaluate on two standard benchmarks and our newly collected traffic dataset. Experimental results show that our model outperforms the state-of-the-art models for both data imputation and forecasting tasks. We visualize the learned attention weights which validate the capability of CDSA to capture important cross-dimensional relationships.

heading: 2 RELATED WORK
text: Statistical data imputation methods. Statistical methods (Ansley & Kohn (1984); Zhang (2003); Shumway & Stoffer (1982); Nelwamondo et al. (2007); Buuren & Groothuis-Oudshoorn (2010)) often impose assumptions over data and reconstruct the missed value by fitting a smooth curve to the available values. For instance, Kriging variogram model (Stein (2012)) was proposed to capture the variance in data w.r.t. the geodesic distance. Matrix completion methods (Acuna & Rodriguez (2004); Yu et al. (2016); Friedman et al. (2001); Cai et al. (2010); Ji & Ye (2009); Ma et al. (2011)) usually enforce low-rank constraint.
𝑉”
𝑉
”
RNN-based data imputation methods. Li et al. (2018) proposed DCGRU for seq-to-seq by adopting graph convolution (Chung & Graham (1997); Shi (2009); Shuman et al. (2012)) to model spatial-temporal relationship. Luo et al. (2018a) built GRUI by incorporating RNN into a Generative Adversarial Network (GAN). Nevertheless, the spatiotemporal and measurements correlation are mixed and indistinguishable. so that the mediate back propagation from loss of available observation can contribute to the missing value updating. Nevertheless, these RNN-based models fundamentally suffer from the constraint of sequential processing, which leads to long training time and prohibits the direct modeling of the relationship between two distant data values.
Self-attention. Recently, Vaswani et al. (2017) introduced the Transformer framework which relies on self-attention, learning the association between each two words in a sentence. Then self-attention has been widely applied in seq-to-seq tasks such as machine translation, image generation (Yang et al. (2016); Zhang et al. (2018a)) and graph-structured data (Veličković et al. (2017)). In this paper, we are the first to apply self-attention for multi-dimensional data imputation and specifically we investigate several choices of modeling self-attention crossing different data dimensions.

heading: 3 APPROACH
text: In Sec. 3.1, we first review the conventional self-attention mechanism in NLP. In Sec. 3.2, we propose three methods for computing attention map cross different dimension. In Sec. 3.3 and 3.4, we present details of using CDSA for missing data imputation.

heading: 3.1 CONVENTIONAL SELF-ATTENTION
text: As shown in Fig. 1(b), for language translation task in NLP, given an input sentence, each word xi is mapped into a Query vector qi of d-dim, a Key vector ki of d-dim, and a Value vector vi of v-dim. The attention from word xj to word xi is effectively the scaled dot-product of qi and kj after Softmax, which is defined as A(i, j) = exp(S(i, j)) (∑T j=1 exp(S(q, j)) )−1 where Then, vi is updated to v′i as a weighted sum of all the Value vectors, defined as v ′ i = ∑T j=1 A(i, j)vj , after which each v′i is mapped to the layer output x ′ i of the same size as xi. In order to adapt the self-attention from NLP to our multivariate, geo-tagged time series data, a straightforward way is to view all data in a time stamp as one word embedding and model the self-attention over time.

heading: 3.2 CROSS-DIMENSIONAL SELF-ATTENTION
text: In order to model Cross-Dimensional Self-Attention (CDSA), in this section we propose three solutions: (1) model attention within each dimension independently and perform late fusion; (2)
model attention crossing all dimension jointly; (3) model attention crossing all dimension in a joint yet decomposed manner. We assume the input X ∈ RT×L×M has three dimensions corresponding time, location, measurement. X can be reshaped into 2-D matrices (i.e. XT ∈ RT×LM , XL ∈ RL×MT , XM ∈ RM×TL) or an 1-D vector (i.e. X ∈ RTLM×1). Similarly, this subscript may be applied on the Query, Key and Value, e.g., Q ∈ RT×L×M×d, QL ∈ RL×MTd and Q ∈ RTLM×d.

heading: 3.2.1 INDEPENDENT
text: As shown in Fig. 2(a), the input X is reshaped into three input matrices XT , XL and XM. Three streams of self-attention layers are built to process each input matrix in parallel. Such as the first layer in stream on XL, each vector XL(l, :) of MT -dim is viewed as a word vector in NLP. Following the steps in Sec. 3.1, XL(l, :) is mapped to QL(l, :) and KL(l, :) of dL-dim, as well as VL(l, :) of vL-dim. The output of every stream’s last layer are fused through element-wise addition, X ′ = αTX ′T + αLX ′L + αMX ′M , where the weights αT , αL and αM are trainable parameters. Besides, the hyper-parameters for each stream such as the number of layers, are set separately.

heading: 3.2.2 JOINT
text: As shown in Fig. 2(b), the three-dimensional input X is reshaped as to X . Each unit X(p) is mapped to Q(p, :) and K(p, :) of d-dim as well as V (p, :) of v-dim, where p = p(t, l,m) denotes the index mapping from the 3-D cube to the vector form. In this way, an attention map of dimension TLM × TLM is built to directly model the cross-dimensional interconnection.

heading: 3.2.3 DECOMPOSED
text: The Independent manner sets multiple attention sub-layers in each stream to model the dimensionspecific attention but fail in modeling cross-dimensional dependency. In contrast, the Joint manner learns the cross-dimensional relationship between units directly but results in huge computation workload. To capture both the dimension-specific and cross-dimensional attention in a distinguishable way, we propose a novel Decomposed manner.
As shown in Fig. 2(c), the input X is reshaped as input matrices XT , XL, XM and X . Each unit X(p) is mapped to vector V (p, :) of v-dim as in the Joint while XT , XL and XM are used for building attention map AT , AL,AM individually as in the Independent. The attention maps are applied on Value vector in order as,
V ′ = AV = ÃMVS,2 = ÃMÃLVS,1 = ÃMÃLÃTV . (1)
The attention map with ˜ is reshaped from the original attention map and consistent with the calculation in (1), e.g., ÃT ∈ RTLM×TLM is reshaped from AT ∈ RT×T . More specifically,
ÃT = AT ⊗ IL ⊗ IM ,
ÃL = IT ⊗AL ⊗ IM ,
ÃM = IT ⊗ IL ⊗AM ,
(2)
where ⊗ denotes tensor product and I is the Identity matrix where the subscript indicates the size, e.g., IT ∈ RT×T . Although the three reshaped attention maps are applied with a certain order, according to (2), we show that each unit in A is effectively calculated as
A(p0, p1) = AT (t0, t1)AL(l0, l1)AM (m0,m1), (3)
where p0 = p(t0, l0,m0), p1 = p(t1, l1,m1). Following the associativity of tensor product, we demonstrate Ãσ(T )Ãσ(L)Ãσ(M) = AT ⊗AL ⊗AM , (4) where σ = σ(T,L,M) denotes the arbitrary arrangement of sequence (T,L,M), e.g., σ =(L,T,M) and σ(T) = L. Effectively, the arrangement σ is the order of attention maps to update V . As (3)- (4) shows that the weight in A is decomposed as the product of weights in three dimension-specific attention maps, the output and gradient back propagation are order-independent. Furthermore, we show in Supp. B that the cross-dimensional attention map has the following property:
TLM∑ p1=1 A(p0, p1) = T∑ t1=1 L∑ l1=1 M∑ m1=1 AT (t0, t1)AL(l0, l1)AM (m0,m1) = 1. (5)
In summary, the Independent builds attention stream for each dimension while the Joint directly model the attention map among all the units. Our proposed CDSA is based on the Decomposed, which forms a cross-dimensional attention map, out of three dimension-specific maps. As an alternative of the Decomposed, the Shared maps unit X(p) to Q(p, :) and K(p, :) of d-dim and calculates all three dimension-specific attention map, e.g., AL = Softmax(QLK>L / √ MTd). As shown in Table 1, by using Tensorflow profile and fixing the hyper-parameters with detailed explanations in Supp., the Decomposed significantly decreases the FLoating point OPerations (FLOPs) compared to the Joint and requires less variables than the Independent. Detailed comparisons are reported in Sec. 4.3.

heading: 3.3 FRAMEWORK
text: Imputation: As shown in Fig. 3(a), we apply our CDSA mechanism in a Transformer Encoder, a stack of N = 8 identical layers with residual connection (He et al. (2016)) and normalization (Lei Ba et al. (2016)) as employed by Vaswani et al. (2017). To reconstruct the missing (along with other) values of the input, we apply a fully connected Feed Forward network on the final Value tensor, which is trained jointly with the rest of the model.
Forecasting: As shown in Fig. 3(b), we apply our CDSA mechanism in Transformer framework where we set N = 9 for both encoder and decoder. Similar to imputation, we use a fully connected feed forward network to generate the predicted values.

heading: 3.4 IMPLEMENTATION DETAILS
text: We normalize each measurement of the input by subtracting the mean and dividing by standard deviation across training data. Then the entries with missed value are set 0. We use the Adam optimizer (Kingma & Ba (2014)) to minimize the Root Mean Square Error (RMSE) between the prediction and ground truth. The model is trained on a single NVIDIA GTX 1080 Ti GPU. More details (e.g., network hyper-parameters, learning rate and batch size) can be found in Supp.

heading: 4 EXPERIMENTS

heading: 4.1 DATASETS, TASKS, EVALUATION METRICS
text: NYC-Traffic. New York City Department of Transportation has set up various street cameras1. Each camera keeps taking a snapshot every a few seconds. The is collected around 1-month data
1https://nyctmc.org/
for 186 cameras on Manhattan from 12/03/2015 to 12/26/2015. For each snapshot, we apply our trained faster-RCNN (Ren et al. (2015)) vehicle detection model to detect the number of vehicles (#vehicle) in each snapshot. To aggregate such raw data into time series, for every non-overlapping 5-minute window, we averaged #vehicle from each snapshot to obtain the average #vehicle as the only measurement. Finally, we obtained 186 time series and the gap between two consecutive time stamps is 5 minutes.
The natural missing rate of the whole dataset is 8.43%. In order to simulate experiments for imputation, we further remove some entries and hold them as ground truth for evaluation. The imputation task is to estimate values of these removed entries. To mimic the natural data missing pattern, we model our manual removal as a Burst Loss, which means at certain location the data is continuously missed for a certain period. More details about vehicle detection and burst loss are be found in Supp. To simulate various data missing extents, we vary the final missing rate after removal from 20% to 90%. For each missing rate, we randomly select 432 consecutive time slots to train our model and evaluate the average RMSE of 5 trials. The dataset will be released publicly.
KDD-2018 (Cup (2018)) is an Air Quality and Meteorology dataset recorded hourly. As indicated in Luo et al. (2018a), 11 locations and 12 measurements are selected. The natural missing rate is 6.83%. In order to simulate experiments for imputation, we follow Luo et al. (2018a) to split the data to every 48 hours, randomly hold values of some available entries and vary the missing rate from 20% to 90%. Mean Squared Error (MSE) is used for evaluation.
METR-LA (Jagadish et al. (2014)). We follow Li et al. (2018) to use this dataset for traffic speed forecasting. This dataset contains traffic speed at 207 locations recorded every 5 minutes for 4 months. Following Li et al. (2018), 80% of data at the beginning of these 4 months is used for training and the remaining 20% is for testing. In order to simulate the forecasting scenario, within either training or testing set, every time series of consecutive 2 hours are enumerated. For each time series, data in the first hour is treated as input and data in the second hour is to be predicted. We respectively evaluate the forecasting results at 15-th, 30-th, 60-th minutes in the second 1 hour and also evalaute the average evaluation results within the total 1 hour. We use RMSE, Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE) as evaluation metrics.

heading: 4.2 COMPARISONS WITH STATE-OF-THE-ART
text: Imputation (NYC-Traffic, KDD-2018) In Table 2 , our CDSA consistently outperforms traditional methods (i.e., Auto Regressive, Kriging expo, Kriging linear) and recent RNN-based methods (i.e. MTSI, BRITS, DCRNN) over a wide range of missing rate. Because CDSA leverages the self-attention mechanism to avoid sequential processing of RNN and directly model the relationship
between distant data. Table 3 shows that our method again achieves significant improvements on cross-dimensional data imputation task. Detailed overview of baselines can be found in Supp.
Forecasting (METR-LA). Table 4 shows that for the forecasting task, our CDSA method outperforms previous methods in most cases. In particular, our method demonstrates clear improvement at long-term forecasting such as 60 min. This again confirms that our CDSA cthe effectiveness of directly modeling the relationship between every two data values (could from different dimensions and of far distance). But RNN-based methods and methods that sequentially conduct spatial conv and temporal conv fail to model the distant spatio-temporal relationship explicitly.

heading: 4.3 DISCUSSIONS
text: The effects of different training losses: For the forecasting task in METR-LA, we compare the performance by setting different training loss in Table 5 and we can see the performance with RMSE as loss metric achieves the best performance.
Ablation study of different cross-dimensional self-attention manners: We compare the performance for different solutions in CDSA mechanism on the three datasets listed above. 1) The way of attention modeling determines the computational complexity. As shown in Table 1, since the Independent calculates dimension-specific Value vectors in parallel, the number of variables and FLOPs are larger than those of the Decomposed. As the Joint and the Shared both share the variables for each dimension, the number of variables is small and basically equals with each other. As the Joint builds a huge attention map, its FLOPs is much larger than others. Since the Decomposed draws attention maps like the Independent but shares Value like the Joint, it reduces the computational complexity significantly. 2) As shown in Table 6 - 8, we evaluate these methods on three datasets and the Decomposed always achieves the best performance thanks to the better learning ability compared to the Joint and Shared. More discussions can be found in Supp.
Study of using the imputed time series for forecasting. On NYC-Traffic of missing rate 50%, we impute missing values in historical data (using statistical methods and our CDSA respectively) and
then feed the imputed historical data into traffic prediction model ARIMA. We compare performances in terms of RMSE: Mean Filling 1.953, Kriging expo 1.681, Kriging linear 1.733, MTSI 1.675, DCRNN 1.666, BRITS 1.579, CDSA 1.536. This indicate that when using the imputed time series for forecasting, our CDSA can achieve significant gains in the downstream forecasting task as well compared to traditional imputation methods. More details can be found in Supp. E.3.
Table 7: Comparisons of different manners to implement CDSA on dataset KDD-2018.
Model \Missing Rate 20% 30% 40% 50% 60% 70% 80% 90% CDSA(Independent) 0.482 0.523 0.351 0.366 0.484 0.573 .608 0.721
CDSA(Joint) 0.451 0.497 0.317 0.336 0.404 0.520 0.558 0.677 CDSA(Shared) 0.783 0.799 0.672 0.692 0.784 0.793 0.791 0.832 CDSA(Decomposed) 0.373 0.393 0.287 0.291 0.387 0.495 0.521 0.631
Table 8: Comparisons of different manners to implement CDSA on dataset METR-LA.
Model \ Dataset 60 min MeanMAE RMSE MAPE MAE RMSE MAPE CDSA(Independent) 3.54 7.02 10.29% 3.25 6.29 8.81%
CDSA(Joint) 3.63 7.62 10.54% 3.30 6.83 9.43% CDSA(Shared) 3.92 7.93 11.17% 3.53 7.33 10.26% CDSA(Decomposed) 3.40 6.27 9.76% 3.16 5.48 8.50%
Attention Map Visualization: Fig. 4 shows an PM10 imputation example in location fangshan at t2. Since the pattern of PM2.5 around t2 is similar to that at t1, the attention in orange box is high. As we can see that PM2.5 and PM10 are strongly correlated , in order to impute PM10 at t2, our model utilizes PM10 at t1 (green arrow) and PM2.5 at t1 (blue arrow), which crosses dimensions. More visualization examples can be found in Supp.
PM2.5 PM10
PM2.5
PM10
4t1t 2t 3t 4t1t 2t 3t
4t
1t 2t 3t
4t
1t 2t 3t

heading: PM2.5
text: PM2.5 PM10
PM2.5
PM10
4t
1t 2t 3t
4t
1t 2t 3t

heading: 5 CONCLUSION
text: In this paper, we have proposed a cross-dimensional self-attention mechanism to impute the missing values in multivariate, geo-tagged time series data. We have proposed and investigated three methods to model the cross-dimensional self-attention. Experiments show that our proposed model achieves superior results to the state-of-the-art methods on both imputation and forecasting tasks. Given the encouraging results, in the future we plan to extend our CDSA mechanism from multivariate, geo-tagged time series to the input that has higher dimension and involves multiple data modalities.

heading: A MODEL ARCHITECTURE

heading: A.1 NORMALIZATION LAYER
text: Under the NLP scenario, each word is embedded as a vector and normalized individually. However, in the Cross-Dimensional scenario, the normalization applied on each individual unit will always lead to a zero-output. As shown in Fig. 5, different measurements may exhibit different correlation, i.e., PM2.5 and PM 10 are significantly positively correlated (ρPM2.5, PM 10 = 0.8278) while NO2 and O3 are negatively correlated (ρNO2,O3 = −0.5117). As discussed in Fig. 4(b) in the paper, different measurement may be used as reference for imputation of other measurements. As such, the normalization cross multiple measurements is unreasonable and we choose to apply normalization for each measurement in parallel which presumes that the time series inside the spatial network is essentially drawn from a standard normal distribution. We subsequently add the trainable scalar and bias to scale the normalized value and the scalars (biases) for different measurements are trained individually.

heading: A.2 UNIT-WISE FEED-FORWARD NETWORK
text: Making use of the approximation property of multi-linear layer Hornik et al. (1989), a fully connected feed-forward network (FFN) is applied to each unit separately and identically. This FFN consists of three fully connected layer while RuLU is set as the activation function.
FFN(x) = max(max(0,xW1 + b1)W2 + b2)W3 + b3 (6)
During experiment, since the FFN is applied on each unit individually, we found the improvement by simply increasing the size of weight and bias of each layer is not obvious while increasing the depth of FFN will lead to obvious improvement.
A.3 IMPUTATION MASK
For the self-attention sub-layer in imputation task NYC and KDD-2018, we modify the attention map with mask in (7) to prevent unit of available observation from contributing to the estimation of itself.
S(i, j) = { −∞ i = j qik > j / √ d otherwise
(7)
where qi and kj are d-dim vectors. This masking, combined with fact that there is no offset between input and output, ensures that the estimation of unit X(p(t, l,m)) depends on all the units except for itself, including both available and complemented units. In this way, the gradient back-propagation can be used to update the missed value effectively. The Table 9 shows the performance improvement of imputation mask applied in our model and demonstrate that mask prevents the estimation of itself and improve the inference ability of the model.
A.4 ATTENTION MAP CALCULATION
Available ObservationMissing Value
Query
Value
(b)
Time (T)
1 1 𝑑
𝑄 𝑝 , : ∈ 𝑅
Joint As shown in Fig. 6, when we build the attention among different units in the Joint, two different kernels will be used to map each input unit X(p) = X (t, l,m) to an 1-D Query vector and an 1-D Key vector individually. As attention map is a scaled dot-product between Query and Key (Fig. 6 Left) after Softmax, each value of attention map in Joint is essentially the scaled numerical multiplication between each two units of input (Fig. 6 Right) after Softmax. As such, the multiple parameters inside that kernels only perform as a single scalar and the learning ability/relationship representation in Joint is limited.
Decomposed According to Sec. 3.2.3, to calculate the dimension-specific attention map, e.g., the attention map of Time AT , the input X will be reshaped into matrix XT . Thus, the units corresponding to the same timestamp, reshaped into one vector XT (t, :), will be mapped into dimension-specific Query vector QT (t, :) and Key vector KT (t, :). As more parameters will be introduced into such vector-vector mapping, each dimensional-specific map can learn the intra-correlation and the crossdimensional attention map can model the relationship among each input units effectively.
Shared Same with the Decomposed, the Shared will calculate the dimension-specific attention maps individually. Like Joint, the Shared will map each unit X(p) into vector. To calculate the dimensionspecific attention map (e.g., AT ), the “vector-vector mapping” is essentially the summation of the units corresponding to the same timestamp while the multiple parameters introduced in the mapping still perform as a single scalar. As a result, the learning ability of the intra-correlation is limited so that the cross-dimensional attention map cannot model the relationship among input units effectively.
Under review as a conference paper at ICLR 2020
B ATTENTION MAP RESHAPE IN THE Decomposed
11t 12t
21t 22t
11t 11t
11t 11t
12t 12t
12t 12t
21t 22t 21t 22t
21t 22t 21t 22t
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11l 12l 21l 22l 11l 12l
21l 22l 11l 12l
21l 22l 11l 12l
21l 22l
11l 12l 21l 22l 11m 12m
21m 22m (a) Attention Map of Time, Location and Measurement
11t 12t
21t 22t
11t 11t
11t 11t
12t 12t
12t 12t
21t 22t 21t 22t
21t 22t 21t 22t
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11l 12l 21l 22l 11l 12l
21l 22l 11l 12l
21l 22l 11l 12l
21l 22l
11l 12l
21l 22l
11m 12m
21m 22m
(b) Time 11t 12t
21t 22t
11t 11t
11t 11t
12t 12t
12t 12t
21t 22t 21t 22t
21t 22t 21t 22t
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11l 12l 21l 22l 11l 12l
21l 22l 11l 12l
21l 22l 11l 12l
21l 22l
11l 12l
21l 22l
1m 12m
1m 22m
(c) Location 11t 12t
21t 22t
11t 11t
11t 11t
12t 12t
12t 12t
21t 22t 21t 22t
21t 22t 21t 22t
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11m 12m 21m 22m
11l 12l 21l 22l 11l 12l
21l 22l 11l 12l
21l 22l 11l 12l
21l 22l
11l 2l
21l 22l
1m 12m
21m 22m
(d) Measurement
Figure 7: Reshape Attention Map: The original attention map in (a) is reshaped into (b-d) respectively. The variable (t, l,m) denotes the attention map units for Time, Location and Measurement where the subscript labels the attention between two units, i.e., t12 = Softmax(q1k>2 / √ d) where q1 and k2 are d-dim vectors. Besides, the empty entry indicates 0.
As described in (2) in the paper, the original attention maps AT , AL and AM are reshaped to A′T , A′L and A ′ M . By setting T = L = M = 2 as an example, we draw the attention maps before and after reshape in Fig. 7. Making use of the matrix structure, we have
TLM∑ p1=1 A(p0, p1) = T∑ t1=1 L∑ l1=1 M∑ m1=1 AT (t0, t1)AL(l0, l1)AM (m0,m1)
= T∑ t1=1
( AT (t0, t1)
L∑ l1=1
( AL(l0, l1)
M∑ m1=1 AM (m0,m1)
)) = 1,
(8)
where ∑T t1=1 AT (t0, t1) = 1, ∑L l1=1 AL(l0, l1) = 1 and ∑M m1=1
AM (m0,m1) = 1. Besides, the reshape operation for attention map on different dimension is only determined by the index mapping function, p = p(t, l,m) = LMt+Ml +m, from the 3-dimensional cube to the vector form.

heading: C CDSA FRAMEWORK FOR DATA FORECASTING
text: Different from time series imputation task, the time series forecasting task use the current observation to estimate the time series in the future. To begin with, We first introduce the framework for time series forecasting and we compare the performance for two different types of input.

heading: C.1 PREDICTION FRAMEWORK
text: As shown in Fig. 8, we apply our CDSA mechanism in Transformer framework and use the same Feed Forward structure as in Sec. 3.3 in the paper. Notably, we set N = 9 for both encoder and decoder and no CDSA module is used to derive a complement input for prediction task where the missing value is replaced with global mean. The architecture detail is shown in Fig. 9 Completed Value Available Observation Missing Value (L) (T)
(M)
Transformer Encoder w/ CDSA(L)
(T)
(M)
(L)
(T)
(M) Loss
(a) Imputation
Transformer Encoder w/ CDSA Transformer Decoder w/ CDSA(L)
(M)
1 T+1
2T
(L)
(M)
Feed Forward
T+1
2T
(L)
(M) Loss
Feed Forward
Figure 8: The framework of using Crossing-Dimensional Self-Attention (CDSA) for data forecasting.
T+1
2T
(L)
(M)
PredictionImputation
(L)
(T)
(M)
(L)
(T)
(M)
1
T
(L)
(M)
Multi-Head CMSA
Feed Forward
Add & Norm
Add & Norm
Inputs
Multi-Head CMSA
Feed Forward
×Ν
Ν
Multi-Head CMSA (Mask L&M)
Feed Forward
Multi-Head CMSA (Mask TLM)
Add & Norm
Add & Norm
Add & Norm
Completed Value
Available Observation
Missing Value
Multi-Head CMSA
Feed Forward
Multi-Head CMSA
Feed Forward
Add & Norm
Add & Norm
Inputs
Outputs
Outputs
Ν
Figure 9: Model Architecture

heading: C.2 DISCUSSION FOR INPUT FORM
text: The decoder in NLP task originally sets the shifted output as input. Take the German-to-English translation scenario as an example where the embedded word vectors of German are set as the Encoder input, the model will first send a [GO] vector into the decoder and generate the first word vector of the translated English sequence, then the predicted vector will be sent into the decoder to predict the next word vector and the decoder will complete the sentence translation by repeating this operation until the end.
Mapping directly from this model setting in NLP task to our series forecasting scenario, we can also use the shifted ground truth as the decoder input, i.e., to forecast the speed of the next T time stamps given the speed of the first T time stamps, the data of T ≤ t ≤ 2T − 1 are sent into the decoder. Consequently, the Casual Mask in Vaswani et al. (2017) need to be modified to make sure that the leftward information flow is prevented.
For data forecasting by CDSA in the Decomposed, the masking on Attention Map on Time is simply masking out (setting to −∞) the values in the input of Softmax which corresponds to the leftward information flow. Same with Vaswani et al. (2017), the masking is only adopted in the Multi-head CDSA layer labeled as (Mask TLM) in Fig. 9. However, as shown in Fig. 10, to calculate the Attention Map of Location and Measurement for data forecasting at 2T − 1, all illegal units of input corresponding to t ≥ 2T − 1 have to be masked out (setting as 0). Then, the Masked Input are mapped to Query, Key and Value to build the Attention Map and calculate the Updated Value.
# Identity # Measurement
Time
H
T
1t 2 t 3t
4t H
T
1t 2 t 3t
4t W
W
3t
Original Input
Masked Out
wid wid Input Attention MapIdentity Sample
(L)
(M)
(L)
(M)
2T‐1
Besides, the decoder generates predictions given previous ground truth observations during training while the ground truth observations are replaced by predictions generated by the model itself during testing. As, the discrepancy between the input distributions of training and testing can cause degraded performance, We adopt the integrated sampling Bengio et al. (2015) as in Li et al. (2018) to mitigate this impact while this method is very time-consuming for the Transformer framework. During testing,
In summary, by setting shifted output as the Decoder input, multiple Attention Map are calculated for forecasting value of different time stamps which requires huge memory usage. Still, integrated
sampling makes this framework suffer from an exhausted training time, since we need to send the predicted output back to decoder (Run) and repeat this Run for T times. During testing, we can use the output corresponding to its own Run (Step) as the predicted result, as well as the output of the last run (Final). As shown in the first 2 columns in Table 10, the performance of outputs in the last run (Final) is better than that of Step mode, which means the leftward information flow still exists to break the auto-regressive property in data forecasting even though the mask is adopted on the input data. For fair comparison, the models for testing are trained in one GPU and the training time are all less than 50 hours.
Typically, missing value still exists in the original dataset. During experiment, we use the global mean to replace the missing value (Mean). We also compare the prediction performance between the input with Mean Filling and Complemented Input of Sec. 3.3. and the results in Table 10 shows the Complemented Input does not lead to performance improvement but increase the training workload. Consequently, we make encoder and decoder share the same input to reduce the memory usage and training time while our model achieves better performance for long-term prediction.

heading: D DATASET DESCRIPTION
text: D.1 DATA AGGREGATION ON KDD 2018 DATASET
The original KDD 2018 dataset consists of an Air Quality data of 35 locations and an Meteorology dataset of 18 locations. each dataset contains 6 different measurements. During experiment, Luo et al. (2018a) select 11 common locations between the two datasets and the measurements of paired locations are concatenated. The location pairs are described in Table 11. Since the unit of some measurements are label-based, e.g., the measurement weather denotes the types of weather including sunny, rainy and etc, these label are replaced with value such as 1, 2, ..., 9. As the range of different measurements varies, Luo et al. (2018a) first apply Z-score normalization for each measurement and the MSE is calculated based on normalized data while the metrics calculation for other dataset are based on the original data.
D.2 NYC TRAFFIC DATASET
Traffic volume extraction: We extract the traffic volume from images using a Faster R-CNN with VGG16 backbone, trained on the MIO-TCD Luo et al. (2018b) dataset. The dataset contains 110k training images, with bounding box annotation for 11 vehicle categories (articulated truck, bicycle, bus, car, motorcycle, motorized vehicle, non-motorized vehicle, pedestrian, pick-up truck, single-unit truck, work van). On the NYC-traffic dataset, we manually annotated bounding boxes for a portion of images to evaluate the vehicle detector. Our model achieves 73% precision and 54% recall with an IoU of 0.5. To construct the NYC traffic time series, we use the model to extract and then sum up the #cars of 11 different types (articulated truck, bicycle, bus, car, motorcycle, motorized vehicle, non-motorized vehicle, pedestrian, pick-up truck, single-unit truck, work van), in non-overlapping 5-minute intervals.
Burst loss simulation: As shown in Fig. 11, we term the loss area marked as blue as the Burst loss area where for a certain camera, the data is continuously missed for a δ time slots. After statistics and analysis, we found the the length of time slots 2 ≤ δ ≤ 134. Then, for those time slots of burst loss, we calculate the mean µ = 6.350773 and standard deviation δ = 9.809643. With the mean and standard deviation, We model the generation of burst loss as Gaussian process.

heading: E DISCUSSION
text: E.1 VISUALIZATION OF DIMENSION-SPECIFIC ATTENTION MAP
We provide the attention map examples extracted from the last CDSA layer. As shown in Fig. 4(b), correlation exists between different measurement, i.e., PM2.5 and PM10 are highly correlated. As shown in Fig. 12(a), the estimation of PM2.5 and PM10 is also highly relied on each other, i.e., for the estimation of PM2.5, the color in second unit, representing the weight of PM10, is darker than the rest in the first row.
As shown in Fig. 12(b) and Fig. 13, the arrow/unit with deeper color indicates a higher weight and the index of location can be found in Table 11. According to the map in Fig. 13, in most cases, neighbouring locations often share higher attention weight, e.g., the estimation at location 1 is mainly relied on the available data from location 2, location 3 and location 4. However, the estimation of location 11 is not relied on its neighbor (location 6), instead, it is mainly relied on the location 8. We think this relation is induced since both location 11 and location 8 are the center of express way while they are away from the urban area. Thus, the air condition from those two location my highly correlated. E.2 VISUALIZATION OF ATTENTION MAP FOR CROSS-DIMENSIONAL IMPUTATIONTimeslots (hour) 0 30 50
500
300
100
10 20 40
PM 10 PM2.5
M ea
su re m en
ts
1t 2t
PM2.5 PM10
PM2.5
PM10
41 5527 41 5527
27
41 55
27
41 55
?
Besides the sample in Fig. 4 where the missing value can be estimated from the cross-dimensional available data, Fig. 14 visualizes another example and further shows that when predicting missing value A, our model pays strong attention to available values C and D while also some attention to another missing value B.

heading: E.3 RUNNING TIME COMPARISON
text: Following the model hyper-parameter setting in Tabls. 1, we further compare the average running time for one segment during testing. As the way of attention modeling determines the computational efficiency, computation method with higher FLOPs also leads to longer running time. As shown in Table. 12, the running time of Joint is much higher than the rest 3 methods. Since the computation schemes of Shared is similar with the Decomposed, while the number of trainable variables of Shared is much less that of Decomposed, the average processing time of Shared is a bit smaller than the running time of the Decomposed.
E.4 DETAIL OF FORECASTING BY USING THE IMPUTED TIME SERIES
1.79676 1.49652 1.66936
1.597 1.602 1.451 1.39012
As described in the main paper, we use the 23-day data of NYC-Traffic for further forecasting. We split the data into two segments, one segment contained the data of the first 20 days and the other contained the data of the rest 3 days. We used the imputed data from the first segment to forecast the value of the second segment. To provide comprehensive comparison, according to different missing ratio (i.e., 30%, 50% 70%), we remove the value of some units in the first segment according to burst loss and then feed the segment with missing value into the data statistical imputation model (i.e., Mean filling, Kriging expo, Kriging linear) and deep learning methods (MTSI, DCRNN, BRITS, CDSA(ours)). Then, we feed the imputed data in to prediction model (ARIMA, Random Forest Friedman et al. (2001)) and evaluate the forecasting performance in terms of RMSE. According to the Fig.15, we can find our proposed model always outperforms than other method.

heading: F EXPERIMENT SETTING

heading: F.1 BASELINES
text: We compare our method with both deep learning based methods and statistical methods while all statistical methods are adopted for each measurement individually.
• Mean Filling: Replace the missing data with global mean. • Auto Regressive Akaike (1969); Orfanidis (1988); Haddadi et al. (1998): Aggregate both
forward and backward auto regressive on each time series data by weighted average and replace the missed value, implemented in MATLAB. • Kriging: Fit linear (Exponential) function between data variance and geodesic distance and
replace the missing value w.r.t. the geo-location. Applied on each time slot and implemented in PyKrige2. • Multi-Imputation by Chained Equation (MICE) Buuren & Groothuis-Oudshoorn (2010):
Replace missing data by creating multiple imputations with chained equations. • k-nearest neighbor (KNN) Speed (2003); Hastie et al. (1999); Troyanskaya et al. (2001):
Use normalized Euclidean distance to find similar samples and impute the missing values with weighted average of its neighbors. • Matrix Factorization (MF) Yu et al. (2016): Factorizes the incomplete matrix into two
low-rank matrices and fill the missing values by l1 sparsity and l2 penalty. • MTSI Luo et al. (2018a): Treat each measurement at different location as variables identi-
cally, align the time series of all the variables into one matrix and impute the missed data through RNN-based GAN. • ST-MVL Yi et al. (2016): Replace the missed value by using geo-location information • BRITS Cao et al. (2018): Treat the missing data as variable of the bidirectional RNN and
impute by getting the delayed gradients for missing values in both forward and backward directions. • IIN Zhou & Huang (2018): Connect three bi-directional RNN in turn while the output of
each RNN are adopted in loss calculation. • FC-LSTM Sutskever et al. (2014): Apply fully connected LSTM hidden units in Encoder-
Decoder RNN structure to model the temporal dependency and forecast the traffic speed. • DCRNN Li et al. (2018): Model traffic speed as signals diffused over bidirectional graph
and modify the gated recurrent unit in encoder-decoder RNN structure to capture the spatiotemporal dependency. • DST-GCNN Wang et al. (2018): Build a two-stream framework for model the condition
stream and flow evolution while the model use the predicted condition to forecast the flow. • GaAN Zhang et al. (2018b): Use attention mechanism for each location and then graph
aggregator to assemble the neighbor nodes and impute the missing value.
2https://github.com/bsmurphy/PyKrige
F.2 IMPLEMENTATION DETAILS
The dataset in METR-LA also has missing data while the missing rate is 91%. Thus, the segment sample whose all units are zero, i.e., all-zero sample, exists. During training, the all-zero sample (in training set) essentially has no contribution for the model training. During testing and validation, the evaluation metric will of such samples will not be counted.
Data Pre-processing We apply Z-score normalization on each measurement as (9) respectively and fill the missing value with 0.
X(t, l)′ = X(t, l)− µ
σ ,
µ =
∑T t=1 ∑L l=1 X(t, l)
TL ,
φ =
√∑T t=1 ∑L l=1(X(t, l)− µ)2 TL .
(9)
Optimizer We use the Adam optimizer Kingma & Ba (2014) while the initial learning rate in each epoch is set as
lr(e) = r0 × αceil(max(0,e−d)/i). (10)
F.3 KDD-2015
For KDD-2018, Luo et al. (2018a) adopts content loss in a GAN-based model to train the random noise and then estimate the missing value, i.e., for one data segment, according to the specified missing rate, some available data will be held to evaluate the imputation performance, while the remaining available data will be set as groundtruth to calculate the content loss. Our experiment on NYC-Traffic follows the same experiment setup as KDD-2018 while the noise is replaced with the remaining available data and the model parameter is trained according to the content loss. Thus, there is no division of training, validation or testing since the training loss is not calculated from the held available data.
To comprehensively develop our experiment, we also adopting our method on KDD-2015 and follow the experiment setup in Yi et al. (2016); Cao et al. (2018) while the available data will be trained to predict the held data directly.
KDD-2015 (Zheng et al. (2015)). This dataset focuses on air quality and meteorology. It contains data recorded hourly, ending up with totally 8,759 time stamps. PM2.5 measurement is recorded at 36 locations and Temperature and Humidity are recorded at 16 locations in Beijing from 05/01/2014 to 04/30/2015, with natural missing rate 13.3%, 21.1% and 28.9% respectively. We treat those two subsets as two separate tasks and evaluate our method on each task separately. Following Yi et al. (2016), data in 3rd, 6th, 9th and 12th months are for testing and the remaining months are for training. We randomly select 36 consecutive time slots to train our model and evaluate Mean Absolute Error (MAE) as well as Mean Relative Error (MRE).
In order to simulate experiments for imputation, besides the natural missing data, for PM2.5 we follow the strategy used in (Yi et al. (2016); Cao et al. (2018); Zhou & Huang (2018)) to further manually remove entries and hold the corresponding value as ground truth. The imputation task is to predict values of these manually removed entries. For Temperature and Humidity, we follow Zhou & Huang (2018) to randomly hold 20% of available data.
KDD-2015. Table 14 shows that for PM2.5, our method outperforms the traditional methods significantly and achieves comparable MAE as IIN (Zhou & Huang (2018)) while better MRE than IIN (Zhou & Huang (2018)). For Temperature and Humidity, our method consistently outperforms state-of-the-art methods.
Since the Decomposed draws attention maps as the Independent but shares Value as the Joint, it reduces the computational complexity significantly. As shown in Table 15, we also evaluate these methods on KDD-2015 datasets and the Decomposed achieves the best performance.

heading: F.4 METRICS
text: Suppose x = [x1, x2, ...xN ] ∈ RN represents the ground truth and x̂ = [x̂1, x̂2, ...x̂N ] ∈ RN represents the predicted value.
Root Mean Square Error (RMSE)
RMSE(x, x̂) = √√√√ 1 N N∑ n=1 (xn − x̂n)2
Mean Squared Error (MSE)
MSE(x, x̂) = 1
N N∑ n=1 (xn − x̂n)2
Mean Absolute Percentage Error (MAPE)
MAPE(x, x̂) = 1
N N∑ n=1 ∣∣∣∣xn − x̂nxn ∣∣∣∣
Mean Absolute Error (MAE)
MAE(x, x̂) = 1
N N∑ n=1 |xn − x̂n|
Mean Relative Error (MRE)
MRE(x, x̂) = ∑N n=1 |xn − x̂n|∑N
n=1 xn

heading: F.5 MODEL HYPER-PARAMETER
text: Since there are missing data (Naturally missing data) in the original dataset, to evaluate the model performance, we manually remove some of the available observation (Manually removed data) and hold those entries’ value as ground truth for evaluation. The rest data are termed as Available data. Thus, as a counterpart of Naturally missing data in the original dataset, the Naturally available data consists of Manually removed data and Available data.
KDD 2015: Following the setting in Yi et al. (2016), we split the data into training set and testing set. During training, we send the Available data into the model to estimate the missing data while the loss is calculated based on the Manually removed data. During testing, we send the Available data into the model and the metric is calculated based on the Manually removed data. For model structure on Air Quality, we set (dL, dT , v) as (6, 12, 3) and 12 heads in each layer. We set (r0, α, d, batch size) as (0.003, 0.2, 60, 23) for at most 100 epoch. The model hyper-parameter on Meteorology is same with that for Air quality except for dM = 15.
KDD-2018: Following the setting in Luo et al. (2018a), we don’t split the data for testing or validation while the imputation task assumes the completed data has no label. Thus, we send the Available data into the model while the loss is calculated based on Available data while the evaluation metric is calculated based on the Manually removed data. For the model structure, we set (dT , dL, dM , v) as (30, 6, 14, 3) and 12 heads in each layer. We set (r0, α, d, batch size) as (0.003, 0.2, 60, 20) for at most 100 epoch.
NYC: Like KDD-2018, the imputation task assumes the completed data has no label. For the estimated data, the loss is calculated on the Available data part while the evaluation metric is calculated on the Manually removed data part.
METR-LA: Following the setting in Li et al. (2018), we split the data into training set and testing set where this prediction task assumes the predicted data has labels. For model structure, we set (dT , dL, v) as (14, 6, 3) and there are 16 heads in each layer. During training, we set (r0, α, d, batch size) as (0.008, 0.5, 40, 16).

Review summaries:
This paper proposes a Transformer-based model with cross-dimensional self-attention for multivariate time series imputation and forecasting. The authors consider time series with observations collected at different locations (L), for different measurements (M), and at different timestamps (T). The authors describe 4 different self-attention mechanisms based on how the three dimensions, and it turns out the proposed Decomposed approach achieves the best performance and has moderate model complexity among the four. Experiments on several traffic and air quality datasets show the superiority of the proposed model.