Low-rank spectral estimation algorithm of learning Markov model

Abstract

This paper proposes a low-rank spectral estimation algorithm of learning Markov model. First, an approximate projection algorithm for the rank-constrained frequency matrix set is proposed, and thereafter its local Lipschitzian error bound established. Then, we propose a low-rank spectral estimation algorithm for estimating the state transition frequency matrix and the probability matrix of Markov model by applying the approximate projection algorithm to correct the maximum likelihood estimation of the frequency matrix, and prove that there is only a multiplying constant difference in estimation errors between the low-rank spectral estimation and the maximum likelihood estimation under appropriate conditions. Finally, numerical comparisons with the prevailing maximum likelihood estimation, spectral estimation, and rank-constrained maximum likelihood estimation show that the low-rank spectral estimation algorithm is effective.

Graphical abstract

Keywords

Markov model / low-rank spectral estimation / error bound / approximate projection

1 Introduction

The learning of Markov model often arose from fields such as statistics, systems engineering, data science, and management [1-3, 14, 16]. For example, the potential urban transportation pattern could be revealed to further facilitate taxi resources allocation, by modeling the urban taxi movement trajectory into a Markov process, and conducting state clustering by estimating its state transition frequency matrix and probability matrix [3, 9]. Similarly, learning corresponding Markov model could also assist in solving the problem of designing the operation route of urban buses [8], personalized pagerank [2], and data sorting in e-commerce [10]. It should be pointed out Markov process originated from the above problems normally possessed a huge state space, but their state transition frequency matrix and probability matrix had been proven to be low-rank or approximately low-rank [3, 9, 16]. Although Markov process had a huge state space, but when the state transition frequency matrix satisfied low rankness, it was possible to accurately learn the process based on a short empirical trajectory [12, 16, 18]. For examples and applications of other low-rank Markov processes, readers, if interested, can refer to references [3, 5, 12, 15, 16] and others.

In recent years, many scholars have conducted research on the estimation of the state transition frequency matrix and probability matrix of low-rank Markov model [5, 16, 18]. As far as we know, existing estimation algorithms with polynomial time complex bound cannot guarantee that the estimation of frequency matrix and probability matrix satisfies the low-rank condition. For example, Zhang and Wang[16] proposed a spectral estimation method for learning low-rank Markov model by combining low-rank truncated singular value decomposition in maximum likelihood estimation of frequency matrix with non-negative projection, established the statistical error bound for the estimation, and proved that the difference between the statistical error bound and the lower bound of minimax error was one logarithmic factor of the Markov chain trajectory length. However, as mentioned in [16], due to the non-negative projection required by the spectral estimation algorithm, the final frequency and probability estimation matrices might not necessarily be low-rank. We have also verified this through numerical experiments, and found out that the spectral estimation matrices are hardly low-rank, but instead full-rank in some cases. This is also the main reason why the error of spectral estimation fails to be near optimal. To overcome this defect, [18] proposed the nuclear norm regularized maximum likelihood estimation and rank-constrained maximum likelihood estimation methods for state transition matrices, established the statistical error bounds for them, and proved that the statistical error bound and the lower bound of the minimax error only differed by one logarithmic factor of the Markov chain state space dimension, which was theoretically near optimal. It is noted that the optimal solution for the nuclear norm regularized problems may not necessarily satisfy the low-rank condition; the optimal solution of rank-constrained optimization problems can otherwise, but it is generally NP-hard to seek the optimal solution in computation. Although [18] designed a class of DC (difference of convex functions) programming algorithm to approximate the solution of rank-constrained maximum likelihood estimation problems. Anyway, nuclear norm regularized maximum likelihood estimation and rank-constrained maximum likelihood estimation require a large amount of computation, because singular value decomposition is needed in each step of the algorithm. Therefore, the estimation method raised in [18] is not suitable for large-scale Markov model learning. In addition, [5] designed a low-rank recovery algorithm for state transition matrices, but the estimation matrix of this method might not necessarily be non-negative, so it couldn’t be always a low-rank state transition matrix. Inspired by the above researches, this paper attempts to seek a method that can quickly obtain low-rank frequency matrix and probability matrix to learn a low-rank Markov model.

Another area related to this paper lies in the study of set error bounds. Error bound plays a core role in the theory and algorithm research of optimization, and the research on it always remained a focus and challenge in the field of optimization [4, 7, 11, 13, 17]. Reference [11] proved that the convex polyhedron set had a global Lipschitz error bound; for general convex feasible systems, Lipschitz error bounds only existed under certain constraint qualifications. Specifically, [11] pointed out that for general non-convex feasible systems, the global error bound, even of Hölder type, could hardly be established. Existing results mainly established local Hölder error bounds for subanalytic systems and polynomial systems [7,11]. As a special class of non-convex feasible system, the research on the error bounds of rank-constrained sets has been limited. For a feasibility system constituting the intersection of a matrix closed convex set and rank-constrained set, when its multivalued function satisfied the condition of calmness at the origin, [4] established the local Lipschitz error bound for estimation of distance from any matrix in the closed convex matrix to the rank-constrained feasibility system. However, the difficulty in verifying the calmness of multivalued function is almost the same with establishing error bound. Therefore, no solution has been found for the error bound of general matrix rank-constrained feasibility system.

The main work of this paper is to propose a low-rank spectral estimation algorithm for learning Markov model based on the approximate projection algorithm of the rank-constrained frequency matrix set. Given a positive integer r\in [1,d-1], this paper defines the rank-constrained frequency matrix set as follows:

\mathcal{F}:=\{Z\in {\mathbb{R}}^{d\times d}:\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(Z)\le r,e^{\mathrm{T}}Ze=1,Z\ge 0\},

where rank(Z) stands for the rank of Z, e\in {\mathbb{R}}^d represents the d-dimensional column vector with all elements as 1, matrix Z⩾0 indicates that its elements Z_{ij}⩾0, 1⩽i,j⩽d. Obviously, \frac{1}{d^2}e{e}^{\mathrm{T}}\in \mathcal{F}, so set \mathcal{F} is non-empty and bounded.

