聚类

scikit-learn中的聚类算法的比较

方法名称	参数	可扩展性	使用场景	几何图形(公制使用)
K-Means	聚类形成的簇的个数	非常大的n_samples,中等的n_clusters使用MiniBatch代码	通用,均匀的cluster size,平面几何,不是太多的clusters	点之间的距离

一、K-means

k-means算法源于信号处理中的一种向量量化方法,现在则更多地作为一种聚类分析方法流行于数据挖掘领域。k-means聚类的目的是:把$n$个点(可以是样本的一次观察或一个实例)划分到$k$个聚类中,使得每个点都属于离他最近的均值(此即聚类中心)对应的聚类,以之作为聚类的标准。这个问题将归结为一个把数据空间划分为Voronoi cells的问题。该算法需要指定簇的数量。

这个问题在计算上是NP困难的,不过存在高效的启发式算法。一般情况下,都使用效率比较高的启发式算法,它们能够快速收敛于一个局部最优解。这些算法通常类似于通过迭代优化方法处理高斯混合分布的最大期望算法(EM算法)。而且,它们都使用聚类中心来为数据建模;然而k-means聚类倾向于在可比较的空间范围内寻找聚类,期望-最大化技术却允许聚类有不同的形状。k-means聚类与k-近邻之间没有任何关系。

K-means相关问题描述:

已知观测集$(x_1,x_2,\cdots,x_n)$,其中每个观测都是一个$d$维实向量,k-means聚类要把这$n$个观测划分到$k$个集合中$(k\le n)$,使得组内平方和(WCSS within-cluster sum of squares)最小。换句话说,它的目标是找到使得下式满足的聚类$S_{i}$，

其中$\mu_{i}$是$S_{i}$中所有点的均值。

• 标准算法

最常用的算法使用了迭代优化的技术。它被称为k-means算法而广为使用,有时也被称为Lloyd算法(尤其在计算机科学领域)。已知初始的$k$个均值点$m_1^{(1)},\cdots,m_k^{(1)}$,算法的按照下面两个步骤交替进行:

1. 分配(Assignment):将每个观测分配到聚类中,使得组内平方和(WCSS)达到最小。因为这一平方和就是平方后的欧氏距离,所以很直观地把观测分配到离它最近的均值点即可。(数学上,这意味依照由这些均值点生成的Voronoi图来划分上述观测)。

$$S_{i}^{(t)}=\left\{x_{p}:\left\|x_{p}-m_{i}^{(t)}\right\|^{2}\leq \left\|x_{p}-m_{j}^{(t)}\right\|^{2}\forall j,1\leq j\leq k\right\}$$

其中每个$x_{p}$都只被分配到一个确定的聚类$S^{{t}}$中，尽管在理论上它可能被分配到2个或者更多的聚类。

1. 更新(Update):对于上一步得到的每一个聚类,以聚类中观测值的图心,作为新的均值点。

$$m_{i}^{(t+1)}={\frac {1}{\left|S_{i}^{(t)}\right|}}\sum_{x_{j}\in S_{i}^{{(t)}}}x_{j}$$

因为算术平均是最小二乘估计，所以这一步同样减小了目标函数组内平方和（WCSS）的值。

这一算法将在对于观测的分配不再变化时收敛。由于交替进行的两个步骤都会减小目标函数WCSS的值,并且分配方案只有有限种,所以算法一定会收敛于某一(局部)最优解。注意:使用这一算法无法保证得到全局最优解。

这一算法经常被描述为"把观测按照距离分配到最近的聚类"。标准算法的目标函数是组内平方和(WCSS),而且按照"最小二乘和"来分配观测,确实是等价于按照最小欧氏距离来分配观测的。如果使用不同的距离函数来代替(平方)欧氏距离,可能使得算法无法收敛。然而,使用不同的距离函数,也能得到k-means聚类的其他变体,如球体k-均值算法和k-中心点算法。

• 初始化方法

通常使用的初始化方法有Forgy和随机划分(Random Partition)方法。Forgy方法随机地从数据集中选择k个观测作为初始的均值点;而随机划分方法则随机地为每一观测指定聚类,然后运行"更新(Update)"步骤,即计算随机分配的各聚类的图心,作为初始的均值点。Forgy方法易于使得初始均值点散开,随机划分方法则把均值点都放到靠近数据集中心的地方。此随机划分方法一般更适用于k-调和均值和模糊k-均值算法。对于期望-最大化(EM)算法和标准k-均值算法,Forgy方法作为初始化方法的表现会更好一些。

这是一个启发式算法,无法保证收敛到全局最优解,并且聚类的结果会依赖于初始的聚类。又因为算法的运行速度通常很快,所以一般都以不同的起始状态运行多次来得到更好的结果。不过,在最差的情况下,k-均值算法会收敛地特别慢:尤其是已经证明了存在这一的点集(甚至在2维空间中),使得k-均值算法收敛的时间达到指数级($2^{\Omega (n)}$)。好在在现实中,这样的点集几乎不会出现:因为k-均值算法的平滑运行时间是多项式时间的。

注:把“分配”步骤视为“期望”步骤,把“更新”步骤视为“最大化步骤”,可以看到,这一算法实际上是广义期望-最大化算法(GEM)的一个变体。

• 复杂度

在$d$维空间中找到k-means聚类问题的最优解的计算复杂度:

NP-hard:一般欧式空间中,即使目标聚类数仅为2

NP困难:平面中,不对聚类数目k作限制

