(二) 三维点云课程---Spectral Clustering(谱聚类）

三维点云课程—Spectral Clustering

由于GMM需要人工设置聚类数且工作在欧式距离上，对形状有一定的假设(数据点分布近似为椭圆形状)等缺点，Spectral Clustering 谱聚类可以克服这些缺点，对谱聚类的原理进行解析

1.前言

谱聚类分为非归一化和归一化的，分别近似于图论中的RatioCut和NormalizedCut，表达形式如下
$$\begin{gathered} RatioCut(A_1,...,A_k)=\sum _{i=1}^k\frac{cut(A_i,\overline{A_i})}{|A_i|} \\ Ncut(A_1,...,A_k)=\sum _{i=1}^k\frac{cut(A_i,\overline{A_i})}{vol(A_i)} \end{gathered}$$
其中$\overline{A_i}$是$A_i$的补集，对于$cut(A_i,\overline{A_i}),|A_i|,vol(A_i)$解释如下

1.1$cut(A_i,\overline{A_i})$的解释

对于两个子集A，B，子图$G=(V,E)$，使得连接A，B的边的权重($w_{ij}$)最小，权重可以自己定义，可以通过连接性，距离等
$$cut(A,B)=\sum _{i\in A,j\in B}{w}_{ij}$$
对于k个子集，定义如下
$$cut(A_1,...,A_k)=\sum _{i=1}^{k}cut(A_i,\overline{A_i})$$

1.2 $|A_i|,vol(A_i)$的解释

如果仅仅使用上面的公式，可能会出现以下的问题：理论上是想在红色黑色的中间的位置切一刀，但根据上面的公式，可能会在蓝色的点的附件切一刀，因为该位置权重最小。

因此需要增加一些限制，使得切割的$A_i$不要太小，那么怎么定义$A_i$的大小呢？对于非归一化和归一化，追求的目标不同,$A_i$的大小自然不同。
$$\begin{gathered} Size(A):UnNormalized\sim RatioCut\to |A|:A中顶点的数量 \\ Size(A):Normalized\sim NCut\to vol(A)=\sum _{i\in A}{d}_i,d_i=\sum _{j=1}^n{w}_{ij}:点i的权重之和 \end{gathered}$$
在非归一化的谱聚类中，GraphCut追求的目标是不同类之间的点不同，更倾向于一个类中的数据点是均等的；归一化的谱聚类中，GraphCut追求的目标是同一个类之间相似，更倾向于每一个类中的粒度是相等的。在实际工程项目中，多采用归一化的谱聚类，但在这里为了方便推导，采用非归一化的谱聚类进行推导。

下图中黑色的竖线表示非归一化的结果，红色的竖线表示归一化的结果

---

2.原理推导

2.1知识铺垫

其实谱聚类需要三个矩阵，相似矩阵W，对角矩阵D，拉普拉斯矩阵L。现在对这三个矩阵进行解释。

2.1.1相似矩阵W三个建立方法

1.Radius 领域搜索建立。$w_{ij}=d(v_i,v_j)$

2.KNN搜索。$w_{ij}=d(v_i,v_j)$

​ a)建立一个边，如果$v_i$是$v_j$KNN领域的一个邻居，或者$v_j$是$v_i$KNN领域的一个邻居

​ b)建立一个边，如果$v_i$是$v_j$KNN领域的一个邻居，并且$v_j$是$v_i$KNN领域的一个邻居

3.全连接图

2.1.2 对角矩阵D

矩阵D是一个对角矩阵，$d_i=\sum _{j=1}^n{w}_{ij}$,表示，对角线的每一个元素都是相似矩阵W的每一行之和

2.1.2 拉普拉斯矩阵L

非归一化拉普拉斯矩阵$L=D-W$,归一化拉普拉斯矩阵
$$\begin{gathered} L_{sym}=D^{-1/2}L{D}^{1/2}=I-D^{-1/2}W{D}^{-1/2} \\ {L}_{rw}=D^{-1}L=I-D^{-1}W \end{gathered}$$

介绍几个关于拉普拉斯矩阵L的性质

