宋素華，喻高航

(杭州電子科技大學理學院，浙江 杭州 310018)

0 引 言

馬爾科夫鏈(Markov Chain)[1]是狀態(tài)空間中從一個狀態(tài)轉(zhuǎn)移到另一個狀態(tài)的隨機過程，廣泛應用于DNA序列分析[2]、多重線性Pagerank[3]、風力機設計[4]、排隊系統(tǒng)[5]等領域，是分析隨時間變化的各種隨機(概率)過程的有力工具。高階馬爾科夫鏈是一階馬爾科夫鏈的高階推廣，近年來，關于高階馬爾可夫鏈的極限概率分布問題的相關研究取得了不少成果。文獻[6]提出一種快速求解稀疏馬爾可夫鏈極限概率分布的算法。文獻[7]研究了用高階冪法(Higher Order Power Method，HOPM)求解馬爾科夫鏈的極限概率分布問題，給出了算法的收斂性分析，并在一定條件下建立了高階馬爾可夫鏈的極限概率分布向量的存在性和唯一性。文獻[8]在高階冪法的基礎上添加松弛項，提出一種松弛高階冪法(Relaxation Higher Order Power Method，RHOPM)，提高了原有冪法的計算效率。本文主要針對由高階馬爾科夫鏈產(chǎn)生的極限概率分布問題，結合動量梯度法，提出一種新的帶動量項的松弛算法(簡稱HOPM-NAG)，并對所提算法進行收斂性分析。

1 預備知識

一個由nm個實數(shù)組成的m階n維實張量，

其中[n]={1,2,…,n},R為實數(shù)集，如果pi1i2…im≥0，則稱為非負張量。當m=2時，上述張量就退化為矩陣。

任意一個n維列向量x∈Rn，張量-向量積定義如下：

(xm-1)i

對于任意2個n維向量x,y∈Rn，

1.1 一階馬爾科夫鏈

馬爾科夫鏈是一組具有馬爾科夫性質(zhì)的離散隨機變量的集合，若集合中某一個時刻的狀態(tài)只依賴于它前一個時刻的狀態(tài)，則稱其為一階馬爾科夫鏈[9]。給定一個隨機過程{Xt,t=0,1,2,…}其取值屬于一個代表系統(tǒng)狀態(tài)的有限集{1,2,…,n}=〈n〉。

定義1假設有一個和時間不相關的固定概率pi,j，使得：

Prob(Xt+1=i|Xt=j,Xt-1=it-1,…,X0=i0)=Prob(Xt+1=i|Xt=j)=pi,j

其中i,j,i0,i1,…,it-1∈〈n〉。這時，稱{Xt}為馬爾科夫鏈過程。

pi,j表示從當前狀態(tài)j轉(zhuǎn)移到狀態(tài)i的概率，有：

一階馬爾科夫鏈在對未來狀態(tài)預測時，只用到了前一個時刻的信息，忽略了其他時刻的信息，進而影響了預測的準確性。因此，文獻[4]提出一種可以充分利用歷史信息的高階馬爾科夫鏈模型。

1.2 高階馬爾科夫鏈

定義2假設存在一個與時間無關的固定概率pi1,i2,…,im，

0≤pi1,i2,…,im=Prob(Xt+1=i1|Xt=i2,…,Xt-m+2=im)≤1

(1)

文獻[7]將高階馬爾科夫鏈的極限概率分布問題轉(zhuǎn)化為張量-向量積形式：

x=

(2)

初始化:給定轉(zhuǎn)移概率張量,最大執(zhí)行步數(shù)kmax,初始向量x0,滿足x01=1,容忍誤差ε。 (1)令k=1; (2) xk=xm-1k-1; (3)令δ=xk-xk-11,如果δ≤ε,算法停止,輸出xk;否則,令k=k+1,返回步驟2。

文獻[8]在高階冪法的基礎上增加了松弛項，提出松弛高階冪法，主要步驟如下。

初始化:給定轉(zhuǎn)移概率張量,最大執(zhí)行步數(shù)kmax,初始向量x0,滿足x01=1,松弛參數(shù)β,容忍誤差ε。 (1)令k=1; (2) yk=xm-1k-1,^xk=βyk+(1-β)xk-1,xk=proj(^xk); (3)令δ=xk-xk-11,如果δ≤ε,算法停止,輸出xk;否則,令k=k+1,返回步驟2。

2 帶動量項的松弛算法及收斂性分析

高階馬爾科夫鏈的極限概率分布問題在一定程度上可以轉(zhuǎn)化為一個優(yōu)化問題：

