Multi-version Tensor Completion for Time-delayed Spatio-temporal Data

1 Introduction

Data completion is the process of estimating missing elements of an N-dimensional array using the observed elements and some presumed structural properties of the data, such as non-negativity and/or low rank. In many realworld spatio-temporal data, missing values are often due to delays or failures in the data acquisition or reporting process. Low-rank tensor factorization and completion is a popular approach to data denoising and completion, with applications ranging from image restoration [Liu et al., 2012; Xu and Yin, 2013] to healthcare data completion [Acar et al., 2011;

∗Equal contribution.

Wang et al., 2015], recommendation systems [Almutairi et al., 2017; Song et al., 2017], and link prediction [Lacroix et al., 2018]. Low-rank tensor factorization models have also been proposed for analyzing and making predictions based on spatio-temporal data [Bahadori et al., 2014; Araujo et al., 2017; Kargas et al., 2021].

Existing methods usually assume that missing elements are randomly distributed, and the noise for each tensor element is zero-mean i.i.d. (independent and identically distributed). However, in many real-world tensor data, such as medical claims data and public health reporting data, both assumptions are violated. For example, insurance claim processing centers receive daily data from third-party providers, where there is often a time lag between data generation and data loading. There are time delays when hospitals report COVID case counts to the government. Such data can be modeled as multiple input tensors with different time-dependent noise levels, where more recent data are more unreliable due to the data acquisition delays. Periodically, updates on tensor elements arrive to produce a new version of the tensor. Over time, we observe multiple versions of the underlying tensor. To model such multi-version tensors, we have to deal with the following challenges.

• Atypical noise model: In our context, the main source of noise is under-reporting. Thus noise is non-positive and negatively correlated in the update mode (if one update is unusually large, the next is likely to be smaller). Underreporting also exhibits predictable patterns (e.g., holidays and weekends).
• Devising efficient schemes for incremental updates: As data updates arrive incrementally over time for both new and historical tensor slices, it is challenging to update the tensor factorization efficiently. Existing works [Sun et al., 2006; Zhou et al., 2016] in dynamic tensor factorization assume that the historical data slices remain unchanged, and only new slices are added.

To address these challenges, we propose a Multi-version Tensor Completion (MTC) framework. Our key contributions are as follows:

• Multi-version tensor model: Instead of modeling the different versions of the data tensor per se, we introduce an

extra tensor mode capturing the data updates across multiple versions of the data tensor. This allows us to explicitly model the update dynamics in latent space.

• Online factorization: We propose a low-rank tensor model to accurately track the updates over time. To estimate the tensor values at any given time point, we perform marginalization over the extra tensor mode. We show how we can extend our approach to track the tensor data over time, leading to significant computational savings.
• We compare MTC against several tensor and time series baselines on multiple real-world spatio-temporal datasets. We observe that MTC can accurately predict the groundtruth values of many real-world tensors and achieve up to 27.2% lower root mean-squared-error than the best baseline method.
• Finally, we extend our method to track the tensor data over time, leading to significant computational savings.

2 RelatedWork

CANDECOMP/PARAFAC decomposition (CPD) [Harshman, 1970] and Tucker decomposition [Tucker, 1966] have been used for imputing missing data for image restoration [Liu et al., 2012; Xu and Yin, 2013], grade prediction [Almutairi et al., 2017], link prediction [Lacroix et al., 2018], and healthcare data completion [Acar et al., 2011; Wang et al., 2015]. Under certain low-rank assumptions, it is possible to recover the missing entries [Yuan and Zhang, 2016]. Tensor models have also been used for joint tracking and imputation [Sun et al., 2006; Sun et al., 2008; Nion and Sidiropoulos, 2009; Mardani et al., 2015; Song et al., 2017] which enables dynamic tensor factorization with streaming data via multi-linear decomposition. Recently, several nonlinear tensor factorization models have been proposed. NeuralCP [Liu et al., 2018] employs two neural networks where the first one predicts the entry value and the second one estimates the noise variance. COSTCO [Liu et al., 2019] is a convolutional tensor framework that models the nonlinear relationship using the local embedding features extracted by two convolutional layers. SPIDER [Fang et al., 2020] is a probabilistic deep tensor factorization method suitable for streaming data.

These methods assume that past slices remain unchanged and the noise in the tensor elements is i.i.d. zero-mean. These assumptions are not appropriate for modeling noise due to under-reporting, as explained earlier, which motivates our Multi-version Tensor Completion (MTC) framework.

3 Problem Statement

