漫谈 Clustering (1): k-means（转）

时间 2019-11-30

标签漫谈 clustering means 繁體版

原文原文链接

转：http://blog.pluskid.org/?p=39python

很久没有写 blog 了，一来是 blog 下线一段时间，而租 DreamHost 的事情又一直没弄好；二来是没有太多时间，每天都跑去实验室。如今主要折腾 Machine Learning 相关的东西，由于不少东西都不懂，因此平时也找一些资料来看。按照我之前的更新速度的话，这么长时间不写 blog 确定是要被闷坏的，因此我也以为仍是不按期地整理一下本身了解到的东西，放在 blog 上，一来梳理老是有助于加深理解的，二来也算共享一下知识了。那么，仍是从 clustering 提及吧。算法

Clustering 中文翻译做“聚类”，简单地说就是把类似的东西分到一组，同 Classification (分类)不一样，对于一个 classifier ，一般须要你告诉它“这个东西被分为某某类”这样一些例子，理想状况下，一个 classifier 会从它获得的训练集中进行“学习”，从而具有对未知数据进行分类的能力，这种提供训练数据的过程一般叫作 supervised learning (监督学习)，而在聚类的时候，咱们并不关心某一类是什么，咱们须要实现的目标只是把类似的东西聚到一块儿，所以，一个聚类算法一般只须要知道如何计算类似度就能够开始工做了，所以 clustering 一般并不须要使用训练数据进行学习，这在 Machine Learning 中被称做 unsupervised learning (无监督学习)。dom

举一个简单的例子：如今有一群小学生，你要把他们分红几组，让组内的成员之间尽可能类似一些，而组之间则差异大一些。最后分出怎样的结果，就取决于你对于“类似”的定义了，好比，你决定男生和男生是类似的，女生和女生也是类似的，而男生和女生之间则差异很大”，这样，你其实是用一个可能取两个值“男”和“女”的离散变量来表明了原来的一个小学生，咱们一般把这样的变量叫作“特征”。实际上，在这种状况下，全部的小学生都被映射到了两个点的其中一个上，已经很天然地造成了两个组，不须要专门再作聚类了。另外一种多是使用“身高”这个特征。我在读小学候，每周五在操场开会训话的时候会按照你们住的地方的地域和距离远近来列队，这样结束以后就能够结队回家了。除了让事物映射到一个单独的特征以外，一种常见的作法是同时提取 N 种特征，将它们放在一块儿组成一个 N 维向量，从而获得一个从原始数据集合到 N 维向量空间的映射——你老是须要显式地或者隐式地完成这样一个过程，由于许多机器学习的算法都须要工做在一个向量空间中。机器学习

那么让咱们再回到 clustering 的问题上，暂且抛开原始数据是什么形式，假设咱们已经将其映射到了一个欧几里德空间上，为了方便展现，就使用二维空间吧，以下图所示：函数

从数据点的大体形状能够看出它们大体聚为三个 cluster ，其中两个紧凑一些，剩下那个松散一些。咱们的目的是为这些数据分组，以便能区分出属于不一样的簇的数据，若是按照分组给它们标上不一样的颜色，就是这个样子：学习

那么计算机要如何来完成这个任务呢？固然，计算机尚未高级到可以“经过形状大体看出来”，不过，对于这样的 N 维欧氏空间中的点进行聚类，有一个很是简单的经典算法，也就是本文标题中提到的 k-means 。在介绍 k-means 的具体步骤以前，让咱们先来看看它对于须要进行聚类的数据的一个基本假设吧：对于每个 cluster ，咱们能够选出一个中心点 (center) ，使得该 cluster 中的全部的点到该中心点的距离小于到其余 cluster 的中心的距离。虽然实际状况中获得的数据并不能保证老是知足这样的约束，但这一般已是咱们所能达到的最好的结果，而那些偏差一般是固有存在的或者问题自己的不可分性形成的。例以下图所示的两个高斯分布，从两个分布中随机地抽取一些数据点出来，混杂到一块儿，如今要让你将这些混杂在一块儿的数据点按照它们被生成的那个分布分开来：优化

因为这两个分布自己有很大一部分重叠在一块儿了，例如，对于数据点 2.5 来讲，它由两个分布产生的几率都是相等的，你所作的只能是一个猜想；稍微好一点的状况是 2 ，一般咱们会将它归类为左边的那个分布，由于几率大一些，然而此时它由右边的分布生成的几率仍然是比较大的，咱们仍然有不小的概率会猜错。而整个阴影部分是咱们所能达到的最小的猜错的几率，这来自于问题自己的不可分性，没法避免。所以，咱们将 k-means 所依赖的这个假设看做是合理的。spa

基于这样一个假设，咱们再来导出 k-means 所要优化的目标函数：设咱们一共有 N 个数据点须要分为 K 个 cluster ，k-means 要作的就是最小化.net

$\displaystyle J = \sum_{n=1}^N\sum_{k=1}^K r_{nk} \|x_n-\mu_k\|^2$

这个函数，其中 $r_{nk}$ 在数据点 n 被归类到 cluster k 的时候为 1 ，不然为 0 。直接寻找 $r_{nk}$ 和 $\mu_k$ 来最小化并不容易，不过咱们能够采起迭代的办法：先固定 $\mu_k$ ，选择最优的 $r_{nk}$ ，很容易看出，只要将数据点归类到离他最近的那个中心就能保证最小。下一步则固定 $r_{nk}$ ，再求最优的 $\mu_k$ 。将对 $\mu_k$ 求导并令导数等于零，很容易获得最小的时候 $\mu_k$ 应该知足：翻译

