張悅++王新梅1++伍恒

摘 要：馬爾可夫過程是一個具備了馬爾可夫性質(zhì)的隨機過程。具備離散狀態(tài)的馬爾可夫過程，稱為馬爾可夫鏈。馬爾科夫鏈通常被用來建模排隊理論與統(tǒng)計學(xué)中的建模，同時還能夠運用到信號模型的熵編碼技術(shù)中。在文中則主要是將其運用到院子軌跡的預(yù)測算法之中。本文描述的是基于馬爾可夫鏈圍繞預(yù)測化學(xué)原子運動軌跡設(shè)計的算法。這種預(yù)測算法旨在節(jié)省分子動力學(xué)中高昂的運算成本。它通過學(xué)習(xí)回溯并糾正的機制對狀態(tài)轉(zhuǎn)移概率矩陣進行更新，并進行原子運動軌跡的預(yù)測。最后用實驗驗證了該算法的有效性。

關(guān)鍵詞：馬爾可夫鏈 馬爾可夫模型 原子運動軌跡 分子動力學(xué)

中圖分類號：O211.62 文獻(xiàn)標(biāo)識碼：A 文章編號：1674-098X（2016）7（b）-0000-00

Abstract：The Markov process meant a random process with Markov properties. A Markov process with discrete state was called Markov chain. This paper was to forecast the algorithm of chemical atomic trajectory design on the basis of Markov chain. This kind of forecast method aimed to save the expensive operation cost in molecular dynamics. A state transition probability matrix， which tends to be stable， was obtained by learning backtracking and correction mechanism to forecast the atomic trajectory. Finally， the experiment proved the effectiveness of the algorithm.

Key words： markov chain； markov model； atomic trajectory； molecular dynamics

軌跡分析作為對物體的行為進行學(xué)習(xí)或描述的基本問題現(xiàn)今已經(jīng)成為了一個十分熱門的話題。本文設(shè)計的算法用于預(yù)測原子的運動軌跡，在這一領(lǐng)域，還未有使用馬爾可夫鏈進行預(yù)測的研究成果。但在普遍的運動軌跡預(yù)測領(lǐng)域，對馬爾可夫鏈的使用是十分廣泛的。

分子動力學(xué)依靠計算機來模擬分子、原子的運動。它是研究凝聚態(tài)系統(tǒng)的有力工具。該技術(shù)不僅可以得到原子的運動軌跡，還可以觀察到原子運動過程中的能量、速度、加速度等。但是其高昂的計算成本和不斷增長的用戶需求之間的矛盾逾演逾烈。因此分子動力學(xué)仿真的優(yōu)化技術(shù)研究一直是數(shù)學(xué)界、計算化學(xué)界以及計算機界的一個研究熱點[1]。在分子動力學(xué)中，通常，分子、原子的軌跡是通過數(shù)值求解牛頓運動方程得到的。正因為在分子動力學(xué)中假定牛頓運動方程決定原子的運動，這也意味著原子的運動可能會產(chǎn)生特定的軌道，因此我們可以考慮能否利用歷史數(shù)據(jù)對原子運動軌跡進行預(yù)測，來節(jié)省大量運算產(chǎn)生的開支。

本文提出了一個基于齊次馬爾可夫鏈C-K方程，通過回溯糾正的機制來更新馬爾可夫鏈狀態(tài)轉(zhuǎn)移概率矩陣的算法。通過糾正更新來使?fàn)顟B(tài)轉(zhuǎn)移概率矩陣適應(yīng)數(shù)據(jù)源，并趨于穩(wěn)定。以穩(wěn)定的狀態(tài)轉(zhuǎn)移概率矩陣來預(yù)測數(shù)據(jù)，盡可能減少分子動力學(xué)仿真中的運算量。

本文第二章將對馬爾可夫鏈的相關(guān)概念和公式進行介紹；第三章將對整個算法進行詳細(xì)描述和分析；第四章將對算法進行數(shù)據(jù)的測試并且分析；第五章對整個研究進行總結(jié)。

1 馬爾可夫鏈

1.1 馬爾可夫鏈與馬爾可夫過程

在一個事件的發(fā)展過程中，若事件的每次狀態(tài)的轉(zhuǎn)移都只和前一時刻的狀態(tài)有關(guān)，和過去的狀態(tài)無關(guān)，或者說事件的狀態(tài)轉(zhuǎn)移過程是無后效性的，則這樣的狀態(tài)轉(zhuǎn)移過程就稱為馬爾可夫過程。舉例來說，設(shè)某個隨機過程的狀態(tài)X可取到一個離散集合中的值，這個值隨時間t變化，我們將該值表示為X（t）。此時，時間變量離散或是連續(xù)并不影響結(jié)果。假設(shè)我們將“過去的時間”集合表示為（…，p2，p1），任何“當(dāng)前時間”n，以及任何“未來時間”t，這些所有時間都在X的取值范圍內(nèi)，有公式（2.1）所示如下。

