Free Mind ? 漫谈 Clustering (1): k

漫谈 Clustering (1): k-means

本文是“漫谈 Clustering 系列”中的第 1 篇，参见本系列的其他文章。

好久没有写 blog 了，一来是 blog 下线一段时间，而租 DreamHost 的事情又一直没弄好；二来是没有太多时间，天天都跑去实验室。现在主要折腾 Machine Learning 相关的东西，因为很多东西都不懂，所以平时也找一些资料来看。按照我以前的更新速度的话，这么长时间不写 blog 肯定是要被闷坏的，所以我也觉得还是不定期地整理一下自己了解到的东西，放在 blog 上，一来梳理总是有助于加深理解的，二来也算共享一下知识了。那么，还是从 clustering 说起吧。

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

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

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

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

那么计算机要如何来完成这个任务呢？当然，计算机还没有高级到能够“通过形状大致看出来”，不过，对于这样的 N 维欧氏空间中的点进行聚类，有一个非常简单的经典算法，也就是本文标题中提到的 k-means 。在介绍 k-means 的具体步骤之前，让我们先来看看它对于需要进行聚类的数据的一个基本假设吧：对于每一个 cluster ，我们可以选出一个中心点 (center) ，使得该 cluster 中的所有的点到该中心点的距离小于到其他 cluster 的中心的距离。虽然实际情况中得到的数据并不能保证总是满足这样的约束，但这通常已经是我们所能达到的最好的结果，而那些误差通常是固有存在的或者问题本身的不可分性造成的。例如下图所示的两个高斯分布，从两个分布中随机地抽取一些数据点出来，混杂到一起，现在要让你将这些混杂在一起的数据点按照它们被生成的那个分布分开来：

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

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

这个函数，其中 在数据点 n 被归类到 cluster k 的时候为 1 ，否则为 0 。直接寻找 和 来最小化 并不容易，不过我们可以采取迭代的办法：先固定 ，选择最优的 ，很容易看出，只要将数据点归类到离他最近的那个中心就能保证 最小。下一步则固定 ，再求最优的 。将 对 求导并令导数等于零，很容易得到 最小的时候 应该满足：

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

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

1. 选定 K 个中心 的初值。这个过程通常是针对具体的问题有一些启发式的选取方法，或者大多数情况下采用随机选取的办法。因为前面说过 k-means 并不能保证全局最优，而是否能收敛到全局最优解其实和初值的选取有很大的关系，所以有时候我们会多次选取初值跑 k-means ，并取其中最好的一次结果。
2. 将每个数据点归类到离它最近的那个中心点所代表的 cluster 中。
3. 用公式 计算出每个 cluster 的新的中心点。
4. 重复第二步，一直到迭代了最大的步数或者前后的 的值相差小于一个阈值为止。

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

1            2            3            4            5            6            7            8            9            10            11            12            13            14            15            16            17            18            19            20            21            22            23            24            25            26            27            28            29            30            31            32            33            34            35            36            37            38            39            40            41            42            43            44            45            46            47            48            49            50            51            52            53            54            55            56            57            58            59            60            61            62            63            64            65            66            67            68            69            70            71            72            73            74            75            76            77            78            79            80            81            82            83            84            85            86            87            88            89            90            91            

#!/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 方法。

Filed under:
Machine Learning by pluskid

13 Responses to “漫谈 Clustering (1): k-means”

1. winsty
   February 5th, 2009 at 9:10 pm | Reply

   如果用模拟退火类似的算法综合k-means 是不是可以走出局部极值

2. pluskid
   February 6th, 2009 at 5:15 am | Reply

   @winsty
   恩，不知道模拟退火算法这个东西，刚才查了一下，发现看起来好像很牛的样子：

   • 初值无关
   • 几乎以概率 1 收敛于全局最优解
   • 具有并行性

   不过大致看了一下它的求解过程，看起来和求解 PageRank 那里用的方法差不多，就是有一个概率能跳出局部最优，而不是死陷在那里。真的要用在这里的话，就是直接去求 J 的极值了，和 K-means 基本上无关了。

   不过 K-means 的真正目的实际上是进行聚类，也就是标 label ，求得 J 的最小值只是其中一个附加产品，就算用退火能求出全局最小的 J ，却只是得到了一个 bound 而已，要从这个 J 的数值推导出所对应的 cluster 形态还是没有办法的事情啊。

