賁樹軍， 翁藝鴻

(華南理工大學(xué)數(shù)學(xué)學(xué)院， 廣州 510641)

馬爾可夫過程的估計(jì)問題是計(jì)算機(jī)科學(xué)、系統(tǒng)工程和數(shù)據(jù)科學(xué)等領(lǐng)域中的一個核心問題[1]。計(jì)算個性化網(wǎng)頁排序的狀態(tài)轉(zhuǎn)移矩陣問題[2]、解決電子商務(wù)中的排序問題[3]和分析城市出租車或公交車的運(yùn)行軌跡問題[4-5]等都可歸結(jié)為馬爾可夫過程的估計(jì)問題。源于上述問題的馬爾可夫過程往往擁有很大的狀態(tài)空間，但是它們的狀態(tài)轉(zhuǎn)移矩陣卻被證明了是低秩或者近似低秩的矩陣[1]。因此，學(xué)者們對低秩馬爾可夫過程的狀態(tài)轉(zhuǎn)移矩陣的估計(jì)及其應(yīng)用問題開展了研究[1,6-10]。

據(jù)我們所知，現(xiàn)有的估計(jì)方法都不能保證得到低秩的轉(zhuǎn)移矩陣估計(jì)。譬如，ZHANG和WANG[1]利用頻率矩陣的經(jīng)驗(yàn)估計(jì)的截?cái)嗥娈愔捣纸饨Y(jié)合非負(fù)投影，提出了低秩馬爾可夫過程的譜估計(jì)方法，并建立了估計(jì)的統(tǒng)計(jì)誤差界，證明了估計(jì)誤差與極小極大誤差的下界相差一個馬爾可夫鏈軌跡長度的對數(shù)因子。但是，由于譜估計(jì)方法利用了非負(fù)投影，導(dǎo)致該文獻(xiàn)最后得到的估計(jì)矩陣不是低秩的。ZHU等[8]提出了狀態(tài)轉(zhuǎn)移矩陣的核范數(shù)正則罰極大似然估計(jì)模型和秩約束極大似然估計(jì)模型，并建立了2種模型的統(tǒng)計(jì)誤差界，證明了估計(jì)誤差與極小極大誤差的下界相差一個馬爾可夫鏈狀態(tài)空間維數(shù)的對數(shù)因子。然而，核范數(shù)正則優(yōu)化問題的最優(yōu)解不一定滿足低秩條件，秩約束優(yōu)化問題的最優(yōu)解雖能滿足秩約束條件，但其求解一般都是NP難。盡管文獻(xiàn)[8]設(shè)計(jì)了一類DC (凸函數(shù)的差) 規(guī)劃算法來近似求解秩約束極大似然估計(jì)模型，但不能保證算法的輸出是一個低秩矩陣。特別地，該DC規(guī)劃算法的每一步都需要進(jìn)行奇異值分解，計(jì)算量非常大，因此不適用于大規(guī)模的馬爾可夫過程估計(jì)問題。

另一方面，誤差界研究一直以來都是最優(yōu)化領(lǐng)域中的重點(diǎn)和難題[11-15]。PANG[11]證明了：凸多面體集合具有全局Lipschitz型誤差界；在一定的約束規(guī)范下，一般凸不等式系統(tǒng)具有Lipschitz型誤差界；對一般非凸不等式系統(tǒng)，全局誤差界(即使是H?lder型的)都很難成立。目前已有的研究主要是對次解析不等式系統(tǒng)和多項(xiàng)式不等式系統(tǒng)建立了局部H?lder型誤差界[14]，對矩陣秩約束系統(tǒng)的誤差界研究還很少。對于一個有界且其多值函數(shù)在原點(diǎn)滿足calmness條件的矩陣秩約束系統(tǒng)，BI和PAN[15]得到了該系統(tǒng)的局部和全局Lipschitz型誤差界。但是，驗(yàn)證多值函數(shù)的calmness條件與建立誤差界的難度基本一樣。因此，我們需要尋求新的工具來研究矩陣秩約束系統(tǒng)的誤差界。

受上述啟發(fā)，本文試圖尋求一個能夠快速獲得低秩轉(zhuǎn)移矩陣的方法，以估計(jì)大規(guī)模的低秩馬爾可夫過程。首先，建立秩約束狀態(tài)轉(zhuǎn)移矩陣集合的局部Lipschitz型誤差界，并尋求該集合高質(zhì)量的近似投影方法；然后，提出一種低秩馬爾可夫過程狀態(tài)轉(zhuǎn)移矩陣的低秩譜估計(jì)算法(LRSEA)，并進(jìn)行數(shù)值實(shí)驗(yàn)。

