葉緯材

(中山大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院//廣東省計算科學(xué)重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510275)

稀疏矩陣向量乘法是大量迭代計算方法的核心。這些方法常用于求解包括線性方程組、奇異值分解、以及數(shù)據(jù)搜索引擎的網(wǎng)頁排序等問題。在迭代過程中，稀疏矩陣保持不變，而矩陣向量乘法被反復(fù)計算。提高并行稀疏矩陣向量乘法的效率就要解決稀疏矩陣的非零元素和輸入輸出向量元素的分布。并行化過程中出現(xiàn)的矩陣元素和輸入輸出向量的分劃問題被稱為并行稀疏矩陣向量乘法所對應(yīng)的數(shù)據(jù)分劃問題。這是一個重要的組合數(shù)學(xué)問題，以下簡稱數(shù)據(jù)分劃問題。

對數(shù)據(jù)分劃問題，我們只關(guān)心矩陣的稀疏模式(sparsity pattern)，即矩陣的元素aij值非 0 即 1。記Zn:={0,1,.…n-1}為指標(biāo)集，n是個正整數(shù)。令矩陣等同于一個指標(biāo)對集合，記為A:={(i,j):i∈Zm,j∈Zn}。記非零元素個數(shù)為nnz(A)。一個矩陣A的p部分劃(p-way partition)是一個以p個A的非空子集合為元素的集合，記為 {A0,A1,…,Ap-1}.這p個子集合是互不相交的，而且滿足∪i∈ZpAi=A。每個Ai稱為A的一部分(part)。記所有矩陣A的非零元素對應(yīng)的指標(biāo)對組成的集合為A1。由于數(shù)據(jù)分劃問題只與矩陣的非零元素有關(guān)，并與矩陣的分劃一一對應(yīng)，下面就用集合A1的數(shù)據(jù)分劃定義矩陣A的數(shù)據(jù)分劃。

下面介紹并行稀疏矩陣向量乘法的計算過程。假設(shè)共有p個計算單元，輸入矩陣A為一個m×n的矩陣。每個計算單元擁有矩陣A的非零元素以及輸入輸出向量元素的一部分。擁有矩陣非零元素aij的計算單元求出乘積aijxj。若aij、xj和yi屬于同一個計算單元，則計算過程為本地的，否則就需要通信。一個自然的并行稀疏矩陣向量乘法由以下四個步驟組成：

1)每個計算單元發(fā)送其擁有的向量元素xj到其他擁有 第j列非零元素的計算單元；

2)yik←∑aij∈Akaijxj，其中yik為第i行在計算單元k上的的部分和，k是計算單元的序號，從 0 到p-1 ；

3)第k號計算單元發(fā)送非零的yik到擁有yi的計算單元上；

任意計算架構(gòu)下和任意數(shù)據(jù)分劃條件下的并行稀疏矩陣向量乘法都可以由以上的算法演變而成。

計算過程中，所有計算單元需要的非本地向量元素個數(shù)的總和稱為通信量。通信量與矩陣和輸入輸出向量的分劃相關(guān)。我們假設(shè)向量的分劃采用非對稱方式，并由矩陣分劃導(dǎo)出。這樣，數(shù)據(jù)分劃問題就只和矩陣的數(shù)據(jù)分劃有關(guān)。

本文討論的數(shù)據(jù)分劃問題定義如下：給定一個稀疏矩陣A和一個正整數(shù)p，尋找A1的一個p部分劃，不但要滿足負(fù)載平衡條件

(1)

而且通信量要低。常數(shù)ε∈[0,1) 稱作負(fù)載不平衡因子，表示最大可接受的計算負(fù)載不平衡程度。

通信量在分布式集群的節(jié)點(diǎn)之間表現(xiàn)為網(wǎng)絡(luò)帶寬的消耗量，而在節(jié)點(diǎn)內(nèi)部表現(xiàn)為內(nèi)存帶寬的消耗量。在多核計算架構(gòu)下，通過數(shù)據(jù)分劃，不但減少內(nèi)存帶寬的消耗量，而且降低內(nèi)存爭用，從而進(jìn)一步提高并行矩陣向量乘法的效率。

文章的結(jié)構(gòu)如下：第一節(jié)以粗化函數(shù)的觀點(diǎn)統(tǒng)一現(xiàn)有的數(shù)據(jù)分劃方法。第二節(jié)給出一種新的自適應(yīng)粗化函數(shù)構(gòu)造方法，并給出該方法對應(yīng)解空間的理論性質(zhì)。第三節(jié)是分劃實(shí)驗(yàn)。最后我們做一總結(jié)。

1 粗化函數(shù)的選擇與一維、二維矩陣分劃

