Spatio-Temporal Tensor Sketching via Adaptive Sampling

Mining massive spatio-temporal data can help a variety of real-world applications such as city capacity planning, event management, and social network analysis. The tensor representation can be used to capture the correlation between space and time and simultaneously exploit the latent structure of the spatial and temporal patterns in an unsupervised fashion. However, the increasing volume of spatio-temporal data has made it prohibitively expensive to store and analyze using tensor factorization. In this paper, we propose SkeTenSmooth, a novel tensor factorization framework that uses adaptive sampling to compress the tensor in a temporally streaming fashion and preserves the underlying global structure. SkeTenSmooth adaptively samples incoming tensor slices according to the detected data dynamics. Thus, the sketches are more representative and informative of the tensor dynamic patterns. In addition, we propose a robust tensor factorization method that can deal with the sketched tensor and recover the original patterns. Experiments on the New York City Yellow Taxi data show that SkeTenSmooth greatly reduces the memory cost and outperforms random sampling and fixed rate sampling method in terms of retaining the underlying patterns.


Introduction
The increasing availability of spatio-temporal data has brought new opportunities in application domains including urban planning, event management, and informed driving [9,26]. Unfortunately, the rapid growth in these data streams can be prohibitively expensive to store, communicate and analyze. In addition, the high-dimensional, multi-aspect spatio-temporal data poses analytic challenges due to the correlations in the measurements from both time and space. Moreover, human-intensive and domain-specific supervised models are not tractable due to the constant and evolving deluge of measurements.
Given the high-dimensional, multi-aspect nature of spatio-temporal data, a tensor serves as an efficient way to represent and model such data [6,9,26]. As an example, each element of the tensor can represent the occurrences of an event at a specific location (encoded as latitude and longitude) within a specific time interval. Compared with the low-dimensional matrix-based methods such as [16], tensors can capture the correlation between each mode. Furthermore, tensor factorization offers an unsupervised, data-driven approach to identify the global structure of the data via a high-order decomposition. It is also more interpretable compared to deep learning methods [32]. Unfortunately, existing models require the full tensor information (i.e., all the data samples must be stored) and do not readily scale to extremely large tensors.
The computation and storage limitations of large tensors motivate the need to approximate such tensors with relatively small "sketches" of the original tensors. Not only are these manageable-sized tensor sketches more readily stored on a single machine, the computationally expensive tasks such as tensor decomposition can be performed on the smaller tensors while still preserving the underlying structure. Although sketching has been proposed as a linear algebra tool for reducing the memory cost, traditional methods involve the linear transformation of the original tensor with a "fat" random projection matrix to reduce the dimensions [23,29]. However, the transformation has inherent shortcomings: 1) the design of the random projection matrix may not capture the evolving patterns without a priori information; and 2) continuously sketching is computationally expensive since it involves the inversion of a fat random projection matrix.
An alternative sketching approach is based on tensor sparsification -subsampling the original tensor while preserving the original tensor structure. Compared with the random projection method, sampling incurs negligible online complexity. Existing tensor sparsification methods include random sampling based on the tensor spectral norm [20], sampling according to the entry values [30], and sampling according to the pre-computed tensor distribution [5]. However, these tensor sparsification algorithms suffer from the following limitations: 1) they cannot deal with streaming data where measurements are not available a priori and the tensor is incrementally updated; 2) there is no formal mechanism to reconstruct the original streams [8]. While tensor reconstruction is often a byproduct of dealing with missing data, reconstruction of the original streams can help track dynamic and abrupt changes at particular locations.
In this paper, we propose SkeTenSmooth, a factorization framework that uses adaptive sampling to generate tensor sketches on-the-fly and preserves the underlying global structure. We explore the problem of serially acquired time slices (measurements appear once they are available). We introduce SkeTen, a tensor sketching method that adaptively adjusts the sampling intervals and samples time slices according to the feedback error between the prior estimate and the true value. Thus it can capture the time slices that are not well modeled by the prior forecasting model, while avoiding the storage of the time slices that contain redundant information (i.e., patterns that are already captured). Furthermore, we propose a novel method SkeSmooth to decompose the tensor sketches and reconstruct the underlying temporal trends. Unlike previous tensor factorization algorithms that deal with randomly missing entries, the sketched tensor has missing time slices (i.e., no data from a specific time point) which poses additional challenges. Thus, we introduce a robust tensor factorization algorithm that incorporates temporal smoothness constraints using auxiliary information about the data gained from SkeTen. Figure 1 shows the overall flow of the proposed framework. We briefly summarize our contributions as: 1) Tensor sketching with adaptive sampling. We propose SkeTen, a tensor sketching technique that helps process large volumes of data with low memory requirements and preserves the underlying temporal dynamics.
2) Tensor factorization with smoothness constraint. We propose a tensor factorization framework, SkeSmooth, that decomposes the smaller tensor "sketches" with smoothness constraints to achieve robust recovery of the underlying latent structure and the missing entries.
3) Case study on New York taxi data. We illustrate the ability to produce small tensor "sketches" and generate smooth temporal factors on real data.

