# Simultaneous Incomplete Traffic Data Imputation and Similarity Pattern Discovery with Bayesian Nonparametric Tensor Decomposition

Authors: Yaxiong Han; Zhaocheng He
Journal: Journal of Advanced Transportation
 (2020)

---

## Abstract

A crucial task in traffic data analysis is similarity pattern discovery, which is of great importance to urban mobility understanding and traffic management. Recently, a wide range of methods for similarities discovery have been proposed and the basic assumption of them is that traffic data is complete. However, missing data problem is inevitable in traffic data collection process due to a variety of reasons. In this paper, we propose the Bayesian nonparametric tensor decomposition (BNPTD) to achieve incomplete traffic data imputation and similarity pattern discovery simultaneously. BNPTD is a hierarchical probabilistic model, which is comprised of Bayesian tensor decomposition and Dirichlet process mixture model. Furthermore, we develop an efficient variational inference algorithm to learn the model. Extensive experiments were conducted on a smart card dataset collected in Guangzhou, China, demonstrating the effectiveness of our methods. It should be noted that the proposed BNPTD is universal and can also be applied to other spatiotemporal traffic data.


---

## Body

## 1. Introduction

Recent advances in data acquisition technologies and mobile computing lead to a collection of large quantities of urban traffic data from various sources, such as loop detectors data, GPS data, and smart card data. These datasets can capture rich spatial-temporal information of the whole transportation system and enable some traffic analysis. A crucial task in a data-driven transportation system is similarity pattern discovery. For example, as shown in Figure1(a), two vehicles have different trajectories but share the same trip origin and trip destination. Figure 1(b) shows that the passenger inflow time series of subway station A and subway station B are similar in a week, which both have peaks in the morning and the peak values of weekends are relatively small to weekdays.Figure 1
The graph description of similar patterns in the transportation system. (a) The same trip origin and trip destination. (b) The similar passenger flow time series.
(a)(b)These similarities are beneficial for urban mobility pattern understanding and the authorities' policy-making. For example, for aggregate-level, the classification management can be adopted in metro systems and the managers should pay more attention to station A and station B to prevent congestion during the morning peak. For individual-level, the travelers that have similar trip rules can be found, that is, familiar stranger [1] so that the authorities can perform more precise demand control and provide personalized travel service. In addition, the similarities can be used for anomaly detection and improving the traffic prediction as a prior knowledge [2–4].In general, the similarity pattern could be extracted by clustering methods, such as K-means algorithm and the density-based spatial clustering of applications with noise (DBSCAN) algorithm [5]. However, the conventional clustering algorithms always require the complete data, while the missing data problem is inevitable in traffic data collection process. For example, some observations will not be recorded if the sensors were broken. Traditionally, the missing data should be filled firstly; then, a clustering approach is implemented for the imputed data [6]. The missing data imputation and clustering are conducted separately, which is complex and time-consuming. Besides, the frequently used clustering methods such as K-means algorithm must set the number of clusters in advance, and the tedious sensitivity analysis is necessary.Aiming at the missing data imputation, a variety of methods have been proposed such as multioutput Gaussian processes [7], deep generative models [8], and Bayesian tensor decomposition [9], among which Bayesian tensor decomposition is proved to be more effective and efficient than the other methods. Tensor (i.e., multiway arrays) is a generalization of matrix, which is widely used for modeling multidimensional traffic data and can capture the multidimensional structural dependencies; for example, the passenger flow of subway stations can be represented with a three-order tensor (station × day × time of the day). Subsequently, the tensor decomposition can be employed to project station, day, and time into a latent space, also known as factor matrices. Finally, the missing data can be filled by the interaction of factor matrices. Recently, many studies have solved tensor decomposition in the Bayesian formula [9,10], which can provide more robust estimation and uncertainty measures for those missing data.Inspired by the recent work of Bayesian tensor decomposition, we propose a novel framework named Bayesian nonparametric tensor decomposition (BNPTD) to achieve incomplete traffic data imputation and similarity pattern discovery simultaneously. The BNPTD consists of two components: (1) Bayesian tensor decomposition and (2) Dirichlet process mixture model (DPMM). DPMM is a nonparametric Bayesian mixture model that has shown great promise for data clustering and allows for the automatic determination of an appropriate number of clusters. These two components are combined with a hierarchical probabilistic model and a variational inference algorithm is presented to derive the posterior distributions of all the model parameters and hyperparameters. It should be noted that the combination not only finds groups of similar objects in the case of missing data but also offers adaptive prior to Bayesian tensor decomposition model which can further improve better imputation performance.In summary, our contributions are summarized as follows:(i)
We proposed a Bayesian nonparametric tensor decomposition model to achieve incomplete traffic data imputation and similarity pattern discovery simultaneously.(ii)
We presented a variational inference algorithm to learn the BNPTD model. Variational inference tends to be faster and easier to scale to large-scale datasets than classical methods, such as Markov chain Monte Carlo sampling.(iii)
Extensive experiments were conducted on a smart card dataset from Guangzhou, China. The results demonstrate that our approach successfully finds interpretable similarity pattern and recovers the missing values well in the case of incomplete traffic data.The rest of this paper is structured as follows. Section2 presents a literature review. Section 3 describes the materials and methods. The results and discussion are presented in Section 4. Finally, Section 5 summarizes the main conclusions of this work.


