基于馬爾可夫鏈改進的原子軌跡預測算法

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

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

關鍵詞：馬爾可夫鏈 馬爾可夫模型 原子運動軌跡 分子動力學

中圖分類號：O211.62 文獻標識碼：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

軌跡分析作為對物體的行為進行學習或描述的基本問題現今已經成為了一個十分熱門的話題。本文設計的算法用于預測原子的運動軌跡，在這一領域，還未有使用馬爾可夫鏈進行預測的研究成果。但在普遍的運動軌跡預測領域，對馬爾可夫鏈的使用是十分廣泛的。

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

本文提出了一個基于齊次馬爾可夫鏈C-K方程，通過回溯糾正的機制來更新馬爾可夫鏈狀態轉移概率矩陣的算法。通過糾正更新來使狀態轉移概率矩陣適應數據源，并趨于穩定。以穩定的狀態轉移概率矩陣來預測數據，盡可能減少分子動力學仿真中的運算量。

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

1 馬爾可夫鏈

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

在一個事件的發展過程中，若事件的每次狀態的轉移都只和前一時刻的狀態有關，和過去的狀態無關，或者說事件的狀態轉移過程是無后效性的，則這樣的狀態轉移過程就稱為馬爾可夫過程。舉例來說，設某個隨機過程的狀態X可取到一個離散集合中的值，這個值隨時間t變化，我們將該值表示為X（t）。此時，時間變量離散或是連續并不影響結果。假設我們將“過去的時間”集合表示為（…，p2，p1），任何“當前時間”n，以及任何“未來時間”t，這些所有時間都在X的取值范圍內，有公式（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，…}為馬爾可夫鏈。在本文中我們討論的是連續時間馬爾可夫鏈，其狀態發生變化的時刻是任意時刻，是連續值。

1.2 齊次馬爾可夫鏈

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

1.3 k步狀態轉移矩陣

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

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

2 算法詳解及思路

2.1 數據預處理與狀態劃分

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

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

原子運動的軌跡由位置決定，而原子運動的位置序列是按時間順序記錄下來的。基于馬爾可夫鏈的預測主要是根據狀態轉移概率矩陣來計算狀態之間直接轉換的概率。

2.2.1 建立類馬爾可夫鏈

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

P = （pij） = （3.1）

2.2.2 確定狀態轉移概率矩陣

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

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

這是一步狀態轉移概率矩陣，根據公式（2.3）可得k步轉移概率矩陣如公式（3.3）所示如下。

P（k）=Pk （3.3）

2.3 預測算法的思路

2.3.1 算法思路

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

2.3.2 主要算法描述及分析

整個算法與普通的馬爾可夫鏈預測算法比較改進的重點在于對狀態轉移概率矩陣的更新。因為我們的目的在于節省計算成本，那么我們可以在較少學習數據的基礎上進行預測，預測3-5步，重新進行數據學習。這樣既可以保證預測的準確度，亦可以保證控制程序運行的時間。

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

3 實驗分析

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

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

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

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

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

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

4 結語

本文提出了一種基于馬爾可夫鏈，改進更新狀態轉移概率矩陣的預測算法，并詳細闡述了預測的方法和步驟。這種方法為得到原子的軌跡提出了一種新的思路，利用歷史原子軌跡基于改進的馬爾可夫鏈進行預測，在學習數據較少的基礎上得到一個可靠的預測結果，為節省計算的開支提供了一種可行方式。

但原子的軌跡的影響因素眾多，在單一的由C-K方程計算來得到下一步的狀態轉移概率矩陣，并由此來計算出下一步的狀態概率的這種方法，對于外界因素的考慮是不夠的。對于狀態的劃分使用簡單的數據預處理也一定程度影響了預測效果，后續的研究可以考慮引入聚類劃分狀態。我們仍需要進一步的研究，改進預測的方法來得到更加可靠的預測結果。

參考文獻

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

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

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

[4] 何麗.基于Markov鏈的用戶瀏覽行為預測方法[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] 夏莉，黃正洪.馬爾可夫鏈在股票價格預測中的應用[J].商業研究.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.