(3)

一定程度上f(x)的梯度可以表述為：f(x)=x-xm-1，針對優(yōu)化問題(3)的梯度下降法具有如下迭代格式：

xk=xk-1-αf(xk-1)=xk-1-α(xk-1-

(4)

當α=1時，求解優(yōu)化問題(3)的梯度法即為高階冪法。

動量梯度法[10]是在原始梯度下降法中加入“動量項”，可以充分利用歷史梯度信息來加快梯度法的效率。針對矩陣情形，Xu等[11]提出一種帶動量項的加速冪法。

Nesterov動量加速法[12]是深度學習和人工智能中常用的加速一階算法，其迭代格式為：

Kim等[13]在Nesterov動量梯度法的基礎上增加了一個松弛項γ(yk+1-xk)。受文獻[11-13]的啟發(fā)，本文提出一種帶NAG型動量項的松弛高階冪法HOPM-NAG。

初始化:給定轉(zhuǎn)移概率張量,最大執(zhí)行步數(shù)kmax,動量參數(shù)β,γ,初始向量x0,y0=x0滿足x01=1,容忍誤差ε。令k=1,當k

在HOPM-NAG算法中，步驟1是式(4)中當α=1時的梯度下降法。步驟2中，β(yk-yk-1)稱為動量項,γ(yk+1-xk)稱為松弛項。

文獻[7-8]給出了如下結果。

引理1[7]如果是m階n維的非負張量，存在非零非負向量滿足張量方程(2)，當是不可約張量，向量是正的，如果張量方程(2)將有唯一解，δm定義如下：

其中，S?〈n〉，S′是S在〈n〉中的補集。

引理2[8]如果是m階n維的非負張量，x,y∈Rn，是n維隨機向量，則有：

其中ηm=(1-δm)(m-1)，顯然0≤ηm<1。

基于以上引理，給出HOPM-NAG算法的收斂性分析。

定理1如果是m階n維的非負張量且記是張量方程(2)的解，當時，由HOPM-NAG算法產(chǎn)生的迭代序列{xk}收斂到并且滿足：

當mod(k,3)≠0時，HOPM-NAG算法執(zhí)行高階冪法迭代步，產(chǎn)生的迭代序列{xk}滿足：

(5)

當mod(k,3)=0時，通過HOPM-NAG算法得到：

(6)

式中，β>0，γ>0，將yk=和代入式(6)得：

γ

(7)

由引理2得：

(8)

將式(8)帶入式(7)得：

通過引理3，有：

(9)

當0<κm<1時，由式(5)，式(9)可以進一步證明，對?k≥1有：

3 數(shù)值實驗

實驗在Windows 10操作系統(tǒng)下的MATLAB R2016b中完成，硬件配置為Intel(R)Xeon(R)W-2145CPU@3.70GHz。

使用文獻[7]給出的4個經(jīng)典算例進行數(shù)值實驗。4個算例為：(1)DNA序列{A/G,C,T}；(2)DNA序列{A,C/T,G}；(3)物理學家職業(yè)流動數(shù)據(jù)；(4)風力數(shù)據(jù)。

表1 不同算法在算例中的參數(shù)

從算法的迭代步數(shù)、運行時間及相對殘差3個方面進行比較分析，實驗結果如表2所示。

表2 不同算法求解算例的結果

從表2可以看出，對于同一個概率轉(zhuǎn)移張量，與HOPM算法和R-HOPM算法相比，本文提出的HOPM-NAG算法的運行時間明顯減少。

在容忍誤差相同的情況下，分別采用HOPM算法、R-HOPM算法和HOPM-NAG算法進行求解，不同算法的迭代步數(shù)與相對誤差的關系如圖1所示。

圖1 不同算法的相對殘差與迭代步數(shù)的關系

從圖1可以看出，同一個算例中，在容忍誤差相同的情況下，HOPM-NAG算法求解極限概率分布所用的迭代步數(shù)比其他2種算法更少，說明HOPM-NAG算法的收斂速度最快，算法效果最優(yōu)。

4 結束語

本文提出一種帶NAG型動量項的松弛高階冪法，采用周期間隔的動量外推，通過選取合適的參數(shù)來求解極限概率分布問題。與高階冪法、松弛高階冪法相比，本文所提算法的運行時間更短，迭代步數(shù)更少，顯示了一定的優(yōu)越性。但是，本文算法沒有給出動量參數(shù)和松弛參數(shù)的選取規(guī)則，后續(xù)將針對這個問題展開進一步研究。