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

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

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

中圖分類(lèi)號(hào)：O211.62 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：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

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

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

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

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

1 馬爾可夫鏈

1.1 馬爾可夫鏈與馬爾可夫過(guò)程

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

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

則該過(guò)程為馬爾可夫過(guò)程。如果對(duì)于所有的取值（…，x（p2），x（p1），x（n），x（t）），滿(mǎn)足公式（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）

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

1.2 齊次馬爾可夫鏈

馬爾可夫鏈在n時(shí)刻的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)，不依賴(lài)于時(shí)刻n，則將此馬爾可夫鏈稱(chēng)為齊次馬爾可夫鏈。

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

當(dāng)k為1時(shí)，Pij（1）稱(chēng)為一步轉(zhuǎn)移概率，記為Pij。所有的k步轉(zhuǎn)移概率為Pij（k），其組成的矩陣P（n）就稱(chēng)為馬爾可夫鏈的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)過(guò)計(jì)算得出的數(shù)據(jù)的精確度是十分高的，但在實(shí)際研究過(guò)程中我們得到精確度在小數(shù)點(diǎn)后五位的數(shù)據(jù)便可。我們的算法并不是將一個(gè)位置作為一個(gè)狀態(tài)，而是將兩個(gè)位置之間的相對(duì)位移作為一個(gè)狀態(tài)。因此在我們實(shí)驗(yàn)的過(guò)程中，狀態(tài)是有正負(fù)之分的。位移區(qū)間的基數(shù)為0.000005，因此狀態(tài)區(qū)間以0.000005分段，四舍五入。如[0.000005，0.000015]之間的相對(duì)位移省略為0.00001。此時(shí)的狀態(tài)劃分為1狀態(tài)。在數(shù)據(jù)處理并進(jìn)行狀態(tài)劃分之后，為了避免給無(wú)用的狀態(tài)分配空間，狀態(tài)區(qū)間需要分為異邊、同邊的情況。在這里我們使用一個(gè)靈活分配空間的二維數(shù)組bound來(lái)儲(chǔ)存狀態(tài)，bound的一維固定長(zhǎng)度為2。并且儲(chǔ)存狀態(tài)的0和1數(shù)組其首位值都用于儲(chǔ)存狀態(tài)的個(gè)數(shù)。如果狀態(tài)區(qū)間異邊，則0代表的一維數(shù)組儲(chǔ)存正狀態(tài)，1代表的一維數(shù)組儲(chǔ)存負(fù)狀態(tài)；如果狀態(tài)區(qū)間同邊且為正，狀態(tài)（包括零狀態(tài)）都儲(chǔ)存在0代表的一維數(shù)組里，1代表的一維數(shù)組的首位值為0（無(wú)負(fù)狀態(tài)）；如果狀態(tài)區(qū)間同邊且為負(fù)，狀態(tài)（包括零狀態(tài)）都儲(chǔ)存在1代表的一維數(shù)組里，0代表的一維數(shù)組首位值為0（無(wú)正狀態(tài)）。

2.2 原子運(yùn)動(dòng)軌跡模型的建立

原子運(yùn)動(dòng)的軌跡由位置決定，而原子運(yùn)動(dòng)的位置序列是按時(shí)間順序記錄下來(lái)的。基于馬爾可夫鏈的預(yù)測(cè)主要是根據(jù)狀態(tài)轉(zhuǎn)移概率矩陣來(lái)計(jì)算狀態(tài)之間直接轉(zhuǎn)換的概率。

2.2.1 建立類(lèi)馬爾可夫鏈

根據(jù)馬爾可夫鏈的定義，將每個(gè)原子運(yùn)動(dòng)的過(guò)程用一個(gè)類(lèi)馬爾可夫鏈來(lái)描述。一個(gè)原子的類(lèi)馬爾可夫鏈預(yù)測(cè)模型我們用一個(gè)三元組表示。其中，X為離散隨機(jī)變量，其取值范圍為{x1，x2，…，xm}。每一個(gè)xi對(duì)應(yīng)一個(gè)原子在兩個(gè)時(shí)間點(diǎn)之間的位移，亦是單個(gè)的狀態(tài)點(diǎn)；P為原子的類(lèi)馬爾可夫鏈預(yù)測(cè)模型的狀態(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)點(diǎn)為m個(gè)。狀態(tài)轉(zhuǎn)移概率矩陣P為一個(gè)m×m的矩陣。pij表示從i相對(duì)位移到j(luò)相對(duì)位移的概率。我們通過(guò)統(tǒng)計(jì)來(lái)得到pij。假設(shè)從i相對(duì)位移到j(luò)相對(duì)位移的次數(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ù)測(cè)算法的思路

2.3.1 算法思路

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

2.3.2 主要算法描述及分析

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

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

3 實(shí)驗(yàn)分析

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

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

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

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

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

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

4 結(jié)語(yǔ)

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

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

參考文獻(xiàn)

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

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

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

[4] 何麗.基于Markov鏈的用戶(hù)瀏覽行為預(yù)測(cè)方法[J].計(jì)算機(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] 夏莉，黃正洪.馬爾可夫鏈在股票價(jià)格預(yù)測(cè)中的應(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.