## 2. Literature Review

In this section, we review the previous approaches regarding missing data imputation and similarity pattern discovery in the traffic data analysis.


### 2.1. Missing Data Imputation

Missing data problem is inevitable in traffic data collection process, and the reasons are manifold. For example, some observations will be lost if the sensors were broken. Besides, when we use mobile sensors such as floating car, some locations will not be covered at some times. So missing data imputation is always a hot topic in academic research and predecessors have solved this problem using a wide range of methods. The traditional approach is to use time series method such as autoregressive integrated moving average (ARIMA) and its variants [11]. Recently, machine learning methods that considering spatial-temporal correlations have achieved the state-of-the-art performance. For example, Rodrigues et al. [7] proposed the use of multioutput Gaussian processes to model the complex spatial and temporal patterns for crowdsourced traffic data imputation. Yoon et al. [12] adapted the well-known Generative Adversarial Nets (GAN) framework and proposed Generative Adversarial Imputation Nets.The collaborative filtering methods such as tensor decomposition have been proved more effective and efficient than other methods. Tensor (i.e., multiway arrays) is a generalization of matrix. The two most popular tensor decomposition frameworks are Tucker and CANDECOMP/PARAFAC (CP), which both could capture the underlying multilinear factors. Kolda et al. [13] gave a comprehensive overview of tensor decomposition and its applications. Tan et al. [14] introduced the tensor to preserve the multiway nature of traffic data and developed a tensor decomposition based imputation method for missing traffic data completion. Zhao et al. [10] formulated CP decomposition using a hierarchical probabilistic model and incorporated a sparsity-inducing prior over factor matrices, resulting in automatic rank determination. Chen et al. [9] extended the Bayesian probabilistic matrix factorization to higher-order tensor to recover missing values. Zhang et al. [15] proposed an iterative tensor decomposition approach and proved its effectiveness under different missing cases.


### 2.2. Similarity Pattern Discovery

A crucial task in traffic data analysis is similarity pattern discovery and the similarities can help to demand control, personalized travel service, and anomaly detection. The similarity pattern can be extracted by clustering methods. Zhao et al. [2] proposed an unsupervised clustering-based approach to classify passengers in terms of the similarity of their travel patterns and find out the abnormal passengers in metro systems. Gan et al. [16] used the k-means algorithm for clustering metro stations and identified urban mobility patterns from a spatiotemporal perspective.Besides, the nonparametric clustering such as Dirichlet process mixture model (DPMM) is widely used for similarity pattern discovery. DPMM is a nonparametric Bayesian mixture model which allows for the automatic determination of an appropriate number of clusters. The most commonly used description of DPMM is Chinese Restaurant Process that can be solved by MCMC sampling [17]. Blei et al. [18] presented an alternative stick-breaking representation and adopted a variational inference algorithm for DPMM. It is universally acknowledged that variational inference is more efficient than MCMC sampling and adapted to large-scale datasets. In traffic data analysis, Ngan et al. [19] and Ngan et al. [20] employed DPMM to outlier detection and validated its effectiveness via numerical experiments. However, DPMM only can apply to complete data while the real-world traffic data is incomplete inevitably.In this paper, we propose Bayesian nonparametric tensor decomposition that is comprised of Bayesian tensor decomposition and Dirichlet process mixture model via a hierarchical probabilistic model to achieve incomplete traffic data imputation and similarity pattern discovery simultaneously. Compared to the previous studies, DPMM is employed to low-rank feature extracted from Bayesian tensor decomposition rather than raw data, which can solve similarity pattern discovery in the case of incomplete traffic data effectively. Moreover, DPMM is beneficial for missing data imputation as an adaptive prior distribution.


## 2.1. Missing Data Imputation

Missing data problem is inevitable in traffic data collection process, and the reasons are manifold. For example, some observations will be lost if the sensors were broken. Besides, when we use mobile sensors such as floating car, some locations will not be covered at some times. So missing data imputation is always a hot topic in academic research and predecessors have solved this problem using a wide range of methods. The traditional approach is to use time series method such as autoregressive integrated moving average (ARIMA) and its variants [11]. Recently, machine learning methods that considering spatial-temporal correlations have achieved the state-of-the-art performance. For example, Rodrigues et al. [7] proposed the use of multioutput Gaussian processes to model the complex spatial and temporal patterns for crowdsourced traffic data imputation. Yoon et al. [12] adapted the well-known Generative Adversarial Nets (GAN) framework and proposed Generative Adversarial Imputation Nets.The collaborative filtering methods such as tensor decomposition have been proved more effective and efficient than other methods. Tensor (i.e., multiway arrays) is a generalization of matrix. The two most popular tensor decomposition frameworks are Tucker and CANDECOMP/PARAFAC (CP), which both could capture the underlying multilinear factors. Kolda et al. [13] gave a comprehensive overview of tensor decomposition and its applications. Tan et al. [14] introduced the tensor to preserve the multiway nature of traffic data and developed a tensor decomposition based imputation method for missing traffic data completion. Zhao et al. [10] formulated CP decomposition using a hierarchical probabilistic model and incorporated a sparsity-inducing prior over factor matrices, resulting in automatic rank determination. Chen et al. [9] extended the Bayesian probabilistic matrix factorization to higher-order tensor to recover missing values. Zhang et al. [15] proposed an iterative tensor decomposition approach and proved its effectiveness under different missing cases.


