【机器学习】谱聚类原理推导

一、前言

本文将介绍一种较为突出的聚类算法——谱聚类。这是一种基于图模型的算法。有过图论或者数据结构基础的人估计会相对来说觉得容易。能力有限，推导并不严谨，请见谅。代码谱聚类代码实现

二、相关知识

2.1、无向图

无向图，由两个要素构成——节点和边。以上面的无向图为例

节点，我们定义为 V = { A , B , C , D , E } V=\{A,B,C,D,E\} V={A,B,C,D,E}。边定义为 E = { w 1 , w 2 , w 3 , w 4 , w 5 } E=\{w_1,w_2,w_3,w_4,w_5\} E={w1​,w2​,w3​,w4​,w5​}对应图中的w1…。

所谓边，就是两个点之间的关联性，一般情况下，关联性越高，所连接的边$w$就越高。

在谱聚类中，我们的一个样本点就是一个节点，那么如果我们有n个样本，那就有n个节点。样本之间的相关性就是用$w$来衡量。既然是聚类，那自然就是相关性越高，则越能聚成一类。

2.2、邻接矩阵

前面我们讲到，某个节点与某个节点之间连接的值用$w$来表示。那有没有一种东西，能够充分表达所有点之间的关系？有，就是邻接矩阵。比如图中就可以表示为（无连接我们暂时用0表示）
$$W={\left[\begin{array}{cccccc}& A& B& C& D& E\\ A& 0& w_1& w_2& 0& 0\\ B& w_1& 0& 0& w_3& 0\\ C& w_2& 0& 0& 0& w_4\\ D& 0& w_3& 0& 0& w_5\\ E& 0& 0& 0& w_4& w_5\end{array}\right]}_{5\times 5}$$
比如节点A，它与B，C有连接，所以在B，C列有值，其余都为0。以此类推。

如果有n个样本，那么它就是一个$n\times n$的方阵。所以我们可以表示成这样
$$W={\left[\begin{array}{cccc}{w}_{11}& w_{12}& \cdots & w_{1n}\\ w_{21}& w_{22}& \cdots & w_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ w_{n1}& w_{n2}& \cdots & w_{nn}\end{array}\right]}_{n\times n}$$

2.3、度

度，意为某个节点所有相连的边的总和。比如图中的节点A，它的度$d=w_1+w_2$。如果有n个样本，那么第$i$个样本的度可以表示为
$$d_i=\sum _{j=1}^n{w}_{ij}$$
其实就等于邻接矩阵的第$i$行求和

那我们又该如何表示所有的节点的度呢？我们用$D$表示
$$D={\left[\begin{array}{cccc}{d}_1& 0& 0& \cdots \\ 0& d_2& 0& \cdots \\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & d_n\end{array}\right]}_{n\times n}$$
其实就是一个对角矩阵（主对角线有值，其余都为0）。所以每个样本的都对应矩阵$D$中的一个度。

2.4、核函数

我们如何去构造邻接矩阵$W$，构造点与点之间的联系？

$\boxed{我们仅介绍一种较为普遍的方法}$——全连接法，即两个点中间的关系表示为
$$w_{ij}=\exp \left\{-\frac{||x_i-x_j|{|}_2^2}{2{\sigma }^2}\right\}$$
分子处为L2范数，其实这就是高斯核函数的定义。同样的，你也可以选择其他核函数。

通过此法，我们便可以构造出邻接矩阵$W$，由$W$可得出矩阵$D$。

2.5、拉普拉斯矩阵

定义$L=D-W$。由这个算出来的矩阵$L$被称为拉普拉斯矩阵。

拉普拉斯矩阵的性质：

①拉普拉斯矩阵是半正定矩阵。

②特征值中0出现的次数就是图连通区域的个数。

③最小特征值是0，因为拉普拉斯矩阵每一行的和均为0。

并且，拉普拉斯矩阵还有一个性质——$\boxed{对任意一个n维y向量，有y^{T}Ly=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_i-y_j{)}^2}$