如果k和$d$都是固定的,时间复杂度为$O(n^{{dk+1}}\log n$),其中$n$为待聚类的观测点数目。相比之下,Lloyds算法的运行时间通常为$O(nkdi)$,$k$和$d$定义如上,$i$为直到收敛时的迭代次数。如果数据本身就有一定的聚类结构,那么收敛所需的迭代数目通常是很少的,并且进行少数迭代之后,再进行迭代的话,对于结果的改善效果很小。鉴于上述原因,Lloyds算法在实践中通常被认为几乎是线性复杂度的。

下面有几个关于这一算法复杂度的近期研究：

Lloyd's k-means算法具有多项式平滑运行时间。对于落在空间$[0,1]^{d}$任意的$n$点集合,如果每一个点都独立地受一个均值为$0$,标准差为$\sigma$的正态分布所影响,那么k-means算法的期望运行时间上界为$O(n^{{34}}k^{{34}}d^{{8}}log^{4}(n)/\sigma ^{6})$,即对于n,k,i,d和$1/\sigma$都是多项式时间的。

在更简单的情况下,有更好的上界。例如在整数网格$\left\{1,...,M\right\}^{d}$中,k-means算法运行时间的上界为$O(dn^{4}M^{2})$。

使得k-means算法效率很高的两个关键特征同时也被经常被视为它最大的缺陷:

聚类数目k是一个输入参数。选择不恰当的k值可能会导致糟糕的聚类结果,这也是为什么要进行特征检查来决定数据集的聚类数目了。

收敛到局部最优解,可能导致"反直观"的错误结果。

k-means算法的一个重要的局限性即在于它的聚类模型。这一模型的基本思想在于:得到相互分离的球状聚类,在这些聚类中,均值点趋向收敛于聚类中心。一般会希望得到的聚类大小大致相当,这样把每个观测都分配到离它最近的聚类中心(即均值点)就是比较正确的分配方案。k-means聚类的结果也能理解为由均值点生成的Voronoi cells。

• 相关应用

