谱聚类

无须事先指定簇数量的聚类

说到聚类，最常见的模型当然是 KMeans。不过如果使用 KMeans 的话，需要在算法运行前指定 $k$ 的值——也就是要在训练前指定最后的结果被分为几簇。

现实中有相当多的聚类问题，无法事先指定簇的数量。KMeans 就无法完成这类任务。

好在聚类方法有很多，有一种算法，不仅不需要事先指定 $k$ 值，还可以在结果中保证每个簇中的个体数量低于某个量值，这就是基于图切割的 谱聚类（Spectral Clustering） 。

算法实现

基于图切割的谱聚类算法过程分为两个大的步骤：

1. 图切割
2. 谱聚类

具体步骤如下。

Step 1：生成一张图 $G = <V,E>$，其中每个顶点（Vertex）对应一个样本对象，每两个顶点之间的边则代表这两个样本之间的距离。

此处的距离可以是欧氏距离、余弦距离，或者任何一种距离，我们用 $c_{ij}$ 表示顶点 $i$ 和顶点 $j$ 之间的距离，那么这张图就可以用矩阵 $C$ 来表示了：

$C = (c_{ij}) $

Step 2：确定距离阈值 $threshold_C$，将所有 $c_{ij} > threshold_C$ 的顶点对 ${ i, j }$ 视为断开。

据此将完整的 $G$ 分割为若干连通图 ${G_1, G_2, … , G_n}$。

计算每一个子图的 Radius (最远边缘节点到中心节点的距离) 和 Size（包含顶点数）：

• 如果 $( cluster_{radius}\leqslant threshold_{radius}) $ && $ (cluster_{size} \leqslant threshold_{size})$， 则该连通图本身就是一个独立的簇；

• 否则，对该簇进行下一个步骤，即 Step 3。

Step 3：图切割，主要包括以下两个步骤。

Step 3.1：将待切割图 $G$ 切割为两个子图 $G_{s_1}$ 和 $G_{s_2}$，使得 $G_{s_1}$ 和 $G_{s_2}$ 之间距离尽量大，而两个子图内部节点间距离尽量小。

具体切割过程如下：

Step 3.1.1：构造一个和 $C$ 同等大小的矩阵：

$W = (w_{ij})$

$w_{ij} = \exp{(-\frac{c_{ij}^2}{2\sigma^2})}$

$w_{ij} = \exp{(-\frac{||x_i – x_j||^2}{2\sigma^2})} $

这里用到了高斯相似度函数（Gaussian Similarity Function），其中 $\sigma$ 是一个用来控制变换速度的参数。

其中，$||x_i-x_j||$ 是样本 $i$ 到样本 $j$ 的距离，也就是说 $||x_i-x_j|| = c_{ij}$。

Step 3.1.2：构造一个对角矩阵 $D$：

对角项 $d_i = \sum_{j}w_{ij}$。

Step 3.1.3：令 $L = D - W$，构造等式：

$Lv = (D-W)v = λv$，其中 $λ$ 为常数，$v$ 为向量。

Step 3.1.4：计算上述等式的第二小的特征值所对应的特征向量 $f$。

注意：为什么要取第二小的特征值对应的特征向量，理由见后面描述。

设被分割图 $G$ 一共包含 $n$ 个顶点，则其由一个 $n \times n$ 矩阵表达，因此得出的特征向量 $f$ 也是 $n$ 维向量。$f$ 中的每一维代表一个顶点（即一个样本）：

$f = (f_1, f_2, …, f_n)$;

如果 $ f_i \geqslant 0 $，那么对应的顶点 $i$ 属于 $G_{s_1}$。

如果 $f_i < 0$，则对应的顶点$i$属于 $G_{s_2}$。

这样就把 $G$ 分成了两部分：$G_{s_1}$ 和 $G_{s_2}$。

Step 3.2：计算 $G_{s_1}$ 和 $G_{s_2}$ 的大小和半径，然后：

IF $((cluster_{radius} > threshold_{radius})$ && $(cluster_{size} > threshold_{size})) $

THEN

重复Step 1，直到所有被分割的结果满足上述条件为止。

Step 4：将 Step 3 运用到 Step 2 中所有连通图 $(G_1, G_2, … , G_n)$ 上。

算法原理

基本原理

