sparse matrix 的分布式存储和计算
矩阵乘法
我们先来补充一下矩阵乘法的数学知识：
- 矩阵乘法的意义: 对一个矩阵进行左乘一个矩阵的运算，相当于对该矩阵的每一列元素做线性变换；对一个矩阵进行右乘矩阵运算相当于对该矩阵的行进行线性变换。
- 矩阵乘法的计算公式： 若矩阵A为m * n的矩阵，矩阵B为n * p的矩阵，那么两矩阵相乘的公式如下：
A * B.png
用向量来解释为下图：
向量示意图.png
分布式计算
- 概念：把需要进行大量计算的工程数据分区成小块，由多台计算机分别计算，再上传运算结果后，将结果统一合并得出数据结论的科学
- 实现工具：目前最流行的分布式处理平台是Hadoop，其流程大致如下图：
Paste_Image.png
数据存储在HDFS上通过指令分发到多个Mapper（多个计算机）上进行处理，随后Mapper输出的结果（key - value）按照key进行shuffle最后传递到不同的reduce（统一合并）统计处理结果输出。
关键点：首先我们要将庞大的问题拆分成运算规则相同的子问题；随后通过设定Mapper输出的key将该合并的问题分类到同一个reduce下实现分布式计算。
稀疏矩阵的存储方式
- Coordinate(COO)
![COO存储方式.png]
如上图所示。我们对于矩阵的存储我们需要直到三个值即可：columnid, rowid, value.于是对于一个稀疏矩阵我们可以存储为三行，三个值对应一个三元组存储为一列。如右图所示。对于等于0的元素我们可以选择不存储。因为在做矩阵运算时，如上面运算公式可以看出，这些元素对结果并没有什么贡献。这样我们对于矩阵的存储可以变为O(3K)的空间复杂度。（K为非零元素的个数）,对于上面的存储方式，我们可以进一步对行索引进行压缩。就是下面所说的行压缩稀疏存储。
- Compressed Sparse Row (CSR)
![CSR存储.png]
如上图所示，这时对行压缩优化的一个存储方式。其中三元组变为：行偏移量、列索引、value。
行向量偏移量：每一行第一个非零元素在所有非零元素中的索引。也就是将所有的非零元素按矩阵的顺序排列在数组中后的索引。如左图：所有非零元素组成数组:[1,7,2,8,5,3,9,6,4],所以每行的偏移量就是该行第一个元素在该数组中的索引。(1,2,5,6) --(0,2,4,7),最后加上所有非零元素个数。
可行性分析： 我们所要做的无非是用上述三个数组的信息是否可以恢复出原有的矩阵。我们发现有几个偏移量就有几行。并且每个元素的列是记录了的，所以可以恢复出原有的矩阵。
空间复杂度: O(2K + m), m为矩阵行数
进一步优化：用两个数组来表示一个矩阵，如下面的方法
- ELLPACK (ELL)
![ELL存储.png]
该方法只是在CSR上做了进一步改进，因为我们发现只要记录了非零元素的列，我们只需直到非零元素在第几行就可以。于是我们将每一行的非零元素末尾加上 #来表示该位置为这一行的结尾。于是我们可以在上面的CSR中拿出列和值的数组，在其中插入#来区分每一行即可，上图矩阵记录如下：
column id [0,1,#,1,2,#,0,2,3,#,1,3]
value [1,7,#,2,8,#,5,3,9,#,6,4]
通过#的个数，我们可以直到该#为第几行的结尾，从而直到上一个#之前的元素是第几行元素。
空间复杂度：降到了二维的存储，并且可以格式很整齐，对后面分布式计算有很大的帮助。
矩阵乘法的分布式计算
-
可行性：我们从矩阵乘法的公式来看：
A * B.png
我们发现最终乘机得到的矩阵是一个m x p的矩阵，相当于每个元素都要进行上述公式的运算 O(n); 求出整个矩阵复杂度就为O( n m p) = n O(m p); 我们思考可以把这个问题拆分成n个子问题，每个子问题需要求出m * p 个值，而这里的每一个值都是对最终结果对应位置的贡献，也就是我们把上面加法的每一个元素放在一个map里去处理，以i#j作为key， 以![ ]为value输出作为索引为(i, j)的元素贡献值。随后在reduce中进对相同key的value加和作为最终结果中（i, j）的结果。这样我们就实现了矩阵在k个机器上的分布式计算。这也是最小粒度相乘的一个思路，将矩阵的乘积转化为数的乘积，随后叠加。如下图：
01.png
-
实现详解：
(1) 数据切分：
首先在hdfs上的文件存储我们采用ELL的存储方式，并将这个二维数组按照#切分成不同的split-part，例如上面的例子我们切分成了：
[0,1] [1,2] [0,2,3] [1,3]
[1,7] [2,8] [5,3,9] [6,4]
上下组成columnid 与 value的对应关系,一共4组，我们对于A矩阵按照列压缩存储(存行号)如下，B矩阵按照行压缩存储(存列号)也就是上面所示。
[0,2] [0,1,3] [1,2] [2,3]
[1,5] [7,2,6] [8,3] [9,4]
随后map中的工作就是将A与B中对应位置的分片值逐一相乘以(i,j)作为key输出；例如A的第一个分片与B的第一个分片相乘放在第一个mapper里去做结果如下：
(0,0) 1 ; (0, 2) 5; (1,0) 7; (1,2) 35这些都是对最终结果(i, j)位置上的值的一个分量。对于A的第二个分片和B的第二个分片我们继续做上述工作得到：(1,0) 14; (1,1) 4;(1,3) 12; (2, 0) 56; (2, 1) 16; (2, 3) 48。A的第三个分片与B的第三个分片：(0,1) 40; (0,2) 15;(2,1) 24; (2, 2) 9; (3，1) 72; (3, 2) 27。A的第四个分片与B的第四个分片：(1, 2) 54; (1, 3) 24; (3, 2) 36; (3, 3) 16;
随后每个mapper将key 设定为 (i # j); value 为计算所得到的数值；将所有的mapper值shuffle后进入reducer。这样，相同的i, j会在一起。相加后得到(1,j) 元素的值。完成计算。
转载请著名出处：http://www.jianshu.com/p/b335ad456990
[ ]: http://latex.codecogs.com/png.latex?a_i_k*b_k_j
[ELL存储.png]:https://img.haomeiwen.com/i5666544/f794e5f9814a9612.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240
[CSR存储.png]:https://img.haomeiwen.com/i5666544/414e9610a52e3a5e.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240
[COO存储方式.png]:https://img.haomeiwen.com/i5666544/1ccdc71896ffb315.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240