First, this paper will establish the local Lipschitz error bound of set \mathcal{F}, that is to give estimation of distance from any frequency matrix F\in \{Z\in {\mathbb{R}}^{d\times d}:e^{\mathrm{T}}Ze=1,Ze\ge \alpha ,Z\ge 0\} to set \mathcal{F}, where \alpha \in (0,\frac{1}{d}) is constant, Ze⩾\alpha indicates (Ze{)}_i⩾\alpha ,i=1,2,\dots ,d. To be specific, a low-rank spectral estimation algorithm for the state transition frequency matrix and probability matrix of the Markov model is proposed. The state transition frequency matrix and probability matrix obtained by our algorithm can satisfy the low-rank condition, thus overcoming the limitation that the spectral estimation in [16] is not a low-rank state transition matrix. Besides, our method does not require solving the optimization problem, so it has a low computational complexity compared to the rank-constrained maximum likelihood estimation [18]. And we further prove that the estimation error between low-rank spectral estimation and maximum likelihood estimation only differs in one multiplying constant under the condition that the probability of each state is equal to or greater than α. Finally, the proposed method is numerically compared with the maximum likelihood estimation [1], spectral estimation [16], and rank-constrained maximum likelihood estimation [18], showing that our method is effective; in addition, we combine the low-rank spectral estimation algorithm and k-means clustering algorithm into the analysis of taxi trajectory in Manhattan Island, New York City, USA.

At the end of this section, we present the main notations of this paper.

(1) Denote E\in {\mathbb{R}}^{d\times d} as a matrix with all elements as 1; e\in {\mathbb{R}}^d as a column vector with all elements as 1.

(2) For any matrix Z\in {\mathbb{R}}^{d\times d}, define \parallel Z{\parallel }_{\mathrm{\infty }}:={\mathrm{m}\mathrm{a}\mathrm{x}}_{1\le i,j\le d} \left|Z_{ij}\right| as an infinite norm, \parallel Z{\parallel }_1:=\sum _{i,j=1}^d \left|Z_{ij}\right| is l₁-norm, \parallel Z{\parallel }_F:=\sqrt{\sum _{i,j=1}^d Z_{ij}^2} is F(robenius) norm, lth row vector of Z is written as \begin{array}{c}Z(l,:)\in {\mathbb{R}}^d\end{array}, given a\in \mathbb{R}, Z⩾a represents Z_{ij}\ge a, 1⩽i,j⩽d. Define the distance from Z\in {\mathbb{R}}^{d\times d} to l₁-norm of \mathcal{F} as

\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(Z,\mathcal{F}):=\mathrm{m}\mathrm{i}\mathrm{n}\{\parallel Z-X{\parallel }_1:X\in \mathcal{F}\}.

(3) Given x\in {\mathbb{R}}^d, define l₁-norm as \parallel x{\parallel }_1:=\sum _{i=1}^d \left|x_i\right|, given a\in \mathbb{R}, x⩾a represents x_i⩾a,1⩽i⩽d.

(4) For any matrix Z\in {\mathbb{R}}^{d\times d}, define {\mathrm{\Pi }}_r\left(Z\right):=\sum _{i=1}^r {\sigma }_i\left(Z\right)u_i\left(Z\right)v_i(Z{)}^{\mathrm{T}}, where {\sigma }_i\left(Z\right) is the ith maximum singular value of Z, \begin{array}{c}{u}_i\left(Z\right)\in {\mathbb{R}}^d\end{array} and \begin{array}{c}{v}_i\left(Z\right)\in {\mathbb{R}}^d\end{array} corresponds to left and right singular vector of {\sigma }_i\left(Z\right), where i=1,2,\dots ,r. As is known to all, {\mathrm{\Pi }}_r\left(Z\right)\in \mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\{\parallel Z-X{\parallel }_F:\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(X)\le r\}.

(5) Define {\mathrm{\Omega }}_{\alpha }:=\{Z\in {\mathbb{R}}^{d\times d}:e^{\mathrm{T}}Ze=1,Ze\ge \alpha ,Z\ge 0\}, where \alpha \in (0,\frac{1}{d}) is constant.

2 Algorithm and its properties

In this section, we will propose a low-rank spectral estimation algorithm for the state transition frequency matrix and probability matrix of a Markov model. We consider a discrete-time Markov process with d(⩾2) states \{S_1,S_2,\dots ,S_d\}, and assume that its state transition probability matrix is \overline{P}\in {\mathbb{R}}^{d\times d} and frequency matrix is \overline{F}\in {\mathbb{R}}^{d\times d}, with satisfied low-rank condition \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\overline{P})=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\overline{F})=r\ll d. It is aimed to estimate the matrices \overline{P} and \overline{F} through the Markov chain \left\{X_0,X_1,\dots ,X_n\right\} with a trajectory length of n + 1. Firstly, we present the assumptions in this section.

Assumption 2.1　 \left\{X_0,X_1,\dots ,X_n\right\} is an ergodic Markov chain and there exists a constant \alpha \in (0,\frac{1}{d}), such that \overline{F}e\ge \alpha.

Assumption 2.1 is also required for the theoretical analysis in [16]. Let’s mark the stationary distribution of Markov chain as \pi \in {\mathbb{R}}^d. It follows from [16] that \overline{F}=\mathrm{D}\mathrm{i}\mathrm{a}\mathrm{g}\left(\pi \right)\overline{P}, where \mathrm{D}\mathrm{i}\mathrm{a}\mathrm{g}\left(\pi \right)\in {\mathbb{R}}^{d\times d} indicates the diagonal matrix of the ith diagonal element π_i(i=1,2,\dots ,d). Since {\overline{P}}_e=e, \pi =\overline{F}e, so Assumption 1.1 is equivalent to \pi ⩾\alpha. When \begin{array}{c}{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d} {\pi }_i\end{array} is relatively small, the probability of some states appearing in X_0,X_1,\dots ,X_n is extremely low, posing a challenge to the estimation of state transition frequency matrix and probability matrix.

2.1 Approximate projection algorithm of rank-constrained frequency matrix set

As is widely known, for any matrix F\in {\mathbb{R}}^{d\times d}, it was NP-hard to compute its projection in \mathcal{F} [6]. In the following, we will propose a fast approximate projection algorithm for any F\in {\mathrm{\Omega }}_{\alpha } in \mathcal{F}, and establish its local Lipschitz error bound.

Algorithm 2.1 (approximate projection algorithm of F\in {\mathrm{\Omega }}_{\alpha } in \mathcal{F})

Input: F\in {\mathrm{\Omega }}_a, r\in [1,d-1], \kappa \in (0,0.5), \mu \ge \eta >0, iteration step T.

