GaBP算法优化与实现

第33卷第2期龙岩学院学报Vol.33 No.22015年4月JOURNAL OF LONGYAN UNIVERSITYApril 2015数学，物理GaBP算法优化与实现郑汉垣(龙岩学院福建龙岩364000)摘要:通过研究经典GaBP算法,实现了同步和异步GaBP算法程序设计和计算实验,并对结果进行了系统的分析。实验表明GaBP优化算法---异 步GaBP算法比经典GaBP算法有更好的计算效率。关键词:大规模稀疏线性方程组;GaBP算法;算法优化中图分类号:TP301.6文献标识码:A文章编号: 1673-4629( 2015 )02-0001-07在对自然科学与社会科学中诸多实际问题进行数值模拟时,最终大都是归结为稀疏线性方程组的求解问题。例如:在结构设计、数值天气预报、油气资源探测、数值风洞、恒星大气分析与核爆模拟等领域,常利用偏微分方程作为数学模型,而在偏微分方程的离散求解中,稀疏线性方程组扮演着十分重要的角色。此外,稀疏线性方程组在数学规划、网络分析、经济分折、离散Markov链等领域中也有着重要的应用。同时,伴随着模拟问题规模的不断增大,稀疏线性方程组求解时间所占的比重也越来越大,在有些应用领域中，其比重已占到了近80%之多。[1-7]正是由于稀疏线性方程组的求解既特别重要,计算又非常耗时,因此众多的科研机构与科学工作者投入了大量的人力和物力来进行研究。GaBP算法是-种针对对称对角占优线性方程组的迭代算法,它是基于递归更新的概率推理算法,具有低计算复杂性和高并行性。正是由于GaBP算法的这两个性质,它很适合用来处理大规模稀疏线性方程组。GaBP 算法不同于经典的迭代算法,也不同于Krylov 子空间算法，GaBP算法对于对称对角占优线性方程组的求解具有良好的收敛性。[8-9]本文首先给出经典GaBP算法的理论，然后类比于经典迭代法,将经典迭代方法的思想引人GaBP算法中,给出了同步和异步GaBP算法和超松驰GaBP算法,同时给出了多种算法的计算效率比较图,试验结果表明了GaBP优化算法一-异步GaBP算法和超松驰GaBP算法同比经典GaBP算法有更好的计算效果。1经典 GaBP算法本节我们回顾文献[8]中关于GaBP算法的具体内容。考虑大规模稀疏线性方程组Ax=b(1)其中,A是已知且非奇异的n阶实数矩阵,b是巳知的n阶向量,x是未知的待求解的n阶向量。收稿日期:2014-05-28中国煤化工作者简介:郑汉垣,男,福建永定人,龙岩学院信息工程学院副教授,主要研究MYHCNMHG理。1GaBP算法优化与实现数学.物理由于本文讨论的是对称稀疏线性方程组,为了方便,后文中出现的系数矩阵A都是对称的。1.1对称线性方程组与概率推理模型首先我们将对称线性方程组与无向图模型对应起来,定义一个无向图{G}=(V,E),其中v是图G中所有顶点组成的集合, E是图G中所有边组成的集合。顶点集V与线性方程组的变量集x= {x.,xo}"相对应;边集E与线性方程组的系数矩阵A中非零元相对应。下面我们给出GaBP算法中需要用到的相关术语和概念。给定系数矩阵A和右端向量b,我们可以给出对应的高斯密度函数:P(x) ~ exp(--x"Ax + b"x)，(2)以及给出与方程组(1)相对应的，由边缘势场( edge potentials) 函数ψv和自势场( self potentials)丽.数φ;构成的无向图G。这里的势函数是由高斯函数式(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等于具有高斯概率密度函数p(x) ~ N(u,A")的图G中的均值向量u A {u,,un,} 。因此根据命题1,求解线性方程组(1)的问题转换为求解图的推理问题，大致的转换过程见图1。Ax =bA.x -b=0min(二x"Ax - b"x)max( --'"Ax +b"x)2max exp( -x"Ax + b"x)图1线性方程组转换为概率推断的过程下面我们就介绍常用于求解概率问题的BP( Belief Propagation) 算法。这里我们先定义一些记号, N(i)为第i个节点对应邻居节点的集合(不含i),N(i)\j为N(i)中除.去节点j。1.2BP算法.BP算法等价于应用Pearl的局部消息传递算法。[10]该方法最初是中国煤化工BP算法最大的优势就是能够利用图的稀疏性。YHCNM HGGaBP算法优化与实现数学.物理BP算法通过图中的边传递一个实值消息，并由“加乘”和“乘”两种计算规则构成。根据1.1节中的式(3) ,传统的加乘规则变成了积分乘规则，"消息m,;(x)表示从节点i沿着相互连接的边传递到节点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算法的特例,因为这里的概率密度分布是高斯分布。下面我们通过1.2节中的连续BP更新规则(式(4)和式(5)取代高斯分布,从而给出GaBP的更新规则。①xφmi| | ψj=ψjixi, φimhi,.-图2节点i 上的消息传递图图2中给出了一个图G中的一部分,主要考虑节点i与相邻节点的传递消息的示意图。从图2中可以看出，每一一个节点上都有-一个自势函数φ; ,每条边上都有- -对(对称的)边缘势函数ψ。消息传递是沿着边上的箭头方向进行的,在传递中只需要计算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算法优化与实现数学.物理|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)对于稠密矩阵而言,无向图G中的一次消息传递需要O(n2)。如果n很大时，传递的消息量是很大的。Bickson等[9将消息量从O(n2 )降低到{O}(n),其主要的思想是:取消节点i向节点j传递消息μ;和P;，而是采用广播的方式，每个节点将相加之后的值(即式(12)和(13))进行广播,P=P+ . Pm，(12)μ;=P7'(Paμ。+ Ewn PuHn),(13)之后根据下面的两个式子，可以重新得到P2y (6)和μvj (7),这样就将传递量从0( n2 )降低到{O}(n)。Pvy=P;-P;,μny=μ; - Pr+;(P: μa)具体详细的GaBP的伪代码可以参看文献[ 8]中的Algorithm 1。2 GaBP算法优化与实现通过对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算法与雅可比迭代算法和高斯-塞德尔迭代算法进行比较后，可知:一般的经典GaBP算法也正如雅可比迭代算法一样是同步更新算法,没有将更新后的消息及时对后继消息更新施加影响。同样的也可以如高斯-塞德尔迭代算法--样对GaBP算法进行优化,让更新后的消息及时地对后继消息的更新施加影响。对应于GaBP算法的“消息累加”和“消息更新”两过程的不同处理方法,也就产生类比雅可比迭代算法和高斯-塞德尔迭代算法的两种不同的实现算法:同步GaBP算法、异步GaBP算法。2.1 同步与异步GaBP算法同步GaBP算法是首先将所有结点进行“消息累加”过程,求得消息聚集和,然后再逐一对结点进行消息传播与更新。其优点是对每一一次的迭代,每-一结点的“消息累加”过程只做- -次 ,计算次数少,一次的迭代步的计算时间会比较小。缺点是没有将更新后的消息及时进中国煤化工息更新施加影响。这样将会带来迭代次数的增加,算法总体的计算时间也就随MHCNM HG4GaBP算法优化与实现数学.物理异步GaBP算法,每次执行“消息传播及更新”处理后,随后也再进行“消息累加”处理,虽说消息累加计算次数增加了,但是这将会使得前面更新后的消息，能立即对后续的消息更新施加影响。也就是“消息累加”处理时,都使用了最新化的消息，缺点是增加了--定量的消息累加计算,却能带来迭代次数的减小，并加快了迭代计算进程,算法的总体计算时间也随之减小了,获得了明显的迭代加速效果。下面同时给出同步和异步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:写= pu/Pr2: 1 end for2: end for28 until convergence: all可convergedOutput: x* = [x]Output: xr = [x]两算法的输入参数是A,b,n,返回解向量x。其中异步算法中说明如下:(1)行号:4- 10是“消息传播及更新”处理过程;(2)行号:11-17是“消息累加”处理过程;中国煤化工(3)行号:20-22是“求解向量"处理过程。MHCNMHG5GaBP算法优化与实现数学.物理3数值试验3.1实验环境本章的实验环境如下:2个Intel Xeon E5-2650的CPU, 96GB内存, Redhat Linux 6. 3 X64操作系统和gcc-4.8.1编译环境。3.2实 验结果及分析我们以美国佛罗里达大学的稀疏矩阵集(UFget)[12]作为实验数据,测试了其中几个稀疏矩阵,给出了两种GaBP算法求解效率比，下面先给出测试矩阵列表(表1):图3是两种GaBP算法迭代次数对比图，其横坐标是测试算法矩阵列表中的矩阵序号，纵坐标是算法求解的迭代次数。图4是两种GaBP算法计算时间对比图,其横坐标是测试算法矩阵列表中的矩阵序号,纵坐标是算法求解的计算时间。可以看出异步GaBP算法比同步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算法计算时间对比4结论本文通过研究经典GaBP算法,实现了同步和异步GaBP算法程序设计和计算实验,通过对两种算法进行数值试验,给出了两种算法的迭代次数与计算时间结果比较表,两结果图表明了异步GaBP算法比同步GaBP算法有更中国煤化空明通过算法优化的异步GaBP算法比经典GaBP算法有更好的计算效率。TY HCNMHG6GaBP算法优化与实现数学.物理参考文献:[1]李开泰,黄艾香.有限元方法及其应用[ M].北京科学出版社,2006.[2]约翰D,安德森.计算流体力学基础及其应用[ M].北京:机械工业出版社, 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].长沙:湖南科学技术出版社,2004.[S]Yousef Saad.系数线性系统的迭代方法[ 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