稀疏矩阵计算中的结果矩阵大小预测

第一个问题是在真实执行SpGEMM之前结果矩阵C的非零元数量未知。稀疏矩阵的非零元数量是决定其内存开销的关键指标，而在SpGEMM中，结果矩阵的稀疏性在真实执行程序之前是未知的。因此，很难提前对结果矩阵进行精确的内存分配。致力于结果矩阵大小预测的研究主要可以分为以下四种方法。

（1）精确的方法：SpGEMM算法使用精确的内存分配方法解决结果矩阵的大小预测问题，通常包括两个阶段：符号阶段和数值阶段。通过符号阶段计算出准确的结果矩阵非零元数，而通过数值阶段得到SpGEMM最终的计算结果。

SpGEMM在KokkosKernels、cuSPARSE、MKL和RMerge中的实现是该方法的典型代表。此外，Demouth、Nagasaka等人和Demirci等人也采用了这种方法。Deveci等人设计了一种类似于颜色压缩思想的图压缩技术，通过将矩阵B的列打包为二进制位来压缩矩阵B，以实现对符号阶段的加速。

提前计算出准确的结果矩阵非零元数不仅可以减少内存开销，还可以在具有相同稀疏结构的输入矩阵参与的不同SpGEMM中重用结果矩阵C的稀疏结构。而且，该方法在图分析中具有显著的优势，因为许多图分析的计算只需要符号阶段的执行，便可得到最终的结果。但是，两个阶段的计算同时也意味着它需要对输入矩阵进行两次迭代计算，相比于以下方法具有更高的计算开销。

（2）概率方法：通过对输入矩阵进行随机采样和概率分析来估计结果矩阵的稀疏结构。

Cohen 认为，可以将估计结果矩阵中非零元数量的问题转换为估计有向图中可达集的大小的问题。然后，他提出了一种基于蒙特卡洛的算法，该算法以高置信度和较小的相对误差来估计可达集的大小。因为了解稀疏矩阵非零元的结构不仅有助于为输出矩阵有效分配内存空间，而且在计算三个或更多矩阵的乘积时可以用来确定矩阵相乘的最佳顺序，其中最佳顺序是指满足总计算数最小化的矩阵相乘顺序。因此，Cohen 基于先前的可达集大小预测算法提出了一种结果矩阵非零元结构的预测方法，该方法是一种线性时间复杂度的算法。后来，Amossen等人将该方法改进为预期的线性时间复杂度O（n），其中相对误差满足（其中z为结果矩阵中非零元的数量）。Pagh等人还提出了一种基于输入矩阵的非零元来估计结果矩阵的列/行大小的方法。

基于概率的预测方法不需要执行类似于SpGEMM的计算，因此该方法的时间成本较低。但是，考虑到此种方法通常给出的预测结果为近似估计，因此预测失败时就必须分配额外的内存空间。

（3）上限方法：通过计算结果矩阵非零元的上限来分配相应的内存空间。最常用的计算结果矩阵C的上限的方法是对于矩阵A中每一个非零元，计算矩阵B中相应行的非零元数。假设输入矩阵A和B均采用CSR格式存储，I和J分别是行指针数组和列索引数组，算法4−4给出了用于计算结果矩阵每一行中非零元数量上限的伪代码，计算结果存储在数组U={u0，u1，…，up−1}中。(www.xing528.com)

Bell等人提出的经典的ESC（扩展，排序，压缩）方法就采用了上限法来为结果矩阵分配存储空间，且这种方法在开源的稀疏线性代数库CUSP中得到了应用。Nagesaka等人还计算了结果矩阵每行中标量非零元相乘数量的最大值。然后，每个线程根据该最大值来创建相应的哈希表，在整个计算过程中通过重新初始化每行对应的哈希表来实现对哈希表的重用。

基于上限的方法高效且易于实现，但因为该方法忽略了计算过程中相同索引项的合并，所以会有过度分配内存空间的缺点。特别是对于某些特定稀疏模式的矩阵，使用上限法估计的结果矩阵的大小有很大误差。

（4）渐进法：根据实时需要来动态分配内存。该方法首先分配适当大小的内存空间，然后进行SpGEMM计算。计算过程中如果出现内存空间不足，则重新分配更大的内存。

Matlab中SpGEMM的实现就采用了这种方法。它首先粗略地猜测结果矩阵的大小，并申请相应大小的内存空间。在某个需要更多内存空间的时刻申请一块新的更大的内存，新申请的内存块的大小往往是之前内存大小的1.5倍。同时，还需要将已经计算出的列拷贝到新的内存块中。

首先，如果此种方法首次猜测成功，则SpGEMM的计算就不会有其他的时间开销。但是，如果第一次估算失败，则这种方法的效率极低，尤其在GPU中新内存的申请和数据的拷贝非常耗时。此外，在线程级并行的SpGEMM中，如果所有线程使用的是各自的局部内存，那么还需要昂贵的时间成本将所有独立的局部内存汇总到共享的结果矩阵中。

此外，还有研究人员采用混合的预测方法来解决该挑战性问题。刘伟峰教授和Brian Vinter提出了一种混合的结果矩阵预分配方法。他们计算结果矩阵每一行中非零元数的上限，并根据该上限将所有行分为多个组。然后，他们采用上限法为短行（即非零元数的上限较小的行，反之即为长行）分配内存，采用渐进法为长行分配内存。