(S1) Compute {\mathrm{\Pi }}_r\left(F\right). If , output \hat{F}=\frac{1}{d^2}, end; or turn to (S2).

(S2) If \begin{array}{c}{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d} \left({\mathrm{\Pi }}_r\right(F){)}_{ij}\ge 0\end{array}, output \hat{F}=\frac{{\mathrm{\Pi }}_r\left(F\right)}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}, end; or turn to (S3).

(S3) Set F_0={\mathrm{\Pi }}_r\left(F\right), k = 1, w\in {\mathbb{R}}^d, where w_i=\mathrm{m}\mathrm{i}\mathrm{n}(\mu ,\mathrm{m}\mathrm{a}\mathrm{x}((F^{\mathrm{T}}e{)}_i,\eta )), t_j=0,1\le j\le T.

Compute

while

　　if k\le T-1

　　 [I,J]=\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}(F_0={\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d} (F_0{)}_{ij});t_k={\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d,j\in J}\frac{F_0(i,j)}{\left({\mathrm{\Pi }}_r\right(F)e{)}_i};

　　 F_0=F_0-t_k{\mathrm{\Pi }}_r\left(F\right)e\frac{w^{\mathrm{T}}}{\underset{j\in J}{min} w_j};

　　else

　　 t_k={\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d}\frac{F_0(i,j)}{\left({\mathrm{\Pi }}_r\right(F)e{)}_i};F_0=F_0-t_k{\mathrm{\Pi }}_r\left(F\right)e{e}^{\mathrm{T}};

　　end if

　　 k = k + 1;

end while

Output \hat{F}=\frac{F_0}{\parallel F_0{\parallel }_1}, end.

Lemma 2.1　 For any F\in {\mathrm{\Omega }}_{\alpha }, output matrix \hat{F} in Algorithm 2.1 satisfies \hat{F}\in \mathcal{F} and \hat{F}e>0.

Proof　If , then \hat{F}=\frac{1}{d^2}E\in \mathcal{F} and \hat{F}e=\frac{1}{d}, so the conclusion holds. If {\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d}\left({\mathrm{\Pi }}_r\right(F\left)e\right)\ge \kappa \alpha且{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d}\left({\mathrm{\Pi }}_r\right(F){)}_{ij}\ge 0, then \hat{F}=\frac{{\mathrm{\Pi }}_r\left(F\right)}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}\in \mathcal{F} and \hat{F}e=\frac{{\mathrm{\Pi }}_r\left(F\right)e}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}\ge \frac{\kappa \alpha }{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}>0, so the conclusion also holds.

Next, let’s assume {\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d}\left({\mathrm{\Pi }}_r\right(F\left)e\right)\ge \kappa \alpha and . It follows from (S3) in Algorithm 2.1 that F_0⩾0 is the output in the end and can be expressed as F_0={\mathrm{\Pi }}_r\left(F\right)-{\mathrm{\Pi }}_r\left(F\right)e{q}^{\mathrm{T}}, where q\in {\mathbb{R}}^d and q⩽0. Since {\mathrm{\Pi }}_r\left(F\right)e>0, \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(F_0\right)\le r, and \parallel F_0{\parallel }_1\ge \parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1\ge \kappa \alpha. Hence, \hat{F}=\frac{F_0}{\parallel F_0{\parallel }_1}\in \mathcal{F}. It is noted that all t_k⩽0 in (S3), so we can get \hat{F}e=\frac{F_{0}e}{\parallel F_0{\parallel }_1}\ge \frac{\kappa \alpha }{\parallel F_0{\parallel }_1}>0.

The following proposition establishes the local Lipschitz error bound of rank-constrained frequency matrix set, namely estimation of distance from any matrix F\in {\mathrm{\Omega }}_{\alpha } to l₁-norm of set \mathcal{F}.

Proposition 2.1　 For any F\in {\mathrm{\Omega }}_{\alpha }, let’s assume the output of Algorithm 2.1 is \hat{F}. The error bound inequation below is obtained:

(2.1)\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F})\le \parallel F-\hat{F}{\parallel }_1\le \frac{2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}}{\alpha \kappa }.

Proof　Given any F\in {\mathrm{\Omega }}_{\alpha }. It follows from Lemma 2.1 that \hat{F}\in \mathcal{F}. \mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F}) is l₁-norm distance, so \mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F})\le \parallel F-\hat{F}{\parallel }_1\le \parallel F{\parallel }_1+\parallel \hat{F}{\parallel }_1=2. We will try to demonstrate from three cases as below.

Case 1　 . Since F\in {\mathrm{\Omega }}_{\alpha }, Fe⩾\alpha. So

\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1\ge \parallel Fe-{\mathrm{\Pi }}_r\left(F\right)e{\parallel }_{\mathrm{\infty }}\ge \alpha -\underset{1\le i\le d}{\mathrm{m}\mathrm{i}\mathrm{n}} \left({\mathrm{\Pi }}_r\right(F\left)e\right)\ge \alpha -\kappa \alpha .

Thus \mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F})\le 2\le \frac{2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1}{\alpha (1-\kappa )}, and inequation (2.1) holds.

Case 2　 {\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d}\left({\mathrm{\Pi }}_r\right(F\left)e\right)\ge \kappa \alpha and 且{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d} \left({\mathrm{\Pi }}_r\right(F){)}_{ij}\ge 0. Then \hat{F}=\frac{{\mathrm{\Pi }}_r\left(F\right)}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}. Based on \parallel F{\parallel }_1=e^{\mathrm{T}}Fe=1, we can get

\begin{array}{rl}\parallel F-\hat{F}{\parallel }_1& ={\Vert F-{\mathrm{\Pi }}_r\left(F\right)+{\mathrm{\Pi }}_r\left(F\right)-\frac{{\mathrm{\Pi }}_r\left(F\right)}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}\Vert }_1\\ & \le {\Vert F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1+\parallel {\mathrm{\Pi }}_r\left(F\right)-\frac{{\mathrm{\Pi }}_r\left(F\right)}{\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1}\Vert }_1\\ & =\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1+|\parallel {\mathrm{\Pi }}_r(F){\parallel }_1-\parallel F{\parallel }_1|\\ & \le 2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1.\end{array}

Thus, \mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F})\le \parallel F-\hat{F}{\parallel }_1\le 2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1\le \frac{2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1}{\alpha (1-\kappa )}, and inequation (1.1) holds.

