学习记录发表在 Science 上的一种新聚类算法及本系列聚类算法的其他链接

最新推荐文章于 2024-04-27 16:11:04 发布

An_Wu

最新推荐文章于 2024-04-27 16:11:04 发布

阅读量989

点赞数

分类专栏： MATLAB，聚类算法学习记录，社团识别

学习记录，社团识别同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

MATLAB，聚类算法

1 篇文章 0 订阅

订阅专栏

作者: peghoty 
出处: http://blog.csdn.net/itplus/article/details/38926837
欢迎转载/分享, 但请务必声明文章出处.
    
  
    2014年 6 月份，Alex Rodriguez 和 Alessandro Laio 在 Science 上发表了一篇名为《Clustering by fast search and find of density peaks》的文章，为聚类算法的设计提供了一种新的思路。虽然文章出来后遭到了众多读者的质疑，但整体而言，新聚类算法的基本思想很新颖，且简单明快，值得学习。这个新聚类算法的核心思想在于对聚类中心的刻画上，本文将对该算法的原理进行详细介绍，并对其中的若干细节展开讨论。

最后，附上作者在补充材料里提供的 Matlab 示例程序 （加了适当的代码注释）。

[plain]view plaincopy 
   
 clear all  
 close all  
 disp('The only input needed is a distance matrix file')  
 disp('The format of this file should be: ')  
 disp('Column 1: id of element i')  
 disp('Column 2: id of element j')  
 disp('Column 3: dist(i,j)')  
   
 %% 从文件中读取数据  
 mdist=input('name of the distance matrix file (with single quotes)?\n');  
 disp('Reading input distance matrix')  
 xx=load(mdist);  
 ND=max(xx(:,2));  
 NL=max(xx(:,1));  
 if (NL>ND)  
   ND=NL;  %% 确保 DN 取为第一二列最大值中的较大者，并将其作为数据点总数  
 end  
   
 N=size(xx,1); %% xx 第一个维度的长度，相当于文件的行数（即距离的总个数）  
   
 %% 初始化为零  
 for i=1:ND  
   for j=1:ND  
     dist(i,j)=0;  
   end  
 end  
   
 %% 利用 xx 为 dist 数组赋值，注意输入只存了 0.5*DN(DN-1) 个值，这里将其补成了满矩阵  
 %% 这里不考虑对角线元素  
 for i=1:N  
   ii=xx(i,1);  
   jj=xx(i,2);  
   dist(ii,jj)=xx(i,3);  
   dist(jj,ii)=xx(i,3);  
 end  
   
 %% 确定 dc  
   
 percent=2.0;  
 fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent);  
   
 position=round(N*percent/100); %% round 是一个四舍五入函数  
 sda=sort(xx(:,3)); %% 对所有距离值作升序排列  
 dc=sda(position);  
   
 %% 计算局部密度 rho (利用 Gaussian 核)  
   
 fprintf('Computing Rho with gaussian kernel of radius: %12.6f\n', dc);  
   
 %% 将每个数据点的 rho 值初始化为零  
 for i=1:ND  
   rho(i)=0.;  
 end  
   
 % Gaussian kernel  
 for i=1:ND-1  
   for j=i+1:ND  
      rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));  
      rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));  
   end  
 end  
   
 % "Cut off" kernel  
 %for i=1:ND-1  
 %  for j=i+1:ND  
 %    if (dist(i,j)<dc)  
 %       rho(i)=rho(i)+1.;  
 %       rho(j)=rho(j)+1.;  
 %    end  
 %  end  
 %end  
   
 %% 先求矩阵列最大值，再求最大值，最后得到所有距离值中的最大值  
 maxd=max(max(dist));   
   
 %% 将 rho 按降序排列，ordrho 保持序  
 [rho_sorted,ordrho]=sort(rho,'descend');  
    
 %% 处理 rho 值最大的数据点  
 delta(ordrho(1))=-1.;  
 nneigh(ordrho(1))=0;  
   
 %% 生成 delta 和 nneigh 数组  
 for ii=2:ND  
    delta(ordrho(ii))=maxd;  
    for jj=1:ii-1  
      if(dist(ordrho(ii),ordrho(jj))<delta(ordrho(ii)))  
         delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj));  
         nneigh(ordrho(ii))=ordrho(jj);   
         %% 记录 rho 值更大的数据点中与 ordrho(ii) 距离最近的点的编号 ordrho(jj)  
      end  
    end  
 end  
   
 %% 生成 rho 值最大数据点的 delta 值  
 delta(ordrho(1))=max(delta(:));  
   
 %% 决策图  
   
 disp('Generated file:DECISION GRAPH')   
 disp('column 1:Density')  
 disp('column 2:Delta')  
   
 fid = fopen('DECISION_GRAPH', 'w');  
 for i=1:ND  
    fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));  
 end  
   
 %% 选择一个围住类中心的矩形  
 disp('Select a rectangle enclosing cluster centers')  
   
 %% 每台计算机，句柄的根对象只有一个，就是屏幕，它的句柄总是 0  
 %% >> scrsz = get(0,'ScreenSize')  
 %% scrsz =  
 %%            1           1        1280         800  
 %% 1280 和 800 就是你设置的计算机的分辨率，scrsz(4) 就是 800，scrsz(3) 就是 1280  
 scrsz = get(0,'ScreenSize');  
   
 %% 人为指定一个位置，感觉就没有那么 auto 了 :-)  
 figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);  
   
 %% ind 和 gamma 在后面并没有用到  
 for i=1:ND  
   ind(i)=i;   
   gamma(i)=rho(i)*delta(i);  
 end  
   
 %% 利用 rho 和 delta 画出一个所谓的“决策图”  
   
 subplot(2,1,1)  
 tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');  
 title ('Decision Graph','FontSize',15.0)  
 xlabel ('\rho')  
 ylabel ('\delta')  
   
 subplot(2,1,1)  
 rect = getrect(1);   
 %% getrect 从图中用鼠标截取一个矩形区域， rect 中存放的是  
 %% 矩形左下角的坐标 (x,y) 以及所截矩形的宽度和高度  
 rhomin=rect(1);  
 deltamin=rect(2); %% 作者承认这是个 error，已由 4 改为 2 了!  
   
 %% 初始化 cluster 个数  
 NCLUST=0;  
   
 %% cl 为归属标志数组，cl(i)=j 表示第 i 号数据点归属于第 j 个 cluster  
 %% 先统一将 cl 初始化为 -1  
 for i=1:ND  
   cl(i)=-1;  
 end  
   
 %% 在矩形区域内统计数据点（即聚类中心）的个数  
 for i=1:ND  
   if ( (rho(i)>rhomin) && (delta(i)>deltamin))  
      NCLUST=NCLUST+1;  
      cl(i)=NCLUST; %% 第 i 号数据点属于第 NCLUST 个 cluster  
      icl(NCLUST)=i;%% 逆映射,第 NCLUST 个 cluster 的中心为第 i 号数据点  
   end  
 end  
   
 fprintf('NUMBER OF CLUSTERS: %i \n', NCLUST);  
   
 disp('Performing assignation')  
   
 %% 将其他数据点归类 (assignation)  
 for i=1:ND  
   if (cl(ordrho(i))==-1)  
     cl(ordrho(i))=cl(nneigh(ordrho(i)));  
   end  
 end  
 %% 由于是按照 rho 值从大到小的顺序遍历,循环结束后, cl 应该都变成正的值了.   
   
 %% 处理光晕点，halo这段代码应该移到 if (NCLUST>1) 内去比较好吧  
 for i=1:ND  
   halo(i)=cl(i);  
 end  
   
 if (NCLUST>1)  
   
   % 初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值  
   for i=1:NCLUST  
     bord_rho(i)=0.;  
   end  
   
   % 获取每一个 cluster 中平均密度的一个界 bord_rho  
   for i=1:ND-1  
     for j=i+1:ND  
       %% 距离足够小但不属于同一个 cluster 的 i 和 j  
       if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))  
         rho_aver=(rho(i)+rho(j))/2.; %% 取 i,j 两点的平均局部密度  
         if (rho_aver>bord_rho(cl(i)))   
           bord_rho(cl(i))=rho_aver;  
         end  
         if (rho_aver>bord_rho(cl(j)))   
           bord_rho(cl(j))=rho_aver;  
         end  
       end  
     end  
   end  
   
   %% halo 值为 0 表示为 outlier  
   for i=1:ND  
     if (rho(i)<bord_rho(cl(i)))  
       halo(i)=0;  
     end  
   end  
   
 end  
   
 %% 逐一处理每个 cluster  
 for i=1:NCLUST  
   nc=0; %% 用于累计当前 cluster 中数据点的个数  
   nh=0; %% 用于累计当前 cluster 中核心数据点的个数  
   for j=1:ND  
     if (cl(j)==i)   
       nc=nc+1;  
     end  
     if (halo(j)==i)   
       nh=nh+1;  
     end  
   end  
   
   fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(i),nc,nh,nc-nh);  
   
 end  
   
 cmap=colormap;  
 for i=1:NCLUST  
    ic=int8((i*64.)/(NCLUST*1.));  
    subplot(2,1,1)  
    hold on  
    plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
 end  
 subplot(2,1,2)  
 disp('Performing 2D nonclassical multidimensional scaling')  
 Y1 = mdscale(dist, 2, 'criterion','metricstress');  
 plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','MarkerEdgeColor','k');  
 title ('2D Nonclassical multidimensional scaling','FontSize',15.0)  
 xlabel ('X')  
 ylabel ('Y')  
 for i=1:ND  
  A(i,1)=0.;  
  A(i,2)=0.;  
 end  
 for i=1:NCLUST  
   nn=0;  
   ic=int8((i*64.)/(NCLUST*1.));  
   for j=1:ND  
     if (halo(j)==i)  
       nn=nn+1;  
       A(nn,1)=Y1(j,1);  
       A(nn,2)=Y1(j,2);  
     end  
   end  
   hold on  
   plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
 end  
   
 %for i=1:ND  
 %   if (halo(i)>0)  
 %      ic=int8((halo(i)*64.)/(NCLUST*1.));  
 %      hold on  
 %      plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
 %   end  
 %end  
 faa = fopen('CLUSTER_ASSIGNATION', 'w');  
 disp('Generated file:CLUSTER_ASSIGNATION')  
 disp('column 1:element id')  
 disp('column 2:cluster assignation without halo control')  
 disp('column 3:cluster assignation with halo control')  
 for i=1:ND  
    fprintf(faa, '%i %i %i\n',i,cl(i),halo(i));  
 end  