证明：
$$\begin{array}{rl}{y}^{T}Ly=& y^{T}Dy-y^{T}Wy\\ =& \sum _{i=1}^n{d}_i{y}_i^2-\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{y}_i{y}_j\\ =& \frac{1}{2}\left(\sum _{i=1}^n{d}_i{y}_i^2-2\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{y}_i{y}_j+\sum _{j=1}^n{d}_j{y}_j^2\right)\\ =& \frac{1}{2}\left(\sum _{i=1}^n\sum _{j=1}^n{w}_{{}_{ij}}{y}_i^2-2\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}{y}_i{y}_j+\sum _{j=1}^n\sum _{i=1}^n{w}_{{}_{ji}}{y}_j^2\right)\\ =& \frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_i-y_j{)}^2\end{array}\text{\hspace{1em}(1)(2)(3)(4)(5)}$$
式（3）到式（4）用到$d_i=\sum _{j=1}^n{w}_{ij}$。式（4）到式（5）用到平方差公式。

三、原理推导

3.1、目标

一个节点就是一个样本。我们要在这张无向图中切一刀，然后得到两部分。比如上图中，用虚线切开，得到两部分。我们就认为上半部分为一类。下半部分为一类。

显然，如何砍这一刀是至关重要的。我们定义砍一刀的代价$\boxed{目标是砍一刀的代价最小}$
$$W(A,B)=\sum _{i\in A}\sum _{j\notin B}{w}_{ij}$$
比如按照上面的无向图，上半部分记为A集合，下半部分记为B集合。那么代价就是a节点与B集合所有点的w（仅有w2，其余为0）加上b节点与B集合所有点的w（仅有w3，其余为0）。所以代价就是$w_2+w_3$。

以此类推，如果要聚成k类，定义代价函数（$\frac{1}{2}$是因为不同类别之间会重复计算代价，故如此）
$$Cut(A_1,A_2,\cdots ,A_k)=\frac{1}{2}\sum _{i=1}^{k}W(A_i,\overline{A_i})$$
$\overline{A_i}$表示不属于第i类的节点集合。

然而，这种形式的代价有问题，求解的时候算法会趋向于将权重w小的单个节点作为一类，这是我们很不想要的。为此我们必须做出整改。

有两种方案：$\boxed{①RatioCut}$、$\boxed{NCut}$。

3.2、RatioCut

定义$|A_i|$表示$i$集合的节点个数。比如上图中的上半部分集合A，只有节点$a,b$，我们就得到$|A|=2$。

