周玉科 劉建文 王 妍

1(中國科學(xué)院地理科學(xué)與資源研究所生態(tài)系統(tǒng)網(wǎng)絡(luò)觀測與模擬重點實驗室 北京 100101) 2(福州大學(xué)空間信息工程研究中心數(shù)據(jù)挖掘與信息共享教育部重點實驗室 福建 福州 350116) 3(中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室 北京 100434)

馬爾科夫鏈蒙特卡洛算法并行化設(shè)計與性能分析

周玉科1劉建文2王 妍3

1(中國科學(xué)院地理科學(xué)與資源研究所生態(tài)系統(tǒng)網(wǎng)絡(luò)觀測與模擬重點實驗室 北京 100101)2(福州大學(xué)空間信息工程研究中心數(shù)據(jù)挖掘與信息共享教育部重點實驗室 福建 福州 350116)3(中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室 北京 100434)

馬爾科夫鏈蒙特卡洛MCMC(Markov Chain Monte Carlo)算法廣泛應(yīng)用于地球系統(tǒng)模型中參數(shù)不確定性分析和模擬。由于地球環(huán)境科學(xué)數(shù)據(jù)的高維度、大容量特性，迫切需求高性能的MCMC算法滿足應(yīng)用需求。采用數(shù)據(jù)分治法實現(xiàn)該算法的多核并行化。利用靜態(tài)和動態(tài)分配策略將算法中的多個輸入鏈分配到各CPU；獨立計算并通過共享內(nèi)存實現(xiàn)進程間通信；主進程回收各單元計算結(jié)果，合成最終的馬爾可夫鏈輸出矩陣。采用控制變量法分析不同樣本和馬爾可夫鏈數(shù)量下的算法加速情況。結(jié)果表明在計算規(guī)模較大、動態(tài)負(fù)載均衡的條件下易于獲得較好的加速比，在4個CPU以內(nèi)時效果顯著，之后隨著CPU增加加速效果出現(xiàn)波動或趨于穩(wěn)定。研究表明并行化MCMC能夠利用多核CPU硬件設(shè)施獲得加速效果，更多核數(shù)的加速性能存在進一步優(yōu)化的空間。

馬爾可夫鏈蒙特卡洛算法 分治法則 多核計算 共享內(nèi)存 加速性能

0 引 言

在實際生產(chǎn)和科研應(yīng)用中，經(jīng)典貝葉斯推斷和分析容易受到未知參數(shù)后驗分布為高維、復(fù)雜分布情況的限制，而基于貝葉斯原理的馬爾科夫鏈蒙特卡洛算法MCMC方法能夠綜合先驗和條件信息，并且給出明確的后驗信息，可以針對高維積分進行運算。MCMC已經(jīng)被廣泛地應(yīng)用于統(tǒng)計分析[1-2]、圖像處理[3-4]、信號處理[5-6]、金融分析等領(lǐng)域[7]。在地球環(huán)境科學(xué)領(lǐng)域，MCMC算法主要應(yīng)用于地球系統(tǒng)模型推斷與參數(shù)不確定性分析[8-9]、地質(zhì)化學(xué)過程[10-13]等研究。衛(wèi)星遙感領(lǐng)域，Posselt等利用MCMC算法分析了在劈窗方法下MODIS云數(shù)據(jù)反演的不確定性[14]。MCMC算法收斂性特征為全局性最小，對應(yīng)的線性最小二乘法則容易在局部收斂情況下停滯或者由于非線性問題收斂緩慢，同時大量的應(yīng)用表明 MCMC方法明顯優(yōu)于線性模型[15]。

