谱聚类的R解释

作者: 飘舞的鼻涕 | 来源:发表于2017-11-02 14:00 被阅读0次

谱聚类的R解释
谱聚类算法总结
待完成：scikit 聚类方法可用源代码
14 聚类算法 - 代码案例六- 谱聚类(SC)算法案例
Clustering
Day 684：机器学习笔记（13）
谱聚类
谱聚类
谱聚类
谱聚类

注，有疑问加QQ群..[174225475].. 共同探讨进步
有偿求助请出门左转 door , 合作愉快

1. 相关概念科普

度矩阵

度矩阵为对角矩阵，对角线上的值表示与该点有来往的其他点的个数即度为与节点相连的个数

邻接矩阵

邻接矩阵表示图上各点之间是否有联系或者联系程度, 可以是 1/0 也可以是具体权值

laplas1

laplas2

拉普拉斯矩阵

图论中的的拉普拉斯矩阵(也被称为导纳矩阵/吉尔霍夫矩阵/离散拉普拉斯)是图的矩阵表示. 它是度矩阵和邻接矩阵的差. 拉普拉斯矩阵结合吉尔霍夫理论可以用来计算图的最小生成树的个数。
拉普拉斯矩阵还可用来寻找图的其他属性: 谱图理论(spectral graph theory)
黎曼几何的Cheeger不等式有涉及了拉普拉斯矩阵的离散模拟. 这或许是谱图理论中最重要的定理也是在算法应用中最有用的facts.它通过拉普拉斯矩阵的第二特征值来近似图的最小割
谱的定义 就是：
方阵作为线性算子，它的所有特征值的全体统称方阵的谱；
方阵的谱半径为最大的特征值；
矩阵A的谱半径：( $A^{T} \cdot A$ )的最大特征值。
其实，这里谱的本质是伪逆，是 $SVD$ 中奇异值的平方

2. 谱聚类的基本思想

2.1 聚类与图论

所谓聚类(Clustering), 就是要把一堆样本合理地分成两份或者 $K$ 份. 从图论的角度来说,聚类的问题就相当于一个图的分割问题. 即给定一个图 $G = (V, E)$ , 顶点集 $V$ 表示各个样本,带权的边表示各个样本之间的相似度,谱聚类的目的便是要找到一种合理的分割图的方法,使得分割后形成若干个子图,连接不同子图的边的权重(相似度)尽可能低,同子图内的边的权重(相似度)尽可能高.被割掉各边的权值和越小代表被它们连接的子图之间的相似度越小,隔得越远. 物以类聚,人以群分,相似的在一块儿,不相似的彼此远离
至于如何把图的顶点集分割/切割为不相交的子图有多种办法,如

cut/Ratio Cut
Normalized Cut
不基于图，而是转换成SVD能解的问题

2.2 最小割

在将数据从点集的结构转换为了图的结构后，我们可以引入更多图论中的概念来帮助理解聚类的本质. 比如说最小割的概念(Minimum cut)
最小割是指去掉图中的一些带权边，在使得图从联通图变为不联通图的前提下，尽可能的让去掉的边权值之和尽可能小。对于数据的相似性图来说，最小割就是要去除一些很弱的相似性，把数据点从一个互相连接的整体分为两个或者多个独立的部分，这其实就是一个聚类的过程。去除的权值就是聚类后不同组别间的相似程度，我们希望这个值尽可能的小，也就是希望不同组别之间的差距尽可能的大
不过，仅仅是用最小割来聚类存在着一些缺陷，因为我们只考虑了不同类的点之间相似度要小，而一个好的聚类还要具有同一类点相似度大的特征。比如下图就是一个例子，一个好的聚类应该是在绿线处分割整个图，但使用最小割却会割在红色的位置

specc3

如下图所示,每个顶点是一个样本,共有7个样本(实际中一个样本是特征向量). 边的权值就是样本之间的相似度. 然后我们需要分割,分割之后要使得连接不同组之间的边的权重尽可能低(组间相似度小),组内部的边的权重尽可能高(组内相似度高)

specc1

比如我们把中间的边去掉,就满足上面的簇间相似性最小，簇内相似性最大。如下图

specc2

3. 算法步骤

谱聚类算法主要有下面几部分：

3.1 未正则拉普拉斯矩阵的谱聚类算法

