选择粗化函数优化并行稀疏矩阵向量乘法*

叶纬材

(中山大学数学与计算科学学院//广东省计算科学重点实验室，广东 广州 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)每个计算单元发送其拥有的向量元素xj到其他拥有 第j列非零元素的计算单元；

2)yik←∑aij∈Akaijxj，其中yik为第i行在计算单元k上的的部分和，k是计算单元的序号，从 0 到p-1 ；

3)第k号计算单元发送非零的yik到拥有yi的计算单元上；

任意计算架构下和任意数据分划条件下的并行稀疏矩阵向量乘法都可以由以上的算法演变而成。

计算过程中，所有计算单元需要的非本地向量元素个数的总和称为通信量。通信量与矩阵和输入输出向量的分划相关。我们假设向量的分划采用非对称方式，并由矩阵分划导出。这样，数据分划问题就只和矩阵的数据分划有关。

本文讨论的数据分划问题定义如下：给定一个稀疏矩阵A和一个正整数p，寻找A1的一个p部分划，不但要满足负载平衡条件

(1)

而且通信量要低。常数ε∈[0,1) 称作负载不平衡因子，表示最大可接受的计算负载不平衡程度。

通信量在分布式集群的节点之间表现为网络带宽的消耗量，而在节点内部表现为内存带宽的消耗量。在多核计算架构下，通过数据分划，不但减少内存带宽的消耗量，而且降低内存争用，从而进一步提高并行矩阵向量乘法的效率。

文章的结构如下：第一节以粗化函数的观点统一现有的数据分划方法。第二节给出一种新的自适应粗化函数构造方法，并给出该方法对应解空间的理论性质。第三节是分划实验。最后我们做一总结。

1 粗化函数的选择与一维、二维矩阵分划

在本节里，我们用粗化函数的观点统一回顾现有的数据分划方法。首先通过分析矩阵分划解的结构，导出求解数据分划问题的过程就是寻找粗化函数的过程；然后提出一类显式粗化函数选择的方法、作用、以及与现有一维、二维矩阵分划方法的关系。

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]，所以求解以上的问题的算法多数是启发式的。文[1]提出的近似算法是现存文献中平均性能最优的，其多层迭代方法计算复杂度[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 细粒度分划方案要灵活的多，但是对大型矩阵来说，计算的规模就增加不少。

近来出现了新的粗化函数选择方法。2-D Mondriaan 方法基于分部 1-D 分划[3]，往往得到比传统 1-D 方法优胜的结果。角型分划方法基于粗化函数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

对箭头型矩阵等例子，角型方法得到的结果是可以与传统的 2-D 方法的结果媲美；但也有部分例子，得到的结果比 1-D 方法的结果更差。以上的两种方法的计算速度都比 2-D 细粒度方法要快得多。

2 粗化函数的选择与分片一维分划

在本节里，我们提出一种新的自适应粗化函数构造方法。目标是提出一种比 1-D 方法更灵活但计算复杂度相若的方法。为了简化讨论，我们在本节中只关注2部分划问题。

(8)

其中νnew是一个序列化函数，将每个j∈CA*进行编号，从n到 |CA*|+n-1。因为CA*是Zn的非空子集合，所以|CA*|≤n。因此，cnew导出的解空间的规模与 1-D 方法相近。

下面的引理和定理给出了新方法在理论上的优势。

引理2 所有的列分划解都属于cnew导出的解空间，即

F(A,ccol∘ι,2)⊂F(A,cnew∘ι,2)

因此，cnew导出的分划映射解本质是一种分片一维分划，并且有以下优势：

(i) 选取最优的行分划初解，质量不比 1-D 方法差；

(ii) 选择性粗化方法的粗化函数是动态自适应选择的，这与“角方法”不同，灵活度提高；

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

(iv) 计算复杂度与 1-D 方法相若。

3 实 验

在本节中，我们将比较各种算法的分划质量。参与比较的算法包括本文提出的选择性粗化分划算法(Selective)、 1-D 超图行分划方法[10]、 2-D Mondriaan方法[3]和 2-D 细粒度超图分划方法[1]。 负载不平衡因子设定为0.03。所有的实验都运行20次，取最佳结果。

实验在Intel Core T2050处理器下运行，内存容量为1GB。实验操作系统： Ubuntu Linux 8.04.每一个程序都用单进程和单核心，独占系统运行。

测试用的矩阵都是在网络上可以取得的，选自不同的问题领域。具体的矩阵信息，请参阅表1。

表1 测试用矩阵的性质

由表2所列的数据表明，2-D 细粒度超图分划方法的平均质量是最好的。这是由于其分划空间理论上是各种方法中最大的。本文提出的选择性粗化分划方法也得到了很不错的效果。对矩阵cage10、Dubcova2和 garon2,选择性粗化分划方法得到的结果胜过所有其他方法的结果，这就验证了选择性粗化分划方法其分划空间在理论上要比1-D 分划方法要好的结论。表2中(*)表示结果不满足负载平衡条件。

表2 分划结果的通信量

4 总 结

本章提出以粗化函数的观点统一现有的数据分划方法，并且提出一种新的粗化函数选择方法。通过理论分析与实验证明，本文提出的选择性粗化方法比 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.