一种外积累加求解稀疏矩阵与稠密矩阵内积的方法

技术领域

本发明属于计算机应用技术领域，尤其涉及一种外积累加求解稀疏矩阵与稠密矩阵内积的方法。

背景技术

矩阵内积乘法是数值计算、机器学习、工程模拟等众多领域算法的基础，而稀疏矩阵A(MxK)与稠密矩阵B(KxN)内积乘法是一类常见的计算场景。稀疏矩阵A(MxK)与稠密矩阵B(KxN)，其特点在于稀疏矩阵A其中存在一定量的零元素。这部分零元素是可以不参与计算依然保证运算结果的正确性，同时对于零元素可以直接放弃存储来提升内存的利用效率。

但是在计算和存储的过程中对稀疏矩阵A之中零元素的忽略导致内存访问的不连续，不利于利用现代处理器架构中的各级高速缓存。因而涉及稀疏矩阵内积的计算效率低于稠密矩阵A与稠密矩阵B的内积。

此处以现有常规方法外积累加求解稀疏矩阵A(4x4)与稠密矩阵B(4x3)的内积举例，稀疏矩阵A(4x4)是一般稀疏矩阵，其特点是零元素的分布没有规律。图1是稀疏矩阵A(4x4)与稠密矩阵B(4x3)内积计算按照外积之和累加方式求解的第一步。本步中首先取出稀疏矩阵A(4x4)的第一列a1(4x1)与稠密矩阵B(4x3)第一行b1(1x3)执行外积运算得到结果c1，以此类推稀疏矩阵A的第2列a2与周密矩阵B的第2行b2也执行外积运算得到结果c2。稀疏矩阵A(4x4)有四列，稠密矩阵B(4x3)有四行，一一对应外积得到c1(4x3)、c2(4x3)、c3(4x3)、c4(4x3)。c1、c2、c3、c4四个外积结果对应元素相加即得到内积结果。而在一般的矩阵内积乘法中，以稀疏矩阵A(4x4)的第i行内积稠密矩阵B(4x3)第j列得到目标矩阵(i,j)位置的值。完成整个内积的计算，稠密矩阵B(4x3)需要被稀疏矩阵A(4x4)各行重复读取计算四次。

图2将稀疏矩阵A(4x4)第一列a1(4x1)与稠密矩阵B(4x3)第一行b1(1x3)外积计算逻辑做个说明。列向量a1(4x1)与行向量b1(1x3)外积结果构成4行3列的矩阵c1。其计算的过程等价于取列向量a1的各个元素分别乘以行向量b1。稀疏矩阵A(4x4)第i列ai(4x1)与稠密矩阵B(4x3)第i行bi(1x3)外积结果ci的计算同c1计算的逻辑。

图3说明了上述外积结果c1、c2、c3、c4四个矩阵相加的过程，矩阵相加遵循数学中对应元素相加的过程。如此也就求解出了稀疏矩阵A(4x4)与稠密矩阵B(4x3)内积结果。可以看到上述外积累加求解稀疏矩阵A(4x4)与稠密矩阵B(4x3)内积过程中，并没有专门考虑稀疏矩阵A(4x4)的零元素与稠密矩阵B的某一行外积相乘后，生成一个零向量。该零向量若是不计算，同时也不参与后面累加过程，不影响内积结果的正确性。但是在现有的外积累加计算中，稀疏矩阵中的零元素参与了外积计算，其不仅降低了外积运算的计算效率也增加了外积计算后零值向内存写入的需求，进而减少了内存写入的有效带宽。

发明内容

本发明针对上述的稀疏矩阵A(MxK)与稠密矩阵B(KxN)内积乘法所存在的技术问题，提出一种方法简单且能够有效避免无效运算、节省内存访问的一种外积累加求解稀疏矩阵与稠密矩阵内积的方法。

为了达到上述目的，本发明采用的技术方案为，一种外积累加求解稀疏矩阵与稠密矩阵内积的方法，包括以下有效步骤：

a、首先完成对稀疏矩阵A(MxK)与稠密矩阵B(KxN)原始数据的读取；

b、对稀疏矩阵A(MxK)进行列向数据压缩，将二维布局的稀疏矩阵A(MxK)列向压缩为一维布局的列向压缩稀疏矩阵A(MxK)；

c、读取列向压缩稀疏矩阵A(MxK)第i列以及读取稠密矩阵B(KxN)第i行外积运算生成ci；

d、对所得到的外积结果进行累加，即可得到稀疏矩阵A(MxK)与稠密矩阵B(KxN)的内积。

作为优选，还包括e步骤，e、判断是否完成到稀疏矩阵A(MxK)的最后一列，如果是，则完成计算；如果否，则i自增1，重复c～e步骤。