k-means聚类(尤其是使用如Lloyd's算法的启发式方法的聚类)即使是在巨大的数据集上也非常容易部署实施。正因为如此,它在很多领域都得到的成功的应用,如市场划分、机器视觉、 地质统计学、天文学和农业等。它经常作为其他算法的预处理步骤,比如要找到一个初始设置。

1. 向量的量化 :k-means起源于信号处理领域,并且现在也能在这一领域找到应用。例如在计算机图形学中,色彩量化的任务,就是要把一张图像的色彩范围减少到一个固定的数目k上来。k-means算法就能很容易地被用来处理这一任务,并得到不错的结果。其它得向量量化的例子有非随机抽样,在这里,为了进一步的分析,使用k-means算法能很容易的从大规模数据集中选出k个合适的不同观测。
2. 聚类分析 :在聚类分析中,k-means算法被用来将输入数据划分到k个部分(聚类)中。然而,纯粹的k-means算法并不是非常灵活,同样地,在使用上有一定局限(不过上面说到得向量量化,确实是一个理想的应用场景)。特别是,当没有额外的限制条件时,参数k是很难选择的(真如上面讨论过的一样)。算法的另一个限制就是它不能和任意的距离函数一起使用、不能处理非数值数据。而正是为了满足这些使用条件,许多其他的算法才被发展起来。
3. 特征学习 :在(半)监督学习或无监督学习中，k-均值聚类被用来进行特征学习(或字典学习)步骤。基本方法是，首先使用输入数据训练出一个k-均值聚类表示,然后把任意的输入数据投射到这一新的特征空间。 k-均值的这一应用能成功地与自然语言处理和计算机视觉中半监督学习的简单线性分类器结合起来。在对象识别任务中，它能展现出与其他复杂特征学习方法（如自动编码器、受限Boltzmann机等）相当的效果。然而，相比复杂方法，它需要更多的数据来达到相同的效果，因为每个数据点都只贡献了一个特征（而不是多重特征）。
4. 与其他统计机器学习方法的关系 :k-means聚类以及它与EM算法的联系,是高斯混合模型的一个特例。很容易能把k-means问题一般化为高斯混合模型。另一个k-means算法的推广则是k-SVD算法,后者把数据点视为“编码本向量”的稀疏线性组合。而k-means对应于使用单编码本向量的特殊情形(其权重为1)。
5. Mean Shift聚类 :基本的Mean Shift聚类要维护一个与输入数据集规模大小相同的数据点集。初始时,这一集合就是输入集的副本。然后对于每一个点,用一定距离范围内的所有点的均值来迭代地替换它。与之对比,k-means把这样的迭代更新限制在(通常比输入数据集小得多的)K个点上,而更新这些点时,则利用了输入集中与之相近的所有点的均值(亦即在每个点的Voronoi划分内)。还有一种与k-means类似的Mean shift算法,即似然Mean shift,对于迭代变化的集合,用一定距离内在输入集中所有点的均值来更新集合里的点。Mean Shift聚类与k-means聚类相比,有一个优点就是不用指定聚类数目,因为Mean shift倾向于找到尽可能少的聚类数目。然而Mean shift会比k-means慢得多,并且同样需要选择一个"宽度"参数。和k-means一样,Mean shift算法有许多变体。
6. 主成分分析(PCA) :有一些研究表明,k-means的放松形式解(由聚类指示向量表示),可由主成分分析中的主成分给出,并且主成分分析由主方向张成的子空间与聚类图心空间是等价的。不过,主成分分析是k-means聚类的有效放松形式并不是一个新的结果,并且还有的研究结果直接揭示了关于“聚类图心子空间是由主成分方向张成的”这一论述的反例。
7. 独立成分分析(ICA) :有研究表明,在稀疏假设以及输入数据经过白化的预处理后,k-means得到的解就是独立成分分析的解。这一结果对于解释k-means在特征学习方面的成功应用很有帮助。
8. 双向过滤 :k-means算法隐含地假设输入数据的顺序不影响结果。双向过滤与k-means算法和Mean shift算法类似之处在于它同样维护着一个迭代更新的数据集(亦是被均值更新)。然而,双向过滤限制了均值的计算只包含了在输入数据中顺序相近的点,这使得双向过滤能够被应用在图像去噪等数据点的空间安排是非常重要的问题中。

• k-means++算法

k-means++算法在聚类中心的初始化过程中的基本原则是使得初始的聚类中心之间的相互距离尽可能远,修改因为初始化聚类中心过短的问题。k-means++算法的初始化过程如下所示:

• 在数据集中随机选择一个样本点作为第一个初始化的聚类中心$c_1$
• 选择出其余的聚类中心:
• 计算样本中的每一个样本点与已经初始化的聚类中心之间的距离,并选择其中最短的距离,记为$D(x)$
• 以概率$\frac{D(x)}{\sum_{x\in X}D(x)}$选择距离最大的样本作为新的聚类中心,重复上述过程,直到$k$个聚类中心都被确定
• 对k个初始化的聚类中心,利用k-means算法计算最终的聚类中心。

• ISODATA算法

正如之前所述,k-means和k-means++的聚类中心数k是固定不变的。而ISODATA算法在运行过程中能够根据各个类别的实际情况进行两种操作来调整聚类中心数k:分裂操作,对应着增加聚类中心数;合并操作,对应着减少聚类中心数。

下面首先给出ISODATA算法的输入(输入的数据和迭代次数不再单独介绍):

第一步 :输入$N$个模式样本$\{x_i,i=1,2,\dots,N\}$,预选$N_c$个初始聚类中心$\{z_1,z_2,\dots,z_{N_c}\}$,它可以不等于所要求的聚类中心的数目,其初始位置可以从样本中任意选取。预选参数:

$k_0$:预期的聚类中心数目;

$\theta_N$,每一聚类域中最少的样本数目,若少于此数即不作为一个独立的聚类;

$\theta_S$,一个聚类域中样本距离分布的标准差;

$\theta_c$:两个聚类中心间的最小距离,若小于此数,两个聚类需进行合并;

$L$:在一次迭代运算中可以合并的聚类中心的最多对数;

$I$:迭代运算的次数。

第二步 :将$N$个样本分给最近的聚类$S_j$,假若$D_j=\min\{\|x-z_i\|,i=1,2,\dots,N_c\}$,即$\|x-z_j\|$的距离最小,则$x\in S_j$。

第三步 :如果$S_j$中的样本数目$S_j<\theta_N$,则取消该样本子集,此时$N_c$减去1。

第四步 :修正各聚类中心

第五步 :计算各聚类域$S_j$中样本与各聚类中心间的平均距离

第六步 :计算全部模式样本和其对应聚类中心的总平均距离

第七步 :判别分裂、合并及迭代运算

若迭代运算次数已达到$I$次,即最后一次迭代,则置$\theta_c=0$,转至第十一步

若$N_c\le \frac{K}{2}$,即聚类中心的数目小于或等于规定值的一半,则转至第八步,对已有聚类进行分裂处理。

若迭代运算的次数是偶数次,或$N_c\ge 2K$,不进行分裂处理,转至第十一步;否则(即既不是偶数次迭代,又不满足$N_c\ge 2K$),转至第八步,进行分裂处理。

第八步 :计算每个聚类中样本距离的标准差向量$\sigma_j=(\sigma_1,\sigma_2,\dots,\sigma_{n_j})^T$,其中向量的各个分量为:

式中$i=1,2,\dots,n$为样本特征向量的维数,$j=1,2,\dots,N_c$为聚类数,$N_j$为$S_j$中的样本个数。

第九步 :求每一标准差向量$\{\sigma_j,j=1,2,\dots,N_c\}$中的最大分量,以$\{\sigma_{j_{max}},j=1,2,\dots,N_c\}$代表。

第十步:在任一最大分量集$\{\sigma_{j_{max}},j=1,2,\dots,N_c\}$中,若有$\sigma_{j_{max}}>\theta_S$,同时又满足如下两个条件之一:

• $\bar{D}_j>\bar{D}$和$N_j>2(\theta_N+1)$,即$S_j$中样本总数超过规定值一倍以上
• $N_c\le \frac{K}{2}$

则将$z_j$分裂为两个新的聚类中心且$N_c$加1。$z_j^+=z_j+\sigma_{j_{max}},z_j^-=z_j-\sigma_{j_{max}}$。如果本步骤完成了分裂运算,则转至第二步,否则继续。

第十一步 :计算全部聚类中心的距离:

第十二步 :比较$D_{ij}$与$\theta_c$的值,将$D_{ij}<\theta_c$的值按最小距离次序递增排列,即

式中$D_{i_1j_1}<D_{i_2j_2}<\cdots D_{i_Lj_L}$.

第十三步 :将距离为$D_{i_kj_k}$的两个聚类中心$Z_{ik}$和$Z_{jk}$合并,得新的中心为:

式中,被合并的两个聚类中心向量分别以其聚类域内的样本数加权,使$Z^*_k$为真正的平均向量。

第十四步 :如果是最后一次迭代运算(即第$I$次),则算法结束;否则,若需要操作者改变输入参数,转至第一步;若输入参数不变,转至第二步。在本步运算中,迭代运算的次数每次应加1。

• scikit-learn代码

class sklearn.cluster.KMeans(n_clusters=8, *, init='k-means++', n_init=10, max_iter=300, tol=0.0001, precompute_distances='deprecated', verbose=0, random_state=None, copy_x=True, n_jobs='deprecated', algorithm='auto')
    """
    Parameters
        n_clustersint, default=8
            The number of clusters to form as well as the number of centroids to generate.
        init{‘k-means++’, ‘random’, ndarray, callable}, default=’k-means++’
            Method for initialization:
            ‘k-means++’ : selects initial cluster centers for k-mean clustering in a smart way to speed up convergence. 
            ‘random’: choose n_clusters observations (rows) at random from data for the initial centroids.
            If an ndarray is passed, it should be of shape (n_clusters, n_features) and gives the initial centers.
            If a callable is passed, it should take arguments X, n_clusters and a random state and return an initialization.
        n_initint, default=10
            Number of time the k-means algorithm will be run with different centroid seeds. The final results will be the best output of n_init consecutive runs in terms of inertia.
        max_iterint, default=300
            Maximum number of iterations of the k-means algorithm for a single run.
        tolfloat, default=1e-4
            Relative tolerance with regards to Frobenius norm of the difference in the cluster centers of two consecutive iterations to declare convergence. It’s not advised to set tol=0 since convergence might never be declared due to rounding errors. Use a very small number instead.
        precompute_distances{‘auto’, True, False}, default=’auto’
            Precompute distances (faster but takes more memory).
            ‘auto’ : do not precompute distances if n_samples * n_clusters > 12 million. This corresponds to about 100MB overhead per job using double precision.
            True : always precompute distances.
            False : never precompute distances.
        verboseint, default=0
            Verbosity mode.
        random_stateint, RandomState instance, default=None
            Determines random number generation for centroid initialization. Use an int to make the randomness deterministic. 
        copy_xbool, default=True
            When pre-computing distances it is more numerically accurate to center the data first. If copy_x is True (default), then the original data is not modified. If False, the original data is modified, and put back before the function returns, but small numerical differences may be introduced by subtracting and then adding the data mean. Note that if the original data is not C-contiguous, a copy will be made even if copy_x is False. If the original data is sparse, but not in CSR format, a copy will be made even if copy_x is False.
        n_jobsint, default=None
            The number of OpenMP threads to use for the computation. Parallelism is sample-wise on the main cython loop which assigns each sample to its closest center.
            None or -1 means using all processors.
    Attributes
        cluster_centers_ndarray of shape (n_clusters, n_features)
            Coordinates of cluster centers. If the algorithm stops before fully converging (see tol and max_iter), these will not be consistent with labels_.
        labels_ndarray of shape (n_samples,)
            Labels of each point
        inertia_float
            Sum of squared distances of samples to their closest cluster center.
        n_iter_int
            Number of iterations run.
    """

参考使用样例:

>>> from sklearn.cluster import KMeans
>>> import numpy as np
>>> X = np.array([[1, 2], [1, 4], [1, 0],
...               [10, 2], [10, 4], [10, 0]])
>>> kmeans = KMeans(n_clusters=2, random_state=0).fit(X)
>>> kmeans.labels_
array([1, 1, 1, 0, 0, 0], dtype=int32)
>>> kmeans.predict([[0, 0], [12, 3]])
array([1, 0], dtype=int32)
>>> kmeans.cluster_centers_
array([[10.,  2.],
       [ 1.,  2.]])

# Vector Quantization Example
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

from sklearn import cluster


try:  # SciPy >= 0.16 have face in misc
    from scipy.misc import face
    face = face(gray=True)
except ImportError:
    face = sp.face(gray=True)

n_clusters = 5
np.random.seed(0)

X = face.reshape((-1, 1))  # We need an (n_sample, n_feature) array
k_means = cluster.KMeans(n_clusters=n_clusters, n_init=4)
k_means.fit(X)
values = k_means.cluster_centers_.squeeze()
labels = k_means.labels_

# create an array from labels and values
face_compressed = np.choose(labels, values)
face_compressed.shape = face.shape

vmin = face.min()
vmax = face.max()

# original face
plt.figure(1, figsize=(3, 2.2))
plt.imshow(face, cmap=plt.cm.gray, vmin=vmin, vmax=256)

# compressed face
plt.figure(2, figsize=(3, 2.2))
plt.imshow(face_compressed, cmap=plt.cm.gray, vmin=vmin, vmax=vmax)

# equal bins face
regular_values = np.linspace(0, 256, n_clusters + 1)
regular_labels = np.searchsorted(regular_values, face) - 1
regular_values = .5 * (regular_values[1:] + regular_values[:-1])  # mean
regular_face = np.choose(regular_labels.ravel(), regular_values, mode="clip")
regular_face.shape = face.shape
plt.figure(3, figsize=(3, 2.2))
plt.imshow(regular_face, cmap=plt.cm.gray, vmin=vmin, vmax=vmax)

# histogram
plt.figure(4, figsize=(3, 2.2))
plt.clf()
plt.axes([.01, .01, .98, .98])
plt.hist(X, bins=256, color='.5', edgecolor='.5')
plt.yticks(())
plt.xticks(regular_values)
values = np.sort(values)
for center_1, center_2 in zip(values[:-1], values[1:]):
    plt.axvline(.5 * (center_1 + center_2), color='b')

for center_1, center_2 in zip(regular_values[:-1], regular_values[1:]):
    plt.axvline(.5 * (center_1 + center_2), color='b', linestyle='--')

plt.show()

二、AP聚类

Affinity Propagation Clustering(简称AP算法)是2007提出的,当时发表在Science上《single-exemplar-based》。特别适合高维、多类数据快速聚类,相比传统的聚类算法,该算法算是比较新的,从聚类性能和效率方面都有大幅度的提升。

AP算法的基本思想:将全部样本看作网络的节点,然后通过网络中各条边的消息传递计算出各样本的聚类中心。聚类过程中,共有两种消息在各节点间传递,分别是吸引度(responsibility)和归属度(availability)。AP算法通过迭代过程不断更新每一个点的吸引度和归属度值,直到产生m个高质量的Exemplar(类似于质心),同时将其余的数据点分配到相应的聚类中。

名词介绍:

Exemplar:指的是聚类中心,K-Means中的质心,AP算法不需要事先指定聚类数目,相反它将所有的数据点都作为潜在的聚类中心。

Similarity(相似度):数据点i和点j的相似度记为s(i,j),是指点j作为点i的聚类中心的相似度。一般使用欧氏距离来计算,一般点与点的相似度值全部取为负值;因此,相似度值越大说明点与点的距离越近,便于后面的比较计算。

Preference:数据点i的参考度称为p(i)或s(i,i),是指点i作为聚类中心的参考度,以S矩阵的对角线上的数值s(k,k)作为k点能否成为聚类中心的评判标准,这意味着该值越大,这个点成为聚类中心的可能性也就越大。一般取s相似度值的中值(scikit-learn中默认为中位数)。聚类的数量受到参考度p的影响,如果认为每个数据点都有可能作为聚类中心,那么p就应取相同的值。如果取输入的相似度的均值作为p的值,得到聚类数量是中等的。如果取最小值,得到类数较少的聚类。

吸引度Responsibility:r(i,k)用来描述点k适合作为数据点i的聚类中心的程度。

归属度Availability:a(i,k)用来描述点i选择点k作为其聚类中心的适合程度。

Damping factor(阻尼系数):主要是起收敛作用的。

在实际计算应用中,最重要的两个参数(也是需要手动指定)是Preference和Damping factor。前者定了聚类数量的多少,值越大聚类数量越多;后者控制算法收敛效果。

AP算法流程:

步骤1:算法初始,将吸引度矩阵R和归属度矩阵初始化为0矩阵;

步骤2:更新吸引度矩阵

$$r_{t+1}(i,k)=\{_{S(i,k)-max_{j\ne k}{\{S(i,j)\}},i= k}^{S(i,k)-max_{j\ne k}{\{a_{t}(i,j)+r_{t}(i,j)\}},i\ne k}$$

步骤3:更新归属度矩阵

$$a_{t+1}(i,k)=\{_{\sum_{j\ne k}{max\{r_{t+1}(j,k),0\}},i=k}^{min\{0,r_{t+1}(k,k)+\sum_{j\ne i,k}{max\{r_{t+1}{(j,k)},0\}}\},i\ne k}$$

步骤4:根据衰减系数$\lambda$对两个公式进行衰减

$$r_{t+1}(i,k)=\lambda*r_{t}(i,k)+(1-\lambda)*r_{t+1}(i,k)\\ a_{t+1}(i,k)=\lambda*a_{t}(i,k)+(1-\lambda)*a_{t+1}(i,k)$$

重复步骤2,3,4直至矩阵稳定或者达到最大迭代次数,算法结束。最终取a+r最大的k作为聚类中心。

AP聚类算法的特点:

无需指定聚类"数量"参数。AP聚类不需要指定K(经典的K-Means)或者是其他描述聚类个数(SOM中的网络结构和规模)的参数,这使得先验经验成为应用的非必需条件,人群应用范围增加。

明确的质心(聚类中心点)。样本中的所有数据点都可能成为AP算法中的质心,叫做Examplar,而不是由多个数据点求平均而得到的聚类中心(如K-Means)。

对距离矩阵的对称性没要求。AP通过输入相似度矩阵来启动算法,因此允许数据呈非对称,数据适用范围非常大。

初始值不敏感。多次执行AP聚类算法,得到的结果是完全一样的,即不需要进行随机选取初值步骤(还是对比K-Means的随机初始值)。

算法复杂度较高,为$O(N^2T)$,$N$为样本数,$T$为迭代次数,而K-Means只是$O(N*K)$的复杂度。因此当N比较大时(N>3000),AP聚类算法往往需要算很久。

若以误差平方和来衡量算法间的优劣,AP聚类比其他方法的误差平方和都要低。(无论k-center clustering重复多少次,都达不到AP那么低的误差平方和)

AP算法相对K-Means鲁棒性强且准确度较高,但没有任何一个算法是完美的,AP聚类算法也不例外:

AP聚类应用中需要手动指定Preference和Damping factor,这其实是原有的聚类“数量”控制的变体。

算法较慢。由于AP算法复杂度较高,运行时间相对K-Means长,这会使得尤其在海量数据下运行时耗费的时间很多。

• scikit-learn代码

class sklearn.cluster.AffinityPropagation(*, damping=0.5, max_iter=200, convergence_iter=15, copy=True, preference=None, affinity='euclidean', verbose=False, random_state='warn')[source]
    """
    Perform Affinity Propagation Clustering of data.
    Parameters
        dampingfloat, default=0.5
            Damping factor (between 0.5 and 1) is the extent to which the current value is maintained relative to incoming values (weighted 1 - damping). This in order to avoid numerical oscillations when updating these values (messages).
        max_iterint, default=200
            Maximum number of iterations.
        convergence_iterint, default=15
            Number of iterations with no change in the number of estimated clusters that stops the convergence.
        copybool, default=True
            Make a copy of input data.
        preferencearray-like of shape (n_samples,) or float, default=None
            Preferences for each point - points with larger values of preferences are more likely to be chosen as exemplars. The number of exemplars, ie of clusters, is influenced by the input preferences value. If the preferences are not passed as arguments, they will be set to the median of the input similarities.
        affinity{‘euclidean’, ‘precomputed’}, default=’euclidean’
            Which affinity to use. At the moment ‘precomputed’ and euclidean are supported. ‘euclidean’ uses the negative squared euclidean distance between points.
        verbosebool, default=False
            Whether to be verbose.
        random_stateint or np.random.RandomStateInstance, default: 0
            Pseudo-random number generator to control the starting state. Use an int for reproducible results across function calls. See the Glossary.

    Attributes
        cluster_centers_indices_ndarray of shape (n_clusters,)  
            Indices of cluster centers
        cluster_centers_ndarray of shape (n_clusters, n_features)
            Cluster centers (if affinity != precomputed).
        labels_ndarray of shape (n_samples,)
            Labels of each point
        affinity_matrix_ndarray of shape (n_samples, n_samples)
            Stores the affinity matrix used in fit.
        n_iter_int
            Number of iterations taken to converge.
    """

方法	说明
fit(self, X[, y])	Fit the clustering from features, or affinity matrix.
fit_predict(self, X[, y])	Fit the clustering from features or affinity matrix, and return cluster labels.
get_params(self[, deep])	Get parameters for this estimator.
predict(self, X)	Predict the closest cluster each sample in X belongs to.
set_params(self, **params)	Set the parameters of this estimator.

样例

>>> from sklearn.cluster import AffinityPropagation
>>> import numpy as np
>>> X = np.array([[1, 2], [1, 4], [1, 0],
...               [4, 2], [4, 4], [4, 0]])
>>> clustering = AffinityPropagation(random_state=5).fit(X)
>>> clustering
AffinityPropagation(random_state=5)
>>> clustering.labels_
array([0, 0, 0, 1, 1, 1])
>>> clustering.predict([[0, 0], [4, 4]])
array([0, 1])
>>> clustering.cluster_centers_
array([[1, 2],
       [4, 2]])

# Demo of affinity propagation clustering algorithm
# #############################################################################
# Generate sample data
centers = [[1, 1], [-1, -1], [1, -1]]
X, labels_true = make_blobs(n_samples=300, centers=centers, cluster_std=0.5,
                            random_state=0)

# #############################################################################
# Compute Affinity Propagation
af = AffinityPropagation(preference=-50).fit(X)
cluster_centers_indices = af.cluster_centers_indices_
labels = af.labels_

n_clusters_ = len(cluster_centers_indices)

print('Estimated number of clusters: %d' % n_clusters_)
print("Homogeneity: %0.3f" % metrics.homogeneity_score(labels_true, labels))
print("Completeness: %0.3f" % metrics.completeness_score(labels_true, labels))
print("V-measure: %0.3f" % metrics.v_measure_score(labels_true, labels))
print("Adjusted Rand Index: %0.3f" % metrics.adjusted_rand_score(labels_true, labels))
print("Adjusted Mutual Information: %0.3f" % metrics.adjusted_mutual_info_score(labels_true, labels))
print("Silhouette Coefficient: %0.3f" % metrics.silhouette_score(X, labels, metric='sqeuclidean'))

# #############################################################################
# Plot result
import matplotlib.pyplot as plt
from itertools import cycle

plt.close('all')
plt.figure(1)
plt.clf()

colors = cycle('bgrcmykbgrcmykbgrcmykbgrcmyk')
for k, col in zip(range(n_clusters_), colors):
    class_members = labels == k
    cluster_center = X[cluster_centers_indices[k]]
    plt.plot(X[class_members, 0], X[class_members, 1], col + '.')
    plt.plot(cluster_center[0], cluster_center[1], 'o', markerfacecolor=col, markeredgecolor='k', markersize=14)
    for x in X[class_members]:
        plt.plot([cluster_center[0], x[0]], [cluster_center[1], x[1]], col)

plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()

https://sklearn.apachecn.org/docs/master/22.html#mini-batch-kmeans

2.2 Mean Shift

在K-Means算法中,最终的聚类效果受初始的聚类中心的影响,K-Means++算法的提出,为选择较好的初始聚类中心提供了依据,但是算法中,聚类的类别个数k仍需事先制定,对于类别个数事先未知的数据集,K-Means和K-Means++将很难对其精确求解,对此,有一些改进的算法被提出来处理聚类个数k未知的情形。Mean Shift算法,又被称为均值漂移算法,与K-Means算法一样,都是基于聚类中心的聚类算法,不同的是,Mean Shift算法不需要事先制定类别个数k。

Mean Shift的概念最早是由Fukunage在1975年提出的,在后来由Yizong Cheng对其进行扩充,主要提出了两点的改进:定义了核函数,增加了权重系数。核函数的定义使得偏移值对偏移向量的贡献随之样本与被偏移点的距离的不同而不同。权重系数使得不同样本的权重不同。

Mean Shift算法在很多领域都有成功应用,例如图像平滑、图像分割、物体跟踪等,这些属于人工智能里面模式识别或计算机视觉的部分;另外也包括常规的聚类应用。

图像平滑:图像最大质量下的像素压缩;
图像分割:跟图像平滑类似的应用,但最终是将可以平滑的图像进行分离已达到前后景或固定物理分割的目的;
目标跟踪:例如针对监控视频中某个人物的动态跟踪;
常规聚类,如用户聚类等。

• Mean Shift算法理论

Mean Shift向量

对于给定的d维空间$R^d$中的n个样本点$x_i,i=1,\cdots,n$,则对于$x$点,其Mean Shift向量的基本形式为:

其中,$S_h$指的是一个半径为h的高维球区域,如上图中的圆形区域。$S_h$的定义为:

里面所有点与圆心为起点形成的向量相加的结果就是Mean shift向量。下图黄色箭头就是$M_h$（Mean Shift向量）。

对于Mean Shift算法,是一个迭代的步骤,即先算出当前点的偏移均值,将该点移动到此偏移均值,然后以此为新的起始点,继续移动,直到满足最终的条件。

Mean Shift聚类就是对于集合中的每一个元素,对它执行下面的操作:把该元素移动到它邻域中所有元素的特征值的均值的位置,不断重复直到收敛。准确的说,不是真正移动元素,而是把该元素与它的收敛位置的元素标记为同一类。

如上的均值漂移向量的求解方法存在一个问题,即在$S_h$的区域内,每一个样本点$x$对样本$X$的共享是一样的。而实际中,每一个样本点$x$对样本$X$的贡献是不一样的,这样的共享可以通过核函数进行度量。

核函数

在Mean Shift算法中引入核函数的目的是使得随着样本与被偏移点的距离不同,其偏移量对均值偏移向量的贡献也不同。核函数是机器学习中常用的一种方式。核函数的定义如下所示:

$X$表示一个$d$维的欧式空间,$x$是该空间中的一个点$x = \{x_1,x_2,x_3\cdots ,x_d\}$,其中,x的模$\| x \|^2=xx^T$,R表示实数域,如果一个函数$K:X\rightarrow R$存在一个剖面函数$k:[0,\infty]\rightarrow R$,即

并且满足:

k是非负的
k是非增的
k是分段连续的

那么,函数$K(x)$就称为核函数。核函数有很多,下图中表示的每一个曲线都为一个核函数。

常用的核函数有高斯核函数。高斯核函数如下所示:

其中,$h$称为带宽(bandwidth),不同带宽的核函数如下图所示:

从高斯函数的图像可以看出,当带宽$h$一定时,样本点之间的距离越近,其核函数的值越大,当样本点之间的距离相等时,随着高斯函数的带宽$h$的增加,核函数的值在减小。

高斯核函数的python实现:

import numpy as np
import math

def gaussian_kernel(distance, bandwidth):
    """
    高斯核函数
        :param distance: 欧氏距离计算函数
        :param bandwidth: 核函数的带宽
        :return: 高斯函数值
    """
    m = np.shape(distance)[0]           # 样本个数
    right = np.mat(np.zeros((m, 1)))    # m*1矩阵
    for i in range(m):
        right[i, 0] = (-0.5 * distance[i] * distance[i].T) / (bandwidth * bandwidth)
        right[i, 0] = np.exp(right[i, 0])
    left = 1 / (bandwidth * math.sqrt(2 * math.pi))
    gaussian_val = left * right
    return gaussian_val

• 引入核函数的Mean Shift向量

假设在半径为$h$的范围$S_h$范围内,为了使得每一个样本点$x$对于样本$X$的共享不一样,向基本的Mean Shift向量形式中增加核函数,得到如下改进的Mean Shift向量形式:

其中,$G(\frac{x_i-x}{h_i})$为核函数。通常,可以取$S_h$为整个数据集范围。

计算$M_h$时考虑距离的影响,同时也可以认为在所有的样本点$X$中,重要性并不一样,因此对每个样本还引入一个权重系数。如此以来就可以把Mean Shift形式扩展为:

其中,$w(x_i)$是一个赋给采样点的权重。

• Mean Shift的代码实现

算法的python实现

import numpy as np
import math

MIN_DISTANCE = 0.00001      # 最小误差

def euclidean_dist(pointA, pointB):
    # 计算pointA和pointB之间的欧式距离
    total = (pointA - pointB) * (pointA - pointB).T
    return math.sqrt(total)

def gaussian_kernel(distance, bandwidth):
    """
    高斯核函数
        :param distance: 欧氏距离计算函数
        :param bandwidth: 核函数的带宽
        :return: 高斯函数值
    """
    m = np.shape(distance)[0]  # 样本个数
    right = np.mat(np.zeros((m, 1)))
    for i in range(m):
        right[i, 0] = (-0.5 * distance[i] * distance[i].T) / (bandwidth * bandwidth)
        right[i, 0] = np.exp(right[i, 0])
    left = 1 / (bandwidth * math.sqrt(2 * math.pi))
    gaussian_val = left * right
    return gaussian_val


def shift_point(point, points, kernel_bandwidth):
    """
    计算均值漂移点
        :param point: 需要计算的点
        :param points: 所有的样本点
        :param kernel_bandwidth: 核函数的带宽
        :return point_shifted:漂移后的点
    """
    points = np.mat(points)
    m = np.shape(points)[0]  # 样本个数
    # 计算距离
    point_distances = np.mat(np.zeros((m, 1)))
    for i in range(m):
        point_distances[i, 0] = euclidean_dist(point, points[i])

    # 计算高斯核
    point_weights = gaussian_kernel(point_distances, kernel_bandwidth)

    # 计算分母
    all = 0.0
    for i in range(m):
        all += point_weights[i, 0]

    # 均值偏移
    point_shifted = point_weights.T * points / all
    return point_shifted


def group_points(mean_shift_points):
    """
    计算所属的类别
        :param mean_shift_points:漂移向量
        :return: group_assignment:所属类别
    """
    group_assignment = []
    m, n = np.shape(mean_shift_points)
    index = 0
    index_dict = {}
    for i in range(m):
        item = []
        for j in range(n):
            item.append(str(("%5.2f" % mean_shift_points[i, j])))

        item_1 = "_".join(item)
        if item_1 not in index_dict:
            index_dict[item_1] = index
            index += 1

    for i in range(m):
        item = []
        for j in range(n):
            item.append(str(("%5.2f" % mean_shift_points[i, j])))

        item_1 = "_".join(item)
        group_assignment.append(index_dict[item_1])
    return group_assignment


def train_mean_shift(points, kernel_bandwidth=2):
    """
    训练Mean Shift模型
        :param points: 特征数据
        :param kernel_bandwidth: 核函数带宽
        :return:
            points:特征点
            mean_shift_points:均值漂移点
            group:类别
    """
    mean_shift_points = np.mat(points)
    max_min_dist = 1
    iteration = 0
    m = np.shape(mean_shift_points)[0]  # 样本的个数
    need_shift = [True] * m             # 标记是否需要漂移

    # 计算均值漂移向量
    while max_min_dist > MIN_DISTANCE:
        max_min_dist = 0
        iteration += 1
        print("iteration : " + str(iteration))
        for i in range(0, m):
            # 判断每一个样本点是否需要计算偏置均值
            if not need_shift[i]:
                continue
            p_new = mean_shift_points[i]
            p_new_start = p_new
            p_new = shift_point(p_new, points, kernel_bandwidth)  # 对样本点进行偏移
            dist = euclidean_dist(p_new, p_new_start)  # 计算该点与漂移后的点之间的距离

            if dist > max_min_dist:  # 记录是有点的最大距离
                max_min_dist = dist
            if dist < MIN_DISTANCE:  # 不需要移动
                need_shift[i] = False

            mean_shift_points[i] = p_new
    # 计算最终的group
    group = group_points(mean_shift_points)  # 计算所属的类别
    return np.mat(points), mean_shift_points, group

以上代码实现了基本的流程,但是执行效率很慢,正式使用时建议使用scikit-learn库中的MeanShift。scikit-learn MeanShift演示

import numpy as np
from sklearn.cluster import MeanShift, estimate_bandwidth

data = []
f = open("k_means_sample_data.txt", 'r')
for line in f:
    data.append([float(line.split(',')[0]), float(line.split(',')[1])])
data = np.array(data)

# 通过下列代码可自动检测bandwidth值
# 从data中随机选取1000个样本,计算每一对样本的距离,然后选取这些距离的0.2分位数作为返回值,
# 当n_samples很大时,这个函数的计算量是很大的。
bandwidth = estimate_bandwidth(data, quantile=0.2, n_samples=1000)
print(bandwidth)
# bin_seeding设置为True就不会把所有的点初始化为核心位置,从而加速算法
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(data)
labels = ms.labels_
cluster_centers = ms.cluster_centers_
# 计算类别个数
labels_unique = np.unique(labels)
n_clusters = len(labels_unique)
print("number of estimated clusters : %d" % n_clusters)

# 画图
import matplotlib.pyplot as plt
from itertools import cycle

plt.figure(1)
plt.clf()  # 清楚上面的旧图形

# cycle把一个序列无限重复下去
colors = cycle('bgrcmyk')
for k, color in zip(range(n_clusters), colors):
    # current_member表示标签为k的记为true 反之false
    current_member = labels == k
    cluster_center = cluster_centers[k]
    # 画点
    plt.plot(data[current_member, 0], data[current_member, 1], color + '.')
    #画圈
    plt.plot(cluster_center[0], cluster_center[1], 'o',
             markerfacecolor=color,  #圈内颜色
             markeredgecolor='k',  #圈边颜色
             markersize=14)  #圈大小
plt.title('Estimated number of clusters: %d' % n_clusters)
plt.show()

执行效果:

scikit-learn MeanShift源码:https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/cluster/_mean_shift.py

https://sklearn.apachecn.org/docs/master/24.html#25-%E5%88%86%E8%A7%A3%E6%88%90%E5%88%86%E4%B8%AD%E7%9A%84%E4%BF%A1%E5%8F%B7%EF%BC%88%E7%9F%A9%E9%98%B5%E5%88%86%E8%A7%A3%E9%97%AE%E9%A2%98%EF%BC%89