## 2.2. Similarity Pattern Discovery

A crucial task in traffic data analysis is similarity pattern discovery and the similarities can help to demand control, personalized travel service, and anomaly detection. The similarity pattern can be extracted by clustering methods. Zhao et al. [2] proposed an unsupervised clustering-based approach to classify passengers in terms of the similarity of their travel patterns and find out the abnormal passengers in metro systems. Gan et al. [16] used the k-means algorithm for clustering metro stations and identified urban mobility patterns from a spatiotemporal perspective.Besides, the nonparametric clustering such as Dirichlet process mixture model (DPMM) is widely used for similarity pattern discovery. DPMM is a nonparametric Bayesian mixture model which allows for the automatic determination of an appropriate number of clusters. The most commonly used description of DPMM is Chinese Restaurant Process that can be solved by MCMC sampling [17]. Blei et al. [18] presented an alternative stick-breaking representation and adopted a variational inference algorithm for DPMM. It is universally acknowledged that variational inference is more efficient than MCMC sampling and adapted to large-scale datasets. In traffic data analysis, Ngan et al. [19] and Ngan et al. [20] employed DPMM to outlier detection and validated its effectiveness via numerical experiments. However, DPMM only can apply to complete data while the real-world traffic data is incomplete inevitably.In this paper, we propose Bayesian nonparametric tensor decomposition that is comprised of Bayesian tensor decomposition and Dirichlet process mixture model via a hierarchical probabilistic model to achieve incomplete traffic data imputation and similarity pattern discovery simultaneously. Compared to the previous studies, DPMM is employed to low-rank feature extracted from Bayesian tensor decomposition rather than raw data, which can solve similarity pattern discovery in the case of incomplete traffic data effectively. Moreover, DPMM is beneficial for missing data imputation as an adaptive prior distribution.


## 3. Materials and Methods

The proposed BNPTD consists of two components: Bayesian tensor decomposition model and Dirichlet process mixture model. In the following subsections, we first give detailed descriptions of these two components. Secondly, we provide an overview of our Bayesian nonparametric tensor decomposition model. Finally, we derive variational approximations to learn the proposed BNPTD. The description below is taken from smart card data as an example and can be applied to other spatiotemporal traffic data.


### 3.1. Bayesian Tensor Decomposition Model

A smart card data record usually has multidimensional attributes such as origin station, destination station, day, and time of day. We organize the smart card data into a three-order tensorX∈ℝI×J×K, where I, J, and K represent the size of station dimension, day dimension, and time dimension, respectively. Each entry in the tensor X is denoted as xijk to represent the inflow (or outflow) of ith station in kth time of jth day.X could be an incomplete tensor when the traffic data is missing. We assume the element xijk is observed if i,j,k∈Ω, where Ω denotes a set of indices. Furthermore, we define a binary tensor O as an indicator of the observed entries, which has the same size as X.Then, CP decomposition factorizes the tensor into a sum of rank-one tensors; that is,(1)X≈U,V,W≡∑r=1Rur∘vr∘wr,where U,V, and W represent the factor matrices of each dimension, respectively, and ur,vr, and wr are the rth column vectors of factor matrices. The symbol ∘ denotes the outer product of the vectors and R is the tensor rank. Figure 2 is a visualization of the CP decomposition of a third-order tensor.Figure 2
CP decomposition of a third-order tensor.Elementwise, (1) is written as(2)xijk≈x^ijk=∑r=1Ruirvjrwkr.CP decomposition can be solved by the alternating least squares method [13]. In order to obtain a more robust solution and avoid overfitting, we resort to Bayesian formula and assume Gaussian observation as follows:(3)xijk∼Nx^ijk,τ−1,where N⋅ denotes the Gaussian distribution and τ is precision.The maximum likelihood estimation of (2) is equivalent to the minimum of square loss function X−X^Ω2. To avoid overfitting, we place conjugate prior distributions over factor matrices and hyperparameters, which correspond to regularization terms:(4)ui:∼Nμu,Λu−1μu,Λu∼Gaussain−Wishartμ0,β0,W0,v0,where ui: is the ith row vector of factor matrix U and V,W have similar prior distributions.We also place a conjugate prior distribution over the precision; that is,(5)τ∼Gammaa0,b0.The graphical model of Bayesian tensor decomposition is illustrated in Figure3.Figure 3
The graphical model of Bayesian tensor decomposition.


### 3.2. Dirichlet Process Mixture Model

