葉緯材
(中山大學數學與計算科學學院//廣東省計算科學重點實驗室,廣東 廣州 510275)
稀疏矩陣向量乘法是大量迭代計算方法的核心。這些方法常用于求解包括線性方程組、奇異值分解、以及數據搜索引擎的網頁排序等問題。在迭代過程中,稀疏矩陣保持不變,而矩陣向量乘法被反復計算。提高并行稀疏矩陣向量乘法的效率就要解決稀疏矩陣的非零元素和輸入輸出向量元素的分布。并行化過程中出現的矩陣元素和輸入輸出向量的分劃問題被稱為并行稀疏矩陣向量乘法所對應的數據分劃問題。這是一個重要的組合數學問題,以下簡稱數據分劃問題。
對數據分劃問題,我們只關心矩陣的稀疏模式(sparsity pattern),即矩陣的元素aij值非 0 即 1。記Zn:={0,1,.…n-1}為指標集,n是個正整數。令矩陣等同于一個指標對集合,記為A:={(i,j):i∈Zm,j∈Zn}。記非零元素個數為nnz(A)。一個矩陣A的p部分劃(p-way partition)是一個以p個A的非空子集合為元素的集合,記為 {A0,A1,…,Ap-1}.這p個子集合是互不相交的,而且滿足∪i∈ZpAi=A。每個Ai稱為A的一部分(part)。記所有矩陣A的非零元素對應的指標對組成的集合為A1。由于數據分劃問題只與矩陣的非零元素有關,并與矩陣的分劃一一對應,下面就用集合A1的數據分劃定義矩陣A的數據分劃。
下面介紹并行稀疏矩陣向量乘法的計算過程。假設共有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的計算單元上;
任意計算架構下和任意數據分劃條件下的并行稀疏矩陣向量乘法都可以由以上的算法演變而成。
計算過程中,所有計算單元需要的非本地向量元素個數的總和稱為通信量。通信量與矩陣和輸入輸出向量的分劃相關。我們假設向量的分劃采用非對稱方式,并由矩陣分劃導出。這樣,數據分劃問題就只和矩陣的數據分劃有關。
本文討論的數據分劃問題定義如下:給定一個稀疏矩陣A和一個正整數p,尋找A1的一個p部分劃,不但要滿足負載平衡條件
(1)
而且通信量要低。常數ε∈[0,1) 稱作負載不平衡因子,表示最大可接受的計算負載不平衡程度。
通信量在分布式集群的節(jié)點之間表現為網絡帶寬的消耗量,而在節(jié)點內部表現為內存帶寬的消耗量。在多核計算架構下,通過數據分劃,不但減少內存帶寬的消耗量,而且降低內存爭用,從而進一步提高并行矩陣向量乘法的效率。
文章的結構如下:第一節(jié)以粗化函數的觀點統(tǒng)一現有的數據分劃方法。第二節(jié)給出一種新的自適應粗化函數構造方法,并給出該方法對應解空間的理論性質。第三節(jié)是分劃實驗。最后我們做一總結。
在本節(jié)里,我們用粗化函數的觀點統(tǒng)一回顧現有的數據分劃方法。首先通過分析矩陣分劃解的結構,導出求解數據分劃問題的過程就是尋找粗化函數的過程;然后提出一類顯式粗化函數選擇的方法、作用、以及與現有一維、二維矩陣分劃方法的關系。
Ak:={(i,j):f(A,i,j)=k,(i,j)∈A1},k∈Zp
矩陣分劃映射函數f都有以下的結構
f(A,i,j)=c°ι(i,j),(i,j)∈A1
(2)
這種將每個非零元素獨立分劃的方法稱為二維細粒度矩陣分劃。該問題可以建模為細粒度超圖分劃問題(fine-grained hypergraph partitioning problem),由文[1]提出。
文[2-7]等都作了不同方面的改進。由于一般超圖分劃問題是NP-hard問題[8],所以求解以上的問題的算法多數是啟發(fā)式的。文[1]提出的近似算法是現存文獻中平均性能最優(yōu)的,其多層迭代方法計算復雜度[9]約為 (nnz(A)log(nnz(A))。為了加速計算過程,最直接的方式就是對粗化函數的選擇加上約束,以降低求解問題的復雜度。
下面討論一類通過特殊粗化函數導出約束分劃問題解的方法。令c0為一個給定的、將Znnz(A)映射到Zk0上的滿射,其中k0>p.如果對應于矩陣A的p部數據分劃問題的解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)
記函數空間為F(A,c0°ι,p),其包含所有關于矩陣A的p部分劃函數f,而且這個空間是由函數c0°ι所引進的。該函數空間具體定義成以下形式
F(A,c0°ι,p)∶={f∶f(A,i,j)=
c°c0°ι(i,j),(i,j)∈A1}
(5)
引理1 令f為一個關于矩陣A的p部分劃問題的分劃映射函數。f∈F(A,c0°ι,p) 當且僅當f滿足條件(4)。
例如文[10]中討論的矩陣的一維(1-D)行分劃問題加入了粗化函數crow,保證同一行的非零元素將分屬同一部分。這個函數的形式為
crow°ι(i,j)=i,(i,j)∈A1
(6)
同理,矩陣的 1-D 列分劃問題加入了粗化函數ccol.這個函數的形式為
ccol°ι(i,j)=j,(i,j)∈A1
(7)
我們稱函數crow和ccol為 1-D 粗化函數。由于引入了約束, 1-D 分劃解空間就比二維( 2-D )細粒度分劃問題的要小的多;但同時分劃的質量有可能變差。盡管 2-D 細粒度分劃方案要靈活的多,但是對大型矩陣來說,計算的規(guī)模就增加不少。
近來出現了新的粗化函數選擇方法。2-D Mondriaan 方法基于分部 1-D 分劃[3],往往得到比傳統(tǒng) 1-D 方法優(yōu)勝的結果。角型分劃方法基于粗化函數ccorner-row和ccorner-col[4]。兩個函數的具體形式為
ccorner-row°ι(i,j)=min(i,j),(i,j)∈A1和
ccorner-col°ι(i,j)=max(i,j),(i,j)∈A1
對箭頭型矩陣等例子,角型方法得到的結果是可以與傳統(tǒng)的 2-D 方法的結果媲美;但也有部分例子,得到的結果比 1-D 方法的結果更差。以上的兩種方法的計算速度都比 2-D 細粒度方法要快得多。
在本節(jié)里,我們提出一種新的自適應粗化函數構造方法。目標是提出一種比 1-D 方法更靈活但計算復雜度相若的方法。為了簡化討論,我們在本節(jié)中只關注2部分劃問題。
(8)
其中νnew是一個序列化函數,將每個j∈CA*進行編號,從n到 |CA*|+n-1。因為CA*是Zn的非空子集合,所以|CA*|≤n。因此,cnew導出的解空間的規(guī)模與 1-D 方法相近。
下面的引理和定理給出了新方法在理論上的優(yōu)勢。
引理2 所有的列分劃解都屬于cnew導出的解空間,即
F(A,ccol°ι,2)?F(A,cnew°ι,2)
因此,cnew導出的分劃映射解本質是一種分片一維分劃,并且有以下優(yōu)勢:
(i) 選取最優(yōu)的行分劃初解,質量不比 1-D 方法差;
(ii) 選擇性粗化方法的粗化函數是動態(tài)自適應選擇的,這與“角方法”不同,靈活度提高;
(iii) 分劃的粒度比“Mordiaan 方法”要小得多,有利于找到更好的解;
(iv) 計算復雜度與 1-D 方法相若。
在本節(jié)中,我們將比較各種算法的分劃質量。參與比較的算法包括本文提出的選擇性粗化分劃算法(Selective)、 1-D 超圖行分劃方法[10]、 2-D Mondriaan方法[3]和 2-D 細粒度超圖分劃方法[1]。 負載不平衡因子設定為0.03。所有的實驗都運行20次,取最佳結果。
實驗在Intel Core T2050處理器下運行,內存容量為1GB。實驗操作系統(tǒng): Ubuntu Linux 8.04.每一個程序都用單進程和單核心,獨占系統(tǒng)運行。
測試用的矩陣都是在網絡上可以取得的,選自不同的問題領域。具體的矩陣信息,請參閱表1。
表1 測試用矩陣的性質
由表2所列的數據表明,2-D 細粒度超圖分劃方法的平均質量是最好的。這是由于其分劃空間理論上是各種方法中最大的。本文提出的選擇性粗化分劃方法也得到了很不錯的效果。對矩陣cage10、Dubcova2和 garon2,選擇性粗化分劃方法得到的結果勝過所有其他方法的結果,這就驗證了選擇性粗化分劃方法其分劃空間在理論上要比1-D 分劃方法要好的結論。表2中(*)表示結果不滿足負載平衡條件。
表2 分劃結果的通信量
本章提出以粗化函數的觀點統(tǒng)一現有的數據分劃方法,并且提出一種新的粗化函數選擇方法。通過理論分析與實驗證明,本文提出的選擇性粗化方法比 1-D 分劃方法好,接近甚至超過 2-D 細粒度方法的結果。
參考文獻:
[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.