… < p2 < p1 < n < t （2.1）

則該過程為馬爾可夫過程。如果對于所有的取值（…，x（p2），x（p1），x（n），x（t）），滿足公式（2.2）所示如下。

P{X（t）=x（t）|X（n）=x（n），X（p1）=x（p1），X（p2）=x（p2），...}

= P{X（t）=x（t）|X（n）=x（n）} （2.2）

則稱{X（n），n=0，1，2，…}為馬爾可夫鏈。在本文中我們討論的是連續(xù)時間馬爾可夫鏈，其狀態(tài)發(fā)生變化的時刻是任意時刻，是連續(xù)值。

1.2 齊次馬爾可夫鏈

馬爾可夫鏈在n時刻的k步轉(zhuǎn)移概率記為，Pij（n，n+k）[2]。若轉(zhuǎn)移概率Pij（n，n+k）只與出發(fā)、到達(dá)狀態(tài)i，j以及轉(zhuǎn)移步數(shù)k有關(guān)，不依賴于時刻n，則將此馬爾可夫鏈稱為齊次馬爾可夫鏈。

1.3 k步狀態(tài)轉(zhuǎn)移矩陣

當(dāng)k為1時，Pij（1）稱為一步轉(zhuǎn)移概率，記為Pij。所有的k步轉(zhuǎn)移概率為Pij（k），其組成的矩陣P（n）就稱為馬爾可夫鏈的k步狀態(tài)轉(zhuǎn)移概率矩陣。根據(jù)C-K方程[3]可得到公式（2.3）如下。

P（k）=PP（k-1）=Pk （2.3）

2 算法詳解及思路

2.1 數(shù)據(jù)預(yù)處理與狀態(tài)劃分

經(jīng)過計算得出的數(shù)據(jù)的精確度是十分高的，但在實際研究過程中我們得到精確度在小數(shù)點后五位的數(shù)據(jù)便可。我們的算法并不是將一個位置作為一個狀態(tài)，而是將兩個位置之間的相對位移作為一個狀態(tài)。因此在我們實驗的過程中，狀態(tài)是有正負(fù)之分的。位移區(qū)間的基數(shù)為0.000005，因此狀態(tài)區(qū)間以0.000005分段，四舍五入。如[0.000005，0.000015]之間的相對位移省略為0.00001。此時的狀態(tài)劃分為1狀態(tài)。在數(shù)據(jù)處理并進行狀態(tài)劃分之后，為了避免給無用的狀態(tài)分配空間，狀態(tài)區(qū)間需要分為異邊、同邊的情況。在這里我們使用一個靈活分配空間的二維數(shù)組bound來儲存狀態(tài)，bound的一維固定長度為2。并且儲存狀態(tài)的0和1數(shù)組其首位值都用于儲存狀態(tài)的個數(shù)。如果狀態(tài)區(qū)間異邊，則0代表的一維數(shù)組儲存正狀態(tài)，1代表的一維數(shù)組儲存負(fù)狀態(tài)；如果狀態(tài)區(qū)間同邊且為正，狀態(tài)（包括零狀態(tài)）都儲存在0代表的一維數(shù)組里，1代表的一維數(shù)組的首位值為0（無負(fù)狀態(tài)）；如果狀態(tài)區(qū)間同邊且為負(fù)，狀態(tài)（包括零狀態(tài)）都儲存在1代表的一維數(shù)組里，0代表的一維數(shù)組首位值為0（無正狀態(tài)）。

2.2 原子運動軌跡模型的建立

原子運動的軌跡由位置決定，而原子運動的位置序列是按時間順序記錄下來的?；隈R爾可夫鏈的預(yù)測主要是根據(jù)狀態(tài)轉(zhuǎn)移概率矩陣來計算狀態(tài)之間直接轉(zhuǎn)換的概率。

2.2.1 建立類馬爾可夫鏈

根據(jù)馬爾可夫鏈的定義，將每個原子運動的過程用一個類馬爾可夫鏈來描述。一個原子的類馬爾可夫鏈預(yù)測模型我們用一個三元組表示。其中，X為離散隨機變量，其取值范圍為{x1，x2，…，xm}。每一個xi對應(yīng)一個原子在兩個時間點之間的位移，亦是單個的狀態(tài)點；P為原子的類馬爾可夫鏈預(yù)測模型的狀態(tài)轉(zhuǎn)移概率矩陣；為原子初始的狀態(tài)概率分布，可表示為（p1，p2，…，pm）[4]。狀態(tài)轉(zhuǎn)移矩陣P公式（3.1）所示如下。