Given a spatio-temporal tensor of I locations and J features over time, we first introduce two time concepts:

Generation date (GD) is the time when data items are generated.

Loading date (LD) is when we receive the data items. On an LD t, we can receive data items related to multiple locations and signals that have been generated in different GDs.

Next, we will introduce several tensors at loading date t: the observed tensor Zt , the ground-truth tensor ˜Zt , and the update tensor Xt all up to time t.

On each LD t, we receive data updates {xijs(t)}s∈St for location i, feature j and GDs St which are earlier than t i.e., s ≤ t ∀s ∈ St. The total number of data that are generated on GD s and are available on LD t is computed by summing over all data received in the interval [s, t], i.e., zijs(t) = ∑tt′=s xijs(t′ ). Stacking all {zijs(t)} into a 3-way tensor, yields Zt ∈ RI×J×St + , where Zt (i, j, s) = zijs(t). Tensor Zt contains the aggregated number of updates which have been received at time t. Because of delays in the data acquisition process, data corresponding to the most recent GDs are usually under-reported. Let ˜Zt be the ground truth tensor at time t. The signal model can then be described by

where Nt ∈ RI×J×St denotes a noise tensor. Our goal is to estimate the groundtruth tensor ˜Zt . We identify three main challenges in estimating ˜Zt .

1. The values corresponding to the latest GDs in Zt , i.e., the latest frontal slabs of Zt are under-reported and thus very noisy. We need to unveil the update patterns from past data and correct these highly corrupted slabs.
2. The noise distribution is unknown in practice. The number of noise changes over time as more corrections to the tensor is introduced. Also, there are patterns in the noise. The method needs to be flexible enough so that it can capture different noise distributions.
3. The dimension corresponding to the GDs in Zt is gradually growing as t increases and more data are introduced. Note that new updates arrive continuously. Therefore, the method needs to be adaptive such that it cannot only correct the under-reported numbers but also refine itself efficiently as soon as new data are received.

To tackle these challenges, the key idea is to track the update tensor. We may assume that data corresponding to a given GD is updated at most K times. This implies that at time t we receive data for GDs (t -K + 1) to t. An example is shown in Figure 1, where we assume K = 3 and Xt(k) denotes the kth update of the data for I location and J features generated on GD t. Each column in Figure 1 shows the updates for the data on GD t. In this example, GD (t - 3) has received three updates at time t-3, t-2 and t-1 respectively. By summing over the three updates, we have the ground truth values. Therefore, we have transformed the problem into an equivalent 4-way tensor completion problem. The final estimate of the target tensor ˜Zt can be estimated by first imputing the missing values of Xt , and then marginalizing over the third mode

where Xt (: , : , k, : ) is the MATLAB notation that selects data from a tensor.

4 Multi-version Tensor Completion (MTC)

Formulation In this work, we approximate X using a lowrank CPD model [Harshman, 1970]. For the ease of notation, w drop the underscript t. CPD expresses tensor X as a sum of F rank-1 tensors,

where A ∈ RI×F , B ∈ RJ×F ,C ∈ RK×F and D ∈ RS×F are the factor matrices for location, feature, loading date (LD) and and generation date (GD), respectively. This kind of lowrank model is intuitive in many applications: Counties in the same state can share the same data loading schedule, and thus they have similar loading patterns. This is also true for updates data collected and sent through the same insurance company. Therefore, there is a high correlation among the data points for similar locations and features, indicating lowrank structures in the tensor. Suppose the optimal low-rank approximation of X is X ≈ [[A, B, C, D]].

Recall that the received data can be split into two parts - a fully observed tensor corresponding to GDs from 1 to (S - K +1) and an incomplete tensor corresponding to GDs from (S -K + 2) to S. The loss function is then expressed as

where

where ‖ · ‖F is the Frobenius norm, ds denotes the s-th row of D, Ωs denotes the index set of the observed entries of X(:, :, :, s), and PΩs (·) is the operation for tensor completion by only keeping the entries from index set Ωs and zero-out the remaining entries, and θ = {A, B, C, D}.

In the context of spatio-temporal analysis, we would like to keep the latent factors interpretable. Towards this end, we employ a graph regularizer which regularizes the latent subspace of locations according to its graph Laplacian L, such that two locations are sufficiently close to each other in the latent subspace if they are connected in the graph and have similar loading patterns. This regularization term is

where T is the transpose and ρA is the regularization parameter. Furthermore, for the LD factor C and GD factor D,we impose smoothness regularization, i.e.,

