CN113794709A - 一种用于二值稀疏矩阵的混合编码方法 - Google Patents

Abstract

本发明公开了一种用于二值稀疏矩阵的混合编码方法，包括：计算稀疏矩阵缩略图中连通域的数量，初步确定可能的稠密区域具体位置；根据每个稠密区域中黑色像素数占该区域边界矩形的比值判定稠密区域形状；将稀疏矩阵缩略图映射回原稀疏矩阵中，根据每个特定形状的稠密区域子矩阵中非零元占比，得到最终的稠密子矩阵；对稠密子矩阵进行反向编码；将原稀疏矩阵中其余零散分布的非零元组成一个超稀疏剩余子矩阵，并采用正向编码；将反向编码和正向编码结果求和，得到稀疏矩阵向量乘计算结果。采用混合编码方式充分压缩了稀疏矩阵的行列索引，减少了主机端和设备端的数据传输量。

Description

一种用于二值稀疏矩阵的混合编码方法

技术领域

本发明涉及二值稀疏矩阵的高效压缩技术领域，更具体的说是涉及一种用于二值稀疏矩阵的混合编码方法。

背景技术

稀疏矩阵最直观的一个特点是其非零元的数量要远远小于其元素总数，通常认为非零元数和矩阵元素总数的比值小于等于0.05的矩阵为稀疏矩阵，但这并非是一个严格的标准。稀疏矩阵几乎涉及到所有的大型科学工程计算领域中，比如计算流体动力学、统计物理、电路模拟、土木和机械工程、金融规划、气候模拟、医学图像、信息检索等。

二值稀疏矩阵是稀疏矩阵的一个特例，其非零元的取值都为1，这类稀疏矩阵常见于图计算相关的应用中，另外在稀疏矩阵和稀疏矩阵相乘的结果矩阵非零元数量预测中也会经常涉及到二值稀疏矩阵的计算。充分利用其非零元取值都为1的特点，往往能得到比通用稀疏矩阵压缩格式更高的压缩比和更优的计算性能。

随着通用图形处理器GPGPU(General-purpose computing on GraphicsProcessing Units)的快速发展和广泛应用，越来越多的研究人员聚焦于稀疏矩阵计算在GPU平台上的性能优化。稀疏矩阵计算在GPU上的实现基本由三个步骤组成：一是将稀疏编码的矩阵从主机端(CPU内存)传输到设备端(GPU内存)，二是利用GPU众多的计算核心完成稀疏矩阵计算，三是将计算得到的结果从设备端传输到主机端。

(1)数据传输的时间开销

目前，许多研究人员提出的稀疏矩阵压缩格式都会考虑稀疏矩阵计算在并行体系架构上的任务划分和负载均衡问题，即目标导向为尽可能减少稀疏矩阵计算部分的时间开销。这是因为稀疏矩阵计算往往出现在稀疏线性系统的迭代求解中，且成百上千次迭代计算中的稀疏矩阵计算往往基于同一个稀疏的系数矩阵，所以从主机端到设备端的数据传输往往只需执行一次。但事实上，在基于GPU的稀疏矩阵计算中，主机端和设备端之间的数据传输的时间开销往往占据了总时间的较大比值。研究表明大部分稀疏矩阵对应的SpMV的数据传输时间占比分布在80％～100％之间，充分证明了压缩二值稀疏矩阵的必要性。

(2)二值稀疏矩阵的特点

二值稀疏矩阵所有非零元的取值均为1，但目前的稀疏矩阵编码格式中并没有充分利用这一特征。比如判断输入的稀疏矩阵为二值稀疏矩阵时，取消值数组从主机端到设备端的传输，或者利用该特征来压缩矩阵的行列索引数组。

因此，如何提供一种用于二值稀疏矩阵的混合编码方法是本领域技术人员亟需解决的问题。

发明内容

