CN111767941A - 一种基于对称非负矩阵分解的改进谱聚类及并行化方法 - Google Patents

Abstract

本发明公开了一种基于对称非负矩阵分解的改进谱聚类及并行化方法。通过单向循环多轮迭代的方法进行不同分区中样本点间相似度的并行计算，并采用t近邻的方法稀疏化相似度矩阵。通过交替方向乘子法来迭代求出与正交约束松弛的谱聚类目标函数等价的对称非负矩阵分解的最优解，避免使用拉普拉斯矩阵进行SVD特征分解。利用改进的K‑means算法对通过对称非负矩阵分解得到的特征向量子集进行聚类。本发明对改进谱聚类算法的计算步骤基于Spark大数据计算框架进行了并行设计及实现，使得算法不仅在聚类效果上优于传统算法，进一步解决了传统谱聚类算法在大规模数据集中计算耗时过长甚至无法完成计算的问题。

Description

一种基于对称非负矩阵分解的改进谱聚类及并行化方法

技术领域

本发明提出了一种改进谱聚类及并行化方法，涉及机器学习、大数据聚类和并行计算领域。

背景技术

谱聚类算法基于谱图理论，将聚类问题转化为图的最优划分问题，由于其能够实现对非凸等任意形状的样本空间的聚类、可有效避免局部最优解、并且可以应用于高维数据的聚类等优点，成为机器学习领域聚类算法中的研究热点。然而，随着大规模数据的普及，传统谱聚类算法由于计算相似度矩阵时空间存储代价大、特征分解的时间复杂度高，存在计算耗时过长甚至无法完成计算的问题，这在实际的大数据应用中是致命的，限制了其在很多领域上的应用。

近年来对于谱聚类算法的研究主要集中在两方面。一方面，通过在算法的实现层面上进行优化，提升算法的执行效率。针对存储相似度矩阵空间复杂度过高的问题，一类解决方法是使相似度矩阵中一些无关紧要的元素归零，将矩阵稀疏化。另一类方法是按照行或列对相似性矩阵进行采样，得到原始矩阵的低秩近似。

针对拉普拉斯矩阵特征分解时间复杂度高的问题，有研究者提出基于Spark平台采用Lanczos分解方法将拉普拉斯矩阵分解成实对称对角矩阵，再进行QR分解来提高算法的运行效率。对于最后的聚类方式，针对普遍使用的K-means聚类方法由于初始聚类中心的随机选取，存在受离群点的影响大并且聚类结果差异大的问题，有研究者使用K-means++预采样过程来确定初始聚类中心，但由于其聚类中心点选择过程中的内在有序性，在扩展方面存在着性能方面的问题。

另一方面，随着MPI、MapReduce并行计算模型的应用以及Hadoop、Spark等分布式并行框架的兴起，实现算法的并行化是提高大数据分析算法运行效率的不错的选择。Song等人利用MPI并行环境设计出并行谱聚类算法并使用大量的数据进行聚类实验，解决了传统的谱聚类算法中存在的计算性能瓶颈的问题，但存在通信开销较大、对控制的要求比较高的缺点。Fei Gao等人提出了一种分布式近似谱聚类算法。这种算法同样也是基于MapReduce编程模型进行设计，并在Hadoop平台之上实现。与MPI和Hadoop系统相比，当前应用广泛的Spark并行计算框架具有良好的优越性。MPI编程模型比较低层次，需要用户理解处理数据流机制和底层架构。Spark提供抽象化编程模型，用户只需要专注于算法的逻辑实现，并不用关心节点之间的通信、失效和恢复等问题。Hadoop只基于map和reduce这两种抽象实现为用户提供高层次的MapReduce编程模型，而map和reduce操作会产生很多中间数据，频繁的磁盘I/O读写限制了任务的高效运行。Spark编程模型将所有数据都抽象成具有丰富的并行操作算子的RDD，基于内存进行RDD的迭代计算，减少了中间结果在磁盘上的读写操作，优化了迭代算法的工作负载。

发明内容

为解决传统谱聚类算法由于相似度矩阵计算和拉普拉斯矩阵特征分解空间复杂度、时间复杂度高而无法应用于大规模数据集的问题，本发明提出了一种基于对称非负矩阵分解的改进谱聚类及并行化方法。对改进谱聚类算法的计算步骤基于Spark大数据计算框架进行了并行设计及实现，使得算法不仅在聚类效果上优于传统算法，进一步解决了传统谱聚类算法在大规模数据集中计算耗时过长甚至无法完成计算的问题。

本发明采用如下技术方案：

基于对称非负矩阵分解的改进谱聚类及并行化方法，包括以下步骤：

1)将数据集上传到HDFS文件系统中，创建Spark上下文环境，从分布式存储系统HDFS上获取数据，将其转化为弹性分布式数据集RDD；根据用户程序中设定的RDD分区数目来启动一定数目的并发线程读取数据；

2)并行化构造相似度矩阵W，通过高斯核函数计算待处理数据集样本之间的相似度值，并采用稀疏向量集方式存储；

3)采用t近邻稀疏化相似度矩阵的方法，将所述相似度矩阵W变换为稀疏化矩阵W₁，所述t为近邻参数；并将上述稀疏化相似度矩阵W₁对称化为矩阵W′；

