聚类算法

聚类是将样本集合中相似的样本(实例)分配到相同的类，不相似的样本分配到不同的类。聚类时，样本通常是欧式空间中的向量，类别不是事先给定，而是从数据中自动发现，但类别的个数通常要预先给定。样本之间的相似度或距离由度量决定。如果一个样本只能属于一个类，则称为硬聚类；如果一个样本可以属于多个类，则称为软聚类。

在这一部分主要介绍四种聚类算法:

• $k$均值聚类
• 层次聚类
• 谱聚类
• DBSCAN

$k$均值聚类

算法思想

$k$均值聚类的思想非常直观:

如果一团样本分布的足够紧凑，那么就可以认为这些样本属于一类，而一团样本的紧凑程度可以通过这些样本距离该团样本中心的距离和来衡量，距离和越小则说明该团样本越紧凑。

算法推导

假设对于样本集$\mathcal{X}$, 初始时随机将样本分成$c$类，那么根据$k$均值的算法思想可以定义损失函数为:

$$J(e) = \sum_{i=1}^c \sum_{y \in \Gamma_i} ||y - m_i||^2$$

公式中$m_i$是目前划归第$i$类样本的均值。想要求得使$J(e)$最小的最优样本划分并不容易，是一个NP难的问题，因此在优化时采用贪心策略，通过迭代优化的形式来进行求解。在迭代过程中，每次移动一个样本点来使$J(e)$减小，为此我们需要计算将一个样本点$y$从$\Gamma_i$类转移到$\Gamma_k$类对损失函数的影响，因为损失函数是依类别进行计算然后累加的，不妨将$\Gamma_i$的类的损失记做$J_i$。在将$y$移动之后，$\Gamma_i$和$\Gamma_k$类均发生了变化，均值分别变为:

$$\begin{aligned} m'_i &= m_i + \frac{m_i - y}{N_i - 1} \\ m'_k &= m_k + \frac{y - m_k}{N_k + 1} \end{aligned}$$

据此，可以分别计算出$J’_i$和$J’_k$并用$J_i, J_k$表示:

$$\begin{aligned} J'_i &= J_i - \frac{N_i}{N_i - 1} ||y - m_i||^2 \\ J'_k &= J_k + \frac{N_k}{N_k + 1} ||y - m_k||^2 \end{aligned}$$

如果$J’_i > J’_k$，就将样本从$i$类移动到$k$类，由此便可以得到$k$均值聚类步骤:

$k$均值:

• Step1: 把样本初始划分成$C$类，并计算$J_e$
• Step2: 选择一个样本$y$，假设$y \in \Gamma_i$
• Step3: 若$N_i = 1$, 则转Step2
• Step4: 计算 $$\rho_j = \begin{cases} \frac{N_j}{N_j + 1} ||y - m_j||^2 & j \neq i \\ \frac{N_i}{N_i - 1} ||y - m_i||^2 & j = i \end{cases}$$
• Step5: 若对$\forall j, \rho_k \leq \rho_j$, 则将$y$从$\Gamma_i$移动到$\Gamma_k$中去。
• Step6: 修正$m_i,m_k$和$J_e$
• Step2: 若连续迭代$N$次(将样本遍历一遍)不变，则算法终止，否则回到Step2

算法分析

算法时间复杂度

假设算法迭代$t$次后收敛，样本个数为$n$，样本维数为$d$， 类别数目为$k$，则算法时间复杂度为:

$$O(tndk)$$

算法超参数

算法要求预先制定聚类类别个数$k$

局部极小

因为$k$均值算法本质上是贪心算法，因此很有可能陷入局部最优，可以通过更新聚类个数以及不同的初始化来一定程度上减轻局部最优，来获得较好的聚类效果。

层次聚类

算法思想

层次聚类首先将各个样本点划分为1类，然后将距离最近的两类进行合并，建立一个新的类，重复这个步骤直到聚类结果满足要求。
从定义来看，层次聚类算法是自下而上进行聚类，在应用该聚类算法时需要明确三个问题:

• 相似性度量的选择
• 合并规则，即如何定义两个类别的相似程度
• 停止条件

