【Scala-spark.mlib】本地矩阵乘法计算效率比较（稠密稀疏哪家强？）

矩阵乘法效率比较

1. 矩阵乘法
2. 效率比较
- 2.1. DenseMatrix(50% zeros) X DenseMatrix
- 2.2. SparseMatrix X DenseMatrix
- - 2.2.1. SparseMatrix(50% zeros) X DenseMatrix
  - 2.2.2. SparseMatrix(0% zeros) X DenseMatrix
3. 底层探究
- 3.1. gemm方法传参
- 3.2. gemm方法实现

1. 矩阵乘法

在很多的科学计算中，矩阵运算是十分基础而重要的计算步骤，有时甚至因为大量迭代运算或大规模矩阵计算而带来十分巨大的资源开销和时间消耗。spark作为分布式计算的框架，在解决大规模计算时有得天独厚的优势，而作为spark中主流的科学计算库mlib，其中的矩阵运算是我们经常会使用到的。之前的博文详细介绍了本地存储的两种矩阵形式——稠密矩阵DenseMatrix和稀疏矩阵SparseMatrix，而mlib也提供了多种分布式矩阵存储的方式和计算方法，这个将在之后的博文中一一探究，本文将着重探寻本地矩阵乘法的效率问题，也同样是分布式矩阵计算的效率基础。
矩阵乘法相信大家都已经十分熟悉了，通过行与列的对应相乘累加，获得结果矩阵中的一个个元素：

而当相乘的两个矩阵中有一个换成了稀疏矩阵，则矩阵乘法所需要计算的元素数量仅为稀疏矩阵的元素数量，通过稀疏矩阵所存储的非0元素的索引，可以跳过0元素的乘法运算，减少不必要的计算开销。

所以从理论上说，利用稀疏矩阵进行乘法计算，则计算所需时间随含0元素比率增加而减少，所以本文正是要讨论在spark的mlib中，矩阵乘法运算利用哪种矩阵存储形式更快。

2. 效率比较

为了使得矩阵乘法运算的时间差异表征不同运算方式带来的时间差异，计算矩阵所需的CPU时间必须远大于其它外界软硬件因素导致的计时误差，故本测试选用100x100含有10000个元素的矩阵进行测试。再考虑到不同计算方法下，尽可能保证计算的公平性和稳定性，故计算所用矩阵的非0值采用固定方法生成，生成测试矩阵的代码如下：

val numSize: Int = 10000 //设置矩阵元素总数
val rowSize: Int = 100 //设置矩阵行数
val colSize: Int = numSize / rowSize //计算矩阵列数
val zerosRate: Int = 50 //设置B矩阵产生0元素的概率（0-100）
val array1 = new Array[Double](numSize)
val array2 = new Array[Double](numSize)
for (i <- 1 to numSize) {array1(i-1) = i //按顺序从小到大填充A矩阵if (Random.nextInt(100) < zerosRate) {array2(i-1) = 0 //概率填充0元素} else {array2(i-1) = numSize - i + 1 //按顺序从大到小填充B矩阵}
}
//创建A稠密矩阵
val matrix1 = new DenseMatrix(rowSize, colSize, array1, true)
//创建B稠密矩阵
val matrix2 = new DenseMatrix(rowSize, colSize, array2, true)
//创建B稀疏矩阵
val matrix2SP: SparseMatrix = matrix2.toSparse

由于mlib仅提供DenseMatrix和SparseMatrix对DenseMatrix的矩阵乘法方法，所以为了尽可能保证计算的公平性，所有计算将从B矩阵调用multiply方法，为了防止页帧高速缓存命中而使得多次重复计算时的性能提升，每种计算都独立运行10次。

2.1. DenseMatrix(50% zeros) X DenseMatrix

// Dense x Dense
val start1: Long = new Date().getTime
val res1: DenseMatrix = matrix2.multiply(matrix1)
val end1: Long = new Date().getTime
val lastTime1: Long =  end1 - start1
println("Dense x Dense : " + lastTime1.toString + " ms")

10次运算时间(ms)：

75	78	73	73	66	66	69	78	77	74

10次均值：72.9ms

2.2. SparseMatrix X DenseMatrix

// Sparse x Dense
val start2: Long = new Date().getTime
val res2: DenseMatrix = matrix2SP.multiply(matrix1)
val end2: Long = new Date().getTime
val lastTime2: Long = end2 - start2
println("Sparse x Dense : " + lastTime2.toString + " ms")

2.2.1. SparseMatrix(50% zeros) X DenseMatrix

先对与第一次实验相同的含有50%的0元素的稀疏矩阵进行10次测试，这10次的运算时间(ms)：

30	27	28	29	23	28	28	24	26	26

10次均值：26.9ms
可以看出，在相同比例0元素的情况下，利用SparseMatrix代替DenseMatrix进行存储和矩阵乘法运算能够节省接近2/3的时间，计算效率增加2.71倍。但我突然意识到，SparseMatrix矩阵的创建时间开销比DenseMatrix大，这样对比并不公平，于是我修改代码，将原有的DenseMatrix转换为SparseMatrix的过程时间也计算在内，代码如下：