3.1.1 计算得到图的邻接矩阵 $W$ ，以及拉普拉斯矩阵 $L$ ；

给定样本的原始特征,我们需要计算两两样本之间的相似度值,才能构造出邻接矩阵 $W$ . 我们一般使用高斯核函数来计算相似度。公式如下：

$S(x,y)=e^{-\frac{||x-y||^{2}}{2\sigma^{2}}}$
即
$S(x,y)=e^{-{\frac{\sqrt{\sum_{i=1}^{k}\left(x_{i}-y_{i}\right)^{2}}}{2\sigma^2}}}$

如果样本 $x,y$ 是对应数据集 $dt$ 的第 $i,j$ 行,则上述公式中的

$\begin {aligned} \sqrt{\sum_{i=1}^{k}\left(x_{i}-y_{i}\right)^{2}}&= \sqrt{\sum_{i,j=1}^{n}(dt[x,]-dt[y,])^2} \\\ &= \sqrt{\sum_{i,j=1}^{n}(dt[i,]-dt[j,]) * t(dt[i,]-dt[j,])} \\\ &=dist(dt) \end {aligned}$

高斯核相似度公式用 R 来表达:

S(dt) = exp(-(dist(dt)^2)/(2*sigma^2))

其中 $x$ 和 $y$ 是两个样本(行向量), $x_{i}$ $y_{i}$ 是样本 $x,y$ 对应的第 $i$ 个原始特征向量，我们可以计算出它们二者之间的相关性(高斯核). 注意,邻接矩阵的对角线元素一致等于零. 这样我们就得到了邻接矩阵 $W$

拉普拉斯矩阵 $L = D - W$ , 其中 $D$ 为上面图的度矩阵，度是图论中的概念，也就是矩阵E行或列的元素之和。后面我们就要基于矩阵 $L$ 来进行研究

3.1.2 聚类划分准则

关于划分准则的问题是聚类技术中的核心。谱聚类被看作属于图论领域的问题，那个其划分也会和边的权重有关

准则1：最小化进行分割时去掉的边的权重之和
$cut (A_{1} \ldots A_{k}) = \frac {1}{2} \cdot \sum_{i=1}^{k} W (A_{i},\overline{A_{i})}$
这个准则直观地让分割之后簇之间的相似度最小了。但是有个问题，就是这个准则没有考虑到簇的大小，容易造成一个样本就能分为一个簇。为了避免这个问题，有了下面的准则
准则2：考虑簇的大小
$RationCut\left(A_{1},\ldots ,A_{k} \right)=\frac{1}{2} \cdot \sum_{i=1}^{k} \frac{W\left(A_{i},\overline{A_{i}}\right)} {||A_{i}||}$