1.对于任意的向量$f\in R^n$,都有
$$f^{T}Lf=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(f_i-f_j{)}^2$$
证明：
$$\begin{array}{l}{f}^{T}Lf=f^{T}Df-f^{T}Wf\\ =\sum _{i=1}^n{f_i}^2{d}_i-\sum _{i=1}^n\sum _{j=1}^n{f}_i{f}_j{w}_{ij}\\ =\frac{1}{2}(\sum _{i=1}^n{f_i}^2{d}_i-2\sum _{i=1}^n\sum _{j=1}^n{f}_i{f}_j{w}_{ij}+\sum _{j=1}^n{f_j}^2{d}_j)\\ =\frac{1}{2}(\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{f}_i^2-2\sum _{i=1}^n\sum _{j=1}^n{f}_i{f}_j{w}_{ij}+\sum _{j=1}^n\sum _{i=1}^n{w}_{ji}{f}_j^2)\\ =\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{(f_i-f_j)}^2\end{array}$$
其中$d_i=\sum _{j=1}^n{w}_{ij}$

2.$L$是对称半正定的矩阵

证明：

对称性：$L^T=(D-W{)}^T=(D^T-W^T)=(D-W)=L$

半正定：$f^{T}Lf=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{(f_i-f_j)}^2\ge =0$

3.L的最小特征值为0，其对应的特征向量为常数向量，方便起见，这里用$1$表示。

证明：
$$\begin{gathered} Lf=(D-W)f=Df-Wf \\ =\left[\begin{array}{c}{d}_1{f}_1\\ ...\\ d_i-f_i\\ ...\\ d_n-f_n\end{array}\right]-\left[\begin{array}{c}{w}_{11}{f}_1+...w_{1n}{f}_n\\ ...\\ w_{i1}{f}_1+...w_{in}{f}_n\\ ...\\ w_{n1}{f}_1+...w_{nn}{f}_n\end{array}\right]=\left[\begin{array}{c}{d}_1{f}_1\\ ...\\ d_i-f_i\\ ...\\ d_n-f_n\end{array}\right]-\left[\begin{array}{c}\sum _{j=1}^n{w}_{1j}{f}_j\\ ...\\ \sum _{j=1}^n{w}_{ij}{f}_j\\ ...\\ \sum _{j=1}^n{w}_{nj}{f}_j\end{array}\right] \\ ={\left[\begin{array}{ccc}...,& d_i{f}_i-\sum _{j=1}^n{w}_{ij}{f}_j,& ...\end{array}\right]}^T=0\bullet f \end{gathered}$$
上式成立的条件是当且仅当f为常量向量，此时L的特征值为0

4.L的特征值是非负的，即$0={\lambda }_1\le {\lambda }_2\le ...\le {\lambda }_n$

证明：通过2，3显然成立

2.2直观理解

情况一：节点0与节点5不相连

左边是10个顶点的连接情况，中间是拉普拉斯矩阵L的特征值分布情况，存在两个0的特征值，右边是两个特征值0对应的特征向量。

情况二：节点0与节点5相连

当节点0与节点5相连时，相似矩阵W，拉普拉斯特征值和特征向量变化情况如下

由于连接了节点0和节点5，原来的零特征值由2个变为1个，但同时也存在一个接近0的特征值0.2984，因为加入的权重是1，试想一下，如果加入的权重比较小，那么特征值就会趋近于0，此时对应的权重就比较小，那么“这条线就是切割分类的线”。

通过以上两种可以得知，有多少个连通域，就有多少个0特征值。一个连通图的特征值至少有一个为0，即至少存在一个独立分区，就是自己本身。如果有两个特征值为0，则存在两个独立分区。并且可以通过0对应的特征向量得知这个类有多少个点属于一类，其中有一个特征向量为常量向量，这里为1，其实为2，3都可以。

2.3 数学推导

2.3.1 k=2的情况

对于两个不相交的子集$A,B\in V$,那么
$$cut(A,B)=\sum _{i\in A,j\in B}{w}_{ij}$$
那么对于非归一化和归一化
$$\begin{gathered} RatioCut(A_1,...,A_k)=\sum _{i=1}^k\frac{cut(A_i,\overline{A_i})}{|A_i|} \\ Ncut(A_1,...,A_k)=\sum _{i=1}^k\frac{cut(A_i,\overline{A_i})}{vol(A_i)} \end{gathered}$$
方便起见，所有推导基于非归一化。谱聚类的问题就是最小化RatioCut函数，即
$$\underset{A\subset V}{\min }RatioCut(A,\overline{A})=\underset{A\subset V}{\min }(\frac{cut(A,\overline{A})}{|A|}+\frac{cut(\overline{A},A)}{|\overline{A}|})$$