Preliminaries and notations
This section summarizes the notations used in this paper. Note that we use X to denote a tensor, X to denote a matrix, and x to denote a vector. The mode-d matricization of a tensor is denoted by X (n) . The row and column vectors are represented by x i: , x :r respectively.

Tensor Factorization
Definition 1. (Khatri-Rao product). Khatri-Rao product is the "columnwise" Kronecker product of two matrices A ∈ R I×R and B ∈ R J×R . The result is a matrix of size (IJ × R) and defined by A B = [a 1 ⊗ b 1 · · · a R ⊗ b R ]. ⊗ denotes the Kronecker product. The Kronecker product of two vectors a ∈ R I , b ∈ R J is (CANDECOMP-PARAFAC Decomposition). The CANDECOMP-PARAFAC (CP) decomposition is to approximate the original tensor Y by the sum of R rank-one tensors where R is the rank of tensor Y. For a three-mode tensor Y ∈ R I×J×K , the CP decomposition can be represented as where a :r ∈ R I , b :r ∈ R J , c :r ∈ R K are the r-th column vectors within the three factor matrices A ∈ R I×R , B ∈ R J×R , C ∈ R K×R , • denotes the outer product.
In this paper, the spatio-temporal tensor has three dimensions. The first dimension is the temporal dimension, while the other two are longitude and latitude dimensions representing the spatial modes.

SkeTen: Tensor Sketching via Adaptive Sampling
The main idea of SkeTen is to discard slices along the temporal mode that are "predictable". Intuitively, this adaptive sampling strategy will capture sudden temporal changes of the incoming streams, which is difficult for other methods such as random sampling and fixed rate sampling to capture. Moreover, SkeTen does not require a priori knowledge of the data and can adjust the sampling interval in real-time. SkeTen consists of two main components. The adaptive sampling module uses a prediction model to measure the "predictability" of the data. The second component selects random time series projections to reduce the training time of the prediction models.

Adaptive Sampling
SkeTen introduces an adaptive sampling method that can adjust according to the detected temporal dynamics and thus can better capture the underlying patterns even with the same amount of data stored. In SkeTen, individual timeseries prediction models are developed for each spatial location. The adaptive sampling strategy then analyzes the slice feedback error, or the feedback error across all the locations. j,k represents the true value at a particular time t, where 0 ≤ t ≤ T , then the slice feedback error at time t, E t is defined as: j,k is the estimated value based on the prediction model, and δ is a small sanity bound.
The slice feedback error reflects how well the current time series models fit the current trends for each time series fiber. If there is a sudden increase across multiple spatial locations (i.e., an increase in the slice feedback error), then the sampling interval should be shortened to better capture the evolving trends.

