Spark中的PIC分类及试验

网络或关联的现象在现代社会中非常普遍，比如电话网络，社交网络，网页的链接，文献的引用，基因分析等等。用我们能够理解的方式来思考关系，最直接也最直观的方式就是关系图：

       但是这个方式只能针对简单关系，节点多了，联系复杂了也就没法继续玩下去了。像下面这种关系网络该怎么用手工的方式画出来和分析呢？而这种复杂的关系图在电话、微博、微信中却是再常见不过的了。

       2、300年前的数学家们想了一个办法出来，他们用矩阵来表示这种关系，这样就可以通过分析矩阵来分析网络的关系了。当然缺点就是与我们普通人的距离一下就拉远了。好在随着理论和技术的发展，现在我们也可以试着来看看这些方头巴脑的矩阵到底能把关系网络搞成什么样子。

完整的示意关系表：

1,2,10.0；1,3,8.0；1,4,12.0；2,3,7.0；2,4,15.0；3,4,13.0；4,5,0.2

5,6,20.0；5,7,15.0；5,8,6.0；6,7,8.0；7,8,22.0；7,9,1.0；8,9,0.6

9,10,15.0；9,11,6.0；9,12,10.0；9,13,8.0；9,14,0.05；10,11,7.0

10,12,30.0；10,13,5.0；11,12,6.0；11,13,20.0；12,13,4.0

以往对图的划分，普遍的方式是基于拉普拉斯矩阵的使用。其过程简洁、分类准确，从60年代末到现在，对其研究与使用一直经久不衰，目前似乎仍然占据主流地位。

但是拉普拉斯矩阵分析关系图必须求解矩阵的特征值，对小矩阵还没有太多的资源要求，不过对巨型矩阵的分析就需要相当多的资源才能操作，所以在通常的硬软件环境下，难以对大规模的关系网络做学习和划分试验。

2014年3月推出的Spark-1.3.0版本，增加了PIC分类算法，大幅度的简化了对图的关系划分，这是spark社区基于2010年卡内基.梅隆大学Frank.Lin的一篇论文实现的。PIC算法的全称是:Power Iteration Clustering，字面的意思是幂迭代分类，也有资料译为快速迭代法。这个算法的核心内容相当简单，就是矩阵和向量的乘法迭代计算。完全不用去求解矩阵的特征值，而且还很方便做分布计算，因而是个非常完美的算法。

有几个步骤或关键的地方：

1、先将关系图转为矩阵表示，识别号需要转换为从0开始的索引数

2、初始向量为每行初值之和除以矩阵总和，即初始向量的构成和归一化

3、对矩阵的每行做归一化处理，使其转变为随机矩阵【或称概率转移矩阵，总之是这类矩阵的名字而已】

4、其中(2)的处理比较特别。按理可以是任意非零向量，但按(2)的方式得到的初始向量收敛到稳态的速度很快

5、迭代后得到的向量再对其做归一化处理

对上图做PIC计算后得到的分类结果和收敛过程的示意图如下：

看得出来PIC算法的计算过程还是比较简单的，结果也很可靠。

不过个人觉得算法的论文和spark的实现似乎还有可以优化的地方，特别是后者，初期版本存在bug。本着先试用再思考的学习过程，就从spark的优化说起吧。

Spark的实现利用了自身的graph包组件。但是源代码初期的版本有缺陷，直到最近1.6.1版本才去除了这个bug。以1.5.1版源码为例：

第380行：curG = GraphImpl.fromExistingRDDs(VertexRDD(v1),g.edges)

上述代码虽然能够通过编译，但是不能实现点属性值【迭代向量值】的更新，去掉斜体字部分才行：

curG = GraphImpl.fromExistingRDDs(v1,g.edges)

或者改为：

curG = GraphImpl.fromExistingRDDs(vertices = v1,edges = g.edges)

在1.6.1版中已经将上述图的生成代码全部更新为

Graph(VertexRDD(v1),g.edges)