Case 3　 {\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d}\left({\mathrm{\Pi }}_r\right(F\left)e\right)\ge \kappa \alpha and . If \parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1>2, then

\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,\mathcal{F})\le 2\le 2(\parallel {\mathrm{\Pi }}_r(F){\parallel }_1-\parallel F{\parallel }_1)\le 2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1\le \frac{2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1}{\alpha (1-\kappa )}.

Therefore, inequation (2.1) is established. Then we assume \parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1\le 2. It follows through Algorithm 2.1 (S3) that t_k⩽0, k=1,2,\dots ,T, and combining \begin{array}{c}{\mathrm{\Pi }}_r\left(F\right)e\frac{w^{\mathrm{T}}}{{\mathrm{m}\mathrm{i}\mathrm{n}}_{{j\in }^J} w_j}\ge 0\end{array}, it can be seen that element of F₀ increases with k. Hence, \left\{t_k\right\} is a non-descending sequence while \left\{\left|t_k\right|\right\} a non-increasing sequence, and

\left|t_k\right|\le \left|t_1\right|\le \left|\underset{1\le i,j\le d}{\mathrm{m}\mathrm{i}\mathrm{n}} \frac{{\mathrm{\Pi }}_r\left(F\right)(i,j)}{\left({\mathrm{\Pi }}_r\right(F)e{)}_i}\right|\le \frac{|{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i,j\le d} {\mathrm{\Pi }}_r(F\left)\right(i,j\left)\right|}{\kappa \alpha },k=1,2,\dots ,T.

Note that F⩾0, , so F⩾0, . Thus

(2.2)\parallel F_0-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1\le \frac{d\mu }{\eta }\parallel {\mathrm{\Pi }}_r\left(F\right){\parallel }_1\sum _{i=1}^T \left|t_i\right|\le \frac{2Td\mu }{\eta }\frac{\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}}{\kappa \alpha },

where F₀ is the final output of Algorithm 2.1 (S3). By using \hat{F}=\frac{F_0}{\parallel F_0{\parallel }_1}, \parallel F{\parallel }_1=e^{\mathrm{T}}Fe=1, and inequation (2.2), it follows

\begin{array}{rl}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(F,F)& \le \parallel F-F_0+F_0-\hat{F}{\parallel }_1\\ & \le \parallel F-F_0{\parallel }_1+{\Vert F_0-\frac{F_0}{\parallel F_0{\parallel }_1}\Vert }_1\\ & =\parallel F-F_0{\parallel }_1+|\parallel F_0{\parallel }_1-1|=\parallel F-F_0{\parallel }_1+|\parallel F_0{\parallel }_1-\parallel F{\parallel }_1|\\ & \le 2\parallel F-F_0{\parallel }_1\le 2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1+2\parallel F_0-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1\\ & \le \frac{2\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_1}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}}{\alpha \kappa }.\end{array}

Therefore, inequation (2.1) is established.

Note 2.1　(i) Although the right side of the error bound inequation in Proposition 2.1 contains dimensional factor d, namely d\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}, d\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}=O(\parallel F-{\mathrm{\Pi }}_r(F\left){\parallel }_1\right) when F-{\mathrm{\Pi }}_r\left(F\right) is a dense matrix with relatively uniform elements. In the future, we will conduct researches on d\parallel F-{\mathrm{\Pi }}_r\left(F\right){\parallel }_{\mathrm{\infty }}=O(\parallel F-{\mathrm{\Pi }}_r(F\left){\parallel }_1\right) derived from frequency matrix estimation of low-rank Markov model.

(ii) It is noted that when μ = n, T = 1, the right side of the error bound inequation in Proposition 2.1 is minimal. This is because we scaled \frac{w}{{\mathrm{m}\mathrm{i}\mathrm{n}}_{{j\in }^J} w_j} up to \frac{\mu }{\eta }e in inequation (2). In fact, when \mu \ge F^{\mathrm{T}}e\ge \eta, w=F^{\mathrm{T}}e, so w is used to depict the weights of each column in matrix F, and it can uplift the quality of approximate projection \hat{F} if applied in correcting F₀.

2.2 Low-rank spectral estimation algorithm

Define n_{ij}:=\left|\right\{k:X_{k-1}=S_i,X_k=S_j,k=1,2,\dots ,n\left\}\right|, 1⩽i, j⩽d, namely the sum of transition times from state S_i to S_jis n_ij. Define matrices \tilde{F}\in {\mathbb{R}}^{d\times d} and \begin{array}{c}\tilde{P}\in {\mathbb{R}}^{d\times d}\end{array} as follows:

\tilde{F}(i,j)=\frac{n_{ij}}{n},\tilde{P}(i,j)=\{\begin{array}{l}\frac{n_{ij}}{{{\displaystyle \sum }}_{j=1}^d{n}_{ij}},\mathrm{i}\mathrm{f}{{\displaystyle \sum }}_{j=1}^d{n}_{ij}>0;\\ {\displaystyle \frac{1}{d}},\mathrm{i}\mathrm{f}{{\displaystyle \sum }}_{j=1}^d{n}_{ij}=0,\end{array}1\le i,j\le d.

[1] has been proved that \tilde{F} and \tilde{P} are the maximum likelihood estimation for \overline{F} and \overline{P}, respectively, and featured by strong consistency. Yet, the maximum likelihood estimation hasn’t utilized the low-rank information in \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\overline{P})=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\overline{F})=r\ll d. Considering the low-rankness of matrices \overline{F} and \overline{P} to be estimated, we present the low-rank spectral estimation algorithm of Markov process.

Algorithm 2.2 (Low-rank spectral estimation algorithm)

Input: \tilde{F}, \alpha \in (0,\frac{1}{d}), r\in [1,d-1], \kappa \in (0,0.5), \mu ⩾\eta >0, integer T.

(S1) Compute the corrected maximum likelihood estimation F_1\in {\mathbb{R}}^{d\times d}:

m_{+}=\sum _{i=1}^d \mathrm{m}\mathrm{a}\mathrm{x}\left(\right(\tilde{F}e{)}_i-\alpha ,0);m_{-}=\sum _{i=1}^d \mathrm{m}\mathrm{a}\mathrm{x}(\alpha -(\tilde{F}e{)}_i,0);

　　　　　For i = 1 to d

　　　　　if (\tilde{F}e{)}_i=0