作为优选，所述b步骤中，列向压缩稀疏矩阵A(MxK)包括列起点偏移量数组、行索引数组以及非零元素值数组，其中，所述列起点偏移量数组用于表示稀疏矩阵A(MxK)每一列中第一个非零元素在非零元素值一维存储结构中相对起始地址的偏移量以及稀疏矩阵A(MxK)转成一维连续存储方式所占据的内存空间的长度的组合；所述行索引数组用于表示原非零元素值数所在稀疏矩阵A(MxK)的行索引值；所述非零元素值数组行用于按照列顺序依次填入非零元素数值。

与现有技术相比，本发明的优点和积极效果在于，

1、本发明提供一种外积累加求解稀疏矩阵与稠密矩阵内积的方法，通过将原有的二维布局的稀疏矩阵A(MxK)列向压缩为一维布局的列向压缩稀疏矩阵A(MxK)，除去了稀疏矩阵A(MxK)的零元素值，进而避免了无效数值的计算，达到节省内存访问和无效运算的目的，从而提高了运算效率，为数值计算、机器学习、工程模拟相关技术的快速运算提供保障。

附图说明

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

图1，稀疏矩阵A(4x4)与稠密矩阵B(4x3)的外积累加求解内积乘法示例；

图2，稀疏矩阵A(4x4)第一列与稠密矩阵B(4x3)第一行执行外积计算示例；

图3，稀疏矩阵A(4x4)与稠密矩阵B(4x3)的外积结果累加示例；

图4，稀疏矩阵A(4x4)依据列压缩存储方式存储的转换过程；

图5，稀疏矩阵A(4x4)CSC压缩存储列起点偏移量数组逻辑示例；

图6，稀疏矩阵A(4x4)CSC压缩存储行索引数组逻辑示例；

图7，稀疏矩阵A(4x4)列向压缩形成一维数组与稠密矩阵B(4x3)外积运算选取的数据示例；

图8，稀疏矩阵A(4x4)列向压缩形成一维数组与稠密矩阵B(4x3)外积累加过程示例；

图9，列向压缩稀疏矩阵A(MxK)与稠密矩阵B(KxN)外积累加流程图。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和实施例对本发明做进一步说明。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互组合。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用不同于在此描述的其他方式来实施，因此，本发明并不限于下面公开说明书的具体实施例的限制。

实施例1，本实施例提供一种外积累加求解稀疏矩阵与稠密矩阵内积的方法

和背景技术一样，本实施例用稀疏矩阵A(4x4)与稠密矩阵B(4x3)来描述本实施例所提供的技术方案，以便更好的理解本发明所提供的技术。

首先完成对原始数据的稀疏矩阵A(4x4)和稠密矩阵B(4x3)读取，其具体数值为图1所示，在此处，对于稀疏矩阵A(4x4)按照图4中401所展示的矩阵形式进行表示，零元素则用黑色方块，行索引与列索引本例从0开始，从1开始亦可，不影响本方法的实施。

然后对稀疏矩阵A(4x4)进行列压缩，去除其中的零元素，如402所示，402中包括列起点偏移量数组403、行索引数组404和非零元素值数组405三个成分构成稀疏矩阵A(4x4)列压缩存储的所有成分，内存地址406只是用来辅助解释所用，实际算法实现并不需要。400中虚线箭头用来指示401中的非零元素经过列向压缩存储后，在402中映射的位置。

列起点偏移量数组403表示用于表示稀疏矩阵A(4x4)每一列中第一个非零元素在非零元素值一维存储结构中相对起始地址的偏移量。非零元素值数组405中在连续的内存地址上依次存储着稀疏矩阵A(4x4)中的非零元素值。行索引数组404某个地址偏移量位置addrN存储的数值，行索引数组用于表示原非零元素值数所在稀疏矩阵A(MxK)的行数。比如，a31处于第一列第三行，那么在行索引数组404中，a31的上方显示的数值，如果从0开始记行数则为2，从1开始记行数则为3。在图4中，采用从0开始记行数，为此，其上方的数值为2。

列起点偏移量数组中除最好一位的数值表示为稀疏矩阵A(MxK)每一列中第一个非零元素在非零元素值一维存储结构中相对起始地址的偏移量，而列起点偏移量数组的最后一位数值则标志着稀疏矩阵A(MxK)转成一维连续存储方式所占据的内存空间的长度，在图5中，用cvalue4来表示，这样，列起点偏移量数组503的第一个地址位置addr0存储的数值是cvalue0，本例中cvalue0实际数值是0，即是第一列的第一个非零元素相对非零元素值数据505的起始地址偏移0位；列起点偏移量数组503的第二个地址位置addr1存储的数值是cvalue1，本例中cvalue1实际数值是3，即是第二列的第一个非零元素相对非零元素值数据505的起始地址偏移3位；依次类推，第三列第一个非零元素相对非零元素值数据505的起始地址偏移5位，第四列第一个非零元素相对非零元素值数据505的起始地址偏移6位，而非零元素值数组的所有非零元素数值的个数为9，那么cvalue4就等于9，所以列起点偏移量数组503的数值就0、3、5、6、9，当然，此类算法是基于列起点偏移量数组503起始地址以零为起点时，如果基于1为起点，其数值就是1、4、6、7、10。

