专栏名称: 酷酷的群

个人公众号：酷酷的群

谱聚类|机器学习推导系列（二十）

酷酷的群 · 简书 · · 2021-01-31 13:40

正文

一、概述

对于下图所示的数据进行聚类，可以采用GMM或者K-Means的方法：

数据

然而对于下图所示的数据，单纯的GMM和K-Means就无效了，可以通过核方法对数据进行转换，然后再进行聚类：

数据

对于上图所示的数据进行聚类可以考虑采用谱聚类（spectral clustering）的方法。

总结来说，聚类算法可以分为两种思路：

①Compactness，这类有 K-means，GMM 等，但是这类算法只能处理凸集，为了处理非凸的样本集，必须引⼊核技巧。
②Connectivity，这类以谱聚类为代表。

二、基础知识

无向权重图

谱聚类的方法基于带权重的无向图，图的每个节点是一个样本点，图的边有权重，权重代表两个样本点的相似度。构建图时可以选择将每个样本点与其相似度最高的 $k$ 个点连接起来。

假设总共 $N$ 个样本点，这些样本点构成的图可以用 $G=(V,E)$ 表示，其中 $V=\left \{v_1,v_2,\cdots ,v_N\right \}$ ，图中的每个点 $v_i$ 也就代表了一个样本 $x_i$ ， $E$ 是边，用邻接矩阵（也是相似度矩阵） $W_{N\times N}$ 来表示， $W=[w_{ij}],1\leq i,j\leq N$ ，由于是无向图，因此 $w_{ij}=w_{ji}$ 。

另外还有度的概念，这里可以类比有向图中的出度和入度的概念，不过图中的点 $v_i$ 的度 $d_i$ 并不是和该点相连的点的数量，而是和其相连的边的权重之和，也就是邻接矩阵的每一行的值加起来，即：

$d_{i}=\sum_{j=1}^{N}w_{ij}$

而图的度矩阵（对角矩阵） $D_{N\times N}$ 可以表示如下：

$D=\begin{bmatrix} d_{1} & & & \\ & d_{2} & & \\ & & & \\ & & & d_{N} \end{bmatrix}$

另外我们定义，对于点集 $V$ 的一个子集 $A\subset V$ ，我们定义：

$|A|:=子集A中点的个数\\ vol(A):=\sum _{i\in A}d_{i}$

邻接矩阵

构建邻接矩阵 $W$ 一共有三种方法，分别是 $\epsilon$ -近邻法、 $k$ 近邻法和全连接法。

$\epsilon$ -近邻法

首先需要设置一个阈值 $\epsilon$ ，比较任意两点 $x_i$ 与 $x_j$ 之间的距离 $s_{ij}=||x_{i}-x_{j}||_{2}^{2}$ 与 $\epsilon$ 的大小，定义邻接矩阵如下：

