Learning Gaussian Mixture Representations for Tensor Time Series Forecasting

1 Introduction

Tensor time series (TTS) data occur as sequences of multidimensional arrays [Rogers et al., 2013; Jing et al., 2021]. For example, urban public transportation management involves data coming frommultiple times, locations, and transportation modes, i.e., Taxi/Bike Inflow/Outflow [Zhang et al., 2016; Zhang et al., 2017], based on which compelling predictions are crucial to subsequent decision-making, especially for those urban-sensitive cases. Compared with multivariate time series (MTS), TTS involves an additional variable apart from the initial time and location variables, namely the source variable. The most significant difference between MTS and TTS is that MTS involves data in a matrix, where the rows represent the time variable and the columns define the location variable. In contrast, TTS concerns tensor data, which usually takes on more substantial heterogeneity and dynamics. Although identifying and extracting similar patterns is reasonable and typical for time series prediction, the components over variables can be dynamic in many real-world applications. Failure to explicitly obtain the quantitative form of patterns may mistakenly cause the model to capture the pattern change in TTS. We use an example of NYC traffic demand data in Figure 1 to illustrate it. We collect the historical data of the three frames/locations/sources and plot their data distribution in the upper part of Figure 1. Based on the characteristics (the mean and variance) of the data distributions, we can trivially draw a conclusion that the instance at 6 am is more similar to 12 am than 8 pm. This is consistent with the evolving patterns of observation in the lower part of Figure 1. Such distinguishability over the time variable also holds water over the location and source variables. Given the complex patterns, capturing heterogeneities arising from multiple variables is an essential problem in TTS forecasting.

∗Corresponding Author

Deep forecasting models have achieved great success in more effective time series forecasting. Graph neural network (GNN) [Wu et al., 2020b; Shao et al., 2022], temporal convolution networks (TCN) [Wu et al., 2020b], recurrent neural networks (RNN) [Lai et al., 2018], and attention mechanisms [Wu et al., 2021; Zheng et al., 2020] are well-used for prediction problems. Each individual or a combination of these manipulations performs well in modeling the characteristics of a single domain or spatial-temporal regularities. Such methods are more prone to focus on low-level hierarchical information and are blind to the rich information and challenges brought by source variables. An increasing number of recent works are drawing forth the evolution from other modes as an auxiliary input, adding ancillary information to borrow similar historical patterns [Wu et al., 2020a; Fang et al., 2021; Han et al., 2021]. These practices, to a certain degree, bring mode-wise awareness. Nevertheless, they mainly help with the spatial-temporal parts, and the heterogeneity or the interactions brought by variables (i.e., time, location, and source) are still underexplored. Recently, cluster-wise learning breaks the independence between instances by exploiting the latent cluster information among samples. The learned representations are expected to retain higher-level heterogeneous information by taking advantage of additional prior information brought by clustering. However, existing approaches [Zhu et al., 2021; Duan et al., 2022] (i) neglect the dynamic nature of time series, assuming that both the representations and the members of the clusters will remain the same over time and (ii) the interactions over multiple variables as well, which are not applicable for intricate and complex TTS. Thus far, while evolution regularities in time series have been studied intensively, dynamic heterogeneity and interactions of multiple variables in TTS have yet to be adequately tackled.

Therefore, we are motivated to propose a novel TTS learning framework for accurate forecasting, namely Gaussian Mixture Representation Learning (GMRL). It consists of three steps: (i) incorporating the multiple variables information into the model through a TTS embedding (TTSE); (ii) designing a dynamic Gaussian mixture representation extractor to explicitly represent the complicated heterogeneous components among multiple variables; (iii) leveraging techniques of memory that enhances the learned representation to learn to distinguish and generalize to diverse scenarios. The contributions of this work are summarized as follows:

• We design a Gaussian Mixture Representation Extractor (GMRE), to explicitly disentangle and quantify the heterogeneity components in the time, the location, and the source variables.
• We propose a Hidden Representation Augmenter (HRA) to equip the learned representations with sequencespecific global patterns on the fly, intrinsically improving model adaptability to TTS.
• We conduct thorough experiments on two real-world datasets, one is traffic dataset containing four sources of transportation demands (i.e., Taxi/Bike Inflow/Outflow), another is air quality dataset containing three sources of pollutants (i.e., PM2.5, PM10, SO2). The results verify the superiority of our method over the state-of-the-arts.

2 Preliminaries