4)并行化计算度矩阵D、邻接矩阵

5)基于ADMM优化算法迭代求解

得到对称非负矩阵分解降维后的基矩阵M；

6)行标准化矩阵M，每一行表示一个样本，对该n个样本使用基于最大最小距离法和快速距离计算的改进并行K-means算法进行聚类，得到聚类结果A₁，A₂，A₃，…，A_k。

进一步的，所述步骤2的具体步骤如下：

步骤2.1)在每个分区内，每个样本只和本分区中序号比自己大的样本使用map算子进行相似度计算，并行地计算每个分区内样本点间的距离；

步骤2.2)先使用mapPartitionsWithlndex算子对每个分区进行标记，再使用单向循环迭代法通过map算子并行地进行不同分区中样本间的相似度计算，总的迭代次数为m＝(n-1)/2，在每轮迭代中，分别计算分区标号为i和j的分区中的样本相似度，其中j＝(i+λ)％n，n是分区个数，λ为当前迭代次数，当λ＞m，终止迭代。

进一步的，所述步骤3的具体步骤如下：

步骤3.1)对每个数据样本点x_i，维护一个大小为t的最大堆H_i，每对样本点计算出的距离将被两次插入到不同的最大堆中。最终，使用filter算子过滤出每个样本点对应的堆中保留的与该样本距离最近的t个其他样本，即t个最近邻居并记为NearestData；

步骤3.2)为了之后计算的相似度矩阵是对称的，根据NearestData中的结果信息，计算出DistData。即通过flatmap算子对NearestData中每一条值进行细分，将原本行号、列号、距离值的元组一组信息，多加一条列号、行号、距离值的元组信息，变成两组信息。再通过groupBykey算子，以行号进行shuffle。shuffle后重新进行分区，保证同一个key中的数据在一个分区当中。这里得到的结果中，key值代表样本向量行的索引，key中所对应的value则是距离key所对应的样本向量最近的t个样本向量的索引和它们之间距离组成的元组。通过上述处理，即可得到对称的稀疏相似度矩阵。

步骤3.3)由于相似度的计算消耗大量资源，故采用persist算子将相似度矩阵RDD缓存起来，避免后续计算中因数据丢失而产生的重复计算。

进一步的，所述步骤4的具体步骤如下：

步骤4.1)读入稀疏对称矩阵W′，将矩阵W′以行向量的稀疏形式存储在各个节点上，使用reduceByKey算子将矩阵每行元素相加，每行元素的和即为度矩阵D对应行的对角元素，由于对角矩阵属于稀疏矩阵，为了方便存储，将其所有对角元素汇总到Driver节点构成一个数组变量，使用map算子将其取倒数并开平方，即可得到对角矩阵D^-1/2；

步骤4.2)将以数组形式存放的对角矩阵D^-1/2的对角线元素分发到各个计算节点，使用map算子从数组中取得行号索引的对角线上元素的值，与存储在本地的矩阵W′行向量中的每个元素相乘得到中间结果记为W₂。从存储对角矩阵D^-1/2的数组中获取与W₂中每行对应位置上的元素值进行相乘，得到对称非负邻接矩阵

进一步的，所述步骤5的具体步骤如下：

步骤5.1)根据目标函数

构造其增广拉格朗日函数的等价形式

其中，X和Y为原始变量，Z为辅助变量，Γ和Λ为拉格朗日乘子，ρ为惩罚项参数。读入对称非负矩阵

惩罚项参数ρ＞0，停止阈值∈，初始化Y⁰，Z⁰，Γ⁰，Λ⁰。使用ADMM优化算法进行对称非负矩阵分解，在每轮的迭代求解过程中，在集群的每个子节点上，每个分区会并行的执行原始变量、辅助变量、对偶变量的更新，Z更新后要判断是否满足停止准则，若不满足，则进入下一轮迭代。此处原始变量X、Y的更新使用大规模矩阵运算库Marlin中实现的Cholesky分解API完成，该接口使用了Spark中的分布式块矩阵(BlockMatrix)数据类型进行设计并实现了分布式的Cholesky分解操作，而对称非负邻接矩阵的计算过程使用的数据类型是Spark中的带有行索引的行矩阵(IndexedRowMatrix)的分布式矩阵类型，因此，首先需要调用IndexedRowMatrix类中的toBlockMatrix()方法进行分布式矩阵数据类型转换。

步骤5.2)按照公式

更新原始变量X，由于在高维数据下，矩阵求逆的时间过长，根据

为对称正定矩阵的性质，因此此处采用Cholesky分解以及顺代、回代求解线性方程组的方式来避免高维矩阵的求逆，具体求解步骤及公式如下；

步骤5.2.1)使用公式(Y^k-1)^TY^k-1+ρI→CC^T进行Cholesky分解得到C、C^T；

步骤5.2.2)使用公式

按行依次顺代求解j₁、j₂、…、j_n获得J；

步骤5.2.3)使用公式C^T(X^k)^T＝J通过按行依次回代的方式求解x_n、x_n-1、…、x₁来更新X^k；

步骤5.3)按照公式

更新原始变量Y，具体的求解步骤及公式如下；

步骤5.3.1)使用公式(X^k)^TX^k+ρI→DD^T进行Cholesky分解得到D、D^T；

