K-means与高斯混合模型

K-means

http://blog.pluskid.org/?p=17

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

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

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

cluster

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

cluster

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

gaussian

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

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

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

亦即 $\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 和 http://matplotlib.org/）：

#!/usr/bin/pythonfrom __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 randomdef 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 = 0def calc_J():sum = 0for i in xrange(N):sum += norm(X[i]-centers[labels[i]])return sumdef 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*xyJprev = calc_J()while True:# notify the observerif 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 centerlabels = dist.argmin(axis=1)# re-calculate each centerfor j in range(k):idx_j = (labels == j).nonzero()centers[j] = X[idx_j].mean(axis=0)J = calc_J()iter += 1if Jprev-J < threshold:breakJprev = Jif iter >= maxiter:break# final notificationif observer is not None:observer(iter, labels, centers)if __name__ == '__main__':# load previously generated pointswith open('cluster.pkl') as inf:samples = pickle.load(inf)N = 0for smp in samples:N += len(smp[0])X = zeros((N, 2))idxfrm = 0for 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 = idxtodef observer(iter, labels, centers):print "iter %d." % itercolors = array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])pyplot.plot(hold=False)  # clear previous plotpyplot.hold(True)# draw pointsdata_colors=[colors[lbl] for lbl in labels]pyplot.scatter(X[:, 0], X[:, 1], c=data_colors, alpha=0.5)# draw centerspyplot.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 个中心点被随机初始化，所有的数据点都还没有进行聚类，默认全部都标记为红色，如下图所示：

iter_00

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

iter_01

可以看到，由于初始的中心点是随机选的，这样得出来的结果并不是很好，接下来是下一次迭代的结果：

iter_02

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

iter_04

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

iter_00_bad

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

iter_03_bad

不得不承认这并不是很好的结果。不过其实大多数情况下 k-means 给出的结果都还是很令人满意的，算是一种简单高效应用广泛的 clustering 方法。

高斯混合模型http://blog.pluskid.org/?p=39

GMM 和 k-means 很像，不过 GMM 是学习出一些概率密度函数来（所以 GMM 除了用在 clustering 上之外，还经常被用于 density estimation ），简单地说，k-means 的结果是每个数据点被 assign 到其中某一个 cluster 了，而 GMM 则给出这些数据点被 assign 到每个 cluster 的概率，又称作 soft assignment 。

得出一个概率有很多好处，因为它的信息量比简单的一个结果要多，比如，我可以把这个概率转换为一个 score ，表示算法对自己得出的这个结果的把握。也许我可以对同一个任务，用多个方法得到结果，最后选取“把握”最大的那个结果；另一个很常见的方法是在诸如疾病诊断之类的场所，机器对于那些很容易分辨的情况（患病或者不患病的概率很高）可以自动区分，而对于那种很难分辨的情况，比如，49% 的概率患病，51% 的概率正常，如果仅仅简单地使用 50% 的阈值将患者诊断为“正常”的话，风险是非常大的，因此，在机器对自己的结果把握很小的情况下，会“拒绝发表评论”，而把这个任务留给有经验的医生去解决。

废话说了一堆，不过，在回到 GMM 之前，我们再稍微扯几句。我们知道，不管是机器还是人，学习的过程都可以看作是一种“归纳”的过程，在归纳的时候你需要有一些假设的前提条件，例如，当你被告知水里游的那个家伙是鱼之后，你使用“在同样的地方生活的是同一种东西”这类似的假设，归纳出“在水里游的都是鱼”这样一个结论。当然这个过程是完全“本能”的，如果不仔细去想，你也不会了解自己是怎样“认识鱼”的。另一个值得注意的地方是这样的假设并不总是完全正确的，甚至可以说总是会有这样那样的缺陷的，因此你有可能会把虾、龟、甚至是潜水员当做鱼。也许你觉得可以通过修改前提假设来解决这个问题，例如，基于“生活在同样的地方并且穿着同样衣服的是同一种东西”这个假设，你得出结论：在水里有并且身上长有鳞片的是鱼。可是这样还是有问题，因为有些没有长鳞片的鱼现在又被你排除在外了。