The smart card data can be represented with a three-order tensor (station× day × time of the day); then Bayesian tensor decomposition factorizes the tensor into factor matrices U and V,W. The factor matrices also denote the features of each attribute; for example, ui: reflects passenger departure (or arrival) demand patterns of ith station in the specific time interval. In order to find groups of stations that have similar demand pattern, we place a Dirichlet process mixture prior over the factor matrix U.Dirichlet process mixture model is the Bayesian nonparametric technique which allows for the automatic determination of an appropriate number of mixture components. In this paper, we adopt the stick-breaking construction of DPMM [18]. Specifically, let zi be a categorical variable that indicates demand pattern category of ith station. T is truncation level that could be chosen large relative to the number of stations. The generative process of DPMM is given as follows:(1)
For each modet∈1,2,…,T(a)
Drawβt∼Beta1,α(b)
Computeπtβ=βt∏m=1t−11−βm(c)
Drawμt∼G0(2)
For each stationi∈1,2,…,I(a)
Drawzi∼Categoricalπβ(b)
Drawui:∼Nμzi,Iwhere α is hyperparameter that could be set in advance. πtβ is the probability value of mode t. G0 is base distribution, which could be chosen Gaussian distribution N0,λI in this work. I is identity matrix and λ is hyperparameter. The graphical model of DPMM is illustrated in Figure 4.Figure 4
The graphical model of DPMM.


### 3.3. The Proposed Bayesian Nonparametric Tensor Decomposition Model

The BNPTD is comprised of Bayesian tensor decomposition and DPMM, and these components are coupled via a hierarchical probabilistic model. The graphical model of BNPTD is illustrated in Figure5.Figure 5
The graphical model of BNPTD.The advantages of the combination are as follows: (1) missing data imputation and similarity pattern discovery can be achieved simultaneously, which is more efficient and avoids errors accumulation, and (2) DPMM can offer adaptive prior distribution to Bayesian tensor decomposition, corresponding to adaptive regularization term and further better imputation performance.


### 3.4. Model Learning via Variational Inference

The goal of BNPTD model learning is to derive the posterior distributions of parameters and hyperparameters. Variational inference is the method that approximates the posterior distributions through optimization [21]. Specifically, for a group of random variables Θ, variational inference seeks a variational distribution qΘ to approximate the true posterior distribution pΘ|XΩ. The approximation is measured by Kullback–Leibler divergence, and a Kullback–Leibler divergence of 0 indicates that the two distributions are identical:(6)KLqΘ∥pΘ|XΩ=∫qΘlnqΘpΘ|XΩdΘ=ElnqΘ−ElnpΘ|XΩ,where all expectations are taken with respect to qΘ. Because the KL divergence is difficult to compute, variational inference optimizes an alternative objective that the minimization of KL is equivalent to the maximization of the evidence lower bound (i.e., ELBO).(7)KLqΘpΘXΩ=ElnqΘ−ElnpΘ,XΩ+lnpXΩ=ELBOqΘ+lnpXΩ.We further focus on the mean-field variational family that the variational distributionqΘ can be factorized with regard to each variable Θj; that is,(8)qΘ=∏j=1mqjΘj.On the basis of the mean-field variational family, the optimal variational posterior distribution of eachΘj can be derived in turn via maximizing the ELBO, and the optimal form is given by(9)lnqjΘj=EqΘ\ΘjlnpXΩ,Θ+const.where notation EqΘ\Θj⋅. denotes an expectation with respect to distributions qΘ over all variables except for Θj. The deduction process of optimal variational posterior distributions with regard to parameters and hyperparameters can be found in Supplementary File. Variational inference updates the model parameters and hyperparameters alternatively. Algorithm 1 outlines the training process of BNPTD.Algorithm 1: Bayesian nonparametric tensor decomposition (BNPTD).
Input: incomplete tensor XΩ∈ℝI×J×K, indicator tensor O∈ℝI×J×K, CP rank R, maximum clusters T.Output: estimated tensor X^∈ℝI×J×K, the results of similarity pattern discovery. Initialize factor matrices and hyperparameters.(1)repeat(2) fort=1 to T do(3) Update the variational posterior distributionqμt,qβt.(4) Update the variational posterior distributionqμv,Λv and qμw,Λw.(5) fori=1 to I do(6) Update the variational posterior distributionqzi,qui.(7)forj=1 to J do(8) Update the variational posterior distributionqvj.(9) fork=1 to K do(10) Update the variational posterior distributionqwk.(11) Update the variational posterior distributionqτ.(12) until maximum iterations exhausted


## 3.1. Bayesian Tensor Decomposition Model

A smart card data record usually has multidimensional attributes such as origin station, destination station, day, and time of day. We organize the smart card data into a three-order tensorX∈ℝI×J×K, where I, J, and K represent the size of station dimension, day dimension, and time dimension, respectively. Each entry in the tensor X is denoted as xijk to represent the inflow (or outflow) of ith station in kth time of jth day.X could be an incomplete tensor when the traffic data is missing. We assume the element xijk is observed if i,j,k∈Ω, where Ω denotes a set of indices. Furthermore, we define a binary tensor O as an indicator of the observed entries, which has the same size as X.Then, CP decomposition factorizes the tensor into a sum of rank-one tensors; that is,(1)X≈U,V,W≡∑r=1Rur∘vr∘wr,where U,V, and W represent the factor matrices of each dimension, respectively, and ur,vr, and wr are the rth column vectors of factor matrices. The symbol ∘ denotes the outer product of the vectors and R is the tensor rank. Figure 2 is a visualization of the CP decomposition of a third-order tensor.Figure 2
CP decomposition of a third-order tensor.Elementwise, (1) is written as(2)xijk≈x^ijk=∑r=1Ruirvjrwkr.CP decomposition can be solved by the alternating least squares method [13]. In order to obtain a more robust solution and avoid overfitting, we resort to Bayesian formula and assume Gaussian observation as follows:(3)xijk∼Nx^ijk,τ−1,where N⋅ denotes the Gaussian distribution and τ is precision.The maximum likelihood estimation of (2) is equivalent to the minimum of square loss function X−X^Ω2. To avoid overfitting, we place conjugate prior distributions over factor matrices and hyperparameters, which correspond to regularization terms:(4)ui:∼Nμu,Λu−1μu,Λu∼Gaussain−Wishartμ0,β0,W0,v0,where ui: is the ith row vector of factor matrix U and V,W have similar prior distributions.We also place a conjugate prior distribution over the precision; that is,(5)τ∼Gammaa0,b0.The graphical model of Bayesian tensor decomposition is illustrated in Figure3.Figure 3
The graphical model of Bayesian tensor decomposition.