步骤5.3.2)使用公式

顺代求解H；

步骤5.3.3)使用公式D^T(Y^k)^T＝H回代求解Y^k；

步骤5.4)在完成一次原始变量的更新后使用公式

更新辅助变量Z；

步骤5.5)在更新完一次原始变量及辅助变量之后，按照公式Λ^k＝Λ^k-1+ρ(Z^k-X^k)、Γ^k＝Γ^k-1+ρ(Z^k-Y^k)更新对偶变量Λ、Γ；

步骤5.6)重复步骤5.2)至5.5)，直到满足

时，终止迭代；

进一步的，所述步骤6的具体步骤如下：

步骤6.1)行标准化由对称非负矩阵分解降维得到的基矩阵M并初始化为K-means的聚类对象，并初始化聚类簇数k与用于终止迭代的距离阈值T；

步骤6.2)从基矩阵M中使用sample算子随机抽取一个较小的数据集S，数据抽取比例需要根据原始数据的数据量来确定，通常为1％至20％之间；

步骤6.3)使用最大最小距离法计算出抽样数据集的聚类中心作为原始聚类对象的初始聚类中心；

进一步的，所述步骤6.3)的具体步骤如下：

步骤6.3.1)从抽样得到的数据集5中使用takeSample算子随机抽取一个数据点s₁，将该点作为抽样数据集聚类中心集合C中的第一个聚类中心点c₁；

步骤6.3.2)计算抽样数据集S中其余数据点到s₁的欧式距离，将距离值最大的数据点作为第二个聚类中心点c₂；

步骤6.3.3)计算数据集S中的数据点s_i分别到聚类中心集合C中每个点c₁，c₂的欧式距离并记为d_i1，d_i2，计算max(min(d_i1，d_i2))对应的数据点s_i作为C中的第三个聚类中心点c₃；

步骤6.3.4)若集合C中有k′个聚类中心点，k′的值小于k，即可通过max(min(d_i1，d_i2，…d_ik′))计算出第k′+1个聚类中心c_k′+1；

步骤6.3.5)重复步骤6.3.4)，直到获得k个聚类中心点，即C＝{c₁，c₂，c₃，…c_k}；

步骤6.4)将k个聚类中心点通过广播的形式分发到各个节点上；在每个节点中采用一种快速计算距离的方法计算本地的样本数据到中心点的距离，使用mapToPair方法将样本点划分到距离最近的中心点所属的类中；

进一步的，所述步骤6.4)的中的快速计算距离具体步骤如下：

步骤6.4.1)矩阵M中的每一行代表一个样本点，计算每个样本点到每个聚类中心的距离；并获得每个点到最近的聚类中心的距离并记为closestCenterDist，

步骤6.4.2)计算样本中任意两点间的二范数之差的平方并记为boundNormDist；将boundNormDist和closestCenterDist的平方值进行比较，当boundNormDist的值大于closestCenterDist的平方时，无需计算欧式距离；

步骤6.4.3)当boundNormDist的值小于closestCenterDist的平方时，则开始计算两个样本点之间的欧式距离。对欧式距离进行解开计算，需要先考虑精度。按照公式precisionBound＝2.0*EPSILON*sumSquaredNorm/(normDiff*normDiff+EPSILON)求解精度。其中，sumSquaredNorm＝norm1*norm1+norm2*norm2，normDiff＝norm1-norm2，EPSILON代表最小中心距离的阈值。

步骤6.4.4)如果精度达到要求，即precisionBound＜precision，其中，precision代表Double类型的精度阈值，值一般设定为1e-6，则按照公式realDist＝sumSquaredNorm-2.0*dot(v1，v2)求解欧式距离。如果精度不满足要求，则需要进行原始的未分解的欧式距离计算。

步骤6.5)根据所有样本新的聚类情况，使用reduceBykey算子重新计算每个类簇的聚类中心；

步骤6.6)重复步骤6.4)至6.5)，直到每个样本点到其所属的中心的距离平方和最小或者当所有的样本点所属类不再发生变化的时候，聚类操作结束。

本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点：

1.本发明是一种基于对称非负矩阵分解的改进谱聚类及并行化方法。首先，通过单向循环多轮迭代的方法进行不同分区中样本点间相似度的并行计算来避免不同分区中样本间相似度的重复计算，并采用t近邻的方法将稠密的相似度矩阵稀疏化来最大化地降低相似度矩阵计算的存储开销和计算开销。

2.针对拉普拉斯矩阵特征分解时间复杂度过高的问题，本发明基于对称非负矩阵分解等价于正交约束松弛的谱聚类的理论，通过交替方向乘子法(ADMM)来迭代求出对称非负矩阵分解的目标函数的最优解，不仅降低了传统谱聚类算法中拉普拉斯矩阵使用SVD特征分解时较高的时间复杂度，同时由于对称非负矩阵分解保留了近似正交性，并且引入了非负约束，相较于传统谱聚类算法，提高了聚类的精度，增强了矩阵分解结果用于聚类的可解释性。