雖然MCMC算法在理論上容易構(gòu)建，但在實際操作中比較難以實現(xiàn)[11]，尤其是在地球科學(xué)領(lǐng)域數(shù)據(jù)集通常維度較高，同時大氣、海洋、水文等過程模型復(fù)雜，這些因素對MCMC的實現(xiàn)也產(chǎn)生嚴(yán)重挑戰(zhàn)。因此不斷改進MCMC算法的實現(xiàn)模式，增強對地球系統(tǒng)科學(xué)數(shù)據(jù)和模型的適用性是地學(xué)領(lǐng)域MCMC算法應(yīng)用的重要發(fā)展方向。隨著共享內(nèi)存式多核計算機的出現(xiàn)，并行計算的格局也從多機集群模式擴展到了單機并行模式，這也使得利用單個高性能計算機實現(xiàn)復(fù)雜模型下的MCMC算法成為可能。

MCMC算法輸入初始樣本和馬爾科夫鏈數(shù)量越多，獲得的分布函數(shù)越逼近于真實的分布函數(shù)，因此基于大量樣本和馬爾科夫鏈設(shè)計分而治之的MCMC算法，充分利用現(xiàn)代高性能計算機的計算能力，有助于提高大規(guī)?？蒲心Ｐ秃蛙浖男?，并獲得較好的MCMC模擬效果。本文從數(shù)據(jù)并行的角度出發(fā)，根據(jù)MCMC算法的特點以馬爾科夫鏈為數(shù)據(jù)分治對象，將這些對象置于共享內(nèi)存塊中，設(shè)計實現(xiàn)多核CPU架構(gòu)下的并行加速方法。

1 MCMC方法簡介

MCMC是貝葉斯推理的標(biāo)準(zhǔn)抽樣方法，其是利用Markov鏈機制探索狀態(tài)空間以生成樣本的方法。通過產(chǎn)生一個待分析概率密度函數(shù)為穩(wěn)定分布的不可分的各態(tài)歷經(jīng)的馬爾科夫鏈，并以該鏈達(dá)到穩(wěn)態(tài)的樣本作為待分析分布的樣本[16]，其重要特性是無后效性，即事物當(dāng)前狀態(tài)只與上一狀態(tài)有關(guān)，而與以前狀態(tài)無關(guān)。它突破了傳統(tǒng)思維，用相對簡單有效的數(shù)值方法近似求解貝葉斯推斷和分析問題，保持了貝葉斯方法的理論最優(yōu)性和穩(wěn)健性。目前MCMC有很多實現(xiàn)算法，如：Metropolis Hastings算法、Gibbs采樣等。Metroplis-Hastings(M-H)算法是在概率密度歸一化常數(shù)(如后驗概率中的證據(jù)因子)未知條件下的典型MCMC算法，具有適應(yīng)面廣的特點[17-18]。

本文并行化設(shè)計方案應(yīng)用M-H算法(見圖1)，其基本思路為：任意選擇一個不可約的轉(zhuǎn)移概率q(x,y)以及一個轉(zhuǎn)移概率α(x,y)(0α(x,y)1) ，對任一組合(x,x′)，定義：

p(x,x′)=q(x,x′)α(x,x′)x≠x′

如果鏈在時刻t處于狀態(tài)x，即Xt=x，則首先由q(?|x)產(chǎn)生一個潛在的轉(zhuǎn)移x→x′，然后以概率α(x,x′)接受x′作為鏈下一時刻的狀態(tài)值，而以概率1-α(x,x′)拒絕轉(zhuǎn)移到x′，從而鏈在下一時刻仍然處于狀態(tài)x。

圖1 馬爾科夫鏈轉(zhuǎn)移概率示意圖

2 MCMC并行化方法

2.1 并行方案

該MCMC多核并行算法遵循分流-合并(Fork-Join)的并行計算模型(圖2)，以數(shù)據(jù)拆分為基礎(chǔ)，各計算節(jié)點通過共享內(nèi)存實現(xiàn)參數(shù)傳遞。圖2中虛線箭頭表示數(shù)據(jù)到計算節(jié)點的流程示意，并非實際的數(shù)據(jù)分配方案。實線箭頭表示算法過程中真實的運算步驟流程。