给予一个子集$A\subset V$，构建一个向量$f=[f_1,...,f_n{]}^T\in R^n$,构建如下(假设$f$已经成功构建出来的，实际上$f$是不知道长什么的)
$$f_i=\{\begin{array}{l}\sqrt{|\overline{A}|/|A|}\mathrm{i}\mathrm{f}{v}_i\in A\\ \sqrt{|A|/|\overline{A}|}\mathrm{i}\mathrm{f}{v}_i\in \overline{A}\end{array}$$
其实此时就可以通过$f_i$的正负就可以判断$v_i$是否属于A类了
$$\{\begin{array}{l}{v}_i\in A\mathrm{i}\mathrm{f}{f}_i\ge 0\\ v_i\in \overline{A}\mathrm{i}\mathrm{f}{f}_i<0\end{array}$$

那么

其中$|V|$表示原来没有切的图中还有多少个点。


另外，$f$是垂直于常量向量，且$||f||=\sqrt{n}$

证明如下
$$\begin{gathered} f^{T}1=\sum _{i=1}^n{f}_i=\sum _{i\in A}\sqrt{\frac{|\overline{A}|}{|A|}}-\sum _{i\in \overline{A}}\sqrt{\frac{|A|}{|\overline{A}|}}=|A|\sqrt{\frac{|\overline{A}|}{|A|}}-|\overline{A}|\sqrt{\frac{|A|}{|\overline{A}|}}=0 \\ ||f|{|}^2=\sum _{i=1}^n{f_i}^2=|A|\frac{|\overline{A}|}{|A|}+|\overline{A}|\frac{|A|}{|\overline{A}|}=|A|+|\overline{A}|=n \end{gathered}$$
因此原先最小化图切问题转化为
$$\underset{A\subset V}{\min }{f}^{T}Lf,s.t.,f\perp 1,||f||=\sqrt{n},f_i=\{\begin{array}{l}\sqrt{|\overline{A}|/|A|}\mathrm{i}\mathrm{f}{v}_i\in A\\ \sqrt{|A|/|\overline{A}|}\mathrm{i}\mathrm{f}{v}_i\in \overline{A}\end{array}$$
由于$f_i$是假设的，因此上式近似为
$$\underset{A\subset V}{\min }{f}^{T}Lf,s.t.,f\perp 1,||f||=\sqrt{n}$$
遇见对称矩阵$A$，求解类似$x^{T}Ax$的最大最小值，采用Rayleigh商进行求解。由于Rayleigh商定理本身没有条件，但是对于上式存在$f\perp 1,||f||=\sqrt{n}$限制条件，那该怎么处理呢？对于$L$的最小特征值就是$min{f}^{T}Lf$的最小值，但是$f$是有条件的，好在$L$的最小特征值0对应的特征向量为常数，又因为常数不垂直于1，取$L$矩阵的第二小特征值作为$f^{T}Lf$的最小值，对应的特征向量就是$f$。

现在存在一个问题，就是怎么通过上述得到的$f$进行分类呢，这时肯定就有人说使用
$$f_i=\{\begin{array}{l}\sqrt{|\overline{A}|/|A|}\mathrm{i}\mathrm{f}{v}_i\in A\\ \sqrt{|A|/|\overline{A}|}\mathrm{i}\mathrm{f}{v}_i\in \overline{A}\end{array}$$
通过$f_i$的正负来进行判断，但是上式假设出来的。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Fm5mXbsP-1635251405573)(E:\资料\三维重建课程\第四课时\图片\理想和非理想.png)]

左图是已知$f_i$的分布，右图$f_i$是未知的，且黄色的点是上述得到的特征向量$f_2$的点。那么此时需要将得到的$f_2$和常量向量$f_1$组合起来，采用KMeans方法进行分类。那么为什么需要结合$f_1$,只采用$f_2$不就可以了吗？如果只存在两个孤立的连通域，只用$f_1$就可以了，但是如果因为两个连通域之间有连线，只用$f_1$就不可以了。其实$f_1$表示的是连通域的信息，$f_2$是图的构造。

2.3.2 k$\ge 2$的情况

对于k个分类，$A_1,...,A_k$，那么
$$cut(A_1,...,A_k)=\sum _{i=1}^{k}cut(A_i,{\overline{A}}_i)$$
构造一个Indication矩阵$H\in R^{n\times k}$
$$h_{ij}=\{\begin{array}{l}1/\sqrt{|A_i|}\mathrm{i}\mathrm{f}i\in A_j\\ 0otherwise\end{array}$$
举个例子

其中 { 1 , 2 , 3 , 4 } ∈ A 1 ; 5 ∈ A 2 ; 6 ∈ A 3 \{1,2,3,4\} \in A_1;5 \in A_2;6\in A_3 {1,2,3,4}∈A1​;5∈A2​;6∈A3​。

类比k=2的情况，我们有
$${h_i}^{T}L{h}_i=\frac{cut(|A_i|,\overline{|A_i|})}{|A_i|},{h_i}^{T}L{h}_i=(H^{T}LH{)}_{ii}$$
那么
$$\begin{gathered} RatioCut(A_1,...,A_k)=\sum _{i=1}^k\frac{cut(|A_i|,\overline{|A_i|})}{|A_i|} \\ =\sum _{i=1}^k{h}_i^{T}L{h}_i=\sum _{i=1}^k(H^{T}LH{)}_{ii}=Tr(H^{T}LH) \end{gathered}$$
现在问题转化为
$$\underset{A_1,...,A_k}{\min }Tr(H^{T}LH)s.t,H^{T}H=I,h_{ij}=\{\begin{array}{l}1/\sqrt{|A_i|}\mathrm{i}\mathrm{f}i\in A_j\\ 0otherwise\end{array}$$
近似为
$$\underset{A_1,...,A_k}{\min }Tr(H^{T}LH).s.t,H^{T}H=I$$
对于矩阵的Rayleigh商，现在给出结论$L$前k个较小的特征值对应的特征向量组合起来就是$H$，具体推导过程可以参考之前PCA推导.

---

3.谱聚类的步骤及效果

3.1步骤

3.1.1非归一化步骤

• 通过图建立一个相似矩阵$W\in R^{n\times n}$,(相似矩阵建立三选一)
• 计算对角矩阵$D$,其中$d_{ij}=\sum _{j=1}^n{w}_{ij}$
• 计算非归一化的拉普拉斯矩阵$L=D-W$
• 计算$L$最小的前k个特征值对应的特征向量$v_1,...,v_k$
• 将上述的特征向量按照行的方向进行排列，组合成$V\in R^{n\times k}$
• 将$V$的每一行数据，记为 Y = { y i , . . . , y n } Y=\{y_i ,...,y_n\} Y={yi​,...,yn​}
• 对上述的$Y$进行KMeans操作，得到 C = { C 1 , . . . , C k } C=\{C_1,...,C_k\} C={C1​,...,Ck​}的聚类结果
• 通过上述$C$对应的索引，将原始数据进行分类

3.1.1归一化步骤

• 通过图建立一个相似矩阵$W\in R^{n\times n}$,(相似矩阵建立三选一)
• 计算对角矩阵$D$,其中$d_{ij}=\sum _{j=1}^n{w}_{ij}$
• 计算归一化的拉普拉斯矩阵$L_{rw}=D^{-1}(D-W)$
• 计算$L_{rw}$最小的前k个特征值对应的特征向量$v_1,...,v_k$
• 将上述的特征向量按照行的方向进行排列，组合成$V\in R^{n\times k}$
• 将$V$的每一行数据，记为 Y = { y i , . . . , y n } Y=\{y_i ,...,y_n\} Y={yi​,...,yn​}
• 对上述的$Y$进行KMeans操作，得到 C = { C 1 , . . . , C k } C=\{C_1,...,C_k\} C={C1​,...,Ck​}的聚类结果
• 通过上述$C$对应的索引，将原始数据进行分类

3.2效果

下图即谱聚类的分类的结果，不同GMM和KMeans，谱聚类可以轻易的将圆分开，因为它是工作在图上的，而不是欧式聚类

3.3 谱聚类自动判断聚类数

将拉普拉斯矩阵L的特征值进行从小到大排布，并且记${\Delta }_k=|{\lambda }_k-{\lambda }_{k-1}|$,如果${\Delta }_k$突变，表示特征值由小到大突变，对应于下图的渐变的过程，此时只需要关心突变前k个特征值，此时k=聚类的个数。

---

4.谱聚类的优、缺点

4.1 优点

• 不会对任何性质进行假设
• 自动的发觉右多少类
• 可以工作在高维空间上
• 工作在图上，并不是欧式距离

4.2 缺点

算法的时间复杂度为$O(n^3)$，效率较低。

注意具体代码参考下一篇，嘻嘻嘻嘻