Time Series Prediction Model SkeTen uses an Autoregressive Integrated
Moving Average (ARIMA) model to predict the values of the temporal slices. Given a time-series sequence {x i }, i = 1, · · · , t, where i is the time index, ARIMA(p, d, q) is defined as where p is the order of lags of the autoregressive model, d is the degree of differencing, and q is the order of the moving average model. α j , j = 1, · · · , p and β j , j = 1, · · · , q represent the coefficients of AR and MA, respectively, ε t is the white noise at time t, and c is the bias term. An ARIMA model is trained on each tensor fiber for the initial sampled time slices. ARIMA has several convenient properties that make it more suitable than other time series prediction models in our setting. (1) The auto-regressive part of ARIMA ({α 1 , · · · , α p }) generates the coefficients that are later exploited in SkeSmooth for the smoothness constraint in Sec. 4.1. This is essential for recovering the original temporal patterns with missing time slices. (2) It doesnt require a large number of training samples so the compression can occur much earlier than more complex models. (3) The coefficients can be used to understand the adaptive sampling for the PID controller which we will explain next.
In preliminary experiments, we tried several state-of-the-art deep learning models such as ConvLSTM [31], but the compression rates and sampled slices only improved marginally. Thus, ARIMA is an optimal choice for our setting.
Feedback Control SkeTen uses a controller system to detect rapid changes in the slice feedback error and adaptively adjusts the sampling rate. We adopt a PID controller similar to [7,14] to change the sampling interval over time. The PID controls the Proportional, Integral and Derivative errors. Proportional error is defined as Γ p = γ p E t to control the sampling intervals by keeping the controller proportional to the current slice feedback error. Integral error is where M t represents how many errors have been taken until time t. Thus the integral error considers the past errors to eliminate offset. Derivative error is defined as Γ d = γ d Et−Et−1 tm−tm−1 to prevent large errors in the future. Therefore, the full PID controller is defined as The PID errors can thus be interpreted as control actions based on the past, the present and the future [3]. We set the PID parameters according to the common practice that proportional > integral > derivative [7]. With the PID error, the sampling interval is adjusted according to where θ and ξ are predefined user-specific parameters. The adaptive sampling process is presented in Algorithm 1.

Random Projection of ARIMA Coefficients
For large scale spatio-temporal data with many spatial locations (countless values of longitude and latitude), training the prediction model (i.e., ARIMA) may be computationally infeasible. To reduce the computational cost and retain the temporal fidelity of the data, we propose a random projection algorithm to train a limited number of ARIMA models, while preserving the competitive performance of the trained models. Inspired by the Locality Sensitive Hashing (LSH) algorithm, which has demonstrated its success in Approximate Nearest Neighbor search, our proposed projection method has the following steps: (2) construct L buckets from the time series set T , each bucket i corresponds to an ARIMA model g i ; (3) map each time series ts j , j = 1, ..., M from T into bucket i according to the RMSE between the true and the predicted values using the ARIMA model g i . In this way, we treat the fitting of each time series to the randomly selected time series as a special "random projection", and then map each time series according to the least RMSE to each bucket i. This process (detailed in Algorithm 2) allows locations that may not have been spatially close to be clustered to the same model and can decrease the training time, which is the bottleneck of the adaptive sampling process.
We compare our proposed random projection algorithm for the ARIMA coefficients with the k-means algorithm where similar time series can also be clustered and thereby reducing training time. Fig. 2 (c), (d) demonstrate that despite the high computation cost of the k-means algorithm, it provided limited improvement over the random projection algorithm.

SkeSmooth: Smooth Tensor Factorization
We then consider analyzing the sketched tensor using the CP decomposition to capture the underlying multi-linear latent structure of the data and reconstruct the unsampled time slices. Built on the success of the Quadratic Variation (QV) matrix constraint [32,35], we consider the delay effect of the time series and utilize the coefficients from the pre-trained ARIMA model from SkeTen to make the temporal pattern smooth and robust to missing entire time slices.