P = （pij） = （3.1）

2.2.2 確定狀態(tài)轉(zhuǎn)移概率矩陣

假設(shè)我們的原子位置之間的位移既馬爾可夫鏈中的狀態(tài)點為m個。狀態(tài)轉(zhuǎn)移概率矩陣P為一個m×m的矩陣。pij表示從i相對位移到j(luò)相對位移的概率。我們通過統(tǒng)計來得到pij。假設(shè)從i相對位移到j(luò)相對位移的次數(shù)用Nij表示。如公式（3.2）所示如下。

Pij= （i≥1，j≤m） （3.2）

這是一步狀態(tài)轉(zhuǎn)移概率矩陣，根據(jù)公式（2.3）可得k步轉(zhuǎn)移概率矩陣如公式（3.3）所示如下。

P（k）=Pk （3.3）

2.3 預(yù)測算法的思路

2.3.1 算法思路

對于原子下一步位置的預(yù)測是基于一個數(shù)據(jù)源的學(xué)習(xí)來進行的。我們的思路是假設(shè)首先取20步數(shù)據(jù)進行學(xué)習(xí)，利用這20步數(shù)據(jù)可以得到一個初始概率以及狀態(tài)轉(zhuǎn)移概率矩陣。當(dāng)預(yù)測第21步的相對位移時，由初始概率以及狀態(tài)轉(zhuǎn)移概率矩陣決定2個可能的相對位移取值，第一次返回的預(yù)測值如果超出容錯范圍，則回溯返回第二個預(yù)測值，如果預(yù)測仍然不準(zhǔn)確，則第21步繼續(xù)學(xué)習(xí)，初始概率以及狀態(tài)轉(zhuǎn)移概率矩陣隨之改變。如果預(yù)測值在容錯范圍以內(nèi)則可將此預(yù)測值加入現(xiàn)有的數(shù)據(jù)源中，更新初始概率以及狀態(tài)轉(zhuǎn)移矩陣。繼續(xù)預(yù)測下一步的值。

2.3.2 主要算法描述及分析

整個算法與普通的馬爾可夫鏈預(yù)測算法比較改進的重點在于對狀態(tài)轉(zhuǎn)移概率矩陣的更新。因為我們的目的在于節(jié)省計算成本，那么我們可以在較少學(xué)習(xí)數(shù)據(jù)的基礎(chǔ)上進行預(yù)測，預(yù)測3-5步，重新進行數(shù)據(jù)學(xué)習(xí)。這樣既可以保證預(yù)測的準(zhǔn)確度，亦可以保證控制程序運行的時間。

在普遍的運用齊次馬爾可夫鏈進行預(yù)測的應(yīng)用中，因為要保證狀態(tài)轉(zhuǎn)移矩陣必須具有一定的穩(wěn)定性。因此必須以足夠多的統(tǒng)計數(shù)據(jù)來保證預(yù)測的精確。這也就決定了馬爾可夫預(yù)測模型必須具有足夠統(tǒng)計數(shù)據(jù)的局限性。我們改進的分子動力學(xué)仿真預(yù)測算法，在基于馬爾可夫鏈C-K方程的基礎(chǔ)上，在數(shù)據(jù)的預(yù)測中，采用糾正并回溯的方法更新狀態(tài)轉(zhuǎn)移矩陣，來縮小所需的數(shù)據(jù)源，節(jié)省需要的計算開支。

3 實驗分析

我們進行實驗來驗證改進算法的有效性。硬件環(huán)境：Intel Core i5-3210M CPU 2.49GHz，內(nèi)存為12GB。軟件環(huán)境：Micros-oft windows 8.1 Professional，Visual Studio 2013。

本文的實驗數(shù)據(jù)是CO2分子使用Hessian-Based Predictor-Correct算法計算得到，計算環(huán)境為NWChem-5.1.1[5]。

實驗數(shù)據(jù)分為三份，三份分別是三次計算的結(jié)果。每一份有300組數(shù)據(jù)，每一組將有9個位置數(shù)據(jù)，分別為碳原子的三維數(shù)據(jù)和兩個氧原子的三維數(shù)據(jù)。

