六种常见聚类算法

 目录

Kmeans

 DBSCAN-基于密度的空间聚类算法

谱聚类

ＧＭＭ-高斯混合模型

 MeanShift-均值迁移

层次聚类

 代码

Kmeans

聚类原则：以空间中k个点为中心进行聚类，对最靠近他们的对象归类。逐次计算各簇中心的值为新的中心值，迭代更新，直至簇中心位置不再改变或者达到最大迭代次数。

 Kmeans的目标函数

 定义为：各簇成员到其簇首的距离的平方和最小，如下所示

式中，C为簇首(聚类中心)集合，共有K个簇首。计算目标函数梯度，令梯度为0，计算簇首C，

 式中l(x)表示簇成员个数。通过迭代优化目标函数来计算最佳参数C。由上式得，在每次迭代中需更新聚类中心为每个簇的簇心即簇成员的均值。

算法流程：

1. 适当选择k个类的初始中心；
2. 在第n次迭代中，对任意一个样本，求其到k个中心的距离，将该样本归到距离最短的中心所在的类/簇；
3. 利用均值等方法更新该类的中心值；

4. 对于所有的k个聚类中心，如果利用（2）（3）的迭代法更新后，聚类中心的位置保持不变，则迭代结束；否则，则继续迭代。

优点：

速度快，简单

缺点：

1. 适合聚类球状类簇，不能发现一些混合度较高，非球状类簇
2. 需指定簇个数
3. 算法结果不稳定，最终结果跟初始点选择相关，容易陷入局部最优

4. 对噪声或离群点比较敏感：无法区分出哪些是噪声或者离群点，只能给每个数据点都判断出一个类别来，这样会导致样本质心偏移，导致误判或者聚类紧密程度降低。

• kmeans算法结果易受初始点影响
  由于kmeans算法结果易受初始点影响，得到的结果是局部最优，为次，我们可以运行n次kmeans算法，每次运行的初始点均为随机生成。然后在n次运行结果中，选取目标函数(代价函数)最小的聚类结果。
• 聚类数量k如何确定?
  通常在数据集有多少个聚类是不清楚的。对于低维的数据，对其可视化，我们可以通过肉眼观察到有多少个聚类，但这样并不准确，而且大部分数据也无法可视化显示。方法1：我们可以通过构建一个代价函数与聚类数量k的关系图，如下图所示，横坐标是聚类数量，每个聚类数量对应的代价函数J均是n次运行的最优结果。观察下图，代价函数随聚类数量增大而减小，但并不是说聚类数量越多越好，比如：对于m个样本，m个聚类的代价函数自然是0，但这样根本不合理。观察下图，代价函数在没有到拐点之前下降很快，拐点之后下降变缓，我们就可以考虑选择这个拐点对应的数量作为聚类数量，这是一种较为合理的方法，但受限于能否得到一个清晰的拐点，如右下图所示，
•                            
  方法2：结果导向。我们使用kmeans算法的目的是什么？来确定k值能更好的服务于后续的目的。根据业务需求，多少个类别能够满足需求，效果更好来确定聚类数量。

 DBSCAN-基于密度的空间聚类算法

 聚类原则：聚类距离簇边界最近的点

算法流程：

核心点：核心点的半径范围内的样本个数≥最少点数

边界点：边界点的半径范围内的样本个数小于最少点数大于0

噪声点：噪声点的半径范围的样本个数为0

1. 根据半径、最少点数区分核心点、边界点、噪声点
2.先对核心点归类：
     while:
           随机选取一核心点作为当前簇
           寻找距离当前簇最近的核心点:与簇边缘最近的核心点，，
           if 若该核心点距离当前簇的边界核心点的距离小于半径:
              则将该核心点归类到当前簇
           if 若剩余的未归类的核心点距离当前簇边界距离均大于半径:
              说明：距离第numC个簇最近的核心点不在该簇邻域内，说明第numC个簇已全部分类完毕，可以 
              生成新簇来归类核心点,则在剩余的未归类的核心点随机选取一核心点归类到新的簇中
           if 核心点已全部归类：
              停止迭代
3. 再根据半径和已归类的核心点来归类边界点

优点：

能够识别任意形状的样本. 该算法将具有足够密度的区域划分为簇，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。

缺点：

需指定最少点个数，半径(或自动计算半径)

谱聚类

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

 算法流程：