我们将代价函数变成这样
$$RatioCut(A_1,A_2,\cdots ,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}$$
也就是要最小化这个代价函数
$$\arg \underset{A_i}{\min }RatioCut(A_1,A_2,\cdots ,A_k)=\arg \underset{A_i}{\min }\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}$$
这种用集合来表达的形式我们是没办法去求解的。因此，我们换一种表达方式。对于每一个样本，我们都定义一个k维的指示向量。
y = { 1 ∣ A i ∣ j ∈ A i 0 j ∉ A i j ∈ { 1 , 2 , ⋯   , k } y=\left\{\begin{matrix} \frac{1}{\sqrt|A_i|} & j\in A_i \\ 0 & j\notin A_i &j\in\{1,2,\cdots,k\} \end{matrix}\right. y={∣ ​Ai​∣1​0​j∈Ai​j∈/Ai​​j∈{1,2,⋯,k}​
意思是如果这个样本属于第$i$类，那么对应位置就有值，否则为0。

比如对于第1个样本，假设它属于第二类：第二个样本属于第三类。则有
$$y_1={\left(\begin{array}{c}0\\ \frac{1}{\sqrt{|}{A}_2|}\\ 0\\ ⋮\\ 0\end{array}\right)}_{k\times 1};y_2={\left(\begin{array}{c}0\\ 0\\ \frac{1}{\sqrt{|}{A}_3|}\\ ⋮\\ 0\end{array}\right)}_{k\times 1};$$
以此类推，所有的样本，用$Y$表示
$$Y=\left(\begin{array}{c}{y}_1^T\\ y_2^T\\ ⋮\\ y_n^T\end{array}\right)={\left[\begin{array}{cccc}{y}_{11}& y_{12}& \cdots & y_{1k}\\ y_{21}& y_{22}& \cdots & y_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ y_{n1}& y_{n2}& \cdots & y_{nk}\end{array}\right]}_{n\times k}$$
所以我们的目标就是
$$\arg \underset{Y}{\min }\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}$$
同样的，里面的目标函数我们也要换成Y的表达形式

对于矩阵$Y$，第$m$列中不为0的部分，就是属于第$m$类的样本。我们用$y_m$表示$Y$的第$m$列
$$\begin{array}{rl}{y}_m^{T}L{y}_m=& \frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_{mi}-y_{mj}{)}^2\\ =& \frac{1}{2}\left(\sum _{i\in A_m}\sum _{j\notin A_m}{w}_{ij}(\frac{1}{\sqrt{|A_m|}}-0{)}^2+\sum _{i\notin A_m}\sum _{j\in A_m}{w}_{ij}(0-\frac{1}{\sqrt{|A_m|}}{)}^2+\sum _{i\in A_m}\sum _{j\in A_m}{w}_{ij}(0-0{)}^2\right)\\ =& \frac{1}{2}\left(\sum _{i\in A_m}\sum _{j\notin A_m}{w}_{ij}\frac{1}{|A_m|}+\sum _{i\notin A_m}\sum _{j\in A_m}{w}_{ij}\frac{1}{|A_m|}\right)\\ =& \frac{1}{2}\left(\frac{Cut(A_m,{\bar{A}}_m)}{|A_m|}+\frac{Cut(A_m,{\bar{A}}_m)}{|A_m|}\right)\\ =& \frac{Cut(A_m,{\bar{A}}_m)}{|A_m|}\end{array}\text{\hspace{1em}(1)(2)(3)(4)}$$
式（1）用到了上面拉普拉斯矩阵的性质。式（2）用到了指示向量。式（3）到式（4）用到了代价函数

发现了吗，它刚好等于我们第$m$类的代价。

我们再将原目标函数转化成矩阵的形式
$$\begin{array}{rl}\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}=& \sum _{i=1}^k\frac{Cut(A_i,{\bar{A}}_i)}{|A_i|}\\ =& \mathbb{T}r{\left[\begin{array}{cccc}\frac{Cut(A_1,\overline{A_1})}{|A_1|}& 0& \cdots & 0\\ 0& \frac{Cut(A_2,\overline{A_2})}{|A_2|}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{Cut(A_k,\overline{A_k})}{|A_k|}\end{array}\right]}_{k\times k}\\ =& \mathbb{T}r{\left[\begin{array}{cccc}{y}_1^{T}L{y}_1& 0& \cdots & 0\\ 0& y_2^{T}L{y}_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & y_k^{T}L{y}_k\end{array}\right]}_{k\times k}\end{array}\text{\hspace{1em}(5)}$$
$\mathbb{T}r$表示求矩阵的迹（也就是对角线求和）。里面是一个对角矩阵。

又因为（对这里不熟悉的，读者可自行尝试，在这里就不展开了）
$$\begin{array}{rl}{Y}^{T}LY=& {\left[\begin{array}{cccc}{y}_1^{T}L{y}_1& y_1^{T}L{y}_2& \cdots & y_1^{T}L{y}_k\\ y_2^{T}L{y}_1& y_2^{T}L{y}_2& \cdots & y_2^{T}L{y}_k\\ ⋮& ⋮& \ddots & ⋮\\ y_k^{T}L{y}_1& y_k^{T}L{y}_2& \cdots & y_k^{T}L{y}_k\end{array}\right]}_{k\times k}\end{array}\text{\hspace{1em}(6)}$$

所以，很容易看到，式（5）和式（6）的对角线是一致的，$\boxed{那么我们就得到目标函数:}$
$$\begin{array}{r}\arg \underset{Y}{\min }\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}=\arg \underset{Y}{\min }\mathbb{T}r(Y^{T}LY)\end{array}$$

3.3、求解