We choose Γ to be a Toeplitz matrix holding (-1, 2, -1) as its principal three diagonals and zeros elsewhere. Γ promotes smoothness on the latent factors so that data corresponding to adjacent time points has a smooth transition in LD and GD dimensions. Let R(θ) = R(A)+R(C)+R(D). We propose solving

where the constraint θ ≥ 0 is selected to ensure that the imputed values are non-negative.

Initialization Due to the special structure of the tensor completion problem in Equation (4) we can obtain a good initialization for our algorithm by first solving

for s = 1, . . . , S -K + 1. This is a standard tensor factorization problem which can be solved efficiently [Xu and Yin, 2013; Fu et al., 2020]. By solving the above problem we obtain initial estimates of A,B,C. Then, F2(θ) reduces to a non-negative least squares problem, i.e.,

for s = S - K + 2, . . . , S, a convex optimization problem which can be solved optimally, where vec(·) is the vectorization operator.

Optimization Having obtained initial estimates for A,B,C,D we turn our attention to optimization Problem (5). A common way of tackling Problem (5) is through alternating optimization. We introduce an auxiliary variable Y ∈ RI×J×K×S of the same size as X and alternatively estimate the factor matrices while imputing the missing data leading to an expectation maximization (EM)-like approach. Now let us express optimization problem (5) as

where

Optimization problems (5) and (6) are equivalent [Xu and Yin, 2013]. Problem (6) can be readily handled through alternating optimization. At each iteration we fix all variables except for one. At τ-th iteration, the subproblemwith respect to A is given by

To proceed, we define ¯Y(τ) =Y(τ) ×4W, ¯D(τ) =WD(τ) , with W = diag(√α, · · · , √α, √1-α, · · · , √1-α). Using the mode-1 unfolding matrix unfolding [Sidiropoulos et al., 2017] we equivalently express F(A) as

where H(τ) A = B(τ) ⊙ C(τ) ⊙ ¯D(τ) and matrix ¯Y(τ) (1) ∈ RI×JKS + is the mode-1 matrix unfolding which stores the mode-1 fibers of the tensor as its columns.

Let f(A) = F(A) + R(A) which can be approximated by the following quadratic function locally around ˆA(τ)

where 〈·〉 is the inner product, γ (τ) A is a parameter associated with the step-size. The gradient ∇A(f( ˆA(τ) )) is given by ∇A(f( ˆA(τ) )) = ( ˆA(τ) (H(τ) A )T - Y(τ) (1) )H(τ ) A + λLˆA(τ) and

denotes an extrapolated point, ν (τ-1) ≥ 0 is the extrapolation weight which can be updated as [Beck and Teboulle, 2009] ν (τ) = 1-e(τ) e(τ+1) where e(0) = 0 and e(τ) = 1+√4(e(τ-1) )2+1 2 . Now the subproblem w.r.t. A reduces to

which admits a closed-form expression as

Following the above analysis, we can easily write down the updates for B,C, and D, i.e.,

where ˆB(τ) , ˆC(τ) and ˆD(τ) have the same definition as ˆA(τ) in (10) 1 . Finally, the update for Y is given by

1 See Supplementary 1 for the analytic expressions of the gradients in (13)-(15).

where ˆX(τ+1) = [[A(τ+1) ,B(τ+1) ,C(τ+1) ,D(τ+1) ]], Ω denotes the index set of the observed entries of X and Ωc is the complementary set of Ω.

Summary of MTC Note that for the factor matrices, we perform local prox-linear approximation to update them. The benefits of doing this are that each subproblem admits a closed-form solution, and the algorithm is guaranteed to converge [Xu and Yin, 2013].

As an example, consider the update for A. Since the original subproblem is a least squares problem, the prox-linear mapping is known to have a Lipschitz continuous gradient. The smallest Lipschitz constant is the maximum eigenvalue of (H(τ)T A H(τ) A + λL). This is also true for the update of B,C,D, where their respective smallest Lipschitz constants are the largest eigenvalues of H(τ)T B H(τ) B , (H(τ)T C H(τ) C + ρΓT Γ) and (H(τ )T D H(τ) D +ρΓT Γ), respectively. In the above, H(τ ) B =A(τ+1) ⊙C(τ) ⊙ ¯D(τ) ,H(τ) C =A(τ+1) ⊙B(τ+1) ⊙ ¯D(τ) andH(τ) D =A(τ+1) ⊙B(τ+1) ⊙C(τ+1) . Furthermore, the subproblem w.r.t. Yτ is also convex with closed-form solution. It follows from [Xu and Yin, 2013](Theorem 2.8):