谱聚类的 目的 就是要找到一种合理的分割，使得分割后形成若干子图，连接不同子图的边的权重尽可能低，同一子图内边的权重尽可能高。

推导过程

具体过程可通过以下操作来实现。

构造矩阵

Step 3.1.1中根据对称阵 $C$ 构造的矩阵 $W$，也是一个对称阵。它描述了 G 中各节点间的相似度。

注意：在 Step 1 中构造的矩阵 $C = (c_{ij})$ 中，$c_{ij}$ 表示顶点 $i$ 到顶点 $j$ 的距离，$c_{ij}$ 值越大，说明距离越远。

但是到了矩阵 $W$ 中的节点：$w_{ij} = \exp{(-\frac{c_{ij}^2}{2\sigma^2})}$。$c_{ij}$ 越大，则 $w_{ij}$ 越小。

也就是说 $W$ 中的节点 $w_{ij}$ 的数值越小，它所表示的对应的两个点之间的距离也就越大。

Step 3.1.2则是构造了 $W$ 的对角矩阵 $D$。

Step 3.1.3中，由相似度矩阵 $W$ 和其对角矩阵 $D$，我们构造了一个新的矩阵：$L= D-W$。

$L$ 是一个拉普拉斯（Laplacian）矩阵，称作非规范化的拉普拉斯矩阵（Unnormalized Laplacian Matrix)。

拉普拉斯矩阵性质

因拉普拉斯矩阵性质得知：

（i）$L$ 是对称半正定矩阵；

（ii）Laplacian 矩阵 $L$ 的最小特征值是 $0$，相应的特征向量是 $I$；

（iii）Laplacian 矩阵 $L$ 有 $n$ 个非负实特征值：$0 \leqslant \lambda_1 \leqslant \lambda_2 \leqslant … \leqslant \lambda_n$。

又因为 $L = D - W$，对于任一实向量 $f$，都可以做如下计算：

$f’Lf = f’Df – f’Wf = \sum_{i=1}^{n}d_i f_i^2 - \sum_{i=1}^{n}\sum_{j=1}^{n} f_i f_j w_{ij} = \frac{1}{2}(\sum_{i=1}^{n}d_i f_i^2 – 2\sum_{i=1}^{n}\sum_{j=1}^{n} f_i f_j w_{ij} + \sum_{j=1}^{n}d_j f_j^2) = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij} (f_i – f_j)^2$

图分割和矩阵运算的关系

现在我们回过头来，看图切割这件事情。

将图切割成两个子图

假设我们把 $L$ 所对应的原图进行图切割，成两个新的图：$A$ 和 $ \overline{A}$。

也就是说，之前 $n\times n$ 矩阵 $L$ 所对应的 $n$ 个顶点被分为了两部分，一部分属于 $A$，另一部分属于它的补 $\overline{A}$。

到底哪些点被分给了$A$，哪些点被分给了它的补呢？

我们可以用一个向量来表示——假设存在一个向量 $f = (f_1, f_2, ..., f_n)^T$，其中不同维度的值可以指示该维度对应的顶点属于新分出来的哪个子图。

具体如下：

$f_i = \begin{cases}\begin{matrix}

\sqrt{\frac{|\overline{A}|}{|A|}}, & if \ vi \in A \\

-\sqrt{\frac{|A|}{|\overline{A}|}}, & if \ vi \in \overline{A}

\end{matrix}\end{cases}$

将 $f$ 带入到上节（iii）中的公式：$f’Lf = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij} (f_i – f_j)^2$；

又因为当 $i$、$j$ 同属于 $A$ 或者 $\overline{A}$ 时，$f_i – f_j = 0$。

则 $f'Lf$ 就可以被转化为如下形式：

$\frac{1}{2}\sum_{i \in A, j \in \overline{A}} w_{ij}(\sqrt{\frac{|\overline{A}|}{|A|}} + \sqrt{\frac{|A|}{|\overline{A}|}})^2 + \frac{1}{2}\sum_{i \in \overline{A}, j \in A} w_{ij}(- \sqrt{\frac{|\overline{A}|}{|A|}} - \sqrt{\frac{|A|}{|\overline{A}|}})^2$