In this section, we formulate the tensor time series forecasting problem, and introduce the used notations (in Table 1). Problem 1 (Tensor Time Series Forecasting). Given T consecutive time steps, L discretized locations, and S sources, the multi-source spatio-temporal data can be aggregated into a Tensor Time Series X ∈ RT×L×S . Given X, the tensor time series forecasting aims to forecast the next O steps as ˆ Y ∈ RO×L×S , which can be expressed by the following conditional distribution:

3 Methodology

In this section, we elaborate on the technical details of the proposed framework Gaussian Mixture Representation Learning (GMRL), as demonstrated in Figure 2.

3.1 Tensor Time Series Embedding

Since the evolution patterns of TTS are constrained by the location and the source variables, it is crucial to incorporate the variables' information into the predictive model. Specifically, we generate location embedding and source embedding to encode corresponding variables into vectors that preserve static information, respectively expressed as eLl ∈ Rdz and eSs ∈ Rdz . The above embedding vectors only provide static representations, which could not represent the dynamic correlations among locations and sources. Thus, we further utilize temporal embedding eTt ∈ Rdz to encode each historical time step into a vector. To obtain the variable representations, we fuse the aforementioned embedding vectors in an additive manner. For location l and source s at time step t, the embedding is defined as: et,l,s = eTt + eLl + eSs . Then we get the TTSE (E ∈ RT×L×S×dz ) of L locations and S sources in T time steps. Next, the initial representation fed into subsequent modules is H = [f(X), E], where f(·) represents linear transformation and [, ] denotes concatenation operation.

3.2 GaussianMixture Representation Extractor

The canonical cluster analysis method is used to reveal the underlying nature and patterns of data, enabling the classification of the internal correlation structure of data. However, existing time series modeling methods rely on cluster distinctions that do not necessarily indicate an optimal solution. Additionally, time series data feature a dynamic context, implying that the characteristics of the clusters, which describe real-time time series data, will change over time as a function of the dynamic context.

To overcome the above limitations, we develop a dynamic Gaussian Mixture Representation Extractor (GMRE) to explicitly capture the patterns of time series data based on the underlying context. Moreover, we expect different channels to hold information about distinct contextual attributes, i.e., the factors of variation, exhibiting mutual independence and complementarity. For instance, some channels correspond to spatial identities, remaining invariant over time, while others correspond to temporal identities, displaying variations, potentially with different frequencies, over time. Motivated by this, we perform channel-wise clustering instead of the conventional sample-wise clustering, in order to separate the clusters for each channel. A practical advantage of injecting inter-channel independence lies in reducing the number of learnable parameters involved in the model, as modeling interactions among the channels is no longer necessary.

Given a representation H ∈ RT×L×S×dk derived from an input tensor X, GMRE aims to learn a mixture of Gaussian distributions over each scalar value within H. Technically, the probability distribution is conditioned on the observation (i.e., input tensor). However, for simplicity, we omit this condition in the following probability formulas. Firstly, we introduce a latent variable z[i] ∈ [1,K] to identify the membership (belonging to which of the K clusters) of each scalar value for each channel i ∈ [1, dk]. Given h ∈ Rdk (one of the T × L × S small cubes in Figure 2) from the representation H, GMRE infers the probability distribution for each scalar value in h[i] (i.e., i-th channel of h). Formally, the probability of h[i] can be decomposed as follows:

Then, we can specify the expressions for P(z[i] = k) and P(h[i]|z[i] = k). Formally, we obtain the prior distribution of z as follows:

where the membership score associated with the ith channel belonging to the kth cluster is parameterized by a linear mapping with W α k,i ∈ RTLS . For P(h[i]|z[i] = k), we assume it follows a Gaussian distribution N(h[i]|µk,i, σ2 k,i ), where the mean µk,i ∈ R and variance σ2 k,i ∈ R are generated by:

where W µ k,i , W σ k,i ∈ RTLS are weights to project representation H[..., i] to µk,i and σ2 k,i , and bµk,i , bσk,i ∈ R are the corresponding biases. exp(·) ensures that σ2 k,i remains greater than zero.

With P(z[i] = k) and P(h|z[i] = k) available, we can derive the posterior probability P(z[i] = k|h[i]) as follows, using Bayes' theorem:

Following Eq. 5, we can assign the ith channel of h to cluster k = arg max k P (z[i] = k|h[i]).

Cluster Normalization. Inspired by the simplex and efficiency of normalization technology, we devise Cluster Norm to increase the distance between clusters while reducing the distance between samples within a cluster. Cluster-wise normalization is applied to each sample belonging to the kth cluster to get a set of representations as:

where ε is a small constant to preserve numerical stability. Then, we can obtain the Gaussian mixture representation H(gm) = [H, ˆH] ∈ RT×L×S×2dk .