MCMC并行化過程首先需要進行并行環(huán)境初始化，包括檢測算法運行平臺的可用CPU核數(shù)，保證輸入多核數(shù)量參數(shù)不大于可用CPU數(shù)量，MCMC標(biāo)準(zhǔn)輸入?yún)?shù)的設(shè)置(概率函數(shù)指針、是否自適應(yīng)模式、接受率等)。計算環(huán)境與算法參數(shù)檢驗通過后進入并行計算過程，輸入的馬爾科夫鏈矩陣以行形式存儲每條鏈，根據(jù)靜態(tài)或者動態(tài)調(diào)度策略，將每行數(shù)據(jù)分發(fā)到各計算單元。實際執(zhí)行的過程是各計算單元通過內(nèi)存標(biāo)記在共享內(nèi)存中訪問對應(yīng)行的馬爾科夫鏈，然后各子進程獨立計算各馬爾科夫鏈。

原始主進程將MCMC概率密度計算表達(dá)式分發(fā)給各子進程，各計算單元獨立計算各自分配的馬爾科夫鏈，最后主進程打開一個管道負(fù)責(zé)回收各計算單元的計算結(jié)果，組合成最終的馬爾可夫鏈矩陣作為MCMC算法輸出。其中馬爾科夫鏈的元素按照分布函數(shù)隨機化生成，因此加載并行程序后還需要初始化隨機數(shù)序列，為確保每個計算單元采用一致的隨機數(shù)，主進程將隨機數(shù)種子發(fā)送到各計算單元，初始化時獨立產(chǎn)生隨機數(shù)。

圖2 MCMC多核并行計算流程圖

2.2 關(guān)鍵技術(shù)

2.2.1 概率密度函數(shù)分發(fā)

本文設(shè)計MCMC算法通過構(gòu)建馬爾科夫鏈從一種banana形式的概率密度分布函數(shù)中進行采樣，采樣中以Log概率密度函數(shù)作為MCMC-MH算法中生成樣本的來源，Log概率密度函數(shù)設(shè)置如下(偽代碼)：

functionlog．pdf(arrayx){ B<-0．04 -x[1]＾2/200-1/2?(x[2]+B?x[1]＾2-100?B)＾2}

圖3為MCMC算法試驗的一個串行運行示例，分別采用控制采樣數(shù)量和自適應(yīng)模式的方法進行驗證。數(shù)據(jù)的目標(biāo)分布采用banana形狀的概率密度，顏色越紅表示分布密度越高，灰色線為等密度線，黑色圈點表示MCMC算法生成樣本。從圖中可以看出，在采樣數(shù)量相同的情況下(圖3(a)和(b)，圖3(c)和(d))，開啟自適應(yīng)模式的MCMC算法生成的結(jié)果分布更接近目標(biāo)分布。比如，無自適應(yīng)模式的生成樣本聚集程度比較高，分布范圍相對較小，在某些高密度和邊緣區(qū)域缺少采樣，而自適應(yīng)模式的采樣分布較均勻，覆蓋樣本區(qū)域范圍較廣。從不同樣本數(shù)量上分析，增加采樣數(shù)量不僅能夠增加采樣密度，而且更加貼近實際分布形狀。根據(jù)串行運行試驗中自適應(yīng)模式的良好效果，以下的并行版本試驗均基于自適應(yīng)模式進行。

圖3 MCMC算法實現(xiàn)示例(不同樣本、不同自適應(yīng)策略)

2.2.2 靜態(tài)與動態(tài)負(fù)載均衡策略

本文MCMC并行方法首先采用靜態(tài)調(diào)度策略，主進程根據(jù)輸入的初始馬爾科夫鏈矩陣行數(shù)和設(shè)定的CPU核數(shù)，將馬爾科夫鏈平均分配到各計算單元進行獨立計算，最后主進程收集計算結(jié)果。