## 3.2. Dirichlet Process Mixture Model

The smart card data can be represented with a three-order tensor (station× day × time of the day); then Bayesian tensor decomposition factorizes the tensor into factor matrices U and V,W. The factor matrices also denote the features of each attribute; for example, ui: reflects passenger departure (or arrival) demand patterns of ith station in the specific time interval. In order to find groups of stations that have similar demand pattern, we place a Dirichlet process mixture prior over the factor matrix U.Dirichlet process mixture model is the Bayesian nonparametric technique which allows for the automatic determination of an appropriate number of mixture components. In this paper, we adopt the stick-breaking construction of DPMM [18]. Specifically, let zi be a categorical variable that indicates demand pattern category of ith station. T is truncation level that could be chosen large relative to the number of stations. The generative process of DPMM is given as follows:(1)
For each modet∈1,2,…,T(a)
Drawβt∼Beta1,α(b)
Computeπtβ=βt∏m=1t−11−βm(c)
Drawμt∼G0(2)
For each stationi∈1,2,…,I(a)
Drawzi∼Categoricalπβ(b)
Drawui:∼Nμzi,Iwhere α is hyperparameter that could be set in advance. πtβ is the probability value of mode t. G0 is base distribution, which could be chosen Gaussian distribution N0,λI in this work. I is identity matrix and λ is hyperparameter. The graphical model of DPMM is illustrated in Figure 4.Figure 4
The graphical model of DPMM.


## 3.3. The Proposed Bayesian Nonparametric Tensor Decomposition Model

The BNPTD is comprised of Bayesian tensor decomposition and DPMM, and these components are coupled via a hierarchical probabilistic model. The graphical model of BNPTD is illustrated in Figure5.Figure 5
The graphical model of BNPTD.The advantages of the combination are as follows: (1) missing data imputation and similarity pattern discovery can be achieved simultaneously, which is more efficient and avoids errors accumulation, and (2) DPMM can offer adaptive prior distribution to Bayesian tensor decomposition, corresponding to adaptive regularization term and further better imputation performance.


## 3.4. Model Learning via Variational Inference

The goal of BNPTD model learning is to derive the posterior distributions of parameters and hyperparameters. Variational inference is the method that approximates the posterior distributions through optimization [21]. Specifically, for a group of random variables Θ, variational inference seeks a variational distribution qΘ to approximate the true posterior distribution pΘ|XΩ. The approximation is measured by Kullback–Leibler divergence, and a Kullback–Leibler divergence of 0 indicates that the two distributions are identical:(6)KLqΘ∥pΘ|XΩ=∫qΘlnqΘpΘ|XΩdΘ=ElnqΘ−ElnpΘ|XΩ,where all expectations are taken with respect to qΘ. Because the KL divergence is difficult to compute, variational inference optimizes an alternative objective that the minimization of KL is equivalent to the maximization of the evidence lower bound (i.e., ELBO).(7)KLqΘpΘXΩ=ElnqΘ−ElnpΘ,XΩ+lnpXΩ=ELBOqΘ+lnpXΩ.We further focus on the mean-field variational family that the variational distributionqΘ can be factorized with regard to each variable Θj; that is,(8)qΘ=∏j=1mqjΘj.On the basis of the mean-field variational family, the optimal variational posterior distribution of eachΘj can be derived in turn via maximizing the ELBO, and the optimal form is given by(9)lnqjΘj=EqΘ\ΘjlnpXΩ,Θ+const.where notation EqΘ\Θj⋅. denotes an expectation with respect to distributions qΘ over all variables except for Θj. The deduction process of optimal variational posterior distributions with regard to parameters and hyperparameters can be found in Supplementary File. Variational inference updates the model parameters and hyperparameters alternatively. Algorithm 1 outlines the training process of BNPTD.Algorithm 1: Bayesian nonparametric tensor decomposition (BNPTD).
Input: incomplete tensor XΩ∈ℝI×J×K, indicator tensor O∈ℝI×J×K, CP rank R, maximum clusters T.Output: estimated tensor X^∈ℝI×J×K, the results of similarity pattern discovery. Initialize factor matrices and hyperparameters.(1)repeat(2) fort=1 to T do(3) Update the variational posterior distributionqμt,qβt.(4) Update the variational posterior distributionqμv,Λv and qμw,Λw.(5) fori=1 to I do(6) Update the variational posterior distributionqzi,qui.(7)forj=1 to J do(8) Update the variational posterior distributionqvj.(9) fork=1 to K do(10) Update the variational posterior distributionqwk.(11) Update the variational posterior distributionqτ.(12) until maximum iterations exhausted


