馬爾可夫過程的低秩譜估計

賁樹軍， 翁藝鴻

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

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

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

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

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

1 局部Lipschitz型誤差界

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

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

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

(3)秩約束狀態(tài)轉移矩陣集合定義為：Π∶={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)對應的左、右奇異向量。眾所周知，在Frobenius范數(shù)距離意義下，Υ(P)是P在上的一個投影矩陣。

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

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

引理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)

證明設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，即ΠΠ。根據(jù)假設(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存在常數(shù)β[1,d),α(0,1]，使得:

(2){X0,X1,…,Xn}是遍歷馬爾可夫鏈，平穩(wěn)分布為π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表示數(shù)學期望。

且

由上述討論，可得

同理可以證明

所以，

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

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

(10)

由于SΩ，由命題1可知

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

(11)

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

3 數(shù)值實驗

3.1 合成數(shù)據(jù)實驗

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

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

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

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

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

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

3.2 分析曼哈頓出租車運行軌跡

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

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

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

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

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

4 結論

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