$=> \frac{1}{2}\sum_{i \in A, j \in \overline{A}} w_{ij}(\frac{|\overline{A}|}{|A|} + 2 + \frac{|A|}{|\overline{A}|}) + \frac{1}{2}\sum_{i \in \overline{A}, j \in A} w_{ij}(\frac{|\overline{A}|}{|A|} + 2 + \frac{|A|}{|\overline{A}|})$

$=> (\frac{|\overline{A}|}{|A|} + 2 + \frac{|A|}{|\overline{A}|}) [\frac{1}{2}\sum_{i \in A, j \in \overline{A}} w_{ij} + \frac{1}{2}\sum_{i \in \overline{A}, j \in A} w_{ij}]$

式子1

$Cut(\cdot)$ 函数

取出上面 式子1 的后一部分：

$ [\frac{1}{2}\sum_{i \in A, j \in \overline{A}} w_{ij} + \frac{1}{2}\sum_{i \in \overline{A}, j \in A} w_{ij}] = \frac{1}{2}W(A,\overline{A}) + \frac{1}{2}W(\overline{A}, A) = \frac{1}{2}\sum_{i=1}^{k}W(A_i, \overline{A}_i)$

其中，$k$ 表示不同类别的个数，这里 $k=2$。

令 $ Cut(A, \overline{A}) = \frac{1}{2}\sum_{i=1}^{k}W(A_i, \overline{A}_i) $，这里的 $W(A，, \overline{A})$ 表示子图 $A$ 和 $\overline{A}$ 之间连通边的权重。

此处定义的 $Cut(\cdot)$ 函数，又可以被称为“ 截函数 ”。

当一个图被划分成为两个子图时，“ 截 ”指子图间的连接密度，即被切割后的子图之间，原本是连通状态（但在切割时被截断）的边的值加权和。

我们要找到一种分割，使得分割后，连接被分割出来的两个子图的边的权重尽可能低，即“ 截最小 ”。

因此， $Cut(\cdot)$ 函数 就是我们求取 图切割 方法的 目标函数 。

求解目标函数

$Cut(\cdot)$ 函数中的值就是 $wij$（顶点 $i$ 位于 $A$，顶点 $j$ 位于 $\overline{A}$）。$wij$ 越小，则对应的两点间的距离越大。

我们既然要让切割出来的结果使两个子图之间的加权距离尽量大，那么自然，我们就要求：

$ \min{Cut(A, \overline{A})} => \min{[\frac{1}{2}\sum_{i=1}^{2}W(A_i, \overline{A}_i) =\frac{1}{2}W(A,\overline{A}) + \frac{1}{2} W(\overline{A}, A)]} $

我们将 $Cut(\cdot)$ 函数带回到 式子1 中，得到结果如下：

$Cut(A,\overline{A})(\frac{|A| + |\overline{A}|}{|A|} + \frac{|A| + |\overline{A}|}{|\overline{A}|})$ $= (|A| + |\overline{A}|)( (\frac{cut(A,\overline{A})}{|A|} + \frac{cut(A,\overline{A})}{|\overline{A}|}))$

其中：