$w_{ij}=\left\{\begin{matrix} 0,s_{ij}>\epsilon \\ \epsilon ,s_{ij}\leq \epsilon \end{matrix}\right.$

这种方法表示如果两个样本点之间的欧氏距离的平方小于阈值 $\epsilon$ ，则它们之间是有边的。

使用这种方法，两点相似度只有 $\epsilon$ 和 $0$ 两个值，这种度量很不精确，因此在实际应用中很少使用 $\epsilon$ -近邻法。

$k$ 近邻法

使用KNN算法遍历所有样本点，取每个样本点最近的 $k$ 个点作为近邻。这种方法会造成构造的邻接矩阵不对称，而谱聚类算法需要一个对称的邻接矩阵。因此有以下两种方法来构造一个对称的邻接矩阵：

①只要一个点在另一个点的 $k$ 近邻内，则 $w_{ij}>0$ ，否则为 $0$ ，相似度 $w_{ij}$ 可以使用径向基函数来度量：

$w_{ij}=w_{ji}=\left\{\begin{matrix} exp\left \{-\frac{||x_{i}-x_{j}||_{2}^{2}}{2\sigma ^{2}}\right \},x_{i}\in KNN(x_{j})\; or\; x_{j}\in KNN(x_{i})\\ 0,x_{i}\notin KNN(x_{j})\; and\; x_{j}\notin KNN(x_{i}) \end{matrix}\right.$

②只有两个点互为 $k$ 近邻，才会有 $w_{ij}>0$ ，否则为 $0$ ：

$w_{ij}=w_{ji}=\left\{\begin{matrix} exp\left \{-\frac{||x_{i}-x_{j}||_{2}^{2}}{2\sigma ^{2}}\right \},x_{i}\in KNN(x_{j})\; and\; x_{j}\in KNN(x_{i})\\ 0,x_{i}\notin KNN(x_{j})\; or\; x_{j}\notin KNN(x_{i}) \end{matrix}\right.$

上述方法是不用先建立图而直接获得邻接矩阵，在编程实现时能够更加简便，构建的邻接矩阵也就表明了哪些样本点之间有边连接。也可以采用先建立图然后再在图上有边的数据点上保留权重获得邻接矩阵的方法。

全连接法

这种方法会使所有的 $w_{ij}$ 都大于 $0$ ，可以选择不用的核函数来度量相似度，比如多项式核函数、径向基核函数和 $sigmoid$ 核函数。最常用的是径向基核函数：

$w_{ij}=exp\left \{-\frac{||x_{i}-x_{j}||_{2}^{2}}{2\sigma ^{2}}\right \}$

在实际应用时选择全连接法建立邻接矩阵是最普遍的，在选择相似度度量时径向基核函数是最普遍的。

拉普拉斯矩阵

图的拉普拉斯矩阵（Graph Laplacian） $L_{N\times N}$ 是一个对称矩阵，用度矩阵减去邻接矩阵得到的矩阵就被定义为拉普拉斯矩阵， $L=D-W$ 。拉普拉斯矩阵有一些性质如下：
①对称性。
②由于其对称性，则它的所有特征值都是实数。
③对于任意向量 $f$ ，有：

$f^{T}Lf=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}(f_{i}-f_{j})^{2}$

这一性质利用拉普拉斯矩阵的性质很容易可以得到：

$f^{T}Lf=f^{T}Df-f^{T}Wf \\ =\sum _{i=1}^{N}d_{i}f_{i}^{2}-\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}f_{i}f_{j}\\ =\frac{1}{2}(\sum _{i=1}^{N}d_{i}f_{i}^{2}-2\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}f_{i}f_{j}+\sum _{j=1}^{N}d_{j}f_{j}^{2})\\ =\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}w_{ij}f_{i}f_{j}+\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}f_{j}^{2})\\ =\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}(f_{i}-f_{j})^{2}$

④拉普拉斯矩阵是半正定的，则其所有特征值非负，这个性质由性质③很容易得出。并且其最小的特征值为 $0$ ，这是因为 $L$ 的每一行和为 $0$ ，对于全 $1$ 向量 $1_{N}=\begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}^{T}$ ，有 $L\cdot 1_{N}=0=0\cdot 1_{N}$ 。

无向图切图

对于无向图 $G$ 的切图，我们的目标是将 $G=(V,E)$ 切成相互没有连接的 $k$ 个子图，每个子图节点的集合为 $A_{1},A_{2},\cdots ,A_{k}$ ，它们满足 $A_{i}\cap A_{j}=\phi$ ，且 $A_{1}\cup A_{2}\cup \cdots \cup A_{k}=V$ 。

对于任意两个子图点的集合 $A,B\subset V$ ，定义 $A$ 和 $B$ 之间的切图权重为：

$W(A,B)=\sum _{i\in A,j\in B}w_{ij}$

对于 $k$ 个子图的集合，定义切图 $cut$ 为：

$cut(A_{1},A_{2},\cdots ,A_{k})=\sum_{i=1}^{k}W(A_{i},\bar{A}_{ i })$

上式中 $\bar{A}_{ i }$ 是 $A_{i}$ 的补集，意为除 $A_{i}$ 子集以外的其他子集的并集。

每个子图就相当于聚类的一个类，找到子图内点的权重之和最高，子图间的点的权重之和最低的切图就相当于找到了最佳的聚类。实现这一点的一个很自然的想法是最小化 $cut$ 。然而这种方法存在问题，也就是最小化的 $cut$ 对应的切图不一定就是符合要求的最优的切图，如下图：

举例

在上面的例子中，我们选择在最小的权重上进行切图，比如在 $C$ 和 $H$ 之间进行切图，这样可以使得 $cut$ 最小，但并不是最优的切图。

接下介绍谱聚类使用的切图方法。

三、谱聚类之切图聚类

为了避免上述最小化 $cut$ 存在的问题，需要对每个子图的规模做出限制，接下来介绍两种切图方法，分别是 $RatioCut$ 和 $NCut$ 。

$RatioCut$ 切图

$RatioCut$ 切图为了避免上述最小切图，对于每个切图，不只考虑最小化 $cut$ ，还考虑最大化每个子图点的个数，即：