　　　　　 F_1(i,:)=\frac{\alpha }{d}e;;

　　　　　elseif

　　　　　 F_1(i,:)=\tilde{F}(i,:)+(\alpha -(\tilde{F}e{)}_i)\frac{\tilde{F}\left(i,:\right)}{(\tilde{F}e{)}_i};

　　　　　else

　　　　　 F_1(i,:)=\tilde{F}(i,:)-\frac{m_{-}}{m_{+}}\left(\right(\tilde{F}e{)}_i-\alpha )\frac{\tilde{F}(i,:)}{(\tilde{F}e{)}_i};

　　　　　end if

　　　　　end for

(S2) Use Algorithm 2.1 to compute the estimation matrix {\hat{F}}_1 of \overline{F} with F₁, r, κ, μ, η, T as the input.

(S3) Computer the estimation matrix {\hat{P}}_1\in {\mathbb{R}}^{d\times d} of \overline{P}, where {\hat{P}}_1(i,:)=\frac{{\hat{F}}_1(i,:)}{\parallel {\hat{F}}_1(i,:){\parallel }_1},i=1,2,\dots ,d.

The following Theorem 2.1(a) shows that Algorithm 2.2 (S1) results in F_1\in {\mathrm{\Omega }}_{\alpha }. It follows through Lemma 2.1 that {\hat{F}}_1\in \mathcal{F} is a low-rank frequency matrix and {\hat{F}}_{1}e>0, thereby the {\hat{P}}_1 generated by Algorithm 2.2 (S3) is meaningful as a low-rank state transition frequency matrix, namely {\hat{P}}_1\in \{P\in {\mathbb{R}}^{d\times d}:\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(P)\le r,Pe=e,P\ge 0\}. Therefore, the estimation matrices {\hat{F}}_1 and {\hat{F}}_1 in Algorithm 2.2 meet the low-rank condition, overcoming the defect of spectral estimation [16]. In particular, Algorithm 2.2 doesn’t need to solve matrix optimization, so there is little computational load. The following is the main conclusions of this section.

Theorem 2.1　 Let’s assume that F₁, {\hat{F}}_1, {\hat{P}}_1 are generated by Algorithm 2.2. If Assumption 2.1 holds, then the following conclusions are tenable.

(a) F_1\in {\mathrm{\Omega }}_{\alpha } and \parallel F_1-\overline{F}{\parallel }_1\le 3\parallel \tilde{F}-\overline{F}{\parallel }_1;

(b) \parallel {\hat{F}}_1-\overline{F}{\parallel }_1\le 3(1+\frac{2d}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta \kappa \alpha })\parallel \tilde{F}-\overline{F}{\parallel }_1;

(c) \parallel {\hat{P}}_1-\overline{P}{\parallel }_1\le \frac{6}{\alpha }(1+\frac{2d}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta \kappa \alpha })\parallel \tilde{F}-\overline{F}{\parallel }_1.

Proof　(a) Fix any index i\in \{1,2,\dots ,d\}. If , obviously (F_{1}e{)}_i=\alpha and F_1(i,j)\ge 0, j=1,2,\dots ,d. If (\tilde{F}e{)}_i\ge \alpha, then (F_{1}e{)}_i=(\tilde{F}e{)}_i-\frac{m_{-}}{m_{+}}\left(\right(\tilde{F}e{)}_i-\alpha ). By using 1>d\alpha and

m₊ > m₋. Hence, if (\tilde{F}e{)}_i\ge \alpha, (F_{1}e{)}_i\ge \alpha, and F_1(i,j)\ge 0, j=1,2,\dots ,d. Through direct verification it can be obtained that \parallel F_1{\parallel }_1=e^{\mathrm{T}}{F}_{1}e=1. So F_1\in {\mathrm{\Omega }}_{\alpha }. Based on Assumption 2.1 and definition of m₊, m₋, it follows

\begin{array}{rl}& \sum _{(\tilde{F}e{)}_i\ge \alpha }\sum _{j=1}^d\left|F_1\right(i,j)-\overline{F}(i,j\left)\right|\\ & =\sum _{(\tilde{F}e{)}_i\ge \alpha }\sum _{j=1}^d\left|\tilde{F}(i,j)-\overline{F}(i,j)-\frac{m_{-}}{m_{+}}\left(\right(\tilde{F}e{)}_i-\alpha )\frac{\tilde{F}(i,j)}{(\tilde{F}e{)}_i}\right|\\ & \le \sum _{(\tilde{F}e{)}_i\ge \alpha } \left(\parallel \tilde{F}(i,:)-\overline{F}(i,:){\parallel }_1+\frac{m_{-}}{m_{+}}\left(\right(\tilde{F}e{)}_i-\alpha )\right)\\ & =\sum _{(\tilde{F}e{)}_i\ge \alpha }\parallel \tilde{F}(i,:)-\overline{F}(i,:){\parallel }_1+m_{-}.\end{array}

Similarly, we can prove . With Assumption 2.1, m_{-}=\sum _{i=1}^d\mathrm{m}\mathrm{a}\mathrm{x}(\alpha -(\tilde{F}e{)}_i,0)\le \sum _{i=1}^d\mathrm{m}\mathrm{a}\mathrm{x}((\overline{F}e{)}_i-(\tilde{F}e{)}_i,0)\le \parallel \tilde{F}-\overline{F}{\parallel }_1. Based on the above,

(b) Based on the definition of {\mathrm{\Pi }}_r\left(F_1\right) and \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\overline{F}\right)\le r, it follows that \parallel F_1-{\mathrm{\Pi }}_r\left(F_1\right){\parallel }_F\le \parallel F_1-\overline{F}{\parallel }_F. Therefore, we can get after calculation based on Proposition 2.1 and result of (a) as well as relationship between l₁-norm and F(robenius) norm.

