马尔可夫过程的低秩谱估计

贲树军， 翁艺鸿

(华南理工大学数学学院， 广州 510641)

马尔可夫过程的估计问题是计算机科学、系统工程和数据科学等领域中的一个核心问题[1]。计算个性化网页排序的状态转移矩阵问题[2]、解决电子商务中的排序问题[3]和分析城市出租车或公交车的运行轨迹问题[4-5]等都可归结为马尔可夫过程的估计问题。源于上述问题的马尔可夫过程往往拥有很大的状态空间，但是它们的状态转移矩阵却被证明了是低秩或者近似低秩的矩阵[1]。因此，学者们对低秩马尔可夫过程的状态转移矩阵的估计及其应用问题开展了研究[1,6-10]。

据我们所知，现有的估计方法都不能保证得到低秩的转移矩阵估计。譬如，ZHANG和WANG[1]利用频率矩阵的经验估计的截断奇异值分解结合非负投影，提出了低秩马尔可夫过程的谱估计方法，并建立了估计的统计误差界，证明了估计误差与极小极大误差的下界相差一个马尔可夫链轨迹长度的对数因子。但是，由于谱估计方法利用了非负投影，导致该文献最后得到的估计矩阵不是低秩的。ZHU等[8]提出了状态转移矩阵的核范数正则罚极大似然估计模型和秩约束极大似然估计模型，并建立了2种模型的统计误差界，证明了估计误差与极小极大误差的下界相差一个马尔可夫链状态空间维数的对数因子。然而，核范数正则优化问题的最优解不一定满足低秩条件，秩约束优化问题的最优解虽能满足秩约束条件，但其求解一般都是NP难。尽管文献[8]设计了一类DC (凸函数的差) 规划算法来近似求解秩约束极大似然估计模型，但不能保证算法的输出是一个低秩矩阵。特别地，该DC规划算法的每一步都需要进行奇异值分解，计算量非常大，因此不适用于大规模的马尔可夫过程估计问题。

另一方面，误差界研究一直以来都是最优化领域中的重点和难题[11-15]。PANG[11]证明了：凸多面体集合具有全局Lipschitz型误差界；在一定的约束规范下，一般凸不等式系统具有Lipschitz型误差界；对一般非凸不等式系统，全局误差界(即使是Hölder型的)都很难成立。目前已有的研究主要是对次解析不等式系统和多项式不等式系统建立了局部Hölder型误差界[14]，对矩阵秩约束系统的误差界研究还很少。对于一个有界且其多值函数在原点满足calmness条件的矩阵秩约束系统，BI和PAN[15]得到了该系统的局部和全局Lipschitz型误差界。但是，验证多值函数的calmness条件与建立误差界的难度基本一样。因此，我们需要寻求新的工具来研究矩阵秩约束系统的误差界。

受上述启发，本文试图寻求一个能够快速获得低秩转移矩阵的方法，以估计大规模的低秩马尔可夫过程。首先，建立秩约束状态转移矩阵集合的局部Lipschitz型误差界，并寻求该集合高质量的近似投影方法；然后，提出一种低秩马尔可夫过程状态转移矩阵的低秩谱估计算法(LRSEA)，并进行数值实验。

1 局部Lipschitz型误差界

本节为秩约束状态转移矩阵集合建立局部Lipschitz型误差界。首先，给出本文的记号。

(1)记Id×d、Ed×d、ed分别为单位矩阵、元素全为1的矩阵、元素全为1的向量。

(2)对于Pd×d，定义‖P‖∞为无穷范数，为Frobenius范数。

(3)秩约束状态转移矩阵集合定义为：Π∶={Pd×d:rank(P)≤r,Pe=e,Pij≥0,1≤i,j≤d}，其中r[1,d-1]是一个给定整数。对于Zd×d，定义Z到集合Π的距离为:

dist(Z,Π)∶=min{‖Z-P‖F:PΠ}。

(4)给定xd，定义l1范数为无穷范数为‖x‖∞记diag(x)d×d是第i个对角元为xi(i=1,2,…,d)的对角矩阵。

(5)定义秩r约束矩阵集合为∶={Zd×d:rank(Z)≤r}。对任意Pd×d，定义

Υ