有鉴于此，本发明提供了一种用于二值稀疏矩阵的混合编码方法，大幅度减少了稀疏矩阵的数据存储开销和主机端与设备端之间的数据传输开销，这在基于GPU的稀疏矩阵计算中具有绝对性优势。另一方面，通过此编码方案，能够大大简化稀疏矩阵计算在GPU端的实现，减少浮点计算的复杂度，从而加速稀疏矩阵计算及其广泛的科学计算应用。

为了实现上述目的，本发明采用如下技术方案：

一种用于二值稀疏矩阵的混合编码方法，包括：

计算稀疏矩阵缩略图中连通域的数量，初步确定可能的稠密区域具体位置；

根据每个稠密区域中黑色像素数占该区域边界矩形的比值判定稠密区域形状；

将稀疏矩阵缩略图映射回原稀疏矩阵中，根据每个特定形状的稠密区域子矩阵中非零元占比，得到最终的稠密子矩阵；

对稠密子矩阵进行反向编码；

将原稀疏矩阵中其余零散分布的非零元组成一个超稀疏剩余子矩阵，并采用正向编码；

将反向编码和正向编码结果传输到GPU上，根据输入向量X分别完成稀疏矩阵向量乘之后求和，得到稀疏矩阵向量乘计算结果。

优选的，稠密区域形状包括对角块状、矩形和三角形。

优选的，稠密子矩阵包括对角块稠密子矩阵、矩形块稠密子矩阵和三角块稠密子矩阵。

优选的，对稠密子矩阵进行反向编码具体包括：

假设稠密子矩阵中全是非零元，计算稠密子矩阵中所有元素均为非零元时稀疏矩阵向量乘结果；

将稠密子矩阵中所有零元采用CSR或其它压缩格式存储，并计算零元看作-1后对应的向量乘结果；

将非零元向量乘结果与零元向量乘结果相加即为稠密子矩阵的编码结果。

优选的，矩形块稠密子矩阵具体反向编码方法包括：

记录矩形块稠密子矩阵左上角位置(x_r，y_r)，矩形块大小a_r×b_r，以及矩形块中每一零元的行列索引信息；

假设矩形块稠密子矩阵中所有元素均为非零元，根据矩形块的位置、大小信息，以及输入向量x，计算矩形块稠密子矩阵所有元素均为非零元时对应的向量乘结果y_rect，具体计算公式为：

其中，y_rect[i]表示向量y_rect的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值，是编程语言中的标记符号；

三角块稠密子矩阵具体反向编码方法包括：

记录三角块稠密子矩阵左上角位置(x_t，y_t)，三角形块大小a_t×b_t，以及三角形块中每一零元的行列索引信息；

假设三角块稠密子矩阵中所有元素均为非零元，根据三角形块的位置、大小信息，以及输入向量x，计算三角块稠密子矩阵所有元素均为非零元时对应的向量乘结果y_tri，具体计算公式为：

其中，y_tri[i]表示向量y_tri的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值；

对角块稠密子矩阵具体反向编码方法包括：

记录对角块稠密子矩阵左上角的位置(x_d，y_d)，对角块边界矩形的大小a_d，b_d，对角块的宽度w_d，对角块中每一零元的行列索引信息；

假设对角块稠密子矩阵中所有元素均为非零元，根据对角块的位置和对角块的宽度计算出对角块中每一条对角线的偏移量，其偏移量的取值范围为[x_d-(y_d+w_d)，(x_d+w_d)-y_d]，计算每一条对角线与输入向量x的相乘结果，计算公式为：

其中，

表示该对角块稠密子矩阵第j条稠密对角线与向量x的相乘结果，

表示向量

的第i个元素，x[i+j]表示向量x的第i+j个元素；

将每一条对角线与向量x的相乘结果规约求和：

其中，y_dia表示该对角块稠密子矩阵所有稠密对角线，y_dia[i]表示该对角块稠密子矩阵和向量x相乘的结果向量的第i个元素；

将对角块稠密子矩阵、矩形块稠密子矩阵和三角块稠密子矩阵中所有零元像素采CSR或其它压缩格式存储，计算零元对应的向量乘结果为y_zero；

所有稠密子矩阵的编码结果为：y_dense＝y_rect+y_tri+y_dia-y_zero。

