Level-Based Blocking for Sparse Matrices ：Sparse Matrix-Power-Vector Multiplication 论文笔记

相关知识

稀疏矩阵向量乘（SpMV）和矩阵幂次（MPK）

稀疏矩阵：稀疏矩阵是指具有大部分元素为零的矩阵。相对于稠密矩阵，稀疏矩阵只存储那些非零元素，而将零元素省略掉，从而节约了存储空间。本文使用压缩行存储（CSR）格式。

CSR包含三个数组（V，C，R）：

1. Value，用来存储矩阵中非零元素的值
2. COL_INDEX，第i个元素记录了V[i]元素的列数
3. ROW_INDEX, 第i个元素记录了前i-1行包含的非零元素的数量。

例：

​ 假如我们想得到“30”的坐标位置，通过V得到index为2，通过COL_INDEX得到列数为1；

​ 又因为index=2∈[ROW_COL[1]=2,ROW_COL[2]=4)，所以行数为1。

稀疏矩阵向量乘(SpMV)是科学、工程和数据分析中使用的各种计算算法的关键构建块。但SpMV内核在现代计算设备上的性能很差，主要原因如下：

1. 非常规的内存访问（irregular memory accesses）；非零元素的位置是随机的
2. 间接的内存查询（indirect memory referencing）；CSR的存储格式限制，查询是间接的
3. 低算术强度（low arithmetic intensity）；计算受限于内存的读取，而不需要大量的算力。

MPK（matrix power kernel）被如下定义：对于一个给定的稀疏矩阵A和一个稠密矩阵x，计算(A^p)*x直到Pm（P= 1，……，Pm）最大，并存储所有结果进行后续计算。我们定义Y(p) = (A^p)*x，且Y(0)=x。

本文不进行SpMv算法的优化，但将视角看向了稀疏矩阵幂次（MPK）中多次SpMv算法调用的并行优化问题。通过扩展了递归的染色算法（RACE）框架，解决了多个SpMv调用之间的依赖关系。

设置一个MPK实验基准

使用标准的SpMv实现或库调用，来计算Pm序列。并以其最小值为基准与本文提出的优化MPK算法作比较。

以下是基于CSR的MPK算法计算A^(Pm)*x。这个代码用C++复现应当不难。

double :: val[Nnz]	#存储A的非零元素值
int :: col [Nnz],rowPtr[Nr+1]	#列索引和行指针
#参考CSR存储格式
double :: y[Nr,0:pm]	#存储输入和输出向量，y(0)是输入向量，y(pm)是输出向量

for p = 1:pm do
	y[:,p] = SpMV(y[:,p-1],0,Nr-1)
	#计算y(p) = A * y（p-1）
end for


function SpMV(double :: in_rhs[Nr], int :: row_s,int :: row_e)
	double :: out_lhs[Nr]
	#以下算法使用OpenMP静态并行处理
	for rouw=row_s:row_e do
		double :: tmp = 0
		for idx=rowPtr[row]:(rowPtr[row＋1]-1) do	#遍历第idx行的所有非零元素
			tmp＋= val[idac] * in_rhs[col[idc]]		#计算改行非零元素与输入向量的乘积
		end for
		out_lhs[row] = tmp							#更新向量该行的值
    end for
	return out_lhs
end function

基准MPK算法存储Pm+1个向量{y(p)} ，然后连续调用Pm进行SpMv算法，如果缓存太小以至于装不下整个矩阵，就必须从主存中读取Pm，就会造成内存限制性能。该基准的最小主存平衡即Bc = 6 byte/flop，Bc是用来衡量主存与代码平衡的指标，是用数据流量除以最小的浮点运算量。改值越小说明没进行相同的计算量下访问更少的主存。

MPK执行时的分块策略

由于重复应用相同的稀疏矩阵A，通过重用缓存中的矩阵块来进行多次幂次的连续计算，从而减少数据传输，从而有很大的性能优化潜力。

MPK算法优化的基础思想在于：为适合缓存的A的部分块计算SpMv，然后在其子集上继续计算下一个幂集的结果，这样就在连续更新的过程中确定了空间结构的依赖。为了演示该依赖关系，使用一个简单的带状稀疏矩阵来进行演示。

(a)是一个简单带状稀疏矩阵，其中黄色区域为适合缓存的部分块，用以计算Ax的相应部分。

(b)当我们计算下一个幂级的部分，即(A^2)x= A*(Ax)的蓝色区域，为了保证多级SpMv步骤之间具有依赖关系，蓝色区域必须是黄色区域的子集。

推广到任意稀疏矩阵，I为A矩阵是行索引的部分子集（对应带状矩阵的黄色区域），C(I)为A矩阵I行的非零元素的列索引集合。

当进行完当前幂级的计算后，即：

计算下一幂级，可以应用下一个SpMv操作在K集合上，K为I的子集。