## 4. Results and Discussion

In order to evaluate the performance of our methods, we conduct extensive experiments on a real-world smart card dataset collected from subway stations in Guangzhou, China. This dataset contains a large number of tap-in/tap-out records from 3/7/2017 to 16/7/2017. Because the subway agency does not provide services all day, we focus on the smart card data from 6 a.m. to 10 p.m., which contain the main trip time of the whole day. We construct the dataset into a three-order tensor (station× day × time of the day). Each entry of this tensor denotes the passenger inflow. The length of each time interval is set as 10 minutes, and the size of the tensor is 148×14×96.


### 4.1. The Missing Data Imputation Performance

In order to show the effectiveness of the adaptive prior distribution, we compare BNPTD with some baselines including DA (daily average), BTD (Bayesian tensor decomposition), and GAIN (deep generative models [8]) in terms of missing data imputation. The BTD is also trained via variational inference algorithm and two methods have the same initialization. We implement experiments in the two missing scenarios, which are random missing scenario and continuous missing scenario. Besides, we set T=12 (the maximum clusters), R=10,20,30(CP rank).In our experiments, we evaluate the imputation accuracy with two widely applied metrics, namely, Root Mean Square Error (RMSE) and Mean Absolute Error (MAE):(10)RMSE=1N∑i=1Nyi−y^i2,(11)MAE=1N∑i=1Nyi−y^i.The missing ratio ranged from 10% to 50%. Table1 shows the performance of our proposed method as compared to other baselines in the random missing scenario, and the result of the continuous missing scenario is shown in Table 2. Our proposed BNPTD outperforms other baselines by achieving the lowest RMSE and MAE in most of the missing ratio cases. In some cases, the performance of BNPTD in Tables 1 and 2 is close to BTD, and the reason may be that the adaptive prior distribution is close to single prior distribution in numerical results. But the adaptive prior distribution is conducive to similarity pattern discovery compared to single prior distribution.Table 1
Comparison of BNPTD with other baselines in the random missing scenario (the numbers 10, 20, 30 after the two methods denote CP rank and the best models are highlighted in bold).
Method10%20%30%40%50%RMSEMAERMSEMAERMSEMAERMSEMAERMSEMAEBTD (10)57.3734.0456.5833.6162.3536.2065.2337.6465.5137.83BNPTD (10)56.3533.5156.2533.4558.9734.1162.9436.5465.4137.79GAIN56.9134.7756.29344662.2237.2965.2438.9165.6639.27BTD (20)46.5228.2647.9328.6448.8628.9653.6630.7450.0129.42BNPTD (20)46.6128.2347.6528.5648.1128.6148.9729.0250.1329.43GAIN46.7429.2448.2629.7149.3130.1854.1532.0450.5130.76BTD (30)43.7926.5944.9626.9045.7927.2146.8027.6649.2128.77BNPTD (30)43.4626.544.726.8645.2727.0945.9627.4848.4328.25GAIN44.1927.6645.4628.0246.3928.4447.4128.9249.9230.20DA74.2944.0772.5743.0079.2345.8081.8146.9380.8246.33Table 2
Comparison of BNPTD with other baselines in the continuous missing scenario (the numbers 10, 20, 30 after the two methods denote CP rank and the best models are highlighted in bold).
Method10%20%30%40%50%RMSEMAERMSEMAERMSEMAERMSEMAERMSEMAEBTD (10)60.9236.7359.6735.7665.7938.3869.1039.7669.7140.03BNPTD (10)60.1935.6659.1035.0765.0137.7066.7638.1866.0837.56GAIN59.2936.8858.3336.0864.5938.9468.0240.7668.9741.36BTD (20)49.8530.7051.2930.7452.2230.9757.6232.7453.7731.35BNPTD (20)49.2629.8150.8030.1451.6030.4155.6731.4450.9729.42GAIN48.6931.0150.0431.1151.1831.5256.4633.5653.0632.39BTD (30)49.2930.6050.7030.8752.0831.3553.7831.9458.0534.02BNPTD (30)48.7129.7150.2130.2751.4630.7951.9630.675.0231.92GAIN48.7731.3150.4231.7152.0132.4854.2233.7358.4136.00DA73.3643.8971.8142.6178.3645.1580.9346.1280.6745.95


### 4.2. The Similarity Pattern Discovery

