【技术分享】快速迭代聚类

本文原作者:尹迪,经授权发布。

1 谱聚类算法的原理

  在分析快速迭代聚类之前,我们先来了解一下谱聚类算法。谱聚类算法是建立在谱图理论的基础上的算法,与传统的聚类算法相比,它能在任意形状的样本空间上聚类且能够收敛到全局最优解。 谱聚类算法的主要思想是将聚类问题转换为无向图的划分问题。

  • 首先,数据点被看做一个图的顶点v,两数据的相似度看做图的边,边的集合由E=AijE=Aij表示,由此构造样本数据集的相似度矩阵A,并求出拉普拉斯矩阵L
  • 其次,根据划分准则使子图内部相似度尽量大,子图之间的相似度尽量小,计算出L的特征值和特征向量
  • 最后,选择k个不同的特征向量对数据点聚类

  那么如何求拉普拉斯矩阵呢?

  将相似度矩阵A的每行元素相加就可以得到该顶点的度,我们定义以度为对角元素的对角矩阵称为度矩阵D。可以通过AD来确定拉普拉斯矩阵。拉普拉斯矩阵分为规范和非规范两种,规范的拉普拉斯矩阵 表示为L=D-A,非规范的拉普拉斯矩阵表示为L=I−D−1AL=I−D−1A 。

  谱聚类算法的一般过程如下:

  • (1)输入待聚类的数据点集以及聚类数k
  • (2)根据相似性度量构造数据点集的拉普拉斯矩阵L
  • (3)选取L的前k个(默认从小到大,这里的k和聚类数可以不一样)特征值和特征向量,构造特征向量空间(这实际上是一个降维的过程);
  • (4)使用传统方法对特征向量聚类,并对应于原始数据的聚类。

  谱聚类算法和传统的聚类方法(例如K-means)比起来有不少优点:

  • K-medoids类似,谱聚类只需要数据之间的相似度矩阵就可以了,而不必像K-means那样要求数据必须是N维欧氏空间中的向量。
  • 由于抓住了主要矛盾,忽略了次要的东西,因此比传统的聚类算法更加健壮一些,对于不规则的误差数据不是那么敏感,而且性能也要好一些。
  • 计算复杂度比K-means要小,特别是在像文本数据或者平凡的图像数据这样维度非常高的数据上运行的时候。

  快速迭代算法和谱聚类算法都是将数据点嵌入到由相似矩阵推导出来的低维子空间中,然后直接或者通过k-means算法产生聚类结果,但是快速迭代算法有不同的地方。下面重点了解快速迭代算法的原理。

2 快速迭代算法的原理

  在快速迭代算法中,我们构造另外一个矩阵W=D−1AW=D−1A ,同第一章做比对,我们可以知道W的最大特征向量就是拉普拉斯矩阵L的最小特征向量。 我们知道拉普拉斯矩阵有一个特性:第二小特征向量(即第二小特征值对应的特征向量)定义了图最佳划分的一个解,它可以近似最大化划分准则。更一般的,k个最小的特征向量所定义的子空间很适合去划分图。 因此拉普拉斯矩阵第二小、第三小直到第k小的特征向量可以很好的将图W划分为k个部分。

  注意,矩阵Lk个最小特征向量也是矩阵Wk个最大特征向量。计算一个矩阵最大的特征向量可以通过一个简单的方法来求得,那就是快速迭代(即PI)。 PI是一个迭代方法,它以任意的向量v0v0作为起始,依照下面的公式循环进行更新。

  在上面的公式中,c是标准化常量,是为了避免vtvt产生过大的值,这里c=||Wvt||1c=||Wvt||1 。在大多数情况下,我们只关心第k(k不为1)大的特征向量,而不关注最大的特征向量。 这是因为最大的特征向量是一个常向量:因为W每一行的和都为1。

  快速迭代的收敛性在文献【1】中有详细的证明,这里不再推导。

  快速迭代算法的一般步骤如下:

  在上面的公式中,输入矩阵W根据W=D−1AW=D−1A来计算。

3 快速迭代算法的源码实现

  在spark中,文件org.apache.spark.mllib.clustering.PowerIterationClustering实现了快速迭代算法。我们从官方给出的例子出发来分析快速迭代算法的实现。

import org.apache.spark.mllib.clustering.{PowerIterationClustering, PowerIterationClusteringModel}
import org.apache.spark.mllib.linalg.Vectors
// 加载和切分数据
val data = sc.textFile("data/mllib/pic_data.txt")
val similarities = data.map { line =>
  val parts = line.split(' ')
  (parts(0).toLong, parts(1).toLong, parts(2).toDouble)
}
// 使用快速迭代算法将数据分为两类
val pic = new PowerIterationClustering()
  .setK(2)
  .setMaxIterations(10)