實驗中，其中20步數(shù)據(jù)作為學(xué)習(xí)數(shù)據(jù)。為了保障狀態(tài)轉(zhuǎn)移概率矩陣的有效性，我們每一100步重新學(xué)習(xí)一次。因此每一份300組數(shù)據(jù)都成3次進行實驗。第一次學(xué)習(xí)1-20步，分別預(yù)測至50，80，100步；第二次學(xué)習(xí)101-120步，分別預(yù)測至150，180，200步；第三次學(xué)習(xí)201-220步，分別預(yù)測至250，280，300步。在實驗中，預(yù)測準(zhǔn)確的標(biāo)準(zhǔn)是：與實際數(shù)據(jù)相差小于0.00005，每一步的9維數(shù)據(jù)預(yù)測準(zhǔn)確的個數(shù)大于等于7判定為這一步預(yù)測準(zhǔn)確。常規(guī)的馬爾可夫鏈預(yù)測方法沒有更新這一環(huán)節(jié)，在表中顯示了常規(guī)方法與改進方法之間預(yù)測效果的對比。因三份數(shù)據(jù)的情況基本類似，在此列出一份數(shù)據(jù)的測試結(jié)果。數(shù)據(jù)的預(yù)測情況如表4.1所示如下。

通過實驗可以看出，預(yù)測的次數(shù)在60次540個數(shù)據(jù)時，次數(shù)準(zhǔn)確率與步數(shù)準(zhǔn)確率整體較高。預(yù)測的次數(shù)在80次750個數(shù)據(jù)時，次數(shù)準(zhǔn)確率仍可達(dá)到65%以上，而步數(shù)預(yù)測準(zhǔn)確率可達(dá)到80%以上。可見基于改進的馬爾可夫鏈預(yù)測出的數(shù)據(jù)是可以考慮使用的。

盡管影響原子運動軌跡的因素我們無法考慮完全，并且數(shù)據(jù)的預(yù)測有一定隨機性。但在超過80%的步數(shù)預(yù)測準(zhǔn)確率的基礎(chǔ)上，對計算的數(shù)據(jù)進行一定程度的預(yù)測來節(jié)省高昂的運算成本是可行的。

4 結(jié)語

本文提出了一種基于馬爾可夫鏈，改進更新狀態(tài)轉(zhuǎn)移概率矩陣的預(yù)測算法，并詳細(xì)闡述了預(yù)測的方法和步驟。這種方法為得到原子的軌跡提出了一種新的思路，利用歷史原子軌跡基于改進的馬爾可夫鏈進行預(yù)測，在學(xué)習(xí)數(shù)據(jù)較少的基礎(chǔ)上得到一個可靠的預(yù)測結(jié)果，為節(jié)省計算的開支提供了一種可行方式。

但原子的軌跡的影響因素眾多，在單一的由C-K方程計算來得到下一步的狀態(tài)轉(zhuǎn)移概率矩陣，并由此來計算出下一步的狀態(tài)概率的這種方法，對于外界因素的考慮是不夠的。對于狀態(tài)的劃分使用簡單的數(shù)據(jù)預(yù)處理也一定程度影響了預(yù)測效果，后續(xù)的研究可以考慮引入聚類劃分狀態(tài)。我們?nèi)孕枰M一步的研究，改進預(yù)測的方法來得到更加可靠的預(yù)測結(jié)果。

參考文獻(xiàn)

[1] 伍恒.分子動力學(xué)仿真化研究[D].湖南：中南大學(xué).2015：10-11.

[2] 彭曲，丁治明，郭黎敏.基于馬爾可夫鏈的軌跡預(yù)測[J].計算機科學(xué).2010.37（8）：189-191.

[3] 汪榮鑫.隨機過程[M].西安：西安交通大學(xué)出版社，2006.

[4] 何麗.基于Markov鏈的用戶瀏覽行為預(yù)測方法[J].計算機工程.2008.34（22）：32-33.

[5] Wu H， Lu S， Zhu N， et al. A high order predictor–corrector integration algorithm for first principle chemical dynamics simulations[J].Journal of Theoretical and Computational Chemistry.2016.15（01）：1650003.

[6] 夏莉，黃正洪.馬爾可夫鏈在股票價格預(yù)測中的應(yīng)用[J].商業(yè)研究.2003（10）：62-65.

[7] Gilks W R. Markov chain monte carlo[M]. New York： John Wiley & Sons，2005.

[8] Lourderaj U， Song K， Windus T L， et al. Direct dynamics simulations using Hessian-based predictor-corrector integration algorithms[J]. The Journal of chemical physics.2007.126（4）.

[9] Doob J L.Stochastic Processes[M].New York： John Wiley and Sons，1953.