Formulation
CP Decomposition with Missing Data. We consider the sketched tensor as an incomplete data problem where the goal is to both learn the underlying patterns of the data and reconstruct the unsampled values. Previous CP-based tensor completion algorithms [2,25] have formulated the completed tensor X as: where * denotes the element-wise product, the observed tensor is factorized as X = [[A, B, C]] and the binary weight tensor W is of the same size as X such ARIMA-based Regularization Matrix. Standard CP decomposition with missing data assumes the entries are missing at random, whereas SkeTen yields a sketched tensor with entire time slices missing. Moreover, they do not account for temporal information which can improve interpretability and robustness. In spatio-temporal datasets, successive observations at the same spatial location (i.e., x j,k ) are unlikely to change significantly. Thus, adjacent time slots are generally smooth. SkeSmooth incorporates the coefficients generated by the various ARIMA models as auxiliary information into the regularization matrix L. Based on Eq. (3), we observe that a i can be approximated as: Thus by utilizing the learned predictive models in SkeTen (the autoregressive coefficients and the order of the number of time lags, p), the algorithm can better reconstruct the missing time slices and learn the underlying patterns. Although each ARIMA model may have different order of time lags, p j,k , we use the largest order p = max p j,k and set the coefficients α p j,k +1 , · · · , α p , to zero if p j,k < p. Since each tensor fiber is clustered to the ARIMA model that fits best, we can thus get the weight for each ARIMA model as the number of tensor fibers clustered to it. We then compute the weighted average of the autoregressive coefficients to produce the regularization matrix, L, defined as: SkeSmooth Optimization Problem. The objective function for SkeSmooth after adding the smoothness constraint is: where Y denotes the sketched tensor with missing entries. Thus higher values of the regularization parameter, ρ, will yield smoother temporal factors. However, this may potentially reduce the overall fit of the tensor.

Algorithm
We use a first-order method to solve the optimization problem, similar to [2]. The gradient-based method has been shown to be robust to overspecification of

Algorithm 3: SkeSmooth
Input: Sparse tensor sketch Y, weight tensor W ∈ R I×J×K , smooth parameter ρ. 1 Randomly initialize the factor matrices A ∈ R I×R , B ∈ R J×R , C ∈ R K×R 2 while A, B, C not converge do the rank. As the regularization matrix is only enforced on the temporal mode (first mode), only the gradient for the temporal mode involves the gradient of the regularization matrix. Thus the updates for all three modes are: where Y (n) is the matricization on the n-th mode of the tensor sketch. Our gradient-based algorithm uses the limited-memory BFGS [21] for the first order optimization. The optimization process is shown in Algorithm 3.

Complexity Analysis
The time complexity of SkeSmooth is dominated by the matricized tensor times Khatri-Rao product (MTTKRP) operations, which is computed in Eq. (10) as the matricized residual tensor Y (n) −Z (n) times the Khatri-Rao product between two factor matrices. For simplicity purpose, we denote an N-th order tensor X ∈ R I×···×I as the representation of the residual tensor Y (n) − Z (n) . For each mode, computing MTTKRP takes O(nnz(X )R), where R is the approximate rank of the tensor, nnz(X ) represents the number of non-zero elements of X . While computing the term ρL T LA for the first mode takes O(I 2 R 2 ). Thus the overall time complexity for SkeSmooth is O(nnz(X )N R + I 2 R 2 ).

Experiments
We evaluate SkeTenSmooth on a large real-world dataset, the New York City (NYC) Yellow Taxi data (see Supplementary for experiment settings). The goal of our evaluation is to assess both SkeTen and SkeSmooth from three aspects: 1. Effectiveness of SkeTen: Analyze memory compression without hurting the performance. Use SkeTen to evaluate the latent structure preservation.

Effectiveness of SkeSmooth:
Evaluate with the sketched tensor on the effectiveness of the smoothness constraints. 3. Sketching result on prediction tasks: Analyze how the decomposed factor matrices perform on a downstream prediction task, to indirectly evaluate whether the sketching results are meaningful.

Dataset
We evaluate SkeTenSmooth on the NYC Yellow Taxi

Evaluation Metrics
The tensor decomposition performance is evaluated by two criteria: Factor Match Scores (FMS) and Tensor Completion Scores (TCS) [2]. FMS is defined as: x r is the r th column of factor matrices. FMS measures the similarity between two tensor decomposition solutions. It ranges from 0 to 1, and the best possible FMS is 1. We choose the factorization result of the complete tensor using CP-OPT [2] as a true factor to compare with. TCS indicates the relative error of the missing entries, and is defined as The best possible TCS is 0.