val model = pic.run(similarities)
//打印出所有的簇
model.assignments.foreach { a =>
  println(s"${a.id} -> ${a.cluster}")
}

  在上面的例子中,我们知道数据分为三列,分别是起始id,目标id,以及两者的相似度,这里的similarities代表前面章节提到的矩阵A。有了数据之后,我们通过PowerIterationClusteringrun方法来训练模型。 PowerIterationClustering类有三个参数:

  • k:聚类数
  • maxIterations:最大迭代数
  • initMode:初始化模式。初始化模式分为RandomDegree两种,针对不同的模式对数据做不同的初始化操作

  下面分步骤介绍run方法的实现。

  • (1)标准化相似度矩阵A到矩阵W
def normalize(similarities: RDD[(Long, Long, Double)]): Graph[Double, Double] = {
    //获得所有的边
    val edges = similarities.flatMap { case (i, j, s) =>
      //相似度值必须非负
      if (s < 0.0) {
        throw new SparkException("Similarity must be nonnegative but found s($i, $j) = $s.")
      }
      if (i != j) {
        Seq(Edge(i, j, s), Edge(j, i, s))
      } else {
        None
      }
    }
    //根据edges信息构造图,顶点的特征值默认为0
    val gA = Graph.fromEdges(edges, 0.0)
    //计算从顶点的出发的边的相似度之和,在这里称为度
    val vD = gA.aggregateMessages[Double](
      sendMsg = ctx => {
        ctx.sendToSrc(ctx.attr)
      },
      mergeMsg = _ + _,
      TripletFields.EdgeOnly)
    //计算得到W , W=A/D
    GraphImpl.fromExistingRDDs(vD, gA.edges)
      .mapTriplets(
        //gAi/vDi
        //使用边的权重除以起始点的度
        e => e.attr / math.max(e.srcAttr, MLUtils.EPSILON),
        TripletFields.Src)
  }

  上面的代码首先通过边集合构造图gA,然后使用aggregateMessages计算每个顶点的度(即所有从该顶点出发的边的相似度之和),构造出VertexRDD。最后使用现有的VertexRDDEdgeRDD, 相继通过fromExistingRDDsmapTriplets方法计算得到最终的图W。在mapTriplets方法中,对每一个EdgeTriplet,使用相似度除以出发顶点的度(为什么相除?对角矩阵的逆矩阵是各元素取倒数,W=D−1AW=D−1A就可以通过元素相除得到)。

  下面举个例子来说明这个步骤。假设有v1,v2,v3,v4四个点,它们之间的关系如下图所示,并且假设点与点之间的相似度均设为1。

  通过该图,我们可以得到相似度矩阵A和度矩阵D,他们分别如下所示。

  通过mapTriplets的计算,我们可以得到从点v1v2,v3,v4的边的权重分别为1/3,1/3,1/3;从点v2v1,v3,v4的权重分别为1/3,1/3,1/3;从点v3v1,v2的权重分别为1/2,1/2;从点v4v1,v2的权重分别为1/2,1/2。 将这个图转换为矩阵的形式,可以得到如下矩阵W

  通过代码计算的结果和通过矩阵运算得到的结果一致。因此该代码实现了W=D−1AW=D−1A 。

  • (2)初始化v0v0

  根据选择的初始化模式的不同,我们可以使用不同的方法初始化v0v0 。一种方式是随机初始化,一种方式是度(degree)初始化,下面分别来介绍这两种方式。

  • 随机初始化
 def randomInit(g: Graph[Double, Double]): Graph[Double, Double] = {
    //给每个顶点指定一个随机数
    val r = g.vertices.mapPartitionsWithIndex(
      (part, iter) => {
        val random = new XORShiftRandom(part)
        iter.map { case (id, _) =>
          (id, random.nextGaussian())
        }
      }, preservesPartitioning = true).cache()
    //所有顶点的随机值的绝对值之和
    val sum = r.values.map(math.abs).sum()
    //取平均值
    val v0 = r.mapValues(x => x / sum)
    GraphImpl.fromExistingRDDs(VertexRDD(v0), g.edges)
  }
  • 度初始化
 def initDegreeVector(g: Graph[Double, Double]): Graph[Double, Double] = {
    //所有顶点的度之和
    val sum = g.vertices.values.sum()
    //取度的平均值
    val v0 = g.vertices.mapValues(_ / sum)
    GraphImpl.fromExistingRDDs(VertexRDD(v0), g.edges)
  }

  通过初始化之后,我们获得了向量v0v0 。它包含所有的顶点,但是顶点特征值发生了改变。随机初始化后,特征值为随机值;度初始化后,特征为度的平均值。

  在这里,度初始化的向量我们称为“度向量”。度向量会给图中度大的节点分配更多的初始化权重,使其值可以更平均和快速的分布,从而更快的局部收敛。详细情况请参考文献【1】。

  • (3)快速迭代求最终的v