作為對比和并行優(yōu)化，同時實現(xiàn)了MCMC多核并行的動態(tài)負(fù)載均衡策略：如果馬爾科夫鏈的任務(wù)隊列數(shù)量p小于多計算節(jié)點(CPU核心)數(shù)量n，則計算任務(wù)分配個p個計算節(jié)點。如果p大于n，那么前n個任務(wù)按照順序依次分配給前n個節(jié)點。當(dāng)某個計算節(jié)點第一個任務(wù)完成后，該節(jié)點成為空閑節(jié)點，則第n+1個任務(wù)被分配到該節(jié)點。按照此順序一直進行，直到p個任務(wù)執(zhí)行完畢。該策略可以很好地利用集群性能，減少子節(jié)點的空轉(zhuǎn)等待時間，但同時也增加了計算節(jié)點通信的代價，可能降低并行計算性能。另外，為防止增加更多的計算節(jié)點間通信開銷，本文方法的負(fù)載均衡策略沒有設(shè)計將指定的任務(wù)分配到指定計算節(jié)點的方案。

3 性能測試與評價

MCMC多核并行試驗使用c語言和OpenMP并行工具庫作為開發(fā)工具，核心算法采用Metropolis-Hastings算法。實驗硬件平臺為Dell PowerEdge R720機架服務(wù)器，物理8核CPU，主頻為2.1 GHz，內(nèi)存為16 GB，操作系統(tǒng)為Windows Server 2014。實驗首先通過控制MCMC算法中關(guān)鍵參數(shù)的不同取值，分析其對并行效果的影響；然后在不同的并行調(diào)度策略下測試加速效果。

3.1 參數(shù)控制并行試驗

MCMC 方法依賴于模擬的收斂性，但至今仍沒有完全可靠的收斂性診斷方法，這使得收斂性的診斷問題成為MCMC方法實施的難點[4]。為控制不同收斂率帶來的并行效果不確定性，本研究MCMC實驗均采用相同的收斂率(0.23)。

對于MCMC算法的主要參數(shù)配置，并行實驗采用1～8個CPU核數(shù)遞增的測試方案(基于靜態(tài)調(diào)度策略)，對每種CPU核數(shù)配置情況進行1 000次模擬實驗，取1 000次的平均耗時作為該CPU核數(shù)下的并行算法耗時，以降低單次運算帶來的誤差和不確定性。首先在設(shè)定10個馬爾科夫鏈的情況下，采取不同樣本數(shù)量(1 000，2 000，3 000，4 000個樣本)測試并行加速情況，加速時間和加速比結(jié)果如圖4所示。從加速時間曲線層次排列的情況來看，在相同的CPU核數(shù)下樣本數(shù)量越多運行時間約長。在四種樣本數(shù)量下試驗，運行時間均隨著CPU核數(shù)的增加而減少，同時可以發(fā)現(xiàn)運行時間從1核到3核時急劇減少，此后加速時間微弱減少并趨于平緩，說明在較少CPU核數(shù)時加速效果較好，核數(shù)量增多加速效果不明顯。同時從加速比觀察，不同的樣本數(shù)量加速比呈現(xiàn)多樣的波動情況。比如，在3 000和4 000個樣本情況下，加速比在5個核數(shù)之內(nèi)呈現(xiàn)顯著的上升趨勢，之后加速比呈現(xiàn)下降趨勢。1 000和2 000個樣本下加速比出現(xiàn)波動，加速比在3個核數(shù)以后開始下降，而且樣本數(shù)量越少加速比下降趨勢越明顯。此部分試驗說明增加計算任務(wù)量有利于獲得更好的加速比，而且可以避免計算任務(wù)量不足引起的加速波動。

圖4 不同樣本數(shù)量下的多核并行加速情況

同樣在固定樣本數(shù)量(3 000個)情況下測試多個(6,10,14,18)馬爾科夫鏈的并行加速情況，由圖5可以看出加速效果和控制樣本數(shù)量實驗的效果一致。兩次試驗共同說明增加每個計算單元的任務(wù)量有利于持續(xù)獲取加速效果，但是隨著使用CPU核數(shù)的增加，各核訪問共享內(nèi)存中數(shù)據(jù)沖突情況加劇，從而引起加速效果的下降。