Parameter Selection
In SkeTen, parameters involved in achieving the adaptive sampling algorithm include: the PID controller parameters {γ p , γ i , γ d }, θ and ξ for the interval adjustment. The PID parameters are set as {γ p , γ i , γ d } = {0.7, 0.2, 0.1} 3 . Parameter θ determines the magnitude of the sampling interval adaptation, its impact is shown in fig. 2(a) and (b). From fig. 2(b), we observe that as θ increases, the reconstruction performance measured by FMS and the TCS decreases due to the enlarged average sampling intervals (shown in fig. 2(a)). The optimal choice shown in the figure is θ = 1. When θ = 0.1, the reconstruction performance (TCS and FMS) is worse than when θ = 1 due to an insufficient interval adjustment. ξ is determined as 3.5 for the tolerance of PID error. For SkeSmooth, the smoothness constraint ρ is chosen as 600 after grid search through {500, 600, 700, 800, 900, 1000}, and the rank of the tensor is determined as 5 after grid search in {5, 10, 15, 20, 25} (see Supplementary for more details). The number of lags is determined as p = 3 based on the majority of all the trained ARIMA models. We use the number of tensor fibers been mapped into the same bucket as weights, and compute the final coefficients for the regularization by weighted average of the coefficients of the randomly trained ARIMA models and is determined as α = [0.55, −0.19, 0.04] 4 .

Tensor Sketching
Here we consider a streaming setup where the tensor sketching will be applied to a succession of incoming time slices.
Baselines We compare SkeTen with two baseline sampling techniques: 3 We experimented with several PID settings, and it hardly affected the sampling rate as long as we set the controller ID according to proportional > integral > derivative. We thus consider our system robust to the control parameters and choose the optimal setting in empirical studies. -Fixed rate sampling samples the data with a predefined, fixed interval.
-Random sampling has been proposed as a way to sparsify the tensor [20].
By evaluating the performance of fixed rate sampling and random sampling, we can gain a better understanding of the benefits of adaptive sampling under the setting of time-evolving streams and how it may provide more informative representations of the streams at the same storage cost.
Evaluation We evaluate the three sampling methods for different levels of dropped data by using the proposed SkeSmooth algorithm. The performance is measured by both FMS and TCS. We explore the sketching level from 50% to 90%. Fig. 2 (c) and (d) show the FMS and TCS of different sampling methods when applied to SkeSmooth algorithm. As more data are sampled, FMS gradually approaches 1, which means the extracted temporal patterns are equivalent to the true patterns. Notice that when the level of sampling is 50%, adaptive sampling and fixed-rate sampling method show FMS as high as 0.9886 and 0.9803, which means that with only 50% of the original data, SkeSmooth can successfully recover the original temporal pattern. Adaptive sampling achieves slightly better reconstruction performance than fixed-rate sampling, because the daily taxi data are fairly regular and periodical so the interval adjustments for adaptive sampling helps marginally. Both adaptive sampling and fixed-rate sampling are better than the random sampling, illustrating that carefully designing the sampling strategy is critical for pattern and data reconstruction.

Tensor Decomposition with Sketches
To evaluate the performance of tensor sketching, we first need to look into how it can preserve the temporal trend and the latent structures on the temporal dimension. We decompose the sketched results to generate the factor matrices as latent patterns, and use these factor matrices to reconstruct the original tensor. We compare SkeSmooth with the following state-of-the-art methods: CP-WOPT [2], HaLRTC [17], SPC [32], BCPF [34], t-SVD [33], t-TNN [13].
Evaluation Since not all baseline algorithms are CP-based tensor completions, it is infeasible to compare the FMS. Therefore, we use TCS as the evaluation metric. Fig. 3 illustrates the TCS as a function of the level of data dropped for different sketching methods. We note that SkeSmooth outperforms the baseline tensor completion algorithms in achieving a lower TCS (relative error), verifying its robustness in dealing with the missed time slices in the sketched tensor. In particular, we observe the noticeable improvement using the ARIMA-based coefficients in comparison with the SPC algorithm which adopts the QV smoothness constraint. This observation is consistent for all sampling mechanisms.

