漫谈 Clustering

时间 2019-11-30

标签漫谈 clustering 繁體版

原文原文链接

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

k-means
很久没有写 blog 了，一来是 blog 下线一段时间，而租 DreamHost 的事情又一直没弄好；二来是没有太多时间，每天都跑去实验室。如今主要折腾 Machine Learning 相关的东西，由于不少东西都不懂，因此平时也找一些资料来看。按照我之前的更新速度的话，这么长时间不写 blog 确定是要被闷坏的，因此我也以为仍是不按期地整理一下本身了解到的东西，放在 blog 上，一来梳理老是有助于加深理解的，二来也算共享一下知识了。那么，仍是从 clustering 提及吧。 node
Clustering 中文翻译做"聚类"，简单地说就是把类似的东西分到一组，同 Classification（分类）不一样，对于一个classifier ，一般须要你告诉它"这个东西被分为某某类"这样一些例子，理想状况下，一个 classifier 会从它获得的训练集中进行"学习"，从而具有对未知数据进行分类的能力，这种提供训练数据的过程一般叫作 supervised learning (监督学习)，而在聚类的时候，咱们并不关心某一类是什么，咱们须要实现的目标只是把类似的东西聚到一块儿，所以，一个聚类算法一般只须要知道如何计算类似度就能够开始工做了，所以 clustering 一般并不须要使用训练数据进行学习，这在 Machine Learning 中被称做 unsupervised learning (无监督学习)。 python
举一个简单的例子：如今有一群小学生，你要把他们分红几组，让组内的成员之间尽可能类似一些，而组之间则差异大一些。最后分出怎样的结果，就取决于你对于"类似"的定义了，好比，你决定男生和男生是类似的，女生和女生也是类似的，而男生和女生之间则差异很大，这样，你其实是用一个可能取两个值"男"和"女"的离散变量来表明了原来的一个小学生，咱们一般把这样的变量叫作"特征"。实际上，在这种状况下，全部的小学生都被映射到了两个点的其中一个上，已经很天然地造成了两个组，不须要专门再作聚类了。另外一种多是使用"身高"这个特征。我在读小学候，每周五在操场开会训话的时候会按照你们住的地方的地域和距离远近来列队，这样结束以后就能够结队回家了。除了让事物映射到一个单独的特征以外，一种常见的作法是同时提取 N 种特征，将它们放在一块儿组成一个 N 维向量，从而获得一个从原始数据集合到 N 维向量空间的映射——你老是须要显式地或者隐式地完成这样一个过程，由于许多机器学习的算法都须要工做在一个向量空间中。 git
那么让咱们再回到 clustering 的问题上，暂且抛开原始数据是什么形式，假设咱们已经将其映射到了一个欧几里德空间上，为了方便展现，就使用二维空间吧，以下图所示：程序员

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

那么计算机要如何来完成这个任务呢？固然，计算机尚未高级到可以"经过形状大体看出来"，不过，对于这样的 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 算法的具体步骤：
选定 K个中心的初值。这个过程一般是针对具体的问题有一些启发式的选取方法，或者大多数状况下采用随机选取的办法。由于前面说过 k-means 并不能保证全局最优，而是否能收敛到全局最优解其实和初值的选取有很大的关系，因此有时候咱们会屡次选取初值跑 k-means ，并取其中最好的一次结果。
将每一个数据点归类到离它最近的那个中心点所表明的 cluster 中。
用公式计算出每一个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 ，我干脆把它上传上来吧，实际上是很容易本身生成的，点击这里下载

k-medoids

上一次咱们了解了一个最基本的 clustering 办法 k-means ，此次要说的 k-medoids 算法，其实从名字上就能够看出来，和 k-means 确定是很是类似的。事实也确实如此，k-medoids 能够算是 k-means 的一个变种。

k-medoids 和 k-means 不同的地方在于中心点的选取，在 k-means 中，咱们将中心点取为当前 cluster 中全部数据点的平均值：

而且咱们已经证实在固定了各个数据点的 assignment 的状况下，这样选取的中心点可以把目标函数最小化。然而在 k-medoids 中，咱们将中心点的选取限制在当前 cluster 所包含的数据点的集合中。换句话说，在 k-medoids 算法中，咱们将从当前 cluster 中选取这样一个点——它到其余全部（当前 cluster 中的）点的距离之和最小——做为中心点。k-means 和 k-medoids 之间的差别就相似于一个数据样本的均值 (mean) 和中位数 (median) 之间的差别：前者的取值范围能够是连续空间中的任意值，然后者只能在给样本给定的那些点里面选。那么，这样作的好处是什么呢？

一个最直接的理由就是k-means 对数据的要求过高了，它使用欧氏距离描述数据点之间的差别 (dissimilarity) ，从而能够直接经过求均值来计算中心点。这要求数据点处在一个欧氏空间之中。

然而并非全部的数据都能知足这样的要求，对于数值类型的特征，好比身高，能够很天然地用这样的方式来处理，可是类别 (categorical) 类型的特征就不行了。举一个简单的例子，若是我如今要对犬进行聚类，而且但愿直接在全部犬组成的空间中进行，k-means 就无能为力了，由于欧氏距离在这里不能用了：一只Samoyed 减去一只 Rough Collie 而后在平方一下？天知道那是什么！再加上一只 German Shepherd Dog 而后求一下平均值？根本无法算，k-means 在这里步履维艰！

在k-medoids 中，咱们把原来的目标函数中的欧氏距离改成一个任意的dissimilarity measure 函数：

Samoyed Rough Collie

最多见的方式是构造一个 dissimilarity matrix 来表明，其中的元素表示第i只狗和第 j只狗之间的差别程度，例如，两只 Samoyed 之间的差别能够设为 0 ，一只 German Shepherd Dog 和一只 Rough Collie 之间的差别是 0.7，和一只 Miniature Schnauzer 之间的差别是 1 ，等等。

除此以外，因为中心点是在已有的数据点里面选取的，所以相对于 k-means 来讲，不容易受到那些因为偏差之类的缘由产生的 Outlier 的影响，更加 robust 一些。

扯了这么多，仍是直接来看看 k-medoids 的效果好了，因为 k-medoids 对数据的要求比 k-means 要低，因此 k-means 能处理的状况天然 k-medoids 也能处理，为了能先睹为快，咱们偷一下懒，直接在上一篇文章中的 k-means 代码的基础上稍做一点修改，还用一样的例子。将代码的 45 到 47 行改为下面这样：

 for j in range(k):

 idx_j = (labels == j).nonzero()

 distj = distmat(X[idx_j], X[idx_j])

 distsum = ml.sum(distj, axis=1)

 icenter = distsum.argmin()

 centers[j] = X[idx_j[0][icenter]]

能够看到 k-medoids 在这个例子上也能获得很好的结果：

并且，同 k-means 同样，运气很差的时候也会陷入局部最优解中：

若是仔细看上面那段代码的话，就会发现，从 k-means 变到 k-medoids ，时间复杂度陡然增长了许多：在 k-means 中只要求一个平均值便可，而在 k-medoids 中则须要枚举每一个点，并求出它到全部其余点的距离之和，复杂度为。

看完了直观的例子，让咱们再来看一个稍微实际一点的例子好了：Document Clustering ——那个永恒不变的主题，不过咱们这里要作的聚类并非针对文档的主题，而是针对文档的语言。实验数据是从 Europarl 下载的包含 Danish、German、Greek、English、Spanish、Finnish、French、Italian、Dutch、Portuguese 和 Swedish 这些语言的文本数据集。

在 N-gram-based text categorization 这篇 paper 中描述了一种计算由不一样语言写成的文档的类似度的方法。一个（以字符为单位的） N-gram 就至关于长度为 N 的一系列连续子串。例如，由 hello 产生的 3-gram 为：hel、ell 和 llo ，有时候还会在划分 N-gram 以前在开头和末尾加上空格（这里用下划线表示）：_he、hel、ell、llo、lo_ 和 o__ 。按照 Zipf's law ：

The nth most common word in a human language text occurs with a frequency inversely proportional to n.

这里咱们用 N-gram 来代替 word 。这样，咱们从一个文档中能够获得一个 N-gram 的频率分布，按照频率排序一下，只保留频率最高的前 k 个（好比，300）N-gram，咱们把叫作一个"Profile"。正常状况下，某一种语言（至少是西方国家的那些类英语的语言）写成的文档，不论主题或长短，一般得出来的 Profile 都差很少，亦即按照出现的频率排序所获得的各个 N-gram 的序号不会变化太大。这是很是好的一个性质：一般咱们只要各个语言选取一篇（比较正常的，也不须要很长）文档构建出一个 Profile ，在拿到一篇未知文档的时候，只要和各个 Profile 比较一下，差别最小的那个 Profile 所对应的语言就能够认定是这篇未知文档的语言了——准确率很高，更难得的是，所须要的训练数据很是少并且容易得到，训练出来的模型也是很是小的。

不过，咱们这里且撇开分类（Classification）的问题，回到聚类（Clustering）上，按照前面的说法，在 k-medoids 聚类中，只须要定义好两个东西之间的距离（或者 dissimilarity ）就能够了，对于两个 Profile ，它们之间的 dissimilarity 能够很天然地定义为对应的 N-gram 的序号之差的绝对值，在 Python 中用下面这样一个类来表示：

class Profile(object):

 def __init__(self, path, N=3, psize=400):

 self.N = N

 self.psize = psize

 self.build_profile(path)

 sep = re.compile(r'\W+')

 def build_profile(self, path):

 grams = {}

 with open(path) as inf:

 for line in inf:

 for tok in self.sep.split(line):

 for n in range(self.N):

 self.feed_ngram(grams, tok, n+1)

 self.create_profile(grams.items())

 def create_profile(self, grams):

 # keep only the top most psize items

 grams.sort(key=itemgetter(1), reverse=True)

 grams = grams[:self.psize]

 self.profile = dict()

 for i in range(len(grams)):

 self.profile[grams[i][0]] = i

 def __getitem__(self, key):

 idx = self.profile.get(key)

 if idx is None:

 return len(self.profile)

 return idx

 def dissimilarity(self, o):

 dis = 0

 for tok in self.profile.keys():

 dis += abs(self[tok]-o[tok])

 for tok in o.profile.keys():

 dis += abs(self[tok]-o[tok])

 return dis

 def feed_ngram(self, grams, tok, n):

 if n != 0:

 tok = '_' + tok

 tok = tok + '_' * (n-1)

 for i in range(len(tok)-n+1):

 gram = tok[i:i+n]

 grams.setdefault(gram, 0)

 grams[gram] += 1

europarl 数据集共有 11 种语言的文档，每种语言包括大约 600 多个文档。我为这七千多个文档创建了 Profile 并构造出一个 7038×7038 的 dissimilarity matrix ，最后在这上面用 k-medoids 进行聚类。构造 dissimilarity matrix 的过程很慢，在我这里花了将近 10 个小时。相比之下，k-medoids 的过程在内存容许的状况下，采用向量化的方法来作实际上仍是很快的，而且一般只要数次迭代就能收敛了。实际的 k-medoids 实现能够在 mltk 中找到，从此若是有时间的话，我会陆续地把一些相关的比较通用的代码放到那里面。

最后，然咱们来看看聚类的结果，因为聚类和分类不一样，只是获得一些 cluster ，而并不知道这些 cluster 应该被打上什么标签，或者说，在这个问题里，包含了哪一种语言的文档。因为咱们的目的是衡量聚类算法的 performance ，所以直接假定这一步能实现最优的对应关系，将每一个 cluster 对应到一种语言上去。一种办法是枚举全部可能的状况并选出最优解，另外，对于这样的问题，咱们还能够用 Hungarian algorithm 来求解。

咱们这里有 11 种语言，全排列有 11! = 39916800 种状况，对于每一种排列，咱们须要遍历一次 label list ，并数出真正的 label （语言）与聚类得出的结果相同的文档的个数，再除以总的文档个数，获得 accuracy 。假设每次遍历并求出 accuracy 只须要 1 毫秒的时间的话，总共也须要 11 个小时才能获得结果。看上去好像也不是特别恐怖，不过相比起来，用 Hungarian algorithm 的话，咱们能够几乎瞬间获得结果。因为文章的篇幅已经很长了，就不在这里介绍具体的算法了，感兴趣的同窗能够参考 Wikipedia ，这里我直接使用了一个现有的 Python 实现。

虽然这个实验很是折腾，不过最后的结果其实就是一个数字：accuracy ——在我这里达到了 88.97% ，证实 k-medoids 聚类和 N-gram Profile 识别语言这两种方法都是挺不错的。最后，若是有感兴趣的同窗，代码能够从这里下载。须要最新版的 scipy， munkres.py 和 mltk 以及 Python 2.6 。

(番外篇): Vector Quantization

在接下去说其余的聚类算法以前，让咱们先插进来讲一说一个有点跑题的东西：Vector Quantization 。这项技术普遍地用在信号处理以及数据压缩等领域。事实上，在 JPEG 和 MPEG-4 等多媒体压缩格式里都有 VQ 这一步。

Vector Quantization 这个名字听起来有些玄乎，其实它自己并无这么高深。你们都知道，模拟信号是连续的值，而计算机只能处理离散的数字信号，在将模拟信号转换为数字信号的时候，咱们能够用区间内的某一个值去代替着一个区间，好比，[0, 1) 上的全部值变为 0 ，[1, 2) 上的全部值变成 1 ，如此类推。其这就是一个 VQ 的过程。一个比较正式一点的定义是：VQ 是将一个向量空间中的点用其中的一个有限子集来进行编码的过程。