1. 生成邻接矩阵(权重矩阵/相似度矩阵)：W:m*m，用高斯函数度量样本之间的距离，距离越远，权重越小

   

2. 生成度量矩阵D：m*m，D=np.diag(np.sum(W,axis=1))、拉普拉斯矩阵L:L=D-W
3. 标准化拉普拉斯矩阵:Ncut切图:L_norm=D^(-0.5)LD^(-0.5),注意：这里是矩阵叉乘

   

4. 生成标准化拉普拉斯矩阵对应的特征值和特征向量，对特征值排序(升序)，选取前k个最小的特征值对应的特征向量，组成m*k矩阵data
5. 对降维后的矩阵data使用kmeans、GMM或其他聚类方法来聚类

拉普拉斯映射

        简单说，谱聚类其实就是使用了一种映射方法，将原有数据映射到新的数据空间，使得原数据中空间位置接近的样本在映射后更加接近，这种映射称为拉普拉斯映射。注意：步骤4中是选取最小的k个特征值对应特征向量，那么问题来了，为什么要选取最小特征值呢？特征值有什么物理意义呢？

        对于一个方阵(L为对称矩阵也是方阵且是半正定矩阵)，可以称之为变换矩阵，当它叉乘某一个向量时，其物理意义为对该向量作旋转和缩放的变换，即将原有向量所在的空间映射到一个新的空间。那么，缩放尺度如何度量？可以用变换矩阵的特征值来度量，如下所示A为变换矩阵，向量α为A的特征向量(对于特征向量，对其变换时，只做缩放变换，其向量方向不变，更多详细内容见看图就懂：线性代数之特征值分解与奇异值分解)，λ为特征值，I为单位向量，向量α左乘A等于对特征向量缩放λ倍。

特征值在物理意义上可以表示为缩放倍数，特征值越大，对应的特征向量放大倍数越大，而特征向量越大，对应的样本点映射后在空间位置上也更分散，而谱聚类或者拉普拉斯特征映射的目标就是将原本比较接近的点在映射后更加接近，使得属于不同簇的样本的可分离性变强，所以拉普拉斯特征映射自然要选择小特征值对应的空间，这样也就实现了原数据中空间位置接近的样本在映射后更加接近。如下图是原始特征空间与拉普拉斯特映射后的样本空间(2d、3d显示)的聚类结果

PCA 

        经过拉普拉斯映射后样本的可分离性比较好，易于聚类。于此相对的是，PCA降维选用的降序特征值即选取最大的k个特征值对应的特征向量，因为PCA的目的是对原数据作旋转变换，将原数据映射到特征基(关于特征基概念见看图就懂：线性代数之特征值分解与奇异值分解)，投影到特征基方向上的样本相对比较分散，如下图所示，二维原数据旋转到二维特征基上和降到1维效果图

 总之，特征值越大对应的向量空间方差越大，即向量方向上的数据越分散。降维数据肯定会丢失一部分信息，PCA适用于处理高维数据，可以在尽可能保留更多信息的同时对高维数据降维处理，如何保留更多数据信息呢？某种意义上，数据分布越分散(方差)或者说越混乱(熵)，能得到的信息越多，PCA通过将数据映射到特征基上(特征基不止一个)，然后从中选取k个方差最大的特征基(k个最大特征值对应的特征向量)组成新的特征空间，从而实现降维+保留大部分信息效果。

拉普拉斯映射与PCA区别

设有原数据X，为m行n列数据、m个样本n个特征，k为维度个数：

1. 拉普拉斯映射是根据样本的相似矩阵作变换得到新的样本空间，这个新样本空间可以直接用来聚类。
2. PCA是根据特征(变量)的协方差矩阵作变换得到一个变换矩阵，然后使用该变换矩阵映射原数据得到：
3. 拉普拉斯映射提取升序特征值对应的特征向量，PCA提取降序特征值对应的特征向量

优点：

谱聚类能够识别任意形状的样本空间且收敛于全局最优解，其基本思想是利用样本数据的相似矩阵(拉普拉斯矩阵)进行特征分解后得到的特征向量进行聚类。

谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好

缺点：

如果最终聚类的维度非常高，则由于降维的幅度不够，算法的运行速度和最终聚类效果均不好

聚类效果依赖于相似度矩阵(权重矩阵/邻接矩阵)，不同的相似度矩阵得到的最终聚类效果可能差别很大

需指定聚类个数，高斯参数。

