杂货店-聚类方法综述

谱聚给你博客园上若干个博客，让你将它们分成K类，你会怎样做？想必有很多方法，本文要聚类的直观解释是根据样本间相似度，将它们分成不同组。谱聚类的思想是将样本看作顶接不同组的边的权重尽可能低（这意味着组间相似度要尽可能低（这意味着组内相似度要尽可能高如下图所示：在上图中，一共有6个顶点（博客），顶点之间的连线表示两个顶点的相似度，现在要将这图分成两半（两个类（去掉哪边条？根据谱聚类的思想，应该去掉的边是用虚线表示的那条。最后，剩下的两半就分别对应两个类了。根据这个思想，可以得到unnormalized谱聚类和normalized谱聚类，由于前者比后者简单，所以本文介绍unnormalized谱聚类的几个步骤（假设要分K个类）：unnormalizedgraphLaplacianmatrixL(LDW,Ddegree计算L的前K把这k个特征向量排列在一起组成一个N*k的矩阵，将其中每一行看作k中的一个向量，并使用K-means算法进行聚类这一节主要从大体上解unnormalized谱聚类的四个步骤是怎么来的，不涉及具体的公谱聚类的思想就是要转化为图分割问题。因此，第一步就是将原问题转化为图。转为图有两个问题要解决：一是两个顶点的边要怎样定义；二是要保留哪些边。对于第一个问题，如果两个点在一定程度上相似，就在两个点之间添加一条边。相似的过常用的Gaussiansimilarityfunction法是建立k-nearestneighborgraph。在这种图中，每个顶点只与K个相似度最高的点连边。unnormalizedgraphLaplacianmatrix(L表示)有很多很好的性质，也正是这个这一组性质将在之后的公式推导中起到决定性作用将原问题转化为图后，接下来的工作就是决定怎样分。图分割问题实际上就是最小割问题(mincutpoblem)。最小割问题可定义为最小化以下目标函数：其中k表示分成k个组，Ai表示第i个组，表示第Ai的补集，W(A,B)表示第A组与第组之间的所有边的权重之和这个式子的直观意义：如果要分成K个组，那么其代价就是进行分割时去掉的边的权重的总和。可惜的是直接最小化这式子通常会导致不好的分割。以分成2类为例，这个式子通常会将图分成这样的两类：一个点为一类，剩下的所有点为另一类。显然，这样的分割是很不好的。因为我们期望着每个类都有合理的大小。所以，要对这个式子进行改进，改进后的公式（称为RatioCut)如下其中|A|表示A组中包含的顶点数目在RatioCut中，如果某一组包含的顶点数越少，那么它的值就越大。在一个最小化问题中，这相当于是惩罚，也就是不鼓励将组分得太小。现在只要将最小化RatioCut解出来，分割就完成了。不幸的是，这是个NP难问题。想要在多项式时间内解出来，就要对这个问题作一个转化了。在转化的过程中，就用到上面提到的L的那一组性质，经过若干推导，最后可以得到这样的一个其中H是一个矩阵，它的元素的定义(Eq.(5))如果H矩阵的元素不为0，则说明第i个点属于第j个类。也就是说，只要得到H矩阵，就能知道要怎样分。可惜的是，这个问题仍然是NP难问题。但是，如果我们让H矩根据Rayleigh-Ritztheorem，这个问题的解是L的前k个最小的特征向量组成的矩阵H，其中特征向量是按列来排H的每一列，均为一个特征向量。在第三步中，我们为了松驰NP难问题，让H矩阵取任意值，因此，解出来的H矩阵不再具有原来的性质——元素值能哪个点属于哪一类。尽管如此，对于k-means来说，将H矩k-meansH矩阵进行聚类作为谱以下是unnormalied谱聚类的 版实现（博客园的代码格式选择中居然没有的。。。这里选个C++的）：functionfunction[C,L,D,Q,V]=SpectralClustering(W,%spectralclustering%input:adjacencymatrixW;numberofcluster%%return:clusterindicatorvectorsascolumnsinC;unnormalizedLaplacianL;degreematrixD; eigenvectorsmatrixQ;eigenvaluesmatrix%calculatedegreematrixdegs=sum(W,2);D=sparse(1:size(W,1),1:size(W,2),%computeunnormalizedLaplacianL=D-W;%computetheeigenvectorscorrespondingtotheksmallest%diagonalmatrixVisNcutL'sksmallestmagnitude%matrixQwhosecolumnsarethecorrespondingeigenvectors.[Q,V]=eigs(L,k,'SA');%usethek-meansalgorithmtoclusterVrow-%Cwillbean-by-1matrixcontainingtheclusternumberforeachdatapointC=kmeans(Q,k);%convertCtoan-by-kmatrixcontainingthekindicatorvectorsascolumnsC=sparse(1:size(D,1),C,1);[1]ATutorialonSpectralClustering[2]ClusteringClustering1k-本文是“漫谈Clustering系列”中的第1篇，参列的其他文章。好久没有写blog了，一来是blog下线一段时间，而租DreamHost的事情又Learning相关的东西，因为很多东西都不懂，所以平时也找一些资料来看。按Clustering中文翻译作“聚类”，简单地说就是把相似的东西分到一组，同Classification(分类)不同，对于一个classifier，通常需要你告诉它“这classifier练数据的过程通常叫做supervisedlearning(监督学习)，而在聚类的时候，因此，一个聚类算法通常只需要知道如何计算相似度就可以开始工作了，因此clustering通常并不需要使用训练数据进行学习，这在MachineLearning中被称作unsupervisedlearning(无监督学习)。映射到一个单独的特征之外，一种常见的做法是同时提取N种特征，将它们放NNclustering大致看出来”，不过，对于这样的N维欧氏空间中的点进行聚类，有一个非常k-meansk-meansclustercenter)cluster所有的点到该中心点的距离小于到其他cluster的中心的距离。虽然实际情况情况是2，通常我们会将它归类为左边的那个分布，因为概率大一些，然而此分性，无法避免。因此，k-means所依赖的这个假设看作是合理的。有N个数据点需要分为K个cluster，k-means要做的就是最小化 在数据点n被归类到clusterk的时候为1，否则0。直接寻找 和来最小化 代的办法：先固定，选择最优的 的。将 对求导并令导数等于零，很容易得到 候应该满足：亦即的值应当是所有clusterk中的数据点的平均值。由于每一次迭代都是取到的最小值，因此只会不断地减小（或者不变），而不会增加，k-meansk-means到全局最优解，但是对于这样的问题，像k-meansk-means选定K个中心的初值。这个过程通常是针对具体的问题有一些启发式的，或者大多数情况下采用随机选取的办法。因为前面k-means并以有时候我们会多次选取初值跑k-means，并取其中最好的一次结果。cluster按照这个步骤写一个k-means实现其实相当容易了，在SciPy或者k-meansk-means123456789123456789fromfuture importcPickleaspicklefrommatplotlibimportfromnumpyimportzeros,array,tilefromscipy.linalgimportnormimportnumpy.matlibasmldefkmeans(X,k,observer=None,threshold=1e-15,maxiter=300):N=len(X)labels=zeros(N,centers=array(random.sample(X,k))iter=0defsum=foriinreturndefdistmat(X,Y):n=len(X)m=xx=ml.sum(X*X,axis=1)yy=ml.sum(Y*Y,axis=1)xy=ml.dot(X,Y.T)returntile(xx,(m,1)).T+tile(yy,(n,1))-Jprev=calc_J()whileTrue:#notifytheifobserverisnotNone:#calculatedistancefromxtoeach#distance_matrixisonlyavailableinscipynewerthan0.7#dist=distance_matrix(X,centers)dist=distmat(X,centers)#assignxtonearst#re-calculateeachcenterforjinrange(k):idx_j=(labels==j).nonzero()iter+=1ifJprev-J<threshold:Jprev=#finalifobserverisnotNone:if =='main#loadpreviouslygeneratedpointswithopen('cluster.pkl')asinf:samples=N=forsmpinX=zeros((N,2))idxfrm=0foriinidxto=idxfrm+len(samples[i][0])X[idxfrm:idxto,0]=samples[i][0]X[idxfrm:idxto,1]=samples[i][1]idxfrm=idxtodefobserver(iter,labels,centers):print"iter%d."%itercolors=array([[1,0,0],[0,1,0],[0,0, #clearpreviousplot#drawdata_colors=[colors[lbl]forlblinpyplot.scatter(X[:,0],X[:,1],c=data_colors,alpha=0.5)#drawcenterskmeans(X,3,代码有些长，不过因为用Python来做这个事情确实不如方便，实际的k-means的代码只是4147341到4345473个中心点，结果不过正如前面所说的那样k-meansk-means都还是很令人满意的，算是一种简单高效应用广泛的clustering方法。Update2010.04.25:很多人都问我要cluster.pkl，脆把它上传上来吧，其实是很容易自己生成的，点击这里。漫谈Clustering(2):k-me本文是“漫谈Clustering系列”中的第2篇，参见本系列的其clusteringk-means事实也确实如此，k-meds可以算是k-means的一个变种。将中心点取为当前cluster中所有数据点的平均值：Rough并且我们已经证明在固定了各个数据点的assignment的情况下这样选取的中 最小化。然而在k-meds中，中心点的选取限制在当前cluster所包含的数据点的集合中。换句话说，在k-meds算法中，从当前cluster中选取这样一个点——它到其他所有（当前cluster中的）点的距离之和最小——作为中心点。k-means和k-meds之间的差异就类似于一个数据样本的均值(mean)和中位数(median)之间的差k-means(dissimilarity)可以很自然地用这样的方式来处理，但是类别(categorical)类型的特征就不组成的空间中进行，k-means就为力了，因为欧氏距离在里不能用了：一只Samoyed减去一只RoughCollie然后在平方一下？天知道那是什么！再加上一只GermanShepherdDog然后求一下平均值？根本没法算，k-means在这里寸步难行！在k-meds中，我们把原来的目标函数中的欧氏距离改为一个任意的dissimilaritymeasure函数最常见的方式是构造一个dissimilaritymatrix 来代表，其中的元素表示第 只狗和第只狗之间的差异程度，例如，两只Samoyed 之间的差异可以设为0，一只GermanShepherdDog和一只RoughCollie之间的差异是0.7，和一只MiniatureSchnauzer之间的差异是1，等等。除此之外，由于中心点是在已有的数据点里面选取的，因此相对于k-means来说，不容易受到那些由于误差之类的原因产生的Outlier robust一些。forjinidx_j=(labels==j).nonzero()distsum=ml.sum(distj,axis=1)icenter=distsum.argmin()centers[j]=扯了这么forjinidx_j=(labels==j).nonzero()distsum=ml.sum(distj,axis=1)icenter=distsum.argmin()centers[j]=可以看到k-meds在这个例子上也能得到很好的结果k-meansk-meansk-meds复杂度陡然增加了许多：在k-means即可，而在k-meds中则需要枚举每个点，并求出它到所有其他点的距离之和，复杂度为看完了直观的例子，让我们再来看一个稍微实际一点的例子好了：Clustering——那个永恒不变的，不过我们这里要做的聚类并不是针对文档的而是针对文档的语言实验数据是从Europarl 和Swedish这些语言的文本数据集。在N-gram-basedtextcategorization这篇paper中描述了一种计算由不同语言写成的文档的相似度的方法。一个（以字符为单位的）N-gram就相当于长度为N的一系列连续子串。例如，由o产生的3-gram为：hel、ell和lloN-gram之前在开头和末尾加上空格（这里用下划线表示）：_he、hel、ell、llo、lo_和o。按照Zipf’slaw：Thenthmostcommonwordinahumanlanguagetextoccurswithafrequencyinverselyproportionalton.N-gramword国家的那些类英语的语言写成的文档，不论或长短，通常得出来的Profile要很长）文档构建出一个Profile，在拿到一篇未知文档的时候，只要和各个Profile比较一下，差异最小的那个Profile所对应的语言就可以认定是这篇123456789123456789classdefinit(self,path,N=3,psize=400):self.N=Nself.psize=sepgrams={}withopen(path)asinf:forlineininf:fortokinself.sep.split(line):forninrange(self.N):#keeponlythetopmostpsizeitemsgrams=file=foriinrange(len(grams)):defgetitem(self,ifidxisNone:returnreturndis=0file.keys():returndeffeed_ngram(self,grams,tok,n):ifn!=0:tok='_'+toktok=tok+'_'*(n-1)foriinrange(len(tok)-n+1):gram=tok[i:i+n]grams[gram]+=europarl数据集共有11600我为这七千多个文档建立了Profile并构造出一个7038×7038的dissimilaritymatrix，最后在这上面用k-meds进行聚类。构造的，并且通常只要数次迭代就能收敛了。实际的k-meds实现可以在mltk 并不知道这些cluster应该被打上什么，或者说，在这个问题里，包含了哪种语言的文档。由于我们的目的是衡量聚类算法的performance，因此直接假定这一步能实现最优的对应关系，将每个cluster对应到一种语言上去。一以用Hungarianalgorithm来求解。1111!=我们需要遍历一次labellist，并数出真正的label（语言）与聚类得出的并求出accuracy只需要1毫秒的时间的话，总共也需要11个小时才能得到具体的算法了，感的同学可以参考Wikipedia，这里我直接使用了一个现有的Python实现。虽然这个实验非常折腾，不过最后的结果其实就是一个数字：accuracy我这里达到了88.97%，证明k-meds聚类和N-gramProfile识别语言这要版的scipy，munkres.py和mltk以及Python2.6。Clustering番外篇Vector本文是“漫谈Clustering系列”中的第3篇，参见本系列的其Vectorzation。这项技术广泛地用在信号处理以及数据压缩等领域。事实上，在JPEG和MPEG-4等多压缩格式里都有VQ这一步。Vectorzation这个名字听起来有些玄乎，其实它本身并没有这么高深。比如，[0,1)上的所有值变为0，[1,2)上的所有值变成1，如此类推。其这就是一个VQ的过程。一个比较正式一点的定义是：VQ是将一个向量空间中一个典型的例子就是图像的编码最简单的情况考虑一个灰度 为黑色1为白色，每个像素的值为[0,1]上的一个实数。现在要把它编码为256阶的灰阶，一个最简单的做法就是将每一个像素值x 数floor(x*255) 现在想要把压缩这个，每个像素只使用4bit（而不是原来的8bit），因此，要将原来的[0,255]区间上的整数值用[0,15]上的整数值来进行编码，一个简单的映射方案是x*15/255 不过这样的映射方案颇有些Naive，虽然能减少颜色数量起到压缩的效果，但例如，如果一个256阶灰阶完全由0和13两种颜色组成，那么通过上面的映射就会得到一个全黑的，因为两个颜色全都被映射到0了。一个更好实际做法就是：将每个像素点当作一个数据，跑一下K-means，得到k个点的像素值。对于彩片来说，也可以用同样的方法来做，例如RGB三色，每一个像素被当作是一个3用本文开头那张RechardStallman大神的来做一下实验好了，VQ2、VQ10和VQ100三张分别显示聚类数目为2、10和100时得到的结果，可VQ100centroids欢的。考虑一种最简单的压缩办法：单独（比如100个）centroids的颜色信息，然后每个像素点centroid的索引而不是颜色信息值，如果一个VQ实现代码很简单直接使用了SciPy 提供的kmeans 写用了PythonImageLibrary：123123456789fromscipy.cluster.vqimportkmeans,vqfromnumpyimportarray,reshape,zerosfrommltkimportimagevqclst=[2,10,100,data=image.read('example.jpg')(height,width,channel)=data.shapeforkinprint'Generatingvq-%d...'%k(centroids,distor)=kmeans(data,k)(code,distor)=vq(data,centroids)print'distor:%.6f'%distor.sum()im_vq=centroids[code,:](height,width,当然，Vectorzation并不一定要用K-means来做，各种能用的聚类方法都可以用，只是K-means通常是最简单的，而且通常都够用了。Clustering3GaussianMixture本文是“漫谈Clustering系列”中的第4篇，参列的其他文章。k-means行的算法：GaussianMixtureModel(GMM)。事实上，GMMk-meansGMM（GMMclusteringdensityestimation），简单地说，k-means数据点被assign到其中某一个cluster了，而GMM则给出这些数据点被assign每个cluster概率，又称作softassignment。以把这个概率转换为一个score情况，比如，49%51%50%废话说了一堆，不过，在回到GMM之前，我们再稍微扯几句。我们知道，不管如果去掉“线性函数”这个归纳偏执，因为对于N个点，我们总是可以构造一个N-1次多项式函数，让它完美地穿过所有的这N个点，或者如果我用任何大大的N，从而得到一个（或者无穷多个）“超级函数”，能够fit这个领域内所有的问题。然而这个（或者这无穷多个）“超级函数”有用吗？只要我们注意没有归纳偏执或者归纳偏执太宽泛会导致Overfitting，然而另一个──望的目的，然而绝大多数模型都会有那么一些“参数”（例如K-means中的GMM，GMMGaussianMixtureModel就是假设数据服从MixtureGaussianDistribution，换句话说，数据可以看作是从数个GaussianDistribution中生成出来的。实际上我们在K-means和 K-meds两篇文章中用到的那个例子就是由三个Gaussian分布从随机选取出来的。实际上，从中心极限定理可以看出，Gaussian分布（也叫做正态(Normal)分布）这个假设其实是比较合理的，除此之外，Gaussian分布在计算上也有一些很好的性质，所以，虽然我们可以用不同的分布来随意地构造XXMixtureModel，但是还是GMM最为流行。另外，MixtureModel本身其实也是可以变得任意复杂的，通过增加Model的个数，每个GMM由 个Gaussian分布组成，每个Gaussian称为一个“Component”，这些Component线性加成在一起就组成了GMM的概率密度函根据上面的式子，如果我们要从GMM以分为两步：首先随机地在这个ComponentComponent被选中的概率实际上就是它的系数，选中了Component之后，再单独地考虑从这个Component的Gaussian分布，转化为了已知的问题。那么如何用GMM来做clustering呢？其实很简单，现在我们有了数据，假定它们是由GMM生成出来的那么我们只要根据数据推出GMM的概率分布来就可以了然后GMM的 个Component实际上就对应了 个cluster了。根据数据来推算概率密度通常被称作densityestimation，特别地，当我们在 个数据点，并假设它们服从某个分布（记作），现在要确定里面的一些参数的值，例如，在GMM中，我们就需要确定 和这些参数。我们的想法是，找到这样一组参数，它所确定的概率于，我们把这个乘积称作似然函数(LikelihoodFunction)。下溢，因此我们通常会对其取对数，把乘积变成加 ，得下面让我们来看一看GMM的log-likelihoodfunction：值。为了解决这个问题，我们采取之前从GMM实际上也就类似于K-means的两步。Component（Component率）：对于每个数据来说，它由第个Component生成的概率为由于式子里的和也是需要我们估计的值，我们采用迭代法，在计算的时候我们假定和均已知取上一次迭代所得估计每个Component的参数：现在我们假设上一步中得到的就是正确的“数据由Component生成的概率”，亦可以当做该Component在生成这个数据上所做的贡献，或者说，我们可以看作这个值其中有这部分是由Component 现在实际上可以看作Component生成了这些点。由于每个Component都是一个标准的Gaussian分布，可以很容易分布求出最其中，并且也顺理成章地可以估计为考PatternRecognitionandMachineLearning这本书的第九章。有了实际的步骤，再实现起来就很简单了。代码如下：矩阵singular的情况，可以参见这篇文章。123123456789K=%randomlypickcentroidsrndp=randperm(N);centroids=X(rndp(1:K),:);K=size(K_or_centroids,1);centroids=K_or_centroids;threshold=1e-[N,D]=%PX:N-by-Kmatrixindicatingtheprobabilityofeachcomponentgeneratingeachpoint.MODEL:astructurecontainingtheparametersforaGMM:MODEL.Miu:aK-by-Dmatrix.MODEL.Pi:a1-by-KX:N-by-DdataK_OR_CENTROIDS:eitherKindicatingthenumberofcomponentsoraK-by-DmatrixindicatingthechoosingoftheinitialK%%%%%%%%%%%%GaussianMixture%%PX=GMM(X,%[PXMODEL]=GMM(X,%%functionvarargout=gmm(X,%%initial[pMiupPipSigma]=whiletruePx=%newvalueforpGammapGamma=Px.*repmat(pPi,N,1);pGamma=pGamma./repmat(sum(pGamma,2),1,%newvalueforparametersofeachComponentNk=sum(pGamma,1);pMiu=diag(1./Nk)*pGamma'*X;pPi=Nk/N;forkk=Xshift=X-repmat(pMiu(kk,:),N,pSigma(:,:,kk)=(Xshift'*(diag(pGamma(:,kk))*Xshift))/%checkforconvergenceL=sum(log(Px*pPi'));ifL-Lprev<thresholdLprev=ifnargout==varargout=model=[];model.Sigma=pSigma;model.Pi=pPi;varargout={Px,model};function[pMiupPipSigma]=init_params()pMiu=centroids;pPi=zeros(1,K);pSigma=zeros(D,D,K);%hardassignxtoeachdistmat=repmat(sum(X.*X,2),1,K)+repmat(sum(pMiu.*pMiu,2)',N,1)-...[dummylabels]=min(distmat,[],forXk=X(labels==k,:);pPi(k)=size(Xk,1)/N;pSigma(:,:,k)=cov(Xk);functionPx=calc_prob()Px=zeros(N,K);fork=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)*Px(:,k)=coef*exp(-函数返回的Px 行中最大的那个概率值所对应的那个Component为 K-means那个cluster有一些点跑得比较远了。当然，因为这个问题原本就是完全有MixtureGaussianDistribution，GMM（如果能求得全局最优解的GMMK-means似（都可以追溯到EM）K-means值，就有可能得到很差的结果。对于K-meansGMMK-means个更流行的做法是先用K-means（已经重复并取最优值了）得到一个粗略的结果，然后将其作为初值（只要将K-means所得的centroids传入gmm函数即可），再用GMM进行细致迭代。如我们最开始所讨论的，GMM所得的结果（Px）不仅仅是数据点的label，而Clustering番外篇Expectation本文是“漫谈Clustering系列”中的第5篇，参列的其他文章。Expectationization(EM)是一种以迭代的方式来解决一类特殊最大似然(umLikelihood)问题的方法，这类问题通常是无法直接求得最优解，我们会看到，上一次说到的GaussianMixtureModel的迭代求解方法可以算是EM算法最典型的应用，而最开始说的K-means其实也可以看作是GaussianMixtureModel的一个变种（固定所有的，并令即可）。然而EM实际上是一种通用的算法（或者说是框架），可GaussianMixtureModel起吧。回顾一下我们之前要解决的问题：求以下Log-likelihoodfunction的 函数里面又有加和没法直接求考虑一下GMM生成sampleGaussianDistribution里进行sample。我们可以很自然地引入一个隐含变 个Component被选中了，那么 个元素置为1，其余的全为0。那么，再来看看，如果除了之前的sample 为。注意到 个元素为1的时候，亦即 Component被选中的时候，这个概率 前是一个和式再替换到最开始的那个Log-likelihoodfunction中得到新的同时关于sample和隐含变量的Log-likelihood:情况瞬间逆转，现在和求和符号换了个位置，直接作用在普通的高斯分布利，完全依赖于一个我们的假设：隐含变量的值已知。然而实际上我们并不0sample的值因此我们可以把这个信息利用起来针对后验概率来取期望。前面，的每一个元素只有0和1两种取值，因此按照期望的中间用贝叶斯定理变了一下形，最后得到的式子正是我们在漫谈GMM中定义的。因此，对于上面那个可以直接求解的Log-likelihoodfunction ，我们只要用其期望亦即代替即可。到这里为止看起来是一个非常完美的方法不过仔细一看发现还有一个：之前的Log-likelihoodfunction可以直接求最大值是建立 是已知况下，现在虽然我们用来代替了 ，但是实际上是一个GMM的方法又推导了一遍。EM现在我们的讨论将不局限于GMM并使用一些稍微紧凑一点的符号用 示所有的sample，同时用 过最大似然的方法估计出中的参数 很，但是要优化却很容易。这就是EM算法能够解决的那一现在我们引入一个关于隐含变量的分布。注意到对于任意的，我们都可以对Log-likelihoodFunction作如下分解：和之间的Kullback-divergence。由于Kullback-Leiblerdivergence个分布完全相同的时候才会取到0。因此我们可以得到关系，亦即是的一个下界。现在考虑EM的迭代过程，记上一次迭代得出的参数为，现在我们要取以令最大，由于并不依赖于因此的上限（在固定的时候）是一个定值，它取到这个最大值的条件就是Kullback-Leiblerdivergence为零，亦即等于后验概率。把它带入到的表达式中可以得到其中const是常量，而是我们之前所得到的“同时包含sampleLog-likelihoodfunction这个对应到EM中的“E”一步。在接下来的“M”一步中，我们讲固定住分布，再选取合适的以将最大化，这次其上界也依赖于变量，并会随着的增大而增大（因为我们有前面的不等式成立）。一般情况下增大的量会比要多一些，这时候Kullback-Leiblerdivergence在新的参数 新回到“E”一步去求新的；另一方面，如果这里Kullback-LeiblerEM因此迭代的过程中得到的likelihood会逐渐近（至少是局部的）最优值。另外，在M一步中除了用最大似然之外，也可以引入先验使用umaPosteriori(MAP)的方法来做。还有一些很的问题，甚至在迭代的过程中GeneralizedEM(GEM)。面我们就来举一个用EM来做中文分词的例子。这个例子来源于论文Self-supervisedChineseWordSegmentation。我在上次MSTC搜索引现在我们有一个字符序列，并希望得到一个模型（依赖于参数）能用来将其词的序列。按照生成模型的方式来考虑，将看成是由生成的序列的话，我们自然希望找到那个生 的可能性最大的，亦即通过最大似然的方式来估计参数 然而我们不知道似然函数应该如何去优化因此我们引入latent 最大化即完成了EM的一次迭代。具体的做法通常是 的集合）按照N-gram分割（N-gram在讲K-meds的那篇文章中介绍过）为，形成一个最初的辞典而模型的参数实际上就描述了各个N-gram的概率，初始值可以直接取为频率值。有了辞典之后对于任意的，我们可以根据辞典枚举出所有可能的分割，而每个分割的后验概率就是其中的单词的概率乘积。其他的就是标准的EM算法的迭代步骤了。就是），并且这个过程是全自动的，主要有两个优点：Wikipedia（繁体中文）上抓下来的一些数据跑过一下小规模的试验，结果还可EM是怎么一回事，此为止了。Clustering4Spectral本文是“漫谈Clustering系列”中的第6篇，参列的其他文章。如果说K-means和 GMM这些聚类的方法是古代流行的算法的话那么这次要讲的SpectralClustering就可以算是现代流行的算法了，中文通常称为SpectralClustering（K-means）和K-meds 类似，SpectralClustering只需要数据之间的相似度矩阵就可以了，而不必像K-means那样要求数据必须是N维欧氏空间中的向量。于不规则的误差数据不是那么敏感，而且performance也要好一些。许多实作为baseline而存在的。K-meansK-meansK-means马，先拉出来溜溜再说。不过，在K-meds那篇文章中曾经实际跑过K-meds算法，最后的结果也就是一个accuracy，一个数字又不能画成图表之类的，看起来实在是没意思，而且K-means结果来自clusteringusinglocalityindexing这篇，这篇实际上是提的另一种聚类方法（下次如果有机会也会讲到），K-meansSpectralClusteringk Reuters-K-K-2349TDT2Reuters-21578K-meansSC(SpectralClustering)抽取了这两个数据集中若干个类别（从210类）的数据进行聚类，得出Clustering这里完胜K-means。以看到SpectralClustering算法的全貌：根据数据构造一个Graph，Graph的每一个节点对应一个数据点，将相似的点连接起来，并且边的权重用于表示数据之间的相似度。把这个Graph用邻接矩 K-meds中用的相似度矩阵作为 果中每一行所属的类别就是原来Graph中的节点亦即最初的 其实，如果你熟悉DimensionalReduction（降维）的话，大概已经看出来了，SpectralClusteringLaplacianEigenmap之后再做K-means的一个过程──听起来土多了。不过，为什么要刚好降到维呢？其实整个模型还可以从另一个角度导出来，所以，让我们不妨先ImageProcessing（我好像之前有听说我对这个领域深恶痛绝？）里有一个问题就是对图像进行Segmentation（区域分割），也就是让相似的像素组成ImageProcessing这个问题，并且有不少方法和Clustering有密切联系。比如我们在谈Vectorzation的时候就曾经用K-means来把颜色相似的像素聚类到一起，不feature（R、G、B一个像素，现在加入x、y两个新的值）即可。另一方面，还有一个经常被研究的问题就是GraphCut，简单地说就是把一个Graph的一些边切断，让他被打散成一些独立的sub-Graph，而这些被切断的边的权值的总和就被称为Cut值。如果用一张中的所有像素来组成一个Graph，并把（比如，颜色和位置上）相似的节点连接起来，边上的权值表示相似程度，那么把分割为几个区域的问题实际上等价于把Graph分sub-GraphCut切断，表示比较相似的点被保留在了同一个sub-Graph中，而彼此之间联系不cut（最小割）本身就是一个被广泛研究的问题，并且有成算法来求解。只多替代的算法提出来，如RatioCut、NormalizedCut等。述清楚。将Graph表示为邻接矩阵的形式，记为 ，其中是节 到节点的权值，如果两个节点不是相连的，权值为零设和 为Graph的两个子集（没有交集），那么两者之间的cut可先考虑最简单的情况，如果把一个GraphMinimumcut就是要最小化（其中表示的补集）。但是由于这样经常会出现孤立节点被分割出来的情况，因此又出现了RatioCut：以及NormalizedCut其中表示中的节点数目，而。两者都可算作的“大小”的一种度量，通过在分母上放置这样的项，就可以有效地防止孤立点的情况出现，达到相对平均一些的分割。事实上，JianboShi的这篇PAMIpaper：NormalizedCutsandImageSegmentation正是把NormalizedCut用在图像分割上了。搬出RatioCutNormalizedCutSpectralClustering是一个NP难问题，不方便求解，为了找到解决办法，让我们先来做做变形。令表示Graph的所有节点的集合，首先定义一个维向量 Usually,everyauthorjustcalls“his”matrixthegraph个有一个性质就是：这个是对任意向量都成立的，很好证明，只要按照定义展开就可以得到了。把我们刚才定义的那个带进去，就可以得到1到和。由于是一个常量，因最小化RatioCut就等价于最小化，当然，要记得加上附加条 以及问题转化到这个样子就好求了，因为有一个叫做Rayleighquotient的东西：征值并且极值 等于对应的特征向量时取到由于是常数因此最小化实际上也就等价于最小化，不过由于的最小的特征值为零，并且对应的特征向量正好为（我们这里仅考虑Graph是的情况），不满足的条件，因此我们取第二个小的特征值，以及对应的特征向量。们耍了一个把戏：之前的问题之所以NP难是因为向量的元素只能取两值和中的一个，是一个离散的问题，而我们求的特征向量其中的元素可以是任意实数，就是说原来的问题限制放宽了。那如何得到原来的解呢？一个最简单的办法就是看零还是小于零，将他们分别对应到离散情况的和不过我们也可以采取稍微复杂一点的办法，用的K-means来将到此为止，已经有SpectralClustering的了：求特征值，再对特征向量K-meansk（数学推导我就不再详细写了），我们就得到了和之前的SpectralClustering一模一样的步骤：求特征值并取前k个最小的，将对应的特征向量排列起来，再按行进行K-means用类似的办法，NormalizedCut也可以等价到SpectralClustering不过这次我就不再讲那么多了感的（还包括其他一些形式的GraphLaplacian以SpectralClustering和Randomwalk的关系TutorialATutorialonSpectralClustering。SpectralClusteringfunctionfunctionidx=spectral_clustering(W,k)D=diag(sum(W));L=D-opt=struct('issym',true,'isreal',true);[Vdummy]=eigs(L,D,k,'SM',idx=kmeans(V,SpectralClustering 是我们构造出来的Graph的邻接矩阵表示，通常我们在构造Graph的时候为了方便进行聚类，更加强到“局部”的连通性，亦抓住了主要，忽略了次要的东西，Performance比传统的K-means要好。SpectralClusteringK-means用LaplacianEigenmap进行降维的后的结果，如果有机会，下次谈DimensionalityReductionk（k很大），在这些低维的数据上做K-means运算量非常小。但是对于原始数据直接做K-means的话，虽然最初的数据是稀疏矩阵，但是K-means中有一个求Centroid的运算，就是求一个平均值：许多稀疏的向量的平均值求出来并不一致普通的K-means巨慢无比，而SpectralClustering等工序的算法则迅Clustering名字来源于Spectraltheory，也就是用特征分解来分析问题的UPDATE2011.11.23:LEKmeansSpectralCluster1thresholdingK类的情况，自然的类比就是降到K-1维，这也是和LDA保持一致。因为Laplacian（0），所以可以K0关于示例代码：注意除非我在这里注明了是发布某个software或者package自己试验一下具体实现和改进，否则可以直接找网上现成的的package，比如SpectralClustering可以考虑这个包。Clustering番外篇Dimensionality本文是“漫谈Clustering系列”中的第7篇，参列的其他文章。例如，一个最简单的将一个文本转化为向量的方法如下：假设待提取的文档是“Tobe,ornottobe:thatisthequestion:”，首先对其进行一些预处理，例如去掉单词的时态后缀、滤掉标点符号等，得到“tobeornottobethatbethequestion”。统计三个维度所对应的单词出现的频率：to2，be3次，the1该文档对应的向量即[2,3,1]不过事实上我们通常会做更复杂一些的处理例如如果你是在做sentimentysis那么你通常会更加关注语气很重的词比如“bad”“terrible”、“awesome”等的重要性就比普通的词要大此时你可以为每一个维度设一个权重，例如，给“bad”设置权重2，那么出现3次的话，向量在该维度对应的值就为2*3=6 是tf-idfInverseFrequency，通常使用如下公式进行计算feature一些，手工设置维度的权重固然是很人力，其实tf-idf也是在假定了原始feature是、term这样的形式（或者类似的模型）的情况下才适用featurefeatureselectiondimensionalityreductiontopic出重要的feature的维度（并抛弃不重要的维度），而后者则是更广泛意义上的值已经不一定是原始的值，也可以把featureselection看作是dimensionalityreductiontf-idffeatureselection（通常）并没有丢弃低权值的维度，并且处理过后一个函数，其输入是一个D维的向量，输出是一个M维的向量。reconstructionerror reconstructionerror另外式是简单地使用variance来衡量所包含信息量，例如，我们要把D1variance 是降维函数。推而广之，如果是降为2维，那么我希望第2维去关注第1维之外的variance，所以要求它在与第一维垂直的情况下也达到variance最大化。以此类推。然而当我们把降维函数限定维线性的时候两种途径会得到同样的结果，就是被广泛使用的PrincipalComponentsysis(PCA)。PCA的降维函数 来表示，因此，一个D维的 经过线性变换之后得到一个M维向量，就是降维的结果。 维的数据矩阵，目标是使其covariance （covariance），当然矩阵不是一个数，不能直接最化，如果我们采用矩阵的Trace 求最大化只需要求出的特征值和特征向量将M个 这也就是PCA的求解过程，得到的降维矩阵 如果熟悉LatentSemanticysis(LSA)的话，大概已经看出PCA和SingularValueposition(SVD)以及LSA之间的关系了。以下这张图（引自《TheElementsofStatisticalLearning》）可以直观地看出来PCA做了什么，对于一些原始为二维的数据，PCAvariancePCA检索等各个领域，其地位类似于聚类中的K-means，在现在关于降维等算法的baseline，PCA个就是PCA降维是线性变换，虽然线性变换计算方便，并且可以很容易地推广其中一个就是KernelPCA，用KernelTrick来将PCA推广到非线性的情况。另外，PCAGaussianlatentvariable模型，它假定数据的mean和variance是重要的特征，并依靠covariance一个典型的问题就是做聚类或分类，回想我们之前谈过的SpectralClustering，就是使用LaplacianeigenmapK-meansPCAreconstructionerror分不同类别的。如果我们有训练数据的类别，则可以用FisherLinearDiscriminantysis来处理这个问题。同PCA一样，FisherLinearDiscriminantysis也是一个线性映射模型VariancevariancevariancePatternClassification》这本书的第三章ComponentysisandDiscriminants。LinearDiscriminantysis就没有办法了。不过，如果我们假设原始的数MDS三维或者二维，用于做visualization。MDSSpectralClusteringLaplacianEigenmapMDS，LE相似度矩阵，不过这里通常要求这个相似度矩阵具有局部性质，亦即只考虑局部领域内的相似性，如果点和距离太远的话，应该为零。1.对每个点选取k个最接近的点作为邻居，与其他的点的相似性则置为零。这里需要注意的是LE要求相似度矩阵具有对称性，因此，我们通常会在 于的k个最接近的邻居且/或反之的时候，就保留的值，否则置为 那么会相对比较大，这样如果映射过后 就会被权重放大，因此最小化目标函数就保证了原来相近的点在映射过后令为将的每一行加起来所得到的对角阵，而 称作是拉斯矩阵，通过求解如下的特征值问题易知最小的那个特征值肯定是0NN类似地推广到M0M的特征向量，转置之后按行排列起来，就是降维后的结果。用代码写出来如下所示（使用了knn来构造相似度矩阵，并且没有用heatkernel）：1212345678910111281%LaplacianEigenmapALGORITHM(usingKnearest%%[Y]=%%X=dataasDxNmatrix(D=dimensionality,N=%K=numberof%dmax=maxembedding%Y=embeddingasdmaxxNfunction[Y]=[D,N]=fprintf(1,'LErunningon%dpointsin%d%STEP1:COMPUTEPAIRWISEDISTANCES&FINDX2=[sorted,index]=sort(distance);neighborhood=index(2:(1+K),:);920920212223242526272829303132333435%STEP2:ConstructsimilaritymatrixWW=zeros(N,N);forii=1:NW(ii,neighborhood(:,ii))=1;W(neighborhood(:,ii),ii)=%STEP3:COMPUTEEMBEDDINGFROMEIGENVECTSOFL=D-W;%CALCULATIONOFoptions.disp=0;options.isreal=1;options.issym=1;[Y,eigenvals]=eigs(L,d+1,0,options);Y=Y(:,2:d+1)'*sqrt(N);%bottomevectis[1,1,1,1...]witheval3363738394044748495051522事实上，LaplacianEigenmap假设数据分布在一个嵌套在高中的低维流LaplacianMatrix则是流形的LaplaceBeltramioperator的LaplacianEigenmap这里做地展开了，下面看一个比较直观的例子SwissRoll。SwissRoll是一个像面包圈一样的结构，可以看作一个嵌套在三中的二SwissRollSwissRollLaplacianEigenmapPCASwissRoll看到LE成功地把这个流形展开把在流形上属于不同位置的点映射到了不同的地方，而PCA的结果则很糟糕，几种颜色都混杂在了一起。LocallyLinearEmbeddingLE 应当等于零。这之后再把 际上，从理论上也可以得出二者之间的联系（LE对应于 而LLE对应于），如果感的话，可以参考LaplacianEigenmapsforDimensionalityReductionandDataRepresentation这篇里的对比。下面是两种方法在SwissRoll数据上的结果，也是非常相似的：影矩阵，这个结果可以保存下来，在之后对任意向量进行投影，而LE和LLEPCALELLE别叫做LocalityPreservingProjection和NeighborhoodPreservingEmbedding。LPPPCA阵来表示，于是LE的目标函数变为得到的按照特征值从小到大排序的特征向量就组成映射矩阵LELE里的矩阵在乘以之后通常就不再稀疏了，计算量会变得比较大，这个问题可以使用SpectralRegression的方法来解决，参见SpectralRegression:AUnifiedApproachforSparseSubspaceLearningpaperKernelTrickLPPLELLENPE另外，虽然LE是unsupervised的，但是如果训练数据确实有可用，也是的决策过程其实人是可以“看到”或者说“理解”的，但是SVM就不那么法发展出“智能”来啊^_^bbClustering5Hierarchical本文是“漫谈Clustering系列”中的第8篇，见本系列的其他文不过觉得HierarchicalClustering这个话题我能说的东西应该不多，所以还是先写了吧（我准备这次一个公式都不贴）。Hierarchicalnaive不过言归正传，这里要说的HierarchicalClustering听上去很，不过其实HierarchicalClustering的想法很简单，主要分为两大类：agglomerative（自底向上）divisive（自顶向下）。首先说前者，看起来很简单吧？其实确实也是比较简单的，不过还是有两个问题需要problemdependentSingleLinkage：又叫做nearest-neighbor的两个点的距离作为这两个集合的距离，容易造成一种叫做Chaining的效离比较近就被合并了，并且这样合并之后Chaining效应会进一步扩大，最后会得到比较松散的cluster。CompleteLinkage：这个则完全是SingleLinkage的，取两个集SingleLinkageCompleteLinkage的方法。整个agglomerativehierarchicalclustering的算法就是这个离最近的两个点，需要一个双重循环，而且GroupAverage计算距离的时候也123456789123456789def#makeacopy,donottouchtheoriginallistnodes