总结起来说，对于稀疏矩阵A(MxK)，其第一列的所有非零元素在非零元素值数组505中所处的地址范围是：大于等于cvalue0，小于cvalue1；第二列的所有非零元素在非零元素值数组505中所处的地址范围是：大于等于cvalue1，小于cvalue2；第K列的所有非零元素在非零元素值数组505中所处的地址范围是：大于等于cvalue(k-1)，小于cvaluek。这样，就由起点偏移量数组以及行索引数组形成对非零元素值数组在原稀疏矩阵A(4x4)的位置的表示。

图6针对行索引数组604做进一步的说明，行索引数组604的第一个地址位置addr0存储的数值是rvalue0，本例中rvalue0的真实值是1，则表示非零元素值数组addr0地址位置的非零值所在的行列位置是第一列第二行。第二个地址位置addr1存储的数值是rvalue1，本例中rvalue1的真实值是2，则表示非零元素值数组addr1地址位置的非零值所在的行列位置是第一列第三行。本例中第一列所有非零元素值是由图5中的cvalue0，cvalue1来指示，具体来说是大于等于cvalue0，小于cvalue1所表明的范围。

图7说明了列向压缩以后的稀疏矩阵A(4x4)与稠密矩阵B(4x3)再依据外积累加方法求解矩阵内积的过程。首先依据列起点偏移量数据703与行索引数组704，计算出稀疏矩阵A(4x4)第一列的所有非零元素在非零元素值数组705中的位置，然后读取相应范围的数据与稠密矩阵B(4x3)第一行706执行外积运算。同样地，计算出稀疏矩阵A(4x4)第二列的所有非零元素在非零元素值数组705中的位置，与稠密矩阵B(4x3)第二行707再次执行外积运算。依次执行稀疏矩阵A(4x4)压缩后的第i列与密矩阵B(4x3)第i行外积运算。

图8展示了稀疏矩阵A(4x4)列向压缩的第一列808与稠密矩阵B(4x3)第一行806外积计算结果813，以及稀疏矩阵A(4x4)列向压缩的第二列809与稠密矩阵B(4x3)第二行807外积计算结果823，经过矩阵加法814，生成部分内积的结果815。外积计算结果813等价于图1中的c1(4x3)，外积计算结果823等价于图1中的c2(4x3)，因此815等价于c1(4x3)+c2(4x3)。图7中的708与709也按照上述外积过程运算，其结果再累加到815之上，即形成完整的内结果C(4x3)＝A(4x4)B(4x3)＝c1(4x3)+c2(4x3)+c3(4x3)+c4(4x3)。

稀疏矩阵A(4x4)列向压缩的第一列808与稠密矩阵B(4x3)第一行806执行外积计算后分别生成结果810、811、812。上述810-812的外积计算过程与图2中所说明的过程相同。由于稀疏矩阵A(4x4)列向压缩的第一列808之中只保留了第二行、第三行、第四行非零的元素，丢弃了第一行的零元素，故而外积计算过程中只用计算出c1(4x3)矩阵的第二行、第三行、第四行，对应生成810、811、812结果写入c1(4x3)矩阵。

通过上述的设置，对比图2之中稀疏矩阵A(4x4)第一列与稠密矩阵B(4x3)第一行，有零元素参与的外积计算，一方面列向压缩后可以只保留有意义的计算，提高了外积运算的计算效率；另一方面放弃零元素的运算，也就放弃外积计算后零值向内存写入的需求，进而提高c1(4x3)矩阵的写入效率，提升了内存写入的有效带宽。

最后，针对矩阵数量较多的稀疏矩阵A(MxK)与稠密矩阵B(KxN)可以描述成如图9所示，具体的说，首先由901模块完成稀疏矩阵A(MxK)与稠密矩阵B(KxN)原始数据的读取。902模块执行稀疏矩阵A(MxK)的列向数据压缩，具体列向数据压缩过程在图4-图6之中举例说明。下面执行外积计算，分别由903模块读取列向压缩的稀疏矩阵A(MxK)第i列以及904模块读取稠密矩阵稠密矩阵B(KxN)第i行生成ci，列与行的外积计算在图2中举例说明。再由906完成外积累加过程，累加一般是在同一块内存空间之中反复累加，将当前结果叠加到该内存空间的值之上。最后判断是否完成到稀疏矩阵A(MxK)的最后一列。如果是，则完成计算；如果否，则i自增1，回到903模块依次执行903-907模块。

需要说明的是，在本实施例中，为了方便对稀疏矩阵A(4x4)和稠密矩阵B(4x3)进行表示，对不同附图引用了不同的标号来进行表示，以方面更好的理解相对应的附图。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。