下面就这三个关键问题进行分析。

算法详述

相似性度量

任何一个聚类算法会面临相似性度量的选择问题，实际中常用的相似性度量有:

• 欧式距离
• 马氏距离
• 测地距离
• ……

具体选择什么样的度量来刻画相似度需要根据实际数据样本的分布情况和待解决的问题来确定。

合并规则

合并规则一般是类间距离最小，若记$\Delta(\Gamma_i, \Gamma_j)$为类间距离，$\delta(x,y)$,为$x$和$y$之间的距离，则类间距离一般有以下几种定义:

• 最近距离 $$\Delta(\Gamma_i, \Gamma_j) = \min_{y \in \Gamma_i, \tilde{y} \in \Gamma_j} \delta(y, \tilde{y})$$
• 最远距离 $$\Delta(\Gamma_i, \Gamma_j) = \max_{y \in \Gamma_i, \tilde{y} \in \Gamma_j} \delta(y, \tilde{y})$$
• 均值距离 $$\Delta(\Gamma_i, \Gamma_j) = \delta(m_i, m_j)$$

停止条件

如果不加停止条件，则最终所有样本都会归于一类，停止条件一般有2种选择:

• 类别个数: 当类别个数到达期望值$k$时，便停止合并，聚类结束。
• 类间距离: 当类间距离大于某个值时便停止合并，聚类结束。

算法流程

层次聚类算法流程总结如下:

层次聚类:
输入: $n$个样本组成的样本集合及样本之间的距离
输出: 对样本集合的一个层次化聚类
算法流程:

• 根据选定的相似性度量方法，计算$n$个样本两两之间的距离$d_{ij}$,得到距离矩阵$D$
• 构造$n$个类，每个类只包含一个样本
• 合并类间距离最小的两个类，构成一个新类，并计算新类与当前各类之间的距离。
• 若满足聚类终止条件，则聚类终止，否则返回第3步

谱聚类

算法思想

谱聚类算法是从图论中衍生出的一种算法，后来在聚类问题中得到了广泛的应用。它的主要思想是把所有的数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。

学习谱聚类算法主要需要理清两个问题:

• 无向权重图的相关定义
• 拉普拉斯矩阵
• 图的切割任务

下面就按这个思路来对谱聚类算法进行介绍。

无向权重图

对于我们的样本集$\mathcal{X} \subseteq \mathbb{R}^d$，可以将每一个样本点看作高维空间中的一个点，若将这些点以某种方式连接起来则构成了一个图，如果对于相连的两个节点$x$和$y$，若并不区分$x \rightarrow y$和$y \rightarrow x$，则构成了一个无向图。

但对于我们的聚类任务，我们期望距离较近的点它们之间的连接权重大一些，距离较远的样本点之间的权重小一些么，从这样的一个想法出发，我们有以下几种方式来构造相似性图:

• $k$近邻图: 若$v_i$是$v_j$的$k$近邻，则$v_i$和$v_j$之间存在一条边，边上权重设置为1。
• 对称$k$近邻图: 若两个点互为$k$近邻，则这两个点之间存在一条边。
• $\epsilon$-近邻图: 任意两个距离小于$\epsilon$的点之间存在一条边
• 全连接图: 相似性大于0的两个点之间均存在一条边，常用高斯核函数来进行相似性度量: $$W_{ij} = S_{ij} = \exp(- \frac{||x_i - x_j||^2}{2 \sigma^2})$$

下面给出一些符号定义:

• $w_{ij}$： 边权重
• 节点$vi \in V$的度: $d_i = \sum{j=1}^n w_{ij}$
• 度矩阵: $D = diag(d_1,d_2, \dots, d_n)$
• 加权邻接矩阵: $W = (w{ij}){i,j= 1,\dots,n}$

拉普拉斯矩阵

拉普拉斯矩阵的定义非常简单：

$$L = D - W$$

拉普拉斯矩阵具有很多很好的性质:

• 拉普拉斯矩阵是对称矩阵，该性质可由$D$和$W$的对称性质得到
• 对于任意的向量$f \in \mathbb{R}^n$,有: $$f^T L f = \frac{1}{2} \sum_{i,j=1}^n w_{ij} (f_i - f_j)^2$$
• $L$是半正定矩阵，该性质由上一条性质可以直接得出

另一种常用形式是拉普拉斯矩阵的归一化形式:

$$L_{rw} = D^{-1} L = I - D^{-1} W$$

为了保证$L_{rw}$也是对称形式，常采用另一种形式:

$$L_{rw} = D^{-\frac{1}{2}} W D^{-\frac{1}{2}}$$

这样，$L_{rw}$也具有原始拉普拉斯矩阵的一些好的性质，比如对称性，半正定等性质。

图的切割任务

在已经定义好权重无向图的情况下，现在我们考虑这样一个任务:

将已经得到的图进行划分，使得不同点集边的权重较小，而同一点集内的边权重较大。以这种思想来完成聚类任务。

直接想到的损失函数是如下形式:

$$\begin{aligned} cut(A_1,A_2,\dots, A_k) &:= \frac{1}{2} \sum_{i=1}^k W(A_i, \bar{A}_i) \\ W(A, B) &:= \sum_{i \in A, j \in B} w_{ij} \end{aligned}$$

但选择这样形式的损失函数，最后往往是有一类包含非常多的点， 其他类都是单独的点。以聚成两类为例，假设有12个点，若两类各6个点，则最终计算损失函数时，需要将$266=72$项求和，若一类11个点，一类1个点，则最终计算时只需要将$2111=22$项求和，这也是为何选用上述损失函数容易得到单独的点。
我们在聚类时，往往期望各个类别的样本点个数比较均衡，因此考虑在损失函数中对样本点较少的类别进行惩罚，由此便可得到两种损失函数定义:

$$\begin{aligned} RatioCut(A_1, A_2, \dots, A_k) &:= \frac{1}{2} \sum_{i=1}^k \frac{W(A_i, \bar{A}_i)}{|A_i|} = \sum_{i=1}^k \frac{cut(A_i, \bar{A}_i)}{|A_i|} \\ Ncut(A_1,A_2, \dots, A_k) &:= \frac{1}{2} \sum_{i=1}^k \frac{W(A_i, \bar{A}_i)}{Vol(A_i)} = \sum_{i=1}^k \frac{cut(A_i, \bar{A}_i)}{Vol(A_i)} \end{aligned}$$

其中$|A_i|$表示集合$A_i$中的元素个数，$vol(A_i)$表示集合$A_i$中元素度的总和。

$RatioCut$切图

若以Ratiocut作为目标函数，则优化问题可以写做:

$$\min_{A_1, A_2, \dots, A_k} RatioCut(A_1,A_2,\dots, A_k)$$

这个问题是很难直接求解的，因此牛人们考虑通过另一种方式来对$RatioCut$进行表示，我们引入指示向量$hj = { h_1, h_2, \dots, h_k}, j = 1,2, \dots, k$,对于任意一个向量$h_j$,它是一个$n$维向量，我们定义$h{ji}$为:

$$h_{ji} = \begin{cases} 0 & v_i \notin A_j \\ \frac{1}{\sqrt{|A_j|}} & v_i \in A_j \end{cases}$$

那么我们对于$h_i^T L h_i$, 有:

$$\begin{aligned} h_i^T L h_i &= \frac{1}{2} \sum_{j=1}^n \sum_{k=1}^n w_{jk}(h_{ij} - h_{ik})^2 \\ &= \frac{1}{2} [\sum_{j \in A_i, k \notin A_i} w_{jk} (\frac{1}{\sqrt{|A_i|}} - 0)^2 + \sum_{j \notin A_i, k \in A_i} w_{jk} (0 - \frac{1}{\sqrt{|A_i|}})^2] \\ &= \frac{1}{2} [\sum_{j \in A_i, k \notin A_i} w_{jk} \frac{1}{|A_i|} + \sum_{j \notin A_i, k \in A_i} w_{jk} \frac{1}{|A_i|}] \\ &= \frac{1}{2} (cut(A_i, \bar{A}_i) \frac{1}{|A_i|} + cut(\bar{A}_i, A_i)\frac{1}{|A_i|} ] \\ &= RatioCut(A_i, \bar{A}_i) \end{aligned}$$