准则2相比于准则1的改进，类似于 $C4.5$ 中的增益比率和 $ID3$ 的信息增益的改进。在 $RatioCut$ 中，如果某一组包含的顶点数越少，那么它的值就越大。在一个最小化问题中，这相当于是惩罚，也就是不鼓励将组分得太小。现在只要将最小化 $RatioCut$ 解出来，分割就完成了。不幸的是，这是个 $NP$ 难问题。想要在多项式时间内解出来，就要对这个问题作一个转化了。在转化的过程中，就用到上面提到的L的那一组性质，经过若干推导，最后可以得到这样的一个问题
$\begin {aligned} \min_{H \in \mathbb{R}^{n \cdot k}} (Tr(H^{T}LH)) \\\ s.t. H^{T}H=I \end {aligned}$
但到了这关键的最后一步，咱们却遇到了一个比较棘手的问题，即由之前得到的拉普拉斯矩阵的性质 $L$ 最小特征值为 $0$ ，对应的特征值为 $1$ .其不满足 $f$ 与 $1$ 正交的条件,那该怎么办呢？根据论文“A Tutorial on Spectral Clustering”中所说的Rayleigh-Ritz 理论，我们可以取第 $2$ 小的特征值，以及对应的特征向量 $v$
更进一步，由于实际中，特征向量 $v$ 里的元素是连续的任意实数，所以可以根据 $v$ 是大于0，还是小于0对应到离散情况下的 $f=(f_{1},\ldots ,f_{n}) \prime \in \mathbb{R}^{n\cdot k}$ ,决定 $f$ 是取 $\sqrt{\frac{|A|} {|\overline{A}|}}$ 还是取 $-\sqrt{\frac{|A|} {|\overline{A}|}}$ .而如果能取 $v$ 的前 $k$ 个特征向量,进行 $Kmeans$ 聚类,得到 $K$ 个簇,便从二聚类扩展到了 $K$ 聚类的问题 ( $\overline{A}$ 表示A的补集)
而所要求的这前 $K$ 个特征向量就是拉普拉斯矩阵的特征向量(计算拉普拉斯矩阵的特征值,特征值按照从小到大顺序排序,特征值对应的特征向量也按照特征值递增的顺序排列,取前 $K$ 个特征向量,便是我们所要求的前 $K$ 个特征向量)!
所以，问题就转换成了：求拉普拉斯矩阵的前 $K$ 个特征值，再对前 $K$ 个特征值对应的特征向量进行 $Kmeans$ 聚类。而两类的问题也很容易推广到 $K$ 类的问题，即求特征值并取前 $K$ 个最小的,将对应的特征向量排列起来,再进行 $Kmeans$ 聚类. 两类分类和多类分类的问题，如出一辙
$RatioCut$ 巧妙地把一个NP难度的问题转换成拉普拉斯矩阵特征值(向量)的问题，将离散的聚类问题松弛为连续的特征向量，最小的系列特征向量对应着图最优的系列划分方法。剩下的仅是将松弛化的问题再离散化，即将特征向量再划分开，便可以得到相应的类别。不能不说妙哉！

3.1.3 求出 $L$ 的前 $k$ 个特征值及对应的特征向量

在本文中,除非特殊说明,否则“前 $k$ 个”指按照特征值的大小从小到大的顺序)以及对应的特征向量
PS： $L=D-W$ 是计算前 $k$ 小的特征向量，如果是 $L=W-D$ 则是计算前 $k$ 大

3.1.4 根据新生成的特征矩阵进行 $kmeans$ 聚类

把这 $k$ 个特征(列)向量排列在一起组成一个 $n \cdot k$ 的矩阵,将其中每一行看作 $k$ 维空间中的一个向量,并使用 $Kmeans$ 算法进行聚类. 聚类的结果中每一行所属的类别就是原来 Graph 中的节点亦即最初的 $n$ 个数据点分别所属的类别

3.1.5 总结一下完整的谱聚类算法的流程

确认输入参数：一个描述数据相似性图的邻接矩阵 $W$ ，和聚类的目标类数目 $k$
计算度矩阵 $D$ 和拉普拉斯矩阵 $L$
计算特征值方程 $L\cdot y = \lambda \cdot D \cdot y$ 的前 $k$ 个特征向量，记为 $y_{1} \ldots y_{k}$
通过 $y_{1} \ldots y_{k}$ 按列排成矩阵 $Y \in \mathbb{R}^{n \cdot k}$
$Y$ 的第 $i$ 行的行向量记为 $x_{i} \prime (i \leq n)$
对 $x_{i} \prime$ 进行 $kmeans$ 聚类
将 $x_{i} \prime$ 的类别分配给 $x_{i}$ , 即 $class(x_{i})=class(x_{i} \prime)$
输出聚类结果
整个过程都只用到了包括求解特征值等基本线性代数运算，因此谱聚类是一种推导有些繁杂，但是实现很简单的聚类算法

3.2 正则的拉普拉斯矩阵的谱聚类算法

3.2.1 对称拉普拉斯矩阵的谱聚类算法

将为正则拉普拉斯中的 $L = D -W$ 替换成 $L=D^{-\frac{1}{2}}\cdot(D-W)\cdot D^{\frac{1}{2}}$
同时在完成 $k$ 维特征向量提取之后,进行 $kmeans$ 聚类之前,添加一步:
将 $k$ 维特征向量对应的各样本进行单位化使得 $|y_{i}|=1$ 即
$mx[i,]=\frac{mx[i,]} {\sqrt{\sum_{i=1}^{k}mx[i,]^2}}$

mx2 <- mx1/sqrt(rowSums(mx1^2))

然后在此基础上继续 $kmeans$ 聚类

下面是在参加某赛事时正则化谱聚类实现步骤

