Individualized passenger travel pattern multi-clustering based on graph regularized tensor latent dirichlet allocation

Individual passenger travel patterns have significant value in understanding passenger’s behavior, such as learning the hidden clusters of locations, time, and passengers. The learned clusters further enable commercially beneficial actions such as customized services, promotions, data-driven urban-use planning, peak hour discovery, and so on. However, the individualized passenger modeling is very challenging for the following reasons: 1) The individual passenger travel data are multi-dimensional spatiotemporal big data, including at least the origin, destination, and time dimensions; 2) Moreover, individualized passenger travel patterns usually depend on the external environment, such as the distances and functions of locations, which are ignored in most current works. This work proposes a multi-clustering model to learn the latent clusters along the multiple dimensions of Origin, Destination, Time, and eventually, Passenger (ODT-P). We develop a graph-regularized tensor Latent Dirichlet Allocation (LDA) model by first extending the traditional LDA model into a tensor version and then applies to individual travel data. Then, the external information of stations is formulated as semantic graphs and incorporated as the Laplacian regularizations; Furthermore, to improve the model scalability when dealing with massive data, an online stochastic learning method based on tensorized variational Expectation-Maximization algorithm is developed. Finally, a case study based on passengers in the Hong Kong metro system is conducted and demonstrates that a better clustering performance is achieved compared to state-of-the-arts with the improvement in point-wise mutual information index and algorithm convergence speed by a factor of two.


Introduction
The public transportation system is the backbone of a city's infrastructure, and the intelligent transportation system (ITS) has been an essential chapter for the smart city blueprint. Most studies for ITS focus on traffic flow prediction (Ren and Xie 2017;Geng et al. 2019;Guo et al. 2019;Shi et al. 2020;Wang et al. 2019;Li et al. 2020). Tensor-based methods, such as tensor decomposition (Ren and Xie 2017), tensor completion (Li et al. 2020), as well as deep learning methods, such as convolutional neural networks (Geng et al. 2019), graph convolutional networks (Yu et al. 2018), and spatiotemporal attentions (Guo et al. 2019), have been developed to predict city-wide traffic flow (Geng et al. 2019), metro station-level passenger flow (Li et al. 2020), origin-destination (OD) flow matrix (Ren and Xie 2017;Wang et al. 2019;Shi et al. 2020), etc. However, the methods mentioned above target traffic flow prediction at a macro level and utilize the traffic data of passengers indiscriminately, consequently neglecting the personalized travel characteristics of individual passengers. For example, to calculate the passenger flow value, only the number of passengers is counted, which abandons individual information (Yi et al. 2019). Thus, those methods could not directly handle individual travel data.
To overcome this issue, we propose to fully utilize the "individual" travel data for an "individualized" travel pattern discovery. Individual travel data preserve the abundant trajectory information, i.e., that passenger u departs from origin o at time t and arrives at destination d at time t . This encourages us to focus on the following individualized analysis tasks due to their high research values (Zhao et al. 2017). We aim at the following two goals: • Clustering of origin, destination, and time (ODT): The latent clusters for origin, destination, and time could be better learned from individual travel pattern data, given the abundant information is well preserved. The intuition is that those origins may belong to the same cluster if they all co-occur in the same type of passengers, as shown in Fig. 1(a). The learned clusters for ODT could guide better urban planning, more suitable station-surrounding facilities, and uncover the peak hour for crowd control. • Clustering of passengers: Passengers will also be clustered into different groups based on their trajectories. For public transport providers, with a better understanding of individual passengers' travel patterns, customized promotions and more suitable operational policies can be designed. For example, the fare surchargereward scheme could be tailored for different passenger groups (Tang et al. 2020).
The two clustering tasks above for Origin, Destination, Time, and Passenger are called for short as "ODT-P Multi-clustering". However, these two tasks are rather challenging due to the multi-mode spatiotemporal big data and the influence from the external environment.
• Challenge 1: Multi-mode spatiotemporal big data. Take Hong Kong as an example, there are 2 million passengers daily: Each passenger has multiple trips, and each trip has multiple modes such as origin, destination, and time. • Challenge 2: external environment. Moreover, passenger behaviors are also affected by the external environment, such as the locations and surroundings of stations. If two stations are geographically adjacent to each other or located in similar functional areas (such as business area, residential area, school), they will attract the same type of passengers.
To tackle the aforementioned challenges, we propose a novel Graph-Regularized Tensor Latent Dirichlet Allocation model (GR-TensorLDA) for ODT-P multiclustering based on individual passenger travel patterns. First, to preserve the multi-mode structure of the high dimensional spatiotemporal data, we focus on the tensor-based methodology (Kolda and Bader 2009), which represents the original data with three-mode tensors, where different modes represent ODT respectively. Secondly, a tensor LDA model is proposed to achieve ODT-P multi-clustering.
The main novelty of our proposed method is that we extend the traditional LDA (Blei et al. 2003) to a tensor version and apply into individual traffic data. An important analogy is made as shown in Fig. 1(b): (1) "Word"-level: A trip is viewed as a three-dimensional word w = (w O , w D , w T ); (2) "Document"-level: A passenger with several trips, i.e., "a bag of words", is treated as a three-dimensional document d u ∈ R O×D×T . Generative processes in the passenger-level and trip-level will be defined along with each mode of ODT; (3) "Topic"-level: Therefore, the latent topic will also be formulated as a tensor, with each element as z = (z O , z D , z T ).
The clusters of ODT-P will be eventually obtained in the following way: (1) ODT clustering: Along each dimension of ODT, the topic is a latent distribution of words, which can be viewed as a cluster containing different words; (2) P clustering: each passenger is represented by the latent distribution of the tensor topics, which will be utilized to cluster passengers.
Our most significant technical contributions are twofold: • Semantic graph structure: To tackle Challenge 2, we incorporate the external environment as graph structures into the model. Precisely, we first formulate the station-related information as two graphs: (1) A geographical graph measures the spatial distance between stations; (2) A contextual graph quantifies whether two stations are located in similar functional areas (Geng et al. 2019;Li et al. 2020).
Then the graph structures are incorporated into the tensor LDA generative process for OD, such that if two stations are close on these two graphs, they are more likely to be in the same topic. We show that by adding such graph regularizations, the interpretability of the learned ODT-P clusters can be significantly improved. • Efficient online algorithm: Since the graph regularization breaks the conjugacy, standard optimization techniques such as Gibbs sampling (Griffiths and Steyvers 2004) are no longer possible, we propose a tensorized variational Expectation-Maximization (EM) algorithm to estimate parameters. Moreover, to tackle Challenge 1, we need an efficient and scalable algorithm to deal with massive passenger data. Therefore, we further propose to conduct the algorithm in an online stochastic learning manner (Hoffman et al. 2010). We show that to reach the same level of performance, the online learning algorithm converges twice faster than the batch learning algorithm.
The remainder of the paper is structured as follows. Section 2 briefly reviews existing tensor methods, individual travel analysis, and topic models. Section 3 formulates the proposed model, and Section 4 proposes an efficient optimization algorithm. Section 5 provides a detailed experiment to demonstrate the improved meaningfulness of the learned clusters and model scalability; Section 6 gives the conclusion and discusses the future work and model generalization.