在本節(jié)里，我們用粗化函數(shù)的觀點(diǎn)統(tǒng)一回顧現(xiàn)有的數(shù)據(jù)分劃方法。首先通過分析矩陣分劃解的結(jié)構(gòu)，導(dǎo)出求解數(shù)據(jù)分劃問題的過程就是尋找粗化函數(shù)的過程；然后提出一類顯式粗化函數(shù)選擇的方法、作用、以及與現(xiàn)有一維、二維矩陣分劃方法的關(guān)系。

Ak:={(i,j):f(A,i,j)=k,(i,j)∈A1},k∈Zp

矩陣分劃映射函數(shù)f都有以下的結(jié)構(gòu)

f(A,i,j)=c°ι(i,j),(i,j)∈A1

(2)

這種將每個非零元素獨(dú)立分劃的方法稱為二維細(xì)粒度矩陣分劃。該問題可以建模為細(xì)粒度超圖分劃問題(fine-grained hypergraph partitioning problem)，由文[1]提出。

文[2-7]等都作了不同方面的改進(jìn)。由于一般超圖分劃問題是NP-hard問題[8]，所以求解以上的問題的算法多數(shù)是啟發(fā)式的。文[1]提出的近似算法是現(xiàn)存文獻(xiàn)中平均性能最優(yōu)的，其多層迭代方法計算復(fù)雜度[9]約為 (nnz(A)log(nnz(A))。為了加速計算過程，最直接的方式就是對粗化函數(shù)的選擇加上約束，以降低求解問題的復(fù)雜度。

下面討論一類通過特殊粗化函數(shù)導(dǎo)出約束分劃問題解的方法。令c0為一個給定的、將Znnz(A)映射到Zk0上的滿射，其中k0>p.如果對應(yīng)于矩陣A的p部數(shù)據(jù)分劃問題的解f都有下面的形式

f(A,i,j)=c°c0°ι(i,j),(i,j)∈A1

(3)

c0°ι(i1,j1)=c0°ι(i2,j2)?f(A,i1,j1)=

f(A,i2,j2),?(i1,j1),(i2,j2)∈A1

(4)

記函數(shù)空間為F(A,c0°ι,p),其包含所有關(guān)于矩陣A的p部分劃函數(shù)f,而且這個空間是由函數(shù)c0°ι所引進(jìn)的。該函數(shù)空間具體定義成以下形式

F(A,c0°ι,p)∶={f∶f(A,i,j)=

c°c0°ι(i,j),(i,j)∈A1}

(5)

引理1 令f為一個關(guān)于矩陣A的p部分劃問題的分劃映射函數(shù)。f∈F(A,c0°ι,p) 當(dāng)且僅當(dāng)f滿足條件(4)。

例如文[10]中討論的矩陣的一維(1-D)行分劃問題加入了粗化函數(shù)crow,保證同一行的非零元素將分屬同一部分。這個函數(shù)的形式為

crow°ι(i,j)=i,(i,j)∈A1

(6)

同理，矩陣的 1-D 列分劃問題加入了粗化函數(shù)ccol.這個函數(shù)的形式為

ccol°ι(i,j)=j,(i,j)∈A1

(7)

我們稱函數(shù)crow和ccol為 1-D 粗化函數(shù)。由于引入了約束， 1-D 分劃解空間就比二維( 2-D )細(xì)粒度分劃問題的要小的多；但同時分劃的質(zhì)量有可能變差。盡管 2-D 細(xì)粒度分劃方案要靈活的多，但是對大型矩陣來說，計算的規(guī)模就增加不少。

近來出現(xiàn)了新的粗化函數(shù)選擇方法。2-D Mondriaan 方法基于分部 1-D 分劃[3]，往往得到比傳統(tǒng) 1-D 方法優(yōu)勝的結(jié)果。角型分劃方法基于粗化函數(shù)ccorner-row和ccorner-col[4]。兩個函數(shù)的具體形式為

ccorner-row°ι(i,j)=min(i,j),(i,j)∈A1和

ccorner-col°ι(i,j)=max(i,j),(i,j)∈A1

對箭頭型矩陣等例子，角型方法得到的結(jié)果是可以與傳統(tǒng)的 2-D 方法的結(jié)果媲美；但也有部分例子，得到的結(jié)果比 1-D 方法的結(jié)果更差。以上的兩種方法的計算速度都比 2-D 細(xì)粒度方法要快得多。

2 粗化函數(shù)的選擇與分片一維分劃

在本節(jié)里，我們提出一種新的自適應(yīng)粗化函數(shù)構(gòu)造方法。目標(biāo)是提出一種比 1-D 方法更靈活但計算復(fù)雜度相若的方法。為了簡化討論，我們在本節(jié)中只關(guān)注2部分劃問題。

(8)

其中νnew是一個序列化函數(shù)，將每個j∈CA*進(jìn)行編號，從n到 |CA*|+n-1。因?yàn)镃A*是Zn的非空子集合，所以|CA*|≤n。因此，cnew導(dǎo)出的解空間的規(guī)模與 1-D 方法相近。