Prediction Task
It is also critical to analyze how meaningful the generated temporal patterns are by performing a downstream prediction task using the decomposed factor  matrices. We fit a multivariate Long short-term memory (LSTM) model, which is well-suited for the traffic demand prediction task [28]. For the prediction task, the taxi demand is predicted on a daily basis where each day's number of pickups are treated as the daily demand with a train-test partition of 70-30. The LSTM model uses the factorized temporal mode factor matrix A of size I × R, where each feature within A (A i; of size 1 × R) represents a latent transportation pattern. The LSTM model is configured as one hidden layer of 30 LSTM units with ReLu activation function, followed by an output dense layer 5 . Table 1 presents the root-mean-square error (RMSE) differences for adaptive sampling (SkeTen), fixed rate sampling, and random sampling with varying percentage of data dropped from 50% to 90%. Results demonstrate that adaptive sampling outperforms other two sampling methods in maintaining the predictive power with 60%-80% data dropped. As more data are dropped, the test error (RMSE) increases, indicating that it becomes harder to fit the LSTM model.

Related Works
We review the previous work regarding the tensor sketching and tensor completion approaches in this section.

Tensor Sketching
Sketching has been proposed and investigated as an indispensable numerical linear algebra tool to help process large volume of data [29]. Thus far, there have been two research directions. Random projection is the predominant form of tensor sketching. Random projection can also be considered as tensor compression via matrix-tensor multiplication. [22] developed an approach to randomly compress a big tensor into a much smaller tensor by multiplying the tensor with a random matrix on each mode. [23] extend this work to improve the permutation matching performances of the resulting factor matrices. [11] expand the work further to the online setting. However, these compression method cannot avoid the inherent bottleneck of the matrix-tensor multiplication operators.
An alternative for tensor sketching is tensor sparsification. The goal of tensor sparsification is to generate a sparse "sketch" of the large tensor. [20] proposed random sampling based on the tensor spectral norm. Entries are sampled with probability proportional to the magnitude of the entries. [30] extended this work with: 1) the extension from cubic tensors to general tensors; 2) the criteria and treatment of "small" entry values. They proposed to keep all large entries, sample proportional to moderate entries, and uniformly sample small entries. [5] proposed a new sampling method that samples entries based on a pre-computed tensor distribution. These tensor sparsification algorithms do not consider the data streaming setting and require prior knowledge about tensor.

Tensor Factorization and Completion
Tensor completion is more considered as a byproduct when dealing with missing data in the tensor factorization process. As a result, CP decomposition and Tucker decomposition are two major types of tensor completion methods [25]. [2] proposed CP-WOPT, a tensor CP decomposition algorithm that can handle missing entries using a binary weight tensor. Tucker-based tensor completion including geomCG [15], HaLRTC/FaLRTC [17], and TNCP [18], etc. There are other approaches besides CP and Tucker decomposition, such as t-SVD [33]. [25] and [24] provide comprehensive surveys on the topic of tensor completion.
Tensor completion has also been widely applied to spatio-temporal data analysis for traffic prediction [12,27], urban mobility pattern mining [26], and urban event detection [6]. Several works have proposed the tensor completion algorithms for spatio-temporal data. [35] applied a similar tensor completion algorithm with smoothness constraints as in [32] and [35] to the imputation of missing internet traffic data, that tries to minimize the difference between adjacent time slots. [10] explored the similarity for each element within one dimension and used this as auxiliary information to better recover the original tensor. Nevertheless, the above algorithms did not consider the tensor streaming problem. [27] proposed a dynamic tensor completion (DTC) algorithm, where traffic data are represented as a dynamic tensor pattern, but still require the entire tensor.

Conclusions
SkeTenSmooth is a tensor factorization framework that generates tensor sketches in a streaming fashion using adaptive sampling mechanism and decomposes the sampled smaller tensor sketch with an ARIMA-based smoothness constraint. It is well suited for the incrementally increased spatio-temporal data analysis that can greatly reduce the storage cost and sample complexity, while preserving the latent global structure. SkeTenSmooth includes two main components: 1) SkeTen as a tensor sketching algorithm that tackles the memory efficiency issue by adaptively adjusting the sampling interval based on the model prediction error; and 2) SkeSmooth incorporates an ARIMA-based smoothness constraint, which better demonstrates the auto correlation between each time point of a time series. Experiments on the large-scale NYC taxi dataset illustrate that with adaptive sampling strategy, SkeTen can greatly reduce the memory cost by 80% and still preserve the latent temporal patterns and the predictive power. Future work will focus on the distributed setting to improve scalability and efficiency.