Tensor and tucker decomposition
We would like to first introduce tensor and tensor decomposition since high dimensional data are usually formulated as a tensor, and tensor decomposition is widely used for clustering (Kolda and Bader 2009;Sun and Axhausen 2016). Tensor is mathematically defined as a multi-dimensional array, which is believed to have sufficient capacity to preserve innate complex correlations across multiple dimensions (Kolda and Bader 2009). One of the most popular techniques is tensor decomposition. Tucker decomposition is a high-order principal component analysis. It decomposes a tensor X O×D×T into a core tensor C J ×K ×L multiplied by a mode matrix along each dimension, U O , U D , U T : X = C × 1 U O × 2 U D × 3 U T : Tensor decomposition has been applied into smart transportation for prediction (Li et al. 2020) and clustering (Sun and Axhausen 2016). However, tensor decomposition might be rather impractical to be applied to our problem. The main reason is: data formulated with each individual passenger are extremely sparse due to curse of dimensionality. According to our preliminary study (Li et al. 2021), the sparsity could reach 99.97%, which paralyzes traditional tensor decomposition methods (Tang et al. 2018). Therefore, a technique that specializes in individual analysis is needed.

Individual travel pattern analysis
The new generation of ITS aims to be more personalized. Recently individual travel data have been utilized for passenger clustering (Briand et al. 2017;Mohamed et al. 2016), station clustering (Mohamed et al. 2016), and personalized services such as travel time estimation (Tang et al. 2018), route recommendation , destination inference (Cheng et al. 2020), driving state recognition (Yi et al. 2019), and activity discovery (Zhao et al. 2020). There are mainly two kinds of approaches as follows.

Spatio-temporal feature engineering
The first kind of approach relies on intense feature engineering to extract useful features such as spatial, temporal, OD pair, transportation mode. It then combines the features with traditional statistical learning models for clustering and prediction, such as Kmean clustering (Zhao et al. 2017), boosting tree , random forest (Yi et al. 2019). However, the feature extraction is rather complicated and differs from one system to another, which does not offer a universal solution. Furthermore, it typically assumes that the feature has the same dimension for each passenger. However, the number of trips of each passenger can be dramatically different. In contrast, our model learns the latent dimension in a data-driven way and can be accommodated to different numbers of trips, which offers a general solution with explainable results.