ＧＭＭ-高斯混合模型

        高斯混合混合模型(GMM)是一种概率生成模型，模型通过学习先验分布后推导后验分布来实现聚类。当我们做聚类任务的时候，可以根据距离最近原则，将样本聚类到距离其最近的聚类中心，如k-means算法，或者将样本聚类到距离某个簇边缘最近的簇，如DBSCAN算法，而GMM是假设样本每个类的特征分布遵循高斯分布，即整个数据集可由不同参数的高斯分布线性组合来描述，GMM算法可以计算出每个样本属于各个簇的概率值并将其聚类到概率值最大的簇。

GMM如何实现？

        GMM是以假设数据各个类别的特征分布是高斯分布为基础，所以首要步骤就是先得到各个类别的概率分布密度函数参数，对于一元变量，需要计算各个类别的均值和期望，而对于多元变量，则需要计算均值和协方差矩阵。有如下数据集X，共有m个样本，n项特征，k个类别

 高斯概率分布密度函数公式如下：

 一元变量：n=1

  多元变量：n>1

 式中，表示协方差矩阵，n行n列矩阵，表示计算第k个类别的协方差矩阵的行列式。

现有如下数据集，如何使用GMM聚类？

在学习GMM聚类算法前，先从大学概率论课本上的一个例子入手。

例如：投掷一颗骰子，实验的样本空间为：S={1,2,3,4,5,6}。

存在以下事件

• 事件X={投掷一个骰子出现i点}，
• 事件={投掷一个骰子出现偶数点}，
• 事件= {投掷一个骰子出现奇数点}，

是样本空间S的一个划分，每做一次实验时，事件必有一个且仅有一个发生。则有

证明过程如下：

则事件X={投掷一个骰子出现i点}的概率p(X)为

  式中的称为先验概率，进一步得到以下公式：贝叶斯公式 

综上我们可以得到，投掷一颗骰子后，发生事件的概率为 

回到聚类任务中，上述 事件X={投掷一个骰子出现i点}其实就是我们的训练集对应X，事件

就是类别，对应y。如下所示：

p(X)是训练集样本出现的概率，公式如下：

GMM代价函数

         GMM首先需要计算每个类别的分布函数的参数，如何估计参数呢？可以采用极大似然估计方法。

        极大似然估计，就是利用已知的样本信息，反推最具有可能(最大概率)导致这些样本出现的模型参数值。换句话说，极大似然估计提供了一种通过给定观察数据来估计模型参数的方法，即：模型已定，参数未知。可以这样理解，既然事情(样本)已经发生了，为什么不让这个结果出现的可能性最大呢？这也就是极大似然估计核心。

        通过极大似然估计方法来使样本出现的概率p(X)最大从而得到分布函数参数，这样我们可以得到代价函数

 由于连乘可能会导致下溢，所以需要取对数。对数函数是一个严格递增的函数，取对数后不会改变数据的相对关系，即对数化后，我们仍然可以获得具有相同临界点的最优化函数，代价函数如下所示

 已知多元变量的高斯分布公式如下

可以得到

注意，式中  ，这里直接用概率分布密度函数表示概率分布可能有一定误导性。为什么这样表达呢？如下 

  是每个类的一个常量乘法因子，在之后的对后验概率进行规范化的时候就抵消掉了。因此可以直接使用来估计类条件概率。

 综上，GMM代价函数的最终公式如下

           通过求解代价函数，就可以得到各个类别的高斯分布函数参数，进而能计算样本属于某个类别的概率，从而实现聚类，公式如下

 以上过程其实就是一个建模过程，建模过程包括如下

1. 分析业务类型(聚类问题)，√
2. 确定模型类型(非参数+概率+生成模型)，√
3. 确定代价函数，√
4. 目标函数：超参数
5. 优化模型：调参

        在得到代价函数后，我们还要确定这个代价函数是否具有凸性，因为凸函数有一个好性质，若代价函数是凸函数，则其局部最优解也是全局最优解。若无法得到凸代价函数，则尝试更换模型或者使用迭代算法来逼近全局最优解。

目标函数一般是代价函数+正则项，用来确定超参数，防止过拟合，来提高模型泛化能力。本文这里先跳过这一步，直接讲如何计算参数。