1 局部Lipschitz型誤差界

本節(jié)為秩約束狀態(tài)轉(zhuǎn)移矩陣集合建立局部Lipschitz型誤差界。首先，給出本文的記號。

(1)記Id×d、Ed×d、ed分別為單位矩陣、元素全為1的矩陣、元素全為1的向量。

(2)對于Pd×d，定義‖P‖∞為無窮范數(shù)，為Frobenius范數(shù)。

(3)秩約束狀態(tài)轉(zhuǎn)移矩陣集合定義為：Π∶={Pd×d:rank(P)≤r,Pe=e,Pij≥0,1≤i,j≤d}，其中r[1,d-1]是一個給定整數(shù)。對于Zd×d，定義Z到集合Π的距離為:

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

(4)給定xd，定義l1范數(shù)為無窮范數(shù)為‖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)對應(yīng)的左、右奇異向量。眾所周知，在Frobenius范數(shù)距離意義下，Υ(P)是P在上的一個投影矩陣。

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

下面給出建立矩陣PΩ與投影矩陣Υ(P)之間關(guān)系的引理。

引理1令γ(0,1)是給定常數(shù)。任取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)為給定常數(shù)。對任意矩陣PΩ,定義如下矩陣

(4)

下面給出建立集合Π的局部Lipschitz型誤差界的命題。

命題1令γ(0,1/2],c是給定常數(shù)。任取PΩ，有ΠΠ，且

(5)

證明設(shè)PΩ,下面分5種情況證明。

(1)假設(shè)‖Υ(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)，進(jìn)而得到定義，由rank()≤r，可得rank(Π)≤r，Π≥0，Πe=e，即ΠΠ。根據(jù)假設(shè)(eTDΥ(P))j≥cα(1≤j≤d)，可得∞。利用假設(shè)β/c≥‖eTDΥ(P)‖∞，可得

(6)

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

(7)

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

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

即不等式(5)成立。

假設(shè)‖Υ(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)假設(shè)‖Υ(P)‖∞≤10‖P‖∞,且由于eTDΥ(P)中存在大于β/c的元素，不妨設(shè)(eTDΥ(P))i>β/c。由(eTP)i≤β，可得

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

(9)

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

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

2 LRSEA算法

首先給出本文的假設(shè)。

假設(shè)1存在常數(shù)β[1,d),α(0,1]，使得:

(2){X0,X1,…,Xn}是遍歷馬爾可夫鏈，平穩(wěn)分布為πd，滿足πi≥α/d(i=1,2,…,d),其中πi是向量π的第i個元素。

算法1 LRSEA算法

第1步：計(jì)算修正的譜估計(jì)

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)得到的低秩譜估計(jì)Π。

(a)SΩ且

其中E表示數(shù)學(xué)期望。

且

由上述討論，可得

同理可以證明

所以，

(b)若假設(shè)1成立，則由文獻(xiàn)[1]的定理1可知：存在一個常數(shù)C，使得

由(a)部分的結(jié)論，可知存在常數(shù)C，使得

(10)

由于SΩ，由命題1可知

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

(11)

由不等式(10)、(11)可得定理結(jié)論。證畢。

3 數(shù)值實(shí)驗(yàn)

3.1 合成數(shù)據(jù)實(shí)驗(yàn)

首先，考慮具有平衡分布的低秩馬爾可夫過程估計(jì)問題。假設(shè)U0d×r,V0d×r，它們的元素由標(biāo)準(zhǔn)正態(tài)分布隨機(jī)生成。定義矩陣d×r,這里[i,:]表示的第i行，[:,j]表示的第j列，⊙表示矩陣Hadamard積。定義真實(shí)狀態(tài)轉(zhuǎn)移矩陣為T。本文利用狀態(tài)轉(zhuǎn)移矩陣生成狀態(tài)數(shù)為d、長度為n=round(qdr(logd)2)的馬爾可夫鏈X0,…,Xn，這里q是常數(shù)。

下面比較文獻(xiàn)[17]的經(jīng)驗(yàn)估計(jì)方法、文獻(xiàn)[1]的譜估計(jì)方法和LRSEA算法的估計(jì)誤差。記P為相應(yīng)的估計(jì)矩陣，本文利用以下2個數(shù)值來衡量其估計(jì)效果：