Proposition 1. MTC converges to a stationary point of (6).

The detailed steps of MTC are summarized in Algorithm 1 in the supplementary material. We note that MTC belongs to a family of block coordinate descent (BCD) algorithms. The computational bottleneck is the matricized tensor times Khatri-Rao product (MTTKRP) in the gradient calculation, which dominates the per-iteration complexity. This can be alleviated by exploiting sparsity [Smith et al., 2015]. The overall complexity per iteration of MTC is O(IJKSF).

5 Adaptive Update

The MTC algorithm leverages all data to complete the missing values, but this is at the expense of high complexity. It requires solving the tensor completion problem exactly at any instant time t. To avoid restarting MTC, we propose an efficient forward-backward propagation based strategy to estimate the missing values and update all factors simultaneously. The resulting algorithm is termed MTC-online.

Let Xt+1 ∈ RI×J×K×St+1 + be obtained after appending newly observed data ˘Xt+1 ∈ RI×J×K + in the fourth dimension, where St+1 = St + 1. In this case, some of the missing values in Xt have been completed, but new missing values are included. Therefore, the index set Ωt changes to Ωt+1 and the factor matrices should be refined accordingly. Assume that the best low-rank approximation of the missing values is

where the dimension of At+1, Bt+1 and Ct+1 are the same as their counterparts obtained at time t while the number of rows in Dt+1 is now expanded to (St + 1).

To develop a linear complexity algorithm, it is necessary to impose a proper approximation for the new data. The choice of such an approximation has a major effect on the performance of the resulting algorithm. One such option is to assume slow variations between data from t to (t + 1). This means that the (i, j, k)-th entry in ˘Xt+1 can be well approximated using the factors from the previous time step. Thus, we have

We then employ a two-step approach to refine the factor matrices, including one forward-propagation (FP) update and one backward-propagation (BP) update, where in each propagation, only one iteration is employed to refine the factor matrices. Specifically, in the FP step, we update the row of the latent factor corresponding to the new GD. The associated optimization problem is described as

which can be solved optimally using projected gradient descent with light-weight operations. After we obtain d, we perform the BP. We first append d to Dt and then complete all missing values via

where Ωt+1 is the index set of the the missing values at time (t+1) and Ωct+1 is its complement. We refine all factor matrix through one iteration of the gradient updates in (12)-(15) in a coordinate manner.

6 Experiments

We evaluate the MTC's performance in terms of prediction accuracy and scalability. We use the following metrics: Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and R2 score. To evaluate scalability, we record the CPU running time. All methods were trained on 2.6 GHz 6-Core Intel Core i7, 16 GB Memory, and 256 GB SSD.

6.1 Experimental Setup

Data We consider the following datasets for evaluation.

• Chicago-Crime dataset [Smith et al., 2017] includes 5,884,473 crime reports in Chicago, ranging from 2001 to 2017. We us biweekly aggregation of the reports resulting in a 77×32×442×10 tensor of density is 0.219. These are 77 locations, 32 crime categories, 442 GDs, and 10 LDs.
• COVID-19 dataset [Dong et al., 2020] summarizes the COVID-19-19 daily reports from Johns Hopkins University. It records 3 categories, i.e., new cases, deaths, and hospitalization of 51 states in the United States between April 15, 2020, to October 31, 2020. The tensor's size is 51×3×200×8, and its density is 0.616, where there are 51 states, 3 features, 200 GDs, and 8 LDs.
• Patient-Claims dataset extract patient claims data from a proprietary claims database from 2018. Each data sample consists of a zip code, diagnosis (i.e., disease), date of diagnosis (i.e., GD), and date that claims loaded to the database (i.e., LD). We transform the zip codes to US counties and map the diagnosis codes to 22 diseases according to the ICD-10-CM category 2 . Finally, we count the claims for a given county, disease, GD, and LD and create a 4-way tensor. The tensor's size is 3027 × 22 × 52 × 12 (over 41m entries) and its density is 0.534, where there are 3027 counties, 22 diseases, 52 GDs and 12 LDs.

2 https://www.icd10data.com/ICD10CM/Codes

Evaluation Strategy We consider both static and dynamic cases. We select the first S GDs from each dataset for the static case while setting the remaining for the dynamic case. Here, S is set to 421, 168, and 46 for ChicagoCrime, COVID-19, and Patient-Claims, respectively. Therefore, their corresponding numbers of GDs for the dynamic case are 21, 32, and 6, respectively. In the static case, the data corresponding to the first (S-K+1) GDs are fully observed, and the last (K - 1) GDs have missing LDs, so the corresponding true values {˜ zijs} in ˜Z are under-reported. Our target is to predict {˜ zijs} corresponding to the last (K-1) GDs.