Probability Regularization. To enable P(h[i]) to describe the empirical data properly, we must regularize the probability of empirical data under P(h[i]). Following the convention, we alternatively minimize minus log P (h[i]), equivalent to maximizing the P (h[i]).

where KL(·) denotes Kullback-Leibler divergence. However, directly maximizing the LB (in Eq. 6) may cause mode collapsing problem, which means the model cannot distinguish features of each cluster in the expected way [Zhu et al., 2021]. We import the mutual information I(z[i], h[i]) into LB [Zhao et al., 2018] to address this issue. Formally, we derive the following cluster objective to minimize:

where Q(z[i]) = EPdata(h[i])P(z[i]|h[i]) is the expectation of posteriors with respect to the empirical data distribution.

3.3 Temporal Encoder

We implement a Temporal Encoder by applying a set of standard extended 1D convolution filters to extract high-level temporal features. Specifically, our Temporal Encoder consists of three parts: (i) two Dilated Causal Convolutions (DCC) [van den Oord et al., 2016] for achieving a larger receptive field for the representation H(gm) , i.e., the receptive field size grows in a linear progression with the depth of the network and the kernel size of the filter; (ii) a Gated Linear Unit (GLU) for introducing the nonlinearity into the network; (iii) a 1×1 convolution function for capturing sequential patterns of time series data. Formally, it works as follows:

where ∗ denotes the dilated causal convolution, Wdc1 and Wdc2 are learnable convolution filters for DCC-1 (filter) and -2 (gate), and ⊙ denotes an element-wise multiplication operator. fconv(·) denotes a 1 × 1 convolution function.

3.4 Hidden Representation Augmenter

Here, we design a Hidden Representation Augmenter (HRA) to further enhance the learned representation to be aware of and adaptive to various scenarios. HRA is motivated by a line of research on memory networks ([Wang et al., 2022; Jiang et al., 2022]), which has shown promising results in capturing global patterns. The basic idea of the memory module is to learn a parameterized memory bank that caches global patterns and projects each variable to the same feature space with the memory. The memory module brings two advantages: (i) learning and storing meaningful patterns from a global perspective; (ii) utilizing the knowledge of patterns to construct global representations.

We thereby utilize the idea of a memory module and further put the final representations conditioned on a plugin memory bank M to encourage the discovery of high-level evolution prototypes, which are representations incorporating the knowledge of the time, the location, and the source variables. Specifically, the memory blank M ∈ Rm×dm starts by taking the skip connection representation H(sc) ∈ RL×S×dk as a query, where m and dm denote the total number of memory records and dimension of each one. The memory latent representation H(me) can be calculated by:

where Q ∈ Rdm indicates the query vector projected from flattened H(sc) , ϕj is the attention score corresponding to j th memory record, and V ∈ Rdm is the reconstructed prototype representation. Then, our HRA generates an augmented representation H(aug) = [H(sc) ,H(me) ] ∈ RL×S×(dk+dm) by concatenating H(me) ∈ RL×S×dm with H(sc) .

3.5 Predictor

The final predictor consists of two fully connected layers with ReLU as an activation function, transforming the dimension of the hidden channel to the desired output dimension:

where Wo1 ∈ R(dk+dm)×dk , Wo2 ∈ Rdk×O , bo1 ∈ Rdk , and bo2 ∈ RO are filter weights and biases. relu(·) denotes ReLU function that is applied element-wise. ˆ Y ∈ RO×L×S indicates the forecasting result of the next O steps for L locations and S sources.

Optimization. In the learning phase, the Mean Squared Error (MSE) is used as the regression loss:

Eq. 7 and Eq. 11 are combined as the loss function for our TTS forecasting model:

where λ ≥ 0 is the balancing parameter.

4 Experiments

4.1 Experiment Setup

Datasets. We conduct experiments on two real-world TTS datasets, namely NYC Traffic Demand and BJ Air Quality as listed in Table 2, the details of which are as follows:

• NYC Traffic Demand dataset 12 is collected from NYC Bike Sharing System, which consists of 98 locations and four sources: Bike Inflow, Bike Outflow, Taxi Inflow, and Taxi Outflow.
• BJ Air Quality dataset 3 is collected from the Beijing Municipal Environmental Monitoring Center, which contains 10 locations and three pollutant sources: PM2.5, PM10, and SO2.

Settings. We implement the network with the Pytorch toolkit. For the model, the number of GMRE-TE layers and cluster components K are set to 4 and 17. The kernel size of each dilated causal convolution component is 2, and the related expansion rate is {2, 4, 8, 16} in each GMRE-TE layer. This enables our model to handle the 16 input steps. The dimension of hidden channels dz is 24. The parameters for memory bank m and dm are set to 8 and 48. The batch size is 8, and the learning rate of the Adam optimizer is 0.0001. In addition, the inputs are normalized by Z-Score.