在这个问题上，机器学习面临着和人一样的问题，在机器学习中，一个学习算法也会有一个前提假设，这里被称作“归纳偏执 (bias)”（bias 这个英文词在机器学习和统计里还有其他许多的意思）。例如线性回归，目的是要找一个函数尽可能好地拟合给定的数据点，它的归纳偏执就是“满足要求的函数必须是线性函数”。一个没有归纳偏执的学习算法从某种意义上来说毫无用处，就像一个完全没有归纳能力的人一样，在第一次看到鱼的时候有人告诉他那是鱼，下次看到另一条鱼了，他并不知道那也是鱼，因为两条鱼总有一些地方不一样的，或者就算是同一条鱼，在河里不同的地方看到，或者只是看到的时间不一样，也会被他认为是不同的，因为他无法归纳，无法提取主要矛盾、忽略次要因素，只好要求所有的条件都完全一样──然而哲学家已经告诉过我们了：世界上不会有任何样东西是完全一样的，所以这个人即使是有无比强悍的记忆力，也绝学不到任何一点知识。

这个问题在机器学习中称作“过拟合 (Overfitting)”，例如前面的回归的问题，如果去掉“线性函数”这个归纳偏执

The Need for Biases in Learning Generalizations

，因为对于 N 个点，我们总是可以构造一个 N-1 次多项式函数，让它完美地穿过所有的这 N 个点，或者如果我用任何大于 N-1 次的多项式函数的话，我甚至可以构造出无穷多个满足条件的函数出来。如果假定特定领域里的问题所给定的数据个数总是有个上限的话，我可以取一个足够大的 N ，从而得到一个（或者无穷多个）“超级函数”，能够 fit 这个领域内所有的问题。然而这个（或者这无穷多个）“超级函数”有用吗？只要我们注意到学习的目的（通常）不是解释现有的事物，而是从中归纳出知识，并能应用到新的事物上，结果就显而易见了。

没有归纳偏执或者归纳偏执太宽泛会导致 Overfitting ，然而另一个极端──限制过大的归纳偏执也是有问题的：如果数据本身并不是线性的，强行用线性函数去做回归通常并不能得到好结果。难点正在于在这之间寻找一个平衡点。不过人在这里相对于（现在的）机器来说有一个很大的优势：人通常不会孤立地用某一个独立的系统和模型去处理问题，一个人每天都会从各个来源获取大量的信息，并且通过各种手段进行整合处理，归纳所得的所有知识最终得以统一地存储起来，并能有机地组合起来去解决特定的问题。这里的“有机”这个词很有意思，搞理论的人总能提出各种各样的模型，并且这些模型都有严格的理论基础保证能达到期望的目的，然而绝大多数模型都会有那么一些“参数”（例如 K-means 中的 k ），通常没有理论来说明参数取哪个值更好，而模型实际的效果却通常和参数是否取到最优值有很大的关系，我觉得，在这里“有机”不妨看作是所有模型的参数已经自动地取到了最优值。另外，虽然进展不大，但是人们也一直都期望在计算机领域也建立起一个统一的知识系统（例如语意网就是这样一个尝试）。https://en.wikipedia.org/wiki/Semantic_Web

废话终于说完了，回到 GMM 。按照我们前面的讨论，作为一个流行的算法，GMM 肯定有它自己的一个相当体面的归纳偏执了。其实它的假设非常简单，顾名思义，Gaussian Mixture Model ，就是假设数据服从 Mixture Gaussian Distribution ，换句话说，数据可以看作是从数个 Gaussian Distribution 中生成出来的。实际上，我们在 K-means 和 K-medoids 两篇文章中用到的那个例子就是由三个 Gaussian 分布从随机选取出来的。实际上，从中心极限定理可以看出，Gaussian 分布（也叫做正态 (Normal) 分布）这个假设其实是比较合理的，除此之外，Gaussian 分布在计算上也有一些很好的性质，所以，虽然我们可以用不同的分布来随意地构造 XX Mixture Model ，但是还是 GMM 最为流行。另外，Mixture Model 本身其实也是可以变得任意复杂的，通过增加 Model 的个数，我们可以任意地逼近任何连续的概率密分布。

每个 GMM 由个 Gaussian 分布组成，每个 Gaussian 称为一个“Component”，这些 Component 线性加成在一起就组成了 GMM 的概率密度函数：

由于在对数函数里面又有加和，我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们采取之前从 GMM 中随机选点的办法：分成两步，实际上也就类似于 K-means 的两步。

估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率）：对于每个数据来说，它由第个 Component 生成的概率为
其中 $N_k = \sum_{i=1}^N \gamma(i, k)$ ，并且 $\pi_k$ 也顺理成章地可以估计为。
重复迭代前面两步，直到似然函数的值收敛为止。

