k-中心點演算法（k-medoids）及Matlab程式碼實現

阿新 • • 發佈：2019-01-07

k-中心點演算法（k-medoids）及Matlab程式碼

1. 提出：

k-means演算法是每次選擇簇的均值作為新的中心點，迭代直到簇中心不再變化（趨於穩定）。其缺點是對離群點特別敏感，因為一個很大的極端值物件會扭曲資料分佈，使簇均值嚴重偏離；於是我們考慮新的簇中心不用均值表示而是選擇簇內的某個物件，只要使總的代價降低就行。。

2. 優化後的演算法：

PAM（圍繞中心點的劃分），具有代表性的k-medoids演算法。

演算法思想： 迭代選出簇中位置最中心的物件，試圖將N個物件給出k個劃分。
具體：

最初隨機選擇k個物件作為中心點，並代表初始簇，然後根據歐氏距離劃分其餘所有物件到各個中心點所代表的簇，得到初始簇劃分。
該演算法反覆利用資料D中所有非代表物件來替換當前代表物件，試圖找出更好的中心點，以改進聚類質量。中心點也叫代表物件，其他物件叫非代表物件。
在每次迭代中，所有可能的物件都被分析，每一對替換中的一個物件是中心點，而另一個是非代表物件。如果一個當前的中心點被一個非代表物件所替換，代價函式將計算平方誤差值所產生的差別；替換的總代價是所有非中心物件去替換所產生的代價之和。
如果總代價為負值，則實際的平方誤差將會減少，則代表物件Oi可被非代表物件Oh替換。

中心點的定義：簇中某點的平均差異性在這一簇中所有點中最小。

3. 演算法描述：

輸入：

k：簇的數目；
N：包含N個物件的資料；

輸出：k個簇，使得所有物件與其距離最近的中心點的相異度總和

最小。

初始化：隨機選擇N個點中的k個點作為初始中心點；
將其餘各點根據歐式距離劃分到這k個類別中；
當損失值減少時：
對於每個中心點o，對於每個非中心點m；
a.)交換o和m，重新計算損失（損失值為，所有點到中心點的距離和）；
b.)如果總的損失增加則不進行交換；

另外一種解釋：
1. 任選k個物件作為初始簇中心；
2. 劃分其餘物件給離它最近的中心點所表示的簇；
3. 選擇一個未被選擇的中心點Oi；
4. 選擇一個未被選擇的非中心點Oh;
5. 計算Oh替換Oi所產生的總代價並記錄在S中；
6. Utill 所有非中心點都被選擇過；
7. Utill 所有中心點都被選擇過；
8. if 在S中的所有替換所產生的總代價有負值存在，then 找出並用該非中心點替換對應的中心點，形成一個新的k箇中心點的集合；
9. Utill 沒有再發生簇的重新分配，即所有的S都大於0.

總的來說：與k均值演算法一樣，初始代表物件任意選取，反覆用一個非代表物件替換一個代表物件，試圖找出更好的中心點，以改進聚類質量。一箇中心點物件被可以產生誤差平方和減少的物件替換，再一次迭代中產生的最佳物件集合成為下次迭代中心點。

為判定一個非代表物件Oh是否是當前一個代表物件Oi的好的替換，對每個非中心點物件Oj，有以下四種情況：

第一種情況：Oj當前隸屬於中心點物件Oi。如果Oi被Oh所代替作為中心點，且Oj離某個中心點Om最近，i≠m，那麼Oj被重新分配給Om。
第二種情況：Oj當前隸屬於中心點物件Oi。如果Oi被Oh所代替作為中心點，且Oj離Oh最近，那麼Oj被重新分配給Oh。
第三種情況：Oj當前隸屬於中心點Om，m≠i。如果Oi被Oh代替作為中心點，而Oj依然離Om最近，那麼物件的隸屬不發生變化。
第四種情況：Oj當前隸屬於中心點Om，m≠i。如果Oi被Oh代替作為一箇中心點，且Oj離Oh最近，那麼Oi被重新分配給Oh。

每當重新分配發生時，絕對誤差E的差會對代價函式有影響。因此，如果一個當前的代表物件被非代表物件所替換時，則代價函式就計算絕對誤差值的差。交換的總代價是所有非代表物件所產生的代價之和。如果總代價為負，則實際的絕對誤差和會減少，Oi可以被Oh所替換作為新中心點。如果代價為正，則本次迭代無變化。

總代價定義如下：
這裡寫圖片描述

其中，TCih表示中心點Oi被非中心點Oh替換後產生的總代價。Cjih表示，Oj在Oi被Oh替換後產生的代價。下面將介紹上述四種情況中代價函式的計算公式，其中引用到的符號有：Oi和Om是兩個原中心點，Oh將替換Oi作為新的中心點。

代價函式計算圖示：
這裡寫圖片描述

Matlab程式碼：

k_medoids.m檔案

clc;  
clear;  
%讀取資料檔案,生成點矩陣  
fileID = fopen('E:\MySoftware\matlabWorks\textKMediods\sample.txt');  
 C=textscan(fileID,'%f %f'); %textscan函式讀取資料
 fclose(fileID); %關閉一個開啟的fileID的檔案
 %顯示陣列結果  
 %celldisp(C);  
 %將cell型別轉換為矩陣型別,這裡只假設原資料為二維屬性，且是二維的座標點  
 CC_init=cat(2,C{1},C{2});%用來儲存初始載入的值  
 CC=CC_init;  %cat函式用於連線兩個矩陣或陣列
  %獲得物件的數量  
 num=length(C{1});  
 %顯示初始分佈圖  
 grid on;%顯示錶格
 scatter(C{1},C{2},'filled');  %filled為實心圓，該函式可以把C中所有座標的點都畫出來。
 %%設定任意k個簇  