优选的，正向编码采用CSR或其它压缩格式编码。

经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种用于二值稀疏矩阵的混合编码方法，提出了针对稠密子矩阵反向编码，稀疏子矩阵正向编码的混合编码方法。一方面充分压缩了稀疏矩阵的行列索引，减少了主机端和设备端的数据传输量；另一方面利用此编码方法可以将SpMV(稀疏矩阵向量乘)计算中的乘累加运算转换为最基本的加减运算，简化并减少了浮点计算，从而加快了核函数的执行。本发明可以应用到所有二值稀疏矩阵相关的应用中，用于加速其对应的稀疏计算。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1附图为本发明二值稀疏矩阵的混合编码方法整体流程图。

图2附图为不同形状子矩阵的边界信息，其中图2(a)表示矩形块稠密子矩阵边界信息，图2(b)表示三角块稠密子矩阵边界信息，图2(c)表示对角块稠密子矩阵边界信息。

图3附图为二值稀疏矩阵混合编码示例。

图4附图为稠密子矩阵与向量相乘的示例。

图5附图为稠密子矩阵的零元与向量相乘的示例。

图6附图为超稀疏剩余子矩阵与向量相乘的示例

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例公开了一种用于二值稀疏矩阵的混合编码方法，如图1所示，包括：

对于给定的稀疏矩阵，将其缩略图作为输入，通过计算稀疏矩阵缩略图中连通域的数量初步得出可能的稠密区域的具体位置；

通过每个稠密区域中黑色像素数占该区域边界矩形的比值判定稠密区域形状：包括对角块状、矩形、三角形；

将稀疏矩阵缩略图映射回原稀疏矩阵，通过稠密区域在缩略图中的位置，得出该区域对应的子矩阵在原始矩阵中的位置，判断稠密区域中非零元数占比得到最终的稠密子矩阵，包括对角块稠密子矩阵、矩形块稠密子矩阵和三角块稠密子矩阵；在本实施例中，具体通过判定稠密区域子矩阵中非零元数是否大于该区域50％的占比得到最终的稠密子矩阵，在实际应用中也可以根据具体情况设计其余数值的具体占比；

对稠密子矩阵进行反向编码，即存储零元(稠密区域中白色像素)的索引信息：假设稠密子矩阵中全是非零元，计算稠密子矩阵中所有元素均为非零元时稀疏矩阵向量乘结果；将稠密子矩阵中所有零元采用CSR(Compressed Sparse Row)或其它压缩格式存储，并计算零元看作-1后对应的向量乘结果；将非零元向量乘结果与零元向量乘结果相加即为稠密子矩阵的编码结果；

稀疏矩阵往往会有少量非零元离散地分布在这些特定的形状之外，将原稀疏矩阵中其余零散分布的非零元组成一个超稀疏剩余子矩阵，并采用正向编码，即存储非零元的行列索引信息；

将反向编码和正向编码结果求和，得到稀疏矩阵向量乘计算结果。

进一步优化技术方案，对稠密子矩阵基于反向编码的SpMV具体包括：

矩形块稠密子矩阵：记录矩形块稠密子矩阵左上角位置(x_r，y_r)，矩形块大小a_r×b_r，以及矩形块中每一零元的行列索引信息，如图2(a)所示；

假设矩形块稠密子矩阵中所有元素均为非零元，根据矩形块的位置、大小信息，以及输入向量x，计算矩形块稠密子矩阵所有元素均为非零元时对应的向量乘结果，具体计算公式为：

其中，y_rect[i]表示向量y_rect的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值；

三角块稠密子矩阵：

记录三角块稠密子矩阵左上角位置(x_t，y_t)，三角形块大小a_t×b_t，以及三角形块中每一零元的行列索引信息，如图2(b)所示；

假设三角块稠密子矩阵中所有元素均为非零元，根据三角形块的位置、大小信息，以及输入向量x，计算三角块稠密子矩阵所有元素均为非零元时对应的向量乘结果，具体计算公式为：