使用EM算法计算GMM参数

        EM(期望最大)算法是一个迭代算法，其思想就是多次迭代来估计参数，直到算法收敛。EM算法就是一个迭代逼近过程，若目标函数为非凸函数，则EM算法不能保证收敛至全局最优解，且EM算法的求解结果与初值的选择有较大关系，该算法对初值敏感。

        EM算法求解参数思路可以将其理解为对两个未知变量的方程求解。首先固定其中一个未知数，求另一个未知数的偏导数，之后再反过来固定后者，求前者的偏导数，反复迭代直到目标函数值不变、收敛。EM算法就是相互迭代，求出一个稳定值，用于在无法直接找到参数的情况下寻找模型的最大似然估计(MLE)。它包括两个步骤:期望步骤和最大化步骤。在聚类问题中，步骤如下

1. 初始化参数：混合高斯分布函数参数
2. 根据估计的参数计算簇成员后验概率：
3. 最大化步骤：代价函数代入，求参数偏导，更新每一个参数
   (分别可以只用一个未知量来表达)
4. 重复1、2步骤直到对数似然值收敛。 

 具体步骤如下

1、初始化参数，根据估计参数计算给定数据点属于类别的概率(后验概率)：

2、最大化步骤：对代价函数求偏导， 根据第一步得到的 来更新每一个参数,

 3、重复1、2步骤，直到算法收敛。

参数推导过程 

        可通过对代价函数求偏导得到，过程比较复杂，所以本文从另一个角度来解释参数的计算过程。

关于均值μ：

 关于协方差矩阵如何计算：给定如下训练集X:m*2，假设有k个类别，则

        GMM的代价函数是最大化数据X出现的概率p(X)。若多次迭代后的p(X)保持不变，则说明程序已收敛。在得到参数，即得到数据集每个类别的概率分布函数，就可以计算每个样本属于各个类别的概率，样本选择最大的概率值对应的聚类中心，从而实现聚类。

GMM算法优缺点

优点：

GMM假设生成数据的是一种混合的高斯分布，为模型找到最大似然估计。与将数据点硬分配到聚类的K-means方法（假设围绕质心的数据呈圆形分布）相比，它使用了将数据点软分配到聚类的方法（即概率性，因此更好）。速度：它是混合模型学习算法中的最快算法;无偏差性：由于此算法仅最大化可能性，因此不会使均值趋于零，也不会使聚类大小具有可能适用或不适用的特定结构。
（A）通过使用软分配捕获属于不同聚类的数据点的不确定性，（B）对圆形聚类没有偏见。即使是非线性数据分布，它也能很好地工作，具有较强的鲁棒性

缺点：

• 奇异性：当每个混合模型的点数不足时，估计协方差矩阵将变得困难，并且除非对协方差进行人为正则化，否则该算法会发散并寻找无穷大似然函数值的解。GMM:需要求逆

• 分量数量：该算法将始终使用它可以使用的所有分量，所以在没有外部提示的情况下，需要留存数据或者信息理论标准来决定用多少个分量。

• 适合聚类球状类簇，不能发现一些混合度较高，非球状类簇

• 需要指定簇个数

 MeanShift-均值迁移

 该算法也叫做均值漂移，在目标追踪中应用广泛。本身其实是一种基于密度的聚类算法。

1、算法思路：

主要思路是：计算某一点A与其周围半径R内的向量距离的平均值M即迁移向量，更新该点下一步的迁移位置（A=M+A）。当该点不再移动时(norm(Mnew-Mold)<Th)，该点与周围半径R内点形成一个类簇，

2、如何计算迁移向量M

step1:随机选择一点为初始簇心点（不完全随机：不能把噪声点/离群点作为初始簇心，以防乒乓效应）

 step2: 根据初始半径R，选择以R为半径的圆的范围内的点分别与簇心点的向量的平均值作为迁移向量M

 表示：属于当前簇c内的点集合，簇c：半径为R，簇心为xc,簇内样本点数量为k。

3、如何迁移

    根据迁移向量M，更新簇心

4、如何收敛

迁移向量M前后变化量小于阈值且样本点(除了离群点)均已归类，则停止运行

5、如何成簇

norm(Mnew-Mold)>TH时:
    判断norm(M)<Th?即当前迁移距离是否大于阈值：
           当迁移距离接近0时，此时根据当前簇内的点扩大簇半径Rnew
           根据扩大后的半径Rnew，计算新的迁移向量M
           判断此时的迁移向量前后变化：
                若norm(Mnew-Mold)<TH:说明 当前簇已稳定，可以判断是否生成下一个簇了。
                      判断是否还有未归类的样本(非离群点)：
                           若还有，则重新初始化半径R，在剩余点内继续分簇；
                           若无，则停止运行。
                若norm(Mnew-Mold)>TH:说明 当前簇未成形，需要继续迁移