Generative models
As the second option, generative models (Briand et al. 2017;Mohamed et al. 2016;Tang et al. 2018;Cheng et al. 2020;Zhao et al. 2020) have been adopted into individual traffic analysis. To cluster passengers' temporal behaviors, Briand et al. (2017) and Mohamed et al. (2016) proposed two-layer generative models with a mixture of Gaussian or unigrams model: The first layer partitions passengers into clusters, and the second layer formulates each cluster's temporal activity. However, the limitation is that passengers are only clustered based on their active or boarding time; Therefore, the passengers' abundant spatiotemporal information is not fully utilized. As a result, no insights about the latent nature of origins and destinations could be obtained.
To capture all dimensions of the spatiotemporal information for individual passengers, researchers have adopted topic models into individual traffic data (Cheng et al. 2020;Zhao et al. 2020), where an individual passenger's travel data is regarded as a document, where each trip is recorded as a word. Specifically, Cheng et al. (2020) and Zhao et al. (2020) proposed a high-dimensional LDA model with a generative process on each dimension, such as location, day, hour, and trip duration. However, the existing methods ignored the underlining spatial correlations in the passengers' travel data, which may lead to a clustering model that could not reflect reality. Compared with them, our most significant advantages are that we incorporate semantic graphs into the LDA generative process. This is inspired by the state-of-the-art topic models (Yao et al. 2017;Li et al. 2019b) that incorporates knowledge graph to improve the interpretability of model output, which we will review them in details in the following section. Such coupling of the "pure-trip" data with external contexts significantly improves topics' interpretability, and that we propose an efficient online stochastic learning algorithm based on a variational EM algorithm.

Graph-based topic models
Incorporating knowledge graphs into topic models could enhance the interpretability of the learned topics (Yao et al. 2017;Li et al. 2019a;Mei et al. 2008;Chen et al. 2016;Li et al. 2019b). In particular, two categories of incorporating methods are considered in state-of-the-art models. The first category embeds words into a continuous space with word relations defined by an external knowledge graph such as DBpedia, WordNet (Yao et al. 2017;Li et al. 2019a). However, in traffic data, such knowledge graph is only applied to a single word representation, which cannot be used for high-dimensional word representation. The second category introduces graph-based regularization (Mei et al. 2008;Chen et al. 2016;Li et al. 2019b), such as graph harmonic function, to encourage entities close on the graph to be more likely to have the same topic. These regularization-based techniques are compatible with our generative model. However, the existing graph-based topic models are formulated only for one-dimensional word, not for high-dimensional data like our passenger travel data. Moreover, the challenges lie in the parameter learning for our corresponding tensor topic model. To this end, we rigorously develop online stochastic learning based on tensorized variational EM algorithm to estimate parameters with higher efficiency and scalability.