圖5 不同馬爾科夫鏈情況下的多核并行加速情況

3.2 動態(tài)負(fù)載均衡的并行實驗

在控制馬爾科夫鏈和樣本數(shù)量兩個變量情況下，基于動態(tài)調(diào)度策略進行實驗獲取兩種情況下(圖6(a)和圖6(b))的加速情況。在馬氏鏈數(shù)量分別設(shè)置為6、10、14、18條的情況下，加速時間曲線與靜態(tài)調(diào)度策略下呈現(xiàn)相似的趨勢，均是在4個CPU核以內(nèi)時間減少比較迅速，之后呈現(xiàn)緩慢的降低趨勢；從加速比對比分析，在4個CPU核以內(nèi)加速比與靜態(tài)策略下基本一致性上升，但是靜態(tài)策略下4個CPU核以上情況下加速比出現(xiàn)波動現(xiàn)象，并且6和10個馬氏鏈情況下加速比較早出現(xiàn)下降。動態(tài)策略下加速比在鏈數(shù)較少和多個CPU核時同樣出現(xiàn)下降(如圖6(b)在10個鏈時)，但是整體呈現(xiàn)上升趨勢。在MCMC樣本數(shù)量設(shè)置為1 000、2 000、3 000、4 000四種情況下(圖7)，動態(tài)策略加速時間曲線在小于4個以內(nèi)CPU核時呈現(xiàn)與靜態(tài)策略時相似的下降趨勢，但是在多于4個核數(shù)時波動情況較小，而且仍然有微弱的加速趨勢；加速比在多于4個核數(shù)時同樣比靜態(tài)調(diào)度下波動較少，只是同樣的有下降趨勢。

圖6 動態(tài)調(diào)度策略下不同馬氏鏈數(shù)加速情況

為定量化地比較分析靜態(tài)調(diào)度與動態(tài)調(diào)度的加速效果，在樣本數(shù)量和馬爾科夫鏈兩個變量控制的基礎(chǔ)上，采用動態(tài)調(diào)度的策略計算運行時間，然后取(Tstatic-Tdynamic)的計算時間差對比分析具體加速效果(圖8)。實驗結(jié)果表明在馬爾科夫鏈數(shù)量控制在3 000條的情況下，使用3個CPU核以內(nèi)時二者的時間差為負(fù)數(shù)，說明靜態(tài)調(diào)度策略優(yōu)于動態(tài)調(diào)度策略，在樣本數(shù)為2 000時出現(xiàn)個別相反情況。在超過3個CPU核時二者時間差為正數(shù)，說明動態(tài)調(diào)度策略優(yōu)于靜態(tài)調(diào)度策略。在兩種參數(shù)控制模式下，3個CPU核內(nèi)的運行時間差整體小于3個以上情況的時間差，說明計算量是影響調(diào)度策略效果的一個重要控制因素。

圖7 動態(tài)調(diào)度策略下不同樣本數(shù)加速情況

3.3 性能與適用性分析

MCMC算法運行過程中并行部分和非并行部分運行時間不可測度，因此上述試驗的運行時間均是算法整體的運行時間。從并行計算的加速效果上分析，MCMC并行實驗效果與理想的線性加速比情況相差較遠(yuǎn)，可以從算法本身特點和并行策略上解釋這種差距。首先，MCMC算法的輸入馬氏鏈內(nèi)部的狀態(tài)矩陣相互依賴性較強，數(shù)據(jù)的耦合性比較高從而降低了各處理器獨立處理數(shù)據(jù)的幾率。后一個矩陣依賴前一個矩陣的狀態(tài)，因此計算過程存在處理器等待上一次計算結(jié)果作為輸入的情況。如果某一個狀態(tài)耗時比較長，會延長整個算法執(zhí)行時間。其次，分布式計算的各條馬氏鏈特征不盡相同，這種差異也會導(dǎo)致各CPU計算耗時不一致從而引起整體計算時間增加。再次，多核CPU的并行策略是基于共享內(nèi)存的方式，雖然內(nèi)存訪問效率比較高，但是在大量的并行循環(huán)過程中也存在內(nèi)存訪問沖突的情況，多個CPU無法及時地獲取輸入數(shù)據(jù)，因此此時并不是并行任務(wù)越多利用效率越高的情況。

