Rotation Invariant Networks for Image Classiﬁcation for HPC and Embedded Systems

: Convolutional Neural Network (CNNs) models’ size reduction has recently gained interest due to several advantages: energy cost reduction, embedded devices, and multi-core interfaces. One possible way to achieve model reduction is the usage of Rotation-invariant Convolutional Neural Networks because of the possibility of avoiding data augmentation techniques. In this work, we present the next step to obtain a general solution to endowing CNN architectures with the capability of classifying rotated objects and predicting the rotation angle without data-augmentation techniques. The principle consists of the concatenation of a representation mapping transforming rotation to translation and a shared weights predictor. This solution has the advantage of admitting different combinations of various basic, existing blocks. We present results obtained using a Gabor-ﬁlter bank and a ResNet feature backbone compared to previous other solutions. We also present the possibility to select between parallelizing the network in several threads for energy-aware High Performance Computing (HPC) applications or reducing the memory footprint for embedded systems. We obtain a competitive error rate on classifying rotated MNIST and outperform existing state-of-the-art results on CIFAR-10 when trained on up-right examples and validated on random orientations.

Introduction
In recent years, Convolutional Neural Networks (CNNs) have the highest accuracy rate on state-of-the-art image classification tasks [1][2][3]. Usually, to obtain these results, the networks are trained for a couple of days or weeks. This training time comes from feeding the network with different image transformations to increase the generalization potential (data-augmentation technique).
As the number of transformations increases, the training time increases consequently. It has been tested and demonstrated that a data augmented model has up to 600% training part than the non-augmented model [4]. Hence, the usage of data augmentation techniques results in increased energy consumption Recently, the increased training time of models in terms of millions of parameters has become a concern. While these models' accuracy is off-charts compared to traditional approaches, better usage of the network's internal resources would result in more energyefficient algorithms. The best scenario would be to achieve the same accuracy level while minimizing parameters, reducing training time and energy costs. This scenario would lead to faster inference times and reduced environmental costs associated with the network's training time.
One possible approach to reduce the number of parameters needed is to endow the network with properties that it usually obtains via data augmentation. It is common knowledge that CNNs are translation invariant up to some extent. This property states that the network can find the features for image classification (i.e., a cat or dog's ear) despite the feature's position on the image. Unfortunately, this invariance property is not present

Related Works
As mentioned previously, a low-footprint network allows reducing the energy spent on training and inference. Thus, the main challenge is to design a network invariant to the image transformations allowing to preserve the accuracy of state-of-the-art networks while keeping the number of trainable parameters low. In this sense, here, we present a state-of-the-art survey of works involved in rotation-invariant networks.
In the literature, two main strategies exist to endow a network with rotation invariance: (i) modify the input image by several images (data-augmentation) (ii) modification of the inner filters of the network. The second one has recently gained attention with different approaches: continuous rotation approaches and discrete angular sampling of the filters. As the main scope of this work is on reducing neural networks, we will present results in the (ii) group, as data-augmentation has demonstrated to increase training time and network size.
Recent works in state-of-the-art converge in providing rotation invariance properties by encoding it in the neural networks' internal filters. Typically, some of these approaches try to find a trade-off between the number of features needed to get a good prediction accuracy and network size.
These approaches modify the network's inner filters in some way (i.e., rotation) to decompose the input image into several orientations. The act of sampling the input allows some of these approaches to obtain rotation-invariant features; hence, the network acquires rotation invariant properties. The main challenge of these works is to find a trade-off between the angular sampling (fine-grained sampling needs more memory and network size) to obtain good accuracy and a reasonable network size.

Continuous Filter Sampling
In the continuous sampling, we can find Worral et al.'s work that achieves a good tradeoff between angular sampling and size in their Harmonic Networks [14] work. The authors use circular harmonics, returning a maximum response, and they achieve continuous angular sampling by constructing any angular filter using the linear combination of a set of base filters.
Other approaches use a similar Harmonic approach to obtain continuous sampling; Cohen et al. [15] use a spherical cross-correlation in their Spherical CNN work. They introduce a set of building blocks to obtain a rotation-equivariant cross-correlation in spherical images. Furthermore, they present results on up-right training and rotated samples validation; however, the computational cost of the spherical space remains very high. Cohen et al. address the reduction of the computational cost in their Icosahedral CNN [16] by sub-sampling the spherical space; the results achieve better computational cost than the previous approach, but the rotation-equivariance is affected.
Continuous sampling approaches reach competitive accuracy and rotation invariant properties. The main drawback of these approaches is the computation involved to obtain continuous sampling. For example, Worral et al. need to add a process to maintain the rotation order's disentanglement that is a product of the combination of complex features [14]. This process leads to increased network size, training time, and network complexity.

Discrete Angular Sampling
The counterpart of the continuous approach is to use a finite number of angular samples. Works in this group sample the input by several angles that become a hyperparameter of the network. Most of these approaches use a fixed angular sampling adequate to the application and the dataset. For example, Cohen et al. [17] find a trade-off between angular sampling and network size. They use a symmetry group composed of 90-degree rotations and pool over the group. Pooling in the first layers discards important angular information, as explained in their paper. Usually, using multiples of 90-degree orientation steps avoid interpolation problems. The main problem is that using only 90 degrees multiples, the network cannot generalize the data orientation and predict smaller angles.
Several previous works use the oriented filter decomposition with group convolutionbased predictors. To mention some, Sifre and Mallat introduced the Scattering transformbased networks in [18]. They use Scattering transforms to generate rotation-equivariant feature extractors. Marcos at al. [19] introduce a CNN architecture encoding rotation equivariance, called Rotation Equivariant Vector Field Networks (RotEqNet). The network applies one filter at different orientations and extracts a vector field feature map, encoding the maximum activation in terms of magnitude and angle. Zhou et al. [20] propose the Oriented Response Networks (ORN) to have Active Rotating Filters. These filters rotate during convolution and produce maps with location and orientation explicitly encoded.
An extension to ORN is the Rotation Invariant Local Binary Networks [21]. The authors replace the first convolutional layer for steerable modules that can be inserted into traditional CNNs. In RotDCF [22], the authors demonstrate the benefits of decomposing the convolution filters leading to implicit regularization and robust representations; they present results on rotated samples validation up to 60°. The main limitation of these approaches is the high number of parameters to learn, the computational complexity, and the fixed number of angular sampling they present.
Follman et al. [23] achieve rotational invariance by rotating the filters and then backrotating the generated feature map by the negative angle of the filter. Weiler et al. [24] introduces the rotation-equivariant CNN with learnable steerable filters. Their SF-CNN allows choosing an arbitrary angular sampling. The idea of learning steerable filters inspired layers based on the Gabor filter bank [25]. Whereas the obtained accuracy is very competitive, the authors present results only for four different orientations, and the number of parameters increases to almost 2 million parameters on the MNIST dataset. While these approaches introduce a lower sampling complexity, they do not provide the roto-translational space to achieve equivariance or angular prediction.
Oyallon and Mallat [26] first present the roto-translational feature space; in their work, they use a non-trainable Scattering transform oriented feature space. Nevertheless, they do not use these oriented feature space properties to obtain an angular prediction. To the best of our knowledge, the first work that uses a translating predictor over an oriented feature space is from Rodriguez et al. [27]. In their work, they use the Scattering transform for the orientation filter decomposition. They apply a predictor to the roto-translational feature space to make the input's angular prediction. The accuracy of both of these works is deeply affected by the non-trainability of the Scattering filters involved in the oriented filter decomposition stage.
Our prediction model uses a set of features generated by one mother filter. These features, oriented and ordered, form a representation space where translation covaries with the input rotation. A subsequent predictor scans over this representation space and predicts the class and angle of the input. As the rotation invariance is the product of the feature space's roto-translational properties and not of the translating predictor, the shared weight predictor keeps a small footprint in the network. Using this approach, the user can select the angular sampling as a network hyper-parameter to adapt the angular resolution to the user's needs.
Furthermore, in the literature we can find references that study the vulnerability of networks based on sampling techniques [28]. In the case of angular sampling techniques, this would mean losing orientation information caused by the attack. The worst case would result in the network's accuracy reduction or the loss of rotation invariance properties. The oriented sampling and shared weight predictor presented in this work make a possible solution against these attacks as the predictor learns from each translation and becomes rotation sensitive. If an oriented wavelet is missing, the predictor has enough information from the complete filter ensemble and the other translations to complete the inference.
In the following section, we present the methodology in which we base our rotation invariant network. We detail how we obtain a roto-translational feature space and how the shared weight predictor is applied to the feature space to obtain the class and angle prediction.

Methodology
In this work, we propose a convolutional network architecture ( Figure 1) with rotationinvariant properties and prediction of the angle without angle labels. With rotationinvariant properties, we refer to the network's capability to correctly predict the class despite the orientation of the object in the input. Furthermore, it can predict the angle by fostering the angle information from the relative position of the network's inner filters.
For easier understanding, we divide the intuition behind the network into two parts. First, the generation of a feature space that translates with respect to the rotation of the input image (roto-translational feature space) (Figure 1 upper row). Second, a predictor that scans the generated roto-translational feature space and outputs a prediction for each of the translations (Figure 1 lower row). The result of these two processes is a probability distribution that contains a prediction for each translation. Then, the class is obtained from the prediction, and the angle from the relative position of the maximum predicted class.

Oriented Representation Mapping
Oriented representation mapping is the result of two operations: oriented components, based on wavelet theory [29], and re-orientation of such components. First, we start by obtaining a set of oriented components from the input image x (Figure 2a). An oriented component is obtained by the product x × g ϕ where ϕ is the orientation of an oriented edge detector g. Furthermore, we have a set of orientations such as i = 0, . . . , N, with N ∈ Z + where N is the number of orientations in which we want to decompose the input image. This means that the set [x × g ϕ i ] contains N oriented components (features) of the input image x. Furthermore, we can define dϕ = 2π/N as the angular sampling magnitude in degrees for a filter with periodicity 2π. Then the orientation of the filter in the position i position is ϕ i = 2πi/N.
Is straightforward to see that the input image x of size [m, n] is decomposed in N oriented components of the same size generating a feature space with size [m, n, N] (Figure 2b). This can be observed in Figure 2b where each feature is an oriented component in the angle ϕ i .
Second, let ρ ϕ be a transformation rotating the support by the angle ϕ. Then we compensate the rotation by rotating the support −ϕ i degrees. This re-orientation results in all the feature components to be aligned to the referential orientation ϕ 0 (Figure 2c).
Then we can define Φ(x) as an oriented feature representation of x: containing oriented components of x, and re-oriented to align with the referential orientation ϕ 0 . As a consequence of the oriented decomposition of the input (x × g ϕ i ) and the reorientation (ρ −ϕ i ) the oriented features in Φ(x) translate i times following idϕ magnitude ( Figure 3). Experimental results demonstrate that the transformation ρ −ϕ i is necessary to obtain a roto-translational feature space. Without this transformation the oriented feature space does not translate and the prediction accuracy is heavily affected.
The oriented components of the input can be obtained with several approaches. Previous works [27,30] used the Scattering transform and Steerable filters. We propose the usage of trainable Gabor filters that can offer a higher degree of freedom to fit the data. The implementation aspects of trainable Gabor filters are explained in the next section. Furthermore, in the result section we compare the Steerable filters approach with the Gabor filters generation of the oriented features.

Gabor Filters
Gabor filters offer a higher degree of flexibility than previous works approaches [27,30], especially in the sense that the angular and frequency selectivity can be learned independently; they can generate an oriented feature space and trained to fit the training data. A 2D Gabor filter can be defined as an oriented sine (imaginary part of the filter) or cosine (real part) modulated by a 2D Gaussian function [31]. The shape of the Gabor function is defined by the frequency, the orientation of the sinusoid, and the scale of the Gaussian function [32].
Consider the rotation matrix in 2D space: Let σ x and σ y divide the rotation matrix by their magnitude: We obtain the respective real and imaginary components of the Gabor function in which σ x corresponds to the filter's angular selectivity and σ y to the frequency selectivity of the filter.
The parameters σ x , σ y , f become trainable parameters of the network and are updated by the back-propagation algorithm ( Figure 4). Recalling ϕ i as a list of N orientations, Gabor filters can generate an oriented filter g separated in real and imaginary components.
The presented Gabor methodology then generates an oriented filter bank g ϕ i with a set of orientations ϕ i that allows the extraction of the oriented components of the input by the convolution operation x × g ϕ (Equation (1)). We implement this operation as a 2D Convolution with as many Gabor filters as orientations N. Each one of these filters is oriented as ϕ i and the shape is determined by the Equation (4) (for the real part filter bank) and Equation (5) (for the imaginary part).
The wavelet re-orientation stage (ρ −ϕ i ) is then applied to these oriented components to obtain the oriented components re-oriented to the referential orientation ϕ 0 . This means that each filter g ϕ of the bank is rotated by a fixed angle −ϕ i calculated by its position i on the bank using bilinear interpolation. As consequence of this Gabor methodology we obtain a roto-translational feature space (Figure 3). The number of oriented components is then defined as a network hyper-parameter, hence the user selects the number of oriented components N that conform the oriented feature space. In this work, we present different N values and results with N = 16.

Feature Extraction
The most typical feature extraction stage in the literature is stacked convolutions with increasing filters over the depth. Usually, these convolutions are followed up by a batch normalization operation and activation. The feature extraction stage has recently moved to architectures that contain connection skipping and the depth of hundreds of layers. Following the state-of-the-art implementations, we connect the feature extraction stage to the representation mapping's oriented feature representation.
This feature extraction stage's main challenge is to preserve the roto-translational feature properties of the oriented feature space. As the translation occurs in the extra dimension of size N, we apply the feature extraction to the feature space dimensions that do not have translation [m, n]. This means that input x gets decomposed in oriented components by a set of oriented filters (previous subsection) and generates a roto-translational feature space with shape(m, n, N, channels) where shape refers to the dimensions of the space. To preserve the translation axis, then the feature extraction is applied to the m, n, and channels axis (Figure 1c). A straightforward way to obtain this is with a 3D Convolution with a window size of one in the translation axis, for example, with shape(3, 3, 1) for a square kernel support of size 3. This implementation way allows the convolutional backbone to extract and refine the feature space's useful features in the same way as the feature extraction stages from the state-of-the-art architectures.
Gabor filters output two channels, the product of the real filter response and the imaginary filters response. Furthermore, it can be possible for the Gabor filters approach to use more than one mother filter and generate as many pairs of real and imaginary filters as necessary. Nevertheless, in this work, we have used a single mother filter in all the experiments.
If observed, this feature extraction convolutional backbone is applied to a 3-axis tensor shape(height, weight, channels). The current state-of-the-art networks use a similar shape as the network's input. This opens the possibility of changing the convolutional backbone to a known state-of-the-art architecture such as ResNet. While the convolutional backbone reaches state-of-the-art results on the MNIST dataset, a more robust backbone is needed for more challenging datasets such as CIFAR-10. To prove this, we propose applying a ResNet backbone to the oriented feature space to extract and refine the features. The translation is then preserved on the additional axis as the backbones' convolution window does not affect this dimension.

State-of-the-Art Architectures as Backbone (ResNet)
The output shape of the Gabor filter approach opened a new possibility to explore: exploring the usage of state-of-the-art backbones (i.e., ResNet). If the connections are managed correctly, it opens the possibility to use a known network as a plug and play backbone of the rotation invariant network. In other words, it is possible to use these networks potential to refine the features ( Figure 5). The first challenge is to adapt the shapes and connections between the custom Gabor layer and the possible backbones. Usually, networks have an input with shape(width, height, channels). This means that the expected input of these networks are not prepared for the dimension N containing the orientations product of the Gabor layer. One way to achieve this shape adaptation is to use the backbone predictor on each one of the filters (apply it on one orientation instead of the complete translation of size N) (Figure 5c). Then, concatenate the outputs and apply a translating convolution hence obtaining the same behavior of having a prediction for each position on the roto-translation feature space.
Several questions could arise when following this methodology. First, the rototranslational properties of the feature space are preserved on the architecture. This is possible because the backbone network is applied to each filter individually and the position of the filter is preserved in the output of this layer.
Another concern comes in the form of the size of the network. As we use a backbone network, the number of trainable parameters increases quickly. For example, using a ResNet20 with 438 k parameters multiplicated by 16 orientations would mean about 7 M trainable parameters. Sharing weights allows the backbone to learn from the different orientations and be implemented once in memory hence keeping the network size constant for any number of orientations N.
This means that the network could be benefited from the usage of any state-of-the-art backbone network instead of the convolutional backbone presented previously. With this the network is able to predict class and the angle without any angular labeling on the training stage. As the network with this backbone is complex enough, we opted for testing ResNet backbones on the CIFAR-10 dataset.

Translating Predictor
As a result of the input decomposition presented in the last section, we obtain an oriented feature representation space. This oriented feature space contains ϕ features ordered in increasing magnitude. As a consequence of this transformation, the feature space translates over this depth axis of size ϕ in a cyclic form.
When we rotate the input x the oriented feature is translated i spaces over the N axis. The translation of the filters is proportional to the magnitude of the input's angle and covariates following idϕ. This means that for any input rotation a translation i exists in the oriented feature space up to ±dϕ. Then is possible to find a translation for an almost unrotated version of the input.
It is possible to find the position i corresponding to the almost unrotated version of the input using a cyclic convolutional predictor over the feature space. This predictor with an attention window of size N reads each translation τ i and obtains a prediction for each position i. Hence, the predictor will output a higher probability on the translation i that corresponds to the unrotated input (the same orientation of the training samples). Consequently, the linear distance of the position i where the maximum probability is obtained to the reference orientation becomes the input angle idϕ.
A straightforward way to implement this cyclic convolutional predictor is to pad the feature space by itself and then applying this convolution over the translation N axis. This is implemented as a 3D Convolution with a kernel size of (3, 3, N) that scans each one of the translations τ i . Then a fully-connected layer with size K (number of classes) generates a class prediction for each one of the translations τ i (Figure 1e). The number of classes K is determined by the problem and dataset used in training (i.e., MNIST contains ten classes K = 10 one for each digit in the dataset).
The output contains N predictions in the form of a probabilistic distribution with a higher value on the translation τ i corresponding to almost the same orientation of the input with slightly lower values at the neighbor positions of the maximum value. A max-pooling operation over the N axis and the K axis outputs the angular and class prediction, respectively (Figure 1f). This max-pooling operation over the output reinforces the convolutional shared predictor on the correct orientation while not updating it at the other orientations. This implementation model allows the proposed prediction model to have rotation-invariant properties angular prediction.
In conclusion, following this method, we obtain a convolutional neural network with rotation invariant properties. This network can predict the input class and angle. The angle is predicted using the relative position i, where the highest prediction is present. As this position is proportional to the rotation of the input (roto-translational feature space), then the angle can be fostered from the translation information (idϕ). This work presents results with two alternatives for the oriented feature space and results with two different backbones: convolutional and different sizes of ResNet. In the following section, we describe the results obtained on the MNIST dataset and CIFAR-10 dataset and compare them with the state-of-the-art approaches.

Memory-Footprint vs. Data Parallelism
One of the main advantages of this approach is the capability of selecting between memory-footprint and data parallelism. Our network uses a set of shared weights convolutions for the feature extraction and for the translating predictor explained previously. This means that the number of parameters is kept low while keeping high accuracy. A closer look at the data dependency between the layers demonstrates that each filter can be processed independently, speeding up the inference process.
For the inference step, instead of a shared weight convolution scanning each filter position, we can implement N instances of the feature extraction and translating predictor and process each position i in parallel. This would mean an increase in the memory size of N times the operation size but speed the inference time N times. An HPC center with several cores could run each orientation in a core as the memory is usually not a limitation in these devices.
On the other hand, the user could select to implement the feature extraction and translating predictor as a single instance scanning each orientation i. This would mean a low-footprint network (10 k learnable parameters) at the expense of augmenting the inference time N times.
The network flexibility between memory footprint and parallelism would allow benefits on several applications. Usually, neural networks require large amounts of onchip memories that occupy significant chip area, making impractical their realization on small footprint devices [33]. Shared weights and the presented methodology allows the implementation of the network in these limited memory embedded devices. For the parallelism, there exist some efforts on deeply pipelined multi-FPGA implementations [11] which could benefit both from the low-footprint size and to run in parallel pipelines.

Results
We evaluate different options for the proposed stages of the architecture. For the oriented feature generators we test Steerable and Gabor filters on the MNIST dataset, and for the feature extractors we test convolutional and ResNet on CIFAR-10 dataset. To test rotational-invariance we train the network with up-right samples (all the samples in the same orientation) and validate with randomly rotated samples. We train the network for 200 epochs using the Adam loss function on a Titan Xp GPU with 16 GB VRAM; unless noted the angular sampling is N = 16. This N value allows a direct comparison between our approach and state-of-the-art approaches but can be changed by the user as a network hyper-parameter.

MNIST
For the MNIST dataset we have trained the network using up-right samples and validated it on randomly rotated samples. This type of validation ensures that the rotation invariance property is adquired by the network. We present results using two different oriented feature filters using Steerable filters and Gabor filters approaches.
We present the rotation-invariant approaches in Table 1. The presented results include the cases when the network is trained on up-right samples and validated on random orientations. Some of the presented approaches obtain the rotation invariant property by increasing the number of parameters over 100 k. Usually, they are group-convolution-based approaches that lose the orientation of the input to ensure the invariance. Instead, we preserve the orientation as a different axis of the network and use it to predict the angle.
The nearest approach in parameter terms is Covariant CNN [27] with lower accuracy; this is mostly due to the static nature of their scattering filters approach instead of trainable ones like the presented on this work.  [15] 6.00% 68 k Icosahedral CNN [16] 30.01% n.c. RI-LBCNNs [21] 25.77% 390 k Covariant CNN [27] 17.21% 7 k RIN (Steerable + Convolutional) (this paper) 2.05% 42 k RIN (Gabor + Convolutional) (this paper) 1.71% 9 k As observed in Table 1 the Gabor filters approach has a lower number of parameters while achieving higher accuracy (best values are depicted in bold); this is mostly caused by the degrees of freedom of this approach. The Steerable filters need to be circular to keep their steerability properties while the Gabor ones have better angular and frequency selectivity. Our results demonstrate the capability of the added functions to endow with rotation invariance and angular prediction properties to neural networks. Furthermore, the number of parameters has been kept lower than 10 k.

CIFAR-10
We tested the capability of our approach on the CIFAR-10 dataset. We opted for Gabor Filters as oriented filters decomposition function because it showed better adaptive capabilities and have higher degree of freedom. We tested two different feature extractors: convolutional and ResNet.
For the first one, the convolutional feature extractor, we have three 2D Convolution stages with filters number increasing with the depth. Each one of the convolutional stages is followed up by a BatchNormalization layer following the state-of-the-art implementations. The convolutional feature extractor shows error rate up to 28% with 238 k learned parameters. Increasing the number of filters, hence the number of trainable parameters, shows no improvement over this point. In conclusion, higher number of filters in the convolutional feature extraction approach tend to overfit the dataset quickly.
For the ResNet feature extractor, we test different depths of the network; with ResNet20 we obtain 21% error rate and preserve the rotation invariance properties. Having more layers on this type of feature extractor did not improve the error rate further than 21%; this means that the network tend to overfit quickly with deeper feature extractors. Table 2 condenses the results with Gabor oriented filters and these two different feature extractors (best values are depicted in bold). We compare these results with previous stateof-the-art implementations on CIFAR-10 dataset when trained on up-right and validated on randomly oriented samples. In general, we observe an accuracy increase over previous methods (Figure 6a) with both approaches. For both feature extractors (Convolutional, Gabor) we find the same error rate with approximately 240 k parameters; more filters in the convolutional feature extractor did not improve the error rate. Nevertheless, ResNet improves the error rate by 6% until 21% error rate is reached. When comparing the ResNet feature extractor with the original ResNet (Figure 6b) we observe that our proposal achieves rotation-invariance having around 80% accuracy; the original ResNet20 accuracy drops when tested with non-upright samples. Increasing the number of orientations N > 16 did not result in an increase on accuracy.

Conclusions
In this paper, we propose a prediction model based on an oriented feature representation mapping and a translating predictor. The concatenation of these two allows the usual convolutional networks to obtain rotational invariance properties. The oriented features have roto-translational properties between the input and an additional axis proportional to the angular sampling of the input. The translational predictor uses the translating feature space to obtain a set of probabilities for each possible discrete orientations of the input. The maximum value of this set contains the prediction for the class and orientation.
We present results obtained with different oriented features mapping as Steerable and Gabor filters for the MNIST dataset. The mapping generated by Steerable filters demonstrates to increase the accuracy results compared to previous state-of-the-art implementations on the MNIST dataset; their main drawback is the circular shape that limits the filter flexibility to adapt to the data. Gabor filters present higher flexibility than Steerable ones; with Gabor filters, one can either prefer a better angular or frequency selectivity. Results on the MNIST dataset show a slight reduction of the error rate (Steerable: 2.05%, Gabor: 1.71%), with the Gabor implementation being almost five times smaller than the first one. These results demonstrate the limited capabilities of the Steerable filters to fit the data. Both of these approaches outperform the previous Scattering approach error rate of 17.21%.
This reduction in the model's size opens the possibility to use these networks in low-memory devices like embedded devices or smartphones. For example, aerial/satellite imagery applications on UAV could be benefited from low-footprint rotation-invariant networks opening the possibility of inference on-site. Furthermore, the reduced model represents an enhancement in energy reduction for reduced training time and less memory access on the inference. The obtention of these results without data augmentation also opens the possibility to focus on developing efficient networks and using the inner filter resources adequately. One example of this is the usage of the relative position of the filters to obtain rotation invariance properties and angular prediction without angle labels.
Furthermore, we tested the proposed prediction model with different feature extractors: convolutional and ResNet. In both cases, we use the headless architecture (feature extraction stage). To validate these feature extractors, we trained them on the CIFAR-10 dataset. Both implementations outperformed (Convolutional: 19.98%, ResNet:11.95%) previous state-of-the-art error rates when trained and validated on up-right samples. Furthermore, these implementations demonstrated to acquire the rotation-invariant properties and angular prediction capabilities described in the paper.
The low footprint of the presented work opens the possibility for implementation in embedded, mobile, and other resource-limited devices. As this model offers the flexibility between memory footprint or multi-thread parallelism, it allows to target the hardware efficiently when more cores are available (e.g., HPC centers) or memory-constrained hardware (e.g., FPGA, single-board computers). In the end, this means better resource usage and energy consumption reduction for inference. In future works, we expect to test the model with other state-of-the-art feature extraction stages and validate the rotationinvariant properties in harder datasets as faces or food. Furthermore, a deeper analysis of the multi-threading advantages and memory constraints are expected in the near future.


extfieldsofstudy:
Computer Science
Mathematics

text:
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.