Baselines. To quantitatively evaluate the prediction accuracy of our model, we implement eight baselines for comparison.

• LSTNet [Lai et al., 2018]. A deep neural network that combines convolutional neural networks and recurrent neural networks.

1 https://ride.citibikenyc.com/system-data

2 https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page

3 https://archive.ics.uci.edu/ml/datasets/Beijing+Multi-Site+

Air-Quality+Data

• MTGNN [Wu et al., 2020b]. A spatial-temporal graph convolutional network that combines graph convolution and dilated convolution.
• Graph Wavenet (GWN) [Wu et al., 2019]. A spatialtemporal network that combines adaptive graph convolutions with dilated casual convolution.
• AGCRN [Bai et al., 2020]. A recurrent neural network with adaptive graph convolution.
• StemGNN [Cao et al., 2020].A spectral temporal graph neural network that maps temp-oral-spatial domain to spectral domain.
• STtrans [Wu et al., 2020a]. A deep neural network based on attention mechanism that integrates dynamic spatial, temporal and semantic dependencies.
• ST-Norm [Deng et al., 2021]. A normalization-based approach that refines the temporal components and the spatial components.
• MiST [Huang et al., 2019]. A co-predictive model for multi-categorical abnormal events forecasting, that models the spatial and categorical embeddings incorporated.

4.2 Overall Performance

We evaluate the performance of our proposed model as well as the above baselines for tensor time series forecasting on all two datasets. We repeat the experiment five times for each model on each dataset and report the average results. Through Table 3, we can find our model outperformed the state-of-the-arts to a large degree. For all four sources of the NYC Traffic Demand dataset, there have improvements of ∆MAE{14.88%, 15.10%, 11.14%, 12.61%}, and ∆RMSE{12.45%, 12.78%, 7.94%, 11.02%} on three horizons average. Among the co-prediction models, MiST performed better than STtrans. However, MiST is second only to MTGNN in terms of bike demand and much inferior in taxi demand. Because it requires reference sources as auxiliary input while outputting the target source prediction, it is not adequate to quantify heterogeneity on three variables (i.e., the time, the location, and the source variable). By contrast, StemGNN and ST-Norm gave a worse performance because ignoring the role of source variables too much makes them not applicable for tensor time series forecasting.

The comparison results with baselines on the BJ Air Quality dataset are shown in Table 4. GMRL has the improvements of ∆MAE{1.87%, 0.90%, 1.45%}, and ∆RMSE{1.05%, 1.66%, 7.94%} for three sources on three horizons average. The BJ Air Quality dataset only has 10 locations and less corrections, and many weather-related features are unavailable. This causes GMRL to perform less well on data with insufficient information, but it fully exploits the correlations among available features. In this case, baseline models are still not as efficient as GMRL.

4.3 Ablation Study

To make a thorough evaluation of key components in our model, we create several variants as follows:

• w/o GMRE. It removes GMRE in each GMRE-TE layer and keeps other components of GMRL, i.e., TTSE is directly used as input to Temporal Encoder.
• w/ KMeans. It applies the KMeans clustering algorithm instead of GMRE to obtain a cluster-wise representation, which randomly initializes K cluster centers and calculates the categories according to the distance between the samples and the centers.
• w/o HRA. It removes HRA in GMRL, i.e., H(aug) = H(sc) , which means the representation H(sc) is directly fed into Predictor to generate the predictions.
• w/o L(cluster) . It removes L(cluster) in GMRE (Eq. 7), i.e., L = L(reg) in the learning phase.

Through Table 5, we can see that compared to w/o GMRE, w/ KMeans will deteriorate the forecasting performance, because it only assembles instances by hard clustering, and the clustering result is affected by the selection of the initial cluster center. Removing HRA and L(cluster) deteriorate performance to a lesser degree. But still, there is a significant gap compared with the full GMRL. All these demonstrate that the proposed GMRL is a complete and indivisible set.

4.4 Case Study

Hidden Representations before/after Cluster Norm

