谱聚类的R解释

注，有疑问加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. 参考文献与推荐阅读

最后编辑于：2018.11.09 15:46:15

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 199,711评论 5赞 468
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 83,932评论 2赞 376
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 146,770评论 0赞 330
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 53,799评论 1赞 271
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 62,697评论 5赞 359
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,069评论 1赞 276
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,535评论 3赞 390
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,200评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,353评论 1赞 294
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,290评论 2赞 317
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,331评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,020评论 3赞 315
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,610评论 3赞 303
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,694评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,927评论 1赞 255
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,330评论 2赞 346
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 41,904评论 2赞 341