(2.3)\begin{array}{rl}\parallel {\hat{F}}_1-\overline{F}{\parallel }_1& =\parallel {\hat{F}}_1-F_1+F_1-\overline{F}{\parallel }_1\le \parallel {\hat{F}}_1-F_1{\parallel }_1+\parallel F_1-\overline{F}{\parallel }_1\\ & \le \frac{2\parallel F_1-{\mathrm{\Pi }}_r\left(F_1\right){\parallel }_1}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F_1-{\mathrm{\Pi }}_r\left(F_1\right){\parallel }_{\mathrm{\infty }}}{\alpha \kappa }+\parallel F_1-\overline{F}{\parallel }_1\\ & \le \frac{2d\parallel F_1-{\mathrm{\Pi }}_r\left(F_1\right){\parallel }_F}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F_1-{\mathrm{\Pi }}_r\left(F_1\right){\parallel }_F}{\alpha \kappa }+\parallel F_1-\overline{F}{\parallel }_1\\ & \le \frac{2d\parallel F_1-\overline{F}{\parallel }_F}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F_1-\overline{F}{\parallel }_F}{\alpha \kappa }+\parallel F_1-\overline{F}{\parallel }_1\\ & \le \frac{2d\parallel F_1-\overline{F}{\parallel }_1}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta }\frac{\parallel F_1-\overline{F}{\parallel }_1}{\alpha \kappa }+\parallel F_1-\overline{F}{\parallel }_1\\ & \le 3\left(1+\frac{2d}{\alpha (1-\kappa )}+\frac{4d\mu T}{\eta \kappa \alpha }\right)\parallel \tilde{F}-\overline{F}{\parallel }_1.\end{array}

(c) It can be seen that \overline{P}(i,:)=\frac{\overline{F}(i,:)}{\parallel \overline{F}(i,:){\parallel }_1}, i=1,2,\dots ,d. With [16, Lemma 2] and \parallel \overline{F}(i,:){\parallel }_1=(\bar{F}e{)}_i\ge \alpha, i=1,2,\dots ,d, we can get

\begin{array}{rl}\parallel {\hat{P}}_1-\overline{P}{\parallel }_1& =\sum _{i=1}^d\parallel {\hat{P}}_1(i,:)-\overline{P}(i,:){\parallel }_1=\sum _{i=1}^d{\Vert \frac{{\hat{F}}_1(i,:)}{\parallel {\hat{F}}_1(i,:){\parallel }_1}-\frac{\overline{F}(i,:)}{\parallel \overline{F}(i,:){\parallel }_1}\Vert }_1\\ & \le \sum _{i=1}^d\frac{2\parallel {\hat{F}}_1(i,:)-\overline{F}(i,:){\parallel }_1}{\parallel \overline{F}(i,:){\parallel }_1}\le \frac{2}{\alpha }\parallel {\hat{F}}_1-\overline{F}{\parallel }_1.\end{array}

With the result of (b), \parallel {\hat{P}}_1-\overline{P}{\parallel }_1\le \frac{6}{\alpha }\left(1+\frac{2d}{\alpha \left(1-\kappa \right)}+\frac{4d\mu T}{\eta \kappa \alpha }\right)\parallel \tilde{F}-\overline{F}{\parallel }_1. □

Note 2.2　(i) Theorem 2.1 shows that the estimation error between Algorithm 2.2 and the maximum likelihood estimation only differs by one multiplying constant. It is noted that the right side of equation in Theorem 2.1 (b) and (c) contains dimensional factor d of state space. That is because we directly scaled \parallel F_1-\overline{F}{\parallel }_F up to \parallel F_1-\overline{F}{\parallel }_1 in inequation (2.3). In fact, F_1-\overline{F} is a dense matrix with relatively uniform elements, so d\parallel F_1-\overline{F}{\parallel }_F=O(\parallel F_1-\overline{F}{\parallel }_1), further improving the conclusion in Theorem 2.1 (b) and (c).

(ii) With [16, Lemma7], it can be seen that there exist constants C>0 and C>1, such that the probability for \parallel \tilde{F}-\overline{F}\parallel \le C\sqrt{\frac{\parallel \overline{F}e{\parallel }_{\mathrm{\infty }}{\tau }_{\ast }{\mathrm{l}\mathrm{o}\mathrm{g}}^2\left(n\right)}{n}} to establish is at least 1-n^{-c}, where \Vert \cdot \Vert stands for the spectral norm of matrix, \parallel \cdot {\parallel }_{\mathrm{\infty }} for the infinite norm of vector, and {\tau }_{\ast } for \frac{1}{4}- mixing time of Markov chain. In this way, we can work out the statistical error bound of low-rank spectral estimation matrix generated by Algorithm 2.2 by using inequation \parallel X{\parallel }_1\le d\parallel X{\parallel }_F\le d\sqrt{d}\parallel X\parallel, \mathrm{\forall }X\in {\mathbb{R}}^{d\times d} in combination with Theorem 2.1.

3 Numerical experiments

In this section, we first compare the numerical performance of low-rank spectral estimation Algorithm 2.2 with spectral estimation [16], maximum likelihood estimation [1] and rank-constrained maximum likelihood estimation [18] through artificially synthesized data experiments. Then, we use Algorithm 2.2 combined with k-means clustering algorithm to analyze the public dataset of taxi trajectory in Manhattan Island, New York City, USA and reveal its potential transportation patterns. All numerical tests in this section have been conducted on laptops configured with Inter (R) Core (TM) i7-8550U CPU @ 1.80GHz and 8.00GB RAM by running Matlab R2019a under the Windows 10 system.

3.1 Synthetic data experiments

Firstly, the true frequency matrix and probability matrix with a rank of r are generated through Examples 3.1 and 3.2 below. Then, based on \overline{F} and \overline{P}, a Markov chain \left\{X_0,X_1,\dots ,X_n\right\}, where k\in \{1,1.5,\dots ,5.5\}, with a trajectory length of n = round(krd log(d)) is simulated through a computer.

In Algorithm 2.2, we assign \alpha =\frac{1}{d^2}, κ = 0.001, \mu =\parallel F_1^{\mathrm{T}}e{\parallel }_{\mathrm{\infty }}, \eta ={\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d} (F_1^{\mathrm{T}}e{)}_i, T = 50. Let’s assume that F and P are estimation of \overline{F} and \overline{P}, respectively, and we measure the estimation performance with the following three values:

{\eta }_F=\parallel F-\overline{F}{\parallel }_1,{\eta }_P=\frac{\parallel P-\overline{P}{\parallel }_1}{d},{\eta }_{U,V}=\mathrm{m}\mathrm{a}\mathrm{x}\left(r-\parallel {\overline{U}}^{\mathrm{T}}U{\parallel }_F^2,r-\parallel {\overline{V}}^{\mathrm{T}}V{\parallel }_F^2\right),