Baselines We compare our method against low-rank tensor factorization methods, a neural network tensor completion method, and deep sequential models.

• Naive: We use the received data until t to benchmark the gap between the received data Zt and the actual data ˜Zt .
• SDF (3-way): This is a 3-way tensor completion approach using the Structured Data Fusion (SDF) framework in Tensorlab [Vervliet et al., 2016]. We complete the missing data of Zt and use the learned model as our estimate for ˜Zt .
• SDF (4-way): We consider a variation of MTC for comparison. This is a 4-way tensor completion approach using SDF/Tensorlab [Vervliet et al., 2016] without smoothness and Laplacian regularization.
• COSTCO [Liu et al., 2019]: This is a convolutional neural network-based model for tensor completion.
• Auto Regressive Integrated Moving Average (ARIMA): We train distinct ARIMA models for different pairs of location and attributes/signals.
• Long Short Term Memory (LSTM) network: We train LSTM to predict the values for the subsequent GDs.

The implementation details are provided in the appendix.

6.2 Results

Table 1 summarizes the results of the static case. The Naive method performs poorly, which indicates that the timedelayed data differ significantly from the ultimate ground truth. The poor performance of SDF (3-way) also verifies such an observation. The 3-way completion approach can complete the missing values of the latest slabs of Zt but cannot correct the under-reported cases. MTC has consistently lower MAE and RMSE than the other methods in the three datasets. Note that in the Patient-Claims dataset, which is the most challenging one, MTC performs the best, and its performance is followed by the SDF (4-way) algorithm. Since both MTC and SDF (4-way) are 4-way tensor completion approaches, the performance improvement of MTC is achieved by the well-designed prox-linear gradient approach and the regularization terms. We empirically observed that COSTCO was more susceptible to overfitting and did not perform well.

Accuracy In Table 2, we calculate the RMSE, MAE, and R2 for each method in the dynamic case. The values were averaged over all GDs in each dataset. It can be seen that the proposed methods achieve the best RMSE and R2 in the three datasets. We observed that COSTCO could easily overfit the data, so it did not perform well. ARIMA, LSTM, and SDF (3way) did not use the LD information, thus failing to achieve satisfactory results. This indicated that using an additional mode - LD can help solve the data delayed issue in spatiotemporal tensors. In the Patient-Claims dataset, MTC-online has a smaller MAE than MTC. In all other cases, MTC has slightly better performance than MTC-online. However, this is at the expense of high complexity, since MTC needs to be restarted once new data received. This is true for the other baselines as well. On the other hand, MTC-online updates all factors in an inexact way with a single proximal gradient update, resulting in very low-complexity.

Scalability Figure 2 shows the CPU running time of MTConline, MTC, ARIMA, and SDF (4-way) on the PatientClaims dataset. We did not include LSTM for comparison because it had a much higher time complexity for training. Figure 2 shows COSTCO and ARIMA has the highest complexity - the CPU time for both methods was close to one hour. ARIMA is slow we need to train a separate model for each location. SDF (4-way) and MTC have similar running time. MTC-online is an order of magnitude faster than the best baseline due to the approximate gradient updates.

Model interpretation Figure 3 shows the latent components corresponding of the GD mode in the COVID-19 dataset. Figure 3 shows that Cluster 2 consisted of the states hit by Covid-19 in the first wave in the United States. The Cluster-1 had about two months delay, exactly for the second wave, including FL, CA, TX, AZ, etc. The five clusters have similar exponential decay trend. Interestingly, Cluster 4 had the lowest loading speed, meaning that data corresponding to this cluster took a longer time to be received than those in the other components. Four out of the top five counties in Cluster 4 were from California, where the laws on processing/collecting claims data are more strict than the other states. More results are reported in the supplement.

7 Conclusion

This paper studies the time-delayed spatio-temporal tensor data, which corresponds to many real-world data aggregation applications. We formulated the problem as a multi-version tensor completion (MTC) problem by introducing an extra mode to capture the data updates. The final estimate of the target tensor is estimated by first imputing the missing values and then marginalizing over the extra mode. We proposed static and online version of MTC algorithms to tackle this problem. The proposed methods employ a graph regularization to reserve the location correlations in latent subspace and smoothness. The experimental results on several real datasets have demonstrated the advantages of the proposed methods.