The CP rankR represents the size of low-rank feature, and large CP rank may lead to redundant feature which is bad for similarity pattern discovery. In this subsection, we set CP rank R=8 and the maximum clusters T=12 empirically. The random missing ratio of incomplete smart card data is 10%. The result of similarity pattern discovery is illustrated in Figure 6, and the passenger flow time series that has similar departure pattern is plotted in a subfigure. There are 7 clusters which represent different departure patterns. Moreover, Figure 7 gives a visualization of stations distribution in Google Maps. The detailed description of these patterns is given as follows:(a)
In the case of Cluster 1, the subway stations have large passenger volume every day, and there is no difference between a working day and nonworking day. From Figure7, we observe that these stations are adjacent to passenger transportation hub such as train station and high-speed railway station.(b)
The second departure pattern is anomalous, and the passenger volume of nonworking days from 8/7/2017 to 9/7/2017 is larger than the working day. With our web search, the reason is that there is an exhibition happening near the subway stations. So, the proposed BNPTD also finds some special events effectively and can be used for anomaly detection.(c)
Aimed at Cluster 3, the passenger volume of the stations is lower than 600 in 10 minutes. These stations are not usually crowded and the transportation potential could be further excavated to service more passengers.(d)
Cluster 4 has a larger passenger volume compared to Cluster 3. The stations belonging to Cluster 4 have small peak values during the morning rush hour or evening rush hour. Besides, these stations are mainly situated in urban fringe area and far from the central business area from Figure7.(e)
Cluster 5 is a morning peak departure pattern, and the peak values of weekends are relatively small to weekdays. The stations are mainly in close proximity to the residential area in Figure7. In the morning, people will leave home and go to work by subway.(f)
The stations with the sixth departure pattern have larger passenger volume in the evening rush hour, and the peak values of weekends are relatively small to weekdays. The stations are mainly located in the work district, and many people working in this area usually take the subway and return home after going off work.(g)
The seventh departure pattern has an extremely large passenger volume in the evening rush hour, and these stations are situated in the middle of the central business district (CBD) from Figure7.Figure 6
The result of similarity pattern discovery. (a) Cluster 1. (b) Cluster 2. (c) Cluster 3. (d) Cluster 4. (e) Cluster 5. (f) Cluster 6. (g) Cluster 7.
(a)(b)(c)(d)(e)(f)(g)Figure 7
The subway stations distribution in Google Maps.


### 4.3. The Clustering Stability with Different Missing Ratios

In this subsection, an evaluation measure called purity is adopted to show the clustering stability of BNPTD with different missing ratios. Purity is an external index that is used to measure the similarity of the formed clusters to external clusters such as ground truth [22], considering G=G1,G2,…,Gm as external clusters and C=C1,C2,…,Gn are the clusters made by the BNPTD. To compute purity, each cluster of C is assigned to the class which is most frequent in the cluster, and the accuracy of this assignment is measured by counting the number of correctly assigned stations and dividing by the number of stations I; that is,(12)purityG,C=1I∑k=1nmaxmCk∩​Gm.It is easy to see that0≤purityG,C≤1, and there is a perfect clustering if purityG,C is equal to 1. We set the clustering result in missing ratio 10% as external clusters and implement experiments from missing ratios 15% to 50% with step size 5%. Figure 8 shows the purity with different missing ratios. We observe that purity is greater than 0.6 in all cases, which indicates the clustering result is stability.Figure 8
The purity with different missing ratios.


## 4.1. The Missing Data Imputation Performance

In order to show the effectiveness of the adaptive prior distribution, we compare BNPTD with some baselines including DA (daily average), BTD (Bayesian tensor decomposition), and GAIN (deep generative models [8]) in terms of missing data imputation. The BTD is also trained via variational inference algorithm and two methods have the same initialization. We implement experiments in the two missing scenarios, which are random missing scenario and continuous missing scenario. Besides, we set T=12 (the maximum clusters), R=10,20,30(CP rank).In our experiments, we evaluate the imputation accuracy with two widely applied metrics, namely, Root Mean Square Error (RMSE) and Mean Absolute Error (MAE):(10)RMSE=1N∑i=1Nyi−y^i2,(11)MAE=1N∑i=1Nyi−y^i.The missing ratio ranged from 10% to 50%. Table1 shows the performance of our proposed method as compared to other baselines in the random missing scenario, and the result of the continuous missing scenario is shown in Table 2. Our proposed BNPTD outperforms other baselines by achieving the lowest RMSE and MAE in most of the missing ratio cases. In some cases, the performance of BNPTD in Tables 1 and 2 is close to BTD, and the reason may be that the adaptive prior distribution is close to single prior distribution in numerical results. But the adaptive prior distribution is conducive to similarity pattern discovery compared to single prior distribution.Table 1
Comparison of BNPTD with other baselines in the random missing scenario (the numbers 10, 20, 30 after the two methods denote CP rank and the best models are highlighted in bold).
Method10%20%30%40%50%RMSEMAERMSEMAERMSEMAERMSEMAERMSEMAEBTD (10)57.3734.0456.5833.6162.3536.2065.2337.6465.5137.83BNPTD (10)56.3533.5156.2533.4558.9734.1162.9436.5465.4137.79GAIN56.9134.7756.29344662.2237.2965.2438.9165.6639.27BTD (20)46.5228.2647.9328.6448.8628.9653.6630.7450.0129.42BNPTD (20)46.6128.2347.6528.5648.1128.6148.9729.0250.1329.43GAIN46.7429.2448.2629.7149.3130.1854.1532.0450.5130.76BTD (30)43.7926.5944.9626.9045.7927.2146.8027.6649.2128.77BNPTD (30)43.4626.544.726.8645.2727.0945.9627.4848.4328.25GAIN44.1927.6645.4628.0246.3928.4447.4128.9249.9230.20DA74.2944.0772.5743.0079.2345.8081.8146.9380.8246.33Table 2
Comparison of BNPTD with other baselines in the continuous missing scenario (the numbers 10, 20, 30 after the two methods denote CP rank and the best models are highlighted in bold).
Method10%20%30%40%50%RMSEMAERMSEMAERMSEMAERMSEMAERMSEMAEBTD (10)60.9236.7359.6735.7665.7938.3869.1039.7669.7140.03BNPTD (10)60.1935.6659.1035.0765.0137.7066.7638.1866.0837.56GAIN59.2936.8858.3336.0864.5938.9468.0240.7668.9741.36BTD (20)49.8530.7051.2930.7452.2230.9757.6232.7453.7731.35BNPTD (20)49.2629.8150.8030.1451.6030.4155.6731.4450.9729.42GAIN48.6931.0150.0431.1151.1831.5256.4633.5653.0632.39BTD (30)49.2930.6050.7030.8752.0831.3553.7831.9458.0534.02BNPTD (30)48.7129.7150.2130.2751.4630.7951.9630.675.0231.92GAIN48.7731.3150.4231.7152.0132.4854.2233.7358.4136.00DA73.3643.8971.8142.6178.3645.1580.9346.1280.6745.95