其中，Ud×r、Vd×r分別是P的前r個最大奇異值對應(yīng)的左奇異向量、右奇異向量，d×r、d×r分別是的前r個最大奇異值對應(yīng)的左奇異向量、右奇異向量。

取d=1 000、r=10、k[1,10]時(shí)，3種方法的估計(jì)效果(圖1)表明：LRSEA算法與譜估計(jì)方法的估計(jì)誤差相差不大，小于經(jīng)驗(yàn)估計(jì)方法的。

圖1 平衡分布下3種估計(jì)方法的比較

圖2 非平衡分布下3種估計(jì)方法的比較

3.2 分析曼哈頓出租車運(yùn)行軌跡

紐約市曼哈頓島于2016年公開的黃色出租車運(yùn)行軌跡數(shù)據(jù)集記錄了1.1×107次乘客的行程(https:∥s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2016 -01.csv)，本研究利用此數(shù)據(jù)集將曼哈頓島分割成幾個區(qū)域，滿足同一區(qū)域中的乘客前往同一個目的地的概率相似。

類似文獻(xiàn)[1]，本文將曼哈頓島細(xì)分為大小一致的小正方形網(wǎng)格，將每個小網(wǎng)格近似地看成馬爾可夫過程的一個狀態(tài)，乘客從一個網(wǎng)格到另一個網(wǎng)格的行程視為該馬爾可夫過程的一次狀態(tài)轉(zhuǎn)移。為了排除干擾項(xiàng)，本文只選取那些作為行程的起點(diǎn)或者終點(diǎn)的總次數(shù)超過2 000的網(wǎng)格作為有效狀態(tài)。然后，利用LRSEA算法來估計(jì)該馬爾可夫過程的狀態(tài)轉(zhuǎn)移矩陣，并利用k-均值聚類方法對估計(jì)矩陣的左奇異子空間進(jìn)行聚類劃分。由k=r分別為4,5,6,7時(shí)的聚類結(jié)果(圖3)可知：當(dāng)聚類數(shù)r增加時(shí)，LRSEA算法可以給曼哈頓島的交通網(wǎng)絡(luò)一個較好的分區(qū)，同一個區(qū)域的乘客前往同一個地點(diǎn)的概率相似。

圖3 LRSEA算法對曼哈頓島交通網(wǎng)絡(luò)的劃分結(jié)果

鑒于在不同時(shí)段，人們出行目的不同，將該出租車運(yùn)行數(shù)據(jù)集化分為3個時(shí)段：早上(06:00～11:59)、中午(12:00～17:59)、晚上(18:00～23:59)，每個時(shí)段的有效狀態(tài)數(shù)分別為769、1 029、1 147個。利用LRSEA算法來估計(jì)每個時(shí)段的馬爾可夫過程的狀態(tài)轉(zhuǎn)移矩陣，并利用k-均值聚類方法對估計(jì)矩陣的左奇異子空間進(jìn)行聚類劃分，其中k=r=5。聚類結(jié)果所展示的在不同時(shí)段下LRSEA算法對曼哈頓島的交通網(wǎng)絡(luò)的分區(qū)結(jié)果(圖4)表明：同一時(shí)段下，同一個區(qū)域的乘客前往同一個目的地的概率相似。

圖4 LRSEA算法在不同時(shí)間段下對曼哈頓島交通網(wǎng)絡(luò)的劃分結(jié)果

4 結(jié)論

針對馬爾可夫過程的估計(jì)問題，利用秩約束狀態(tài)轉(zhuǎn)移矩陣集合的近似投影，本文對現(xiàn)有的譜方法進(jìn)行低秩修正，提出了一個低秩譜估計(jì)算法(LRSEA)，以快速得到滿足秩約束條件的狀態(tài)轉(zhuǎn)移矩陣。此外，通過建立秩約束狀態(tài)轉(zhuǎn)移矩陣集合的局部Lipschitz型誤差界，給出該算法的統(tǒng)計(jì)誤差界，建立了算法的理論保證。數(shù)值實(shí)驗(yàn)結(jié)果表明，對于具有非平衡分布的低秩馬爾可夫過程的估計(jì)問題，LRSEA算法的估計(jì)誤差小于譜估計(jì)方法和經(jīng)驗(yàn)估計(jì)方法的。下一步，將把LRSEA算法應(yīng)用到強(qiáng)化學(xué)習(xí)問題以及系統(tǒng)工程領(lǐng)域中的控制問題。