对于这种方法，I集合的选择对于性能的影响是决定性的：

1. 对于I和C(I)相关联的矩阵元素必须适合缓存。
2. 必须存储为高空间局部性和时间局部性的结构，而无需间接访问
3. MPK涉及的I行索引和C(I)列索引要尽可能相近，我的理解（对于行多，列少的矩阵，可以更好的并行计算MPK）。

解决这些问题的是将SpMv操作视为RACE着色方案中做的图遍历问题BFS。

基于等级分块（LB）的MPK

该方案的图论定义

• G = （V，E），V(G)代表顶点集合，E (G)代表边集合，对于一个稀疏矩阵，V代表所有的行/列索引（对称矩阵，无向图，该方法也适用于非对称矩阵），对于E则是如下定义：
  • 
• 邻域。顶点u的邻域（与u有边相邻的顶点集合）是N(u)，定义如下：
  • 
• 子图H，包含G部分顶点，与这些顶点相对应的边。

现在在这个图的定义下，当前幂级的计算，即更新u∈V(G)的值：

理解：对于一个顶点u，即第u行，所有u 的邻域即第u行第v列有非零元素，与向量v行相乘再相加就得到了下一幂级的该顶点值。

使用一个简单带状稀疏矩阵来说明上述图的结构：

等级

基于BFS进行等级划分，流程如下：

这种排列增加了相邻顶点之间的数据局部性，因此这种排列被广泛用于基于SpMv的算法预处理步骤。

处理结束后，我们可以进一步划分基于等级的邻域划分，即：

这个公式定义了不同幂级下各个等级之间的SpMv计算的依赖关系。

基于等级分块（LB）的MPK

接下来，我们引入LP图来可视化MPK计算中各等级分块之间的依赖关系，LP图中x轴为等级分块（0-14），y轴为幂级（1-5）。

1. 对角线内，即i+p = const，都从下到上遍历
2. 对角线间，从左到右遍历

比如计算p=4时L(6)的值，需要先计算p=3时L(5)、L(6)、L(7)的值，如上图的第38次计算需要用到第27、32、37次的计算结果。

对于上述算法的缓存局部性由所给等级的重用距离来表示，在LP图中即该次计算在下一次何时才会重用（即同一级别上的两个计算之间的执行步数）。

由于幂级的数量级比等级的数量级小得多，所以可以假设最大重用距离为pm+1（最大幂级+1），也就是两次计算之间的数据都可以由缓存提供，除了第一个幂级需要访问主存，之后的计算都在缓存内进行，这将速度提升了一倍。

算法实现

for d=1:Lm+Pm-1 do #遍历每一个对角线 d=i+p
	#以上LG图为例
	P_start = max(1,d-(Lm-1))	#特殊情况，例65、69、72、74
	P_end = min(d,Pm)			#特殊情况，例0、2、5、9
	#其他都是1-5
	#计算对角线从右下到左上
	for P=P_start:P_end do #遍历每一个对角线上的顶点
		#以38点为例，d=10，p=4
		i = (d - p) #因为 d=i+p
		#i=6
		y[p] = SpMv(y[p-1],level_ptr[i],level_ptr[i+1]-1)
		#level_ptr存储了每个等级的第一个顶点条目
		#y[p-1]即上一幂级的运算结果，对于38这个顶点来说就是27、32、37
		#按对角线运算次序来说，这三个顶点一定在计算38之前都运算结束了
	end for
end for

OpenMP并行化是在SpMv函数中使用静态调度完成的。

初始版本的性能分析

对于上述这个版本，已经有了一定的性能改进，但是它往往达不到预期的性能，有百分之五十的矩阵加速不到百分之十，有十个矩阵的性能甚至下降。

主要分析两个方面：

1. 多线程性能。
2. 主存代码平衡。

对于矩阵pwtk，我们可以看到经典的内存带宽饱和模式的性能增长曲线（橙色线和蓝色线）。这两个曲线均达到了Baseline最小内存代码平衡即Bc=6 byte/flop。

注：算第一幂级是需要访问主存，并将计算结果放入缓存中，所以这一部分是应当和基准实验性能相当。

当Pm=4，LB算法的Bc只有1.5 byte/flop，符合预期，优化了接近四倍。但是问题是尽管优化了主存代码平衡，但性能依旧有评价，性能并没有随着核心的增加而上升，甚至出现了倒退的情况。

通过进一步分析，与基准实验对比，LB MPK(Pm=4)的指令在OpenMP的内存屏障中（用于线程同步），有大部分在自旋等待循环中执行，这表明当核心数增加以后，线程之间的同步是一个潜在的性能瓶颈。一个等级的工作负载太低，因此无法忽略同步成本。