$RatioCut(A_{1},A_{2},\cdots ,A_{k})=\sum_{i=1}^{k}\frac{W(A_{i},\bar{ A }_{i})}{|A_{i}|}$

为了最小化 $RatioCut$ 这个函数，我们引入指示向量 $h_{i}\in \left \{h_{1},h_{1},\cdots ,h_{k}\right \},i=1,2,\cdots ,k$ ，对于每一个向量 $h_{i}$ ，它是一个 $N$ 维向量，另外定义 $h_{ij}$ 为：

$h_{ij}=\left\{\begin{matrix} 0,v_{j}\notin A_{i}\\ \frac{1}{\sqrt{|A_{i}|}},v_{j}\in A_{i} \end{matrix}\right.$

那么对于 $h_{i}^{T}Lh_{i}$ ，有：

$h_{i}^{T}Lh_{i}=\frac{1}{2}\sum_{m=1}^{N}\sum_{n=1}^{N}w_{mn}(h_{im}-h_{in})^{2}\\ =\frac{1}{2}(\sum _{m\in A_{i},n\notin A_{i}}w_{mn}(\frac{1}{\sqrt{|A_{i}|}}-0)^{2}+\sum _{m\notin A_{i},n\in A_{i}}w_{mn}(0-\frac{1}{\sqrt{|A_{i}|}})^{2})\\ =\frac{1}{2}(\sum _{m\in A_{i},n\notin A_{i}}w_{mn}\frac{1}{|A_{i}|}+\sum _{m\notin A_{i},n\in A_{i}}w_{mn}\frac{1}{|A_{i}|})\\ =\frac{1}{2}(\frac{1}{|A_{i}|}\sum _{m\in A_{i},n\notin A_{i}}w_{mn}+\frac{1}{|A_{i}|}\sum _{m\notin A_{i},n\in A_{i}}w_{mn})\\ =\frac{1}{2}(\frac{cut(A_{i}, \bar{A}_{i} )}{|A_{i}|}+\frac{cut( \bar{A}_{i} ,A_{i})}{|A_{i}|})\\ =\frac{cut(A_{i},\bar{A}_{i})}{|A_{i}|}$

由上式可知，某一个子图的 $RatioCut$ 也就是 $h_{i}^{T}Lh_{i}$ ，所有的 $k$ 个子图的 $RatioCut$ 表达式也就是：

$RatioCut(A_{1},A_{2},\cdots ,A_{k})=\sum_{i=1}^{k}h_{i}^{T}Lh_{i} \\ =\sum_{i=1}^{k}(H^{T}LH)_{ii}\\ =tr(H^{T}LH)$

上式中 $tr(H^{T}LH)$ 为矩阵 $H^{T}LH$ 的迹， $H=\begin{pmatrix} h_{1} & h_{2} & \cdots & h_{k} \end{pmatrix}$ ，需要注意这里的 $H$ 满足 $H^TH=I$ ，并且 $H$ 的元素只能取 $0$ 或者 $\frac{1}{|A_{i}|}$ 。也就是说我们需要优化以下目标函数：

$\underset{H}{argmin}\; tr(H^{T}LH)\; \; s.t.\; H^{T}H=I$

由于每个元素只能取两个值，因此上面的目标函数是不可求导的。这里每个指示向量都是 $N$ 维的，而且每个元素只有两种取值，所以就有 $2^N$ 种取值方式，一共有 $k$ 个指示向量，因此共有 $k2^N$ 种 $H$ ，因此想要找到满足使目标函数最小的 $H$ 是一个NP难的问题。

由于存在上述问题，所以我们采用降维的思想来考虑解决这个优化问题。我们需要最小化 $tr(H^{T}LH)$ ，也就是需要优化每一个 $h_{i}^{T}Lh_{i}$ ，这里的 $h$ 是单位正交基， $L$ 是对称矩阵，因此 $h_{i}^{T}Lh_{i}$ 的最大值是 $L$ 的最大特征值，最小值是 $L$ 的最小特征值。之所以有这种结论可以参考主成分分析PCA的解法，在PCA中需要找到协方差矩阵（类比此处的拉普拉斯矩阵 $L$ ，它们都是对称的）的最大特征值，而在谱聚类中需要找到最小的 $k$ 个非零特征值，然后得到这些特征值对应的特征向量，通过这个过程我们也就完成了数据的降维，最终 $H_{N\times k}$ 就是降维的结果，使用这个结果来近似解决这个NP难的问题。

一般我们仍然需要对 $H$ 按行做标准化，也就是：

$h_{ij}^{*}=\frac{h_{ij}}{(\sum_{t=1}^{k}h_{it}^{2})^{1/2}}$

由于在降维时损失了少量信息，导致得到的优化后的指示向量 $h$ 对应的 $H$ 现在不能完全指示各样本的归属，因此在得到降维结果 $H$ 后还需要进行一次传统的聚类，比如K-Means。

$NCut$ 切图

$NCut$ 切图的方法与 $RatioCut$ 切图的方法很类似，只是把 $RatioCut$ 的分母 $|A_{i}|$ 换成 $vol(A_{i})$ 。使用 $NCut$ 切图时，由于子图样本个数多不一定权重就大（只有权重大时，子图内样本点的相似度才高），因此切图时基于权重也更符合目标，一般来说 $NCut$ 切图优于 $RatioCut$ 切图：

$NCut(A_{1},A_{2},\cdots ,A_{k})=\sum_{i=1}^{k}\frac{W(A_{i},\bar{ A }_{ i } )}{vol(A_{i})}$

另外需要修改指示向量的表示形式， $RatioCut$ 的指示向量使用 $\frac{1}{\sqrt{|A_{i}|}}$ 来标示样本归属，而 $NCut$ 使用子图权重 $\frac{1}{\sqrt{vol(A_{i})}}$ 来标示指示向量 $h$ ，定义如下：

$h_{ij}=\left\{\begin{matrix} 0,v_{j}\notin A_{i}\\ \frac{1}{\sqrt{vol(A_{i})}},v_{j}\in A_{i} \end{matrix}\right.$

类似的，对于 $h_{i}^{T}Lh_{i}$ ，有：

$h_{i}^{T}Lh_{i}=\frac{1}{2}\sum_{m=1}^{N}\sum_{n=1}^{N}w_{mn}(h_{im}-h_{in})^{2}\\ =\frac{1}{2}(\sum _{m\in A_{i},n\notin A_{i}}w_{mn}(\frac{1}{\sqrt{vol(A_{i})}}-0)^{2}+\sum _{m\notin A_{i},n\in A_{i}}w_{mn}(0-\frac{1}{\sqrt{vol(A_{i})}})^{2})\\ =\frac{1}{2}(\sum _{m\in A_{i},n\notin A_{i}}w_{mn}\frac{1}{vol(A_{i})}+\sum _{m\notin A_{i},n\in A_{i}}w_{mn}\frac{1}{vol(A_{i})})\\ =\frac{1}{2}(\frac{1}{vol(A_{i})}\sum _{m\in A_{i},n\notin A_{i}}w_{mn}+\frac{1}{vol(A_{i})}\sum _{m\notin A_{i},n\in A_{i}}w_{mn})\\ =\frac{1}{2}(\frac{cut(A_{i},\bar{ A }_{ i } )}{vol(A_{i})}+\frac{cut( \bar{ A }_{ i } ,A_{i})}{vol(A_{i})})\\ =\frac{cut(A_{i}, \bar{ A }_{ i } )}{vol(A_{i})}$

同样的优化目标也就是：

$NCut(A_{1},A_{2},\cdots ,A_{k})=\sum_{i=1}^{k}h_{i}^{T}Lh_{i} \\ =\sum_{i=1}^{k}(H^{T}LH)_{ii}\\ =tr(H^{T}LH)$

但是现在的约束条件不再满足 $H^TH=I$ ，而是 $H^TDH=I$ ，证明如下：

$H^{T}DH=\begin{pmatrix} h_{1}^{T}\\ h_{2}^{T}\\ \vdots \\ h_{k}^{T} \end{pmatrix}\begin{pmatrix} d_{1} & & & \\ & d_{2} & & \\ & & \ddots & \\ & & & d_{N} \end{pmatrix}\begin{pmatrix} h_{1} & h_{2} & \cdots & h_{k} \end{pmatrix}\\ =\begin{pmatrix} h_{11}d_{1} & h_{12}d_{2} & \cdots & h_{1N}d_{N}\\ h_{21}d_{1} & h_{22}d_{2} & \cdots & h_{2N}d_{N}\\ \vdots & \vdots & \ddots & \vdots \\ h_{k1}d_{1} & h_{k2}d_{2} & \cdots & h_{kN}d_{N} \end{pmatrix}\begin{pmatrix} h_{1} & h_{2} & \cdots & h_{k} \end{pmatrix}\\ =\begin{pmatrix} \sum_{i=1}^{N}h_{1i}^{2}d_{i} & \sum_{i=1}^{N}h_{1i}h_{2i}d_{i} & \cdots & \sum_{i=1}^{N}h_{1i}h_{ki}d_{i}\\ \sum_{i=1}^{N}h_{2i}h_{1i}d_{i} & \sum_{i=1}^{N}h_{2i}^{2}d_{i} & \cdots & \sum_{i=1}^{N}h_{2i}h_{ki}d_{i}\\ \vdots & \vdots & \ddots & \vdots \\ \sum_{i=1}^{N}h_{ki}h_{1i}d_{i} & \sum_{i=1}^{N}h_{ki}h_{2i}d_{i} & \cdots & \sum_{i=1}^{N}h_{ki}^{2}d_{i} \end{pmatrix}$

对于对角线元素有：

$\sum_{j=1}^{N}h_{ij}^{2}d_{j}=\frac{1}{vol(A_{i})}\sum _{j\in A_{i}}d_{j}=\frac{1}{vol(A_{i})}vol(A_{i})=1$

由于 $h_{mi}$ 和 $h_{ni}$ 不可能同时非零，因此对于非对角线元素有：

$\sum_{i=1}^{N}h_{mi}h_{ni}d_{i}=\sum_{i=1}^{N}0\cdot d_{i}=0$

因此有 $H^TDH=I$ 。我们最终优化的目标函数为：

$\underset{H}{argmin}\; tr(H^{T}LH)\; \; s.t.\; H^{T}DH=I$

此时指示向量 $h$ 并不是标准正交基，所以在 $RatioCut$ 中的降维思想不能直接使用。对于这个问题，只需要将指示向量 $h$ 做一个转化即可，我们令 $H=D^{-1/2}F$ ，则：

$H^{T}LH=F^{T}{\color{Red}{D^{-1/2}LD^{-1/2}}}F\\ H^{T}DH=F^{T}F=I$

也就是说优化的目标变成了：

$\underset{F}{argmin}\; tr(F^{T}{\color{Red}{D^{-1/2}LD^{-1/2}}}F)\; \; s.t.\; F^{T}F=I$

可以发现这个式子和 $RatioCut$ 基本一致，只是中间的 $L$ 变成了 $D^{-1/2}LD^{-1/2}$ 。如此我们就可以按照 $RatioCut$ 的思想，求出 $D^{-1/2}LD^{-1/2}$ 的前 $k$ 个最小非零特征值，然后求对应的特征向量再进行标准化得到最后的特征矩阵 $F$ ，然后再使用K-Means等传统方法进行聚类即可。

一般来说， $D^{-1/2}LD^{-1/2}$ 相当于对拉普拉普斯矩阵做了一次标准化，即 $(D^{-1/2}LD^{-1/2})_{ij}=\frac{L_{ij}}{\sqrt{d_{i}*d_{j}}}$ 。

四、总结

算法流程

以 $NCut$ 切图为例总结一下谱聚类算法的流程：

输入：样本集 $D=(x_{1},x_{2},\cdots ,x_{N})$ ，邻接矩阵的生成方式，降维后的维度 $k_1$ ，聚类方法，聚类的簇个数 $k_2$ 。
输出：簇划分 $C(c_{1},c_{2},\cdots ,c_{k_{2}})$ 。
①根据输入的邻接矩阵生成方式构建样本的邻接矩阵矩阵 $W$ 和度矩阵 $D$ ；
②计算拉普拉斯矩阵 $L$ ；
③构建标准化后的拉普拉斯矩阵 $D^{-1/2}LD^{-1/2}$ ；
④计算 $D^{-1/2}LD^{-1/2}$ 的最小的前 $k_1$ 个非零特征值对应的特征向量 $f$ ；
⑤将各自对应的特征向量 $f$ 组成的矩阵按行标准化，最终得到 $N\times k_1$ 维的特征矩阵 $F$ ；
⑥对 $F$ 的每一行作为一个 $k_1$ 维的样本，共 $N$ 个样本，用输入的聚类方法进行聚类，聚类的簇的个数为 $k_2$ ；
⑦得到簇划分 $C(c_{1},c_{2},\cdots ,c_{k_{2}})$ 。

优缺点

谱聚类的优点有：
①谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到。
②由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。

谱聚类的缺点有：
①如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。
②聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。

参考资料

ref: 谱聚类（spectral clustering）原理总结