【BFS应用】两种层同步并行BFS算法与SpMV和SpMSpV之间的对应关系

Lord King

已于 2023-03-28 20:57:59 修改

阅读量1.1k

点赞数 7

文章标签： bfs 图搜索算法矩阵

于 2023-03-28 20:40:02 首次发布

本文链接：https://blog.csdn.net/m0_63104232/article/details/129775219

版权

BFS：Breadth-First Search，广度优先搜索

SpMV：Sparse-Matrix (Dense-)Vector Multiplication，稀疏矩阵向量乘法

SpMSpV：Sparse-Matrix Sparse-Vector Multiplication，稀疏矩阵稀疏向量乘法

基于CSR的SpMV计算方法

CSR（Compressed Sparse Row，压缩行存储）是存储稀疏矩阵的一种有效方式，避免了使用二维数组方式时存储大量0值的情况。事实上，CSR对矩阵的稀疏性没有要求，是一种适用于任何矩阵的通用存储方法，在诸稀疏矩阵的存储方式中也不见得是最高效的。

Example 1: 图1是来自Parallel Programming for FPGAs书中的举例，a)和b)分别展示了矩阵M的二维数组表示和CSR表示。CSR作为一种存储方式，也作为一种数据结构，由三个数组构成（所有的index默认从0开始）：

values：存储M中的非零元素的值；

columnIndex：存储非零元素在二维数组中水平方向的位置，即columnIndex[i]记录了values[i]元素的列序号；

rowPtr：该数组的长度为矩阵的行数加1，其中rowPtr[0]=0，其后的rowPtr[i]记录了矩阵前i行（即序号0至i-1的行）中包含非零元素的数量。

图1
Reference: Kastner, R., Matai, J., & Neuendorffer, S. (2018). Parallel Programming for FPGAs. *ArXiv, abs/1805.03648*.

我们可以来模拟一下矩阵M（CSR表示）和向量 $gif.latex?x%3D%28x_%7B0%7D%2C%20x_%7B1%7D%2C%20x_%7B2%7D%2C%20x_%7B3%7D%29%5E%7BT%7D$ 的乘法 $gif.latex?Mx%3Dy%3D%28y_%7B0%7D%2C%20y_%7B1%7D%2C%20y_%7B2%7D%2C%20y_%7B3%7D%29%5E%7BT%7D$ 的计算过程：

首先要用M的第一行向量和x相乘，由rowPtr[1]-rowPtr[0]=2知第一行共有2个非零元素；

从columnIndex中找到这两个元素所处的列，即columnIndex[0]=0和columnIndex[1]=1，这就意味着只需要考虑 $gif.latex?x_%7B0%7D$ 和 $gif.latex?x_%7B1%7D$ 分别与该行前两列元素的乘积（因为其余列元素都是0）；

再从values中找到这两个元素的值values[0]=3和values[1]=4，分别计算values[0]× $gif.latex?x_%7B0%7D$ 和values[1]× $gif.latex?x_%7B1%7D$ 并将两个乘积加起来即得到 $gif.latex?y_%7B0%7D$ 的值。

后面 $gif.latex?y_%7Bi%7D$ 的计算都是先从rowPtr[i+1]-rowPtr[i]=k开始，然后查找columnIndex[rowPtr[i]]至columnIndex[rowPtr[i+k-1]]这k个列序号值，再然后根据列序号值在values中查找对应位置的k个元素值，分别与 $gif.latex?x_%7BrowPtr%5Bi%5D%7D$ 至 $gif.latex?x_%7BrowPtr%5Bi+k-1%5D%7D$ 相乘后将k个乘积相加得到 $gif.latex?y_%7Bi%7D$ 的值。

在程序实现中，上述操作需要2个嵌套循环，外层循环访问每一行，内层循环访问当前行中每一列元素。具体的实现代码从略。

尽管上文并没有刻意给出一个十分“稀疏”矩阵的例子，但也容易分析知：相比于二维数组，采用CSR存储方式避免了大量乘零冗余运算，的确能提高SpMV的实现效率。

由SpV切入的SpMSpV计算方法

在前一部分中，我们讨论了基于CSR的矩阵向量乘法，其核心思想是从矩阵的压缩着手，想办法避免由矩阵中的零元素造成的冗余运算。由（前一部分）Mx的计算过程可以看出，矩阵要先“主动宣布”自己在当前行有哪些列上的元素是非零的，然后向量再“被动派出”与这些非零元素对应的分量（“对应”是指非零元素的列序号与分量的行序号相等）与之相乘。当向量也十分稀疏时，它所“派出”的分量也有很大概率为零，这同样会导致相当多的冗余运算。既然由SpM切入可以使SpMV的计算更加高效，那么由SpV切入能否改善SpMSpV的计算效率呢？

