GaBP算法優(yōu)化與實(shí)現

第33卷第2期龍巖學(xué)院學(xué)報Vol.33 No.22015年4月JOURNAL OF LONGYAN UNIVERSITYApril 2015數學(xué)，物理GaBP算法優(yōu)化與實(shí)現鄭漢垣(龍巖學(xué)院福建龍巖364000)摘要:通過(guò)研究經(jīng)典GaBP算法,實(shí)現了同步和異步GaBP算法程序設計和計算實(shí)驗,并對結果進(jìn)行了系統的分析。實(shí)驗表明GaBP優(yōu)化算法---異 步GaBP算法比經(jīng)典GaBP算法有更好的計算效率。關(guān)鍵詞:大規模稀疏線(xiàn)性方程組;GaBP算法;算法優(yōu)化中圖分類(lèi)號:TP301.6文獻標識碼:A文章編號: 1673-4629( 2015 )02-0001-07在對自然科學(xué)與社會(huì )科學(xué)中諸多實(shí)際問(wèn)題進(jìn)行數值模擬時(shí),最終大都是歸結為稀疏線(xiàn)性方程組的求解問(wèn)題。例如:在結構設計、數值天氣預報、油氣資源探測、數值風(fēng)洞、恒星大氣分析與核爆模擬等領(lǐng)域,常利用偏微分方程作為數學(xué)模型,而在偏微分方程的離散求解中,稀疏線(xiàn)性方程組扮演著(zhù)十分重要的角色。此外,稀疏線(xiàn)性方程組在數學(xué)規劃、網(wǎng)絡(luò )分析、經(jīng)濟分折、離散Markov鏈等領(lǐng)域中也有著(zhù)重要的應用。同時(shí),伴隨著(zhù)模擬問(wèn)題規模的不斷增大,稀疏線(xiàn)性方程組求解時(shí)間所占的比重也越來(lái)越大,在有些應用領(lǐng)域中，其比重已占到了近80%之多。[1-7]正是由于稀疏線(xiàn)性方程組的求解既特別重要,計算又非常耗時(shí),因此眾多的科研機構與科學(xué)工作者投入了大量的人力和物力來(lái)進(jìn)行研究。GaBP算法是-種針對對稱(chēng)對角占優(yōu)線(xiàn)性方程組的迭代算法,它是基于遞歸更新的概率推理算法,具有低計算復雜性和高并行性。正是由于GaBP算法的這兩個(gè)性質(zhì),它很適合用來(lái)處理大規模稀疏線(xiàn)性方程組。GaBP 算法不同于經(jīng)典的迭代算法,也不同于Krylov 子空間算法，GaBP算法對于對稱(chēng)對角占優(yōu)線(xiàn)性方程組的求解具有良好的收斂性。[8-9]本文首先給出經(jīng)典GaBP算法的理論，然后類(lèi)比于經(jīng)典迭代法,將經(jīng)典迭代方法的思想引人GaBP算法中,給出了同步和異步GaBP算法和超松馳GaBP算法,同時(shí)給出了多種算法的計算效率比較圖,試驗結果表明了GaBP優(yōu)化算法一-異步GaBP算法和超松馳GaBP算法同比經(jīng)典GaBP算法有更好的計算效果。1經(jīng)典 GaBP算法本節我們回顧文獻[8]中關(guān)于GaBP算法的具體內容?？紤]大規模稀疏線(xiàn)性方程組Ax=b(1)其中,A是已知且非奇異的n階實(shí)數矩陣,b是巳知的n階向量,x是未知的待求解的n階向量。收稿日期:2014-05-28中國煤化工作者簡(jiǎn)介:鄭漢垣,男,福建永定人,龍巖學(xué)院信息工程學(xué)院副教授,主要研究MYHCNMHG理。1GaBP算法優(yōu)化與實(shí)現數學(xué).物理由于本文討論的是對稱(chēng)稀疏線(xiàn)性方程組,為了方便,后文中出現的系數矩陣A都是對稱(chēng)的。1.1對稱(chēng)線(xiàn)性方程組與概率推理模型首先我們將對稱(chēng)線(xiàn)性方程組與無(wú)向圖模型對應起來(lái),定義一個(gè)無(wú)向圖{G}=(V,E),其中v是圖G中所有頂點(diǎn)組成的集合, E是圖G中所有邊組成的集合。頂點(diǎn)集V與線(xiàn)性方程組的變量集x= {x.,xo}"相對應;邊集E與線(xiàn)性方程組的系數矩陣A中非零元相對應。下面我們給出GaBP算法中需要用到的相關(guān)術(shù)語(yǔ)和概念。給定系數矩陣A和右端向量b,我們可以給出對應的高斯密度函數:P(x) ~ exp(--x"Ax + b"x)，(2)以及給出與方程組(1)相對應的，由邊緣勢場(chǎng)( edge potentials) 函數ψv和自勢場(chǎng)( self potentials)麗.數φ;構成的無(wú)向圖G。這里的勢函數是由高斯函數式(2)分解得到:P(x) x I,"=1中(x,)II ψ(x,x). .(3)事實(shí)上:ψ(x,x) A exp(-一xAx),φ(x) A exp(-一Ax? + b.x).命題1 ([8]中 的Proposition 1)線(xiàn)性方程組的解向量x° =A'b等于具有高斯概率密度函數p(x) ~ N(u,A")的圖G中的均值向量u A {u,,un,} 。因此根據命題1,求解線(xiàn)性方程組(1)的問(wèn)題轉換為求解圖的推理問(wèn)題，大致的轉換過(guò)程見(jiàn)圖1。Ax =bA.x -b=0min(二x"Ax - b"x)max( --'"Ax +b"x)2max exp( -x"Ax + b"x)圖1線(xiàn)性方程組轉換為概率推斷的過(guò)程下面我們就介紹常用于求解概率問(wèn)題的BP( Belief Propagation) 算法。這里我們先定義一些記號, N(i)為第i個(gè)節點(diǎn)對應鄰居節點(diǎn)的集合(不含i),N(i)\j為N(i)中除.去節點(diǎn)j。1.2BP算法.BP算法等價(jià)于應用Pearl的局部消息傳遞算法。[10]該方法最初是中國煤化工BP算法最大的優(yōu)勢就是能夠利用圖的稀疏性。YHCNM HGGaBP算法優(yōu)化與實(shí)現數學(xué).物理BP算法通過(guò)圖中的邊傳遞一個(gè)實(shí)值消息，并由“加乘”和“乘”兩種計算規則構成。根據1.1節中的式(3) ,傳統的加乘規則變成了積分乘規則，"消息m,;(x)表示從節點(diǎn)i沿著(zhù)相互連接的邊傳遞到節點(diǎn)j，即m(x)∞」ψ(x,x)φ,(x)Iew) ,m(x)dx,.(4)邊緣的計算通常是由乘規則確定,即p(x)=ad;(x*)Iev( ma(x).(5)這里的x是歸一化常數。1.3 GaBP算法GaBP算法是連續的BP算法的特例,因為這里的概率密度分布是高斯分布。下面我們通過(guò)1.2節中的連續BP更新規則(式(4)和式(5)取代高斯分布,從而給出GaBP的更新規則。①xφmi| | ψj=ψjixi, φimhi,.-圖2節點(diǎn)i 上的消息傳遞圖圖2中給出了一個(gè)圖G中的一部分,主要考慮節點(diǎn)i與相鄰節點(diǎn)的傳遞消息的示意圖。從圖2中可以看出，每一一個(gè)節點(diǎn)上都有-一個(gè)自勢函數φ; ,每條邊上都有- -對(對稱(chēng)的)邊緣勢函數ψ。消息傳遞是沿著(zhù)邊上的箭頭方向進(jìn)行的,在傳遞中只需要計算m;。下面考慮式(4)的右端，首先引人下面的引理1。引理1 ([8] Lemma2) 令f(x)和f2(x) 都有高斯概率分布，即f(x) ~ N(μ,pil) 和f2(x) ~N(μr,p2l) ,則f(x)=f(x)f2(x) ~ N(μ,p-') ,這里μ=ρ '(p1M + P21H2)。令φ(x)IIew m2(x) ~ N(uv,P)。根據式子(3)可知φ(x,) x N(Aa,pi"),ma(x,) ∞N(Mg，pi') ,再根據引理1,即有Pv=P1+ Es keN()\Pu，(6)和μiy=Pr}(P;μ; +二n Piun),(7)其中P: A A: ,μ; A b;/A.中國煤化工根據高斯積分MHCNMHG3GaBP算法優(yōu)化與實(shí)現數學(xué).物理|exp( - ax2 + bx)dx= V(π/a)exp( b2/4a)和式(4),消息m(x;)的分布函數為P; =- APy，μg=- P;Ag Mury，(9)按照上面相同的步驟處理式(5),則得到邊緣的高斯概率密度函數N(μ, ,P7I) ,其中P:=P:+ 2 P,( 10)和μ;=pj"(P:μ; + 2 P&Mx).(11)對于稠密矩陣而言,無(wú)向圖G中的一次消息傳遞需要O(n2)。如果n很大時(shí)，傳遞的消息量是很大的。Bickson等[9將消息量從O(n2 )降低到{O}(n),其主要的思想是:取消節點(diǎn)i向節點(diǎn)j傳遞消息μ;和P;，而是采用廣播的方式，每個(gè)節點(diǎn)將相加之后的值(即式(12)和(13))進(jìn)行廣播,P=P+ . Pm，(12)μ;=P7'(Paμ。+ Ewn PuHn),(13)之后根據下面的兩個(gè)式子，可以重新得到P2y (6)和μvj (7),這樣就將傳遞量從0( n2 )降低到{O}(n)。Pvy=P;-P;,μny=μ; - Pr+;(P: μa)具體詳細的GaBP的偽代碼可以參看文獻[ 8]中的Algorithm 1。2 GaBP算法優(yōu)化與實(shí)現通過(guò)對GaBP計算過(guò)程進(jìn)行簡(jiǎn)化與整理,可以將GaBP算法寫(xiě)成如下主要的初始化和迭代兩部分:1、初始化:P;=0,μq=0(i≠j),且有P:AAuμa≌B,2、迭代的三個(gè)主要計算部分(1)消息累加:P=P。+ 2Puμ=μu+ SwopMu(2)消息傳播及更新:Pvy=P,-P41y=μ.-μ; P;=-APv,uy =- PrVA,uy(3)求解向量:x; =μ;/P:.通過(guò)對GaBP算法與雅可比迭代算法和高斯-塞德?tīng)柕惴ㄟM(jìn)行比較后，可知:一般的經(jīng)典GaBP算法也正如雅可比迭代算法一樣是同步更新算法,沒(méi)有將更新后的消息及時(shí)對后繼消息更新施加影響。同樣的也可以如高斯-塞德?tīng)柕惴?-樣對GaBP算法進(jìn)行優(yōu)化,讓更新后的消息及時(shí)地對后繼消息的更新施加影響。對應于GaBP算法的“消息累加”和“消息更新”兩過(guò)程的不同處理方法,也就產(chǎn)生類(lèi)比雅可比迭代算法和高斯-塞德?tīng)柕惴ǖ膬煞N不同的實(shí)現算法:同步GaBP算法、異步GaBP算法。2.1 同步與異步GaBP算法同步GaBP算法是首先將所有結點(diǎn)進(jìn)行“消息累加”過(guò)程,求得消息聚集和,然后再逐一對結點(diǎn)進(jìn)行消息傳播與更新。其優(yōu)點(diǎn)是對每一一次的迭代,每-一結點(diǎn)的“消息累加”過(guò)程只做- -次 ,計算次數少,一次的迭代步的計算時(shí)間會(huì )比較小。缺點(diǎn)是沒(méi)有將更新后的消息及時(shí)進(jìn)中國煤化工息更新施加影響。這樣將會(huì )帶來(lái)迭代次數的增加,算法總體的計算時(shí)間也就隨MHCNM HG4GaBP算法優(yōu)化與實(shí)現數學(xué).物理異步GaBP算法,每次執行“消息傳播及更新”處理后,隨后也再進(jìn)行“消息累加”處理,雖說(shuō)消息累加計算次數增加了,但是這將會(huì )使得前面更新后的消息，能立即對后續的消息更新施加影響。也就是“消息累加”處理時(shí),都使用了最新化的消息，缺點(diǎn)是增加了--定量的消息累加計算,卻能帶來(lái)迭代次數的減小，并加快了迭代計算進(jìn)程,算法的總體計算時(shí)間也隨之減小了,獲得了明顯的迭代加速效果。下面同時(shí)給出同步和異步GaBP算法,其描述如下:(算法1和算法2)算法1同步GaBP算法.算法2異步 GaBP算法Algorithm 1同步GaBP算法Algorithm 2異步GaBP算法Initialization:1: Piy=0; =0;1: Piy=0; y=0;lteration:上P:=Ag; 山=b;2: repeatx: ()消息累加(行號:4--12)Iteration:for llie 0.ndo3: repeatP:=As: μ=br;4: for alli∈0.n do& forallje0.ndofor allj∈0.n doif(i≠jAnd Anj≠0) thenif(i≠jAnd Ay≠0) thenP=P:+Pp山=山+HpPy=-AE/(P,-P)8Huy= -Ag(4i-up/(P.- Py)end forend it讓endfou13: (2)消息傳播及更新(行號:14--21)forallie 0.ndo1:P,=Ag; μ=b;: :for all jc0.n do2for all je0.ndoit(i≠jAnd Ayj≠)thenPj=-Aj/(P-Pμ)4:P=P+Pp阿= -Ajy(4-pa)(P.-Pp)s:endifμ=內+μμ0l&cend if21: end for22 (3)求解向量(行號:23--27)18: end for2: forallie0.n doP:=Au+ EuempPu19 until convergence: all μ4J converged25:μ=bn+ ZevwnPu20 for alli∈0.n do ;26有=μ/P:21:寫(xiě)= pu/Pr2: 1 end for2: end for28 until convergence: all可convergedOutput: x* = [x]Output: xr = [x]兩算法的輸入參數是A,b,n,返回解向量x。其中異步算法中說(shuō)明如下:(1)行號:4- 10是“消息傳播及更新”處理過(guò)程;(2)行號:11-17是“消息累加”處理過(guò)程;中國煤化工(3)行號:20-22是“求解向量"處理過(guò)程。MHCNMHG5GaBP算法優(yōu)化與實(shí)現數學(xué).物理3數值試驗3.1實(shí)驗環(huán)境本章的實(shí)驗環(huán)境如下:2個(gè)Intel Xeon E5-2650的CPU, 96GB內存, Redhat Linux 6. 3 X64操作系統和gcc-4.8.1編譯環(huán)境。3.2實(shí) 驗結果及分析我們以美國佛羅里達大學(xué)的稀疏矩陣集(UFget)[12]作為實(shí)驗數據,測試了其中幾個(gè)稀疏矩陣,給出了兩種GaBP算法求解效率比，下面先給出測試矩陣列表(表1):圖3是兩種GaBP算法迭代次數對比圖，其橫坐標是測試算法矩陣列表中的矩陣序號，縱坐標是算法求解的迭代次數。圖4是兩種GaBP算法計算時(shí)間對比圖,其橫坐標是測試算法矩陣列表中的矩陣序號,縱坐標是算法求解的計算時(shí)間?？梢钥闯霎惒紾aBP算法比同步GaBP并行有更好的計算效率。表1測試算法的 矩陣列表- +同步GaEP算法0十七異步CaBP算法I0 idGroupNameRows Nonzeros872Pothenmeshlem048306十873meshleml30十meshlem6222HBnos667532555875mesh2el2018圖3同步與異步 GaBP算法迭代次數對比6 766Nemethnemeth029506 394808|40←同步GaBP算法甘異步GaBP算法nemeth033025.8768nemeth04 9506 3948089769nemeth0595063948081510010 770nemeth06|5011887Norrisv19604 85264123456789101112 888Norris .v2980187025圖4同步與異步 GaBP算法計算時(shí)間對比4結論本文通過(guò)研究經(jīng)典GaBP算法,實(shí)現了同步和異步GaBP算法程序設計和計算實(shí)驗,通過(guò)對兩種算法進(jìn)行數值試驗,給出了兩種算法的迭代次數與計算時(shí)間結果比較表,兩結果圖表明了異步GaBP算法比同步GaBP算法有更中國煤化空明通過(guò)算法優(yōu)化的異步GaBP算法比經(jīng)典GaBP算法有更好的計算效率。TY HCNMHG6GaBP算法優(yōu)化與實(shí)現數學(xué).物理參考文獻:[1]李開(kāi)泰,黃艾香.有限元方法及其應用[ M].北京科學(xué)出版社,2006.[2]約翰D,安德森.計算流體力學(xué)基礎及其應用[ M].北京:機械工業(yè)出版社, 2007.[3] Em Karmiadakis G, Sherwin s. Spectral/hp Element Methods for Computational Fluid Dynamics [ M]. Oxford Univ.Press，. Oxford, etc., 2nd edt. 2005.[4]吳建平,王正華,李曉梅.稀疏線(xiàn)性方程組的高效求解與并行計算[ M].長(cháng)沙:湖南科學(xué)技術(shù)出版社,2004.[S]Yousef Saad.系數線(xiàn)性系統的迭代方法[ M].2版.北京:科學(xué)出版社, 2009.[6]Sun X H, Zhang W.A parallel two-level hybrid method for tridiagonal systems and its application to fast poisson solvers[J].IEEE Transactions on Parallel and Distributed Systems , 2004, 15:97- 106.[7]Thomas J. Ashby, Pieter Ghysels, W imHeirman, WimV anroose. The Impact of Global Communication Latency at ExtremeScales on Krylov Methods[ C].12th Intermational Conference on Algorithms and Architerctures for Parallel Processing( ICA3PP -12) , Fukuoka , Japan, 2012: 428-442.[8]Ori Shental, Paul H. Siegel, Jack K. Wolf, Danny Bickson, Dany Dolev: Gaussian belief propagation solver for systemsof linear equations[ C]. ISIT 2008: 1863- 1867.[9]Yousef EI-Kurdi, Warren J Gross, and Dennis Giannacopoulos. Eficient implementation of Gaussian Belief PropagationSolver for Large Sparse Diagonally Dominant linear systems[J]. IEEE Transactions on magnetics, 2012, 48(2) :471-474.[ 10]Pearl J Probabilistic Reasoning in Ielligent Systems: Networks of Plausible Inference [ M]. San Francisco: MorganKaufmann，1988.[11] Weiss Y, Freeman W T. Correctness of belief propagation in Gaussian graphical models of arbitrarytopology[ J]. NeuralComputation, 2001, 13( 10): 2173-2200.[12] DAVIS A，Hu Y. The University of Florida Sparse Matrix Collection [ J]. ACM Transactions on MathematicalSoftware, 2011, 38(1): 1-28.[責任編輯:邱維敦]Optimization and Implementation of GaBP AlgorithmZHENG HanyuanAbstract: By studying the classical GaBP algorithm， the synchronous and asynchronous implementations using variousprogramming languages are given and the test results are systematically analyzed. The results ilustrate that the iterative acceleration;aBP algorithm，asynchronous GaBP algorithm, is faster than the classical counterpart in convergence speed.Key words: large-scale sparse linear systems; GaBP algorithm; algorithm optimization中國煤化工MYHCNMHG7