3.本发明通过对对称非负邻接矩阵分解降维得到的特征向量子集中的样本数据进行抽样，采用最大最小距离法确定K-means算法初始聚类中心，避免K-means聚类由于初始聚类中心点选择不恰当导致迭代次数增加，并将一种快速计算样本间距离的方法应用于改进的并行K-means算法中，进一步提升谱聚类算法的速率和精度。同时基于RDD数据的分布式存储以及丰富的RDD操作算子实现改进谱聚类算法SNMFSC的并行化，提升了算法的运行速率及处理大规模数据的能力。

附图说明

图1是基于Spark的改进谱聚类SNMFSC并行算法的流程图。

图2是改进K-means算法的流程图。

图3是改进K-means算法中样本间距离计算优化的流程图。

图4是基于Spark的改进谱聚类SNMFSC并行算法运行效率测试实验结果图。

具体实施方式

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

以下结合具体的实施例子对上述方案做进一步说明，本发明的优选实施例详述如下：

实施例一

在本实施例中，参见图1-图4，一种基于对称非负矩阵分解的改进谱聚类及并行化方法，其步骤如下：

1)将数据集上传到HDFS文件系统中，创建Spark上下文环境，从分布式存储HDFS上获取数据，将其转化为弹性分布式数据集RDD；根据用户程序中设定的RDD分区数目来启动一定数目的并发线程读取数据；

2)并行化构造相似度矩阵W，通过高斯核函数计算待处理数据集样本之间的相似度值，并采用稀疏向量集方式存储；

3)采用t近邻稀疏化相似度矩阵的方法，将所述相似度矩阵W变换为稀疏矩阵W₁，所述t为近邻参数；并将上述稀疏相似度矩阵W₁对称化为矩阵W′；

4)并行地计算度矩阵D、对称非负邻接矩阵

5)基于ADMM优化算法迭代求解

得到对称非负矩阵分解降维后的基矩阵M；

6)行标准化矩阵M，每一行表示一个样本，对该n个样本使用基于最大最小距离法和快速距离计算的改进并行K-means算法进行聚类，得到聚类结果A₁，A₂，A₃，…，A_k。

本发明对改进谱聚类算法的计算步骤基于Spark大数据计算框架进行了并行设计及实现，使得算法不仅在聚类效果上优于传统算法，进一步解决了传统谱聚类算法在大规模数据集中计算耗时过长甚至无法完成计算的问题。

实施例二

本实施例与实施例一基本相同，特别之处在于：

在本实施例中，一种基于对称非负矩阵分解的改进谱聚类及并行化方法，其步骤如下：

步骤1)将数据集上传到HDFS文件系统中，先进行SparkContext初始化，使用SparkContext中的textFile函数从HDFS读取待处理的数据集，通过textFile函数的第二个参数定义分区个数并将其转换为RDD存储，由于下文中使用单向循环迭代法求解相似度矩阵来避免分区间样本相似度的重复计算，分区个数设定为奇数；

步骤2)并行化构造相似度矩阵W，通过高斯核函数计算待处理数据集样本之间的相似度值，并采用稀疏向量集方式存储；

步骤2.1)在每个分区内，每个样本只和本分区中序号比自己大的样本使用map算子进行相似度计算，并行地计算每个分区内样本点间的距离；

步骤2.2)先使用mapPartitionsWithlndex算子对每个分区进行标记，再使用单向循环迭代法通过map算子并行地进行不同分区中样本间相似度计算。总的迭代次数为m＝(n-1)/2，在每轮迭代中，分别计算分区标号为i和j的分区中的样本相似度，其中j＝(i+λ)％n，n是分区个数，λ为当前迭代次数，当λ＞m，终止迭代。

步骤3)采用t近邻稀疏化相似度矩阵的方法，将所述相似度矩阵W变换为稀疏矩阵W₁，所述t为近邻参数；并将上述稀疏相似度矩阵W₁对称化为矩阵W′：

步骤3.1)对每个数据样本点x_i，维护一个大小为t的最大堆H_i，每对样本点计算出的距离将被两次插入到不同的最大堆中。最终，使用filter算子过滤出每个样本点对应的堆中保留的与该样本距离最近的t个其他样本，即t个最近邻居并记为NearestData；

步骤3.2)为了之后计算的相似度矩阵是对称的，根据NearestData中的结果信息，计算出DistData。即通过flatmap算子对NearestData中每一条值进行细分，将原本行号、列号、距离值的元组一组信息，多加一条列号、行号、距离值的元组信息，变成两组信息。再通过groupBykey算子，以行号进行shuffle。shuffle后重新进行分区，保证同一个key的数据在一个分区当中。这里得到的结果中，key值代表样本向量行的索引，key中所对应的value则是距离key所对应的样本向量最近的t个样本向量的索引和它们之间的距离组成的元组。通过上述处理，即可得到对称的稀疏相似度矩阵。

步骤3.3)由于相似度的计算消耗大量资源，故采用persist算子将相似度矩阵RDD缓存起来，避免后续计算中因数据丢失而产生的重复计算。

步骤4)并行地计算度矩阵D以及对称非负邻接矩阵

步骤4.1)读入稀疏对称矩阵W′，将矩阵W′以行向量的稀疏形式存储在各个节点上，使用reduceByKey算子将矩阵每行元素相加，每行元素的和即为度矩阵D对应行的对角元素，由于对角矩阵属于稀疏矩阵，为了方便存储，将其所有对角元素汇总到Driver节点构成一个数组变量，使用map算子将其取倒数并开平方，即可得到对角矩阵D^-1/2；