Example 2：

$gif.latex?M%3D%28a_%7Bij%7D%29_%7Bm%5Ctimes%20n%7D%280%5Cleqslant%20i%5Cleqslant%20m-1%2C%200%5Cleqslant%20j%5Cleqslant%20n-1%29$ $gif.latex?x%3D%280%2C%20...%2C%200%2C%20x_%7Br%7D%2C%200%2C%20...%2C%200%2C%20x_%7Bs%7D%2C%200%2C%20...%2C%200%29%5E%7BT%7D%280%5Cleqslant%20r%3C%20s%5Cleqslant%20n-1%29$

$gif.latex?Mx%3Dy%3D%28y_%7B0%7D%2C%20y_%7B1%7D%2C%20...%2C%20y_%7Bn-1%7D%29%5E%7BT%7D$

$gif.latex?%5CRightarrow%20y_%7Bk%7D%3Dx_%7Br%7Da_%7Bkr%7D+x_%7Bs%7Da_%7Bks%7D%280%5Cleqslant%20k%5Cleqslant%20m-1%29$

由上例 $gif.latex?y_%7Bk%7D$ 的表达式可以看出，对于仅含两个非零分量 $gif.latex?x_%7Br%7D$ 和 $gif.latex?x_%7Bs%7D$ 的稀疏向量x，想要计算出M和x的乘积y，需要进行的乘法运算仅有 $gif.latex?x_%7Br%7Da_%7B0r%7D%2C%20x_%7Br%7Da_%7B1r%7D%2C%20...%2C%20x_%7Br%7Da_%7Bm-1%2Cr%7D$ 以及 $gif.latex?x_%7Bs%7Da_%7B0s%7D%2C%20x_%7Bs%7Da_%7B1s%7D%2C%20...%2C%20x_%7Bs%7Da_%7Bm-1%2Cs%7D$ ，也就是说，参与Mx计算的M中的元素仅有列序号为r和s（和 $gif.latex?x_%7Br%7D$ 、 $gif.latex?x_%7Bs%7D$ 的行序号对应相等）的两列。这样一来，我们只需要根据向量非零分量的行序号，选出矩阵对应列上的全部非零元素，例如我们假设列r的 $gif.latex?a_%7Bkr%7D%20%28k%3D0%2C%201%2C%20...%2C%20m-1%29$ 中仅有 $gif.latex?a_%7Bpr%7D$ 和 $gif.latex?a_%7Bqr%7D$ 这两个非零元素，那么只需要将二者分别与 $gif.latex?x_%7Br%7D$ 相乘，将两个乘积分别暂时保存在两个列表中，这两个列表分别存储所有形如 $gif.latex?x_%7B*%7Da_%7Bp*%7D$ 和 $gif.latex?x_%7B*%7Da_%7Bq*%7D$ 的结果，这是为了方便最后计算 $gif.latex?y_%7Bp%7D$ 和 $gif.latex?y_%7Bq%7D$ （直接求对应列表中全部元素的和即得）。可以看到，此时SpV成为了“主动”的一方。

从上面所描述的计算过程可以看出，这里我们需要一种存储（或表示）方法，使得按列获取矩阵的非零元素更为便捷。第一部分中提到的CSR是按行优先访问非零元素的，它对行进行压缩；我们这里需要按列优先访问非零元素，那当然也可以对列进行压缩，这便是与CSR对应的CSC（Compressed Sparse Column，压缩列存储）。CSC与CSR的原理和结构极为类似，只不过是将“行”和“列”互换了而已，因此这里就不再详述了。

我们再来考虑一个更为简化的问题：将M和x中的非零元素都置为1，计算Mx。具体如下：

Example 3: Mx=y

