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