## 4.2. The Similarity Pattern Discovery

The CP rankR represents the size of low-rank feature, and large CP rank may lead to redundant feature which is bad for similarity pattern discovery. In this subsection, we set CP rank R=8 and the maximum clusters T=12 empirically. The random missing ratio of incomplete smart card data is 10%. The result of similarity pattern discovery is illustrated in Figure 6, and the passenger flow time series that has similar departure pattern is plotted in a subfigure. There are 7 clusters which represent different departure patterns. Moreover, Figure 7 gives a visualization of stations distribution in Google Maps. The detailed description of these patterns is given as follows:(a)
In the case of Cluster 1, the subway stations have large passenger volume every day, and there is no difference between a working day and nonworking day. From Figure7, we observe that these stations are adjacent to passenger transportation hub such as train station and high-speed railway station.(b)
The second departure pattern is anomalous, and the passenger volume of nonworking days from 8/7/2017 to 9/7/2017 is larger than the working day. With our web search, the reason is that there is an exhibition happening near the subway stations. So, the proposed BNPTD also finds some special events effectively and can be used for anomaly detection.(c)
Aimed at Cluster 3, the passenger volume of the stations is lower than 600 in 10 minutes. These stations are not usually crowded and the transportation potential could be further excavated to service more passengers.(d)
Cluster 4 has a larger passenger volume compared to Cluster 3. The stations belonging to Cluster 4 have small peak values during the morning rush hour or evening rush hour. Besides, these stations are mainly situated in urban fringe area and far from the central business area from Figure7.(e)
Cluster 5 is a morning peak departure pattern, and the peak values of weekends are relatively small to weekdays. The stations are mainly in close proximity to the residential area in Figure7. In the morning, people will leave home and go to work by subway.(f)
The stations with the sixth departure pattern have larger passenger volume in the evening rush hour, and the peak values of weekends are relatively small to weekdays. The stations are mainly located in the work district, and many people working in this area usually take the subway and return home after going off work.(g)
The seventh departure pattern has an extremely large passenger volume in the evening rush hour, and these stations are situated in the middle of the central business district (CBD) from Figure7.Figure 6
The result of similarity pattern discovery. (a) Cluster 1. (b) Cluster 2. (c) Cluster 3. (d) Cluster 4. (e) Cluster 5. (f) Cluster 6. (g) Cluster 7.
(a)(b)(c)(d)(e)(f)(g)Figure 7
The subway stations distribution in Google Maps.


## 4.3. The Clustering Stability with Different Missing Ratios

In this subsection, an evaluation measure called purity is adopted to show the clustering stability of BNPTD with different missing ratios. Purity is an external index that is used to measure the similarity of the formed clusters to external clusters such as ground truth [22], considering G=G1,G2,…,Gm as external clusters and C=C1,C2,…,Gn are the clusters made by the BNPTD. To compute purity, each cluster of C is assigned to the class which is most frequent in the cluster, and the accuracy of this assignment is measured by counting the number of correctly assigned stations and dividing by the number of stations I; that is,(12)purityG,C=1I∑k=1nmaxmCk∩​Gm.It is easy to see that0≤purityG,C≤1, and there is a perfect clustering if purityG,C is equal to 1. We set the clustering result in missing ratio 10% as external clusters and implement experiments from missing ratios 15% to 50% with step size 5%. Figure 8 shows the purity with different missing ratios. We observe that purity is greater than 0.6 in all cases, which indicates the clustering result is stability.Figure 8
The purity with different missing ratios.


## 5. Conclusions

In this paper, a novel Bayesian nonparametric tensor decomposition model called BNPTD is proposed. The model combines the Bayesian tensor decomposition and Dirichlet process mixture model via a hierarchical probabilistic model, which can achieve incomplete traffic data imputation and similarity pattern discovery simultaneously. Moreover, we derive a variational inference algorithm to learn the model efficiently. Experiments on a real-world smart card dataset show the effectiveness of the proposed model.Actually, the similarity pattern extracted in this paper not only helps us understand the urban mobility patterns, but also can improve the traffic prediction as prior knowledge. For future work, we plan to develop multitask prediction based on the above similarities.



---

*Source: 8810753-2020-08-12.xml*