val start2: Long = new Date().getTime
val res2: DenseMatrix = matrix2.toSparse.multiply(matrix1)
val end2: Long = new Date().getTime
val lastTime2: Long = end2 - start2
println("Sparse x Dense : " + lastTime2.toString + " ms")

得到10次的测试结果(ms)：

28	27	28	28	26	29	29	28	31	27

10次均值：28.1ms
结果变化并不大，从这样的测试可以得到两点结论：
1. 矩阵从DenseMatrix调用toSparse转化为SparseMatrix的时间开销有，但并不大。
2. 对比DenseMatrix，SparseMatrix的存储形式在矩阵运算时有很大的效率优势。

2.2.2. SparseMatrix(0% zeros) X DenseMatrix

这时我突然想到，如果不含0的矩阵转化为稀疏矩阵的存储形式进行存储计算，理论上计算效率和速度应当与稠密矩阵的计算相当，由此我改变原代码中的含0率，令矩阵变为稠密矩阵，而采用稀疏矩阵形式进行存储，计算过程仍然采用上一节改进的DenseMatrix先转换SparseMatrix再进行矩阵乘法的流程，进行了10次独立运行测试，得到结果(ms)：

33	33	33	38	33	34	31	39	34	33

10次均值：34.1ms
可以看出，计算时间确实增加了，增幅达到21.3%，但仍然远远不及DenseMatrix的72.9ms计算时间，究竟是什么原因导致了这样的测试结果，我们需要深入mlib的源码来探究。

3. 底层探究

仔细查看class DenseMatrix和class SparseMatrix的源码会发现，两个类中都并没有multiply方法，那么multiply一定在它们共同的父类trait Matrix中被实现了，经过查找，果真在这个“接口”中发现了multiply的实现方法：

def multiply(y: DenseMatrix): DenseMatrix = {val C: DenseMatrix = DenseMatrix.zeros(numRows, y.numCols)BLAS.gemm(1.0, this, y, 0.0, C)C
}

这样也就解释了为什么DenseMatrix和SparseMatrix都只能与DenseMatrix进行乘积，因为它们继承了共同的一个multiply方法，那么这个方法又是怎么实现稠密矩阵和稀疏矩阵的矩阵乘法操作的呢，就需要进入mlib的矩阵向量计算模块中这个神秘的私有静态类BLAS，虽然它是私有类，但它却实现了许多有用的功能，今天我们先探究multiply方法所调用到的这个gemm方法。

3.1. gemm方法传参

先来看看gemm的参数列表：

再对应方法体上的官方注释，不难理解这些参数的含义：

而需要实现矩阵的A*B，则只需要使alpha = 1，beta = 0，C为任意矩阵即可，这时我们再返回主调函数中，看看这些参数是怎么设置的：

可以看出，这些参数正是按照所需要的数值设置，而C矩阵被定义为了一个与输入矩阵等大小的0矩阵，这也是执行效率最高的一种方式。知道了gemm方法是如何调用的，那么它内部如何实现不同的DenseMatrix和SparseMatrix的矩阵计算的呢？

3.2. gemm方法实现

继续看gemm的方法体：

这里对输入的矩阵A，也就是主调的对象类型进行了判断，再次分支进入gemm的两个重写方法中，我们先进入DenseMatrix方法看看gemm如何实现稠密矩阵的矩阵相乘的：

gemm针对DenseMatrix的实现十分简单，在对各个输入矩阵的矩阵主序进行判断和标识后，就输入了nativeBLAS.dgemm方法中，而这个方法是com.github.fommil.netlib向量计算包中的方法，针对顺序存储的元素数组按行列进行拆分相乘并累加。这样的调用也就解释了为什么DenseMatrix的multiply方法运行会如此之慢，因为其跨包调用了计算方法，而这样比较简单的计算方法我相信如果在gemm方法中编码实现将在效率上有不小的进步。
接下来，我们进入gemm针对SparseMatrix的实现方法，看看效率如此之高的稀疏矩阵乘积计算是如何实现的：
（gemm对SparseMatrix的乘积运算是硬编码在方法体内的，但针对输入的矩阵主序不同有许多计算分支，主要差异在行列选取的不同上，这里选取第一个分支进行分析）

A（主调矩阵——SparseMatrx）为行主序，B（DenseMatrix）为列主序进行分析）

忽略C矩阵的运算操作，从代码逻辑中可以看出，通过AcolPtrs(rowCounterForA)和AcolPtrs(rowCounterForA+1)的计算可以得到A稀疏矩阵每行的坐标始末，再通过Bstart + ArowIndices(i)得到B矩阵每个列元素的所在位置，通过循环累加得到一个元素的值，从而不断得到整个结果矩阵，这样的计算只针对稀疏矩阵存储中的非0元素，从而使得计算时间大幅减少，不会对0元素进行无意义的计算，使得整体矩阵乘积的计算速度十分迅速。
所以综上所述，即便是稠密矩阵的计算，都采用利用稀疏矩阵SparseMatrix调用multiply方法得到矩阵乘积能得到最大的计算效率。当然通过对gemm的修改，直接将稠密矩阵乘积的方法写入也可以提高DenseMatrix直接调用multiply方法的效率，但是作为private修饰私有静态类中的私有方法，无法利用继承进行重写，只能对框架源码进行修改，这会带来很大的风险，所以不建议这么做。