机器学习：Clustering by fast search and find of density peaks

本文主要是介绍机器学习：Clustering by fast search and find of density peaks，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

聚类算法，作为机器学习里常用的一种无监督方法，一直以来都受到很大的关注，聚类算法，简单来说就是希望把同一类的样本或者样本聚到一起，比如说常见的图像分类，我们希望猫的图片能聚到一起，狗的图片能聚到一起，不希望猫和狗的图片混在一起。

聚类算法发展到现在，已经有了很多不同的算法，常见的比如 K-means 就是比较经典的一种聚类算法。

今天介绍发表在 science 2014 上的一篇文章，Clustering by fast search and find of density peaks，这篇文章对聚类中心的的性质做了两点归纳，或者说从两个维度来进行刻画：

密度：聚类中心的密度要大，也就是说，围绕着这个中心的点应该要尽可能的多
距离：与其他聚类中心的距离要尽可能的远，这个也比较直观

给定一个数据集， $S=\{ \mathbf{x}_i \}_{i=1}^{N}$ , $d_{ij} = dis(\mathbf{x}_i, \mathbf{x}_j)$ 表示数据点 $\mathbf{x}_i$ 和 $\mathbf{x}_j$ 之间的距离，对于 $S$ 中的任何一个点 $\mathbf{x}_i$ ，都可以求出该点的局部密度 $\rho_i$ 和相对距离 $\sigma_i$ ，首先来看局部密度的计算：

局部密度 $\rho_i$

有两种计算方式：第一种方式：

$\rho_{i} = \sum_{j \in I_S \setminus \{i\}} \mathcal{X} (d_{ij} - d_c)$

其中函数：

$\mathcal{X} (x)= \{ \begin{matrix} 1 \quad x < 0 \\ 0 \quad x \geq 0 \end{matrix}$

第二种方式：

$\rho_{i} = \sum_{j \in I_S \setminus \{i\}} e^{-({\frac{d_{ij}}{d_c})^{2}}}$

其中 $d_c$ 表示截断距离，由用户指定。第一种方式可以看成是一种 hard mapping，第二种方式是一种 soft mapping。这两种方式都是去计算围绕着某个点附近的数据点的数量，类似以某个点为中心，以 $d_c$ 为半径画了一个圆圈，计算落在圈里的点数，第二种方式可能更加soft一点，根据距离来表示数据点与中心点的权重，离的越近，权重越高，离得越远，权重就越低。

除了局部密度之外，另外一个指标就是相对距离，

-相对距离 $\sigma_{i}$

假设 ${ q_i\}_{i=1}^{N}$ 表示 $\{ \rho_i\}_{i=1}^{N}$ 的一个降序排序，即：

$\rho_{q_1} \geq \rho_{q_2} \geq \rho_{q_3} \geq ... \geq \rho_{q_N}$

则可定义，

$\sigma_{q_i}= \{ \begin{matrix} \min_{q_j, j < i} \{ d_{q_i q_j} \} \quad i \geq 2 \\ \max_{j \geq 2} \{ \sigma_{q_j} \} \quad i = 1 \end{matrix}$

相对距离，简单来说就是这样，首先根据上面的局部密度，可以对所有的数据点进行排序，比如说根据密度从大到小进行排序，然后对于密度最大的那个点，其相对距离就取除改点之外的局部密度最大值作为该点的相对距离，而其他的点，就从所有局部密度大于该点的数据点中，离该点最近的距离作为相对距离。

结合这两个指标，局部密度大并且相对距离大的点，就可以看成是一个聚类中心。

最后，附上作者的代码：

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 取为第一二列最大值中的较大者，并将其作为数据点总数
endN=size(xx,1); %% xx 第一个维度的长度，相当于文件的行数（即距离的总个数）%% 初始化为零
for i=1:NDfor j=1:NDdist(i,j)=0;end
end%% 利用 xx 为 dist 数组赋值，注意输入只存了 0.5*DN(DN-1) 个值，这里将其补成了满矩阵
%% 这里不考虑对角线元素
for i=1:Nii=xx(i,1);jj=xx(i,2);dist(ii,jj)=xx(i,3);dist(jj,ii)=xx(i,3);
end%% 确定 dcpercent=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:NDrho(i)=0.;
end% Gaussian kernel
for i=1:ND-1for j=i+1:NDrho(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:NDdelta(ordrho(ii))=maxd;for jj=1:ii-1if(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)endend
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:NDfprintf(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:NDind(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:NDcl(i)=-1;
end%% 在矩形区域内统计数据点（即聚类中心）的个数
for i=1:NDif ( (rho(i)>rhomin) && (delta(i)>deltamin))NCLUST=NCLUST+1;cl(i)=NCLUST; %% 第 i 号数据点属于第 NCLUST 个 clustericl(NCLUST)=i;%% 逆映射,第 NCLUST 个 cluster 的中心为第 i 号数据点end
endfprintf('NUMBER OF CLUSTERS: %i \n', NCLUST);disp('Performing assignation')%% 将其他数据点归类 (assignation)
for i=1:NDif (cl(ordrho(i))==-1)cl(ordrho(i))=cl(nneigh(ordrho(i)));end
end
%% 由于是按照 rho 值从大到小的顺序遍历,循环结束后, cl 应该都变成正的值了. %% 处理光晕点，halo这段代码应该移到 if (NCLUST>1) 内去比较好吧
for i=1:NDhalo(i)=cl(i);
endif (NCLUST>1)% 初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值for i=1:NCLUSTbord_rho(i)=0.;end% 获取每一个 cluster 中平均密度的一个界 bord_rhofor i=1:ND-1for j=i+1:ND%% 距离足够小但不属于同一个 cluster 的 i 和 jif ((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;endif (rho_aver>bord_rho(cl(j))) bord_rho(cl(j))=rho_aver;endendendend%% halo 值为 0 表示为 outlierfor i=1:NDif (rho(i)<bord_rho(cl(i)))halo(i)=0;endendend%% 逐一处理每个 cluster
for i=1:NCLUSTnc=0; %% 用于累计当前 cluster 中数据点的个数nh=0; %% 用于累计当前 cluster 中核心数据点的个数for j=1:NDif (cl(j)==i) nc=nc+1;endif (halo(j)==i) nh=nh+1;endendfprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(i),nc,nh,nc-nh);endcmap=colormap;
for i=1:NCLUSTic=int8((i*64.)/(NCLUST*1.));subplot(2,1,1)hold onplot(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:NDA(i,1)=0.;A(i,2)=0.;
end
for i=1:NCLUSTnn=0;ic=int8((i*64.)/(NCLUST*1.));for j=1:NDif (halo(j)==i)nn=nn+1;A(nn,1)=Y1(j,1);A(nn,2)=Y1(j,2);endendhold onplot(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:NDfprintf(faa, '%i %i %i\n',i,cl(i),halo(i));
end