$gif.latex?%5Cbegin%7Bbmatrix%7D%200%260%20%26%201%20%26%200%20%26%200%20%26%200%20%26%200%20%260%20%5C%5C%201%260%20%260%20%260%20%26%200%20%260%20%26%200%20%260%20%5C%5C%200%26%200%20%260%20%261%20%26%200%20%26%200%20%26%200%260%20%5C%5C%201%260%20%26%200%20%260%20%260%20%260%20%260%20%260%20%5C%5C%200%26%201%20%26%200%260%20%260%20%26%201%26%200%20%260%20%5C%5C%200%26%201%20%26%200%20%26%201%20%26%200%20%26%200%26%200%20%260%20%5C%5C%200%26%200%20%26%201%26%200%26%200%20%26%200%26%200%20%261%20%5C%5C%200%26%200%26%200%20%26%201%260%20%260%20%260%20%260%20%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D%200%5C%5C%201%5C%5C%200%5C%5C%201%5C%5C%200%5C%5C%200%5C%5C%200%5C%5C%200%20%5Cend%7Bbmatrix%7D%3D%5Cbegin%7Bbmatrix%7D%200%5C%5C%200%5C%5C%201%5C%5C%200%5C%5C%201%5C%5C%202%5C%5C%200%5C%5C%201%20%5Cend%7Bbmatrix%7D$

如果我们只关注结果向量有哪些分量（行）非零，那么就根据向量x的第2、4两个分量非零，去矩阵M中查看第2、4两列，发现两列中的非零元素分布于第3、5、6、8四行（不重复列举），于是结果向量中第3、5、6、8四个分量就是非零的。如果要求出这四个分量具体的值，也很简单，第6行第2、4两列都是1，而其余没行都只有一列是1，所以结果向量第6个分量是2，第3、5、8三个分量都是1。可结合上面右边的图示来思考这一计算过程，这张图在后文会再次出现。

两种层同步并行BFS算法

层同步并行BFS有两种算法，分别是自顶向下（top-down）和自底向上（bottom-up）方法。通俗来讲，我们熟知的BFS树扩展方向是自顶向下，也就是parent-to-child的；而还有一种child-to-parent的BFS树扩展方向，就是自底向上。本部分主要介绍这两种算法的操作过程，以及引入自底向上方法的原因。

需要声明的一点是，本文并不对什么是“层同步并行”做解释（笔者本人也不懂这方面的知识），这里只需要了解这两种算法的执行过程，看懂其中的抽象操作即可（对这两种算法的分析也不需要有关“层同步并行”的知识）。下文所写多为简要概括性的内容，如有兴趣了解更多，可以参考Beamer, S., Asanović, K., & Patterson, D. (2013). Direction-optimizing breadth-first search. Scientific Programming, 21(3-4).这篇论文。

层同步并行BFS算法（Level Synchronized BFS）是最常用的并行图遍历方法：

1. 从任意随机顶点开始，作为当前层活跃顶点（Current Frontier）；

2. 由Current Frontier向未访问的邻接顶点（Neighbor）扩展；

3. 新遍历的顶点集合成为下一层遍历迭代的活跃顶点（Next Frontier）。

图2
Reference: https://www.ccita.net/wp-content/uploads/dlm_uploads/2021/01/4.2.4-面向高通量计算机的图算法优化技术-（中科院计算所-曹华伟）.pdf

Top-Down方法：

1. 扫描Current Frontier的所有顶点；

2. 检查其所有的Neighbor；

3. 添加未访问的Neighbor到Next Frontier。

Remark: Current Frontier活跃顶点数量较多时，遍历过程存在大量的冗余检查；并行执行时存在更新冲突。



Bottom-Up方法：

1. 扫描未访问的顶点；

2. 检查其Neighbor是否在Current Frontier，若有，则终止并检查下一顶点；

3. 将该顶点添加到Next Frontier。

Remark: 从未访问的顶点开始遍历，发现Neighbor位于Current Frontier即跳出循环。

Pseudocode Reference: 深入浅出BFS(1)_delibk的博客-CSDN博客

C++ Implementation: 深入浅出BFS(2)：自顶向下(Top-down method) & 自底向上(Bottom-up method)_自底向上bfs_delibk的博客-CSDN博客

BFS与矩阵向量乘法的联系

对于笔者这样的初学者来说，能够发现两个看似毫不相干概念之间的某种内在联系，无疑是一件令人兴奋的事情。让我们回头看看前文的Example 3，例子中给出的矩阵M实际上正是下面这个（有向无权）图G的邻接矩阵（Adjacency Matrix）：

为方便表示，设置该图的顶点表为： $\begin{pmatrix} i & 0 & 1 & 2 & 3 & 4 & 5 & 6 &7 \\ vertices[i]& A & B & C & D & E & F & G & H \end{pmatrix}$ .

