第33卷第2期龍巖學(xué)院學(xué)報(bào)Vol.33 No.22015年4月JOURNAL OF LONGYAN UNIVERSITYApril 2015數(shù)學(xué)，物理GaBP算法優(yōu)化與實(shí)現(xiàn)鄭漢垣(龍巖學(xué)院福建龍巖364000)摘要:通過(guò)研究經(jīng)典GaBP算法,實(shí)現(xiàn)了同步和異步GaBP算法程序設(shè)計(jì)和計(jì)算實(shí)驗(yàn),并對(duì)結(jié)果進(jìn)行了系統(tǒng)的分析。實(shí)驗(yàn)表明GaBP優(yōu)化算法---異 步GaBP算法比經(jīng)典GaBP算法有更好的計(jì)算效率。關(guān)鍵詞:大規(guī)模稀疏線性方程組;GaBP算法;算法優(yōu)化中圖分類號(hào):TP301.6文獻(xiàn)標(biāo)識(shí)碼:A文章編號(hào): 1673-4629( 2015 )02-0001-07在對(duì)自然科學(xué)與社會(huì)科學(xué)中諸多實(shí)際問(wèn)題進(jìn)行數(shù)值模擬時(shí),最終大都是歸結(jié)為稀疏線性方程組的求解問(wèn)題。例如:在結(jié)構(gòu)設(shè)計(jì)、數(shù)值天氣預(yù)報(bào)、油氣資源探測(cè)、數(shù)值風(fēng)洞、恒星大氣分析與核爆模擬等領(lǐng)域,常利用偏微分方程作為數(shù)學(xué)模型,而在偏微分方程的離散求解中,稀疏線性方程組扮演著十分重要的角色。此外,稀疏線性方程組在數(shù)學(xué)規(guī)劃、網(wǎng)絡(luò)分析、經(jīng)濟(jì)分折、離散Markov鏈等領(lǐng)域中也有著重要的應(yīng)用。同時(shí),伴隨著模擬問(wèn)題規(guī)模的不斷增大,稀疏線性方程組求解時(shí)間所占的比重也越來(lái)越大,在有些應(yīng)用領(lǐng)域中，其比重已占到了近80%之多。[1-7]正是由于稀疏線性方程組的求解既特別重要,計(jì)算又非常耗時(shí),因此眾多的科研機(jī)構(gòu)與科學(xué)工作者投入了大量的人力和物力來(lái)進(jìn)行研究。GaBP算法是-種針對(duì)對(duì)稱對(duì)角占優(yōu)線性方程組的迭代算法,它是基于遞歸更新的概率推理算法,具有低計(jì)算復(fù)雜性和高并行性。正是由于GaBP算法的這兩個(gè)性質(zhì),它很適合用來(lái)處理大規(guī)模稀疏線性方程組。GaBP 算法不同于經(jīng)典的迭代算法,也不同于Krylov 子空間算法，GaBP算法對(duì)于對(duì)稱對(duì)角占優(yōu)線性方程組的求解具有良好的收斂性。[8-9]本文首先給出經(jīng)典GaBP算法的理論，然后類比于經(jīng)典迭代法,將經(jīng)典迭代方法的思想引人GaBP算法中,給出了同步和異步GaBP算法和超松馳GaBP算法,同時(shí)給出了多種算法的計(jì)算效率比較圖,試驗(yàn)結(jié)果表明了GaBP優(yōu)化算法一-異步GaBP算法和超松馳GaBP算法同比經(jīng)典GaBP算法有更好的計(jì)算效果。1經(jīng)典 GaBP算法本節(jié)我們回顧文獻(xiàn)[8]中關(guān)于GaBP算法的具體內(nèi)容。考慮大規(guī)模稀疏線性方程組Ax=b(1)其中,A是已知且非奇異的n階實(shí)數(shù)矩陣,b是巳知的n階向量,x是未知的待求解的n階向量。收稿日期:2014-05-28中國(guó)煤化工作者簡(jiǎn)介:鄭漢垣,男,福建永定人,龍巖學(xué)院信息工程學(xué)院副教授,主要研究MYHCNMHG理。1GaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理由于本文討論的是對(duì)稱稀疏線性方程組,為了方便,后文中出現(xiàn)的系數(shù)矩陣A都是對(duì)稱的。1.1對(duì)稱線性方程組與概率推理模型首先我們將對(duì)稱線性方程組與無(wú)向圖模型對(duì)應(yīng)起來(lái),定義一個(gè)無(wú)向圖{G}=(V,E),其中v是圖G中所有頂點(diǎn)組成的集合, E是圖G中所有邊組成的集合。頂點(diǎn)集V與線性方程組的變量集x= {x.,xo}"相對(duì)應(yīng);邊集E與線性方程組的系數(shù)矩陣A中非零元相對(duì)應(yīng)。下面我們給出GaBP算法中需要用到的相關(guān)術(shù)語(yǔ)和概念。給定系數(shù)矩陣A和右端向量b,我們可以給出對(duì)應(yīng)的高斯密度函數(shù):P(x) ~ exp(--x"Ax + b"x)，(2)以及給出與方程組(1)相對(duì)應(yīng)的，由邊緣勢(shì)場(chǎng)( edge potentials) 函數(shù)ψv和自勢(shì)場(chǎng)( self potentials)麗.數(shù)φ;構(gòu)成的無(wú)向圖G。這里的勢(shì)函數(shù)是由高斯函數(shù)式(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)線性方程組的解向量x° =A'b等于具有高斯概率密度函數(shù)p(x) ~ N(u,A")的圖G中的均值向量u A {u,,un,} 。因此根據(jù)命題1,求解線性方程組(1)的問(wèn)題轉(zhuǎn)換為求解圖的推理問(wèn)題，大致的轉(zhuǎ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線性方程組轉(zhuǎn)換為概率推斷的過(guò)程下面我們就介紹常用于求解概率問(wèn)題的BP( Belief Propagation) 算法。這里我們先定義一些記號(hào), N(i)為第i個(gè)節(jié)點(diǎn)對(duì)應(yīng)鄰居節(jié)點(diǎn)的集合(不含i),N(i)\j為N(i)中除.去節(jié)點(diǎn)j。1.2BP算法.BP算法等價(jià)于應(yīng)用Pearl的局部消息傳遞算法。[10]該方法最初是中國(guó)煤化工BP算法最大的優(yōu)勢(shì)就是能夠利用圖的稀疏性。YHCNM HGGaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理BP算法通過(guò)圖中的邊傳遞一個(gè)實(shí)值消息，并由“加乘”和“乘”兩種計(jì)算規(guī)則構(gòu)成。根據(jù)1.1節(jié)中的式(3) ,傳統(tǒng)的加乘規(guī)則變成了積分乘規(guī)則，"消息m,;(x)表示從節(jié)點(diǎn)i沿著相互連接的邊傳遞到節(jié)點(diǎn)j，即m(x)∞」ψ(x,x)φ,(x)Iew) ,m(x)dx,.(4)邊緣的計(jì)算通常是由乘規(guī)則確定,即p(x)=ad;(x*)Iev( ma(x).(5)這里的x是歸一化常數(shù)。1.3 GaBP算法GaBP算法是連續(xù)的BP算法的特例,因?yàn)檫@里的概率密度分布是高斯分布。下面我們通過(guò)1.2節(jié)中的連續(xù)BP更新規(guī)則(式(4)和式(5)取代高斯分布,從而給出GaBP的更新規(guī)則。①xφmi| | ψj=ψjixi, φimhi,.-圖2節(jié)點(diǎn)i 上的消息傳遞圖圖2中給出了一個(gè)圖G中的一部分,主要考慮節(jié)點(diǎn)i與相鄰節(jié)點(diǎn)的傳遞消息的示意圖。從圖2中可以看出，每一一個(gè)節(jié)點(diǎn)上都有-一個(gè)自勢(shì)函數(shù)φ; ,每條邊上都有- -對(duì)(對(duì)稱的)邊緣勢(shì)函數(shù)ψ。消息傳遞是沿著邊上的箭頭方向進(jìn)行的,在傳遞中只需要計(jì)算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)。根據(jù)式子(3)可知φ(x,) x N(Aa,pi"),ma(x,) ∞N(Mg，pi') ,再根據(jù)引理1,即有Pv=P1+ Es keN()\Pu，(6)和μiy=Pr}(P;μ; +二n Piun),(7)其中P: A A: ,μ; A b;/A.中國(guó)煤化工根據(jù)高斯積分MHCNMHG3GaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理|exp( - ax2 + bx)dx= V(π/a)exp( b2/4a)和式(4),消息m(x;)的分布函數(shù)為P; =- APy，μg=- P;Ag Mury，(9)按照上面相同的步驟處理式(5),則得到邊緣的高斯概率密度函數(shù)N(μ, ,P7I) ,其中P:=P:+ 2 P,( 10)和μ;=pj"(P:μ; + 2 P&Mx).(11)對(duì)于稠密矩陣而言,無(wú)向圖G中的一次消息傳遞需要O(n2)。如果n很大時(shí)，傳遞的消息量是很大的。Bickson等[9將消息量從O(n2 )降低到{O}(n),其主要的思想是:取消節(jié)點(diǎn)i向節(jié)點(diǎn)j傳遞消息μ;和P;，而是采用廣播的方式，每個(gè)節(jié)點(diǎn)將相加之后的值(即式(12)和(13))進(jìn)行廣播,P=P+ . Pm，(12)μ;=P7'(Paμ。+ Ewn PuHn),(13)之后根據(jù)下面的兩個(gè)式子，可以重新得到P2y (6)和μvj (7),這樣就將傳遞量從0( n2 )降低到{O}(n)。Pvy=P;-P;,μny=μ; - Pr+;(P: μa)具體詳細(xì)的GaBP的偽代碼可以參看文獻(xiàn)[ 8]中的Algorithm 1。2 GaBP算法優(yōu)化與實(shí)現(xiàn)通過(guò)對(duì)GaBP計(jì)算過(guò)程進(jìn)行簡(jiǎn)化與整理,可以將GaBP算法寫成如下主要的初始化和迭代兩部分:1、初始化:P;=0,μq=0(i≠j),且有P:AAuμa≌B,2、迭代的三個(gè)主要計(jì)算部分(1)消息累加:P=P。+ 2Puμ=μu+ SwopMu(2)消息傳播及更新:Pvy=P,-P41y=μ.-μ; P;=-APv,uy =- PrVA,uy(3)求解向量:x; =μ;/P:.通過(guò)對(duì)GaBP算法與雅可比迭代算法和高斯-塞德?tīng)柕惴ㄟM(jìn)行比較后，可知:一般的經(jīng)典GaBP算法也正如雅可比迭代算法一樣是同步更新算法,沒(méi)有將更新后的消息及時(shí)對(duì)后繼消息更新施加影響。同樣的也可以如高斯-塞德?tīng)柕惴?-樣對(duì)GaBP算法進(jìn)行優(yōu)化,讓更新后的消息及時(shí)地對(duì)后繼消息的更新施加影響。對(duì)應(yīng)于GaBP算法的“消息累加”和“消息更新”兩過(guò)程的不同處理方法,也就產(chǎn)生類比雅可比迭代算法和高斯-塞德?tīng)柕惴ǖ膬煞N不同的實(shí)現(xiàn)算法:同步GaBP算法、異步GaBP算法。2.1 同步與異步GaBP算法同步GaBP算法是首先將所有結(jié)點(diǎn)進(jìn)行“消息累加”過(guò)程,求得消息聚集和,然后再逐一對(duì)結(jié)點(diǎn)進(jìn)行消息傳播與更新。其優(yōu)點(diǎn)是對(duì)每一一次的迭代,每-一結(jié)點(diǎn)的“消息累加”過(guò)程只做- -次 ,計(jì)算次數(shù)少,一次的迭代步的計(jì)算時(shí)間會(huì)比較小。缺點(diǎn)是沒(méi)有將更新后的消息及時(shí)進(jìn)中國(guó)煤化工息更新施加影響。這樣將會(huì)帶來(lái)迭代次數(shù)的增加,算法總體的計(jì)算時(shí)間也就隨MHCNM HG4GaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理異步GaBP算法,每次執(zhí)行“消息傳播及更新”處理后,隨后也再進(jìn)行“消息累加”處理,雖說(shuō)消息累加計(jì)算次數(shù)增加了,但是這將會(huì)使得前面更新后的消息，能立即對(duì)后續(xù)的消息更新施加影響。也就是“消息累加”處理時(shí),都使用了最新化的消息，缺點(diǎn)是增加了--定量的消息累加計(jì)算,卻能帶來(lái)迭代次數(shù)的減小，并加快了迭代計(jì)算進(jìn)程,算法的總體計(jì)算時(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: ()消息累加(行號(hào):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)消息傳播及更新(行號(hào):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μ=內(nèi)+μμ0l&cend if21: end for22 (3)求解向量(行號(hào):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:寫= pu/Pr2: 1 end for2: end for28 until convergence: all可convergedOutput: x* = [x]Output: xr = [x]兩算法的輸入?yún)?shù)是A,b,n,返回解向量x。其中異步算法中說(shuō)明如下:(1)行號(hào):4- 10是“消息傳播及更新”處理過(guò)程;(2)行號(hào):11-17是“消息累加”處理過(guò)程;中國(guó)煤化工(3)行號(hào):20-22是“求解向量"處理過(guò)程。MHCNMHG5GaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理3數(shù)值試驗(yàn)3.1實(shí)驗(yàn)環(huán)境本章的實(shí)驗(yàn)環(huán)境如下:2個(gè)Intel Xeon E5-2650的CPU, 96GB內(nèi)存, Redhat Linux 6. 3 X64操作系統(tǒng)和gcc-4.8.1編譯環(huán)境。3.2實(shí) 驗(yàn)結(jié)果及分析我們以美國(guó)佛羅里達(dá)大學(xué)的稀疏矩陣集(UFget)[12]作為實(shí)驗(yàn)數(shù)據(jù),測(cè)試了其中幾個(gè)稀疏矩陣,給出了兩種GaBP算法求解效率比，下面先給出測(cè)試矩陣列表(表1):圖3是兩種GaBP算法迭代次數(shù)對(duì)比圖，其橫坐標(biāo)是測(cè)試算法矩陣列表中的矩陣序號(hào)，縱坐標(biāo)是算法求解的迭代次數(shù)。圖4是兩種GaBP算法計(jì)算時(shí)間對(duì)比圖,其橫坐標(biāo)是測(cè)試算法矩陣列表中的矩陣序號(hào),縱坐標(biāo)是算法求解的計(jì)算時(shí)間。可以看出異步GaBP算法比同步GaBP并行有更好的計(jì)算效率。表1測(cè)試算法的 矩陣列表- +同步GaEP算法0十七異步CaBP算法I0 idGroupNameRows Nonzeros872Pothenmeshlem048306十873meshleml30十meshlem6222HBnos667532555875mesh2el2018圖3同步與異步 GaBP算法迭代次數(shù)對(duì)比6 766Nemethnemeth029506 394808|40←同步GaBP算法甘異步GaBP算法nemeth033025.8768nemeth04 9506 3948089769nemeth0595063948081510010 770nemeth06|5011887Norrisv19604 85264123456789101112 888Norris .v2980187025圖4同步與異步 GaBP算法計(jì)算時(shí)間對(duì)比4結(jié)論本文通過(guò)研究經(jīng)典GaBP算法,實(shí)現(xiàn)了同步和異步GaBP算法程序設(shè)計(jì)和計(jì)算實(shí)驗(yàn),通過(guò)對(duì)兩種算法進(jìn)行數(shù)值試驗(yàn),給出了兩種算法的迭代次數(shù)與計(jì)算時(shí)間結(jié)果比較表,兩結(jié)果圖表明了異步GaBP算法比同步GaBP算法有更中國(guó)煤化空明通過(guò)算法優(yōu)化的異步GaBP算法比經(jīng)典GaBP算法有更好的計(jì)算效率。TY HCNMHG6GaBP算法優(yōu)化與實(shí)現(xiàn)數(shù)學(xué).物理參考文獻(xiàn):[1]李開(kāi)泰,黃艾香.有限元方法及其應(yīng)用[ M].北京科學(xué)出版社,2006.[2]約翰D,安德森.計(jì)算流體力學(xué)基礎(chǔ)及其應(yīng)用[ M].北京:機(jī)械工業(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]吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行計(jì)算[ M].長(zhǎng)沙:湖南科學(xué)技術(shù)出版社,2004.[S]Yousef Saad.系數(shù)線性系統(tǒng)的迭代方法[ 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.[責(zé)任編輯:邱維敦]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中國(guó)煤化工MYHCNMHG7