谱聚类：原理 + python 实现

这篇文章讲一下谱聚类的推导、实现以及思考（公式太多，可能加载不出来，多刷新几次就行了🥱）

目录展开

1. 简介

谱聚类是一种来自于图论的聚类方法，可区分线性不可分的数据。本质上就是最小化图的割，实际上是将原始样本转换到谱空间再聚类。

2. 问题描述

我们有一张图，希望把图上的点分成 $k$ 类，怎么分呢？谱聚类的想法就是，在图上割 $k$ 刀，分成 $k$ 类，如果割下的边代价最小，那么我们得到的分割就是最好的。

接下来我们就将这个过程转化为数学语言：

定义一个图 $G = {V,E}$，其中 $V = \{1,2,···,N\}$ 为 $N$ 个顶点；

定义矩阵 $W$ 为相似度矩阵（Affinity matrix），其中每一项 $w_{ij}$ 表示样本 $i$ 和 $j$ 的相似度，在图结构中，即为 $i$ 和 $j$ 之间边的权重。 $w_{ij}$ 可以有两种定义方式：

$$w_{ij}=w_{ji}= \begin{cases} 0& {x_i \notin N_j \;and \;x_j \notin N_i}\\ exp(-\dfrac{\|x_i-x_j\|_2^2}{2\sigma^2})& {x_i \in N_j\; or \; x_j \in N_i} \end{cases} $$

如果样本 $i$ 是 $j$ 附近 $k$ 个点的集合 $N_j$，同时 $j$ 也在 $i$ 附近，那么他们的 Affinity 就定义为距离归一再取负指数幂，否则相似度就是 0。

$$w_{ij}=w_{ji}= \begin{cases} 0& {x_i \notin N_j \;or \;x_j \notin N_i}\\ exp(-\dfrac{\|x_i-x_j\|_2^2}{2\sigma^2})& {x_i \in N_j\; and\; x_j \in N_i} \end{cases} $$

另外一种规定，只要一点在另一点附近，即相似。

两种表示方法都可以保证 $W$ 为对称矩阵。

我们再选择图中的两类 $A$ 和 $B$，$A \subset V, B \subset V, A \cap B = \emptyset$，那么 $A$ 和 $B$ 两类之间的相似度

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

下面就进入正题，来定义整张图切割的损失：

$$Ncut(V) = \sum_{i=1}^k\frac{w(A_k,\overline{A_k})}{D(A_k)}$$

其中， $w(A_k,\overline{A_k})$ 表示 $k$ 类与它的补集的相似度，$D(A_k)$ 表示 $k$ 类中所有点两两间边权之和：

$$D(A_k) = \sum_{i\in A_k}d_i =\sum_{i,j\in A_k} w_{ij}$$

这里补充一句，$d_i$ 是点 $i$ 的度，表示和它相连的所有边权之和，用它还能定义一个度矩阵，即：

$$D = \left[ \begin{array}{ccc} d_1 & & &\\ &d_2 & &\\ & &\ddots &\\ & & &d_k\\ \end{array} \right] $$

为什么 $Ncut$ 是这个公式？分子是尽量使切割的每一刀，都切掉图中 $w_{ij}$ 最小的边；分母是使切出的几个部分大小尽量平均，对过大规模的类进行惩罚。

在这张图中，由于加了分母的惩罚，得以选择 Best Cut，而非 Smallest Cut

为了表示“切”这个过程，我们引入一组指示向量 $\{y_i\}_{i=i}^n$，满足

$$y_i\quad s.t.\quad y_i\in \{0,1\}^k, \sum_{j=1}^ky_{ij} = 1$$

即 $y_i$ 为一个长度为 $k$ 的向量，其中只有在 $j$ 号位置是 1，其余位置全是零，$y_{ij} = 1$ 表示样本 $i$ 属于 $j$ 类

那么 $Y = (y_1,y_2,···,y_n)^T_{n\times k}$ 就表示了图分割成 $k$ 类的一种方式