尽管去除了bug，但是不加以优化的话，仍然存在内存耗用过快，速度变慢的情况，比如下面不到20个点的微网络。这个网络纯粹是用来试验的，在一台10G【spark，Excutor的内存】的笔记本上计算，按理是不该有任何问题的。

上图的数据关系：

55,2,15；55,4,10；2,4,40；2,9,1.0；4,5,5；4,6,20；55,6,4

6,7,10；7,8,15；55,8,3；33,55,8；33,11,1；33,12,1；33,13,1

33,10,1；33,14,1；33,15,1；33,16,10；33,77,1；77,18,15

77,19,18；77,20,13；77,21,11；77,22,22

上面这个图的矩阵，在超过200次迭代之后，速度开始变慢，继续下去，能够勉强通过迭代计算，但之后，在接下来的Kmenas算法中就会停滞不前，最后因内存不够退出计算，这个情况很难让人接受。下面有两个改进，能够很大程度上克服图算法对资源的过度消耗。

参考spark的PageRank算法，发现该算法在迭代之初，就对图做了

curG = cache()

缓冲，因此提高了内存使用效率。将该代码引入PIC算法后，的确显著加快了矩阵迭代速度，提速效果极佳。但后期，随着迭代次数的增加，仍会面临内存不足Kmeans退出的问题。

另外，迭代计算最后的图更新，实际上只是顶点属性的更新，因此可以考虑join的方法更新图

curG = curG.joinVertices(v1){(x,y,z)=>z}

这两个改动之后，PIC的图实现方法至少对微网络的试验是没有什么问题了。

Spark的PIC算法，网络上几乎没有什么资料，估计跟算法的粗糙和给的例子过于简单【即使早期存在bug的时候也能“正确区分”出来，但是略微复杂的网络划分就完全是一场稀里糊涂的噩梦】都有关系。

矩阵的幂迭代算法，听上去有点吓人，但是实际上可以理解为数组乘以数组，所以方法的实现和计算过程其实是再简单不过的了，而且scala的数组也可以通过par方法做并行计算，因此单机的多线程计算和多机的分布式计算，都是很容易实现的。原论文作者Frank.Lin也特别强调了该算法利于分布计算。期间还考虑并使用过scala的breeze数学计算包，因为包有相当多的矩阵计算函数，比如计算矩阵的特征值和特征向量。但breeze只有密集矩阵的计算，没有稀疏矩阵的计算函数，因此在生成矩阵时，受内存的限制，最多只能生成比较有限的维度，对于上百万、上亿维度的矩阵则没有办法。所以最终还是放弃了，只是在对小规模网络做对比试验时用一下，比如对PIC和拉普拉斯算法的比较。下面列出想到的一些方法的比较：

限于时间和精力，试验时采用了第二种和最后一种方法。前者采用稀疏矩阵实现PIC算法，后者主要为了使用拉普拉斯算法。两种方式都是单机操作，后面会看到，尽管不是分布计算，但在200万节点规模下SparseMatrix还是能够完成计算的。对于小规模的试验网络，由于迭代实际上是算术运算，所以速度相当快，远远超过graphx包的图算法。作为试验，这种单机算术计算还是完全可行的。

说完spark的优化，再回头看看PIC算法。目前的PIC算法似乎只对纯粹的团状图或二分图比较有效，但是对混合图则比较弱，个人的试验结论是区分不出来，需要改进。而拉普拉斯方法则不存在这个问题，能够很完美的划分。

所谓团状图，指内部关联度强，外部关联度弱或没有关联；

所谓结构图，指单点对多点或没有三角形回路关系的图，比如星状图或二分图

所谓混合图，指上述两种图的混合体

比如对前面的一个图，重新画出来：

上图目视可以分为三个子图，最左边的是团状图，中间和右边是结构图。

PIC和拉普拉斯划分结果：

可见在混合图的划分上，拉氏划分还是占相当的优势的。而PIC算法出现了跳空节点分在一组的情况，这对关系划分来说是完全不能接受的。那么PIC算法能不能做点改进达到后者的效果呢？我试了一个办法，还是相当不错的。微小网络试验时不再出现跳空的节点组，但是孤点仍然存在。这个问题可以通过降低分组数来改善。