$(\frac{cut(A,\overline{A})}{|A|} + \frac{cut(A,\overline{A})}{|\overline{A}|}) = \sum_{i=1}^{k}\frac{cut(A_i, \overline{A_i}}{|A_i|} = RatioCut(A, \overline{A})$

因此：

$(|A| + |\overline{A}|) RatioCut(A, \overline{A}) = |V| RationCut(A, \overline{A})$

其中 $|V|$ 表示的是顶点的数目，对于确定的图来说是个常数。

由上述的推导可知，由 $f’Lf$ 推导出了 $RatioCut(\cdot)$ 函数。到此，我们得出了：

$f’Lf = |V| RatioCut(A, \overline{A})$

因为 $Cut(\cdot)$ 函数和 $RatioCut(\cdot)$ 函数相差的是一个常数，因此求 $Cut(\cdot)$ 的最小值就是求 $RatioCut(\cdot)$ 的最小值。

又因为 $|V|$ 是常数，因此我们求 $RatioCut(\cdot)$ 函数的最小值就是求 $f’Lf$ 的最小值。

到此时，图切割问题，就变成了求 $f’Lf$ 的最小值的问题。

通过求 $f’Lf$ 的最小值来切割图

假设 $λ$ 是 Laplacian 矩阵 $L$ 的特征值，$f$ 是特征值 $λ$ 对应的特征向量，则有：$Lf = λf$。

在上式的两端同时左乘 $f’$，得到：$f’Lf = λf’f$。

已知 $||f|| = n^{\frac{1}{2}}$，则 $f’f = n$，上式可以转化为：$f’Lf = λn$。

既然我们的目标是求 $\min{ f’Lf}$，那么我们只需求得最小特征值 $λ$。

由 Laplacian 矩阵的性质可知，Laplacian 矩阵的最小特征值为 $0$，相应的特征向量是 $I$。

向量 $I$ 中所有维度都为 $1$，无法将对应顶点分为两份。因此我们用 $L$ 第二小的特征值（也就是最小的非零特征值）来近似取 $RatioCut(\cdot)$ 的最小值（此处背后实际的理论依据是 Rayleigh-Ritz 理论。）

我们先求出 $L$ 第二小的特征向量 $f$，再通过如下变换，将 $f$ 转化为一个离散的指示向量。

对于求解出来的特征向量 $f = (f_1,f_2,…, f_n)^T$ 中的每一个分量 $f_i$，根据每个分量的值来判断对应的点所属的类别：

$\begin{cases}\begin{matrix}

vi \in A, & if f_i \geqslant 0 \\

vi \in \overline{A}, & if \ f_i < 0

\end{matrix}\end{cases}$

这也就是 Step 3.1.4 中描述的内容。

从切割成 $2$ 份到切割成 $k$ 份的推演

如果不是要一次将图切成 $2$ 份，而是要切成 $k$ 份，那么就要先求 $L$ 的前 $k$ 小个特征向量。

这 $k$ 个特征向量表示为 $f^{(1)}, f^{(2)}, …, f^{(k)}$。

由特征向量构成如下这样一个 $n \times k$ 的特征向量矩阵：

$\left( \begin{matrix}

f_1^{(1)} & f_1^{(2)} & ... & f_1^{(k)} \\

f_2^{(1)} & f_2^{(2)} & ... & f_2^{(k)} \\

... \\

f_n^{(1)} & f_n^{(2)} & ... & f_n^{(k)}

\end{matrix} \right)$

将特征向量矩阵中的每一行作为一个样本，利用 KMeans 聚类方法对其进行聚类。也就是对 $n$ 个 $k$ 维向量进行聚类，将其聚为 $k$ 个簇。

聚类完成之后，如果特征矩阵中的第 $i$ 个 $k$ 维向量被聚集到了第 $j$ 个簇中，则原本图中的第 $i$ 个点就被聚集到了第 $j$ 个簇中。

以上，就是根据非规范化拉普拉矩阵进行基于图切割的谱聚类的算法原理。

规范化拉普拉斯矩阵

$L$ 也可以被规范化，$D^{-\frac{1}{2}}L D^{-\frac{1}{2}}$ 就是 $L$ 的规范化形式。

$L' = D^{-\frac{1}{2}}L D^{-\frac{1}{2}}$ 又称为规范化的拉普拉斯矩阵（Normalized Laplacian Matrix）。

对于规范化的拉普拉斯矩阵，不能直接求其特征值和特征向量来做图切割。不过大体过程和思路与非规范化拉普拉斯矩阵一致，在此不赘述。

实例

将10个人的身高体重数据用谱聚类进行聚类：

from sklearn.cluster import SpectralClustering
import numpy as np
import math

X = np.array([[185.4, 72.6],
[155.0, 54.4],
[170.2, 99.9],
[172.2, 97.3],
[157.5, 59.0],
[190.5, 81.6],
[188.0, 77.1],
[167.6, 97.3],
[172.7, 93.3],
[154.9, 59.0]])

w, h = 10, 10;

    #构建相似度矩阵，任意两个样本间的相似度= 100 - 两个样本的欧氏距离
Matrix = [[100- math.hypot(X[x][0]- X[y][0], X[x][1]- X[y][1]) for x in range(w)] for y in range(h)]

sc = SpectralClustering(3, affinity='precomputed', n_init=10)
sc.fit(Matrix)

print('spectral clustering')   
print(sc.labels_)

输出为：

spectral clustering
[2 1 0 0 1 2 2 0 0 1]