Multi-view subspace clustering
Last but not the least, subspace clustering is also a popular method for highdimensional clustering (Parsons et al. 2004), which learns data representation in certain low-dimensional subspaces and clusters the data points. Multi-view subspace clustering (Gao et al. 2015;Zhang et al. 2017Zhang et al. , 2018 specifically deals with data represented by multiple distinct feature sets. We would like to emphasize the difference between our model and subspace clustering from the perspectives of data and model: (1) The typical multi-view data are formulated as X v ∈ R d v ×n , where d v and n are the number of features and samples on the v-th view. Our data instead present a hierarchical structure: a passenger has a sequence of a few trip records, and each trip is an instance in a three-dimensional space of ODT. Moreover, our data suffer from high sparsity, which also hinders subspace clustering, e.g., factorization-based methods, from normal functioning. (2) A typical formulation of multi-view subspace clustering is based on the data's selfexpression property (Gao et al. 2015), which is to use the data set to represent itself: is the subspace representation matrix of the v-th view, and the nonzero elements in Z v correspond to the data points from the same subspace. Various methods are proposed to add different regularizations on Z v , such as sparsity (Elhamifar and Vidal 2013), low rank (Liu et al. 2012) and smoothness (Hu et al. 2014). However, self-expression property cannot be applied to our model since our data are already high-dimensional and extremely sparse, which may lead to a even higher-dimensional and more sparse Z.

Proposed methodology
We will introduce the proposed methodology. Section 3.1 gives the data formulation and the notations; Section 3.2 introduces the tensor topic and tucker decomposition; Section 3.3 formulates the generative process along each dimension; Section 3.4 formulates the graph structure for dimension OD; Section 3.5 gives the final loss function. The concepts of "passenger" and "document", "trip" and "word", "topic" and "cluster" are interchangeably used here.

Data representation and notation
Firstly, the notations throughout this paper are as follows: We denote scalars in italics, e.g., n, vectors by lowercase letters in boldface, e.g., β, matrices by uppercase boldface letters, e.g., B, and tensors by boldface script capital W.
Then, we would like to define the data representation. A trip is defined as a threedimensional tuple (i.e., word) w = (w O , w D , w T ), indicating a trip that starts from origin w O at time w T and heads to destination w D . V O , V D , V T are the vocabulary sizes for ODT respectively. A passenger who has traveled several trips is regarded as "a bag of words" (i.e., document): d u = {w 1 , . . . , w i , . . . , w N u }, with i as the i-trip of the passenger u, N u as the number of trips from passenger u ∈ R M , and M as the total number of passengers. All the notations in the paper are summarized in Table 1.

Tensor topic definition and decomposition
The topic for the i-th word w i in passenger u is also formulated as a three-element indicates the i-th word belongs to the j-th 'O' topic, k-th 'D' topic and l-th 'T' topic, respectively (Cheng et al. 2020).
The vocabulary of all unique ODT Variational variable for z O , the probability of the word at i-th position generated from j-th topic.
According to Bayes' theorem, the probability of the i-trip from passenger d u can be written as: (1) (1) can be presented in Tucker decomposition as follows: The essence of the model is revealed as probabilistic tucker decomposition (Kolda and Bader 2009), where the tensor ODT data for each passenger u is W u ∈ R V O ×V D ×V T , which is decomposed into a core tensor C u ∈ R J ×K ×L , and along It is worth mentioning that although the essence of the model is a decomposition, yet since W u is intractable, C u , B O , B D , B T could not be learned by decomposing W u . Instead, we learn the latent parameters first and then W u could be calculated.

Generative process
The generative process for a trip w will be defined along ODT.
Prior: Dirichlet distribution is known as a good conjugate prior for multinomial distribution (Zhou 2018). The tensor topic distribution C u ∈ R J ×K ×L for the uth passenger is generated from the 3D-Dirichlet distribution with parameter A ∈ R J ×K ×L and each element c u, j,k,l ( j k l c u, j,k,l = 1) defines the possibility for the passenger to have trips from topic z j,k,l : (3) Passenger to Tensor Topic: The topic for the i-trip in the u-th passenger is drawn from the multinomial distribution: (4)

Graph structure on origin and destination
If two stations are geographically close to each other or located in the similar functional area, intuitively these two stations are more likely to be in the same topic. This external information will be formulated as a graph and then introduced into the model as the Laplacian regularization. Precisely, two graphs are defined accordingly to capture the inter-relationships of different stations: (1) Geographical graph G net : describes how two stations are geographically close to each other on the network; (2) Functional similarity graph G poi : quantifies how similar the functions of two stations' locations are. These two graphs have effects on both OD dimensions, with definition details in Section 5.2.
The graph regularization term is defined as follows. Take one graph in origin dimension as an example, where the weight between entity o 1 and o 2 , and L is the Laplacian matrix for graph G O (Li et al. 2020;Wang et al. 2015;Yu et al. 2019). The intuition is that two stations that are closer on the graph will be more likely to have the same topic (Mei et al. 2008).
With the corresponding Laplacian matrices as L net and L poi , the final graph Laplacian penalty on OD could be formulated as: where the tuning parameters μ and ν adjust the relative effect of two graphs on OD respectively. The whole generative process is shown in Fig. 3.

Loss function
In order to learn the model parameters A, B O , B D and B T for the GR-TensorLDA model, the log-likelihood function could be formulated as follows: where is the marginal distribution of a passenger which can be defined as: Generative Process for trips from each passenger via latent topics z cannot be computed tractably due to the coupling between A and B O , B D , B T in the summation over latent topics Blei et al. (2003). Luckily, variational inference provides us with a tractable lower bound on the log likelihood, which will be elaborated in detail in Section 4.
After combining the external knowledge, the model parameters are learned by maximizing the regularized likelihood function, with λ tuning the penalty strength:

Parameter estimation
In this section, we describe a tensorized variational EM-algorithm to optimize the model parameters A, B O , B D , B T in Eq. (9), efficiently. Based on the existing variational EM-algorithm (Blei et al. 2003), we mainly emphasize the following significant contributions in our algorithm: (1) For E-step in Section 4.1, the variational E-step is extended from one-dimension words to high-dimension words with tucker decomposition; (2) For M-step in Section 4.2, the gradient ascend method is adopted to address the graph regularizations; (3) Most importantly, in Section 4.3, an online learning algorithm is proposed to handle the big data problem in smart transportation systems. The algorithm is summarized as follows: • Tensorized variational E-step: to approximate posterior, four variational distributions q(·)s are introduced for C, z O , z D , z T with free variational parameters E, O , D , T respectively, as shown in Fig. 4; then the lower bound (L B) for the original log-likelihood is calculated by Jensen's inequality, with optimal variational parameters learned to maximize the L B;

Tensorized variational E-Step
The variational distribution is formulated as Eq. (10) to approximate the posterior distribution of each passenger: where φ O i ∈ R J , with φ O i j interpreting the probability that word at i-th position in current document is generated from origin topic j.
A tight lower-bound is found by minimizing Kullback-Leibler (KL) divergence between the inference distribution Q and posterior P: As shown in Appendix A.2, the optimal variational parameters E * , ( O * , D * , T * ) are learned by computing the derivatives of the KL divergence and setting them to zero, with results shown in Eqs. (12) and (13).
To estimate φ O i j for the u-th passenger by Appendix A.2.1: The parameter of one dimension, for example, φ O i j , is not only related to its own dimension but also other dimensions φ D ik and φ T il . Therefore, we could perform the block coordinate descent algorithm, which iteratively update the parameters for ODT dimensions until convergence.
To estimate j,k,l for the u-th passenger via Appendix A.2.2:  u, j,k,l = α u, j,k,l + N u J ×K ×L for all u, i, j, k, l. 4: for iter = 1 to max_iter I n f er do 5: for u = 1 to M do 6:

M-Step
In the M-step, we aim to maximize the lower bound learned from E-step with respect to B O , B D , B T and A.
As shown in Appendix A.3.1, B O and B D cannot be solved in the closed-form solution due to the graph regularization. Therefore, we propose to use the gradient ascend method to update B O and B D : the gradient with respect to Finally, similar to the original LDA model (Blei et al. 2003), A can be estimated using the Newton-Raphson method: Its detailed derivation is given in Appendix A.3.3. It is worth mentioning that Eqs. (14), (15), and (16) do not show a direct relation between J , K , L and model parameters: this is because J , K , L affect the variational Algorithm 2 Online GR-TensorLDA: Tensorized Variational EM Algorithm Input: passengers, J , K , L, λ, μ, ν, stopping tolerance tol, max_iter I n f er, max I ter β , r , κ, τ 0 . 1: for s = 0 to ∞ do 2: E-Step: 3: if parameter i+1 − parameter i < tol, break 10: end for 11: M-Step: 12: The whole algorithm is shown as Algorithm 1: it is in a batch learning manner which needs to read through the whole document set for each iteration.

Online learning algorithm
In practice, the proposed Algorithm 1 is very computationally intense since it updates parameters in a batch learning manner, which iterates between analyzing each observation and updating dataset-wide variational parameters. Therefore, as shown in Line 5 in Algorithm 1, E-step needs a full pass through the entire corpus each iteration, which is impractical when dealing with a large dataset containing tens of thousands of passengers. For example, Hong Kong, as an international transport hub, the monthly visitor arrivals were recorded to 6 million in Dec-2018, and the batch learning algorithm is not suitable for the situation where new visitors are continually arriving (Hoffman et al. 2010).
To this end, to efficiently implement the proposed model in real traffic systems, we further develop an online stochastic algorithm, which outputs good estimates outstandingly faster than the batch algorithm. To avoid repetitious details, only the critical differences in the algorithm will be explained.
In E-step, the updating equations from Eq. (12) to Eq. (13) remain the same except that variational variables are updated each time a passenger s is read, as shown in Line 6-7 in Algorithm 2.
In M-step, to update model parameters stochastically, once observing the current passenger s, we first assume the optimal model parameters are learned if the entire corpus contained this passenger s repeated M times: under this setting (denoted with the accent symbolβ ),B O andB D are updated by stochastic gradient ascend, with the gradient calculated based on the single observation d s times M: where Similar with Eq. (15),B T has closed-form solution with passenger s repeated M times:β Then the final model parameters B O,D,T are updated by using a weighted average of its previous value andB O,D,T : where ρ s = (τ 0 + s) −κ . κ ∈ (0.5, 1] controls the rate at which old values ofB are forgotten and guarantees convergence; τ 0 ≤ 0 slows down the early iterations. The proposed online learning algorithm of the GR-TensorLDA method is summarized in Algorithm 2. The computational complexity of each iteration in the batch learning algorithm is O (M N (J + K + L)) (M ≫ N , J , K , L), whereas the complexity of each iteration in the online learning algorithm reduces to O(N (J + K + L)) (Teh et al. 2007).
Mini-Batches: To reduce noise, parameters could be updated with a mini-batch containing multiple observations, with mini-batch size as S.B T is updated asβ T .B O,D are updated by stochastic gradient ascend with mini-batches.

Dataset
The individual travel data from 1st-Jan-2017 to 31st-Mar-2017 are chosen for analysis. Each trip has recorded the anonymized passenger ID, entry station, exit station, entry time, and exit time. In this implementation, entry station, exit station, and hour stamp of entry time have been collected for each trip and aggregated over the whole three months to ensure each passenger has enough trips for analysis, with average amount of trips around 134. The Hong Kong MTR system has 98 stations in total and operates in 24 hours. Thus the vocabulary size for origin, destination, and time is 98, 98, 24. We randomly pick 50,000 passengers as the training set, 1000 passengers as the validation set for tuning parameter selection, and another 1000 passengers as the test set. The data information is summarized in Table 2.

Graph definition
As discussed in Section 3.4, the geographical graph and the functional similarity graph affect passengers' travel patterns. Here we would like to define the two graphs in detail. Geographical graph: Two spatially close stations are more likely to be in the same topic. The distance from station i to j in the graph {G net } i, j is usually simplified as an "H -hop" binary graph: if from station i to j less than H -hops are needed, two stations are connected, and the edge between them is '1'; Otherwise, the edge is '0'. We set H = 3 since a survey stated that passengers are willing to travel freely if two stations are only three hops away (Geng et al. 2019;Li et al. 2021).
Functional similarity graph: Two stations located in highly similar functional areas are also prone to be in the same topic. A functional similarity graph is commonly formulated with the point of information (POI) (Li et al. 2020;Zhong et al. 2017). We collect each station's surrounding POIs 1 with the following seven services: hotel, leisure shopping, major building, public facilities, residential, school, and public transport. Each element of the POI vector indicates the amount of the corresponding service nearby the station. The element {G poi } i, j can be defined as cosine similarity between POI vectors of station i and station j, and a higher value in G poi means a higher functional similarity between two stations:

Benchmark methods
We apply the following benchmark methods to passenger travel data and compare the results with the proposed model. However, a relatively limited amount of research targets ODT-P multi-clustering based on individual passenger travel data.
• One-dimensional LDA (1d-LDA): It defines the generative process from document to topic, topic to word in a single dimension (Blei et al. 2003). To apply it, three-dimensional data (w O , w D , w T ) = (o, d, t) are flattened into onedimensional w = odt: each different element creates a new word; thus, the total vocabulary size for the new data format is expanded to 98 × 98 × 24 ∼ 10 5 , and the computational complexity of each iteration is O (M N K ). • Tucker Decomposition (Tucker): It decomposes an ODT flow tensor into a core tensor and mode matrices along each dimension. Each rank vector from a mode matrix can be regarded as a cluster (Sun and Axhausen 2016). However, this method is not applicable to individual travel data due to extreme sparsity. Merely to examine its ODT clustering performance, we feed it with macro-level passenger flow data with the dimension of origin, destination, and time (Sun and Axhausen 2016). • CP Decomposition with Graphs (CP-G): Similar to tucker decomposition, the input is passenger flow data to check its ODT clustering performance. It decomposes an ODT flow tensor into a weight vector and mode matrices along each dimension, with graphs on OD (Li et al. 2020). • Three-dimensional LDA with Gibbs sampling (3d-LDA(Gibbs)): It also defines a generative process for each dimension, however, without any semantic graph structure. Parameters are estimated by Gibbs sampling (Cheng et al. 2020), with the computational complexity of each iteration as O (M N (J + K + L)) (Porteous et al. 2008;Xiao and Stibor 2010).
All the methods are compared by whether it is an individualized analysis (Indiv. for short), tensor-based (Tensor), and graph-structured (Graph) model with efficient online algorithm (Eff.) and low computational complexity (Complexity) as shown in Table 3. Only our model ticks all the boxes.

Evaluation metrics
Traditional topic models only measure the quality of the model via perplexity, but ignore how "interpretable" the learned topics are. For example, a topic containing all words related to covid (e.g., omicron, vaccine, quarantine, etc: these words are all connected in a knowledge graph) is more "interpretable" than a topic containing words from various themes (e.g., covid, solar energy, iPhone: these words are far away to each other in a knowledge graph). Such "interpretability" is usually measured by point-wise mutual information (PMI), known as topic coherence. Besides, given our graph structure on origin and destination dimensions, we also innovatively design distance of graph to measure the "interpretability": a topic with words that are close to each other on a graph is more interpretable than the one that does not. Topic Coherence PMI: PMI is to evaluate how meaningful the learned topics along each dimension are (Yao et al. 2017;Newman et al. 2010). For example, topic j in dimension of the origin stations is calculated as P Distance on Graph (d G ): Based on the Laplacian matrix of a graph, d G measures the distance of the word components from a topic. A smaller value means a more concentrated topic. For example, the distance on graph for origin topic j is defined as (Blei et al. 2003) examines the likelihood of the proposed model in the test set. A lower perplexity means a higher likelihood.

Parameter tuning
The number of topics in each ODT dimension is set as J = 10, K = 10, L = 4 in our dataset. This is chosen by the expert knowledge: J , K = 10 since the POI of each station has seven elements; L = 4 since there are usually at least three timecomponents capturing morning peak, evening peak, and midday trend. Generally, if there is no prior information about the parameters, J , K , L could be determined by a grid search to minimize perplexity (Blei et al. 2003) or maximize topic coherence (Yao et al. 2017). Theoretically, with bigger J , K , L, the dimensions of model parameters A, B O , B D , B T also increase, which means more latent clusters are introduced to describe data pattern: this will naturally increase model's likelihood and decrease perplexity. However, too large J , K , L may cause overfitting. Tuning parameters for graph regularization λ, μ, ν are searched from grids λ, μ, ν ∈ {0.1, 0.2, . . . , 0.9}, and configuration parameters for online algorithm κ ∈ {0.5, 0.6, . . . , 1.0}, τ 0 ∈ {1, 4, 16, 64, 256, 1024} and S ∈ {1, 4, 16, 64, 256, 1024}. The optimal values are chosen to maximize likelihood in the validation set, with λ = 0.5, μ = 0.4, ν = 0.2 and S = 128, κ = 0.5, τ 0 = 256.

Topic matching
It is worth mentioning that before the comparison, the topics learned from different methods need to be matched first. We match two topics from two methods if they have the highest cosine similarity. For example, topic i from method m 1 refers to topicĵ from method m 2 whenĵ = arg max j {cos-similarity(β m 1 i , β m 2 j )}.

PMI and distance on graph
The PMI and d G for the learned origin topics are shown in Table 4, with the best performance highlighted in boldface and the second-best performance highlighted in underline.
From Table 4, we find that: (1) Tucker and CP decomposition have the worst PMI since they are methods targeting macro-level traffic analysis. Thus it ignores the individual passenger information; However, tensor decomposition with graphs considered (i.e., CP-G) still has higher PMI and lower d G than that without graphs; (2) Our proposed method achieves twice higher PMI than 3d-LDA(Gibbs) for most topics, which means the proposed model can discover more meaningful topics. This is due to the external information introduced as graphs, with more than 50% lower d G observed on both graphs.

Perplexity vs interpretability
In Table 5, our model has a 10% higher perplexity in the test set, which means a lower likelihood score for these passengers. When we introduce the regularization  Table 4 term into the loss function, it naturally decreases the likelihood score because regularization terms generally reduce the model fitting accuracy (i.e., likelihood score) in the exchange of a better generalizability on the testing samples. However, in the literature, it has been shown that perplexity is not a good measure compared to the topic coherence PMI score. This is because the perplexity itself does not reflect the meaningfulness of topics, and topics with lower perplexity might even conflict with the real-world knowledge (Yao et al. 2017;Chang et al. 2009). Our experiments show that, by adding the graph regularization, although the model likelihood score (i.e., 10% higher perplexity measure) is lower, the model interpretability and generalization power (i.e., as seen in the twice higher PMI and 50% lower d G measures) are significantly improved. Therefore, we observed such a trade-off between perplexity and interpretability.
1d-LDA has the highest perplexity since it is not a tensor-based model, thus cannot preserve the innate spatiotemporal correlations.

Improved interpretability in station topics
To better demonstrate the enhanced interpretability of the proposed model's topics and check how they reflect the real world, we also visualize: (1) the topic POI features, (2) topic locations, and (3) topic station components based on the real metro map, in comparison to the second-best model 3d-LDA (Gibbs) in origin dimension only. The discovered topics can be used as station clusters.
(1) Topic POI Feature to check the POI feature of each topic: Usually, more distinguishable topics are better (Zhu et al. 2012). The POI features of topics from 3d-LDA(Gibbs) and our model are calculated as j( poi) ∈ R N poi indicates the POI feature of this topic. As shown in Fig. 5. (a2), topics from the proposed model capture more distinct POI groups: origin topic 4, 5, 6, and 8 capture the POI leisure shopping, major building, residential and school respectively; Topics from 3d-LDA(Gibbs) instead, as shown in Fig. 5. (a1), capture the topics with similar and non-distinguishable POI distribution. The distinct POI patterns are observed in our destination topics too.
(2) Topic Location on Map to check each topic's location on map: The top 10 stations with the highest weights from origin topics 4, 5, 6, and 8 are located on the metro map. In Fig. 5. (b), the stations from our topics are concentrated in the same line/region; however, topics without graphs are dispersed among different regions. Therefore, the topics learned from our proposed method have significant improvement in interpretability and reflect external knowledge.
(3) Station Components Analysis to check each topic's station component in terms of the station location and POI: A further detailed study about the top 10 station components with the highest weights inside each topic (here topic 6 is chosen) is conducted to check those stations' exact locations and the surrounding POIs. (1) In Fig. 6. (b), the top 10 stations of our topic 6 are all located in the same metro line, and the POI feature of each station is also mainly residential; (2) On the contrary, in Fig. 6. (a) the top 10 stations of topic 6 without graphs are scattered over different three lines, and those stations also have quite different POI features, such as station 80 in leisure shopping, station 99 in the major building.
To conclude, the identified topics from the proposed model have better physical meanings and interpretability. Applications of OD Clustering: Metro companies usually categorize stations by their expert knowledge, such that different station categories have different operational, marketing, and urban planning strategies. However, this categorization is usually outof-date since a station's feature evolves over time. The learned station clustering is purely data-driven and updated with data, which provides better insights.

Time topics
The time topics from our method are shown in Fig. 7(a), Topic 0 captures the first morning peak (7-8 AM); Topic 1 captures the second morning peak (9-10 AM); Topic 2 captures the mid-day trend, and Topic 3 captures the evening peak (8 PM).
Applications of Time Clustering: The learned time topics could offer clear insights about the peak hours and enable crowd management.

Passenger clustering
The passenger cluster could be learned from each passenger's tensor topic distribution parameter C u . The distance between two passengers' topic distributions could be measured by the Euclidean distance, the Jensen Shannon (JS) divergence  (Cheng et al. 2020) and so on. We choose the JS divergence since it is a symmetric measure for distributions. Then clustering methods such as K-means could be applied to cluster passengers.
Applications of passenger clustering: (1) Customized Services: Passenger clustering could help public transport companies better understand passenger demographics, enabling customized travel reward plans or tailored advertisements for different passenger clusters. (2) Destination Inference: Moreover, the conditional probability based on the latent parameters B O , B D , B T , and C u enables more potential applications. For example, given the passenger u , origin o , and entry time t , the destination could be predicted by: As a result, a destination crowd warning message in the mobile application could then be directed towards each passenger. With B O , B D better estimated with graphs, for passengers who travel between two stations (accounts for 83.8% of the population in our dataset), the destination inference accuracy is improved by 18% compared with 3d-LDA(Gibbs), as shown in Table 6. In Fig. 8, the most popular OD pairs (i.e., o → d ) are selected to visualize destination inference. With the input of passenger u , who has  9 Convergence comparison of (a) log-likelihood evolution and (b) origin perplexity evolution between Online Algorithm 2 in green and Batch Algorithm 1 in red. Each point marker '·' in online version denotes 10 iterations, and each cross marker '×' in batch version denotes 1 iteration these OD pairs, and the time t when he/she enters o , the destination is inferred by our method with higher accuracy (correct d in green); However, the method without graphs makes more mistakes (wrong d in red).

Faster convergence from online algorithm
Last but not least, as shown in Fig. 9, we compare the convergence speed of the algorithm's batch version and its online version in the same computation environment. The proposed online stochastic algorithm needs more iterations to converge but 60% less time: The online version (shown in color green) converges at t ≈ 300, more than twice faster than the batch algorithm (shown in color red), which converges at t ≈ 700. Moreover, the online algorithm also converges with with better parameter estimates: (1) Higher log-likelihood: As shown in Fig. 9. (a), online version converges at loglikelihood ≈ −113 × 10 3 , higher than batch version's convergence at log-likelihood ≈ −118 × 10 3 ; (2) Lower perplexity: as shown in Fig. 9. (b), online version converges at origin perplexity ≈ 60, lower than batch version's convergence at origin perplexity ≈ 89.

Conclusion
In this paper, we studied the ODT-P multi-clustering problem for the individual passenger travel pattern to achieve meaningful clusters on each dimension (origin, destination, and time) and on the individual passenger by incorporating external information on the origin and destination stations. To solve this challenge, we proposed a novel graphregularized tensor Latent Dirichlet Allocation model, which applies to the travel data of each passenger with the consideration of the external information as the graph regularization. We proposed a tensorized variational EM-algorithm to estimate parameters. To improve the scalability, an online learning algorithm is further proposed. In the case study based on the Hong Kong metro system, we demonstrate our superiority over state-of-the-art methods in terms of two times higher topic coherence, 50% lower distance on graph, and better interpretability. Our improvement is also reflected on its 20% more accurate individual destination inference. The proposed online learning method can also converge twice faster with the same good performance as the batch learning method.

Future work
Our model will be extended to cover trip duration, which is a continuous variable. The generative process will be extended to handle continuous distribution correspondingly. Besides, due to the independent assumption of Dirichlet distribution, correlations between passengers will also be further examined.

Generalization
This work could also be applied to non-metro data such as bus and sharing rides if the ODT information is recorded. In the road traffic, the nodes in G net and G poi could be defined as different grids, road segments, or zip code zones, and the edges could be defined similarly if distance and POI are available. The five terms are further expanded using Appendix A.1 result:

A.2.1 Variational multinomial
We keep the terms in lower bound containing φ O i j for example: To maximize it with respect to φ O i j , derivative is calculated and set to be zero: Then we have: Same steps are followed to get φ D ik and φ T il .

A.2.2 Variational dirichlet
The term containing j,k,l is simplified as follows: Similarly, to maximize the lower bound, derivative is calculated: with respect to β O jo is as follows: As mentioned in Section 4.3, B O will be updated in an online stochastic manner, where gradient is calculated with only one observation s repeated M times, thus, B D will be updated same way.

A.3.2 Conditional multinomials B T
There is no regularization in T dimension, so the estimation is similar with the work in (Blei et al. 2003).
A.3.3 Dirichlet A α j,k,l will also be updated with Newton-Raphson method. The gradient g j,k,l with respect to α j,k,l is as below: where c = jkl g j,k,l /h j,k,l z −1 + jkl h j,k,l , h j,k,l = −M (α j,k,l ), and z = M ( jkl α j,k,l ) (Blei et al. 2003).