密度峰值聚類算法(DPC)
凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/
具體實例見：密度峰值聚類算法MATLAB程序 - 凱魯嘎吉 - 博客園
1. 簡介
基於密度峰值的聚類算法全稱為基於快速搜索和發現密度峰值的聚類算法(clustering by fast search and find of density peaks, DPC)。它是2014年在Science上提出的聚類算法，該算法能夠自動地發現簇中心，實現任意形狀數據的高效聚類。
該算法基於兩個基本假設：1）簇中心（密度峰值點）的局部密度大於圍繞它的鄰居的局部密度；2）不同簇中心之間的距離相對較遠。為了找到同時滿足這兩個條件的簇中心，該算法引入了局部密度的定義。
假設數據點
的局部密度為
，數據點
到局部密度比它大且距離最近的數據點
的距離為
，則有如下定義：
式中，
為
和
之間的距離；
為截斷距離；
為邏輯判斷函數，
，否則
。
這里對於局部密度最大的數據點
，它的
。
根據以上定義，通過構造
相對於
的決策圖，進行數據點分配和噪聲點剔除，可以快速得到最終的聚類結果。算法1給出了基於快速搜索和發現密度峰值的聚類算法的具體步驟。首先，基於快速搜索和發現密度峰值的聚類算法對任意兩個數據點計算它們之間的距離，並依據截斷距離計算出任意數據點
的
和
；然后，算法根據
和
，畫出對應的聚類決策圖；接着，算法利用得到的決策圖，將
和
都相對較高的數據點標記為簇的中心，將
相對較低但是
相對較高的點標記為噪聲點；最后，算法將剩余的數據點進行分配，分配的規則為將每個剩余的數據點分配到它的最近鄰且密度比其大的數據點所在的簇。
算法1 基於快速搜索和發現密度峰值的聚類算法
下面舉一個簡單的例子。在圖1中，數據點的分布情況如圖1（左）所示，可以看出數據集包含兩個簇，分別用藍色和紅色標出，噪聲用黑色標出。利用基於快速搜索和發現密度峰值算法，可以得到圖1（右）的決策圖。在決策過程中可以發現點1和點10的
和
都相對較高，因此它們被標記為中心點。點26~點28的
相對較低但是
相對較高，因此它們被標記為噪聲。其他的點將被分配到它的最近鄰且密度比其大的數據點所在的簇中去。
圖1 基於快速搜索和發現密度峰值的聚類算法例子
基於快速搜索和發現密度峰值的聚類算法，想法非常直觀，能夠快速發現密度峰值點，並能夠高效進行樣本分配和發現噪聲點。同時，因為該方法非常適用於大規模數據的聚類分析，因此具有很好的研究價值和應用前景。
2. MATLAB程序
數據請參考文獻[3]
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\n','s');
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');
%% 人為指定一個位置
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')
fig=subplot(2,1,1);
rect = getrect(fig);
%% 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
%% 將其他數據點歸類 (assignation)
for i=1:ND
if (cl(ordrho(i))==-1)
cl(ordrho(i))=cl(nneigh(ordrho(i)));
end
end
%halo
%% 由於是按照 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
3. 結果
>> cluster_dp The only input needed is a distance matrix file The format of this file should be: Column 1: id of element i Column 2: id of element j Column 3: dist(i,j) name of the distance matrix file example_distances.dat Reading input distance matrix average percentage of neighbours (hard coded): 2.000000 Computing Rho with gaussian kernel of radius: 0.033000 Generated file:DECISION GRAPH column 1:Density column 2:Delta Select a rectangle enclosing cluster centers NUMBER OF CLUSTERS: 5 Performing assignation CLUSTER: 1 CENTER: 149 ELEMENTS: 378 CORE: 260 HALO: 118 CLUSTER: 2 CENTER: 451 ELEMENTS: 326 CORE: 250 HALO: 76 CLUSTER: 3 CENTER: 1310 ELEMENTS: 884 CORE: 785 HALO: 99 CLUSTER: 4 CENTER: 1349 ELEMENTS: 297 CORE: 208 HALO: 89 CLUSTER: 5 CENTER: 1579 ELEMENTS: 115 CORE: 115 HALO: 0 Performing 2D nonclassical multidimensional scaling Generated file:CLUSTER_ASSIGNATION column 1:element id column 2:cluster assignation without halo control column 3:cluster assignation with halo control
注：出錯的話，將Y1 = mdscale(dist, 2, 'criterion','metricstress');換一個准則函數，比如改為Y1 = mdscale(dist, 2, 'criterion','sstress');
4. 參考文獻
[1] Rodriguez A, Laio A. Clustering by fast search and find of density peaks [J]. Science, 2014, 344(6191): 1492-1496.
[2] 張憲超. 數據聚類. 北京：科學出版社, 2017.06.
[3] MATLAB程序：Clustering by fast search-and-find of density peaks
[4] 密度峰值聚類算法MATLAB程序