where, U\in {\mathbb{R}}^{d\times r}, V\in {\mathbb{R}}^{d\times r} are the matrices composed of the left and right singular vectors corresponding to the foregoing r maximum singular values of P, \overline{U}\in {\mathbb{R}}^{d\times r}, \overline{V}\in {\mathbb{R}}^{d\times r} are the matrices composed of the left and right singular vectors corresponding to the foregoing r maximum singular values of \overline{P}, and {\eta }_{U,V} can be regarded as the estimation error in the characteristic subspace. Since the rank-constrained maximum likelihood estimation [18] only provided an estimation of the state transition probability matrix \overline{P}, we only compare η_p and {\eta }_{U,V}.

Example 3.1　Generate matrix U_0\in {\mathbb{R}}^{d\times r}, V_0\in {\mathbb{R}}^{d\times r} with independent elements that obey standard normal distribution. Consider

Type I: Z=\left|U_0\right||V_0^{\mathrm{T}}|, where elements of \left|U_0\right|\in {\mathbb{R}}^{d\times r}, \left|V_0\right|\in {\mathbb{R}}^{d\times r} are absolute values of U₀, V₀ elements;

Type II: Z={\tilde{U}}_0{\tilde{V}}_0^{\mathrm{T}}, where {\tilde{U}}_0\in {\mathbb{R}}^{d\times r}, {\tilde{V}}_0\in {\mathbb{R}}^{d\times r} is defined as ({\tilde{U}}_0{)}_{ij}:=\frac{(U_0{)}_{ij}^2}{\parallel U_0(i,:){\parallel }^2}, ({\tilde{V}}_0{)}_{ij}:=\frac{(V_0{)}_{ij}^2}{\parallel V_0(:,j){\parallel }^2}.

Define \overline{F}\in {\mathbb{R}}^{d\times d} and \overline{P}\in {\mathbb{R}}^{d\times d}:\overline{F}:=\frac{Z}{\parallel Z{\parallel }_1}, \overline{P}(i,:):=\frac{\overline{F}(i,:)}{\parallel \overline{F}(i,:){\parallel }_1}, i=1,2,\dots ,d.

Fig.1 and Fig.2 present the computation results of four methods for Example 3.1. In terms of the estimation error of the frequency matrix, Algorithm 2.2 is higher than spectral estimation and maximum likelihood estimation. In terms of the estimation error of the probability matrix, when \overline{F} and \overline{P} are given by Type I, Algorithm 2.2 is optimal, followed by rank-constrained maximum likelihood estimation. When \overline{F} and \overline{P} are given by Type II, Algorithm 2.2 has similar estimation errors to rank-constrained maximum likelihood estimation and spectral estimation, both of which are lower than maximum likelihood estimation. For the estimation error of characteristic subspace, Algorithm 2.2, spectral estimation, and rank-constrained maximum likelihood estimation enjoy advantages over maximum likelihood estimation. When it comes to computational time, Algorithm 2.2, spectral estimation and maximum likelihood estimation are all less than rank-constrained maximum likelihood estimation because the former three methods do not require solving optimization problem while the latter requires using DC programming to solve a matrix rank-constrained optimization problem.

Example 3.2　Generate matrices U_0\in {\mathbb{R}}^{d\times r}, V_0\in {\mathbb{R}}^{d\times r} with independent elements that obey standard normal distribution and diagonal matrix D\in {\mathbb{R}}^{d\times d}, where diagonal elements obey distribution beta(γ, γ). Consider

Type I: Z=\left|U_0\right||V_0^{\mathrm{T}}|D, where elements of \left|U_0\right|\in {\mathbb{R}}^{d\times r}, \left|V_0\right|\in {\mathbb{R}}^{d\times r} are absolute values of U₀, V₀ elements;

Type II: Z={\tilde{U}}_0{\tilde{V}}_0^{\mathrm{T}}D, where {\tilde{U}}_0\in {\mathbb{R}}^{d\times r}, {\tilde{V}}_0\in {\mathbb{R}}^{d\times r} is defined as ({\tilde{U}}_0{)}_{ij}:=\frac{(U_0{)}_{ij}^2}{\parallel U_0(i,:){\parallel }^2}, ({\tilde{V}}_0{)}_{ij}:=\frac{(V_0{)}_{ij}^2}{\parallel V_0(:,j){\parallel }^2} .

Define \overline{F}\in {\mathbb{R}}^{d\times d} and \mid \overline{P}\in {\mathbb{R}}^{d\times d}:\overline{F}:=\frac{Z}{\parallel Z{\parallel }_1}, \overline{P}(i,:):=\frac{\overline{F}(i,:)}{\parallel \overline{F}(i,:){\parallel }_1}, i=1,2,\dots ,d.

It is noted that \rho :=\frac{{\mathrm{m}\mathrm{a}\mathrm{x}}_{1\le i\le d} {\pi }_i}{{\mathrm{m}\mathrm{i}\mathrm{n}}_{1\le i\le d} {\pi }_i}., where \pi =\overline{F}e, which can be used to characterize the equilibrium of the stationary distribution of the Markov model. The smaller the value of ρ, the more equilibrium the stationary distribution of the Markov model has. Otherwise, there is more non-equilibrium. Compared with Example 3.1, Example 3.2 can be seen as a Markov model with a “non-equilibrium” stationary distribution. Fig.3 and Fig.4 show the computational results of the four methods for Example 3.2, where γ = 0.1. For the estimation error of the frequency matrix, Algorithm 2.2 does not significantly differ from spectral estimation and is lower than maximum likelihood estimation; for the estimation error of probability matrix, Algorithm 2.2 and rank-constrained maximum likelihood estimation are lower than spectral estimation and maximum likelihood estimation; for the estimation error of the characteristic subspace, rank-constrained maximum likelihood estimation is optimal. Furthermore, in order to examine the numerical performance of spectral estimation, rank-constrained maximum likelihood estimation and Algorithm 2.2 in the Markov model’s stationary distribution shifting from non-equilibrium to equilibrium, we let γ = 0.1,0.2,\dots ,1 (beta (1, 1) is uniform distribution on [0, 1]) in Example 3.2.

Fig.5 shows the computational results of three three methods. It can be seen that when γ is relatively small, rank-constrained maximum likelihood estimation is optimal; when γ gradually increases, the stationary distribution of the Markov model tends to equilibrium, with the differences between the three methods decreasing and the estimation errors displaying a lowering trend. In fact, as the stationary distribution of the Markov model tends to equilibrium, let’s assume the constant α in Assumption 2.1 increases, it can be inferred that the estimation error of Algorithm 2.2 decreases.