这个方法，对200万节点的网络分析也有满意的结果。

对大规模网络的分类，也许还有两个问题需要提及。最初做分类试验的时候，曾经想一蹴而就，将200万个点一次性分割为平均10个一组的圈子，这样大约需要20万组。但是刚做1万次【32核，spark excutor内存为160G的服务器】时就失败了，后来降到2000才完成计算，再试图上升到3000个分组时，又失败了。而且结果也不是很满意，有很多组只有一个点的孤点划分。这说明两个问题，第一个就是最多不超过3000的分组距希望的分组还差很远，第二个就是感觉不应该出现一个孤点的分组。

第一个问题睡一觉后就反应过来了：如果不能一蹴而就，那么就多蹴以就。比如先做少量分组的粗粒度划分，然后对每一个子图形做再次划分，这样不断下去直到可能的最小分组。经考虑和试验，个人认为粒度由粗到细的划分过程比较合理。粗了可以继续细分，但细了却容易分错，而且不便于检查及重算。

第二个问题比较麻烦，想不出来为什么，考虑到Kmeans本身就是个局部最优的算法，也不能有太多的期望，一时也想不出什么办法。好像没几天，spark-1.6.0发布了，在看有哪些变化的介绍时，发现这个版本新增了一个Bisecting的分类算法，大约是对Kmeans算法的改进。于是又花了些时间学这个新算法。具体原理没细看，如果仅就代码的编写而言，倒是巨简单，可以直接把Kmeans改成Bisecting就行了。Bisectting这个算法比较有意思，首次分类时不会出现一个孤点的分组，这是其与Kmeans最大的不同，对关系划分也是至关重要的。

对于关系网络而言，按常理是不应该出现孤点的情况，至少都应该是两个点一组，所以这就是Bisecting算法最吸引人的地方，因此果断换下Kmeans，改用新人上场。

接下来的试验却又有点意外。Bisecting算法比Kmeans的速度要慢点，看其代码，似乎比Kmeans更为复杂，好在也没慢多少，因此也不算什么。但在后续分类时也出现了孤点。当然，最搞怪的是它的分组能力——越往后，越分不开主集。所谓主集，是指在分类结果有巨量的节点分在一组或两组中，几乎占了总结点数的90%，其余10%的节点则零星分散其它组中。对其中的一个主集做再次分类，结果仍然是90%的点聚在一起，也就是说Bisecting划分不开主集。而且此时也会出现孤立的点了。

此时换Kmeans上场操刀，则能顺利的劈开主集。

没有细看Bisecting的算法代码，猜测跟Bisecting使用缺省的最小区隔距离有关。Kmeans可以设置epsilon值，而Bisecting内部做子分类时调用的是Kmeans算法，故用的缺省值是1e-4，因此无法劈开更小的分类。不知这个猜测对不对，总之Bisecting对主集的拆分是不得力的。

好在刚才两个算法的使用方法几乎完全一样，因此在具体操作时可以使用开关方式来分别使用两个算法。比如在首次分类时用的Bisecting算法，其后的分类全部用的是Kmeans。

在200万节点的分类中，对其中一个节点的关系分类一直做到第8次，才见分晓。

对更大规模的节点分组，应该也是可以如此操作的。

 

参考资料：

1、谱聚类算法：http://www.cnblogs.com/sparkwen/p/3155850.html

2、PIC原论文：https://www./publication/221346512_Power_Iteration_Clustering?enrichId=rgreq-c230346b-458d-4453-b94f-596f6934b371&enrichSource=Y292ZXJQYWdlOzIyMTM0NjUxMjtBUzoxMDIyMDE5MjE5MDA1NDdAMTQwMTM3ODI0NTM4Nw==&el=1_x_2

3、A Tutorial on Spectral Clustering：https://www./publication/234801250_A_Tutorial_on_Spectral_Clustering