因此若以$RatioCut$作为目标函数，则可以表示为:

$$RationCut(A_1,A_2, \dots, A_k) = \sum_{i=1}^k h_i^T L h_i = trace(H^T L H)$$

其中:

$$H = [h_1, h_2, \dots, h_k]$$

因此我们如果要最小化$RatioCut$切图可以近似等价于优化如下问题:

$$\arg \min_{H} tr(H^T L H) \quad s.t. H^TH = I， h_{ji} = \begin{cases} 0 & v_i \notin A_j \\ \frac{1}{\sqrt{|A_j|}} & v_i \in A_j \end{cases}$$

但因为在上面推导过程中我们对矩阵$H$的形式做了限制，这是一个离散优化问题，找到满足优化目标的$H$矩阵是一个NP难的问题，因此我们考虑将离散约束放松，求解近似问题:

$$\arg \min_{H} tr(H^T L H) \quad s.t. H^TH = I$$

该优化问题其实就是一个瑞丽熵优化问题，如果要聚成$k$类，只需要取矩阵$L$的前$k$小的特征值所对应的特征向量即可。同时注意到$L$矩阵有着很好的性质:

拉普拉斯矩阵的的最小特征值为0，特征值0的重数等于图$G$连通分量$A1, \dots, A_k$的个数，且相应的特征空间可由$\mathbf{1{A1}}, \dots, \mathbf{1{Ak}}$张成，其中$\mathbf{1{A_i}}$中与$A_i$相对应的分量为1，其余分量为0。
证明: 特征值0所对应的特征向量$f$应当满足:

$$0 = f^T L f = \frac{1}{2} \sum_{i,j=1}^n w_{ij} (f_i - f_j)^2$$

也就是说对于$w_{ij} \neq 0$(两点之间连通)的情况，应当有$f_i = f_j$。

而该特征向量刚好与上面最小割优化问题中$hi$的要求相吻合，但因为我们根据相似性度量所构造的图可能是一个全连接图或者连通分量的个数与我们的聚类个数并不相等，则最终得到的$H{n \times k}$很可能是如下的形式(以$n=5,k=3$为例,未对特征向量做归一化):

$$H = \begin{bmatrix} 0.95 & 0.01 & 0.01 \\ 0.96 & 0.01 & 0.01 \\ 0.01 & 0.95 & 0.01 \\ 0.01 & 0.95 & 0.01 \\ 0.01 & 0.01 & 0.95 \end{bmatrix}$$

因此在在降维之后，往往会对图分割向量$h_i$再做一次k-means聚类。若选用$NCut$作为优化目标，则最终可以推导出等价于优化如下问题:

$$\arg \min_{F} tr(F^T D^{-\frac{1}{2}} L D^{-\frac{1}{2}} F) \quad s.t. F^TF =I$$

所以最终只需要取归一化对称拉普拉斯矩阵的前$k$小的特征值所对应的特征向量组成$F_{n \times k}$，然后再使用k-means算法进行一次聚类就可以。

谱聚类算法流程总结

下面给出谱聚类算法流程:

谱聚类算法:
输入: 样本集$D = (x_1, x_2, \dots, x_n)$,相似矩阵生成方式，聚类类别$k$，输入聚类方法
输出: 簇划分$C = (A_1, A_2, \dots, A_k)$
步骤:

• 根据输入的相似矩阵生成方式构建样本的相似矩阵$S$
• 根据相似矩阵$S$构建邻接矩阵$W$和度矩阵$D$
• 计算拉普拉斯矩阵$L = D - W$或归一化拉普拉斯矩阵$L_{rw} = D^{-\frac{1}{2}} L D^{-\frac{1}{2}}$
• 计算(归一化)拉普拉斯矩阵最小的$k$个特征值所对应的特征向量
• 利用$k$个特征向量构成矩阵$F_{n\times k}$
• 每一行原始样本在$k$维空间的投影点，利用输入聚类方法(常用k-means)进行聚类
• 得到簇划分$C = (A_1, A_2, \dots, A_k)$