4 結(jié) 語

基于數(shù)據(jù)分治并行加速的思路，利用多核CPU和共享內(nèi)存的方法實現(xiàn)了MCMC算法的并行化加速。實驗結(jié)果表明MCMC算法在不同任務(wù)級別下，利用4個以內(nèi)的CPU核心實現(xiàn)并行化的加速效果比較顯著，在4～8個CPU核心上加速效果不明顯而且個別案例呈現(xiàn)加速波動情況。采用變量控制法測試了不同MCMC算法參數(shù)設(shè)置情況下的加速效果，分別控制樣本和馬氏鏈數(shù)量的實驗呈現(xiàn)相似的加速情況，隨著樣本或者馬氏鏈的增加,運行時間由逐漸增加到趨于平穩(wěn)。隨著CPU核數(shù)的增加，加速效果也從逐漸上升到趨于平穩(wěn)或稍微下降。該分析表明增加MCMC計算規(guī)模有利于充分利用多核CPU獲得加速效果，由于MCMC算法中馬爾科夫鏈內(nèi)轉(zhuǎn)移狀態(tài)的計算過程相互依賴性強(不可并行)，多核并行整體效果跟理想的線性加速比存在一定的差距，但是加速情況已經(jīng)明顯優(yōu)于串行MCMC算法效率。

[1] Berthelsen K K, M?ller J. Likelihood and Non-parametric Bayesian MCMC Inference for Spatial Point Processes Based on Perfect Simulation and Path Sampling[J]. Scandinavian Journal of Statistics, 2003, 30(3): 549-564.

[2] Besag J, Green P J. Spatial statistics and Bayesian computation[J]. Journal of the Royal Statistical Society. Series B (Methodological), 1993,55(1): 25-37.

[3] 成敏,張建州, 于巖. 基于多核的MCMC圖像去噪算法并行實現(xiàn)[J]. 計算機工程與應(yīng)用, 2014, 50(18):152-155.

[4] 李成錄,張永仁. 基于貝葉斯與馬爾科夫鏈蒙特卡洛融合算法的圖像處理探究[J]. 電子測試, 2013(16): 014.

[5] 靳曉艷,周希元,張琬琳.多徑衰落信道中基于自適應(yīng) MCMC的調(diào)制識別[J]. 北京郵電大學(xué)學(xué)報, 2014 (1): 31-34.

[6] 李晶,趙擁軍,李冬海.基于馬爾科夫鏈蒙特卡羅的時延估計算法[J]. 物理學(xué)報, 2014, 63(13): 1-7.

[7] 朱新玲. 馬爾科夫鏈蒙特卡羅方法研究綜述[J]. 統(tǒng)計與決策, 2009(21):151-153.

[8] Gallagher K, Charvin K, Nielsen S, et al. Markov chain Monte Carlo (MCMC) sampling methods to determine optimal models, model resolution and model choice for Earth Science problems[J]. Marine and Petroleum Geology, 2009, 26(4): 525-535.

[9] Hassan A E, Bekhit H M, Chapman J B. Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model[J]. Environmental Modelling & Software, 2009, 24(6): 749-763.

[10] Korenaga J. A method to estimate the composition of the bulk silicate Earth in the presence of a hidden geochemical reservoir[J]. Geochimica et Cosmochimica Acta, 2009, 73(22): 6952-6964.