In this subsection, we conduct multiple studies to qualitatively analyze the representations learned from GMRE of the last GMRE-TE layer. In Figure 3, we visualize variations of clusters in different channels before and after applying Cluster Norm. In Figure 3a, the evolution curves of clusters 2, 13, and 14 in Channel 1 overlap and intertwine with each other, especially for clusters 13 and 14. This is because clusters 13 and 14 both originate from the same source (Bike Inflow), and the overlapping features correspond to several identical locations. As shown in Figure 3b, after Cluster Norm, the evolution curves of clusters in Channel 1 become distinguishable. The clusters corresponding to spatial identifiers maintain a stable behavior over time, which remains unchanged even after Cluster Norm. The performance of clusters in Channel 8 before and after Cluster Norm is similar to that in Channel 1, as shown in Figure 3c and Figure 3d. However, in Channel 8, clusters 2, 9, and 16 represent temporal identifiers, and they exhibit changes with varying frequencies over time compared to Channel 1. Clearly, after applying Cluster Norm, each cluster in different channels exhibits heterogeneous evolution with sufficient discriminability. This study confirms the component clustering capability and time-adaptability of GMRL.

Gaussian Mixture Parameters

In Figure 4, we display multiple Gaussian mixture parameters learned in GMRE. As the GMRE-TE layer goes deeper, Gaussian mixture parameters α, µ, and σ2 can perfectly describe six heterogeneous components after passing through the last GMRE-TE layer. For example, cluster components 4 (purple) and 5 (brown) tend to be similar on parameters α and µ, but their variances are very different. When the models converge, the Gaussian distribution they determine is sufficient to represent heterogeneous components.

4.5 Hyper-parameter Study

We further study the effect of key hyper-parameters manually set in the proposed GMRL, including the number of cluster components K, the dimension of the hidden channel dz, and the balance parameter λ, as shown in Figure 5. We find that as K and λ increase, the prediction error of GMRL first decreases and then increases. This phenomenon follows our understanding that the model can suffer from underfitting or overfitting due to improper settings of cluster components and balance parameter. Surprisingly, the dimension of the hidden channel has little effect on GMRL: setting it to 24 can achieve fairly good performance, and increasing the dimension of the hidden channel hardly changes the performance.

5 RelatedWork

Time series forecasting has been studied for decades and applied to various fields. Some studies utilized the attention mechanism for dealing with lost memories [Li et al., 2019; Vaswani et al., 2017; Zhou et al., 2021]. To extract complex spatial-temporal patterns, existing works applied various operations over temporal and spatial domains. Specifically, the graph neural network is developed to model the spatial relation [Yu et al., 2017; Li et al., 2017; Oreshkin et al., 2021; Shang et al., 2021; Chen et al., 2020]; the attention mechanism and its variants contain spatial attention [Fang et al., 2019; Zheng et al., 2020] and temporal attention [Wu et al., 2021; Zheng et al., 2020; Liu et al., 2021; Zhou et al., 2022]; the convolution operator is leveraged to spatial convolution [Deng et al., 2022b; Guo et al., 2021], temporal convolution [Wu et al., 2020b; Wu et al., 2019], spatial-temporal convolution [Guo et al., 2019; Yang et al., 2021] and adaptive convolution [Pan et al., 2019]. [Zhu et al., 2021] proposes an end2end mixture model to cluster microscopic series for macroscopic time series forecasting. [Jiang et al., 2021] reviews the deep learning models and builds a standard benchmark. [Deng et al., 2021; Deng et al., 2022a] refine the high-frequency and local components from MTS data by using normalization technology. Unlike MTS involved data within a matrix, data points within TTS usually take on more substantial heterogeneity and dynamics. Since such unique properties of TTS compared to MTS, traditional methods for predicting MTS (or MTS with reference source) are not entirely suitable for TTS prediction. The inability to explicitly learn and encode interactions from TTS brings up some intuitive concerns. A new line of research further extends the spatialtemporal dependency modeling by considering another domain (i.e., source or categorical domain) [Ye et al., 2019; Huang et al., 2019; Huang et al., 2021; Han et al., 2021; Wang et al., 2021] or extra data source [Jiang et al., 2023]. Our work distinguishes itself from these methods in the following aspects: (i) our GMRE is dynamic, which could explicitly represent the complicated heterogeneous components among multiple variables in TTS; (ii) our HRA enhances the learned representations with sequence-specific global patterns to distinguish and generalize to diverse scenarios, which could help model adaptability to TTS.

6 Conclusion

In this paper, we develop a dynamic Gaussian mixture representations learning framework GMRL, which seeks to individually model each heterogeneity component implied in multiple variables. In addition, we design an HRA module that enhances the learned representation with sequencespecific global patterns to distinguish and generalize to diverse scenarios. Extensive experiments on two real-world datasets verify our perspective's insightfulness and the proposed approach's effectiveness. In the next step, we plan to improve the speed of adaptation and efficiency of GMRL.

Contribution Statement

Renhe Jiang is the corresponding author for this study.