其中，y_tri[i]表示向量y_tri的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值；

对角块稠密子矩阵：

记录对角块稠密子矩阵左上角的位置(x_d，y_d)，对角块边界矩形的大小a_d×b_d，对角块的宽度w_d，对角块中每一零元的行列索引信息，如图2(c)所示；

假设对角块稠密子矩阵中所有元素均为非零元，根据对角块的位置和对角块的宽度计算出对角块中每一条对角线的偏移量，其偏移量的取值范围为[x_d-(y_d+w_d)，(x_d+w_d)-y_d]，计算每一条对角线与对应的向量x的相乘结果，计算公式为：

其中，

表示该对角块稠密子矩阵第j条稠密对角线与向量x的相乘结果，

表示向量

的第i个元素，x[i+j]表示向量x的第i+j个元素；

将每一条对角线与向量x的相乘结果规约求和：

其中，

表示该对角块稠密子矩阵第j条稠密对角线与向量x的相乘结果，

表示向量

的第i个元素，x[i+j]表示向量x的第i+j个元素。

将稠密子矩阵中的所有零元的行列索引信息采用COO或CSR格式存储，计算其和向量x的相乘结果y_zero，则y_dense＝y_rect+y_tri+y_dia-y_zero即为稠密区域的计算结果。

进一步优化技术方案，对于超稀疏剩余子矩阵，采用COO或CSR格式编码非零元的位置，即正向编码，采用基于COO或CSR格式的SpMV算法计算超稀疏剩余子矩阵对应的计算结果y_sparse。

将反向编码和正向编码结果求和，得到稀疏矩阵向量乘计算结果为：y＝y_dense+y_sparse。

图3展示了一个10×10的二值稀疏矩阵混和编码的示例。该矩阵由3个稠密块组成，分别为矩形块稠密子矩阵(包含10个非零元)、三角块稠密子矩阵(包含9个非零元)和对角块稠密子矩阵(包含15个非零元)，以及一个超稀疏剩余子矩阵(包含5个非零元)，共39个非零元。按照传统的稀疏矩阵压缩方法，即正向编码稀疏矩阵的非零元，则COO格式需要存储的数据量为39×2＝78,CSR格式的数据量为11+39＝50。而本发明提出的反向编码方法的数据量为13+10+10＝33，针对COO格式的压缩比为78/33＝2.36，针对CSR格式的压缩比为50/33＝1.52。事实上，随着矩阵规模的增大，反向编码稠密子矩阵的收益更加明显，从而能够得到更高的压缩比。

以图3中的二值稀疏矩阵为例，图4、图5、图6介绍了该稀疏矩阵基于反向编码的SpMV计算流程。图4展示三个不同形状的稠密子矩阵与向量x相乘。即假设三个稠密子矩阵中的元素均为非零元，计算对应的结果向量y_rect，y_tri，y_dia，三个结果向量的和记为y_one。图5展示了稠密子矩阵中的零元和向量x相乘，结果为y_zero，图6展示了超稀疏剩余子矩阵和向量x相乘，结果为y_sparse。最终的混合编码SpMV计算结果为y＝y_one-y_zero+y_sparse。

相比于现有的稀疏格式压缩编码方案，本发明主要针对于二值稀疏矩阵的行列索引压缩。通过挖掘稀疏中有稠密，稠密中有稀疏的非零元分布信息，本发明提出了混和编码的压缩方案。对于稠密子矩阵，即非零元多于零元的子矩阵，我们对零元进行编码从而保证最少的数据量；利用稠密子矩阵所表现出的矩形、三角形和对角形的形态特征，本发明提出了特征驱动的压缩方案。而对于稀疏子矩阵，仍旧采用现有的正向稀疏编码方案。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (6)

1.一种用于二值稀疏矩阵的混合编码方法，其特征在于，包括：

计算稀疏矩阵缩略图中连通域的数量，初步确定可能的稠密区域具体位置；

根据每个稠密区域中黑色像素数占该区域边界矩形的比值判定稠密区域形状；