一个典型的例子就是图像的编码。最简单的状况，考虑一个灰度图片，0 为黑色，1 为白色，每一个像素的值为 [0, 1] 上的一个实数。如今要把它编码为 256 阶的灰阶图片，一个最简单的作法就是将每个像素值 x 映射为一个整数 floor(x*255) 。固然，原始的数据空间也并不以必定要是连续的。好比，你如今想要把压缩这个图片，每一个像素只使用 4 bit （而不是原来的 8 bit）来存储，所以，要将原来的 [0, 255] 区间上的整数值用 [0, 15] 上的整数值来进行编码，一个简单的映射方案是 x*15/255 。

Rechard Stallman VQ2

VQ100 VQ10

不过这样的映射方案很有些 Naive ，虽然能减小颜色数量起到压缩的效果，可是若是原来的颜色并非均匀分布的，那么的出来的图片质量可能并非很好。例如，若是一个 256 阶灰阶图片彻底由 0 和 13 两种颜色组成，那么经过上面的映射就会获得一个全黑的图片，由于两个颜色全都被映射到 0 了。一个更好的作法是结合聚类来选取表明性的点。

实际作法就是：将每一个像素点看成一个数据，跑一下 K-means ，获得 k 个 centroids ，而后用这些 centroids 的像素值来代替对应的 cluster 里的全部点的像素值。对于彩色图片来讲，也能够用一样的方法来作，例如 RGB 三色的图片，每个像素被看成是一个 3 维向量空间中的点。

用本文开头那张 Rechard Stallman 大神的照片来作一下实验好了，VQ 二、VQ 10 和 VQ 100 三张图片分别显示聚类数目为 2 、10 和 100 时获得的结果，能够看到 VQ 100 已经和原图很是接近了。把原来的许多颜色值用 centroids 代替以后，总的颜色数量减小了，重复的颜色增长了，这种冗余正是压缩算法最喜欢的。考虑一种最简单的压缩办法：单独存储（好比 100 个）centroids 的颜色信息，而后每一个像素点存储 centroid 的索引而不是颜色信息值，若是一个 RGB 颜色值须要 24 bits 来存放的话，每一个（128 之内的）索引值只须要 7 bits 来存放，这样就起到了压缩的效果。

实现代码很简单，直接使用了 SciPy 提供的 kmeans 和 vq 函数，图像读写用了 Python Image Library ：

#!/usr/bin/python

from scipy.cluster.vq import kmeans, vq

from numpy import array, reshape, zeros

from mltk import image

vqclst = [2, 10, 100, 256]

data = image.read('example.jpg')

(height, width, channel) = data.shape

data = reshape(data, (height*width, channel))

for k in vqclst:

 print 'Generating vq-%d...' % k

 (centroids, distor) = kmeans(data, k)

 (code, distor) = vq(data, centroids)

 print 'distor: %.6f' % distor.sum()

 im_vq = centroids[code, :]

 image.write('result-%d.jpg' % k, reshape(im_vq,

 (height, width, channel)))

固然，Vector Quantization 并不必定要用 K-means 来作，各类能用的聚类方法均可以用，只是 K-means 一般是最简单的，并且一般都够用了。

Gaussian Mixture Model

上一次咱们谈到了用 k-means 进行聚类的方法，此次咱们来讲一下另外一个很流行的算法：Gaussian Mixture Model (GMM)。事实上，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)"，例如前面的回归的问题，若是去掉"线性函数"这个概括偏执，由于对于 N 个点，咱们老是能够构造一个 N-1 次多项式函数，让它完美地穿过全部的这 N 个点，或者若是我用任何大于 N-1 次的多项式函数的话，我甚至能够构造出无穷多个知足条件的函数出来。若是假定特定领域里的问题所给定的数据个数老是有个上限的话，我能够取一个足够大的 N ，从而获得一个（或者无穷多个）"超级函数"，可以 fit 这个领域内全部的问题。然而这个（或者这无穷多个）"超级函数"有用吗？只要咱们注意到学习的目的（一般）不是解释现有的事物，而是从中概括出知识，并能应用到新的事物上，结果就显而易见了。

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

废话终于说完了，回到 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个 Component 之中选一个，每一个 Component 被选中的几率实际上就是它的系数，选中了 Component 以后，再单独地考虑从这个 Component 的分布中选取一个点就能够了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

那么如何用 GMM 来作 clustering 呢？其实很简单，如今咱们有了数据，假定它们是由 GMM 生成出来的，那么咱们只要根据数据推出 GMM 的几率分布来就能够了，而后 GMM 的 K个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算几率密度一般被称做 density estimation ，特别地，当咱们在已知（或假定）了几率密度函数的形式，而要估计其中的参数的过程被称做"参数估计"。

如今假设咱们有 N 个数据点，并假设它们服从某个分布（记做），如今要肯定里面的一些参数的值，例如，在 GMM 中，咱们就须要肯定、和这些参数。咱们的想法是，找到这样一组参数，它所肯定的几率分布生成这些给定的数据点的几率最大，而这个几率实际上就等于，咱们把这个乘积称做似然函数 (Likelihood Function)。一般单个点的几率都很小，许多很小的数字相乘起来在计算机里很容易形成浮点数下溢，所以咱们一般会对其取对数，把乘积变成加和，获得 log-likelihood function 。接下来咱们只要将这个函数最大化（一般的作法是求导并令导数等于零，而后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，咱们就认为这是最合适的参数，这样就完成了参数估计的过程。

下面让咱们来看一看 GMM 的 log-likelihood function：

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

估计数据由每一个 Component 生成的几率（并非每一个 Component 被选中的几率）：对于每一个数据来讲，它由第 k 个 Component 生成的几率为

因为式子里的和也是须要咱们估计的值，咱们采用迭代法，在计算的时候咱们假定和均已知，咱们将取上一次迭代所得的值（或者初始值）。

估计每一个 Component 的参数：如今咱们假设上一步中获得的就是正确的"数据由 Component k 生成的概率"，亦能够当作该 Component 在生成这个数据上所作的贡献，或者说，咱们能够看做这个值其中有这部分是由 Component k 所生成的。集中考虑全部的数据点，如今实际上能够看做 Component 生成了这些点。因为每一个 Component 都是一个标准的 Gaussian 分布，能够很容易分布求出最大似然所对应的参数值：

其中，而且也瓜熟蒂落地能够估计为。

重复迭代前面两步，直到似然函数的值收敛为止。

固然，上面给出的只是比较"直观"的解释，想看严格的推到过程的话，能够参考 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 centroids

 rndp = randperm(N);

 centroids = X(rndp(1:K), :);

 else

 K = size(K_or_centroids, 1);

 centroids = K_or_centroids;