k=3;  
%臨時存放k箇中心點的陣列  
C_temp=zeros(k,2);  
%判斷所設定的k值是否小於物件的數量  
if k<num  
    %產生隨機的k個整數  
   randC=randperm(num);  
   randC=randC(1:k);  
   %從原陣列中提出這三個點     
   for i=1:k  
       C_temp(i,:)=CC(randC(1,i),:);  
   end  
   %將原陣列中的這三個點清空  
    for j=1:k  
       CC(randC(1,j),:)=zeros(1,2);  
    end    
    idZero=find(CC(:,1)==0);  
    %刪除為零的行  
    [i1,j1]=find(CC==0);  
    row=unique(i1);  
    CC(row,:)=[];  
   %分配k個二維陣列，用來存放聚類點  
   %分配行為k的儲存單元  
   cluster=cell(k,1,1);   
    %將剔除的三個點加入到對應的三個儲存單元,每個單元的第一行置為0，為了儲存相對應的簇中心  
   for m=1:k  
       cluster{m}=C_temp(m,:);  
   end    
   %計算其他點到這k個點的距離，然後分配這些點，第一次遍歷  
   for i=1:num-k  
       %分別計算到三個點的距離         
       minValue=1000000;%最小值，要根據實際情況設定該值  
       minNum=-1;%最小值序號  
       for j=1:k  
           if minValue>sqrt((CC(i,1)-C_temp(j,1))*(CC(i,1)-C_temp(j,1))+(CC(i,2)-C_temp(j,2))*(CC(i,2)-C_temp(j,2)))  
               minValue=sqrt((CC(i,1)-C_temp(j,1))*(CC(i,1)-C_temp(j,1))+(CC(i,2)-C_temp(j,2))*(CC(i,2)-C_temp(j,2)));  
               minNum=j;  
           end  
       end  
       cluster{minNum}=cat(1,cluster{minNum},CC(i,:));         
   end  
   %隨機選擇p點  
   flag=1;  
   count=0;  
   while flag==1  
       randC=randperm(num-k);  
       randC=randC(1:1);    
       o_random=CC(randC,:);  
       %找出該隨機點所在的簇  
       recordN=0;  
       for i=1:k        
           for j=1:size(cluster{i},1)        
               cc=cluster{i}(j,:);  
               cc=cc-o_random;  
               if cc==0  
                   recordN=i;  
                   break;  
               end  
           end  
       end  
       %將選擇的隨機點從點集中刪除  
       CC(randC,:)=[];  
       %計算替換代價  
       o=cluster{recordN}(1,:);  
       o_rand_sum=0;  
       o_sum=0;  
       for i=1:length(CC)  
           o_rand_sum=o_rand_sum+sqrt((CC(i,1)-o_random(1,1))*(CC(i,1)-o_random(1,1))+(CC(i,2)-o_random(1,2))*(CC(i,2)-o_random(1,2)));  
           o_sum=o_sum+sqrt((CC(i,1)-o(1,1))*(CC(i,1)-o(1,1))+(CC(i,2)-o(1,2))*(CC(i,2)-o(1,2)));  
       end  
       %如果隨機選擇的點的代價小於原始代表點的代價，則替換該代表點，然後重新聚類  
       if o_rand_sum<o_sum  
           cluster{recordN}(1,:)=o_random;  
           %將代表點放入物件集  
           CC=cat(1,CC,o);  
           %對所有物件重新進行聚類  
           %將cluster除第一行之外的資料全部清空  
           for i=1:k  
               c=cluster{i}(1,:);  
               cluster{i}=[];  
               cluster{i}=c;  
           end   
           %重新聚類  
           for i=1:num-k  
               %分別計算到三個點的距離         
               minValue=1000000;%最小值，要根據實際情況設定該值  
               minNum=-1;%最小值序號  
               for j=1:k  
                   if minValue>sqrt((CC(i,1)-C_temp(j,1))*(CC(i,1)-C_temp(j,1))+(CC(i,2)-C_temp(j,2))*(CC(i,2)-C_temp(j,2)))  
                       minValue=sqrt((CC(i,1)-C_temp(j,1))*(CC(i,1)-C_temp(j,1))+(CC(i,2)-C_temp(j,2))*(CC(i,2)-C_temp(j,2)));  
                       minNum=j;  
                   end  
               end  
               cluster{minNum}=cat(1,cluster{minNum},CC(i,:));         
           end             
       else  
           %將隨機點重新放入物件集  
           CC=cat(1,CC,o_random);  
           %終止迴圈  
           flag=0;  
       end  
       count=count+1;  
   end     
   %繪製聚類結果  
   for i=1:k  
       scatter(cluster{i}(:,1),cluster{i}(:,2),'filled');  
       hold on  
       grid on;%顯示錶格
   end     
end

執行結果如下：

當K=3時：
這裡寫圖片描述

當K=4時：

這裡寫圖片描述

參考文獻及例子：

歡迎掃描提問碼交流
CSDN程式碼下載地址：

k-中心點演算法（k-medoids）及Matlab程式碼實現