下面的引理和定理給出了新方法在理論上的優(yōu)勢。

引理2 所有的列分劃解都屬于cnew導(dǎo)出的解空間，即

F(A,ccol°ι,2)?F(A,cnew°ι,2)

因此，cnew導(dǎo)出的分劃映射解本質(zhì)是一種分片一維分劃，并且有以下優(yōu)勢：

(i) 選取最優(yōu)的行分劃初解，質(zhì)量不比 1-D 方法差；

(ii) 選擇性粗化方法的粗化函數(shù)是動態(tài)自適應(yīng)選擇的，這與“角方法”不同，靈活度提高；

(iii) 分劃的粒度比“Mordiaan 方法”要小得多，有利于找到更好的解；

(iv) 計算復(fù)雜度與 1-D 方法相若。

3 實(shí) 驗(yàn)

在本節(jié)中，我們將比較各種算法的分劃質(zhì)量。參與比較的算法包括本文提出的選擇性粗化分劃算法(Selective)、 1-D 超圖行分劃方法[10]、 2-D Mondriaan方法[3]和 2-D 細(xì)粒度超圖分劃方法[1]。 負(fù)載不平衡因子設(shè)定為0.03。所有的實(shí)驗(yàn)都運(yùn)行20次，取最佳結(jié)果。

實(shí)驗(yàn)在Intel Core T2050處理器下運(yùn)行，內(nèi)存容量為1GB。實(shí)驗(yàn)操作系統(tǒng)： Ubuntu Linux 8.04.每一個程序都用單進(jìn)程和單核心，獨(dú)占系統(tǒng)運(yùn)行。

測試用的矩陣都是在網(wǎng)絡(luò)上可以取得的，選自不同的問題領(lǐng)域。具體的矩陣信息，請參閱表1。

表1 測試用矩陣的性質(zhì)

由表2所列的數(shù)據(jù)表明，2-D 細(xì)粒度超圖分劃方法的平均質(zhì)量是最好的。這是由于其分劃空間理論上是各種方法中最大的。本文提出的選擇性粗化分劃方法也得到了很不錯的效果。對矩陣cage10、Dubcova2和 garon2,選擇性粗化分劃方法得到的結(jié)果勝過所有其他方法的結(jié)果，這就驗(yàn)證了選擇性粗化分劃方法其分劃空間在理論上要比1-D 分劃方法要好的結(jié)論。表2中(*)表示結(jié)果不滿足負(fù)載平衡條件。

表2 分劃結(jié)果的通信量

4 總 結(jié)

本章提出以粗化函數(shù)的觀點(diǎn)統(tǒng)一現(xiàn)有的數(shù)據(jù)分劃方法，并且提出一種新的粗化函數(shù)選擇方法。通過理論分析與實(shí)驗(yàn)證明，本文提出的選擇性粗化方法比 1-D 分劃方法好，接近甚至超過 2-D 細(xì)粒度方法的結(jié)果。

參考文獻(xiàn)：

[1]CATALYUREK U,AYKANAT C.A fine-grain hypergraph model for 2D decomposition of sparse matrices [C]//Proceedings of 15th International Parallel and Distributed Processing Symposium(IPDPS),San Francisco,CA,2001.

[2]CATALYUREK U,AYKANAT C.A hypergraph-partitioning approach for coarse-grain decomposition [C]// ACM/IEEE SC2001,Denver,CO,2001.

[3]VASTENHOUW B,BISSELING ROB H.A two-dimensional data distribution method for parallel sparse matrix-vector multiplication [J].SIAM Review,2005,47(1): 67-95.

[4]WOLF M M,BOMAN E G,HENDRICKSON B.Optimizing parallel sparse matrix-vector multiplication by corner partitioning [C]//PARA08,Trondheim,Norway,2008.

[5]CATALYUREK UMIT V,BOMAN ERIK G.Hypergraph-based dynamic load balancing for adaptive scientific computations [C]//Proceedings of IPDPS’07,Best Algorithms Paper Award,2007.

[6]BOMAN ERIK G.A nested dissection approach for sparse matrix partitioning [C]// Proceedings in Applied Mathematics and Mechanics,2008.

[7]HENDRICKSON B.Graph partitioning and parallel solvers: has the emperor no clothes? [J].Lecture Notes in Computer Science,1998,1457: 218-225.

[8]LENGAUER T.Combinatorial algorithms for integrated circuit layout [M].UK: Wiley,1990.

[9]SAAB Y.Combinatorial optimization by dynamic contraction [J].Journal of Heuristics,1997,3: 207-224.

[10]CATALYUREK U,AYKANAT C.Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication [J].IEEE Transaction on Parallel and Distributed Systems,1999,10(7): 673-693.