基于譜方法的無向賦權圖剖分算法
2009-01-01 00:00:00冷明孫凌宇郁松年
計算機應用研究 2009年6期
摘 要：在多水平方法初始剖分階段提出了一種基于譜方法的無向賦權圖剖分算法SPWUG，給出了基于Lanczos迭代計算Laplacian矩陣次小特征值及特征向量的實現細節。SPWUG算法借助Laplacian矩陣次小特征值對應的特征向量，刻畫了節點間相對距離，將基于非賦權無向圖的Laplacian譜理論在圖的剖分應用方面擴展到無向賦權圖上，實現了對最小圖的初始剖分。基于ISPD98電路測試基準的實驗表明，SPWUG算法取得了一定性能的改進。實驗分析反映了在多水平方法中，最小圖上的全局近似最優剖分可能是初始圖的局部最優剖分，需要加強優化階段的遷移優化算法逃離局部最優的能力。
關鍵詞：多水平方法; 剖分； 無向賦權圖； 譜方法
中圖分類號：TP391文獻標志碼：A
文章編號：1001-3695(2009)06-2086-04
doi:10.3969/j.issn.1001-3695.2009.06.025
Spectral-based weighted undirected graph partitioning algorithm
LENG Ming1，2， SUN Ling-yu1， YU Song-nian2
(1. Dept. of Computer Science， Jinggangshan University， Ji’an Jiangxi 343009， China；2.School of Computer Engineering Science，Shanghai University，Shanghai 200072， China)
Abstract:During the initial partitioning phase of multilevel method， this paper proposed the algorithm of spectral partitioning for weighted undirected graph (SPWUG) and gave the Fiedler vector calculation of Laplacian matrix using the Lanczos iteration method. The components of the Fiedler vector were weighted on the corresponding vertices of graph such that differences of the components provide information about the distances between the vertices. According to the distances， the SPWUG algorithm computed the initial partition of the coarsest graph by extending the Laplacian spectral graph theory from the unweighted undirected graph into the weighted undirected graph. The experimental evaluations show that the SPWUG algorithm produces the performance improvement. Meanwhile， the analysis shows that the approximate global minimum partition of the coarsest graph may be the local minimum partition of the finest graph and need to strengthen the global search ability of refinement algorithm in the refinement phase.
Key words:multilevel method; partitioning; weighted undirected graph; spectral method
0 引言
對非規則、非結構的無向圖進行最優剖分，屬于組合優化范疇，即根據輸入的無向圖基本拓撲模型，在得到節點權值之和大致相等的節點子集的同時，使得處于不同子集的節點之間邊的權值之和最小。Kernighan等人[1]提出利用啟發信息迭代找到節點子集，對原有剖分進行組遷移優化的K-L算法，是第一個成功有效的啟發式算法。Fiduccia等人[2]在K-L算法的基礎上，引入節點的收益值和平衡約束條件，提出的F-M算法具備良好的收斂性和應用性。多水平方法是針對超大圖而發展起來的剖分技術[3~6]。該方法基于多水平思想，包含圖的粗化、初始剖分、投影優化三個階段。在粗化階段，將某些圖節點結合在一起，得到下一級粗化圖，重復此過程直到粗化圖足夠小為止，即得到最小圖；在初始剖分階段對該最小圖進行對分，得到初始剖分；最后在投影優化階段，將剖分從最小圖投影回初始圖，在每一水平層的粗化圖中對剖分進行優化。
譜圖論主要研究圖的鄰接矩陣和Laplacian矩陣的特征值和特征向量，建立圖的拓撲結構和圖矩陣表示的置換相似不變量之間的聯系。著名組合學家Nash-Williams在論及譜圖論時講道：“有些特征上是純組合的主要結果有這樣奇妙的特性，即如果不借助涉及圖的鄰接矩陣的特征根的代數方法，是不可能得到現在所知結論的。”Laplacian譜圖論及其在圖剖分中的應用從20世紀70年代以來有很多論述[7~9]。雖然譜方法得到的圖剖分質量比較好，但由于Laplacian矩陣對應的次小特征值計算量比較大，求解該向量占用了主要的執行時間，導致求解問題的時間是其他方法的幾個數量級[10]。基于多水平初始剖分階段的最小圖節點數一般為20~40，其相應的特征向量計算量相當有限。本文在初始剖分階段引入譜方法對最小圖進行初始剖分，提出了一種基于譜方法的無向賦權圖剖分算法SPWUG，為優化階段提供更好的初始剖分。
1 無向賦權圖的Laplacian譜圖論
由于多水平方法在初始剖分階段得到的最小圖為無向賦權圖，需要把基于非賦權無向圖的Laplacian譜圖論[10]在圖剖分的應用擴展到無向賦權圖上。
定義1 對于一無向賦權圖G=(V，E)給出的二剖分P={V1，V2}，該剖分的割切(cut)為屬于不同剖分節點之間的邊權值之和，即剖分割切記為
cut(P)=cut({V1，V2})=(v，u)∈E∧v∈V1∧u∈V2W(v，u) (1)
其中:W(v，u)指的是邊(v，u)的權值。
定義2 對于一無向賦權圖G=(V，E)，稱式(2)給出的A矩陣為圖G的鄰接矩陣，式(3)給出的D對角線矩陣為圖G的邊權值和矩陣，式(4)給出的L矩陣為圖G的Laplacian矩陣。
AG=0 v=u
W(v，u)v≠u(2)
DG=adjwgt(v) v=u
0v≠uwhere adjwgt(v)=u∈N(v)W(v，u)(3)
LG=DG-AG=adjwgt(v)v=u
-W(v，u)v≠u(4)
定義3 對于一無向賦權圖G=(V，E)給出的二剖分P={V1，V2}，對應于一個|V|維的剖分向量X=[X1，X2，…，X|V|]T。其中：Xv記為式(5)且滿足XT#8226;X=n。
Xv=+1 if v∈V1
-1 if v∈V2(5)
定理1 對于一無向賦權圖G=(V，E)給出的二剖分P={V1，V2}，其對應于剖分向量X滿足式(6)。
XT#8226;L#8226;X=(v，u)∈EW(v，u)(Xv-Xu)2(6)
證明 對于一無向賦權圖G=(V，E)中每一條邊，隨機地確定每一條邊的方向，從而對于任意一邊e∈E包含起始節點和終結節點。式(7)給出的|V|×|V|邊權值的對角線矩陣Q和式(8)給出的|V|×|E|定向關聯矩陣C=[Cve]。
Q=W(ei) i=j
0i≠j(7)
Cve=+1 if v is the initial vertex of e
-1if v is the terminal vertex of e
0otherwise(8)
基于式(5)(8)，CT#8226;X為|E|維向量Y=[Y1，Y2，…，Y|E|]T。其中：Ye=Xv-Xu且節點v是邊e的起始節點，節點u是邊e的終結節點。根據式(4)(7)(8)，可得L=C#8226;Q#8226;CT，代入式(6)經過推導，定理1得證于式(9)。
XT#8226;L#8226;X=XT#8226;C#8226;Q#8226;CT#8226;X=(CT#8226;X)T#8226;Q#8226;(CT#8226;X)=
Y#8226;Q#8226;Y=[…，(Xv-Xu)，…]#8226;Q#8226;[…，(Xv-Xu)，…]T=[…，(Xv-Xu)，…]#8226;[…，W(v，u)(Xv-Xu)，…]T=
(v，u)∈EW(v，u)(Xv-Xu)2(9)
定理2 對于一無向賦權圖G=(V，E)給出的二剖分P={V1，V2}，其對應于剖分向量X滿足
cut(P)=1/4XT#8226;L#8226;X(10)
證明 XT#8226;L#8226;X=(v，u)∈EW(v，u)(Xv-Xu)2包含三種情況：
a)邊e的節點(v，u)都處于V1中，即(v，u)∈E∧v∈V1∧u∈V1W(v，u)(Xv-Xu)2；
b)邊e的節點(v，u)都處于V2中，即(v，u)∈E∧v∈V2∧u∈V2W(v，u)(Xv-Xu)2；
c)邊e的節點(v，u)分別處于V1和V2中，即(v，u)∈E∧v∈V1∧u∈V2W(v，u)(Xv-Xu)2。
根據定義3可以得知，情況a)b)的累加和為零，且第三種情況中的|Xv-Xu|=2，因此定理2得證于式(11)。
XT#8226;L#8226;X=(v，u)∈E∧v∈V1∧u∈V2W(v，u)（|Xv-Xu|）2=
4×(v，u)∈E∧v∈V1∧u∈V2W(v，u)=4×cut(P)(11)
根據定義1和定理2，無向圖剖分優化問題可以轉換為式(12)所定義的優化問題。其中:e=(1，1，…，1)T。
minimize cut(P)=1/4XT#8226;L#8226;X
subject to eT#8226;X=0，and X#8226;XT=n(12)
針對式(12)所定義的問題應用Lagrange乘數因子法進行函數條件極值的求解。首先構造輔助函數F(x)=f(x)-ξ1g1(x)-ξ2g2(x)。其中：f(x)=X T#8226;L#8226;X，g1(x)=eT#8226;X，g2(x)=X T#8226;X-n；接著求輔助函數F(x)的一階偏導數并使之為零，即
Δf(x)-ξ1Δg1(x)-ξ2Δg2(x)=0，得到L#8226;X-ξ1e-ξ2X=0，即(L-ξ2I)#8226;X=ξ1e，左右兩邊都左乘上eT，得到eT#8226;L#8226;X-ξ2eT#8226;X=ξ1eT#8226;e。由于eT#8226;L#8226;X=0，eT#8226;X=0，eT#8226;e=n，得到ξ1=0。將ξ1=0代入公式(L-ξ2I)#8226;X=ξ1e，得到(L-ξ2I)#8226;X=0。
(L-ξ2I)#8226;X=0反映出：無向圖剖分問題在轉換為式(12)所定義的優化問題后，再次轉換為矩陣的特征值和特征向量的計算。在問題求解過程中需要對該問題本身的性質進行分析，以判斷所求得的點是否是極值點。根據圖Laplacian矩陣的定義，可以得到矩陣L的最小特征值λ1=0，且特征向量為[1，1，…，1]T，不符合無向圖剖分優化問題的要求。于是將式(12)所定義的優化問題的剖分向量取值范圍由{-1，+1}放松到實數空間，進行式(13)推導后，無向圖剖分問題最終轉換為對Laplacian矩陣次小特征值λ2及特征向量的計算。
瘙 綆 n，eT#8226;x=0，xT#8226;x=n，‖x‖2=n1/4×‖L#8226;X‖1≥1/4n|λ2|(13)
2 基于譜方法的無向賦權圖剖分算法
2.1 算法的描述
SPWUG算法的主要思想是：借助無向賦權圖對應的Laplacian矩陣的次小特征值給出最優剖分對應割切的下界，對應的特征向量作為度量無向賦權圖中各節點的相對距離，將特征向量中的分量進行排序之后對分得到初始剖分。算法的基本步驟如下：
a)清空節點子集V1和V2。
b)基于無向賦權圖G求解其相應的鄰接矩陣A、邊權值和對角線矩陣D；基于矩陣A和D，得到無向賦權圖G對應的Laplacian矩陣L。
c)求解出Laplacian矩陣L的次小特征值λ2和對應的特征向量 χ。
d)標注特征向量χ的分量都處于未訪問狀態，并按照特征向量χ的分量進行非降序排序。
e)選擇χ的非降序的第一個分量，以該分量起點并按照升序依次訪問處于未訪問狀態的分量，執行以下步驟：
(a)按照升序依次訪問處于未訪問狀態的分量，將該分量對應節點v加入節點子集V1，并標注節點v處于訪問狀態；
(b)如果節點子集V1權值之和滿足公式(1-r)S(V)/2≤S(V1)≤(1+r)S(V)/2，則訪問結束；否則重復上一步，直至訪問結束。
f)將處于未訪問狀態的節點加入到節點子集V2。
2.2 算法的關鍵實現
SPWUG算法實現的關鍵部分為第三個步驟，即基于Lanczos迭代[11，12]計算Laplacian矩陣次小特征值及特征向量。首先引入矩陣Q=[q1，q2，…，q|V|]，以及形如式(14)Jacobi矩陣形式的對稱三對角正交矩陣T，且滿足條件
QT#8226;L#8226;Q=T和QT#8226;Q=I。經過推導和展開得到βj+1qj+1=L#8226;qj-αjqj-βjqj-1，可以將矩陣L的特征值和特征向量計算轉換為矩陣T的特征值和特征向量計算。由于T滿足對稱三對角正交矩陣基本特性，可以采用雅可比方法求解出其特征值和特征向量。
基于Lanczos迭代計算矩陣L特征向量的偽代碼如下：
input:graph G′s Laplacian matrix L
output:matrix L′s eigenvalue λ and eigenvector χ
1 begin
2choose a unit-norm vectorq and set β1 =0
3for j=1 to n do begin
4wj+1=L#8226;qj-βj qj-1
5αj=qj#8226;wj+1
6wj+1=wj+1-αj qj
7βj+1=‖wj+1‖
8qj+1=wj+1/βj+1
9end for
10 T=QT#8226;L#8226;Q
11 compute T ′s eigenvalue λ and eigenvector v
12 χ=QT#8226;v
13 end
2.3 算法的復雜度分析
設初始剖分階段的最小圖節點數為n，邊數為m，將兩數算術操作作為復雜度分析的元運算。基于Lanczos迭代計算矩陣L特征向量算法的時間復雜度主要包括四部分：a)第3~9行的循環次數為n，主要運算集中在第4行的矩陣L與向量qj的乘運算，時間復雜度為Θ(n2)；b)第10行為矩陣之間的乘運算，時間復雜度為Θ(n3)；c)第11行采用雅可比方法計算對稱三對角正交矩陣的特征向量，時間復雜度為Θ(n3)；d)第12行為矩陣qj與向量v的乘運算，時間復雜度為Θ(n2)。因此總時間復雜度為Θ(n3)。
SPWUG算法的時間復雜度主要包括四部分：a)步驟b)計算L，時間復雜度為
Θ(max(m，n))；b)步驟c)基于Lanczos迭代計算矩陣特征向量，時間復雜度為Θ(n3)；c)步驟d)采用比較排序得到非降序排序，時間復雜度為Θ(n×log(n));d)步驟e)和f)將所有節點分為兩個節點子集，時間復雜度為Θ(n)。因此總時間復雜度為Θ(n3)，其計算復雜性主要集中在矩陣特征向量的計算上。由于SPWUG算法運行在初始剖分階段，其最小圖節點數一般為20~40，算法的計算量相當有限。
SPWUG算法的空間復雜度主要包括兩部分：a)壓縮存儲格式(compressed storage format，CSR)[13]的圖存儲空間。基于CSR存儲結構定義，所需空間大小為2×n+4×m。b)運算輔助空間。階數為n的矩陣有鄰接矩陣A、對角線矩陣D和Laplacian矩陣L；維數為n的向量有w和q，所需空間大小為3×n2+2×n。通常情況下用于剖分的圖實例滿足m>n，因此總空間復雜度為Ω(n2+m)。
3 實驗設計及結果分析
針對SPWUG算法，本文的實驗選用IBM公司的ISPD98電路測試基準[14]。它基于IBM公司Austin等設計部門給出的一系列內部電子產品經過轉換得到。本文采用了文獻[15]給出的一種ISPD98電路網表到圖的轉換算法，得到用于本實驗的測試基準圖。測試圖的基本特征如表1所示。節點數為12 752~210 613，邊數為109 183~2 235 716。經典算法實驗方案α為RM與SHEM組合粗化策略[3，13]、GGGP初始剖分算法[3，13]和BKL遷移優化算法[3，13]；對比實驗方案β為RM與SHEM組合粗化策略、SPWUG初始剖分算法和BKL遷移優化算法。采用的評估指標為20次實例驗證中得到的剖分割切最小值(minimum cut，Min-Cut)和剖分割切平均值(average cut，Ave-Cut)。該實驗結果如表1所示。可以觀察到，SPWUG算法在20次實例驗證中，相比GGGP初始剖分算法在Min-Cut獲得了7.9%的改進，在Ave-Cut獲得了6.9%的改進。
在ISPD98數據集的個別測試基準上，采用SPWUG算法之后導致剖分割切更差。對多水平方法的特點分析后可知，最小圖在拓撲結構上只是初始圖拓撲結構的一種近似，因此在最小圖上的全局近似最優剖分可能是初始圖的局部最優剖分。此外，在多水平優化階段采用的是BKL算法，它作為一種局部鄰域搜索方法采用的局部優化策略，是在當前解的鄰域中貪婪搜索，導致逃離局部最優的能力是有限的，使得個別測試基準上剖分割切更大。
4 結束語
本文借助輔助索引樹方法提高調整操作元運算的效率來彌補，基于多水平方法最小圖節點數一般為20~40的特點，提出了基于譜方法的無向賦權圖剖分算法，并應用在多水平方法的初始剖分階段。該算法借助無向賦權圖的Laplacian矩陣次小特征值對應的特征向量，刻畫了無向賦權圖中各節點的相對距離，計算出無向賦權圖中節點的相對距離并按序排列，然后將序列對分得到初始剖分。該算法的時間復雜度為Θ(n3)，空間復雜度為Ω(n2+m)，且實驗對比數據表明，SPWUG相比GGGP算法可以找到更優的剖分。針對實驗分析出來的問題，即在最小圖上的全局近似最優剖分可能是初始圖的局部最優剖分，今后將引入智能算法來改進多水平遷移優化算法，展開進一步的研究，增強其逃離局部最優的能力，從而能在優化階段找到初始圖的全局近似最優剖分。
參考文獻：
［1］KERNIGHAN B W，LIN S.An efficient heuristic procedure for partitioning graphs[J].Bell System Technical Journal，1970，49（1）:291-307.
［2］FIDUCCIA C，MATTHEYSES R.A linear-time heuristics for improving network partitions[C]//Proc of the 19th Conference on Design Automation.Piscataway:IEEE Press，1982:175-181.
［3］LENG Ming，YU Song-nian.An effective multi-level algorithm for bisecting graph [C]//Proc of the 2nd International Conference on Advanced Data Mining and Applications.Berlin:Springer， 2006:493-500.
［4］LENG Ming，YU Song-nian.An effective multi-level algorithm based on ant colony optimization for bisecting graph[C]//Proc of the 11th Pacific Asia Conference on Knowledge Discovery and Data Mining.Berlin:Springer，2007:138-149.
［5］SUN Ling-yu，LENG Ming.An effective refinement algorithm based on swarm intelligence for graph bipartitioning[C]//Proc of International Symposium on Combinatorics， Algorithms， Probabilistic and Experimental Methodologies.Berlin:Springer，2007:60-69.
［6］SUN Ling-yu，LENG Ming.An effective multi-level algorithm based on genetic algorithm for bisecting graph[J].Journal of Computational Information Systems，2008，4(4):1347-1354.
［7］PAIGE C C.Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem[J].Linear Algebra Appl，1980(34):235-258.
［8］POTHEN A，SIMON H.Towards a fast implementation of spectral nested dissection [C]//Proc of ACM/IEEE Conference on Super Computing.Berlin:Springer，1992:42-51.
［9］BARNARD S T，SIMON H D.A fast multilevel implementation of recursive spectral bisection for partitioning unstructed problems[C]//Proc of the 6th SIAM Conference on Parallel Processing for Scientific Computing. Philadelphia: SIAM，1993:711-718.
［10］DONGARRA J，FOX G，KENNEDY K，et al.Sourcebook of parallel computing[M].San Francisco: Morgan Kaufmann Publishers，2003.
［11］PARLETT B，SIMON H，STRINGER L.On estimating the Largest eigenvalue with the Lanczos algorithm[J].Mathematics of Computation，1982，38(137):153-165.
［12］STEWART G W.Matrix algorithms，volume Ⅱ:eigensystems[M].Philadelphia: SIAM，2001.
［13］KARYPIS G，KUMAR V.MeTiS 4.0:unstructured graphs partitioning and sparse matrix ordering system[R].Minnesota:University of Minnesota，1998.
［14］ALPERT C J.The ISPD98 circuit benchmark suite[C]//Proc of International Symposium on Physical Design. California:ACM Press，1998:80-85.
［15］孫凌宇，彭宣戈，冷明. 一種ISPD98電路網表到圖的轉換算法[J].井岡山大學學報:自然科學版，2008，29(4):19-21.