有了目标表达式，接下来就是求解，但是还有个问题，就是矩阵$Y$的每一行是一个指示向量，由这个可得
$$Y^{T}Y=\mathbb{I}$$
$\mathbb{I}$是单位矩阵，比如我以开头的图构造矩阵$Y$，我们切的那刀分成了两份，可得
$$Y^{T}Y=\left[\begin{array}{ccccc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0& 0& 0\\ 0& 0& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\end{array}\right]\left[\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ \frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{3}}\\ 0& \frac{1}{\sqrt{3}}\\ 0& \frac{1}{\sqrt{3}}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$
并且，因为指示向量是离散的。我们是没办法去求导，一般情况下，我们会忽略掉这个指示向量这个约束，而保留$Y^{T}Y=\mathbb{I}$

$\boxed{所以我们的问题就变成了}$
$$\begin{array}{rl}\arg \underset{Y}{\min }\mathbb{T}r& (Y^{T}LY)\\ s.t.& Y^{T}Y=\mathbb{I}\end{array}$$

3.3.1、拉格朗日乘数法

由于$Y$难以求解，所以我们对$Y$的每一列进行求导求解，又因为$\mathbb{T}r(Y^{T}LY)=\sum _i^k{y}_i^{T}L{y}_i$，所以只要每一个最小，那么最后求和也就是最小。但是这种情况其实仅是一个近似解罢了，实际上并不一样，只是由于$Y$太难求了，所以退而求其次。

所以最终的目标函数为
$$\begin{array}{rl}\arg \underset{y_i}{\min }& (y_i^{T}L{y}_i)\\ s.t.& y_i^T{y}_i=1\end{array}$$
约束条件由$Y^{T}Y=\mathbb{I}\to y_i^T{y}_i=1$得到

构造拉格朗日函数
$$L(y_i,\lambda )=y_i^{T}L{y}_i+\lambda (1-y_i^T{y}_i)$$
在这里给出两个求导公式，不懂的可直接背，对其原理感兴趣自行百度查阅
$$\begin{gathered} \boxed{\frac{d{x}^{T}Ax}{x}=2Ax，其中A为对称阵，x为列向量} \\ \boxed{\frac{d{x}^{T}x}{x}=2x，x为列向量} \end{gathered}$$
对$y_i$求导
$$\begin{array}{rl}\frac{\partial L(y_i,\lambda )}{\partial y_i}=& 2L{y}_i-2\lambda y_i\\ =& L{y}_i-\lambda y_i\end{array}$$
令其为0得
$$L{y}_i=\lambda y_i$$

如果你对特征值和特征向量敏感，便一眼可以看出来，此时的$\lambda$就是特征值，而$y_i$就是特征向量。

式子左右左乘$y_i^T$得
$$y_i^{T}L{y}_i=\lambda y_i^T{y}_i=\lambda$$
左边不就是我们的目标函数吗？要最小化它，不就是最小化$\lambda$？而$\lambda$又是拉普拉斯矩阵的特征值。

所以，我们是否可以这么想，对于$Y$的每一列，我们都用拉普拉斯矩阵的特征向量组成，而特征向量的组成便是由前k个特征值（升序）所对应的特征向量。

3.3.2、重离散化

前面我们提到，我们是忽略了指示向量那个离散约束，使得y的取值为全体实数。所以最后我们得出来的$Y$值它也是在全体实数当中的，难以判定它属于哪一类。

所以，我们可以直接在矩阵$Y$当中按行作为样本，使用传统的聚类算法（比如k-mean），实现最终的聚类。

3.4、NCut

弄懂了RatioCut，NCut其实也懂了一半了。

同样的，我们定义
$$vol(A_i)=\sum _{i\in A_i}{d}_i$$
也就是对于类别为$A_i$的集合，集合内所有的节点的度的总和。

而我们的代价函数则定义为
$$NCut(A_1,A_2,\cdots ,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{vol(A_i)}=\sum _{i=1}^k\frac{Cut(A_i,\overline{A_i})}{vol(A_i)}$$
对于我们每一个样本的指示向量，定义为
y = { 1 v o l ( A i ) j ∈ A i 0 j ∉ A i j ∈ { 1 , 2 , ⋯   , k } y=\left\{\begin{matrix} \frac{1}{\sqrt{vol(A_i)}} & j\in A_i \\ 0 & j\notin A_i &j\in\{1,2,\cdots,k\} \end{matrix}\right. y={vol(Ai​) ​1​0​j∈Ai​j∈/Ai​​j∈{1,2,⋯,k}​
同样有（推导过程与RatioCut同，不再赘述）
$$y_i^{T}L{y}_i=\frac{Cut(A_i,\overline{A_i})}{vol(A_i)}$$
所以，同样的都有
$$NCut(A_1,A_2,\cdots ,A_k)=\mathbb{T}r(Y^{T}LY)$$
但是，此时$Y^{T}Y\ne \mathbb{I}$，而是
$$Y^{T}DY=\mathbb{I}$$
证明$\boxed{y_i仍然代表Y的第i列}$：
$$\begin{array}{rl}{y}_i^{T}D{y}_i=& \sum _{j=1}^n{y}_{ij}^2{d}_j\\ =& \sum _{j\in A_i}\frac{1}{vol(A_i)}{d}_j\\ =& \frac{1}{vol(A_i)}\sum _{j\in A_i}{d}_i\\ =& \frac{1}{vol(A_i)}vol(A_i)\\ =& 1\end{array}\text{\hspace{1em}(1)(2)(3)(4)(5)}$$
式子（1）将矩阵展开可得，式子（2）利用指示向量，式子（4）利用了前面定义的$vol(A_i)$。

因为$y_i^{T}D{y}_i$和指示向量的定义，可得$Y^{T}DY=\mathbb{I}$。

所以，目标函数就变成了
$$\begin{array}{rl}& \arg \underset{Y}{\min }\mathbb{T}r(Y^{T}LY)\\ & s.t.Y^{T}DY=\mathbb{I}\end{array}$$
由于$y_i^T{y}_i\ne 1$，所以我们没办法像RatioCut那样去解决。所以我们不妨对这个式子进行一下转化。

令$Y=D^{-\frac{1}{2}}F$，所以
$$\begin{array}{rl}\arg \underset{F}{\min }(Y^{T}LY)=& \arg \underset{F}{\min }\mathbb{T}r(F^T{D}^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}F)\\ & s.t.F^{T}F=\mathbb{I}\end{array}$$
所以，我们只需要将$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$当作原来的$L$即可，那么这样一切便都是照常即可。

同样我们一样一列一列地看（假设某一列为f）
$$\begin{array}{rl}& \arg \underset{f}{\min }(f^T{D}^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}f)\\ & s.t.f^{T}f=1\end{array}$$
构造拉格朗日函数
$$L(f,\lambda )=f^T{D}^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}f+\lambda (1-f^{T}f)$$
对f求导
$$\begin{array}{r}\frac{\partial L(f,\lambda )}{\partial f}=2(D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}})f-2\lambda f\end{array}$$
令其等于0，得
$$(D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}})f=\lambda f$$
所以，一样的，我们只需要求出$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$的前k个最小值。然后组成矩阵F。不同的是，$\boxed{NCut需要将得到的F按行按行归一化}$，再用传统聚类方法聚类即可。

3.4.1、标准化（一点小补充）

其实，$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$也被叫做对拉普拉斯矩阵L的归一化。同样的还有其他的归一化方法，感兴趣可自行查阅。

4、流程

以RatioCut为例：假设保留前$k_1$个特征，聚成$k_2$类

①求出邻接矩阵w，度矩阵D

②计算拉普拉斯矩阵L

③求出L的特征值，并取出前$k_1$个特征值所对应的特征向量组成Y。

④以Y作为样本，使用传统的聚类方法进行聚类，比如K_mean聚类，聚成$k_2$类

5、结束

以上，便是谱聚类的原理的简单推导了，由于笔者能力有限，还有很多的细节都并未深入解读。如有问题，还望指出，阿里嘎多。