对于矩阵Flan_1565显示出相反的特征。LB MPK(Pm=1)与基准实验结果一致，LB MPK(Pm=4)则是实现了1.2倍的性能提升，从主存代码平衡Bc可以看出，基准与LB MPK(Pm=4)也基本上是1.2倍的性能提升。这时候，限制性能的是矩阵本身。在划分等级的时候，有的等级可能包含20000行（每行大概75个非零），而有的可能很小，这就导致了LB MPK(Pm=4)在核心增多的情况下会出现不规则的性能收缩。

基于以上两个问题，本文提出了三个改进方法，前两个目标是通过形成更大的“等级组”并通过点对点同步代替昂贵的屏障来减低同步成本。第三种优化方法使用递归地分割这些矩阵来提高缓存利用效率，从而提高矩阵中具有主导或最为庞大的部分的性能。

等级组（LG）

等级组的思想是：连续的级别被聚合成所谓的等级组。并行执行在每个组内进行，并且所有线程在每个组计算之后进行同步。

由上文可知，如果我们假设等级组的大小相似，那么缓存必须能装下相邻的等级组，否则就会需要向主存寻找数据。如下：

Nnz(T(i))表示第i个等级组T(i)所具有的非零元素数，(Pm+1表示最大重用距离，12 bytes为设定的每个非零条目的大小。）

C为缓存大小，f为安全稀疏（设置为0.5，用以考虑来自数据结构的额外流量和缓存替换策略的低效率）

该限制下，T(i)的设置以i=0为例：

j为最大使得该等式成立的值。

明显，该过程使得等级组的数量远小于等级的数量，大大减少了同步的次数。同时由于缓存限制条件的存在，主存代码平衡仅略有增加。

点对点同步

下图展示等级组的划分和同步过程：

图a展示了矩阵图在进行等级分组后的分组情况。可以看到每个等级组具有的定点数相差不大。

图b的Lp图展示了五个等级组的执行次序以及依赖关系。我们可以在这张图中看出同步条件：

1. 该顶点的正下邻居顶点执行结束
2. 该顶点的左下邻居顶点执行结束
3. 该顶点的右下邻居顶点执行结束

需要注意的是，当满足1、2条件时，则3条件一定满足。即图b的两个红色箭头执行完成时，蓝色箭头一定执行结束。例如对于16顶点来说，11和15执行结束时，7一定执行结束，因为在12顶点计算的过程中需要7顶点计算结束。

不一定正确的想法：从理论上说，当12顶点执行结束时，7和11顶点应当都执行结束了，因此16顶点只需要与15顶点同步即可。我的理解是只需要考虑3条件即可，但下述方案考虑了1、3两个条件。猜测可能是T(0)的顶点没有左上父顶点，只有条件1能约束。

为了满足条件1，使用一个易失性整型数组U来进行同步，U的大小等于等级组数目，Nt为线程数。初始化时U[x]为0。当x等级组的一个线程完成了一次工作时，U[x]自动加一。检查y等级组是否完成了p-1幂级的计算即检查是否U[y] = (p-1) x Nt（Nt个线程全部完成一次工作即表示某幂级全部工作进行结束）。如果不满足则需要其他线程自旋等待直到其他线程完成它们的计算。

为了满足条件3，使用一个易失性整型数组V来进行同步，V的大小为Lp图的维度，即Lm x Pm。这里只有在两个相邻级别组的边界上工作的线程需要交互。例如图d，标红圈的线程为需要同步的线程，在T(3)上的1、2线程需要等待T(4)上的0线程工作完毕才能开始工作。为了实现这一点，需要检查下一等级组的计算p-1幂级的线程是否完成计算。设置V(x)(p)，初始值为0，当T(x)上工作的线程在执行完p幂级的计算后自加1，。若V[y+1][p-1] = H[y+1]则满足条件3（H是预先计算的在每个等级组第一层工作的线程数，该等式表明相邻的线程全部计算完成）。

实验证明这种同步方法将性能提升了1.2倍。图见性能分析。

递归

之前通过分析Flan_1565矩阵，我们可以发现如果等级组过于庞大，就会导致并行性能不规则抖动。因此提出了一种基于RACE染色算法的递归算法，通过递归的方式连续生成新等级组或者缩小等级组，直到他们适合缓存。

首先，需要定位不适合缓存的连续等级，并分离这些等级形成的子图。然后通过BFS遍历得到一组新的等级组，这些新的等级组往往比原先等级组小很多。

如图：

对于a图，假设认为超过六行的等级组是不适合缓存的，那么，将这些顶点做构成的子图分离开，得到一个子图。

然后对这个子图进行等级组划分的算法。最终递归进行这个过程，直到每个等级组都满足缓存限制条件。

下面使用s来表示当前的递归次数，sm表示最大递归次数。我们是应当对sm进行一定的限制，因为该递归算法会使得子图的数据局部性丢失，因为子图是内部的BFS，没有考虑和子图外的邻居。

对于每个递归创建的子图都划分了自己的等级组，因此每个递归阶段都应当有自己的Lp图。

这里展示pm=2时s=0到s=1的计算过程。首先划分出递归等级组的区域，即黄色区域。子图共有四个等级组两个幂级共八个顶点放入s=1的Lp子图中，按Lp顶点标号为次序进行计算。

当pm变大，这个过程就会变得复杂。如上图：这是pm=3（上）和pm=5（下）时的Lp图，其中6-8的等级组需要递归计算（矩形框区域）。

可以看出在相关区域（平行四边形区域），有三种依赖关系：

1. 只需要6-8三个等级组提供输出的关系，蓝色区域
2. 只需要给6-8三个等级组提供输入的关系，橙色区域
3. 两个关系都有，即菱形区域

红色箭头表示与递归边界点之间最长的输入/输出关系。

在进行子图递归的时候，也需要将菱形区域内的边界点进行递归计算。也就是对于T(4)和T(5)以及T(9)和T(10)都要参与递归划分和计算。将第三种依赖关系加入到递归处理的好处是，子图外的顶点计算顺序与依赖关系都不用改变。

同步过程有两点改变：

1. 首先，每个递归过程都需呀单独定义和使用U、V、H三个独立的数组来同步子图内部相应递归阶段的计算。
2. 依赖关系1和2分别需要执行在子图计算之后和之前。这可以通过检查U数组的值来实现

通过递归计算的方法，相对于没有使用递归的方法，性能提升了1.4倍。通过计算Bc可以看出，其主存代码平衡减少到了一半。这些数据表明该方法对于主存的访问不再是瓶颈。

性能分析

实验设置了三个平台，Intel Cascade Lake(CLX)、Intel Ice Lake(ICL)、 Intel Ice Lake(ROME)。实验结果如下图:

x轴表示矩阵对应的名称，并且按矩阵大小从左向右排序。虚线表示各自硬件的总缓存大小，可以看到在虚线左侧的矩阵在Baseline上性能好，这是因为矩阵能全部放入缓存中，不受主存读取速度的限制。而当缓存不足以装下整个矩阵时，RACE算法的优越性就体现出来了。

对于a、d两个图，可以看到在最大的缓存内矩阵，RACE算法就已经超过Baseline了。这里有两个原因：

1. 与AMD ROME相比，英特尔架构上L3和主存带宽之间的过渡是渐进的，即使工作集超过缓存大小也会存在较明显的缓存效应。如下图，CS表示cache size：

   

2. 对于a图缓存内矩阵的性能都超过了Baseline，而g图缓存内性能都不如Baseline。这里一方面是L2和L3缓存在英特尔架构上几乎具有相似的大小，而RACE的分块主要针对L3和L2缓存的组合，因此对于适合L3缓存的较小矩阵，与Baseline相比，可以减少L2的流量。另一方面对于ROME来说，L3缓存比L2缓存大得多，因此分块只在L3缓存中执行，因此对于L3缓存中的矩阵匹配没有任何好处。

   ps：不懂在说什么胡话。。。

对于内存驻留的情况，RACE在所有架构上都有这压倒性的优势，与Baseline相比，可以实现2-5倍的性能提升。这与b、e、h三张图展示的L2、L3和主存的流量情况相关。

可以看出，在三个测试平台上，大多数情况下减少了主存的流量，并小于Baseline的最小内存代码平衡即Bc=6 byte/flop。因为主存流量的减少，RACE在ROME上实现了最高的速度，因为较大的L2和L3缓存使得对于主存的流量显著小于英特尔平台。

我们观察到在ROME上平均加速为3.2倍，最大5倍。而在CLX和ICL上，平均加速1.6倍和1.9倍，最大加速为3.2倍和3倍。

预处理成本

已经分析完了RACE的性能，现在研究它的预处理开销。c、f、i中的图像显示了RACE对内存驻留情况下的预处理成本统计。预处理所需要的时间等效于SpMv操作的数量。

在ROME上，对于大多数矩阵，预处理时间远低于40个SpMv的等效时间。而在英特尔系统上却要多不少。这是由于较小的缓存需要更多递归分组和计算。大部分预处理时间（大于百分之九十五）花在BFS的确定等级上。

总结和展望

在本文中，我们开发了一种基于级别的分块算法(RACE MPK)来提高稀疏矩阵幂向量乘法(MPK)的性能。RACE算法使用由BFS生成的等级，通过在连续的幂级计算中重用的矩阵条目来增加矩阵条目的时间局部性。引入了各种面向硬件的算法优化策略，如等级组LG、点对点同步和递归的分块方案，以进一步提高性能。最后对42个代表矩阵进行全面的性能分析表明，本文方法在现代Intel和AMD的cpu上的性能比标准MPK实现平均高出2倍和3.5倍。