[11] Posselt D J. Markov chain Monte Carlo methods: Theory and applications[M]//Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. II). Springer Berlin Heidelberg, 2013: 59-87.

[12] 張廣智,王丹陽,印興耀.利用 MCMC 方法估算地震參數(shù)[J]. 石油地球物理勘探, 2011, 46(4): 605-609.

[13] 孫海,姚軍,張磊,等.基于孔隙結(jié)構(gòu)的頁巖滲透率計算方法[J]. 中國石油大學(xué)學(xué)報: 自然科學(xué)版, 2014, 38(2): 92-98.

[14] Posselt D J, L'Ecuyer T S, Stephens G L. Exploring the error characteristics of thin ice cloud property retrievals using a Markov chain Monte Carlo algorithm[J]. Journal of Geophysical Research: Atmospheres, 2008, 113(D24).

[15] Posselt D J, Bishop C H. Nonlinear parameter estimation: Comparison of an ensemble Kalman smoother with a Markov chain Monte Carlo algorithm[J]. Monthly Weather Review, 2012, 140(6): 1957-1974.

[16] Andrieu Christophe, De Freitas Nando, Doucet Amaud, et al. An introduction to MCMC for Machine Learning[J]. Machine Learning, 2003, 50: 5-43.

[17] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equations of state calculations by fast computing machines[J]. Journal of Chemical Physics, 1953, 21: 1087-1091.

[18] Hastings W K. Monte Carlo sampling methods using Markov chains and their Applications [J]. Biometrika, 1970, 57: 97-109.

PARALLELDESIGNANDPERFORMANCEANALYSISOFMARKOVCHAINMONTECARLOALGORITHM

Zhou Yuke1Liu Jianwen2Wang Yan3

1(EcologyObservingNetworkandModelingLaboratory,InstituteofGeographicandNatureResourcesResearch,CAS,Beijing100101,China)2(KeyLaboratoryofSpatialDataMiningandInformationSharingofMinistryofEducation,SpatialInformationResearchCenterofFujianProvince,FuzhouUniversity,Fuzhou, 350116,FujianChina)3(StateKeyLaboratoryofSimulationandRegulationofWaterCycleinRiverBasin,ChinaInstituteofWaterResourcesandHydropowerResearch,Beijing100434,China)

Markov chain Monte Carlo (MCMC) algorithm is widely used in the uncertain analysis and simulation of model parameters in the domain of earth system science. Due to the high dimension and high volume characteristic of earth environment data, it is urgent to develop a high performance version of MCMC algorithm. In this paper, we used the data ‘divide and conquer’ strategy to convert the classic MCMC to a multi-core parallel method. Firstly, the input Markov chains were assigned to multiple CPU cores using static or dynamic load balance plan. Secondly, every standalone core computed the Markov chains assigned by the main process, and the communication between the cores was realized through the sharing memory in the computer. Finally, the main process collected the output Markov chains from every CPU core and aggregated to a Markov chains matrix. The speedup of this parallel MCMC method was tested by selecting different sample size and input Markov chains number. The results indicated that the parallel MCMC got better speedup performance when given more computing task to the CPU cores, as well as using a dynamic load balance

plan when assigned the task to every CPU core. This study showed that parallel version of MCMC method could fully use the hardware of modern computer to obtain acceleration. In our future work, the speedup in the case of more than 4 CPU cores would be further optimized.

MCMC Divide and conquer Multi-core computing Sharing memory Speedup

2017-02-23。國家自然科學(xué)基金項目(41601478)；國家重點研發(fā)計劃項目(2016YFC0500103);中科院STS項目(KFJ-SW-STS-167)；資源與環(huán)境信息系統(tǒng)國家重點實際室開放基金項目(2016);中科院地理資源所青年啟動基金。周玉科，助理研究員，主研領(lǐng)域：并行空間分析與生態(tài)遙感。劉建文，碩士生。王妍，助理研究員。

TP3

A

10.3969/j.issn.1000-386x.2017.12.048