最终我们就把切图转化成一个优化问题：

$$\hat{Y} = \mathop{\arg\min}\limits_{Y} Ncut(V)$$

3. 优化问题求解

$$\begin{aligned} Ncut(V) &= \sum_{i=1}^k\frac{W(A_i,\overline{A_i})}{\sum_{i\in A_i} d_i} \\ &= tr \left[ \begin{smallmatrix} \dfrac{W(A_1,\overline{A_1})}{\sum_{i\in A_1} d_i} & & &\\ & \dfrac{W(A_2,\overline{A_2})}{\sum_{i\in A_2} d_i}& & \\ & & \ddots & \\ & & &\dfrac{W(A_k,\overline{A_k})}{\sum_{i\in A_k} d_i}\\ \end{smallmatrix} \right] \\ &= tr \left[ \begin{smallmatrix} W(A_1,\overline{A_1})& & & \\ & W(A_2,\overline{A_2})& & \\ & & \ddots & \\ & & & W(A_k,\overline{A_k})\\ \end{smallmatrix} \right] \cdot \left[ \begin{smallmatrix} \sum_{i\in A_1} d_i& & & \\ & \sum_{i\in A_2} d_i& & \\ & & \ddots & \\ & & & \sum_{i\in A_k} d_i\\ \end{smallmatrix} \right]^{-1} \\ &= tr(O_{k\times k}\cdot P_{k\times k}^{-1})\\\end{aligned}$$

待优化问题便转换为：

$$\mathop{\arg\min}\limits_{Y} tr(O\cdot P^{-1})$$

而我们已知的只有 $Y,W,D$，因此接下来尝试用它们表示矩阵 $O,P$

3.1 矩阵 P 的表示

$$ \begin{aligned} Y\cdot Y^T_{k \times k} = \left( \begin{array}{ccc} y_1& y_2 & \cdots & y_n\\ \end{array} \right) \cdot \left( \begin{array}{ccc} y_1^T\\ \vdots\\ y_n^T\\ \end{array} \right) = \sum_{i=1}^ny_i\cdot y_i^T \end{aligned} $$

我们可以看到，当样本 $i$ 属于第 $j$ 类时，由于 $y_i$ 为长度为 $k$ 的向量，因此 $y_i\cdot y_i^T$ 的结果是 $k\times k$ 的矩阵；由于 $y_i$ 只在 $j$ 位置是 1，其余位置全零，因此 $y_i\cdot y_i^T$ 的结果在 $(j,j)$ 处为 1 ，其余所有位置全零。

那么便可以得到： $$ Y\cdot Y^T = \sum_{i=1}^ny_i\cdot y_i^T = \left[ \begin{smallmatrix}N_1 & & &\\\\ &N_2 & &\\\\ & &\ddots &\\\\ & & &N_k\\\\ \end{smallmatrix} \right] = \left[ \begin{smallmatrix}\sum_{i \in A_1}1 & & &\\ &\sum_{i \in A_2}1 & &\\ & &\ddots &\\ & & &\sum_{i \in A_k}1\\ \end{smallmatrix} \right] $$

即 $Y\cdot Y^T$ 结果也是一个 $k\times k$ 的矩阵，在 $(j,j)$ 处为第 $j$ 类包含的元素数量，其余所有位置全零。

此时发现，$Y\cdot Y^T$ 与我们定义的矩阵 $p$ 是那么地相似！只是 $\sum 1$ 与 $\sum d_i$ 的区别！显然：

$$P = Y^T D Y $$

其中 $D$ 是我们定义的度矩阵

3.2 矩阵 O 的表示

