谱聚类:原理 + 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$,对于图拉普拉斯矩阵,有如下性质:

  1. 由于 $D,W$ 均为对称矩阵,因此 $L$ 也是对称矩阵,其特征值为实数
  2. 对于任意的向量 $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,紧接着就有个跳变,因此选择两个特征值。为啥这样做捏?选多了当然不行,因为上文说过,会误差积累。一维又肯定不行——如果选一维就够,那数据就是线性可分的。但是“选跳变”这个方法的根本原因,我还是没能琢磨明白。

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

转换到二维空间中,样本点是上图的分布

将相似度矩阵按照类别排列,可得到右图,是优雅的分块对角化形式。这就表明,类内挺相似,类类之间不相似,整挺好的。

我又用同样的半月形数据重复试验,还发现一个问题:在数据噪声不同的情况下,最好的特征空间维度也可能不同!比如在噪声较小时,二维特征空间是最好的;而噪声较大时,四维特征空间是最好的。

画出他们的最小特征值,发现果然,特征值跳变不明显,聚类效果也不咋地(我为了统计性能,每个 noise 都做了 50 遍,就不上图了)

特征维数随着噪音变化。这是为什么呢?我认为,特征空间转换与拟合函数类似,是一个不适定的问题,其解是敏感的:当输入有微小的变化,输出可能翻天覆地。

[mathjax]

Leave a Comment

电子邮件地址不会被公开。