library(data.table)
library(dplyr)
library(ggplot2)
library(geosphere)
setwd('xxx/spec')
train1 <- fread('train.csv',header = FALSE,encoding = 'UTF-8')
test1 <- fread('test.csv',header = FALSE,encoding = 'UTF-8')
dt_all <- rbind(train1,test1)

dist2 = distm(dt_all[,c('longitude','latitude')])/1000
lamb1 <- sqrt(ncol(dt_all[,c('longitude','latitude')]))/sqrt(2)
W1 <- exp(-dist2^2/(2*lamb1^2))
rm(dist2)
W1[W1<1e-15]=0
W2 <- as(W1,'sparseMatrix')
diag1 <- rowSums(W1)
rm(W1)
D1 <- Matrix::sparseMatrix(i=seq(1,length(diag1)),
                           j=seq(1,length(diag1)),
                           x=diag1) # 度矩阵
D2 <- Matrix::sparseMatrix(i=seq(1,length(diag1)),
                           j=seq(1,length(diag1)),
                           x=diag1^(-1/2))
D3 <- Matrix::sparseMatrix(i=seq(1,length(diag1)),
                           j=seq(1,length(diag1)),
                           x=diag1^(1/2))
L1 <- D1-W2
L2 <- D2%*%L1%*%D3 #拉普拉斯矩阵
rm(L1,D1,D2,D3)
eig1 <- eigen(L2)
mx1 <- eig1$vectors[,-c(1:(ncol(eig1$vectors)-k))]
mx2 <- mx2/sqrt(rowSums(mx1^2))
# find best k value for kmeans
library(factoextra)
fviz_nbclust(area1, kmeans, method = "silhouette",k.max = 50)
# the larger the better for result value(this formu is friendly for bigdata)
kmean1 <- kmean(mx2,49)
# plot
ggplot()+
  geom_point(data=dt_all,aes(x=longitude,y=latitude,color=as.factor(kmean1$cluster)))+
  geom_text(data=NULL,aes(x=kmean1$centers[,1],
                          y=kmean1$centers[,2],
                          label=row.names(kmean1$centers)))+
  theme(legend.position = 'none')

由于样本huge,eigen时cost too much, so kmeans替代，以上代码作为R实现步骤留念，指不定哪天R真正解决了分布式问题,谱聚类就能跟kmeans一样高效了
Kernlab包中有specc函数可实现此功能,但碍于样本量巨大,2个小时没出结果也只能作罢

3.2.2 随机游走拉普拉斯矩阵的谱聚类算法

它的变化就是将 $L$ 矩阵的 $L=D-W$ 计算方式改为 $L=D^{-1}(D-W)$ ,其他步骤与未正则拉普拉斯矩阵的谱聚类算法相同

4. 参考文献与推荐阅读

谱聚类的R解释
注，有疑问加QQ群..[174225475].. 共同探讨进步有偿求助请出门左转 door , 合作愉快 1....
谱聚类算法总结
聚类三种方法：k-means聚类、密度聚类、层次聚类和谱聚类Spectrum Clustering 简述谱聚类是...
待完成：scikit 聚类方法可用源代码
kmeans聚类 spec 谱聚类
14 聚类算法 - 代码案例六- 谱聚类(SC)算法案例
13 聚类算法 - 谱聚类需求使用scikit的相关API创建模拟数据，然后使用谱聚类算法进行数据聚类操作，并...
Clustering
本文结构安排经典聚类算法：线性聚类 Kmeans 经典聚类算法：非线性聚类 DBSCAN、谱聚类新兴聚类算法：...
Day 684：机器学习笔记（13）
谱聚类 KMeans需要事先确定有多少簇，谱聚类可以不需要事先指定。基于图切割的谱聚类算法分两个主要步骤：图切割...
谱聚类
Spectral Clustering 和传统的聚类方法（例如 K-means）比起来有不少优点：和K-medo...
谱聚类
原理：谱聚类(Spectral Clustering, SC)是一种基于图论的聚类方法——将带权无向图划分为两个或...
谱聚类
谱聚类原理介绍： https://www.cnblogs.com/pinard/p/6221564.html 谱聚...
谱聚类
先收藏下，数学不好的我，还要再看看谱聚类（spectral clustering）是广泛使用的聚类算法，比起传统...