将稀疏矩阵缩略图映射回原稀疏矩阵中，根据每个特定形状的稠密区域子矩阵中非零元占比，得到最终的稠密子矩阵；

对稠密子矩阵进行反向编码；

将原稀疏矩阵中其余零散分布的非零元组成一个超稀疏剩余子矩阵，并采用正向编码；

将反向编码和正向编码结果传输到GPU上，根据输入向量X分别完成稀疏矩阵向量乘之后求和，得到稀疏矩阵向量乘计算结果。

2.根据权利要求1所述的一种用于二值稀疏矩阵的混合编码方法，其特征在于，稠密区域形状包括对角块状、矩形和三角形。

3.根据权利要求1所述的一种用于二值稀疏矩阵的混合编码方法，其特征在于，稠密子矩阵包括对角块稠密子矩阵、矩形块稠密子矩阵和三角块稠密子矩阵。

4.根据权利要求3所述的一种用于二值稀疏矩阵的混合编码方法，其特征在于，对稠密子矩阵进行反向编码具体包括：

假设稠密子矩阵中全是非零元，计算稠密子矩阵中所有元素均为非零元时稀疏矩阵向量乘结果；

将稠密子矩阵中所有零元采用CSR格式存储，并计算将零元看作-1后对应的向量乘结果；

将非零元向量乘结果与零元向量乘结果相加即为稠密子矩阵的编码结果。

5.根据权利要求4所述的一种用于二值稀疏矩阵的混合编码方法，其特征在于，矩形块稠密子矩阵具体反向编码方法包括：

记录矩形块稠密子矩阵左上角位置(x_r，y_r)，矩形块大小a_r×b_r，以及矩形块中每一零元的行列索引信息；

假设矩形块稠密子矩阵中所有元素均为非零元，根据矩形块的位置、大小信息，以及输入向量x，计算矩形块稠密子矩阵所有元素均为非零元时对应的向量乘结果y_rect，具体计算公式为：

其中，y_rect[i]表示向量y_rect的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值；

三角块稠密子矩阵具体反向编码方法包括：

记录三角块稠密子矩阵左上角位置(x_t，y_t)，三角形块大小a_t×b_t，以及三角形块中每一零元的行列索引信息；

假设三角块稠密子矩阵中所有元素均为非零元，根据三角形块的位置、大小信息，以及输入向量x，计算三角块稠密子矩阵所有元素均为非零元时对应的向量乘结果y_tri，具体计算公式为：

其中，y_tri[i]表示向量y_tri的第i个元素，x[j]是向量x的第j个元素，i和j是向量的索引值；

对角块稠密子矩阵具体反向编码方法包括：

记录对角块稠密子矩阵左上角的位置(x_d，y_d)，对角块边界矩形的大小a_d×b_d，对角块的宽度w_d，对角块中每一零元的行列索引信息；

假设对角块稠密子矩阵中所有元素均为非零元，根据对角块的位置和对角块的宽度计算出对角块中每一条对角线的偏移量，其偏移量的取值范围为[x_d-(y_d+w_d)，(x_d+w_d)-y_d]，计算每一条对角线与输入向量x的相乘结果，计算公式为：

其中，

表示该对角块稠密子矩阵第j条稠密对角线与向量x的相乘结果，

表示向量

的第i个元素，x[i+j]表示向量x的第i+j个元素；

将每一条对角线与向量x的相乘结果规约求和：

其中，y_dia表示该对角块稠密子矩阵所有稠密对角线，y_dia[i]表示该对角块稠密子矩阵和向量x相乘的结果向量的第i个元素；

将对角块稠密子矩阵、矩形块稠密子矩阵和三角块稠密子矩阵中所有零元像素采用CSR格式存储，计算零元对应的向量乘结果为y_zero；

所有稠密子矩阵和向量x相乘的计算结果为：y_dense＝y_rect+y_tri+y_dia-y_zero。

6.根据权利要求1所述的一种用于二值稀疏矩阵的混合编码方法，其特征在于，正向编码采用CSR格式编码。

Images