$$W(A_k,\overline{A_k}) = W(A_k,V) – W(A_k,A_k) = \sum_{i\in A_k}d_i – \sum_{i,j\in A_k}w_{ij}$$ $$\begin{aligned} O =& \left[ \begin{smallmatrix} W(A_1,\overline{A_1})& & & \\ & W(A_2,\overline{A_2})& & \\ & & \ddots & \\ & & & W(A_k,\overline{A_k})\\ \end{smallmatrix} \right]\\ =& \left[ \begin{smallmatrix} \sum_{i\in A_1} d_i& & & \\ & \sum_{i\in A_2} d_i& & \\ & & \ddots & \\ & & & \sum_{i\in A_k} d_i\\ \end{smallmatrix} \right] – \left[ \begin{smallmatrix} W(A_1,A_1)& & & \\ & W(A_2,A_2)& & \\ & & \ddots & \\ & & & W(A_k,A_k)\\ \end{smallmatrix} \right] \\ =& Y^T D Y – Q \end{aligned} $$

现在只需表示出矩阵 $Q$ 即可。我们可以先算一下 $Y^T W Y$：

$$\begin{aligned} Y^T_{k\times n} W_{n \times n} Y_{n \times k} &= \left( \begin{array}{ccc} y_1& y_2 & \cdots & y_n\\ \end{array} \right) \cdot \left( \begin{array}{ccc} W_{11}&\cdots&W_{n1}\\ \vdots&\ddots&\vdots\\ W_{1n}&\cdots&W_{nn}\\ \end{array} \right) \cdot \left( \begin{array}{ccc} y_1^T\\ \vdots\\ y_n^T\\ \end{array} \right)\\ &= \left( \begin{array}{ccc} \sum_{i=1}^n y_iw_{i1}&\cdots & \sum_{i=1}^n y_iw_{in}\\ \end{array} \right) \cdot \left( \begin{array}{ccc} y_1^T\\ \vdots\\ y_n^T\\ \end{array} \right)\\ &= \sum_{i=1}^n \sum_{j=1}^n y_iy_j^Tw_{ij} \end{aligned}$$

这是我们又看到了熟悉的 $y_iy_j^T$，根据刚才的理解，可以得到

$$\begin{aligned} Y^T_{k\times n} W_{n \times n} Y_{n \times k} &= \left( \begin{array}{ccc} \sum_{i\in A_1}\sum_{j\in A_1}w_{ij}&\cdots&\sum_{i\in A_1}\sum_{j\in A_k}w_{ij}\\ \vdots&\ddots&\vdots\\ \sum_{i\in A_k}\sum_{j\in A_1}w_{ij}&\cdots&\sum_{i\in A_k}\sum_{j\in A_k}w_{ij}\\ \end{array} \right) \end{aligned}$$

此时我们又发现一个巧合：要表示的矩阵 $Q$ 即为 $Y^T W Y$ 对角线上的元素！由于要求的优化问题是几个对角矩阵相乘再求迹，那么如果用 $Y^T W Y$ 代替 $Q$，也不会影响最后的结果！因此，我们也得到了矩阵 $O$ 的表示。

3.3 优化问题的解

$$\hat{Y} = \mathop{\arg\min}\limits_{Y} tr(O\cdot P^{-1}) = \mathop{\arg\min}\limits_{Y} tr(\frac{Y^T(D-W)Y}{Y^TDY})$$

这里的 $D-W$ 就是图拉普拉斯矩阵 $L$，对于图拉普拉斯矩阵，有如下性质：

由于 $D,W$ 均为对称矩阵，因此 $L$ 也是对称矩阵，其特征值为实数
对于任意的向量 $f$，有：$$\begin{aligned} f^TLf &= f^TDf – f^TWf = \sum_{i=1}^nd_if_i^2 – \sum_{i=1}^n\sum_{j=1}^nw_{ij}f_if_j\\ &= \frac{1}{2}(\sum d_if_i^2 – 2 \sum_{i=1}^n\sum_{j=1}^nw_{ij}f_if_j + \sum d_jf_j^2)\\ &= \frac{1}{2}\sum_{i=1}^n\sum_{j=1}^nw_{ij}(f_i-f_j)^2 \ge 0 \end{aligned}$$因此 $L$ 半正定，其特征值均 $\ge$ 0