构造这个邻接矩阵的规则是： $M=(a_{ij})_{8\times 8},\:0\leqslant i,j\leqslant 7;\:a_{ij}=\left\{\begin{matrix} 1 & if <j, i>\in\!E,\\ 0 & if <j, i>\notin\!E. \end{matrix}\right.$

注意，这里列下标代表的是有向边的起点，行下标代表的是有向边的终点。至于为什么要这么设定（通常都应该是行代表起点、列代表终点的，这里反过来了），看到后面自然就会明白。

回到Example 3，我们给向量x变个形，变形规则是： $x=(x_{0},x_{1},...,x_{7})\rightarrow x'=(x_{0}^{'},x_{1}^{'},...,x_{7}^{'}):\;x_{i}^{'}=\left\{\begin{matrix} 0 & if \;x_{i}=0,\\ vertices[i] & if\; x_{i}=1. \end{matrix}\right.$

经这样的变形后，Example 3中的矩阵向量乘法运算式就变成了下面这样：

Example 4: Mx‘=y‘

$\begin{bmatrix} 0& 0 & 1 & 0 &0 &0 &0 &0 \\ 1& 0& 0& 0& 0 & 0 & 0 &0 \\ 0& 0 & 0 & 1 & 0 & 0 & 0 &0 \\ 1& 0 & 0 & 0 & 0 & 0 & 0 &0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 &0 \\ 0& 1 & 0 & 1 & 0 & 0 & 0 & 0\\ 0& 0 & 1 & 0 & 0 & 0 & 0 &1 \\ 0& 0 & 0 & 1 & 0 & 0 & 0 &0 \end{bmatrix}\begin{bmatrix} 0\\ B\\ 0\\ D\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}=\begin{bmatrix} 0\\ 0\\ D\\ 0\\ B\\ B+D\\ 0\\ D \end{bmatrix}\;\;\begin{matrix} A\\ B\\ C\\ D\\ E\\ F\\ G\\ H \end{matrix}$

为什么要对向量做这么一个变形呢？有何用意？我们回顾一下第三部分所讲的层同步并行BFS算法：如果在图G中选择从顶点A开始层同步并行BFS过程，那么下一层遍历迭代的顶点就是B和D；接下来，以B和D为Current Frontier，进行下一轮扩展，其未访问的Neighbor有C、E、F、H（其中E、F是B的Neighbor，C、F、H是D的Neighbor，二者共有的Neighbor是F），于是C、E、F、H四个顶点成为Next Frontier。再看看Example 4中结果向量y'的各分量与一旁标注出来的顶点之间的对应关系，对照着想想前面层同步并行BFS树扩展的过程，是否就有所发现了？

图G以及Example 3、Example 4都是来源于GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs这篇论文中的举例，如图3所示。

图3
Reference: Hu, Y., Du, Y., Ustun, E., & Zhang, Z. (2021). GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs. *2021 IEEE/ACM International Conference On Computer Aided Design (ICCAD)*, (1–9).

这张图已经把层同步并行BFS和矩阵向量乘法之间的联系表现得足够清晰，前文的分析也已经有足够的提示性，因此，为了省时省力，笔者就不在这里给出这种联系的具体的数学语言表达了。

经进一步分析我们还可以发现，层同步并行BFS的两种方法（top-down和bottom-up）还分别与SpMSpV和SpMV具有对应关系。图3(a)右边表现的就是第一部分所讲的按行优先遍历的SpMV计算过程，左边就是这个SpMV所对应的BFS过程，可以看出是按bottom-up方法进行的（其实看右边矩阵中各行的红色箭头所指的方向，可以看作是“由‘行’指向‘列’”，而“行”代表的是终点dst vertex，“列”代表的是起点src vertex，“由终点指向起点”，正是child-to-parent的方向，也就判断出是按bottom-up方法遍历的了；后面对top-down方法也可以像这样去判断）；而图3(b)右边表现的是按列优先遍历的SpMSpV计算过程，左边就是与之对应的按top-down方法进行的BFS过程。这样对照着看下来，本文标题所说的“两种层同步并行BFS算法与SpMV和SpMSpV之间的对应关系”，便十分明晰了。

GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs中还提到：

BFS typically applies SpMSpV at the first few iterations when the frontier is small (i.e., the input vector has a high sparsity), and switches to SpMV as the frontier grows (i.e., the input vector becomes denser).

BFS通常在前几次迭代中当Frontier的规模较小（即输入向量较为稀疏）时应用SpMSpV算法，并在当Frontier的规模增长（即输入向量变得更为稠密）时切换到SpMV算法。

至此，我们知道了不仅两种层同步并行BFS算法与SpMV和SpMSpV之间存在着两个对应关系，它们之间还能随着遍历迭代中由量变引起的质变而相互转化。我们可以将四者之间的联系表现在下面的这张图中，并以此作为本文的结尾：