$\displaystyle \mu_k=\frac{\sum_n r_{nk}x_n}{\sum_n r_{nk}}$

亦即 $\mu_k$ 的值应当是全部 cluster k 中的数据点的平均值。因为每一次迭代都是取到的最小值，所以只会不断地减少（或者不变），而不会增长，这保证了 k-means 最终会到达一个极小值。虽然 k-means 并不能保证老是能获得全局最优解，可是对于这样的问题，像 k-means 这种复杂度的算法，这样的结果已是很不错的了。

下面咱们来总结一下 k-means 算法的具体步骤：

选定 K 个中心 $\mu_k$ 的初值。这个过程一般是针对具体的问题有一些启发式的选取方法，或者大多数状况下采用随机选取的办法。由于前面说过 k-means 并不能保证全局最优，而是否能收敛到全局最优解其实和初值的选取有很大的关系，因此有时候咱们会屡次选取初值跑 k-means ，并取其中最好的一次结果。
将每一个数据点归类到离它最近的那个中心点所表明的 cluster 中。
用公式 $\mu_k = \frac{1}{N_k}\sum_{j\in\text{cluster}_k}x_j$ 计算出每一个 cluster 的新的中心点。
重复第二步，一直到迭代了最大的步数或者先后的的值相差小于一个阈值为止。

按照这个步骤写一个 k-means 实现其实至关容易了，在 SciPy 或者 Matlab 中都已经包含了内置的 k-means 实现，不过为了看看 k-means 每次迭代的具体效果，咱们不妨本身来实现一下，代码以下（须要安装 SciPy 和 matplotlib）：

#!/usr/bin/python
 
from __future__ import with_statement
import cPickle as pickle
from matplotlib import pyplot
from numpy import zeros, array, tile
from scipy.linalg import norm
import numpy.matlib as ml
import random
 
def kmeans(X, k, observer=None, threshold=1e-15, maxiter=300):
    N = len(X)
    labels = zeros(N, dtype=int)
    centers = array(random.sample(X, k))
    iter = 0
 
    def calc_J():
        sum = 0
        for i in xrange(N):
            sum += norm(X[i]-centers[labels[i]])
        return sum
 
    def distmat(X, Y):
        n = len(X)
        m = len(Y)
        xx = ml.sum(X*X, axis=1)
        yy = ml.sum(Y*Y, axis=1)
        xy = ml.dot(X, Y.T)
 
        return tile(xx, (m, 1)).T+tile(yy, (n, 1)) - 2*xy
 
    Jprev = calc_J()
    while True:
        # notify the observer
        if observer is not None:
            observer(iter, labels, centers)
 
        # calculate distance from x to each center
        # distance_matrix is only available in scipy newer than 0.7
        # dist = distance_matrix(X, centers)
        dist = distmat(X, centers)
        # assign x to nearst center
        labels = dist.argmin(axis=1)
        # re-calculate each center
        for j in range(k):
            idx_j = (labels == j).nonzero()
            centers[j] = X[idx_j].mean(axis=0)
 
        J = calc_J()
        iter += 1
 
        if Jprev-J < threshold:
            break
        Jprev = J
        if iter >= maxiter:
            break
 
    # final notification
    if observer is not None:
        observer(iter, labels, centers)
 
if __name__ == '__main__':
    # load previously generated points
    with open('cluster.pkl') as inf:
        samples = pickle.load(inf)
    N = 0
    for smp in samples:
        N += len(smp[0])
    X = zeros((N, 2))
    idxfrm = 0
    for i in range(len(samples)):
        idxto = idxfrm + len(samples[i][0])
        X[idxfrm:idxto, 0] = samples[i][0]
        X[idxfrm:idxto, 1] = samples[i][1]
        idxfrm = idxto
 
    def observer(iter, labels, centers):
        print "iter %d." % iter
        colors = array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
        pyplot.plot(hold=False)  # clear previous plot
        pyplot.hold(True)
 
        # draw points
        data_colors=[colors[lbl] for lbl in labels]
        pyplot.scatter(X[:, 0], X[:, 1], c=data_colors, alpha=0.5)
        # draw centers
        pyplot.scatter(centers[:, 0], centers[:, 1], s=200, c=colors)
 
        pyplot.savefig('kmeans/iter_%02d.png' % iter, format='png')
 
    kmeans(X, 3, observer=observer)

代码有些长，不过由于用 Python 来作这个事情确实不如 Matlab 方便，实际的 k-means 的代码只是 41 到 47 行。首先 3 个中心点被随机初始化，全部的数据点都尚未进行聚类，默认所有都标记为红色，以下图所示：

而后进入第一次迭代：按照初始的中心点位置为每一个数据点着上颜色，这是代码中第 41 到 43 行所作的工做，而后 45 到 47 行从新计算 3 个中心点，结果以下图所示：

能够看到，因为初始的中心点是随机选的，这样得出来的结果并非很好，接下来是下一次迭代的结果：

能够看到大体形状已经出来了。再通过两次迭代以后，基本上就收敛了，最终结果以下：

不过正如前面所说的那样 k-means 也并非万能的，虽然许多时候都能收敛到一个比较好的结果，可是也有运气很差的时候会收敛到一个让人不满意的局部最优解，例如选用下面这几个初始中心点：

最终会收敛到这样的结果：

不得不认可这并非很好的结果。不过其实大多数状况下 k-means 给出的结果都仍是很使人满意的，算是一种简单高效应用普遍的 clustering 方法。

Update 2010.04.25: 不少人都问我要 cluster.pkl ，我干脆把它上传上来吧，实际上是很容易本身生成的，点击这里下载。