K-Means 和K-Medoids算法及其MATLAB实现

[cpp]  view plain copy

1.   

K-Means和K-Medoids算法是学习领域比较普通的聚类算法(无监督学习)，本文介绍原理及Matlab实现代码。

1.问题：

     给定数据点集P，d-by-N，将这些数据点集聚类到K类中去<K是给定的类的数目，可以不给定K，但本文我们处理不那么复杂的问题>

同时要求下式值最小：Sk是聚类形成的数据集合，mk是每个类集合的“中心”——K-Means与K-Medoids唯一不同的地方

                                                       

演示图：

2.K-Means算法：

       1. 将数据分为K个非空子集

       2. 计算每个类中心 

       3. 将每个数据点 xj 到最近的 mk 

       4. 返回2，当聚类结果(如：计算得到的中心m1—mk，本文采用的方法)不再变化的时候stop

3.K-Medoids算法:

       1. 随机选择K个点作为初始medoid

       2.将每个数据点分配到最近的medoid

      3. 更新每个类的medoid ——此步导致比K-Means算法的计算量加大

                                                

       4. 返回2，当各类medoid不再变化的时候stop

4.特点：

       -聚类结果与初始点有关（因为是做steepest descent from a random initial starting oint）

       -是局部最优解

       -在实际做的时候，随机选择多组初始点，最后选择拥有最低TSD（Totoal Squared Distance）的那组

=======================================

用matlab实现上面两个算法：

说明：   输入：P=[x1 x2 ... xN]  N个数据点集，K 类的数目，method可选‘K-Means' 或者' K-Medoids'            输出： best_Label 最终得到的最优聚类标签<对应P中点的顺序>； best_Label 最佳聚类中心。

[cpp]  view plain copy

1. function [best_Label best_Center best_ind label] = KM(P,K,method)  
2. %%%%-----------------------------------------------------------------------  
3. %   Version 1.0   
4. %   Author: feitengli@foxmail.com   from DUT  
5. %   CreateTime: 2012-11-29  
6. %%%------------------------------------------------------------------------  
7. %KM   K-Means Clustering or K-Medoids Clustering  
8. %    P is an d-by-N data matrix   
9. %    K is the clustering number  
10. %    method = KMeans    :K-Means Clustering  
11. %           = KMedoids  :K-Medoids Clustering  
12. %References：  
13. %        1.The Elements of Statistical Learning 2nd Chapter14.3.6&&14.3.10  
14. %%%%-----------------------------------------------------------------------  
15.   
16. [d N] = size(P);   
17. %% 本算法要求数据矩阵P的每列代表一个数据点，如果不是 需要转置矩阵  
18. if d > N  
19.     ButtonName = questdlg('数据维数小于点的个数,是否转置矩阵',...  
20.                           'MATLAB quest','Yes','No','Yes');  
21.     if  strcmp(ButtonName, 'Yes')  
22.         P = P';  
23.         [d N] = size(P);   
24. %     else  
25. %         return  
26.     end  
27. end  
28.       
29. %% 选取初始点 方法2   
30. max_Initial = max(20,N/(5*K));  
31. label = zeros(max_Initial,N);  
32. center = zeros(d,K,max_Initial);  
33. C = zeros(1,N);  
34. %% 主循环  
35. for initial_Case = 1:max_Initial  
36.       
37.     pointK = Initial_center(P,K);      
38.     iter = 0;  
39.     max_iter = 1e+3;  
40.     % xK = pointK;  
41.     disp(['------------KM进行第 ' num2str(initial_Case) ' 次重新选择初始中心-----------'])  
42.     %% 每次初始化K个中心点后，进行的循环  
43.     while iter < max_iter  
44.         iter = iter+1;  
45.         if mod(iter,50)==0  
46.             disp(['  内部循环进行第 ' num2str(iter) ' 次迭代'])  
47.         end  
48.         %%%根据数据矩阵P中每个点到中心点的距离(最小)确定所属分类  
49.         for i = 1:N  
50.             dert = repmat(P(:,i),1,K)-pointK;  
51.             distK = sqrt(diag(dert'*dert));  
52.             [~,j] = min(distK);  
53.             C(i) = j;  
54.         end  
55.         %%%重新计算K个中心点    
56.         xK_ = zeros(d,K);  
57.         for i = 1:K  
58.             Pi = P(:,find(C==i));  
59.             Nk = size(Pi,2);  
60.             % K-Means K-Medoids唯一不同的地方：选择中心点的方式  
61.             switch lower(method)  
62.                 case 'kmeans'    
63.                     xK_(:,i) = sum(Pi,2)/Nk;  
64.                 case 'kmedoids'  
65.                     Dx2 = zeros(1,Nk);  
66.                     for t=1:Nk  
67.                        dx = Pi - Pi(:,t)*ones(1,Nk);  
68.                        Dx2(t) = sum(sqrt(sum(dx.*dx,1)),2);  
69.                     end  
70.                     [~,min_ind] = min(Dx2);  
71.                     xK_(:,i) = Pi(:,min_ind);  
72.                 otherwise  
73.                     errordlg('请输入正确的方法：kmeans-OR-kmedoids','MATLAB error');  
74.             end  
75.         end  
76.           
77.         % 判断是否达到结束条件  
78.         if xK_==pointK   % & iter>50  
79.             disp(['###迭代 ' num2str(iter) ' 次得到收敛的解'])  
80.             label(initial_Case,:) = C;  
81.             center(:,:,initial_Case) = xK_;  
82.           % plot_Graph(C);  
83.             break  
84.         end  
85.           
86.         pointK = xK_;  
87.         %xK = xK_;  
88.     end  
89.     if iter == max_iter  
90.          disp('###达到内部最大迭代次数1000,未得到收敛的解')   
91.          label(initial_Case,:) = C;  
92.          center(:,:,initial_Case) = xK_;  
93.         % plot_Graph(C);  
94.          % break  
95.     end  
96.       
97. end  
98.   
99. %%%%增加对聚类结果最优性的比较    
100. %距离差  
101.  dist_N = zeros(max_Initial,K);  
102.  for initial_Case=1:max_Initial       
103.      for k=1:K  
104.          tem = find(label(initial_Case,:)==k);  
105.          dx = P(:,tem)-center(:,k,initial_Case)*ones(1,size(tem,2));  
106.          dxk = sqrt(sum(dx.*dx,1));  
107.          dist_N(initial_Case,k) = sum(dxk);         
108.          % dist_N(initial_Case,k) = dxk;     
109.      end       
110.  end  
111.    
112.  %%%%对于max_Initial次初始化中心点得到的分类错误  
113.  %%%%取错误最小的情况的Label作为最终分类  
114.  dist_N_sum = sum(dist_N,2); %求K类总的误差  
115.  [~,best_ind] = min(dist_N_sum);  
116.  best_Label = label(best_ind,:);   
117.  best_Center = center(:,:,best_ind);  
118.    

选取初始中心：

[cpp]  view plain copy

1. function center = Initial_center(X,K)    
2.         N = size(X,2);  
3.         rnd_Idx = randperm(N);    
4.         center = X(:,rnd_Idx(1:K));    
5. end     

没有对算法精度进行测量，但是可以肯定的是：对于简单的情形，KM算法几乎都可以得到最优的结果。
下面是两张聚类的图：二维

三维：

来源：

http://blog.csdn.net/lifeitengup/article/details/8245028 

References：
1.酷壳网：http://coolshell.cn/articles/7779.html
2.abcjennifer：http://blog.csdn.net/abcjennifer/article/details/8197072该博客的MachineLearning系列很不错。