步骤4.2)将以数组形式存放的对角矩阵D^-1/2的对角线元素分发到各个计算节点，使用map算子从数组中取得行号索引的对角线上元素的值，与存储在本地的矩阵W′行向量中的每个元素相乘得到中间结果记为W₂。从存储对角矩阵D^-1/2的数组中获取与W₂中每行对应位置上的元素值进行相乘，得到对称非负邻接矩阵

步骤5)根据目标函数

构造其增广拉格朗日函数的等价形式

其中，X和Y为原始变量，Z为辅助变量，Γ和Λ为拉格朗日乘子，ρ为惩罚项参数。基于ADMM优化算法迭代求解

得到对称非负矩阵分解降维后的基矩阵M；

步骤5.1)输入对称非负矩阵

惩罚项参数ρ＞0，停止阈值∈，初始化Y⁰，Z⁰，Γ⁰，Λ⁰；

步骤5.2)按照公式

更新原始变量X，此处采用Cholesky分解以及顺代、回代求解线性方程组的方式来避免高维矩阵的求逆；

步骤5.2.1)使用公式(Y^k-1)^TY^k-1+ρI→CC^T进行Cholesky分解得到C、C^T；

步骤5.2.2)使用公式

按行依次顺代求解j₁、j₂、…、j_n获得J；

步骤5.2.3)使用公式C^T(X^k)^T＝J通过按行依次回代的方式求解x_n、x_n-1、…、x₁来更新X^k；

步骤5.3)按照公式

更新原始变量Y，具体的求解步骤及公式如下；

步骤5.3.1)使用公式(X^k)^TX^k+ρI→DD^T进行Cholesky分解得到D、D^T；

步骤5.3.2)使用公式

顺代求解H；

步骤5.3.3)使用公式D^T(Y^k)^T＝H回代求解Y^k；

步骤5.4)在完成一次原始变量的更新后使用公式

更新辅助变量Z；

步骤5.5)在更新完一次原始变量及辅助变量之后，按照公式Λ^k＝Λ^k-1+ρ(Z^k-X^k)、Γ^k＝Γ^k-1+ρ(Z^k-Y^k)更新对偶变量Λ、Γ；

步骤5.6)重复步骤5.2)至5.5)，直到满足

时，终止迭代；

步骤6)行标准化矩阵M，每一行表示一个样本，对该n个样本使用基于最大最小距离法和快速距离计算的改进K-means算法进行聚类，得到聚类结果A₁，A₂，A₃，…，A_k。

步骤6.1)行标准化由对称非负矩阵分解降维得到的基矩阵M并初始化为K-means的聚类对象，并初始化聚类簇数k与用于终止迭代的距离阈值T；

步骤6.2)从基矩阵M中使用sample算子随机抽取一个较小的数据集S，数据抽取比例需要根据原始数据的数据量来确定，通常为1％至20％之间；

步骤6.3)使用最大最小距离法计算出抽样数据集的聚类中心作为原始聚类对象的初始聚类中心；

步骤6.3.1)从抽样得到的数据集5中使用takeSample算子随机抽取一个数据点s₁，将该点作为抽样数据集聚类中心集合C中的第一个聚类中心点c₁；

步骤6.3.2)计算抽样数据集S中其余数据点到s₁的欧式距离，将距离值最大的数据点作为第二个聚类中心点c₂；

步骤6.3.3)计算数据集S中的数据点s_i分别到聚类中心集合C中每个点c₁，c₂的欧式距离并记为d_i1，d_i2，计算max(min(d_i1，d_i2))对应的数据点s_i作为C中的第三个聚类中心点c₃；

步骤6.3.4)若集合C中有k′个聚类中心点，k′的值小于k，即可通过max(min(d_i1，d_i2，…d_ik′))计算出第k′+1个聚类中心c_k′+1；

步骤6.3.5)重复步骤6.3.4)，直到获得k个聚类中心点，即C＝{c₁，c₂，c₃，…c_k}；

步骤6.4)将k个聚类中心点通过广播的形式分发到各个节点上；在每个节点中通过map算子采用一种快速计算距离的方法计算本地的样本数据到中心点的距离，使用mapToPair方法将样本点划分到距离最近的中心点所属的类中；

步骤6.4.1)矩阵M中的每一行代表一个样本点，计算每个样本点到每个聚类中心的距离；并获得每个点到最近的聚类中心的距离并记为closestCenterDist，步骤6.4.2)计算样本中任意两点间的二范数之差的平方并记为boundNormDist；将boundNormDist和closestCenterDist的平方值进行比较，当boundNormDist的值大于closestCenterDist的平方时，无需计算欧式距离；

步骤6.4.3)当boundNormDist的值小于closestCenterDist的平方时，则开始计算两个样本点之间的欧式距离。对欧式距离进行解开计算，需要先考虑精度。按照公式precisionBound＝2.0*EPSILON*sumSquaredNorm/(normDiff*normDiff+EPSILON)求解精度。其中，sumSquaredNorm＝norm1*norm1+norm2*norm2，normDiff＝norm1-norm2，EPSILON代表最小中心距离的阈值。