3. winsty
   February 6th, 2009 at 6:32 am | Reply

   @pluskid
   不是啊 在最小化J的同时也能够求出对应的参数rnk吧
   这里有个伪代码
   2.1.1步骤就是按照你这篇文章里提到的那个办法，本质没有区别
   只是用SA走出局部极值，避免发生你最后一个图片那样的问题

   Simulated Annealing (SA) Algorithm:
   1 初始化：系统初温T ，初始状态S0 ，马尔可夫链长L，终止条件AIM
   2 while (true)
   2.1 对于k=1..L, 执行2.1.1到2.1.4
   2.1.1 从当前解S,产生新解SN ，他们之间的差值为D.
   2.1.2 若 (D<0 或 满足概率 exp(-D/T)),则 S:=SN.
   2.1.3 若(当前解S<当前最优解 SB),则SB:=S.
   2.1.4 if (T趋于0 或 连续AIM次迭代物最优值), 则可近似认为SB为最优，转3.
   2.2 降温
   2.3 S=SB
   3 输出 SB

4. pluskid
   February 6th, 2009 at 6:53 am | Reply

   @winsty
   哦！我大致明白了，就是每次迭代的时候实际上是有一个概率是否接受新的解了。不过还是不能直接套到 K-means 里面去，因为 K-means 每次产生新的解的方式是固定的，而不是随机的，换句话说，初值确定之后，后面会得到什么样的结果就已经定了。要用在这里的话，还要设计一个产生新解的步骤，就是对应到你那个伪代码的 2.1.1 那一步。

5. winsty
   February 6th, 2009 at 8:42 am | Reply

   @pluskid
   嗯 是的-.-
   不过这个也应该不太困难

6. pluskid
   February 7th, 2009 at 12:54 am | Reply

   @winsty
   恩，回头好好研究下，这个好像和 Markov random walk 有关系。

   ps: blog 上的时间好像和中国时间相差了十几个小时啊，得好好设置一下……

7. winsty
   February 7th, 2009 at 1:16 am | Reply

   @pluskid
   期待后续连载……
   花了一上午看这些文章和相关的wiki链接
   爽死了

8. pluskid
   February 7th, 2009 at 4:00 pm | Reply

   @winsty
   赞一下看链接的人，难得我用心良苦呢。恩，后面的会陆续出来，不过也急不得，每写一篇都得花不少功夫呢。我也是一边学一边写啊。

9. rhythm
   February 7th, 2009 at 4:03 pm | Reply

   看到clustering，第一反应是服务器集群，结果发现完全不是……不过也是很有趣的话题。话说这篇文章里的代码、公式和图是不是用CodeColorer、Latex for WordPress以及gnuplot做的？

10. pluskid
    February 7th, 2009 at 4:09 pm | Reply

    @rhythm
    代码是用 wp-syntax 高亮的，推荐一下这个插件。公式是 LaTeXRender 这个插件吧，好像不太好找，你需要的话我可以拷贝给你。图多是用 Python 的 matplotlib 或者直接用 Matlab 画的。

11. windstorm
    March 8th, 2009 at 3:33 pm | Reply

    写得很清楚，赞一个。特别是图，下次我写考试笔记的时候就直接用了，呵呵。

12. lyslys34
    June 1st, 2009 at 4:38 pm | Reply

    请问使用matlab如何画图以及运行此实例呢？

13. pluskid
    June 1st, 2009 at 4:51 pm | Reply

    @lyslys34
    这里的代码是 Python 的，Matlab 里自带了一个 kmeans 函数可以用的。