In summary, for Markov models with equilibrium stationary distribution, Algorithm 2.2 can quickly obtain high-quality estimation of low-rank frequency and low-rank probability matrices. For Markov models with non-equilibrium stationary distribution, rank-constrained maximum likelihood estimation is more effective, but its computational complexity is too high. Therefore, Algorithm 2.2 is suitable for solving the learning problem of high-dimension and low-rank Markov models.

3.2 Movement trajectory analysis of taxis in Manhattan Island, New York City, USA

The yellow taxi operation trajectory dataset of Manhattan Island, New York City, USA made public in 2016 recorded the itinerary of 1.1×10⁷ passengers. Our goal is to divide Manhattan Island into several regions based on this dataset, so that passengers in the same region have similar preferences for traveling to the same destination. By such partitioning, we can assist in optimizing the balance of urban transportation supply and demand and resource distribution.

Similar to [16], Manhattan Island is firstly subdivided into multiple 0.001° × 0.001° square grids according to longitude and latitude (40.70°−40.82° N, 73.93°−74.02° W). Then, each small grid is approximately taken as a state of Markov process, and the journey of passengers from one grid to another can be regarded as a state transition of the process. To eliminate interference, we only select grids that have been the starting point more than 1000 times as effective states. Then, low-rank spectral estimation Algorithm 2.2 is used to estimate the state transition probability matrix of the Markov process, and k-means clustering method to estimate the state transition probability matrix of the Markov process. The left singular subspace of the matrix is used for clustering-based partitioning. Fig.6-Fig.8 show the clustering results when k = r is 5, 6, and 7, respectively.

Considering that people have different travel purposes at different time periods, we divide the taxi operation dataset into three time periods, namely morning (06:00−11:59), afternoon (12:00−17:59), and evening (18:00−23:59). Under this time period division, the effective number of states for each time period can be obtained as 803, 1036, and 1130. We take k = r = 5 and use Algorithm 2.2 to estimate the probability matrix of state transition for different time periods and cluster them. The clustering results are shown in Fig.9 (morning), Fig.10 (afternoon), and Fig.11 (evening).

4 Summary

This paper proposes a low-rank spectral estimation method for the state transition frequency matrix and probability matrix estimation of Markov models. With the help of approximate projection algorithm for rank-constrained frequency matrix set, our method can obtain low-rank estimation of the real matrix, thus overcoming the limitation of the spectral estimation method proposed in [16] that cannot meet the low-rank condition. In addition, our method does not require matrix optimization, so the computational complexity is very low and can be applied to large low-rank Markov model learning. We have tested synthetic data and Manhattan Island taxi movement trajectory data in New York City, USA, and numerical experiments show that the low-rank spectral estimation method is effective. Next, we will study the relationship between the statistical error bound and the lower bound of the minimax error in low-rank spectral estimation algorithm.

References

[1] Anderson T W, Goodman L A. Statistical inference about Markov chains. Ann Math Statist 1957; 28(1): 89–110

[2] BenczúrA ACsalogányKSarlósT. On the feasibility of low-rank approximation for personalized PageRank. In: Special Interest Tracks and Posters of the 14th International Conference on World Wide Web. New York: ACM, 2005, 972–973

[3] Benson A R, Gleich D F, Lim L-H. The spacey random walk: a stochastic process for higher-order data. SIAM Rev 2017; 59(2): 321–345

[4] Bi S J, Pan S H. Error bounds for rank constrained optimization problems and applications. Oper Res Lett 2016; 44(3): 336–341

[5] HuangQ QKakadeS MKongW HValiantG. Recovering structured probability matrices. In: 9th Innovations in Theoretical Computer Science, LIPIcs Leibniz Int Proc Inform, Vol 94. Wadern: Schloss Dagstuhl Leibniz-Zent Inform, 2018, 46

[6] JainPMekaRDhillonI. Guaranteed rank minimization via singular value projection. In: Advances in Neural Information Processing Systems 23 (24th Annual Conference on Neural Information Processing Systems 2010), Vol 1. La Jolla, CA: NIPS, 2010, 937–945

[7] Li G Y, Mordukhovich B S, Nghia T T A, Pham T S. Error bounds for parametric polynomial systems with applications to higher-order stability analysis and convergence rates. Math Program 2018; 168(1/2): Ser B, 313–346

[8] Li Y H, Ren T Y, Shao P D, Song W P, Li Z Y, Gou Q Z. Development of driving cycle of bus in Xi’an City based on Markov chain. China Sciencepaper 2019; 14(2): 121–128

[9] Liu Y, Kang C G, Gao S, Xiao Y, Tian Y. Understanding intra-urban trip patterns from taxi trajectory data. J Geographical Systems 2012; 14(4): 463–483

[10] Negahban S, Oh S, Shah D. Rank centrality: ranking from pairwise comparisons. Oper Res 2017; 65(1): 266–287

[11] Pang J-S. Error bounds in mathematical programming. Mathematical Programming 1997; 79(1/2/3): Ser B, 299–332

[12] Sanders J, Proutière A, Yun S-Y. Clustering in block Markov chains. Ann Statist 2020; 48(6): 3488–3512

[13] Tseng P. Approximation accuracy, gradient methods, and error bound for structured convex optimization. Math Program 2010; 125(2): Ser. B, 263–295

[14] Xu X B, Li Z P. Research progress and prospects for application of reinforcement learning in operations research. Operations Research and Management Science 2020; 29(5): 227–239

[15] Zhai Y, Liu J, Chen J, Xing X C, Du J, Li H, Zhu J. Analysis of dockless bike-sharing fleet size based on Markov chain. J Beijing Jiaotong Univ 2019; 43(5): 27–36

[16] Zhang A R, Wang M D. Spectral state compression of Markov processes. IEEE Trans Inform Theory 2020; 66(5): 3202–3231

[17] Zhou Z R, So A M-C. A unified approach to error bounds for structured convex optimization problems. Math Program 2017; 165(2): Ser A, 689–728

[18] Zhu Z W, Li X D, Wang M D, Zhang A R. Learning Markov models via low-rank optimization. Oper Res 2022; 70(4): 2384–2398

RIGHTS & PERMISSIONS

Higher Education Press 2024