其中，σi(P)(i=1,2,…,r)是P的第i个最大奇异值，ui(P)d、vi(P)d分别是σi(P)对应的左、右奇异向量。众所周知，在Frobenius范数距离意义下，Υ(P)是P在上的一个投影矩阵。

(6)定义集合Ω∶={Zd×d:Ze=e,β/d≥Zij≥α/d,1≤i,j≤d}，其中α(0,1)和β(1,d)是给定常数。

下面给出建立矩阵PΩ与投影矩阵Υ(P)之间关系的引理。

引理1令γ(0,1)是给定常数。任取PΩ，若(Υ(P)e)i≥γ(i=1,2,…,d)，则有：

(1)

‖P-DΥ(P)‖∞∞,

(2)

(3)

证明由于(Υ(P)e)j≥γ>0 (j=1,2,…,d)且Pe=e，利用〈x,y〉≤‖x‖1‖y‖∞(x,yd)，有

因此，不等式(1)成立。利用不等式(1)、I-D是对角矩阵以及d‖P‖∞≥1>γ，可得

‖P-DΥ(P)‖∞=‖P-Υ(P)+(I-D)Υ(P)‖∞≤

‖P-Υ(P)‖∞+‖I-D‖∞‖Υ(P)‖∞≤

‖P-Υ(P)‖∞∞≤

不等式(3)成立。证毕。

令γ>0,c(0,1)为给定常数。对任意矩阵PΩ,定义如下矩阵

(4)

下面给出建立集合Π的局部Lipschitz型误差界的命题。

命题1令γ(0,1/2],c是给定常数。任取PΩ，有ΠΠ，且

(5)

证明设PΩ,下面分5种情况证明。

(1)假设‖Υ(P)‖∞≤10‖P‖∞,(Υ(P)e)j≥γ且β/c≥(eTDΥ(P))j≥cα(1≤j≤d)。由和ti的定义，可得il=(DΥ(P))il-ti(eTDΥ(P))l≥0(1≤i,l≤d)。注意-ti(eTDΥ(P))l≥0，可得il≥(DΥ(P))il(1≤i,l≤d)，因此(e)i≥(DΥ(P)e)i=1 (i=1,2,…,d)，进而得到定义，由rank()≤r，可得rank(Π)≤r，Π≥0，Πe=e，即ΠΠ。根据假设(eTDΥ(P))j≥cα(1≤j≤d)，可得∞。利用假设β/c≥‖eTDΥ(P)‖∞，可得

(6)

‖Υ(P)‖∞∞‖-DΥ(P)‖∞

(7)

因此，利用d‖P‖∞≤β及不等式(2)、(3)、(7)，可得

‖P-DΥ(P)‖∞∞，

即不等式(5)成立。

假设‖Υ(P)‖∞>10‖P‖∞或(Υ(P)e)<γ或或由式(4)可知Π=E/d。显然ΠΠ且

‖Υ(P)-P‖∞≥9‖P‖∞

不等式(5)成立。

d‖P-Υ(P)‖∞≥‖(P-Υ(P))e‖∞≥

(Pe)i-(Υ(P)e)i≥1-γ≥0.5，

则

2β‖P-Υ(P)‖∞。

由此可得不等式(5)成立。

d‖P-DΥ(P)‖∞≥‖eT(P-DΥ(P))‖∞≥

(eTP)i-(eTDΥ(P))i≥α-cα。

(8)

由不等式(8)、d‖P‖∞≤β及引理1，可得

由c可得不等式(5)成立。

(5)假设‖Υ(P)‖∞≤10‖P‖∞,且由于eTDΥ(P)中存在大于β/c的元素，不妨设(eTDΥ(P))i>β/c。由(eTP)i≤β，可得

d‖P-DΥ(P)‖∞≥‖eT(P-DΥ(P))‖∞≥

(9)

由不等式(9)、d‖P‖∞≤β及引理1,可得

由c可得不等式(5)成立。证毕。

2 LRSEA算法

首先给出本文的假设。

假设1存在常数β[1,d),α(0,1]，使得:

(2){X0,X1,…,Xn}是遍历马尔可夫链，平稳分布为πd，满足πi≥α/d(i=1,2,…,d),其中πi是向量π的第i个元素。

算法1 LRSEA算法

第1步：计算修正的谱估计

fori=1 tod

forj=1 tod

end for

else

end if

end for

fori=1 tod

forj=1 tod

end for

else

end if

end for