norm(Mnew-Mold)<TH时:     
    判断是否还有未归类的样本(非离群点)：
        若还有，则重新初始化半径R，在剩点内继续分簇；
        若无，则停止运行。

6、如何调整簇半径大小？

当迁移向量长度接近0时，需要扩大簇半径即扩大搜索范围，看看扩大范围后，簇内点是否增加，即迁移向量增量是否大于阈值Th，若小于阈值，则没必要继续迁移(当前簇已稳定)，反之继续进行迁移计算。

关于调整簇半径策略：根据当前簇成员调整簇半径

Rold-->当前半斤 Rold范围内簇成员集合：

 这里max()表示：如a=[[3,5],[9,2]],max(a)=[9,5]

7、扩展

计算均值的迁移向量-->计算核函数形式的迁移向量

 优点：

与kmeans相比，均值迁移无需指定簇个数，自动发现簇个数，需指定初始半径(带宽)(或者算法自动计算)。算法相对k-means分类结果比较稳定

缺点：

算法分类结果易受噪声点(离群点)干扰，若初始簇心点为离群点，可能会发生乒乓效应，算法无限循环，需数据预处理(过滤噪声点)或者改变距离计算方法。

适合聚类球状类簇，不能发现一些混合度较高，非球状类簇

计算复杂度高

R半径非自适应的menshift聚类算法聚类结果取决于半径(带宽)的设置，半径太小，收敛太慢，簇类个数过多；半径太大，一些簇类可能会丢失。

层次聚类

训练数据X，

例如有下数据

各个样本点之间的距离：

 算法流程：

1、每个样本都视为一个簇

2、计算各个样本点之间的距离(相似度)，生成m*m的距离矩阵dis

3、寻找最近的两个簇(寻找矩阵dis中的非0元素中的最小值对应的索引，i1，i2)，将它们归为一类(计算两个点/簇的均值means，如何归为同一类？-->X[[i1,i2],:]=menas)，本算法中衡量距离方法：计算样本点到簇心的距离-> average 

a=copy.deepcopy(X)
i,j = np.unravel_index(np.argmin(dis), dis.shape)#获取多维矩阵xx值对应的索引
listIJ = findRow(a,i1,i2)#在X中寻找与a[i1,:],a[i2,:]样本点相同的所有样本点对应的索引
means = np.mean(X[listIJ,:],axis=0)#计算最近的样本点/簇的均值
a[listIJ,:] = means#赋予相同的均值即视为归为同一类
Tree[iter] = copy.deepcopy(a)#生成聚类树

def findRow(inX,i,j):
    listI,listJ= [],[]
    for m in range(inX.shape[0]):
        if sum(abs(inX[m,:]-inX[i,:]))<1e-10:
            listI.append(m)
        elif sum(abs(inX[m, :] - inX[j,:]))<1e-10:
            listJ.append(m)
    listI.extend(listJ)
    return listI

4、重复步骤2、3，直到所有样本归为一类，或者归为指定类别个数k

a=copy.deepcopy(X)
while len(set(a[:,0]))>1：#while len(set(a[:,0]))>k
    i,j = np.unravel_index(np.argmin(dis), dis.shape)#获取多维矩阵xx值对应的索引
    listIJ = findRow(a,i1,i2)#在X中寻找与a[i1,:],a[i2,:]样本点相同的所有样本点对应的索引
    means = np.mean(X[listIJ,:],axis=0)#计算最近的样本点/簇的均值
    a[listIJ,:] = means#赋予相同的均值即视为归为同一类
    Tree[iter] = copy.deepcopy(a)#生成聚类树
#聚类树：每一层都是一个分类结果，簇数量由多到少。

优点：

一次性得到聚类树，每一层均为分类结果。算法相对k-means分类结果比较稳定。

相似度规则(距离衡量)容易定义

可以发现类别的层次关系

缺点：

计算复杂度高，不适合数据量大的

由于类聚类之间的距离(相似度)衡量方法不同，算法分类结果易受噪声点(离群点)干扰，需数据预处理(过滤噪声点)或者改变距离计算方法。

适合聚类球状类簇，不能发现一些混合度较高，非球状类簇