回到优化问题哈：

$$\hat{Y} =\mathop{\arg\min}\limits_{Y} tr(\frac{Y^TLY}{Y^TDY})$$

令 $Y = D^{-\frac{1}{2}}Z$，那么 $Y^TY = Z^TZ, Y^TLY = Z^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}Z$，且 $Z^TZ = I$

优化问题就变成

$$\mathop{\arg\min}\limits_{Z} tr(\frac{Z^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}Z}{Z^TZ})$$

我们要求的是 $trace$，那么可以拆成 $k$ 个子问题，即让矩阵对角线上每个值都最小，即

$$\mathop{\arg\min}\limits_{Z} \frac{z_i^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}z_i}{z_i^Tz_i}$$

其中 $z_i^Tz_i = \textbf{1}$，所以分母去掉，变成带约束的优化问题：

$$\mathop{\arg\min}\limits_{Z} z_i^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}z_i,\qquad s.t.\quad z_i^Tz_i = \textbf{1}$$

利用拉格朗日乘数法，则：

$$\begin{aligned} L(z_i) &= z_i^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}z_i – \lambda(z_i^Tz_i – 1) \\ \frac{\partial L(z_i)}{\partial z_i} &= D^{-\frac{1}{2}}LD^{-\frac{1}{2}}z_i – \lambda z_i = 0 \\ \Rightarrow \lambda \cdot z_i &= D^{-\frac{1}{2}}LD^{-\frac{1}{2}}z_i \end{aligned}$$

这不就是对 $D^{-\frac{1}{2}}LD^{-\frac{1}{2}}$ 进行特征值分解嘛！所以最小的 $z_i$ 就是最小特征值对应的特征向量，最小的 $z_i$ 对应的 $y_i$ 即为把图切一刀的结果。

现在终于理解，为什么谱聚类是选最小的特征值了吧！

但是这里还有一个问题：我们一开始设计的指示向量是 $y_i\in \{0,1\}^k$，是一个冲激信号，不现实。我们特征分解后得到的特征向量，也一定是连续的：每取一个最小的特征值，就意味着有一些误差；切太多刀，误差就会积累。那怎么办呢？

4. 对谱方法的理解

先来理解一下“最小特征值对应的特征向量”，她代表什么呢？

上图的方向是最大特征值对应特征向量所代表的方向。特征值最大，就相当于在这个特征向量代表的线性子空间，是拉得最开的。也就四说，样本点在这个方向上最分散，即方差最大，这就是 PCA 降维的目的。

上图的方向是最小特征值对应特征向量所代表的方向。各点到这个方向距离的差距将最大。我认为，在这个方向切一刀，可以把这张图最好地切开：在这个维度上，坐标为正的属于一类，负属于另一类。每切一刀，就代表在特征空间里，选出一个维度，在这个维度某点处左右分开。选择 $k$ 个，就相当于将原始数据映射到 $k$ 维特征空间。因此，谱聚类实际上是将原始样本转换到谱空间，再聚类。

5. 实现

谱聚类的流程及代码各种博客都有写，这里不再赘述。（建议参考此文的实现）

与 k-means 对比，的确可以分开线性不可分的数据

6. 进一步思考

上文其实遗漏了一个问题，应该取几个最小的特征值呢？

请教了导师，导师给出这样的方法：

画出拉普拉斯矩阵的所有特征值（左图），观察最小的几个特征值（右图），发现前两个特征值非常小，接近于 0，紧接着就有个跳变，因此选择两个特征值。为啥这样做捏？选多了当然不行，因为上文说过，会误差积累。一维又肯定不行——如果选一维就够，那数据就是线性可分的。但是“选跳变”这个方法的根本原因，我还是没能琢磨明白。

后来还发现一个规律：特征值跳变的越明显，那么谱聚类效果越好。