第2步：令P=S,γ(0,1/2],c由式(4)得到的低秩谱估计Π。

(a)SΩ且

其中E表示数学期望。

且

由上述讨论，可得

同理可以证明

所以，

(b)若假设1成立，则由文献[1]的定理1可知：存在一个常数C，使得

由(a)部分的结论，可知存在常数C，使得

(10)

由于SΩ，由命题1可知

C2‖S-Υ(S)‖∞‖1≤C2‖S-Υ(S)‖F+

(11)

由不等式(10)、(11)可得定理结论。证毕。

3 数值实验

3.1 合成数据实验

首先，考虑具有平衡分布的低秩马尔可夫过程估计问题。假设U0d×r,V0d×r，它们的元素由标准正态分布随机生成。定义矩阵d×r,这里[i,:]表示的第i行，[:,j]表示的第j列，⊙表示矩阵Hadamard积。定义真实状态转移矩阵为T。本文利用状态转移矩阵生成状态数为d、长度为n=round(qdr(logd)2)的马尔可夫链X0,…,Xn，这里q是常数。

下面比较文献[17]的经验估计方法、文献[1]的谱估计方法和LRSEA算法的估计误差。记P为相应的估计矩阵，本文利用以下2个数值来衡量其估计效果：

其中，Ud×r、Vd×r分别是P的前r个最大奇异值对应的左奇异向量、右奇异向量，d×r、d×r分别是的前r个最大奇异值对应的左奇异向量、右奇异向量。

取d=1 000、r=10、k[1,10]时，3种方法的估计效果(图1)表明：LRSEA算法与谱估计方法的估计误差相差不大，小于经验估计方法的。

图1 平衡分布下3种估计方法的比较

图2 非平衡分布下3种估计方法的比较

3.2 分析曼哈顿出租车运行轨迹

纽约市曼哈顿岛于2016年公开的黄色出租车运行轨迹数据集记录了1.1×107次乘客的行程(https:∥s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2016 -01.csv)，本研究利用此数据集将曼哈顿岛分割成几个区域，满足同一区域中的乘客前往同一个目的地的概率相似。

类似文献[1]，本文将曼哈顿岛细分为大小一致的小正方形网格，将每个小网格近似地看成马尔可夫过程的一个状态，乘客从一个网格到另一个网格的行程视为该马尔可夫过程的一次状态转移。为了排除干扰项，本文只选取那些作为行程的起点或者终点的总次数超过2 000的网格作为有效状态。然后，利用LRSEA算法来估计该马尔可夫过程的状态转移矩阵，并利用k-均值聚类方法对估计矩阵的左奇异子空间进行聚类划分。由k=r分别为4,5,6,7时的聚类结果(图3)可知：当聚类数r增加时，LRSEA算法可以给曼哈顿岛的交通网络一个较好的分区，同一个区域的乘客前往同一个地点的概率相似。

图3 LRSEA算法对曼哈顿岛交通网络的划分结果

鉴于在不同时段，人们出行目的不同，将该出租车运行数据集化分为3个时段：早上(06:00～11:59)、中午(12:00～17:59)、晚上(18:00～23:59)，每个时段的有效状态数分别为769、1 029、1 147个。利用LRSEA算法来估计每个时段的马尔可夫过程的状态转移矩阵，并利用k-均值聚类方法对估计矩阵的左奇异子空间进行聚类划分，其中k=r=5。聚类结果所展示的在不同时段下LRSEA算法对曼哈顿岛的交通网络的分区结果(图4)表明：同一时段下，同一个区域的乘客前往同一个目的地的概率相似。

图4 LRSEA算法在不同时间段下对曼哈顿岛交通网络的划分结果

4 结论

针对马尔可夫过程的估计问题，利用秩约束状态转移矩阵集合的近似投影，本文对现有的谱方法进行低秩修正，提出了一个低秩谱估计算法(LRSEA)，以快速得到满足秩约束条件的状态转移矩阵。此外，通过建立秩约束状态转移矩阵集合的局部Lipschitz型误差界，给出该算法的统计误差界，建立了算法的理论保证。数值实验结果表明，对于具有非平衡分布的低秩马尔可夫过程的估计问题，LRSEA算法的估计误差小于谱估计方法和经验估计方法的。下一步，将把LRSEA算法应用到强化学习问题以及系统工程领域中的控制问题。