步骤6.4.4)如果精度达到要求，即precisionBound＜precision，其中，precision代表Double类型的精度阈值，值一般设定为1e-6，则按照公式realDist＝sumSquaredNorm-2.0*dot(v1，v2)求解欧式距离。如果精度不满足要求，则需要进行原始的未分解的欧式距离计算。

步骤6.5)根据所有样本新的聚类情况，使用reduceByKey算子重新计算每个类簇的聚类中心；

步骤6.6)重复步骤6.4)至6.5)，直到每个样本点到其所属的中心的距离平方和最小或者当所有的样本点所属类不再发生变化的时候，聚类操作结束；

实验方案及结果分析：

本次实验环境是在四个节点为基础的小型Spark集群上进行的，其中每个节点均配有Intel Corei7 7700 3.6GHz处理器、8GBRAM、硬盘参数为1TB/7200转、核数为四核八线程。每台节点装置Ubuntu 18.04 64位操作系统，Spark版本号为2.3.1。选取DaviesBouldin index(DBI)作为聚类精确度的评估标准。

DBI衡量聚类结果的分散程度，DBI就是计算簇内距离之和与簇外距离之比，它的计算公式为：

其中c是聚类结果簇的个数，c_i是簇i的聚类中心，σ_i是簇i中所有元素到c_i的平均值。d(c_i，c_j)，是两个簇的聚类中心的距离。当DBI值越小的时候，就代表聚类结果越准确。

使用SciKit-Leam机器学习库中的make_blobs函数来生成聚类的测试数据，该函数中可以设置需要生成的数据总数，每个数据的维度，聚类中心的数目，和每个簇之间的方差。因为测试规模的扩大，数据的维度选取为5，以降低数据之间的密集程度。聚类中心的个数选取为N/1000向下取整。具体数据样本的数据规模如表1所示，其中第一条样本数据是由球状结构组合，其他样本数据中混淆了非球状的簇结构。

基于Spark实现的改进谱聚类SNMFSC并行算法性能测试，主要测试谱聚类算法的聚类精度、运行效率及聚类数据规模。在生成的聚类样本上对基于对称非负矩阵分解的改进谱聚类并行算法(SNMFSC)与稀疏化近似谱聚类算法(SASC)、基于局部敏感哈希谱聚类算法(LSHSC)，基于Spark的K-means算法进行对比实验。其中运用到稀疏化的方法时，控制稀疏度为0.01。实验数据取的是多次实验的平均值。

表1是基于Spark的改进谱聚类SNMFSC并行算法与稀疏化近似谱聚类算法(SASC)、基于局部敏感哈希谱聚类算法(LSHSC)，基于Spark的K-means算法进行对比实验，测试聚类精度及聚类数据规模的实验结果。

表1 基于Spark的SNMFSC算法聚类精度及聚类数据规模测试结果表

由表1的实验结果可以看出，当数据样本为球状结构时，四种算法都能够精确聚类，K-means算法相对于其他算法而言能够快速完成聚类。而当遇到带有非球状数据集的数据时，K-means算法的精确度却大大降低。K-means算法对于11种数据集聚类的平均DBI为8.42，SASC对于6种数据集聚类的平均DBI为1.68，LSHSC对于8种数据集的平均DBI为1.595，SNMFSC对于11种数据集的平均DBI为1.48，所以三种谱聚类算法精确度方面较K-means提高了很多，并且在三种谱聚类算法中，SNMFSC的平均DBI最低，即具有最好的聚类效果。

从实验数据规模来看，在本文的实验环境下，SASC、LSHSC这两种谱聚类算法的数据计算规模分别为6万和8万，当数据计算规模进一步扩大时，会出现无法计算或者计算时间过长的问题，存在着性能瓶颈。而基于Spark的SNMFSC算法的数据计算规模可以达到57万，表明改进的并行谱聚类算法具有更好的数据扩展性。

根据图4基于Spark的改进谱聚类SNMFSC并行算法运行效率测试实验结果图，LSHSC相对于SASC执行时间更短，处理的数据量级更大一些。同时，也可以看出相较于另外两种谱聚类算法，SNMFSC算法在图中所有规模的数据集上均使用最短的运行时间，表明改进的算法具有更高的运行速率。

综上所述，本发明是一种基于对称非负矩阵分解的改进谱聚类及并行化方法。通过单向循环多轮迭代的方法进行不同分区中样本点间相似度的并行计算，并采用t近邻的方法稀疏化相似度矩阵。通过交替方向乘子法(ADMM)来迭代求出与正交约束松弛的谱聚类目标函数等价的对称非负矩阵分解的最优解，避免使用拉普拉斯矩阵进行SVD特征分解。引入了最大最小距离法确定初始聚类中心以及一种快速计算样本间距离的方法改进传统K-means算法，使用于对对称非负矩阵分解降维得到的特征向量子集进行聚类的步骤中。本发明对改进谱聚类算法的计算步骤基于Spark大数据计算框架进行了并行设计及实现，使得算法不仅在聚类效果上优于传统算法，进一步解决了传统谱聚类算法在大规模数据集中计算耗时过长甚至无法完成计算的问题。