for (iter <- 0 until maxIterations if math.abs(diffDelta) > tol) {
      val msgPrefix = s"Iteration $iter"
      // 计算w*vt
      val v = curG.aggregateMessages[Double](
        //相似度与目标点的度相乘
        sendMsg = ctx => ctx.sendToSrc(ctx.attr * ctx.dstAttr),
        mergeMsg = _ + _,
        TripletFields.Dst).cache()
      // 计算||Wvt||_1,即第二章公式中的c
      val norm = v.values.map(math.abs).sum()
      val v1 = v.mapValues(x => x / norm)
      // 计算v_t+1和v_t的不同
      val delta = curG.joinVertices(v1) { case (_, x, y) =>
        math.abs(x - y)
      }.vertices.values.sum()
      diffDelta = math.abs(delta - prevDelta)
      // 更新v
      curG = GraphImpl.fromExistingRDDs(VertexRDD(v1), g.edges)
      prevDelta = delta
    }

  在上述代码中,我们通过aggregateMessages方法计算WvtWvt 。我们仍然以第(1)步的举例来说明这个方法。假设我们以度来初始化v0v0 , 在第一次迭代中,我们可以得到v1(注意这里的v1是上面举例的顶点)的特征值为(1/3)*(3/10)+(1/3)*(1/5)+(1/3)*(1/5)=7/30v2的特征值为7/30v3的特征值为3/10,v4的特征值为3/10。即满足下面的公式。

  • (4)使用k-means算法对v进行聚类
def kMeans(v: VertexRDD[Double], k: Int): VertexRDD[Int] = {
    val points = v.mapValues(x => Vectors.dense(x)).cache()
    val model = new KMeans()
      .setK(k)
      .setRuns(5)
      .setSeed(0L)
      .run(points.values)
    points.mapValues(p => model.predict(p)).cache()
  }

  如果对graphX不太了解,可以阅读spark graph使用和源码解析

4 参考文献

【1】Frank Lin,William W. Cohen.Power Iteration Clustering

【2】漫谈 Clustering (4): Spectral Clustering

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏量子位

头像神器!照片一键秒转简笔画,清华刘永进等CVPR 19 Oral研究 | 在线可玩

清华大学和英国卡迪夫大学的研究人员提出了一种名为APDrawingGAN系统,随意输入一张人脸照片,系统输出黑白人物线条画。

41130
来自专栏数据派THU

谷歌提出新型卷积网络EfficientNet:推理速度提升5.1倍,参数减少88%(附论文&代码)

谷歌提出了一种新型CNN网络EfficientNet,该网络具备极高的参数效率和速度。

23930
来自专栏ATYUN订阅号

Facebook发布PyTorch Hub,一行代码简单重现AI模型

可重现性是许多研究领域的基本要求,包括基于机器学习技术的研究领域。然而,许多机器学习出版物要么不可再现,要么难以复制。

13510
来自专栏小小挖掘机

带答案面经分享-面试中最常考的树模型!

树模型可以说在机器学习的面试中,是面试官最喜欢问的一类问题,今天小编就带你一起回顾一下去年校招中我所经历的树模型相关的问题,这次带答案呦~~(答案是小编的理解,...

17640
来自专栏机器学习与统计学

机器学习是最容易得到错误结论的一种解决方案

机器学习是最容易得到错误结论的一种解决方案。和编程、做表格、或者纯粹的数学建模不同,机器学习是由数据驱动,并有很强的黑箱性。因此很多时候容易得出似是而非的结论。...

12050
来自专栏小小挖掘机

R&S | 手把手搞推荐[0]:我的推荐入门小结

去年年末至今,学习推荐系统已经接近半年,在各种事情的忙碌下,依然坚持着完成一些自学,推荐系统的入门流程逐渐到达尾声,自己学的内容也逐步完善,门是基本入了,所以打...

16440
来自专栏数据派THU

人脸照片秒变艺术肖像画:清华大学提出APDrawingGAN CVPR 2019 oral paper

该项工作被CVPR 2019录取为oral paper。CVPR是计算机视觉和人工智能领域内的国际顶级会议,2019共收到投稿5160篇,录取1300篇,其中o...

15740
来自专栏腾讯高校合作

【犀牛鸟·硬核】腾讯-华中科技大学联合实验室最新研究成果入选SIGMOD国际顶级会议研究类长文

? 前言:腾讯与华中科技大学于2018年成立智能云存储技术联合研究中心,联合研究中心旨在通过强强联合建设一流的智能云存储技术创新和人才培养平台,吸引汇聚顶尖专...

18840
来自专栏小小挖掘机

推荐系统遇上深度学习(五十)-使用强化学习优化用户的长期体验

在现有的推荐模型中,往往优化的目标是点击率,而忽略了用户的长期体验。特别是在信息流推荐中,给用户推荐一个标题很吸引人但内容比较无聊的消息,往往点击率很高,但用户...

45430
来自专栏机器学习与统计学

100天搞定机器学习|Day33-34 随机森林

前言: 随机森林是一个非常灵活的机器学习方法,从市场营销到医疗保险有着众多的应用。它可以用于市场营销对客户获取和存留建模或预测病人的疾病风险和易感性。

11220

扫码关注云+社区

领取腾讯云代金券

年度创作总结 领取年终奖励