end

 % initial values

 [pMiu pPi pSigma] = init_params();

 Lprev = -inf;

 while true

 Px = calc_prob();

 % new value for pGamma

 pGamma = Px .* repmat(pPi, N, 1);

 pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);

 % new value for parameters of each Component

 Nk = sum(pGamma, 1);

 pMiu = diag(1./Nk) * pGamma' * X;

 pPi = Nk/N;

 for kk = 1:K

 Xshift = X-repmat(pMiu(kk, :), N, 1);

 pSigma(:, :, kk) = (Xshift' * ...

 (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);

end

 % check for convergence

 L = sum(log(Px*pPi'));

 if L-Lprev < threshold

 break;

end

 Lprev = L;

end

 if nargout == 1

 varargout = {Px};

 else

 model = [];

 model.Miu = pMiu;

 model.Sigma = pSigma;

 model.Pi = pPi;

 varargout = {Px, model};

end

 function [pMiu pPi pSigma] = init_params()

 pMiu = centroids;

 pPi = zeros(1, K);

 pSigma = zeros(D, D, K);

 % hard assign x to each centroids

 distmat = 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:K

 Xk = X(labels == k, :);

 pPi(k) = size(Xk, 1)/N;

 pSigma(:, :, k) = cov(Xk);

end

end

 function Px = calc_prob()

 Px = zeros(N, K);

 for k = 1:K

 Xshift = 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);

end

end

end

函数返回的 Px 是一个的矩阵，对于每个，咱们只要取该矩阵第 i 行中最大的那个几率值所对应的那个 Component 为所属的 cluster 就能够实现一个完整的聚类方法了。对于最开始的那个例子，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 自己只是一个模型，咱们这里给出的迭代的办法并非惟一的求解方法。感兴趣的同窗能够自行查找相关资料。

Regularized Gaussian Covariance Estimation
我以前写过一篇介绍 Gaussian Mixture Model (GMM) 的文章，并在文章里贴了一段 GMM实现的 Matlab 示例代码，而后就不断地有人来问我关于那段代码的问题，问得最多的就是你们常常发如今跑那段代码的时候估计出来的 Covariance Matrix 是 singular 的，因此在第 96 行求逆的时候会挂掉。这是今天要介绍的主要话题，我会讲得罗嗦一点，把关于那篇文章我被问到的一些其余问题也都提到一下，不过，在步入正题以前，不妨先来小小地闲扯一下。
我本身之前就很是在痛恨书里看到伪代码，又不能编译运行，还搞得像模像样的像代码通常，因此我本身在写 blog 的时候就尝试给出直接可运行的代码，可是如今我渐渐理解为何书的做者喜欢用伪代码而不是可运行的代码了（除非是专门讲某编程语言的书）。缘由很简单：示例用的代码和真实项目中的代码每每是差异很大的，直接给出可运行的示例代码，读者就会直接拿去运行（由于包括我在内的大部分人都是很懒的，不是说"懒"是程序员的天性么？），而每每示例代码为告终构清晰并突出算法的主要脉络，对不少琐碎状况都没有处理，都使用最直接的操做，没有作任何优化（例如个人那个 GMM 示例代码里直接用 matlab 的inv 函数来求 Covariance 矩阵的逆），等等，所以直接运行这样的代码——特别是在实际项目中使用，是很是不明智的。
但是那又为何不直接给出实际项目级别的代码呢？第一个缘由固然是工做量，我想程序员都知道从一个算法到一个健壮、高效的系统之间有多长的路要走，并且不少时候都须要根据不一样的项目环境和须要作不一样的修改。因此当一我的其实是在介绍算法的时候让他去弄一套工业级别的代码出来做为示例，实在是太费力了。固然更重要的缘由是：工业级别的代码一般通过大量的优化并充斥着大量错误处理、以及为了处理其余地方的 bug 或者整个系统混乱的 API 等而作的各类 workaround 等等……也许根本就看不懂，基本彻底失去了示例的做用。因此，伪代码就成为了惟一的选择。
罗嗦了半天，如今咱们回到正题，那篇文章讲的是 GMM 的学习，因此关于 GMM 的问题能够直接参考那篇文章，出问题的地方仅在于 GMM 的每一个 component (Gaussian Distribution) 的参数估计上，特别地，在于 Covariance Matrix 的估计上。因此，咱们今天主要来探讨 Gaussian 分布的协方差矩阵的估计问题。首先来看一下多维 Gaussian 分布的几率密度函数：
eq: 1 »
参数估计的过程是给定一堆数据点，并假设他们是由某个真实的 Gaussian 分布独立地 (IID) 采样出来的，而后根据这堆点来估计该 Gaussian 分布的参数。对于一个多维 Gaussian 分布来讲，须要估计的就是均值  和协方差矩阵  两"个"参数——用引号是由于和  都不是标量，假设咱们的数据是  维的，那么  其实是  个参数，而做为一个的矩阵，则是由  个参数组成。根据协方差矩阵的对称性，咱们能够进一步把  所包含的参数个数下降为  个。
值得一提的是，咱们从参数估计的角度并不能一会儿判定  就是对称的，须要稍加说明（这其实是《Pattern Recognition and Machine Learning》书里的习题 2.17）。记，假设咱们的  不是对称的，那么咱们把他写成一个对称矩阵和一个反对称阵之和（任何一个矩阵均可以这样表示）：
将起代入式 (eq: 1) 的指数部分中，注意到对于反对称部分（书写方便起见咱们暂时另）：
因此指数部分只剩下对称阵  了。也就是说，若是咱们找出一个符合条件的非对称阵，用它的对称部分  代替也是能够的，因此咱们就能够"不失通常性"地假设  是对称的（由于对称阵的逆矩阵也是对称的）。注意咱们不用考虑式 (eq: 1) 指数前面的系数（里的那个  的行列式），由于前面的系数是起 normalization 做用的，使得全几率积分等于 1，指数部分变更以后只要从新 normalize 一下便可（能够想像一下，若是  真的是非对称的，那么 normalization 系数那里应该也不会是  那么简单的形式吧？有兴趣的同窗能够本身研究一下）。
好，如今回到参数估计上。对于给定的数据，将其代入式 (eq: 1) 中能够获得其对应的 likelihood ，注意这个并不是  的几率。和离散几率分布不一样的是，这里谈论单个数据点的几率是没有多大意义的，由于单点集的测度为零，它的几率也老是为零。因此，这里的 likelihood 出现大于 1 的值也是很正常的——这也是我常常被问到的一个问题。虽然不是几率值，可是我想从 likelihood 或者几率密度的角度去理解，应该也是能够获得直观的认识的。
参数估计经常使用的方法是最大似然估计（Maximum Likelihood Estimation, MLE）。就是将全部给定数据点的似然所有乘起来，获得一个关于参数（这里是  和）的函数，而后最大化该函数，所对应的参数值就是最大似然估计的值。一般为了不数值下溢，会使用单调函数  将相乘变成相加在计算。
对于最大似然估计，个人理解是寻找最可能生成给定的数据的模型参数。这对于离散型的分布来讲是比较好解释的，由于此时一点的似然实际上就是该点的几率，不过对于像 Gaussian 这样的针对连续数据的分布来讲，就不是那么天然了。不过咱们其实有另外一种解释，这是我在Coursera 的 Probabilistic Graphical Model 课上看到的：
对于一个真实的分布，假设咱们获得了它的一个估计，那么如何来衡量这个估计的好坏呢？比较两个几率分布的一个经常使用的办法是用 KL-Divergence，又叫 Relative Entropy（以及其余若干外号），定义为：
不过因为咱们一般是不知道真实的  的，因此这个值也无法直接算，因而咱们把它简单地变形一下：
其中蓝色部分是分布  本身的 Entropy，因为咱们比较不一样的模型估计好坏的时候这个量是不变的，因此咱们能够忽略它，只考虑红色部分。为了使得 KL-Divergence 最小，咱们须要将红色部分（注意不包括前面的负号）最大化。不过因为仍然依赖于真实分布，因此红色部分仍是不能直接计算，可是这是一个咱们熟悉的形式：某个函数关于真实分布的指望值，而咱们手里面又有从真实分布里 IID 采样出来的  个数据点，因此能够直接使用标准的近似方法：用数据的经验分布指望来代替原始的指望：
因而咱们很天然地就获得了最大似然估计的目标函数。这也从另外一个角度解释了最大似然估计的合理性：同时对离散型和连续型的数据都适用。
接下来咱们再回到 Gaussian 分布，将 Gaussian 的几率密度函数 (eq: 1) 代入最大似然估计的目标函数，能够获得：
eq: 2 »
因为本文主要讨论协方差矩阵的估计，咱们这里就无论均值矩阵  的估计了，方法上是相似的，下面式子中的  将表示已知的或者一样经过最大似然估计计算出来的均值。为了找到使得上式最大化的协方差矩阵，咱们将它对于  求导：
另导数等于零，便可获得最优解：
eq: 3 »
因而  就是协方差矩阵的最大似然估计。到这里问题就出来了，再回到 Gaussian 分布的几率密度函数 (eq: 1)，若是把  带进去的话，能够看到是须要对它进行求逆的，这也是我被问得最多的问题：因为  singular 了，致使求逆失败，因而代码也就运行失败了。
从 (eq: 3) 能够大体看到何时  会 singular：因为该式的形式，能够知道  的 rank 最大不会超过，因此若是，也就是数据的维度大于数据点的个数的话，获得的估计值  确定是不满秩的，因而就 singular 了。固然，即便不是这么极端的状况，若是  比较小或者比较大的时候，仍然会有很大的概率出现不满秩的。
这个问题还能够从更抽象一点的角度来看。咱们以前计算过估计 Covariance 矩阵其实是要估计  这么多个参数，若是  比较大的话，这个数目也会变得很是多，也就意味着模型的复杂度很大，对于很复杂的模型，若是训练数据的量不够多，则没法肯定一个惟一的最优模型出来，这是机器学习中很是常见的一个问题，一般称为 overfitting。解决 overfitting 问题的方法有很多，最直接的就是用更多的训练数据（也就是增大），固然，有时候训练数据只有那么多，因而就只好从反面入手，下降模型复杂度。
在下降模型复杂度方面最经常使用的方法大概是正则化了，而正则化最简单的例子也许是 Ridge Regression 了，经过在目标函数后面添加一项关于参数的 -norm 的项来对模型参数进行限制；此外也有诸如 LASSO 之类的使用 -norm 来实现稀疏性的正则化，这个咱们在以前也曾经简单介绍过。
在介绍估计 Covariance 的正则化方法以前，咱们首先来看一种更直接的方法，直接限制模型的复杂度：特别地，对于咱们这里 Gaussian 的例子，好比，咱们能够要求协方差矩阵是一个对角阵，这样一来，对于协方差矩阵，咱们须要估计的参数个数就变成了  个对角元素，而不是原来的  个。大大减小参数个数的状况下，overfitting 的可能性也极大地下降了。
注意 Covariance Matrix 的非对角线上的元素是对应了某两个维度之间的相关性 (Correlation) 的，为零就表示不相关。而对于 Gaussian 分布来讲，不相关和独立 (Independence) 是等价的，因此说在咱们的假设下，向量  的各个维度之间是相互独立的，换句话说，他们的联合几率密度能够分解为每一个维度各自的几率密度（仍然是 Gaussian 分布）的乘积，这个特性使得咱们能够很方便地对协方差矩阵进行估计：只有每一个维度上当作一维的 Gaussian 分布各自独立地进行参数估计就能够了。注意这个时候除非某个维度上的坐标所有是相等的，不然 Covariance 矩阵对角线上的元素都能保证大于零，也就不会出现 singular 的状况了。
那么直观上将协方差矩阵限制为对角阵会产生什么样的模型呢？咱们不妨看一个二维的简单例子。下面的代码利用了 scikit-learn 里的 GMM 模块来估计一个单个 component 的 Gaussian，它在选项里已经提供了 Covariance 的限制，这里咱们除了原始的 full 和对角阵的 diag 以外，还给出一个 spherical，它假定协方差矩阵是  这样的形式，也就是说，不只是对角阵，并且全部对角元素都相等。下面的代码把三种状况下 fit 出来的模型画出来，在图 (fig: 1)中能够看到。
1. ```
import numpy as np 
```
2. ```
import pylab as pl 
```
3. ```
from sklearn import mixture 
```
5. ```
n_samples = 300 
```
6. ```
c_types = ['full', 'diag', 'spherical'] 
```
8. ```
np.random.seed(0) 
```
9. ```
C = np.array([[0., -0.7], [3.5, 1.7]]) 
```
10. ```
X_train = np.dot(np.random.randn(n_samples, 2), C) 
```
12. ```
pl.figure(dpi=100,figsize=(3,3)) 
```
13. ```
pl.scatter(X_train[:, 0], X_train[:, 1], .8) 
```
14. ```
pl.axis('tight') 
```
15. ```
pl.savefig('GaussianFit-data.svg') 
```
16. ```
pl.close() 
```
18. ```
for c_type in c_types: 
```
19. ```
 clf = mixture.GMM(n_components=1, covariance_type=c_type) 
```
20. ```
 clf.fit(X_train) 
```
22. ```
 x = np.linspace(-15.0, 20.0, num=200) 
```
23. ```
 y = np.linspace(-10.0, 10.0, num=200) 
```
24. ```
 X, Y = np.meshgrid(x, y) 
```
25. ```
 XX = np.c_[X.ravel(), Y.ravel()] 
```
26. ```
 Z = np.log(-clf.eval(XX)[0]) 
```
27. ```
 Z = Z.reshape(X.shape) 
```
29. ```
 pl.figure(dpi=100,figsize=(3,3)) 
```
30. ```
 CS = pl.contour(X, Y, Z) 
```
31. ```
 pl.scatter(X_train[:, 0], X_train[:, 1], .8) 
```
33. ```
 pl.axis('tight') 
```
34. ```
 pl.savefig('GaussianFit-%s.svg' % c_type) 
```
35. ```
 pl.close() 
```
fig: 1 »

Training samples.
Spherical Covariance.
Diagonal Covariance.
Full Covariance.
能够看到，从直观上来说：spherical Covariance 的时候就是各个维度上的方差都是同样的，因此造成的等高线在二维上是正圆。而 diagonal Covariance 则能够 capture 更多的特征，由于它容许造成椭圆的等高线，可是它的限制是这个椭圆是不能旋转的，永远是正放的。Full Covariance 是 fit 最好的，可是前提是数据量足够来估计这么多的参数，虽然对于 2 维这样的低维度一般并非问题，可是维度高了以后就常常会变得很困难。
通常来讲，使用对角线的协方差矩阵是一种不错的下降模型复杂度的方式，在 scikit-learn 的 GMM 模块中这甚至是默认的。不过若是有时候你以为这样的限制实在是和你的真实数据想去甚远的话，也有另外的处理方式。最简单的处理办法是在求出了 Covariance 矩阵的估计值以后直接加上一个，也就是在对角线上统一加上一个很小的正数，一般称做正则化系数。例如，Matlab 自带的统计工具箱里的 gmdistribution.fit 函数就有一个可选参数叫作Regularize，就是咱们这里说的。为何加上  以后就不会 singular 了呢？首先  自己至少确定是半正定的，因此：
因此这样以后确定就是正定的了。此外在 scikit-learn 的 GMM 学习的函数中也有一个相似的参数 min_covar，它的文档里说的是"Floor on the diagonal of the covariance matrix to prevent overfitting. Defaults to 1e-3."感受好像是若是原来的 Covariance 矩阵对角线元素小于该参数的时候就将它设置为该参数，不过这样的作法是否是必定产生正定的协方差矩阵彷佛就不像上面那种那样直观能够得出结论了。
不过这样的作法虽然可以解决 singular 的问题，可是总让人有些难以信服：你莫名其妙地在个人协方差矩阵上加上了一些东西，这究竟是什么意思呢？最简单的解释能够从式 (eq: 1) 那里协方差矩阵你的地方来看，对于矩阵求逆或者解方程的时候出现 singular 的状况，加上一个也算是数值上的一种标准处理方式，叫作 Tikhonov Regularization，Ridge Regression 自己也能够经过这个角度来解释。从数值计算的角度来讲，这样处理可以增长数值稳定性，直观地来说，稳定性是指  元素的微小数值变化，不会使得求逆（或者解方程之类的）以后的解产生巨大的数值变化，若是  是 singular 的话，一般就不具备这样的稳定性了。此外，随着 regularization coefficient  逐渐趋向于零，对应的解也会逐渐趋向于真实的解。
若是这个解释还不能令出信服的话，咱们还能够从 Bayes 推断的角度来看这个问题。不过这自己是一个很大的话题，我在这里只能简略地讲一个思路，想了解更多的同窗能够参考《Pattern Recognition and Machine Learning》第二章或者其余相关的书籍。
总的来讲，咱们这里要经过 Maximum a posteriori (MAP) Estimation 来对协方差矩阵进行估计。作法是首先为  定义一个先验分布 (prior) ，而后对于数据，咱们根据 Bayes 公式，有
eq: 4 »
其中等式左边称做后验 (posterior)，右边红色部分是似然 (likelihood)，蓝色部分是用于将几率分布 normalize 到全几率为 1 的，一般并不须要直接经过  去计算，而只要对求出的几率分布进行归一化就能够了。MAP 的思路是将 posterior 最大的那个  做为最佳估计值，直观上很好理解，就是给定了数据点以后"最可能的"那个。而计算上咱们是经过 (eq: 4) 来实现的，由于分母是与  无关的，因此在计算的时候（由于只要比较相对大小嘛）忽略掉，这样一来，就能够看到，其实 MAP 方法和刚才咱们用的最大似然 (MLE) 方法惟一的区别就在于前面乘了一个先验分布。
虽然从最终的计算公式上来看差异很细微，可是 MAP 却有不少好处。和本文最相关的固然是先验的引入，这在很大程度上和咱们平时用正则化的目的是同样的：加入咱们本身对模型的先验知识或者假设或者要求，这样能够避免在数据不够的时候产生 overfitting 。此外还有另外一个值得一提的好玩特性就是 MAP 能够在线计算：注意到后验分布自己也是关于  的一个合法分布，因此在咱们计算出后验以后，若是又获得了新的训练数据，那么咱们能够将以前算出来的后验又做为先验，代入新的一轮的计算，这样能够一直不断地重复下去。:D
不过，虽然 MAP 框架看上去很美，可是在 general 的状况下计算却可能会变得很复杂，所以，虽然理论上来讲，咱们能够把任何先验知识经过先验分布的方式加入到模型中来，可是在实际中先验的选择每每是根据"哪一个更好计算"而不是根据"哪一个更合理"的准则来作出的。基于这个准则，有了 Conjugate Prior 的概念，简单地来讲，若是一个 prior 和 likelihood 相乘再归一化以后获得的 posterior 形式上是和 prior 同样的话，那么它就是该 likelihood 的一个 conjugate prior。
选择 conjugate prior 的好处是显而易见的，因为先验和后验的形式是同样的，都是某一类型的分布，只是参数的值发生了变化，因此说能够当作是在观察到数据以后对模型参数进行修正（注意这里的模型是关于  的，而  自己又是关于数据的模型——一个 Gaussian 分布——的参数，不要搞混淆了），而且这种修正有时候能够获得很是直观的解释，不过我在这里就不细说了。此外，因为形式没有发生变化，因此在使用在线计算的时候，每一次观察到新的数据所使用的计算公式都仍是同样的。
回到咱们的问题，关于 Gaussian 分布的协方差矩阵的 Conjugate Prior 是一个叫作 Inverse Wishart Distribution 的家伙，它的几率密度函数是这个样子的：
其中  和参数  都是  的正定矩阵，而  则是 Multivariate Gamma Function。是否是看上去有点恐怖？^_^bbbb我也以为有点恐怖……接下来让咱们来看 MAP，首先将这个 prior 乘以咱们的 likelihood，和 MLE 同样的，咱们在  空间中计算，因此对整个式子取对数，而后丢掉和  无关的常数项（由于是关于  最大化嘛），获得下面的式子：
能够看到如今已经变得不是那么复杂了，形式上和咱们以前的 (eq: 2) 是同样的，只是多了红色的一项，接下来对  进行求导并令导数等于零，相似地获得以下方程：
因而获得最优解为
这里  是以前的最大似然估计。接下来若是咱们将咱们 prior 中的参数  选为，就能够获得 MAP 估计为：
就获得了咱们刚才的正则化的形式。有一点细微的区别就是这里前面多了一个全局的缩放系数，不过这样的系数是没有关系的，能够直接丢掉，由于它总也会在最终归一化的时候被 Gaussian 分布前面的系数给抵消掉。这样一来，咱们就获得了这种正则化方法的一个看起来比较靠谱的解释：实际上就是咱们先验地假设协方差矩阵知足一个均值为 spherical 形状（全部对角元素都相等的一个对角阵）的 Inverse Wishart 分布（关于 Inverse Wishart 分布的均值的公式能够参见 wikipedia），而后经过数据来对这个先验进行修正获得后验几率，并从后验中取了几率（密度）最大的那个  做为最终的估计。
这种方式比暴力地直接强制要求协方差矩阵具备 diagonal 或者 spherical 形式要温柔得多，而且随着训练数据的增多，先验的影响也会逐渐减少，从而趋向于获得真实的参数，而强制结构限制的结果则是不论你最后经过什么手段搞到了多少数据，若是真实模型自己不是 diagonal 或者 spherical 的话，那么你的参数估计的 bias 将永远都没法被消除，这个咱们已经在图 (fig: 1)中直观地看到了。固然，在训练数据自己就不足的状况下，两种选择谁好谁好就没有定论了，不要忘了，强制结构限制的计算复杂度要小不少。:)
最后，对于不想看长篇大论的人，我作一下简单的总结：若是碰到估计 Gaussian 的协方差矩阵时结果 singular 的状况，解决方法通常有两种：

直接限制协方差矩阵的结构，好比要求它是一个对角阵，这种状况下数据的各个维度将会是独立的，因此能够把每一个维度当作一个一维的 Gaussian 分布来单独进行估计。
在估计出来的协方差矩阵的对角线上统一加上一个很小的正数，使得它变成正定的。这种方法的合理性能够经过正则化或者 MAP 估计来进行解释。

而后再补充几点不知道该放在哪里的注意事项：

prior、posterior 和 likelihood 比较容易搞混淆，记住 prior 和 posterior 都是关于参数（好比咱们这里就是）的分布，prior 是咱们所假设的参数原本的分布，而 posterior 则是在观察到训练数据以后获得的条件分布。而 likelihood 则彻底不同，一方面它是关于数据的，另外一方面它没有被归一化，因此也并非一个合法的几率分布。
标量值函数关于矩阵变量的求导原理上其实就是和多元标量值函数求导同样的，不过求起来很是繁琐，通常不会本身算，网上能够找到一个叫作《Matrix Cookbook》的小册子，里面有搜集了各类经常使用的矩阵求导的公式，基本上直接从那里查询就能够解决大部分问题了。

(番外篇): Expectation Maximization
Expectation Maximization (EM) 是一种以迭代的方式来解决一类特殊最大似然 (Maximum Likelihood) 问题的方法，这类问题一般是没法直接求得最优解，可是若是引入隐含变量，在已知隐含变量的值的状况下，就能够转化为简单的状况，直接求得最大似然解。
咱们会看到，上一次说到的 Gaussian Mixture Model 的迭代求解方法能够算是 EM 算法最典型的应用，而最开始说的 K-means 其实也能够看做是 Gaussian Mixture Model 的一个变种（固定全部的，并令便可）。然而 EM 其实是一种通用的算法（或者说是框架），能够用来解决不少相似的问题，咱们最后将以一个中文分词的例子来讲明这一点。
为了不问题变得太抽象，咱们仍是先从上一次的 Gaussian Mixture Model 提及吧。回顾一下咱们以前要解决的问题：求如下 Log-likelihood function 的最大值：
```
 
```
可是因为在  函数里面又有加和，无法直接求。考虑一下 GMM 生成 sample 的过程：先选一个 Component ，而后再从这个 Component 所对应的那个普通的 Gaussian Distribution 里进行 sample 。咱们能够很天然地引入一个隐含变量  ，这是一个  维向量，若是第  个 Component 被选中了，那么咱们讲其第  个元素置为 1 ，其他的全为 0 。那么，再来看看，若是除了以前的 sample  的值以外，咱们同时还知道了每一个  所对应的隐含变量  的值，状况又会变成怎么样呢？
由于咱们同时观察到了  和  ，因此咱们如今要最大化的似然函数也变为  。注意到  能够表示为：

而  的几率，当  的第  个元素为 1 的时候，亦即第  个 Component 被选中的时候，这个几率为  ，统一地写出来就是：

带入上面个式子，咱们获得  的几率是一个大的乘积式（对比以前  是一个和式）。再替换到最开始的那个 Log-likelihood function 中，获得新的同时关于 sample  和隐含变量  的 Log-likelihood:

状况瞬间逆转，如今  和求和符号换了个位置，直接做用在普通的高斯分布上了，一会儿就变成了能够直接求解的问题。不过，事情之因此能发展得如此顺利，彻底依赖于一个咱们伪造的假设：隐含变量的值已知。然而实际上咱们并不知道这个值。问题的结症在这里了，若是咱们有了这个值，全部的问题都迎刃而解了。回想一下，在相似的地方，咱们是如何处理这样的状况的呢？一个很相似的地方就是（好比，在数据挖掘中）处理缺失数据的状况，通常有几种办法：
1. 用取值范围类的随机值代替。
2. 用平均值代替。
3. 填 0 或者其余特殊值。
这里咱们采起一种相似于平均值的办法：取指望。由于这里咱们至少有 sample 的值，所以咱们能够把这个信息利用起来，针对后验几率来取指望。前面说过，的每个元素只有 0 和 1 两种取值，所以按照指望的公式写出来就是：
```
 
```
中间用贝叶斯定理变了一下形，最后获得的式子正是咱们在漫谈 GMM 中定义的。所以，对于上面那个能够直接求解的 Log-likelihood function 中的，咱们只要用其指望亦即代替便可。
到这里为止，看起来是一个很是完美的方法，不过仔细一看，发现还有一个漏洞：以前的 Log-likelihood function 能够直接求最大值是创建在是已知状况下，如今虽然咱们用来代替了，可是实际上是一个反过来以很是复杂的关系依赖所要求参数的一个式子，而不是一个"已知的数值"。解决这个问题的办法就是迭代。
到此为止，咱们就勾勒出了整个 EM 算法的结构，同时也把上一次所讲的求解 GMM 的方法又推导了一遍。EM 名字也来源于此：
1. Expectation 一步对应于求关于后验几率的指望亦即；
2. Maximization 一步则对应于接下来的正常的最大似然的方法估计相关的参数亦即、和。
若是你尚未厌烦这些公式的话，让咱们不妨再稍微花一点时间，从偏理论一点的角度来简略地证实一下 EM 这个迭代的过程每一步都会对结果有所改进，除非已经达到了一个（至少是局部的）最优解。
如今咱们的讨论将不局限于 GMM ，并使用一些稍微紧凑一点的符号。用表示全部的 sample ，同时用表示全部对应的隐含变量。咱们的问题是要经过最大似然的方法估计出中的参数。在这里咱们假设这个问题很困难，可是要优化却很容易。这就是 EM 算法可以解决的那一类问题。
如今咱们引入一个关于隐含变量的分布。注意到对于任意的，咱们均可以对 Log-likelihood Function 做以下分解：
```
 
```
其中是分布和之间的 Kullback-Leibler divergence 。因为 Kullback-Leibler divergence 是非负的，而且只有当两个分布彻底相同的时候才会取到 0 。所以咱们能够获得关系，亦即是的一个下界。
如今考虑 EM 的迭代过程，记上一次迭代得出的参数为，如今咱们要选取以令最大，因为并不依赖于，所以的上限（在固定的时候）是一个定值，它取到这个最大值的条件就是 Kullback-Leibler divergence 为零，亦即等于后验几率。把它带入到的表达式中能够获得
```
 
```
其中 const 是常量，而  则正是咱们以前所获得的"同时包含了 sample 和隐含变量的 Log-likelihood function 关于后验几率的指望"，所以这个对应到 EM 中的"E"一步。
在接下来的"M"一步中，咱们讲固定住分布  ，再选取合适的  以将  最大化，此次其上界  也依赖于变量  ，并会随着  的增大而增大（由于咱们有前面的不等式成立）。通常状况下  增大的量会比  要多一些，这时候 Kullback-Leibler divergence 在新的参数  下又不为零了，所以咱们能够进入下一轮迭代，从新回到"E"一步去求新的  ；另外一方面，若是这里 Kullback-Leibler divergence 在新的参数下仍是等于 0 ，那么说明咱们已经达到了一个（至少是局部的）最优解，迭代过程能够结束了。
上面的推导中咱们能够看到每一次迭代 E 和 M 两个步骤都是在对解进行改进，所以迭代的过程当中获得的 likelihood 会逐渐逼近（至少是局部的）最优值。另外，在 M 一步中除了用最大似然以外，也能够引入先验使用 Maximum a Posteriori (MAP) 的方法来作。还有一些很困难的问题，甚至在迭代的过程当中 M 一步也不能直接求出最大值，这里咱们经过把要求放宽──并不必定要求出最大值，只要可以获得比旧的值更好的结果便可，这样的作法一般称做 Generalized EM (GEM)。
固然，一开始咱们就说了，EM 是一个通用的算法，并不仅是用来求解 GMM 。下面咱们就来举一个用 EM 来作中文分词的例子。这个例子来源于论文 Self-supervised Chinese Word Segmentation 。我在上次 MSTC 搜索引擎系列小课堂讲文本处理的时候提到过。这里为了把注意力集中到 EM 上，我略去一些细节的东西，简单地描述一下基本的模型。
如今咱们有一个字符序列，并但愿获得一个模型（依赖于参数）能用来将其词的序列  。按照生成模型的方式来考虑，将  当作是由  生成的序列的话，咱们天然但愿找到那个生成  的可能性最大的  ，亦即经过最大似然的方式来估计参数  。
然而咱们不知道似然函数  应该如何去优化，所以咱们引入 latent variable  ，若是咱们知道  的话，似然值很容易求得：

其中  的值是从模型  中已知的。可是如今咱们不知道  的值，所以咱们转而取其关于后验几率的指望：

而后将这个指望针对  最大化即完成了 EM 的一次迭代。具体的作法一般是先把一个初始文本（好比全部的  的集合）按照 N-gram 分割（N-gram 在讲 K-medoids 的那篇文章中介绍过）为，造成一个最初的辞典，而模型  的参数  实际上就描述了各个 N-gram 的几率  ，初始值能够直接取为频率值。有了辞典以后对于任意的  ，咱们能够根据辞典枚举出全部可能的分割  ，而每一个分割的后验几率  就是其中的单词的几率乘积。其余的就是标准的 EM 算法的迭代步骤了。
实际上在实际产品中咱们会使用一些带了更多启发式规则的分词方法（好比 MMSEG），而这个方法更大的用处在于从一堆文本中"学习"出一个词典来（也就是  ），而且这个过程是全自动的，主要有两个优势：
1. 不须要人参与手工分词、标记等。
2. 能自动从文本中学习，换句话说，你给它一些某个领域的专业文章，它能从中学习到这个领域的专业词汇来。
无论怎样，以上两点看起来都很是诱人。不过理论离实际每每仍是有很多差距的。我不知道实际分词系统中有没有在用这样的方法来作训练的。以前我曾经用 Wikipedia （繁体中文）上抓下来的一些数据跑过一下小规模的试验，结果还能够，可是也并不如想像中的完美。由于当时也并无弄明白 EM 是怎么一回事，加上这个过程自己计算负荷仍是很是大的，也就没有深刻下去。也许之后有时间能够再尝试一下。
总之，EM 是一种应用普遍的算法，虽然讲了点理论也讲了例子，可是一没有贴代码二没有贴图片，彷佛实在是不像个人做风，不过精力和篇幅有限，也只能到此为止了。

Spectral Clustering

若是说 K-means 和 GMM 这些聚类的方法是古代流行的算法的话，那么此次要讲的 Spectral Clustering 就能够算是现代流行的算法了，中文一般称为"谱聚类"。因为使用的矩阵的细微差异，谱聚类实际上能够说是一"类"算法。

Spectral Clustering 和传统的聚类方法（例如 K-means）比起来有很多优势：

和 K-medoids 相似，Spectral Clustering 只须要数据之间的类似度矩阵就能够了，而没必要像 K-means 那样要求数据必须是 N 维欧氏空间中的向量。
因为抓住了主要矛盾，忽略了次要的东西，所以比传统的聚类算法更加健壮一些，对于不规则的偏差数据不是那么敏感，并且 performance 也要好一些。许多实验都证实了这一点。事实上，在各类现代聚类算法的比较中，K-means 一般都是做为 baseline 而存在的。
计算复杂度比 K-means 要小，特别是在像文本数据或者平凡的图像数据这样维度很是高的数据上运行的时候。

忽然冒出这么一个要求比 K-means 要少，结果比 K-means 要好，算得还比 K-means 快的东西，实在是让人不得不怀疑是否是江湖骗子啊。因此，是骡子是马，先拉出来溜溜再说。不过，在K-medoids 那篇文章中曾经实际跑过 K-medoids 算法，最后的结果也就是一个 accuracy ，一个数字又不能画成图表之类的，看起来实在是没意思，并且 K-means 跑起来实在是太慢了，因此这里我仍是稍微偷懒一下，直接引用一下一篇论文里的结果吧。

结果来自论文 Document clustering using locality preserving indexing 这篇论文，这篇论文其实是提的另外一种聚类方法（下次若是有机会也会讲到），不过在它的实验中也有 K-means 和 Spectral Clustering 这两组数据，抽取出来以下所示：

k	TDT2		Reuters-21578
k	K-means	SC	K-means	SC
2	0.989	0.998	0.871	0.923
3	0.974	0.996	0.775	0.816
4	0.959	0.996	0.732	0.793
…
9	0.852	0.984	0.553	0.625
10	0.835	0.979	0.545	0.615

其中 TDT2 和 Reuters-21578 分别是两个被普遍使用的标准文本数据集，虽然在不一样的数据集上得出来的结果并不能直接拿来比较，可是在同一数据集上 K-means 和 SC (Spectral Clustering) 的结果对比就一目了然了。实验中分别抽取了这两个数据集中若干个类别（从 2 类到 10 类）的数据进行聚类，得出的 accuracy 分别在上表中列出（我偷懒没有所有列出来）。能够看到，Spectral Clustering 这里完胜 K-means 。

这么强大的算法，再戴上"谱聚类"这么一个高深莫测的名号，若不是模型无比复杂、包罗宇宙，就确定是某镇山之宝、不传秘籍吧？其实不是这样，Spectral Clustering 无论从模型上仍是从实现上都并不复杂，只须要能求矩阵的特征值和特征向量便可──而这是一个很是基本的运算，任何一个号称提供线性代数运算支持的库都理应有这样的功能。而关于 Spectral Clustering 的秘籍更是满街都是，随便从地摊上找来一本，翻开即可以看到 Spectral Clustering 算法的全貌：

根据数据构造一个 Graph ，Graph 的每个节点对应一个数据点，将类似的点链接起来，而且边的权重用于表示数据之间的类似度。把这个 Graph 用邻接矩阵的形式表示出来，记为。一个最偷懒的办法就是：直接用咱们前面在 K-medoids 中用的类似度矩阵做为。
把的每一列元素加起来获得个数，把它们放在对角线上（其余地方都是零），组成一个的矩阵，记为。并令。
求出的前个特征值（在本文中，除非特殊说明，不然"前个"指按照特征值的大小从小到大的顺序）以及对应的特征向量。
把这个特征（列）向量排列在一块儿组成一个的矩阵，将其中每一行看做维空间中的一个向量，并使用 K-means 算法进行聚类。聚类的结果中每一行所属的类别就是原来 Graph 中的节点亦即最初的个数据点分别所属的类别。

就是这么几步，把数据作了一些诡异的变换，而后还在背后偷偷地调用了 K-means 。到此为止，你已经能够拿着它上街去招摇撞骗了。不过，若是你仍是以为不太靠谱的话，不妨再接着往下看，咱们再来聊一聊 Spectral Clustering 那几个"诡异变换"背后的道理何在。

其实，若是你熟悉 Dimensional Reduction （降维）的话，大概已经看出来了，Spectral Clustering 其实就是经过 Laplacian Eigenmap 的降维方式降维以后再作 K-means 的一个过程──听起来土多了。不过，为何要恰好降到维呢？其实整个模型还能够从另外一个角度导出来，因此，让咱们不妨先跑一下题。

在 Image Processing （我好像以前有据说我对这个领域深恶痛绝？）里有一个问题就是对图像进行 Segmentation （区域分割），也就是让类似的像素组成一个区域，好比，咱们通常但愿一张照片里面的人（前景）和背景被分割到不一样的区域中。在 Image Processing 领域里已经有许多自动或半自动的算法来解决这个问题，而且有很多方法和 Clustering 有密切联系。好比咱们在谈 Vector Quantization 的时候就曾经用 K-means 来把颜色类似的像素聚类到一块儿，不过那还不是真正的 Segmentation ，由于若是仅仅是考虑颜色类似的话，图片上位置离得很远的像素也有可能被聚到同一类中，咱们一般并不会把这样一些"游离"的像素构成的东西称为一个"区域"，但这个问题其实也很好解决：只要在聚类用的 feature 中加入位置信息（例如，原来是使用 R、G、B 三个值来表示一个像素，如今加入 x、y 两个新的值）便可。

另外一方面，还有一个常常被研究的问题就是 Graph Cut ，简单地说就是把一个 Graph 的一些边切断，让他被打散成一些独立联通的 sub-Graph ，而这些被切断的边的权值的总和就被称为 Cut值。若是用一张图片中的全部像素来组成一个 Graph ，并把（好比，颜色和位置上）类似的节点链接起来，边上的权值表示类似程度，那么把图片分割为几个区域的问题实际上等价于把 Graph 分割为几个 sub-Graph 的问题，而且咱们能够要求分割所得的 Cut 值最小，亦即：那些被切断的边的权值之和最小，直观上咱们能够知道，权重比较大的边没有被切断，表示比较类似的点被保留在了同一个 sub-Graph 中，而彼此之间联系不大的点则被分割开来。咱们能够认为这样一种分割方式是比较好的。

实际上，抛开图像分割的问题不谈，在 Graph Cut 相关的一系列问题中，Minimum cut （最小割）自己就是一个被普遍研究的问题，而且有成熟的算法来求解。只是单纯的最小割在这里一般并非特别适用，不少时候只是简单地把和其余像素联系最弱的那一个像素给分割出去了，相反，咱们一般更但愿分割出来的区域（的大小）要相对均匀一些，而不是一些很大的区块和一些几乎是孤立的点。为此，又有许多替代的算法提出来，如 Ratio Cut 、Normalized Cut 等。

不过，在继续讨论以前，咱们仍是先来定义一下符号，由于仅凭文字仍是很难表述清楚。将 Graph 表示为邻接矩阵的形式，记为，其中是节点到节点的权值，若是两个节点不是相连的，权值为零。设和为 Graph 的两个子集（没有交集），那么二者之间的 cut 能够正式定义为：

先考虑最简单的状况，若是把一个 Graph 分割为两个部分的话，那么 Minimum cut 就是要最小化（其中表示的补集）。可是因为这样常常会出现孤立节点被分割出来的状况，所以又出现了 RatioCut ：

以及 NormalizedCut ：

其中表示中的节点数目，而。二者均可以算做的"大小"的一种度量，经过在分母上放置这样的项，就能够有效地防止孤立点的状况出现，达到相对平均一些的分割。事实上，Jianbo Shi 的这篇 PAMI paper：Normalized Cuts and Image Segmentation 正是把 NormalizedCut 用在图像分割上了。

搬出 RatioCut 和 NormalizedCut 是由于它们和这里的 Spectral Clustering 实际上有很是紧密的联系。看看 RatioCut ，式子虽然简单，可是要最小化它倒是一个 NP 难问题，不方便求解，为了找到解决办法，让咱们先来作作变形。

令表示 Graph 的全部节点的集合，首先定义一个维向量：

再回忆一下咱们最开始定义的矩阵，其实它有一个名字，叫作 Graph Laplacian ，不过，咱们后面能够看到，其实有好几个相似的矩阵都叫作这个名字：

Usually, every author just calls "his" matrix the graph Laplacian.

其实也能够理解，就好象如今全部的厂家都说本身的技术是"云计算"同样。这个有一个性质就是：

这个是对任意向量都成立的，很好证实，只要按照定义展开就能够获得了。把咱们刚才定义的那个带进去，就能够获得

另外，若是令为各个元素全为 1 的向量的话，直接展开能够很容易获得和。因为是一个常量，所以最小化 RatioCut 就等价于最小化，固然，要记得加上附加条件以及。

问题转化到这个样子就好求了，由于有一个叫作 Rayleigh quotient 的东西：

他的最大值和最小值分别等于矩阵的最大的那个特征值和最小的那个特征值，而且极值在等于对应的特征向量时取到。因为是常数，所以最小化实际上也就等价于最小化，不过因为的最小的特征值为零，而且对应的特征向量正好为（咱们这里仅考虑 Graph 是联通的状况），不知足的条件，所以咱们取第二个小的特征值，以及对应的特征向量。

到这一步，咱们看起来好像是很容易地解决了前面那个 NP 难问题，其实是咱们耍了一个把戏：以前的问题之因此 NP 难是由于向量的元素只能取两个值和中的一个，是一个离散的问题，而咱们求的的特征向量其中的元素能够是任意实数，就是说咱们将原来的问题限制放宽了。那如何获得原来的解呢？一个最简单的办法就是看的每一个元素是大于零仍是小于零，将他们分别对应到离散状况的和，不过咱们也能够采起稍微复杂一点的办法，用的 K-means 来将的元素聚为两类。

到此为止，已经有 Spectral Clustering 的影子了：求特征值，再对特征向量进行 K-means 聚类。实际上，从两类的问题推广到 k 类的问题（数学推导我就再也不详细写了），咱们就获得了和以前的 Spectral Clustering 如出一辙的步骤：求特征值并取前 k 个最小的，将对应的特征向量排列起来，再按行进行 K-means 聚类。分绝不差！

用相似的办法，NormalizedCut 也能够等价到 Spectral Clustering 不过此次我就再也不讲那么多了，感兴趣的话（还包括其余一些形式的 Graph Laplacian 以及 Spectral Clustering 和 Random walk 的关系），能够去看这篇 Tutorial ：A Tutorial on Spectral Clustering 。

为了缓和一下气氛，我决定贴一下 Spectral Clustering 的一个简单的 Matlab 实现：

function idx = spectral_clustering(W, k)

 D = diag(sum(W));

 L = D-W;

 opt = struct('issym', true, 'isreal', true);

 [V dummy] = eigs(L, D, k, 'SM', opt);

 idx = kmeans(V, k);

end

最后，咱们再来看一下本文一开始说的 Spectral Clustering 的几个优势：

只须要数据的类似度矩阵就能够了。这个是显然的，由于 Spectral Clustering 所须要的全部信息都包含在中。不过通常并不老是等于最初的类似度矩阵——回忆一下，是咱们构造出来的 Graph 的邻接矩阵表示，一般咱们在构造 Graph 的时候为了方便进行聚类，更增强到"局部"的连通性，亦即主要考虑把类似的点链接在一块儿，好比，咱们设置一个阈值，若是两个点的类似度小于这个阈值，就把他们看做是不链接的。另外一种构造 Graph 的方法是将 n 个与节点最类似的点与其链接起来。
抓住了主要矛盾，忽略了次要的东西，Performance 比传统的 K-means 要好。实际上 Spectral Clustering 是在用特征向量的元素来表示原来的数据，并在这种"更好的表示形式"上进行 K-means 。实际上这种"更好的表示形式"是用 Laplacian Eigenmap 进行降维的后的结果，若是有机会，下次谈 Dimensionality Reduction 的时候会详细讲到。而降维的目的正是"抓住主要矛盾，忽略次要的东西"。
计算复杂度比 K-means 要小。这个在高维数据上表现尤其明显。例如文本数据，一般排列起来是维度很是高（好比，几千或者几万）的稀疏矩阵，对稀疏矩阵求特征值和特征向量有很高效的办法，获得的结果是一些 k 维的向量（一般 k 不会很大），在这些低维的数据上作 K-means 运算量很是小。可是对于原始数据直接作 K-means 的话，虽然最初的数据是稀疏矩阵，可是 K-means 中有一个求 Centroid 的运算，就是求一个平均值：许多稀疏的向量的平均值求出来并不必定仍是稀疏向量，事实上，在文本数据里，不少状况下求出来的 Centroid 向量是很是稠密，这时再计算向量之间的距离的时候，运算量就变得很是大，直接致使普通的 K-means 巨慢无比，而 Spectral Clustering 等工序更多的算法则迅速得多的结果。

说了这么多，好像有些乱，不过也只能到此打住了。最后再多嘴一句，Spectral Clustering 名字来源于 Spectral theory ，也就是用特征分解来分析问题的理论了。

UPDATE 2011.11.23: 有很多同窗问我关于代码的问题，这里更新两点主要的问题：

关于 LE 降维的维度和 Kmeans 聚类的类别数的关系：我上面的代码里，取成了同样的，可是这二者并非要求必定要同样的。最初 Spectral Cluster 是分析聚两类的状况，就降到 1 维，而后用 thresholding 的方法来分红两类。对于K 类的状况，天然的类比就是降到 K-1 维，这也是和 LDA 保持一致。由于 Laplacian 矩阵的特征向量有一个全一的向量（对应于特征值 0 ），因此能够求 K 个特征向量而后把特征值 0 对应的那个特征向量去掉。可是在实际中并不必定要保持这二者一致的，也就是说，这个降维的维度能够做为一个参数进行调节的，选择效果好的参数。
关于示例代码：注意除非我在这里注明了是发布某个 software 或者 package 什么的，不然这里贴的代码主要都是为了示例做用，为了只显示出算法的主要部分，所以一般会省略掉一些实现细节，能够当作是"可执行的伪代码"，不推荐直接用这里的代码去跑实验之类的（包括其余 post 里的相关代码）。除非你想本身试验一下具体实现和改进，不然能够直接找网上现成的专用的 package ，好比 Spectral Clustering 能够考虑这个包。

(番外篇): Dimensionality Reduction

因为老是有各类各样的琐事，这个系列的文章居然一会儿拖了好几个月，（实际上其余的日志我也写得比较少），如今决定仍是先把这篇降维的日志写完。我甚至都以及忘记了在这个系列中以前有没有讲过"特征"（feature）的概念了，这里不妨再稍微提一下。机器学习应用到各个领域里，会遇到许多不一样类型的数据要处理：图像、文本、音频视频以及物理、生物、化学等实验还有其余工业、商业以及军事上获得的各类数据，若是要为每一种类型的数据都设计独立的算法，那显然是很是不现实的事，所以，机器学习算法一般会采用一些标准的数据格式，最多见的一种格式就是每个数据对应欧几里德空间里的一个向量。

若是原始的数据格式不兼容，那么就须要首先进行转换，这个过程一般叫作"特征提取"（Feature Extraction），而获得的标准数据格式一般叫作 Feature 。例如，一个最简单的将一个文本 Document 转化为向量的方法以下：

选定特征空间，这里采用三维欧氏空间，三个维度（依次）分别由 to 、be 和 the 表示。
假设待提取的文档是"To be, or not to be: that is the question:"，首先对其进行一些预处理，例如去掉单词的时态后缀、滤掉标点符号等，获得"to be or not to be that be the question"。
统计三个维度所对应的单词出现的频率：to 2 次，be 3 次，the 1 次。
该文档对应的向量即 [2, 3, 1] 。

固然，在实际中咱们几乎不会这样人工设定空间的各个维度所对应的单词，而一般是从一个数据集中统计出全部出现的词，再将其中的一些挑选出来做为维度。怎样挑选呢？最简单的办法是根本不作任何挑选，或者简单地只是把出现频率过低的单词（维度）去掉。

不过，事实上咱们一般会作更复杂一些的处理，例如，若是你是在作 sentiment analysis ，那么你一般会更加关注语气很重的词，好比 "bad"、"terrible"、"awesome" 等的重要性就比普通的词要大，此时你能够为每个维度设一个权重，例如，给 "bad" 设置权重 2 ，那么出现 3 次的话，向量在该维度对应的值就为 2*3 = 6 。固然这样手工指定权重只在小范围内可行，若是要为数百万个维度指定权重，那显然是不可能的，另外一个稍微自动一点的办法是 tf-idf 。

tf 就是 Term Frequency ，就和刚才说的单词出现的次数差很少，而 idf 则是 Inverse Document Frequency ，一般使用以下公式进行计算：

这至关于自动计算出来的每一个单词的权重，其想法很简单：若是在许多文档中都出现的词（例如虚词、语气词等），它所包含的信息量一般会比较小，因此用以上的公式计算出来的权重也会相对较小；另外一方面，若是单词并非在不少文档中都出现的话，颇有可能就是出现的那些文档的重要特征，所以权重会相对大一些。

前面说了一堆，其实主要是想要对 feature 作一些"预"处理，让它更"好"一些，手工设置维度的权重当然是很人力，其实 tf-idf 也是在假定了原始 feature 是 document 、term 这样的形式（或者相似的模型）的状况下才适用的（例如，在门禁之类的系统里，feature 可能有声音、人脸图像以及指纹等数据，就不能这样来处理），所以也并不能说是一种通用的方法。

然而，既然机器学习的算法能够在不考虑原始数据的来源和意义的状况下工做，那么 feature 的处理应该也能够。事实也是如此，包括 feature selection 和 dimensionality reduction 两个 topic 都有许多很经典的算法。前者一般是选出重要的 feature 的维度（并抛弃不重要的维度），然后者则是更普遍意义上的把一个高维的向量映射为一个低维向量，亦即"降维"，获得的结果 feature 的值已经不必定是原始的值，也能够把 feature selection 看做是 dimensionality reduction 的一种特殊状况。举一个例子，tf-idf 实际上不算 feature selection ，由于它（一般）并无丢弃低权值的维度，而且处理事后的特征的每一个维度都被乘上了一个权值，再也不是原来的值了；可是它却能够被看做一种降维，虽然严格意义上来讲维度并无"下降"。简单来讲降维能够看做一个函数，其输入是一个 D 维的向量，输出是一个 M 维的向量。

按照机器学习的方法的一向做风，咱们须要定义目标函数并进行最优化。不一样的目标也就致使了不一样的降维算法，这也正是今天要讲的话题。

然而，咱们的目的到底是什么呢？一个比较直接的问题是原始的数据量维度过高，咱们没法处理，所以须要降维，此时咱们一般但愿在最大限度地下降数据的维度的前提下可以同时保证保留目标的重要的信息，就比如在作有损的数据压缩时但愿压缩比尽可能大可是质量损失不要太多。因而问题又转化为如何衡量对信息的保留程度了。

一个最直接的办法就是衡量 reconstruction error ，即

其中是所对应的低维表示再从新构造出来的高维形式，就至关因而压缩以后解压出来的结果，不过虽然有许多压缩方法都是无损的，就是说这个差值会等于零，可是大部分降维的结果都是有损的。不过咱们仍然但愿把上面的 reconstruction error 最小化。

另一种方式是简单地使用 variance 来衡量所包含信息量，例如，咱们要把一些 D 维的向量降为 1 维，那么咱们但愿这一维的 variance 达到最大化，亦即：

其中是降维函数。推而广之，若是是降为 2 维，那么我但愿第 2 维去关注第 1 维以外的 variance ，因此要求它在与第一维垂直的状况下也达到 variance 最大化。以此类推。

然而，当咱们把降维函数限定维线性的时候，两种途径会获得一样的结果，就是被普遍使用的 Principal Components Analysis(PCA) 。PCA 的降维函数是线性的，能够用一个维的矩阵来表示，所以，一个 D 维的向量通过线性变换以后获得一个 M 维向量，就是降维的结果。把原始数据按行排列为一个维的矩阵，则就是降维后的维的数据矩阵，目标是使其 covariance 矩阵最大。在数据被规则化（即减去其平均值）过的状况下，协方差矩阵（covariance），固然矩阵不是一个数，不能直接最大化，若是咱们采用矩阵的 Trace （亦即其对角线上元素的和）来衡量其大小的话，要对求最大化，只须要求出的特征值和特征向量，将 M 个最大的特征值所对应的特征向量按列排列起来组成线性变换矩阵便可。这也就是 PCA 的求解过程，获得的降维矩阵能够直接用到新的数据上。若是熟悉 Latent Semantic Analysis (LSA) 的话，大概已经看出 PCA 和 Singular Value Decomposition (SVD) 以及 LSA 之间的关系了。如下这张图（引自《The Elements of Statistical Learning》）能够直观地看出来 PCA 作了什么，对于一些原始为二维的数据，PCA 首先找到了 variance 最大的那一个方向：

PCA 做为一种经典的降维方法，被普遍地应用于机器学习、计算机视觉以及信息检索等各个领域，其地位相似于聚类中的 K-means ，在如今关于降维等算法的研究中也常常被做为 baseline 出现。然而，PCA 也有一些比较明显的缺点：一个就是 PCA 降维是线性变换，虽然线性变换计算方便，而且能够很容易地推广到新的数据上，然而有些时候线性变换也许并不合适，为此有许多扩展被提出来，其中一个就是 Kernel PCA ，用 Kernel Trick 来将 PCA 推广到非线性的状况。另外，PCA 实际上能够看做是一个具备 Gaussian 先验和条件几率分布的 latent variable 模型，它假定数据的 mean 和 variance 是重要的特征，并依靠 covariance 最大化来做为优化目标，而事实上这有时候对于解决问题帮助并不大。

一个典型的问题就是作聚类或分类，回想咱们以前谈过的 Spectral Clustering ，就是使用 Laplacian eigenmap 降维以后再作 K-means 聚类，若是问题定下来了要对数据进行区分的话，"目的"就变得明朗了一些，也就是为了可以区分不一样类别的数据，再考虑直观的状况，咱们但愿若是经过降维把高维数据变换到一个二维平面上的话，能够很容易"看"出来不一样类别的数据被映射到了不一样的地方。虽然 PCA 极力下降 reconstruction error ，试图获得能够表明原始数据的 components ，可是却没法保证这些 components 是有助于区分不一样类别的。若是咱们有训练数据的类别标签，则能够用 Fisher Linear Discriminant Analysis 来处理这个问题。

同 PCA 同样，Fisher Linear Discriminant Analysis 也是一个线性映射模型，只不过它的目标函数并非 Variance 最大化，而是有针对性地使投影以后属于同一个类别的数据之间的 variance 最小化，而且同时属于不一样类别的数据之间的 variance 最大化。具体的形式和推导能够参见《Pattern Classification》这本书的第三章 Component Analysis and Discriminants 。

固然，不少时候（好比作聚类）咱们并不知道原始数据是属于哪一个类别的，此时 Linear Discriminant Analysis 就没有办法了。不过，若是咱们假设原始的数据形式就是可区分的的话，则能够经过保持这种可区分度的方式来作降维，MDS 是 PCA 以外的另外一种经典的降维方法，它降维的限制就是要保持数据之间的相对距离。实际上 MDS 甚至不要求原始数据是处在一个何种空间中的，只要给出他们之间的相对"距离"，它就能够将其映射到一个低维欧氏空间中，一般是三维或者二维，用于作 visualization 。

不过我在这里并不打算细讲 MDS ，而是想说一下在 Spectral Clustering 中用到的降维方法 Laplacian Eigenmap 。同 MDS 相似，LE 也只须要有原始数据的类似度矩阵，不过这里一般要求这个类似度矩阵具备局部性质，亦即只考虑局部领域内的类似性，若是点和距离太远的话，应该为零。有两种很直接的办法可让普通的类似度矩阵具备这种局部性质：

经过设置一个阈值，类似度在阈值如下的都直接置为零，这至关于在一个 -领域内考虑局部性。
对每一个点选取 k 个最接近的点做为邻居，与其余的点的类似性则置为零。这里须要注意的是 LE 要求类似度矩阵具备对称性，所以，咱们一般会在属于的 k 个最接近的邻居且/或反之的时候，就保留的值，不然置为零。

构造好以后再来考虑降维，从最简单的状况开始，即降到一维，经过最小化以下目标函数来实现：

从直观上来讲，这样的目标函数的意义在于：若是原来和比较接近，那么会相对比较大，这样若是映射事后和相差比较大的话，就会被权重放大，所以最小化目标函数就保证了原来相近的点在映射事后也不会彼此相差太远。

令为将的每一行加起来所获得的对角阵，而，被称做是拉普拉斯矩阵，经过求解以下的特征值问题

易知最小的那个特征值确定是 0 ，除此以外的最小的特征值所对应的特征向量就是映射后的结果。特征向量是一个 N 维列向量，将其旋转一下，正好是 N 个原始数据降到一维以后的结果。

相似地推广到 M 维的状况，只须要取除去 0 以外的最小的 M 个特征值所对应的特征向量，转置以后按行排列起来，就是降维后的结果。用 Matlab 代码写出来以下所示（使用了 knn 来构造类似度矩阵，而且没有用 heat kernel ）：

% Laplacian Eigenmap ALGORITHM (using K nearest neighbors)

% [Y] = le(X,K,dmax)

% X = data as D x N matrix (D = dimensionality, N = #points)

% K = number of neighbors

% dmax = max embedding dimensionality

% Y = embedding as dmax x N matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Y] = leigs(X,K,d)

[D,N] = size(X);

fprintf(1,'LE running on %d points in %d dimensions\n',N,D);

% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS

fprintf(1,'-->Finding %d nearest neighbours.\n',K);

X2 = sum(X.^2,1);

distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[sorted,index] = sort(distance);

neighborhood = index(2:(1+K),:);

% STEP2: Construct similarity matrix W

fprintf(1,'-->Constructing similarity matrix.\n');

W = zeros(N, N);

for ii=1:N

 W(ii, neighborhood(:, ii)) = 1;

 W(neighborhood(:, ii), ii) = 1;

end

% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF L

fprintf(1,'-->Computing embedding.\n');

D = diag(sum(W));

L = D-W;

% CALCULATION OF EMBEDDING

options.disp = 0; options.isreal = 1; options.issym = 1;

[Y,eigenvals] = eigs(L,d+1,0,options);

Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0

fprintf(1,'Done.\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

事实上，Laplacian Eigenmap 假设数据分布在一个嵌套在高维空间中的低维流形上， Laplacian Matrix 则是流形的 Laplace Beltrami operator 的一个离散近似。关于流形以及 Laplacian Eigenmap 这个模型的理论知识就不在这里作更多地展开了，下面看一个比较直观的例子 Swiss Roll 。

Swiss Roll 是一个像面包圈同样的结构，能够看做一个嵌套在三维空间中的二维流形，以下图所示，左边是 Swiss Roll ，右边是从 Swiss Roll 中随机选出来的一些点，为了标明点在流形上的位置，给它们都标上了颜色。

而下图是 Laplacian Eigenmap 和 PCA 分别对 Swiss Roll 降维的结果，能够看到 LE 成功地把这个流形展开把在流形上属于不一样位置的点映射到了不一样的地方，而 PCA 的结果则很糟糕，几种颜色都混杂在了一块儿。

另外，还有一种叫作 Locally Linear Embedding 的降维方法，它同 LE 同样采用了流形假设，并假定平滑流形在局部具备线性性质，一个点能够经过其局部邻域内的点重构出来。首先它会将下式最小化

以求解出最优的局部线性重构矩阵，对于距离较远的点和，应当等于零。这以后再把看成已知量对下式进行最小化：

这里成了变量，亦即降维后的向量，对这个式子求最小化的意义很明显了，就是要求若是原来的数据能够以矩阵里对应的系数根据其邻域内的点重构出来的话，那么降维事后的数据也应该保持这个性质。

通过一些变换以后获得的求解方法和 LE 相似，也是要求解一个特征值问题，实际上，从理论上也能够得出两者之间的联系（LE 对应于而 LLE 对应于），若是感兴趣的话，能够参考Laplacian Eigenmaps for Dimensionality Reduction and Data Representation 这篇论文里的对比。下面是两种方法在 Swiss Roll 数据上的结果，也是很是类似的：

有一点须要注意的是，LE 和 LLE 都是非线性的方法，PCA 获得的结果是一个投影矩阵，这个结果能够保存下来，在以后对任意向量进行投影，而 LE 和 LLE都是直接得出了数据降维以后的结果，可是对于新的数据，却没有获得一个"降维函数"，没有一个合适的方法获得新的数据的降维结果。因此，在人们努力寻求非线性形式扩展 PCA 的时候，另外一些人则提出了 LE 和 LLE 的线性形式，分别叫作 Locality Preserving Projection 和 Neighborhood Preserving Embedding 。

在 LPP 中，降维函数跟 PCA 中同样被定为是一个线性变换，用一个降维矩阵来表示，因而 LE 的目标函数变为

通过相似的推导，最终要求解的特征值问题以下：

获得的按照特征值从小到大排序的特征向量就组成映射矩阵，和 LE 不一样的是这里不须要去掉第一个特征向量。另外一点是在 LE 中的特征值是一个稀疏的特征值问题，在只须要求解最小的几个特征值的时候能够比较高效地求解，而这里的矩阵在乘以以后一般就再也不稀疏了，计算量会变得比较大，这个问题能够使用 Spectral Regression 的方法来解决，参见 Spectral Regression: A Unified Approach for Sparse Subspace Learning 这篇 paper 。若是采用 Kernel Trick 再把 LPP 非线性化的话，又会回到 LE 。而 LLE 的线性版本 NPE 也是采用了相似的办法来获得的，就不在这里多讲了。

另外，虽然 LE 是 unsupervised 的，可是若是训练数据确实有标签可用，也是能够加以利用的——在构造类似度矩阵的时候，属于同一类别的类似度要大一些，而不一样类别的类似度则会小一些。

固然，除去聚类或分类以外，降维自己也是一种比较通用的数据分析的方法，不过有许多人批评降维，说获得的结果没有意义，不用说非线性，就是最简单的线性降维，除去一些非藏极端的特殊状况的话，一般将原来的份量线性组合一下都不会获得什么有现成的物理意义的量了。然而话也说回来，如今的机器学习几乎都是更 prefer "黑盒子"式的方法吧，好比决策树，各个分支对应与变量的话，它的决策过程其实人是能够"看到"或者说"理解"的，可是 SVM 就不那么"直观"了，若是再加上降维处理，就更加"不透明"了。不过我以为这没什么很差的，若是只是靠能够清晰描诉出来的 rule 的话，彷佛感受神秘感不够，无法发展出"智能"来啊 ^_^bb 最后，所谓有没有物理意义，其实物理量不过也都是人为了描述问题方便而定义出来的吧。

Hierarchical Clustering
本系列不当心又拖了很久，其实正儿八经的 blog 也很久没有写了，由于比较忙嘛，不过以为 Hierarchical Clustering 这个话题我能说的东西应该很少，因此仍是先写了吧（我准备此次一个公式都不贴）。Hierarchical Clustering 正如它字面上的意思那样，是层次化的聚类，得出来的结构是一棵树，如右图所示。在前面咱们介绍过很多聚类方法，可是都是"平坦"型的聚类，然而他们还有一个更大的共同点，或者说是弱点，就是难以肯定类别数。实际上，（在某次不太正式的电话面试里）我曾被问及过这个问题，就是聚类的时候如何肯定类别数。
我能想到的方法都是比较 naive 或者比较不靠谱的，好比：
1. 根据数据的来源使用领域相关的以及一些先验的知识来进行估计——说了等于没有说啊……
2. 降维到二维平面上，而后若是数据形状比较好的话，也许能够直观地看出类别的大体数目。
3. 经过谱分析，找相邻特征值 gap 较大的地方——这个方法我只了解个大概，并且我以为"较大"这样的词也让它变得不能自动化了。
当时对方问"你还有没有什么问题"的时候我居然忘记了问他这个问题到底有没有什么更好的解决办法，过后真是至关后悔啊。不事后来在实验室里询问了一下，获得一些线索，总的来讲复杂度是比较高的，待我下次有机会再细说（先本身研究研究）。
不过言归正传，这里要说的 Hierarchical Clustering 从某种意义上来讲也算是解决了这个问题，由于在作 Clustering 的时候并不须要知道类别数，而获得的结果是一棵树，过后能够在任意的地方横切一刀，获得指定数目的 cluster ，按需取便可。
听上去很诱人，不过其实 Hierarchical Clustering 的想法很简单，主要分为两大类：agglomerative（自底向上）和 divisive（自顶向下）。首先说前者，自底向上，一开始，每一个数据点各自为一个类别，而后每一次迭代选取距离最近的两个类别，把他们合并，直到最后只剩下一个类别为止，至此一棵树构造完成。
看起来很简单吧？其实确实也是比较简单的，不过仍是有两个问题须要先说清除才行：
如何计算两个点的距离？这个一般是 problem dependent 的，通常状况下能够直接用一些比较通用的距离就能够了，好比欧氏距离等。

如何计算两个类别之间的距离？一开始全部的类别都是一个点，计算距离只是计算两个点之间的距离，可是通过后续合并以后，一个类别里就不止一个点了，那距离又要怎样算呢？到这里又有三个变种：

Single Linkage：又叫作 nearest-neighbor ，就是取两个集合中距离最近的两个点的距离做为这两个集合的距离，容易形成一种叫作 Chaining 的效果，两个 cluster 明明从"大局"上离得比较远，可是因为其中个别的点距离比较近就被合并了，而且这样合并以后 Chaining 效应会进一步扩大，最后会获得比较松散的 cluster 。
Complete Linkage：这个则彻底是 Single Linkage 的反面极端，取两个集合中距离最远的两个点的距离做为两个集合的距离。其效果也是恰好相反的，限制很是大，两个 cluster 即便已经很接近了，可是只要有不配合的点存在，就顽固到底，老死不相合并，也是不太好的办法。
Group Average：这种方法看起来相对有道理一些，也就是把两个集合中的点两两的距离所有放在一块儿求一个平均值，相对也能获得合适一点的结果。

总的来讲，通常都不太用 Single Linkage 或者 Complete Linkage 这两种过于极端的方法。整个 agglomerative hierarchical clustering 的算法就是这个样子，描述起来仍是至关简单的，不过计算起来复杂度仍是比较高的，要找出距离最近的两个点，须要一个双重循环，并且 Group Average 计算距离的时候也是一个双重循环。

另外，须要提一下的是本文一开始的那个树状结构图，它有一个专门的称呼，叫作 Dendrogram，其实就是一种二叉树，画的时候让子树的高度和它两个后代合并时相互之间的距离大小成比例，就能够获得一个相对直观的结构概览。不妨再用最开始生成的那个三个 Gaussian Distribution 的数据集来举一个例子，我采用 Group Average 的方式来计算距离，agglomerative clustering 的代码很简单，没有作什么优化，就是直接的双重循环：

def do_clustering(nodes):

 # make a copy, do not touch the original list

 nodes = nodes[:]

 while len(nodes) > 1:

 print "Clustering [%d]..." % len(nodes)

 min_distance = float('inf')

 min_pair = (-1, -1)

 for i in range(len(nodes)):

 for j in range(i+1, len(nodes)):

 distance = nodes[i].distance(nodes[j])

 if distance < min_distance:

 min_distance = distance

 min_pair = (i, j)

 i, j = min_pair

 node1 = nodes[i]

 node2 = nodes[j]

 del nodes[j] # note should del j first (j > i)

 del nodes[i]

 nodes.append(node1.merge(node2, min_distance))

 return nodes[0]

数据点又一千多个，画出来的 dendrogram 很是大，为了让结果看起来更直观一点，我把每一个叶节点用它自己的 label 来染色，而且向上合并的时候按照权重混合一下颜色，最后把图缩放一下获得这样的一个结果（点击查看原图）：

或者能够把全部叶子节点所有拉伸一下看，在右边对齐，彷佛起来更加直观一点：

从这个图上能够很直观地看出来聚类的结果，造成一个层次，并且也在整体上把上个大类分开来了。因为这里我把图横过来画了，因此在须要具体的 flat cluster 划分的时候，直观地从图上能够当作竖着划一条线，打断以后获得一片"森林"，再把每一个子树里的全部元素变成一个"扁平"的集合便可。完整的 Python 代码以下：

from scipy.linalg import norm

from PIL import Image, ImageDraw

def make_list(obj):

 if isinstance(obj, list):

 return obj

 return [obj]

class Node(object):

 def __init__(self, fea, gnd, left=None, right=None, children_dist=1):

 self.__fea = make_list(fea)

 self.__gnd = make_list(gnd)

 self.left = left

 self.right = right

 self.children_dist = children_dist

 self.depth = self.__calc_depth()

 self.height = self.__calc_height()

 def to_dendrogram(self, filename):

 height_factor = 3

 depth_factor = 20

 total_height = int(self.height*height_factor)

 total_depth = int(self.depth*depth_factor) + depth_factor

 im = Image.new('RGBA', (total_depth, total_height))

 draw = ImageDraw.Draw(im)

 self.draw_dendrogram(draw, depth_factor, total_height/2,

 depth_factor, height_factor, total_depth)

 im.save(filename)

 def draw_dendrogram(self,draw,x,y,depth_factor,height_factor,total_depth):

 if self.is_terminal():

 color_self = ((255,0,0), (0,255,0), (0,0,255))[int(self.__gnd[0])]

 draw.line((x, y, total_depth, y), fill=color_self)

 return color_self

 else:

 y1 = int(y-self.right.height*height_factor/2)

 y2 = int(y+self.left.height*height_factor/2)

 xc = int(x + self.children_dist*depth_factor)

 color_left = self.left.draw_dendrogram(draw, xc, y1, depth_factor,

 height_factor, total_depth)

 color_right = self.right.draw_dendrogram(draw, xc, y2, depth_factor,

 height_factor, total_depth)

 left_depth = self.left.depth

 right_depth = self.right.depth

 sum_depth = left_depth + right_depth

 if sum_depth == 0:

 sum_depth = 1

 left_depth = 0.5

 right_depth = 0.5

 color_self = tuple([int((a*left_depth+b*right_depth)/sum_depth)

 for a, b in zip(color_left, color_right)])

 draw.line((xc, y1, xc, y2), fill=color_self)

 draw.line((x, y, xc, y), fill=color_self)

 return color_self

 # use Group Average to calculate distance

 def distance(self, other):

 return sum([norm(x1-x2)

 for x1 in self.__fea

 for x2 in other.__fea]) \

 / (len(self.__fea)*len(other.__fea))

 def is_terminal(self):

 return self.left is None and self.right is None

 def __calc_depth(self):

 if self.is_terminal():

 return 0

 return max(self.left.depth, self.right.depth) + self.children_dist

 def __calc_height(self):

 if self.is_terminal():

 return 1

 return self.left.height + self.right.height

 def merge(self, other, distance):

 return Node(self.__fea + other.__fea,

 self.__gnd + other.__gnd,

 self, other, distance)

def do_clustering(nodes):

 # make a copy, do not touch the original list

 nodes = nodes[:]

 while len(nodes) > 1:

 print "Clustering [%d]..." % len(nodes)

 min_distance = float('inf')

 min_pair = (-1, -1)

 for i in range(len(nodes)):

 for j in range(i+1, len(nodes)):

 distance = nodes[i].distance(nodes[j])

 if distance < min_distance:

 min_distance = distance

 min_pair = (i, j)

 i, j = min_pair

 node1 = nodes[i]

 node2 = nodes[j]

 del nodes[j] # note should del j first (j > i)

 del nodes[i]

 nodes.append(node1.merge(node2, min_distance))

 return nodes[0]

agglomerative clustering 差很少就这样了，再来看 divisive clustering ，也就是自顶向下的层次聚类，这种方法并无 agglomerative clustering 这样受关注，大概由于把一个节点分割为两个并不如把两个节点结合为一个那么简单吧，一般在须要作 hierarchical clustering 但整体的 cluster 数目又不太多的时候能够考虑这种方法，这时能够分割到符合条件为止，而没必要一直分割到每一个数据点一个 cluster 。

总的来讲，divisive clustering 的每一次分割须要关注两个方面：一是选哪个 cluster 来分割；二是如何分割。关于 cluster 的选取，一般采用一些衡量松散程度的度量值来比较，例如 cluster 中距离最远的两个数据点之间的距离，或者 cluster 中全部节点相互距离的平均值等，直接选取最"松散"的一个 cluster 来进行分割。而分割的方法也有多种，好比，直接采用普通的 flat clustering 算法（例如 k-means）来进行二类聚类，不过这样的方法计算量变得很大，并且像 k-means 这样的和初值选取关系很大的算法，会致使结果不稳定。另外一种比较经常使用的分割方法以下：

待分割的 cluster 记为 G ，在 G 中取出一个到其余点的平均距离最远的点 x ，构成新 cluster H；
在 G 中选取这样的点 x' ，x' 到 G 中其余点的平均距离减去 x' 到 H 中全部点的平均距离这个差值最大，将其纳入 H 中；
重复上一个步骤，直到差值为负。
到此为止，个人 hierarchical clustering 介绍就结束了。总的来讲，在我我的看来，hierarchical clustering 算法彷佛都是描述起来很简单，计算起来很困难（计算量很大）。而且，无论是 agglomerative 仍是 divisive 实际上都是贪心算法了，也并不能保证能获得全局最优的。而获得的结果，虽说能够从直观上来获得一个比较形象的大局观，可是彷佛实际用处并不如众多 flat clustering 算法那么普遍。

Deciding the Number of Clusterings

如何肯定聚类的类别个数这个问题常常有人问我，也是一直以来让我本身也比较困惑的问题。不过说到底其实也没什么困惑的，由于这个问题自己就是一个比较 ill posed 的问题呀：给一堆离散的点，要你给出它们属于几个 cluster，这个基本上是没有惟一解或者说是没有合适的标准来衡量的。好比若是简单地用每一个类别里的点到类中心的距离之和来衡量的话，一会儿就会进入到"全部的点都独立成一类"这样的尴尬境界中。

可是 ill posed 也并非一个很好的理由，由于咱们其实大部分时候都是在处理 ill posed 的问题嘛，好比 Computer Vision 整个一个领域基本上就没有啥问题是 well posed 的…… =.=bb，好比下面盗用一下 Bill Freeman 的 Slides 中的一张讲 deblur 的问题的图：

Blurry image Sharp image Blur kernel

		=
		=
		=

从 Stata Center 的模糊照片找出清晰照片和模糊核的这个过程（特别对于计算机来讲）就是很是 ambigous 的。（顺便这个 Slides 自己也是颇有意思的，推荐看一下。）

因此么，仍是让咱们先抛开各类借口，回到问题自己。固然此次的标题取得有点宏大，其实写这篇日志的主要目的不过是想吐槽一下上一次 6.438 课的一道做业题而已……=.=bb

总而言之呢，像 Kmeans 之类的大部分经典的聚类算法，都是须要事先指定一个参数K做为类别数目的。可是不少时候这个K值并非那么好肯定的，更麻烦的是，即便你使用很暴力的方法把某个范围内的整数K所有都试一遍，一般也无法知道哪一个K才是最好的。

因此呢，咱们天然会问：那有没有算法不须要指定类别数呢？固然有！最简单的就是层次聚类，它会将你的数据点用一个树形结构给连起来，但是我想要的是聚类结果啊！这也好办，只要在合适的树深度的地方把树切开，变成一个一个的子树，每一个子树就是一个cluster。不过问题就来了，究竟什么深度才是合适的深度呢？事实上不一样的深度会产生不一样的类别数目，这基本上把原来的选择K的问题转化成了选择树深度的问题。

不要紧，咱们还有其余货，好比著名的 Mean Shift 算法，也是不须要指定类别数目就能够聚类的，并且聚类的结果也不是一棵树那种奇奇怪怪的东西。具体 Mean Shift 算法是什么样的今天就不在这里讲了，不过它其实也有一个参数 Window Size 须要选择——对，正如你们所料，这个 Window Size 其实也是会影响最终聚类的类别数目。换句话说，原来选择K的问题如今被转化为了选择 Window Size 的问题。

图 2 Factor graph demonstration for Affinity Propagation of 3 points.

而后咱们进入贝叶斯推断的领域，以前跟别人聊天的时候就据说过有个叫作 Dirichlet Process 的东西，特别神奇，可以自动地根据数据 adaptive 地肯定合适的类别数目。我感受这背后可能像其余算法同样会隐藏着其余的参数须要调节，不过鉴于我对 DP 彻底不懂，就不在这里说胡话了，感兴趣的同窗能够去研究一下。

而后贝叶斯推断聚类里的另外一个成员就是 Affinity Propagation，实际上是在讲了 Loopy Belief Propagation 算法以后的一道做业题，给了一个 factor graph，让把消息传递的公式推出来而后实现出来。看起来好像蛮容易的样子，实际上折腾了好几天，由于实现的时候碰到各类 tricky 的问题。

AP 算法把聚类问题当作一个 MAP 推断问题，假设咱们有n个数据点，这里咱们要求类别中心其实是这些数据点中的一个子集，聚类的结果能够用这个binary variable 来表示，其中表示点被归到点所表明的那个类别中。

给定个随机变量以后，接下来是要经过 local factor 来定义他们的 (unnormalized) joint distribution。首先咱们要求已知一个 Affinity Matrix S，其中表示点和点之间的类似程度（这里并不必定要求是对称的），一个简单的作法就是令

而后整个 joint distribution 主要由两类 factor 构成，在图 (fig: 2) 中给出了 3 个数据点时候的 factor graph 的例子。其中橙色的 factor 定义为

直观地来说，该 factor 的第一部分是说每一个点只能且必须被 assign 到一个 cluster；第二部分是讲被assign 到的时候 factor 的值为，也就是说会偏向于 assign 到类似度较大的那个 cluster 去。

而绿色的那些 factor 则定义为：

这样的 indicator factor 要求若是有任何点被 assign 到的话，那么必需要被 assign 到它本身。换句话说，每一个 cluster 的中心必须是属于它本身那个类的，这也是很是合理的要求。

好了，模型的定义就是这么简单，这样咱们会获得一个 loopy factor graph，接下来只要在这个 graph 上跑一下 Max-Product Algorithm 就能够获得最优的，从而获得聚类结果了，而且都不用指定类别数目什么的，由于合适的类别数目能够从数据中自动地推断出来！

简直太神奇了，我都不敢相信这是真的！结果，果真童话里都是骗人的……由于隐藏在背后的还有许多细节须要处理。首先有一些 trick 来下降算法的复杂程度，好比说，因为全部的变量都是 binary 的，所以在消息传递的过程当中咱们能够只传递等于 1 和等于 0 两种状态时的消息差；而后因为数值下溢的缘由，通常都会对全部的消息取log ，这样 Max-Product 算法会变成 Max-Sum（或者 Min-Sum）；而后有些消息是不作任何操做就直接 literally 传到接下来的节点的，这样的消息能够从中间步骤中去掉，最后化简以后的结果会只留下两种消息，在经典的 Affinity Propagation 算法中分别被成为 responsibility 消息和 availability 消息，能够有另外的 intuitive 的解释，具体能够参见 (Frey & Dueck, 2007)。

图 3 Clustering result of Affinity Propagation.

不过这些仍是不太够，由于 Loopy BP 的收敛性实际上是没有什么保证的，随便实现出来极可能就不收敛。好比说，若是直接进行 parallel 地 message updating 的话，极可能就因为 update 的幅度过大致使收敛困难，所以须要用 dampen 的方法将 fixed-point 的式子

改成

而后用这个来作 updating。固然也能够用 0.5 之外的其余 dampen factor。不过这样仍是不够的，反正我折腾了很久最后发现必须把 parallel updating 改为 in-place updating 才能收敛，这样一来更新的幅度就更小了一些，不过 dampen 仍然是须要的。另外就是一同讨论的另外的同窗还发现实现中消息更新的顺序也是极端重要……把顺序变一下就立马从彻底不 work 变成结果超好了。-_-!!因而忽然想起来 CV 课上老师拿 Convolutional Network 开玩笑的时候说，这个算法效果超级好，可是就是太复杂了，在发明这个算法的那个 lab 所在地的方圆 50 米以外彻底没人能成功训练 Convolutional Network……

图 4 Clustering result of Affinity Propagation, with large diagonal values for the affinity matrix.

此外收敛性的判断也是有各类方法，好比消息之差的绝对值、或者是连续多少轮迭代产生的聚类中心都没有发生变化等等。即便在收敛以后，如何肯定 assignment 也并非显而易见的。由于 Max-Product 算法获得的是一些 Max-Marginal，若是全局最优的 configuration 不是 unique 的话，general 地是不能直接局部地经过各自的 Max-Marginal 来肯定全局 configuration 的，因而可能须要实现 back-tracking，总之就是动态规划的东西。不过也有简单的 heuristic 方法能够用，好比直接经过的Max-Marginal 来肯定点究竟是不是一个 cluster center，肯定了 center 以后则再能够经过 message 或者直接根据 Affinity 矩阵来肯定最终的类别 Assignment。图 (fig: 3) 展现了一个二维状况下 100 个点经过 AP 聚类的结果。

不过，你有没有发现一点不对劲的地方？一切彷佛太完美了：什么参数都不用给，就能自动肯定类别数？果真仍是有问题的。事实上，若是你直接就去跑这个算法，确定得不到图上的结果，而是会获得一个n个类别的聚类结果：全部的点都成为了中心而且被归类到了本身。这样的结果也是能够理解的，由于咱们计算 Affinity 的方式（距离的相反数）致使每一个点和本身的 Affinity 是最大的……因此呢，为了解决这个问题，咱们须要调整一个参数，就是S矩阵的对角线上的值，这个值反应了每一个点本身成为类别中心的偏好程度，对的，如你们所料，这个就是背后的隐藏参数，直接影响类别数目。

图 5 Clustering result of Affinity Propagation, with small diagonal values for the affinity matrix.

如图 (fig: 4) 和 (fig: 5) 分别是把对角线上的元素值设置得比较大和比较小的状况，能够看到前者产生了更多的类别，由于你们都比较倾向于让本身成为类别中心；然后者则只产生了两个类，由于你们都很不情愿……（我发现这种画一些线连到类别中心的 visualization 方法彷佛会给人产生很强的暗示感受，因而各类聚类结果都看起来好像很"对"的样子）

因此呀，到最后好像仍是竹篮打水一场空的样子，怎么绕怎么绕好像也绕不开类别数目的这个参数。不过好像也并非彻底没有任何进展的，由于虽然各类算法都或多或少地须要指定一些参数，可是这些参数从不一样的角度来阐释对应的问题，好比说在 AP 中，通常能够直接将的对角线元素设置为其余因此元素的中位数，一般就能获得比较合理的结果；有时候对于特定的实际问题来讲，从某一个角度来进行参数选择可能会有比较直观的 heuristic 能够用呢。因此说在是否须要指定参数这个问题上，仍是不用太钻牛角尖了啊。

最后我还想提一下，最近看到的 lab 的两篇 NIPS 文章 (Canas, Poggio & Rosasco, 2012),(Canas & Rosasco, 2012)，里面的视角也是颇有意思的。

咱们知道，好比说，KMeans 算法的目标函数是这样定义的：

其中是一个由K个点组成的集合，而一个点到的距离被定义为到S中全部点的距离中最小的那一个值。这样能够当作是用的这K个点去近似本来的n个点所带来的偏差。然而这样的目标函数对于选择合适的K并无什么指导意义，由于随着K的增大咱们获得的最优的会使得目标函数愈来愈小。

可是咱们能够将这个目标函数推广一下，相似于 Supervised Learning Theory 中那样引入几率空间。具体地来讲，咱们假设数据点是从一个未知的几率分布中采样而获得的，而且目标函数也从对n个数据点的近似推广到整个生成流形的近似：

固然因为是未知的，因此是无法计算的，可是咱们能够从 data sample 获得的来对它进行近似，并获得 bound 之类的。直观地来说，此时对数值K的选择将会影响到模型空间的复杂度，因而会出现一个trade-off，因而从这个角度下去探讨"最优的K值"就变成一个很合理的问题了。感兴趣的同窗能够详细参考 paper 以及里面的参考文献，都在 arXiv 上能够下载到的。

References

Canas, G. D., & Rosasco, L. (2012). Learning Probability Measures with respect to Optimal Transport Metrics. In NIPS.
Canas, G. D., Poggio, T., & Rosasco, L. (2012). Learning Manifolds with K-Means and K-Flats. In NIPS.
Frey, B. J., & Dueck, D. (2007). Clustering by passing messages between data points.Science, 315, 972–977.

1. 漫谈 Clustering (1): k-means（转）
2. [Scikit-learn] 2.3 Clustering - Spectral clustering
3. 漫谈慢慢谈。。。
4. 漫谈并发
5. 漫谈Puppet4
6. 漫谈JSTL
7. 漫谈集成
8. 漫谈HBase Filter
9. 漫谈gRPC
10. 并发漫谈
更多相关文章...
• NoSQL和传统数据库有什么区别？NoSQL能取代传统数据库吗？ - Redis教程
• MyBatis SqlSessionFactory - MyBatis教程
• 漫谈MySQL的锁机制
• Docker容器实战(八) - 漫谈 Kubernetes 的本质