上面对本发明实施例结合附图进行了说明，但本发明不限于上述实施例，还可根据本发明的发明创造的目的做出多种变化，凡依据本发明技术方案的精神实质和原理下做的改变、修饰、替代、组合或简化，均应为等效的置换方式，只要符合本发明的发明目的，只要不背离本发明基于对称非负矩阵分解的改进谱聚类及并行化方法的技术原理和发明构思，都属于本发明的保护范围。

Claims (8)

1.一种基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，其步骤如下：

1)将数据集上传到HDFS文件系统中，创建Spark上下文环境，从分布式存储HDFS上获取数据，将其转化为弹性分布式数据集RDD；根据用户程序中设定的RDD分区数目来启动一定数目的并发线程读取数据；

2)并行化构造相似度矩阵W，通过高斯核函数计算待处理数据集样本之间的相似度值，并采用稀疏向量集方式存储；

3)采用t近邻稀疏化相似度矩阵的方法，将所述相似度矩阵W变换为稀疏矩阵W₁，所述t为近邻参数；并将上述稀疏相似度矩阵W₁对称化为矩阵W′；

4)并行地计算度矩阵D、对称非负邻接矩阵

5)基于ADMM优化算法迭代求解

得到对称非负矩阵分解降维后的基矩阵M；

6)行标准化矩阵M，每一行表示一个样本，对该n个样本使用基于最大最小距离法和快速距离计算的改进并行K-means算法进行聚类，得到聚类结果A₁，A₂，A₃，…，A_k。

2.根据权利要求1所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤2的具体步骤如下：

步骤2.1)在每个分区内，每个样本只和本分区中序号比自己大的样本使用map算子进行相似度计算，并行地计算每个分区内样本点间的距离；

步骤2.2)先使用mapPartitionsWithlndex算子对每个分区进行标记，再使用单向循环迭代法通过map算子并行地进行不同分区中样本间的相似度计算，总的迭代次数为m＝(n-1)/2，在每轮迭代中，分别计算分区标号为i和j的分区中的样本相似度，其中j＝(i+λ)％n，n是分区个数，λ为当前迭代次数，当λ＞m，终止迭代。

3.根据权利要求1所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤3的具体步骤如下：

步骤3.1)对每个数据样本点x_i，维护一个大小为t的最大堆H_i，每对样本点计算出的距离将被两次插入到不同的最大堆中；最终，使用filter算子过滤出每个样本点对应的堆中保留的与该样本距离最近的t个其他样本，即t个最近邻居并记为NearestData；

步骤3.2)为了之后计算的相似度矩阵是对称的，根据NearestData中的结果信息，计算出DistData；即通过flatmap算子对NearestData中每一条值进行细分，将原本行号、列号、距离值的元组一组信息，多加一条列号、行号、距离值的元组信息，变成两组信息；再通过groupBykey算子，以行号进行shuffle；shuffle后重新进行分区，保证同一个key的数据在一个分区当中；这里得到的结果中，key值代表样本向量行的索引，key中所对应的value则是距离key所对应的样本向量最近的t个样本向量的索引和它们之间的距离组成的元组；通过上述处理，即可得到对称的稀疏相似度矩阵；

步骤3.3)由于相似度的计算消耗大量资源，故采用persist算子将相似度矩阵RDD缓存起来，避免后续计算中因数据丢失而产生的重复计算。

4.根据权利要求1所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤4的具体步骤如下：

步骤4.1)读入稀疏对称矩阵W′，将矩阵W′以行向量的稀疏形式存储在各个节点上，使用reduceByKey算子将矩阵每行元素相加，每行元素的和即为度矩阵D对应行的对角元素，由于对角矩阵属于稀疏矩阵，为了方便存储，将其所有对角元素汇总到Driver节点构成一个数组变量，使用map算子将其取倒数并开平方，即可得到对角矩阵D^-1/2；

步骤4.2)将以数组形式存放的对角矩阵D^-1/2的对角线元素分发到各个计算节点，使用map算子从数组中取得行号索引的对角线上元素的值，与存储在本地的矩阵W′行向量中的每个元素相乘得到中间结果记为W₂；从存储对角矩阵D^-1/2的数组中获取与W₂中每行对应位置上的元素值进行相乘，得到对称非负邻接矩阵

5.根据权利要求1所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤5的具体步骤如下：

步骤5.1)根据目标函数

构造其增广拉格朗日函数的等价形式

其中，X和Y为原始变量，Z为辅助变量，Γ和Λ为拉格朗日乘子，ρ为惩罚项参数；读入对称非负矩阵

惩罚项参数ρ＞0，停止阈值∈，初始化Y⁰，Z⁰，Γ⁰，Λ⁰；使用ADMM优化算法进行对称非负矩阵分解，在每轮的迭代求解过程中，在集群的每个子节点上，每个分区会并行的执行原始变量、辅助变量、对偶变量的更新，Z更新后要判断是否满足停止准则，若不满足，则进入下一轮迭代；此处原始变量X、Y的更新使用大规模矩阵运算库Marlin中实现的Cholesky分解API完成，该接口使用了Spark中的分布式块矩阵(BlockMatrix)数据类型进行设计并实现了分布式的Cholesky分解操作，而对称非负邻接矩阵的计算过程使用的数据类型是Spark中的带有行索引的行矩阵(IndexedRowMatrix)的分布式矩阵类型，因此，首先需要调用IndexedRowMatrix类中的toBIockMatrix()方法进行分布式矩阵数据类型转换；