Supplementary
We provide experiment settings details and additional experiment results in the supplementary material.

Experiment Settings
The tensor factorization task is performed using Matlab 2019a on a r5.12xlarge instance of AWS EC2 with Tensor Toolbox Version 2.6 [4] for tensor computation and l − bf gs for optimization. For the prediction task, we build the Long short-term memory (LSTM) model using Keras with python 3.6 on a p2.xlarge instance of AWS EC2. Figure 4 illustrates the adaptive sampling results for 90% data dropped according to the PID errors. To better show the sampling results, we select a portion (200 time slots) from the long streams. From the sampled data points, we observe that for time point at 19-21, the PID error incurs a sudden increase and drop, this reflects on more frequent time slots being sampled at these points. While the error switches to a moderate fluctuation, the sampling interval becomes even. Furthermore, we explore the effect of smoothing in fig. 5. We show one of the extracted temporal pattern under the case of 50% sampled data. Note that we only present 200 time slices to give a better illustration. The upper figure shows the original pattern discovered by CP-OPT [1] with the complete tensor, compared with the pattern extracted by SkeSmooth with tensor sketches acquired by different sketching methods. Combining the lower and upper figures, we observe that with the ARIMA-based smoothness constraints, the dropped data points are smoothed out to form a complete and smooth temporal pattern that resembles the original pattern discovered with the full tensors. With only 50% of data remaining, SkeSmooth is able to recover the original temporal pattern discovered with the full tensor. With more of the time slices dropped, SkeSmooth will still be able to recover the overall trends, as the values of FMS are still above 0.5 for adaptive and fixed-rate sampling while some of the details might be lost. We notice that adaptive sampling better captures the trend for 50% data dropped due   6: Taxi demand prediction obtained with temporal patterns generated from different sampling methods with 50% data dropped, compared with the ground truth and the prediction with the temporal pattern extracted from the full tensor.

Tensor Sketching
to the error feedback system, as it samples more frequently at the time points with higher fluctuations. Fig. 7 and 8 show the FMS and TCS of different sampling methods with and without the smoothness constraints when applied to SkeSmooth algorithm. It is obvious that the smoothness constraints achieves significantly better performance as the FMS with smoothness constraints are overall higher than those without any smoothness constraints; and the TCS with smoothness constraints are overall lower than those without the smoothness constraints.

Prediction Task
Predicting the taxi demand is critical for taxi companies for real-time decision making, monitoring the traffic flow, and helping drivers to determine the most time-saving routes [19]. Fig. 6 shows the taxi demand prediction of 254 timestamps with 50% data dropped by different sampling mechanisms. We observe that with only 50% data from the original data stream, we are still able to extract meaningful temporal patterns that can successfully fit a taxi demand prediction model.

Hourly Basis Anlysis
Besides collecting data on the daily basis, we also show the results based on an hourly collected data. Fig. 9 shows the effect of smoothness constraints. We show one of the extracted temporal pattern under two extreme cases: 50% and 90% of dropped data. Note that we only present 100 time slices to give a better illustration. It can be seen that consistent with Fig. 5, the smoothness constraints play an essential part in recovering smooth patterns. Adaptive sampling (SkeTen) achieves the best temporal pattern recovery, while fixed-rate is slightly worse than SkeTen. Both adaptive and fixed-rate sampling are better than random sampling method, indicating the necessity of carefully design the sampling strategies.  10 shows the FMS and TCS with the change of dropped data amount from 50% to 90%. Consistent with Fig. 2 (c), (d) in the paper which illustrate the FMS and TCS on data collected on daily basis, they develop the same trend as FMS decrease and TCS increase with more data dropped. Adaptive sampling technique also shows advantageous over the fixed rate sampling and random sampling strategies. Fig. 11 compares SkeSmooth with other state-ofthe-art baselines. Results demonstrates that SkeSmooth is more effective and generally has lower recovery errors compared with other baseline methods. The results for both daily basis and hourly basis data collection strategy suggests that our proposed framework SkeTenSmooth is effective with different kinds of data collection granularity, showing the generalizability of SkeTenSmooth.