爱华网漫谈 Clustering (5): Hierarchical Clustering
本文是“漫谈 Clustering 系列”中的第 8 篇，参见本系列的其他文章。
系列不小心又拖了好久，其实正儿八经的 blog 也好久没有写了，因为比较忙嘛，不过觉得 HierarchicalClustering 这个话题我能说的东西应该不多，所以还是先写了吧（我准备这次一个公式都不贴 ）。Hierarchical Clustering正如它字面上的意思那样，是层次化的聚类，得出来的结构是一棵树，如右图所示。在前面我们介绍过不少聚类方法，但是都是“平坦”型的聚类，然而他们还有一个更大的共同点，或者说是弱点，就是难以确定类别数。实际上，（在某次不太正式的电话面试里）我曾被问及过这个问题，就是聚类的时候如何确定类别数。
我能想到的方法都是比较 naive 或者比较不靠谱的，比如：
当时对方问“你还有没有什么问题”的时候我竟然忘记了问他这个问题到底有没有什么更好的解决办法，事后真是相当后悔啊。不过后来在实验室里询问了一下，得到一些线索，总的来说复杂度是比较高的，待我下次有机会再细说（先自己研究研究）。
不过言归正传，这里要说的 Hierarchical Clustering 从某种意义上来说也算是解决了这个问题，因为在做Clustering 的时候并不需要知道类别数，而得到的结果是一棵树，事后可以在任意的地方横切一刀，得到指定数目的 cluster，按需取即可。
听上去很诱人，不过其实 Hierarchical Clustering的想法很简单，主要分为两大类：agglomerative（自底向上）和divisive（自顶向下）。首先说前者，自底向上，一开始，每个数据点各自为一个类别，然后每一次迭代选取距离最近的两个类别，把他们合并，直到最后只剩下一个类别为止，至此一棵树构造完成。
看起来很简单吧？ 其实确实也是比较简单的，不过还是有两个问题需要先说清除才行：
- 如何计算两个点的距离？这个通常是 problem dependent的，一般情况下可以直接用一些比较通用的距离就可以了，比如欧氏距离等。
- 如何计算两个类别之间的距离？一开始所有的类别都是一个点，计算距离只是 计算两个点之间的距离，但是经过后续合并之后，一个类别里就不止一个点了，那距离又要怎样算呢？到这里又有三个变种：
总的来说，一般都不太用 Single Linkage 或者 Complete Linkage 这两种过于极端的方法。整个agglomerative hierarchical clustering的算法就是这个样子，描述起来还是相当简单的，不过计算起来复杂度还是比较高的，要找出距离最近的两个点，需要一个双重循环，而且 GroupAverage 计算距离的时候也是一个双重循环。
另外，需要提一下的是本文一开始的那个树状结构图，它有一个专门的称呼，叫做 Dendrogram，其实就是一种二叉树，画的时候让子树的高度和它两个后代合并时相互之间的距离大小成比例，就可以得到一个相对直观的结构概览。不妨再用最开始生成的那个三个Gaussian Distribution 的数据集来举一个例子，我采用 Group Average的方式来计算距离，agglomerative clustering 的代码很简单，没有做什么优化，就是直接的双重循环：
123456789101112131415161718192021 | 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 代码如下：
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106 | from scipy.linalg import normfrom PIL import Image, ImageDrawdef 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这样受关注，大概因为把一个节点分割为两个并不如把两个节点结合为一个那么简单吧，通常在需要做 hierarchicalclustering 但总体的 cluster数目又不太多的时候可以考虑这种方法，这时可以分割到符合条件为止，而不必一直分割到每个数据点一个 cluster 。
总的来说，divisive clustering 的每一次分割需要关注两个方面：一是选哪一个 cluster来分割；二是如何分割。关于 cluster 的选取，通常采用一些衡量松散程度的度量值来比较，例如 cluster中距离最远的两个数据点之间的距离，或者 cluster 中所有节点相互距离的平均值等，直接选取最“松散”的一个 cluster来进行分割。而分割的方法也有多种，比如，直接采用普通的 flat clustering 算法（例如k-means）来进行二类聚类，不过这样的方法计算量变得很大，而且像 k-means这样的和初值选取关系很大的算法，会导致结果不稳定。另一种比较常用的分割方法如下：
到此为止，我的 hierarchical clustering 介绍就结束了。总的来说，在我个人看来，hierarchicalclustering 算法似乎都是描述起来很简单，计算起来很困难（计算量很大）。并且，不管是 agglomerative 还是divisive实际上都是贪心算法了，也并不能保证能得到全局最优的。而得到的结果，虽然说可以从直观上来得到一个比较形象的大局观，但是似乎实际用处并不如众多flat clustering 算法那么广泛。