步骤5.2)按照公式

更新原始变量X，由于在高维数据下，矩阵求逆的时间过长，根据

为对称正定矩阵的性质，因此此处采用Cholesky分解以及顺代、回代求解线性方程组的方式来避免高维矩阵的求逆，具体求解步骤及公式如下；

步骤5.2.1)使用公式(Y^k-1)^TY^k-1+ρI→CC^T进行Cholesky分解得到C、C^T；

步骤5.2.2)使用公式

按行依次顺代求解j₁、j₂、…、j_n获得J；

步骤5.2.3)使用公式C^T(X^k)^T＝J通过按行依次回代的方式求解x_n、x_n-1、…、x₁来更新X^k；

步骤5.3)按照公式

更新原始变量Y，具体的求解步骤及公式如下；

步骤5.3.1)使用公式(X^k)^TX^k+ρI→DD^T进行Cholesky分解得到D、D^T；

步骤5.3.2)使用公式

顺代求解H；

步骤5.3.3)使用公式D^T(Y^k)^T＝H回代求解Y^k；

步骤5.4)在完成一次原始变量的更新后使用公式

更新辅助变量Z；

步骤5.5)在更新完一次原始变量及辅助变量之后，按照公式Λ^k＝Λ^k-1+ρ(Z^k-X^k)、Γ^k＝Γ^k-1+ρ(Z^k-Y^k)更新对偶变量Λ、Γ；

步骤5.6)重复步骤5.2)至5.5)，直到满足

时，终止迭代。

6.根据权利要求1所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤6的具体步骤如下：

步骤6.1)行标准化由对称非负矩阵分解降维得到的基矩阵M并初始化为K-means的聚类对象，并初始化聚类簇数k与用于终止迭代的距离阈值T；

步骤6.2)从基矩阵M中使用sample算子随机抽取一个较小的数据集S，数据抽取比例需要根据原始数据的数据量来确定，通常为1％至20％之间；

步骤6.3)使用最大最小距离法计算出抽样数据集的聚类中心作为原始聚类对象的初始聚类中心；

步骤6.4)将k个聚类中心点通过广播的形式分发到各个节点上；在每个节点中采用一种快速计算距离的方法计算本地的样本数据到中心点的距离，使用mapToPair方法将样本点划分到距离最近的中心点所属的类中；

步骤6.5)根据所有样本新的聚类情况，使用reduceByKey算子重新计算每个类簇的聚类中心；

步骤6.6)重复步骤6.4)至6.5)，直到每个样本点到其所属的中心的距离平方和最小或者当所有的样本点所属类不再发生变化的时候，聚类操作结束。

7.根据权利要求6所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤6.3)的具体步骤如下：

步骤6.3.1)从抽样得到的数据集S中使用takeSample算子随机抽取一个数据点s₁，将该点作为抽样数据集聚类中心集合C中的第一个聚类中心点c₁；

步骤6.3.2)计算抽样数据集S中其余数据点到s₁的欧式距离，将距离值最大的数据点作为第二个聚类中心点c₂；

步骤6.3.3)计算数据集S中的数据点s_i分别到聚类中心集合C中每个点c₁，c₂的欧式距离并记为d_i1，d_i2，计算max(min(d_i1，d_i2))对应的数据点s_i作为C中的第三个聚类中心点c₃；

步骤6.3.4)若集合C中有k′个聚类中心点，k′的值小于k，即可通过max(min(d_i1，d_i2，…d_ik′))计算出第k′+1个聚类中心c_k′+1；

步骤6.3.5)重复步骤6.3.4)，直到获得k个聚类中心点，即C＝{c₁，c₂，c₃，…c_k}。

8.根据权利要求6所述的基于对称非负矩阵分解的改进谱聚类及并行化方法，其特征在于，所述步骤6.4)的中的快速计算距离方法具体步骤如下：

步骤6.4.1)矩阵M中的每一行代表一个样本点，计算每个样本点到每个聚类中心的距离；并获得每个点到最近的聚类中心的距离并记为closestCenterDist，

步骤6.4.2)计算样本中任意两点间的二范数之差的平方并记为boundNormDist；将boundNormDist和closestCenterDist的平方值进行比较，当boundNormDist的值大于closestCenterDist的平方时，无需计算欧式距离；

步骤6.4.3)当boundNormDist的值小于closestCenterDist的平方时，则开始计算两个样本点之间的欧式距离；对欧式距离进行解开计算，需要先考虑精度；按照公式precisionBound＝2.0*EPSlLON*sumSquaredNorm/(normDiff*normDiff+EPSlLON)求解精度；其中，sumSquaredNorm＝norm1*norm1+norm2*norm2，normDiff＝norm1-norm2，EPSlLON代表最小中心距离的阈值；

步骤6.4.4)如果精度达到要求，即precisionBound＜precision，其中，precision代表Double类型的精度阈值，值一般设定为1e-6，则按照公式reaIDist＝sumSquaredNorm-2.0*dot(v1，v2)求解欧式距离；如果精度不满足要求，则需要进行原始的未分解的欧式距离计算。

Images