今年 6 月份，Alex Rodriguez 和 Alessandro Laio 在 Science 上發表了一篇名為《Clustering by fast search and find of density peaks》的文章，為聚類演算法的設計提供了一種新的思路。雖然文章出來後遭到了眾多讀者的質疑，但整體而言，新聚類演算法的基本思想很新穎，且簡單明快，值得學習。這個新聚類演算法的核心思想在於對聚類中心的刻畫上，本文將對該演算法的原理進行詳細介紹，並對其中的若干細節展開討論。

最後，附上作者在補充材料裡提供的 Matlab 示例程式 （加了適當的程式碼註釋）。

1. clear all  
2. close all  
3. disp('The only input needed is a distance matrix file')  
4. disp('The format of this file should be: ')  
5. disp('Column 1: id of element i')  
6. disp('Column 2: id of element j')  
7. disp('Column 3: dist(i,j)')  
8. %% 從檔案中讀取資料  
9. mdist=input('name of the distance matrix file (with single quotes)?\n');  
10. disp('Reading input distance matrix')  
11. xx=load(mdist);  
12. ND=max(xx(:,2));  
13. NL=max(xx(:,1));  
14. if (NL>ND)  
15.   ND=NL;  %% 確保 DN 取為第一二列最大值中的較大者，並將其作為資料點總數  
16. end  
17. N=size(xx,1); %% xx 第一個維度的長度，相當於檔案的行數（即距離的總個數）  
18. %% 初始化為零  
19. for i=1:ND  
20.   for j=1:ND  
21.     dist(i,j)=0;  
22.   end  
23. end  
24. %% 利用 xx 為 dist 陣列賦值，注意輸入只存了 0.5*DN(DN-1) 個值，這裡將其補成了滿矩陣  
25. %% 這裡不考慮對角線元素  
26. for i=1:N  
27.   ii=xx(i,1);  
28.   jj=xx(i,2);  
29.   dist(ii,jj)=xx(i,3);  
30.   dist(jj,ii)=xx(i,3);  
31. end  
32. %% 確定 dc  
33. percent=2.0;  
34. fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent);  
35. position=round(N*percent/100); %% round 是一個四捨五入函式  
36. sda=sort(xx(:,3)); %% 對所有距離值作升序排列  
37. dc=sda(position);  
38. %% 計算區域性密度 rho (利用 Gaussian 核)  
39. fprintf('Computing Rho with gaussian kernel of radius: %12.6f\n', dc);  
40. %% 將每個資料點的 rho 值初始化為零  
41. for i=1:ND  
42.   rho(i)=0.;  
43. end  
44. % Gaussian kernel  
45. for i=1:ND-1  
46.   for j=i+1:ND  
47.      rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));  
48.      rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));  
49.   end  
50. end  
51. % "Cut off" kernel  
52. %for i=1:ND-1  
53. %  for j=i+1:ND  
54. %    if (dist(i,j)<dc)  
55. %       rho(i)=rho(i)+1.;  
56. %       rho(j)=rho(j)+1.;  
57. %    end  
58. %  end  
59. %end  
60. %% 先求矩陣列最大值，再求最大值，最後得到所有距離值中的最大值  
61. maxd=max(max(dist));   
62. %% 將 rho 按降序排列，ordrho 保持序  
63. [rho_sorted,ordrho]=sort(rho,'descend');  
64. %% 處理 rho 值最大的資料點  
65. delta(ordrho(1))=-1.;  
66. nneigh(ordrho(1))=0;  
67. %% 生成 delta 和 nneigh 陣列  
68. for ii=2:ND  
69.    delta(ordrho(ii))=maxd;  
70.    for jj=1:ii-1  
71.      if(dist(ordrho(ii),ordrho(jj))<delta(ordrho(ii)))  
72.         delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj));  
73.         nneigh(ordrho(ii))=ordrho(jj);   
74.         %% 記錄 rho 值更大的資料點中與 ordrho(ii) 距離最近的點的編號 ordrho(jj)  
75.      end  
76.    end  
77. end  
78. %% 生成 rho 值最大資料點的 delta 值  
79. delta(ordrho(1))=max(delta(:));  
80. %% 決策圖  
81. disp('Generated file:DECISION GRAPH')   
82. disp('column 1:Density')  
83. disp('column 2:Delta')  
84. fid = fopen('DECISION_GRAPH', 'w');  
85. for i=1:ND  
86.    fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));  
87. end  
88. %% 選擇一個圍住類中心的矩形  
89. disp('Select a rectangle enclosing cluster centers')  
90. %% 每臺計算機，控制代碼的根物件只有一個，就是螢幕，它的控制代碼總是 0  
91. %% >> scrsz = get(0,'ScreenSize')  
92. %% scrsz =  
93. %%            1           1        1280         800  
94. %% 1280 和 800 就是你設定的計算機的解析度，scrsz(4) 就是 800，scrsz(3) 就是 1280  
95. scrsz = get(0,'ScreenSize');  
96. %% 人為指定一個位置，感覺就沒有那麼 auto 了 :-)  
97. figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);  
98. %% ind 和 gamma 在後面並沒有用到  
99. for i=1:ND  
100.   ind(i)=i;   
101.   gamma(i)=rho(i)*delta(i);  
102. end  
103. %% 利用 rho 和 delta 畫出一個所謂的“決策圖”  
104. subplot(2,1,1)  
105. tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');  
106. title ('Decision Graph','FontSize',15.0)  
107. xlabel ('\rho')  
108. ylabel ('\delta')  
109. subplot(2,1,1)  
110. rect = getrect(1);   
111. %% getrect 從圖中用滑鼠擷取一個矩形區域， rect 中存放的是  
112. %% 矩形左下角的座標 (x,y) 以及所截矩形的寬度和高度  
113. rhomin=rect(1);  
114. deltamin=rect(2); %% 作者承認這是個 error，已由 4 改為 2 了!  
115. %% 初始化 cluster 個數  
116. NCLUST=0;  
117. %% cl 為歸屬標誌陣列，cl(i)=j 表示第 i 號資料點歸屬於第 j 個 cluster  
118. %% 先統一將 cl 初始化為 -1  
119. for i=1:ND  
120.   cl(i)=-1;  
121. end  
122. %% 在矩形區域內統計資料點（即聚類中心）的個數  
123. for i=1:ND  
124.   if ( (rho(i)>rhomin) && (delta(i)>deltamin))  
125.      NCLUST=NCLUST+1;  
126.      cl(i)=NCLUST; %% 第 i 號資料點屬於第 NCLUST 個 cluster  
127.      icl(NCLUST)=i;%% 逆對映,第 NCLUST 個 cluster 的中心為第 i 號資料點  
128.   end  
129. end  
130. fprintf('NUMBER OF CLUSTERS: %i \n', NCLUST);  
131. disp('Performing assignation')  
132. %% 將其他資料點歸類 (assignation)  
133. for i=1:ND  
134.   if (cl(ordrho(i))==-1)  
135.     cl(ordrho(i))=cl(nneigh(ordrho(i)));  
136.   end  
137. end  
138. %% 由於是按照 rho 值從大到小的順序遍歷,迴圈結束後, cl 應該都變成正的值了.   
139. %% 處理光暈點，halo這段程式碼應該移到 if (NCLUST>1) 內去比較好吧  
140. for i=1:ND  
141.   halo(i)=cl(i);  
142. end  
143. if (NCLUST>1)  
144.   % 初始化陣列 bord_rho 為 0,每個 cluster 定義一個 bord_rho 值  
145.   for i=1:NCLUST  
146.     bord_rho(i)=0.;  
147.   end  
148.   % 獲取每一個 cluster 中平均密度的一個界 bord_rho  
149.   for i=1:ND-1  
150.     for j=i+1:ND  
151.       %% 距離足夠小但不屬於同一個 cluster 的 i 和 j  
152.       if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))  
153.         rho_aver=(rho(i)+rho(j))/2.; %% 取 i,j 兩點的平均區域性密度  
154.         if (rho_aver>bord_rho(cl(i)))   
155.           bord_rho(cl(i))=rho_aver;  
156.         end  
157.         if (rho_aver>bord_rho(cl(j)))   
158.           bord_rho(cl(j))=rho_aver;  
159.         end  
160.       end  
161.     end  
162.   end  
163.   %% halo 值為 0 表示為 outlier  
164.   for i=1:ND  
165.     if (rho(i)<bord_rho(cl(i)))  
166.       halo(i)=0;  
167.     end  
168.   end  
169. end  
170. %% 逐一處理每個 cluster  
171. for i=1:NCLUST  
172.   nc=0; %% 用於累計當前 cluster 中資料點的個數  
173.   nh=0; %% 用於累計當前 cluster 中核心資料點的個數  
174.   for j=1:ND  
175.     if (cl(j)==i)   
176.       nc=nc+1;  
177.     end  
178.     if (halo(j)==i)   
179.       nh=nh+1;  
180.     end  
181.   end  
182.   fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(i),nc,nh,nc-nh);  
183. end  
184. cmap=colormap;  
185. for i=1:NCLUST  
186.    ic=int8((i*64.)/(NCLUST*1.));  
187.    subplot(2,1,1)  
188.    hold on  
189.    plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
190. end  
191. subplot(2,1,2)  
192. disp('Performing 2D nonclassical multidimensional scaling')  
193. Y1 = mdscale(dist, 2, 'criterion','metricstress');  
194. plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','MarkerEdgeColor','k');  
195. title ('2D Nonclassical multidimensional scaling','FontSize',15.0)  
196. xlabel ('X')  
197. ylabel ('Y')  
198. for i=1:ND  
199.  A(i,1)=0.;  
200.  A(i,2)=0.;  
201. end  
202. for i=1:NCLUST  
203.   nn=0;  
204.   ic=int8((i*64.)/(NCLUST*1.));  
205.   for j=1:ND  
206.     if (halo(j)==i)  
207.       nn=nn+1;  
208.       A(nn,1)=Y1(j,1);  
209.       A(nn,2)=Y1(j,2);  
210.     end  
211.   end  
212.   hold on  
213.   plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
214. end  
215. %for i=1:ND  
216. %   if (halo(i)>0)  
217. %      ic=int8((halo(i)*64.)/(NCLUST*1.));  
218. %      hold on  
219. %      plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));  
220. %   end  
221. %end  
222. faa = fopen('CLUSTER_ASSIGNATION', 'w');  
223. disp('Generated file:CLUSTER_ASSIGNATION')  
224. disp('column 1:element id')  <