DBSCAN

前面介绍的几种聚类算法可以归为“基于距离的聚类”，即在某种距离度量下，距离较近的点划归一类。但对于某些数据，使用基于距离的聚类方法可能并不理想，比如下图:

我们人眼直觉上会将外圈(蓝色样本点)划分为一类，将内圈(红色样本点)划分为一类，我们会有这样的直觉是因为外圈和内圈的样本点是分别连续的，为了让计算机能够实现这种情况下样本的聚类，基于密度的聚类算法便诞生了。
密度聚类算法假设聚类结构能够通过样本分布的紧密程度确定。通常情形下，密度聚类算法从样本密度的角度来考察样本之间的可连接性，并基于可连接样本不断扩展聚类簇以获得最终的聚类结果。
DBSCAN是一种著名的密度聚类算法，它基于一组邻域参数($\epsilon, Minpts$)来刻画样本分布的紧密程度，给定数据集$D = (x_1, x_2, \dots, x_n)$, 定义下面几个概念:

• $\epsilon$-邻域: 对$x_j \in D$,其$\epsilon$-邻域中包含样本集$D$中与$x_j$的距离不大于$\epsilon$的样本，即: $$N_\epsilon(x_j) = \{ x_i \in D | dist(x_i, x_j) \leq \epsilon \}$$
• 核心对象: 若$xj$的$\epsilon$-邻域中至少包含$MinPts$个样本，即$|N\epsilon(x_j)| \geq MinPts$, 则$x_j$是一个核心对象。
• 密度直达: 若$x_j$位于$x_i$的$\epsilon$邻域中，且$x_i$是核心对象，则称$x_j$由$x_i$密度直达。
• 密度可达: 对$xi$与$x_j$，若存在样本序列$p_1,p_2,\dots, p_n$，其中$p_1 = x_i, p_n = x_j$,且$p{i+1}$由$p_i$密度直达，则称$x_j$由$x_i$密度可达。
• 密度相连: 对$x_i$与$x_j$,若存在$x_k$使得$x_i$与$x_j$均由$x_k$密度可达，则称$x_i$与$x_j$密度相连。

基于这些概念，DBSCAN将“簇”定义为:

由密度可达关系导出的最大的密度相连样本集合，形式化地说，给定邻域参数$(\epsilon, Minpts)$,簇$C \subseteq D$是满足以下性质的非空样本子集:

• 连接性: $x_i \in C, x_j \in C \Rightarrow$ $x_i$与$x_j$密度相连
• 最大性: $x_i \in C$, $x_j$由$x_i$密度可达 $\Rightarrow x_j \in C$

在给出了聚类簇的定义后，我们需要解决的问题便是如何从数据集$D$中找出满足以上性质的聚类簇。实际上，若$x$为核心对象，则$x$密度可达的所有样本组成的集合记为$X = { x’ \in D | x’由 x 密度可达 }$,容易证明$X$即为满足连接性与最大性的簇。

于是，DBSCAN算法先任选数据集中的一个核心对象作为“种子”，再由此出发确定相应的聚类簇，算法流程如下:

DBSCAN算法
输入: 样本集$D = (x_1,\dots,x_n)$,邻域参数$(\epsilon, MinPts)$
输出: 簇划分$C = (A_1, A_2, \dots, A_k)$
算法流程:

• 初始化核心对象集合$\Omega = \varnothing$
• 遍历所有样本点，将所有核心对象放入集合$\Omega$
• 初始化聚类簇数$k=0$, 初始化未访问样本集合$\Gamma = D$
• 随机从$\Omega$中取一个核心对象，然后基于该核心对象得到聚类簇$C_k$,更新: $$\Omega = \Omega \backslash C_k, \Gamma = \Gamma \backslash C_k$$
• 重复以上步骤直到$\Omega$为空

但需要注意的是，采用$DBSCAN$聚类方法，如果邻域参数$(\epsilon, MinPts)$选择不恰当会导致一些样本点不会被归到任何类别。