当然，上面给出的只是比较“直观”的解释，想看严格的推到过程的话，可以参考 Pattern Recognition and Machine Learning 这本书的第九章。有了实际的步骤，再实现起来就很简单了。Matlab 代码如下：

（Update 2012.07.03：如果你直接把下面的代码拿去运行了，碰到 covariance 矩阵 singular 的情况，可以参见这篇文章。）

function varargout = gmm(X, K_or_centroids)
% ============================================================
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
%  - X: N-by-D data matrix.
%  - K_OR_CENTROIDS: either K indicating the number of
%       components or a K-by-D matrix indicating the
%       choosing of the initial K centroids.
%
%  - PX: N-by-K matrix indicating the probability of each
%       component generating each point.
%  - MODEL: a structure containing the parameters for a GMM:
%       MODEL.Miu: a K-by-D matrix.
%       MODEL.Sigma: a D-by-D-by-K matrix.
%       MODEL.Pi: a 1-by-K vector.
% ============================================================threshold = 1e-15;[N, D] = size(X);if isscalar(K_or_centroids)K = K_or_centroids;% randomly pick centroidsrndp = randperm(N);centroids = X(rndp(1:K), :);elseK = size(K_or_centroids, 1);centroids = K_or_centroids;end% initial values[pMiu pPi pSigma] = init_params();Lprev = -inf;while truePx = calc_prob();% new value for pGammapGamma = Px .* repmat(pPi, N, 1);pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);% new value for parameters of each ComponentNk = sum(pGamma, 1);pMiu = diag(1./Nk) * pGamma' * X;pPi = Nk/N;for kk = 1:KXshift = X-repmat(pMiu(kk, :), N, 1);pSigma(:, :, kk) = (Xshift' * ...(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);end% check for convergenceL = sum(log(Px*pPi'));if L-Lprev < thresholdbreak;endLprev = L;endif nargout == 1varargout = {Px};elsemodel = [];model.Miu = pMiu;model.Sigma = pSigma;model.Pi = pPi;varargout = {Px, model};endfunction [pMiu pPi pSigma] = init_params()pMiu = centroids;pPi = zeros(1, K);pSigma = zeros(D, D, K);% hard assign x to each centroidsdistmat = repmat(sum(X.*X, 2), 1, K) + ...repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...2*X*pMiu';[dummy labels] = min(distmat, [], 2);for k=1:KXk = X(labels == k, :);pPi(k) = size(Xk, 1)/N;pSigma(:, :, k) = cov(Xk);endendfunction Px = calc_prob()Px = zeros(N, K);for k = 1:KXshift = X-repmat(pMiu(k, :), N, 1);inv_pSigma = inv(pSigma(:, :, k));tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));Px(:, k) = coef * exp(-0.5*tmp);endend
end

函数返回的 Px 是一个 $N\times K$ 的矩阵，对于每一个 x_i ，我们只要取该矩阵第行中最大的那个概率值所对应的那个 Component 为 x_i 所属的 cluster 就可以实现一个完整的聚类方法了。对于最开始的那个例子，GMM 给出的结果如下：

gmm

相对于之前 K-means 给出的结果，这里的结果更好一些，左下角的比较稀疏的那个 cluster 有一些点跑得比较远了。当然，因为这个问题原本就是完全有 Mixture Gaussian Distribution 生成的数据，GMM （如果能求得全局最优解的话）显然是可以对这个问题做到的最好的建模。

另外，从上面的分析中我们可以看到 GMM 和 K-means 的迭代求解法其实非常相似（都可以追溯到 EM 算法，下一次会详细介绍），因此也有和 K-means 同样的问题──并不能保证总是能取到全局最优，如果运气比较差，取到不好的初始值，就有可能得到很差的结果。对于 K-means 的情况，我们通常是重复一定次数然后取最好的结果，不过 GMM 每一次迭代的计算量比 K-means 要大许多，一个更流行的做法是先用 K-means （已经重复并取最优值了）得到一个粗略的结果，然后将其作为初值（只要将 K-means 所得的 centroids 传入 gmm 函数即可），再用 GMM 进行细致迭代。

如我们最开始所讨论的，GMM 所得的结果（Px）不仅仅是数据点的 label ，而包含了数据点标记为每个 label 的概率，很多时候这实际上是非常有用的信息。最后，需要指出的是，GMM 本身只是一个模型，我们这里给出的迭代的办法并不是唯一的求解方法。感兴趣的同学可以自行查找相关资料。