機器學習：Python實現聚類算法(二)之AP算法

分類：IT技術 時間:2017-05-27

1.算法簡介

        AP(Affinity Propagation)通常被翻譯為近鄰傳播算法或者親和力傳播算法，是在2007年的Science雜誌上提出的一種新的聚類算法。AP算法的基本思想是將全部數據點都當作潛在的聚類中心(稱之為exemplar)，然後數據點兩兩之間連線構成一個網絡(相似度矩陣)，再通過網絡中各條邊的消息(responsibility和availability)傳遞計算出各樣本的聚類中心。

 

2.相關概念(假如有數據點i和數據點j)

         

                         （圖1）                                     （圖2）                                       （圖3）

     1）相似度： 點j作為點i的聚類中心的能力，記為S(i,j)。一般使用負的歐式距離，所以S(i,j)越大，表示兩個點距離越近，相似度也就越高。使用負的歐式距離，相似度是對稱的，如果采用其他算法，相似度可能就不是對稱的。

     2）相似度矩陣：N個點之間兩兩計算相似度，這些相似度就組成了相似度矩陣。如圖1所示的黃色區域，就是一個5*5的相似度矩陣(N=5)

     3)  preference：指點i作為聚類中心的參考度(不能為0)，取值為S對角線的值(圖1紅色標註部分)，此值越大，最為聚類中心的可能性就越大。但是對角線的值為0，所以需要重新設置對角線的值，既可以根據實際情況設置不同的值，也可以設置成同一值。一般設置為S相似度值的中值。(有的說設置成S的最小值產生的聚類最少，但是在下面的算法中設置成中值產生的聚類是最少的) 

     4）Responsibility(吸引度):指點k適合作為數據點i的聚類中心的程度，記為r(i,k)。如圖2紅色箭頭所示，表示點i給點k發送信息，是一個點i選點k的過程。

     5）Availability(歸屬度):指點i選擇點k作為其聚類中心的適合程度，記為a(i,k)。如圖3紅色箭頭所示，表示點k給點i發送信息，是一個點k選diani的過程。

     6）exemplar：指的是聚類中心。

     7）r (i, k)加a (i, k)越大,則k點作為聚類中心的可能性就越大,並且i點隸屬於以k點為聚類中心的聚類的可能性也越大

 

3.數學公式

    1）吸引度叠代公式:

             （公式一）

        說明1：Rt+1(i,k)表示新的R(i,k)，Rt(i,k)表示舊的R(i,k)，也許這樣說更容易理解。其中λ是阻尼系數，取值[0.5,1)，用於算法的收斂

        說明2：網上還有另外一種數學公式：

             （公式二）

          sklearn官網的公式是：

                 （公式三）

            我試了這兩種公式之後，發現還是公式一的聚類效果最好。同樣的數據都采取S的中值作為參考度，我自己寫的算法聚類中心是5個，sklearn提供的算法聚類中心是十三個，但是如果把參考度設置為p=-50，則我自己寫的算法聚類中心很多，sklearn提供的聚類算法產生標準的3個聚類中心(因為數據是圍繞三個中心點產生的)，目前還不清楚這個p=-50是怎麽得到的。

    2）歸屬度叠代公式

      

       說明：At+1(i,k)表示新的A(i,k)，At(i,k)表示舊的A(i,k)。其中λ是阻尼系數，取值[0.5,1)，用於算法的收斂

 

4.詳細的算法流程

   1）設置實驗數據。使用sklearn包中提供的函數，隨機生成以[1, 1], [-1, -1], [1, -1]三個點為中心的150個數據。    

def init_sample():
    ## 生成的測試數據的中心點
    centers = [[1, 1], [-1, -1], [1, -1]]
    ##生成數據
    Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                            random_state=0)
    #3數據的長度，即：數據點的個數
    dataLen = len(Xn)

    return Xn,dataLen
View Code

   2）計算相似度矩陣，並且設置參考度，這裏使用相似度矩陣的中值

def cal_simi(Xn):
    ##這個數據集的相似度矩陣，最終是二維數組
    simi = []
    for m in Xn:
        ##每個數字與所有數字的相似度列表，即矩陣中的一行
        temp = []
        for n in Xn:
            ##采用負的歐式距離計算相似度
            s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
            temp.append(s)
        simi.append(temp)

    ##設置參考度，即對角線的值，一般為最小值或者中值
    #p = np.min(simi)   ##11個中心
    #p = np.max(simi)  ##14個中心
    p = np.median(simi)  ##5個中心
    for i in range(dataLen):
        simi[i][i] = p
    return simi
View Code

   3）計算吸引度矩陣，即R值。

     如果有細心的同學會發現，在上述求R和求A的公式中，求R需要A，求A需要R，所以R或者A不是一開始就可以求解出的，需要先初始化，然後再更新。(我開始就陷入了這個誤區，總覺得公式有問題，囧)

##初始化R矩陣、A矩陣
def init_R(dataLen):
    R = [[0]*dataLen for j in range(dataLen)] 
    return R

def init_A(dataLen):
    A = [[0]*dataLen for j in range(dataLen)]
    return A

##叠代更新R矩陣
def iter_update_R(dataLen,R,A,simi):
    old_r = 0 ##更新前的某個r值
    lam = 0.5 ##阻尼系數,用於算法收斂
    ##此循環更新R矩陣
    for i in range(dataLen):
        for k in range(dataLen):
            old_r = R[i][k]
            if i != k:
                max1 = A[i][0] + R[i][0]  ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if A[i][j] + R[i][j] > max1 :
                            max1 = A[i][j] + R[i][j]
                ##更新後的R[i][k]值
                R[i][k] = simi[i][k] - max1
                ##帶入阻尼系數重新更新
                R[i][k] = (1-lam)*R[i][k] +lam*old_r
            else:
                max2 = simi[i][0] ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if simi[i][j] > max2:
                            max2 = simi[i][j]
                ##更新後的R[i][k]值
                R[i][k] = simi[i][k] - max2
                ##帶入阻尼系數重新更新
                R[i][k] = (1-lam)*R[i][k] +lam*old_r
    print("max_r:"+str(np.max(R)))
    #print(np.min(R))
    return R
View Code

  4）計算歸屬度矩陣，即A值

##叠代更新A矩陣
def iter_update_A(dataLen,R,A):
    old_a = 0 ##更新前的某個a值
    lam = 0.5 ##阻尼系數,用於算法收斂
    ##此循環更新A矩陣
    for i in range(dataLen):
        for k in range(dataLen):
            old_a = A[i][k]
            if i ==k :
                max3 = R[0][k] ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if R[j][k] > 0:
                            max3 += R[j][k]
                        else :
                            max3 += 0
                A[i][k] = max3
                ##帶入阻尼系數更新A值
                A[i][k] = (1-lam)*A[i][k] +lam*old_a
            else :
                max4 = R[0][k] ##註意初始值的設置
                for j in range(dataLen):
                    ##上圖公式中的i!=k 的求和部分
                    if j != k and j != i:
                        if R[j][k] > 0:
                            max4 += R[j][k]
                        else :
                            max4 += 0

                ##上圖公式中的min部分
                if R[k][k] + max4 > 0:
                    A[i][k] = 0
                else :
                    A[i][k] = R[k][k] + max4
                    
                ##帶入阻尼系數更新A值
                A[i][k] = (1-lam)*A[i][k] +lam*old_a
    print("max_a:"+str(np.max(A)))
    #print(np.min(A))
    return A
View Code

 5）叠代更新R值和A值。終止條件是聚類中心在一定程度上不再更新或者達到最大叠代次數

##計算聚類中心
def cal_cls_center(dataLen,simi,R,A):
    ##進行聚類，不斷叠代直到預設的叠代次數或者判斷comp_cnt次後聚類中心不再變化
    max_iter = 100    ##最大叠代次數
    curr_iter = 0     ##當前叠代次數
    max_comp = 30     ##最大比較次數
    curr_comp = 0     ##當前比較次數
    class_cen = []    ##聚類中心列表，存儲的是數據點在Xn中的索引
    while True:
        ##計算R矩陣
        R = iter_update_R(dataLen,R,A,simi)
        ##計算A矩陣
        A = iter_update_A(dataLen,R,A)
        ##開始計算聚類中心
        for k in range(dataLen):
            if R[k][k] +A[k][k] > 0:
                if k not in class_cen:
                    class_cen.append(k)
                else:
                    curr_comp += 1
        curr_iter += 1
        print(curr_iter)
        if curr_iter >= max_iter or curr_comp > max_comp :
            break
    return class_cen
View Code

 6）根據求出的聚類中心，對數據進行分類

     這個步驟產生的是一個歸類列表，列表中的每個數字對應著樣本數據中對應位置的數據的分類

 ##根據聚類中心劃分數據
    c_list = []
    for m in Xn:
        temp = []
        for j in class_cen:
            n = Xn[j]
            d = -np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
            temp.append(d)
        ##按照是第幾個數字作為聚類中心進行分類標識
        c = class_cen[temp.index(np.max(temp))]
        c_list.append(c)
View Code

  7）完整代碼及效果圖

from sklearn.datasets.samples_generator import make_blobs
import numpy as np
import matplotlib.pyplot as plt
'''
第一步：生成測試數據
    1.生成實際中心為centers的測試樣本300個，
    2.Xn是包含150個(x,y)點的二維數組
    3.labels_true為其對應的真是類別標簽
'''

def init_sample():
    ## 生成的測試數據的中心點
    centers = [[1, 1], [-1, -1], [1, -1]]
    ##生成數據
    Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                            random_state=0)
    #3數據的長度，即：數據點的個數
    dataLen = len(Xn)

    return Xn,dataLen

'''
第二步：計算相似度矩陣
'''
def cal_simi(Xn):
    ##這個數據集的相似度矩陣，最終是二維數組
    simi = []
    for m in Xn:
        ##每個數字與所有數字的相似度列表，即矩陣中的一行
        temp = []
        for n in Xn:
            ##采用負的歐式距離計算相似度
            s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
            temp.append(s)
        simi.append(temp)

    ##設置參考度，即對角線的值，一般為最小值或者中值
    #p = np.min(simi)   ##11個中心
    #p = np.max(simi)  ##14個中心
    p = np.median(simi)  ##5個中心
    for i in range(dataLen):
        simi[i][i] = p
    return simi

'''
第三步：計算吸引度矩陣，即R
       公式1：r(n+1) =s(n)-(s(n)+a(n))-->簡化寫法，具體參見上圖公式
       公式2：r(n+1)=(1-λ)*r(n+1)+λ*r(n)
'''

##初始化R矩陣、A矩陣
def init_R(dataLen):
    R = [[0]*dataLen for j in range(dataLen)] 
    return R

def init_A(dataLen):
    A = [[0]*dataLen for j in range(dataLen)]
    return A

##叠代更新R矩陣
def iter_update_R(dataLen,R,A,simi):
    old_r = 0 ##更新前的某個r值
    lam = 0.5 ##阻尼系數,用於算法收斂
    ##此循環更新R矩陣
    for i in range(dataLen):
        for k in range(dataLen):
            old_r = R[i][k]
            if i != k:
                max1 = A[i][0] + R[i][0]  ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if A[i][j] + R[i][j] > max1 :
                            max1 = A[i][j] + R[i][j]
                ##更新後的R[i][k]值
                R[i][k] = simi[i][k] - max1
                ##帶入阻尼系數重新更新
                R[i][k] = (1-lam)*R[i][k] +lam*old_r
            else:
                max2 = simi[i][0] ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if simi[i][j] > max2:
                            max2 = simi[i][j]
                ##更新後的R[i][k]值
                R[i][k] = simi[i][k] - max2
                ##帶入阻尼系數重新更新
                R[i][k] = (1-lam)*R[i][k] +lam*old_r
    print("max_r:"+str(np.max(R)))
    #print(np.min(R))
    return R
'''
    第四步：計算歸屬度矩陣，即A
'''
##叠代更新A矩陣
def iter_update_A(dataLen,R,A):
    old_a = 0 ##更新前的某個a值
    lam = 0.5 ##阻尼系數,用於算法收斂
    ##此循環更新A矩陣
    for i in range(dataLen):
        for k in range(dataLen):
            old_a = A[i][k]
            if i ==k :
                max3 = R[0][k] ##註意初始值的設置
                for j in range(dataLen):
                    if j != k:
                        if R[j][k] > 0:
                            max3 += R[j][k]
                        else :
                            max3 += 0
                A[i][k] = max3
                ##帶入阻尼系數更新A值
                A[i][k] = (1-lam)*A[i][k] +lam*old_a
            else :
                max4 = R[0][k] ##註意初始值的設置
                for j in range(dataLen):
                    ##上圖公式中的i!=k 的求和部分
                    if j != k and j != i:
                        if R[j][k] > 0:
                            max4 += R[j][k]
                        else :
                            max4 += 0

                ##上圖公式中的min部分
                if R[k][k] + max4 > 0:
                    A[i][k] = 0
                else :
                    A[i][k] = R[k][k] + max4
                    
                ##帶入阻尼系數更新A值
                A[i][k] = (1-lam)*A[i][k] +lam*old_a
    print("max_a:"+str(np.max(A)))
    #print(np.min(A))
    return A

'''
   第5步：計算聚類中心
'''

##計算聚類中心
def cal_cls_center(dataLen,simi,R,A):
    ##進行聚類，不斷叠代直到預設的叠代次數或者判斷comp_cnt次後聚類中心不再變化
    max_iter = 100    ##最大叠代次數
    curr_iter = 0     ##當前叠代次數
    max_comp = 30     ##最大比較次數
    curr_comp = 0     ##當前比較次數
    class_cen = []    ##聚類中心列表，存儲的是數據點在Xn中的索引
    while True:
        ##計算R矩陣
        R = iter_update_R(dataLen,R,A,simi)
        ##計算A矩陣
        A = iter_update_A(dataLen,R,A)
        ##開始計算聚類中心
        for k in range(dataLen):
            if R[k][k] +A[k][k] > 0:
                if k not in class_cen:
                    class_cen.append(k)
                else:
                    curr_comp += 1
        curr_iter += 1
        print(curr_iter)
        if curr_iter >= max_iter or curr_comp > max_comp :
            break
    return class_cen
  
   
if __name__=='__main__':
    ##初始化數據
    Xn,dataLen = init_sample()
    ##初始化R、A矩陣
    R = init_R(dataLen)
    A = init_A(dataLen)
    ##計算相似度
    simi = cal_simi(Xn)   
    ##輸出聚類中心
    class_cen = cal_cls_center(dataLen,simi,R,A)
    #for i in class_cen:
    #    print(str(i)+":"+str(Xn[i]))
    #print(class_cen)

    ##根據聚類中心劃分數據
    c_list = []
    for m in Xn:
        temp = []
        for j in class_cen:
            n = Xn[j]
            d = -np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
            temp.append(d)
        ##按照是第幾個數字作為聚類中心進行分類標識
        c = class_cen[temp.index(np.max(temp))]
        c_list.append(c)
    ##畫圖
    colors = ['red','blue','black','green','yellow']
    plt.figure(figsize=(8,6))
    plt.xlim([-3,3])
    plt.ylim([-3,3])
    for i in range(dataLen):
        d1 = Xn[i]
        d2 = Xn[c_list[i]]
        c = class_cen.index(c_list[i])
        plt.plot([d2[0],d1[0]],[d2[1],d1[1]],color=colors[c],linewidth=1)
        #if i == c_list[i] :
        #    plt.scatter(d1[0],d1[1],color=colors[c],linewidth=3)
        #else :
        #    plt.scatter(d1[0],d1[1],color=colors[c],linewidth=1)
    plt.show()
View Code

 叠代11次出結果：

       補充說明：這個算法重點在講解實現過程，執行效率不是特別高，有優化的空間。以後我會補充進來

 

5.sklearn包中的AP算法

    1）函數：sklearn.cluster.AffinityPropagation

    2）主要參數：

         damping : 阻尼系數，取值[0.5,1)

         convergence_iter ：比較多少次聚類中心不變之後停止叠代，默認15

         max_iter ：最大叠代次數

         preference :參考度

    3）主要屬性

        cluster_centers_indices_ : 存放聚類中心的數組

        labels_ :存放每個點的分類的數組

        n_iter_ : 叠代次數

    4）示例     

        preference(即p值)取不同值時的聚類中心的數目在代碼中註明了。

from sklearn.cluster import AffinityPropagation
from sklearn import metrics
from sklearn.datasets.samples_generator import make_blobs
import numpy as np


## 生成的測試數據的中心點
centers = [[1, 1], [-1, -1], [1, -1]]
##生成數據
Xn, labels_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5,
                            random_state=0)



simi = []
for m in Xn:
    ##每個數字與所有數字的相似度列表，即矩陣中的一行
    temp = []
    for n in Xn:
         ##采用負的歐式距離計算相似度
        s =-np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2)
        temp.append(s)
    simi.append(temp)

p=-50   ##3個中心
#p = np.min(simi)  ##9個中心，
#p = np.median(simi)  ##13個中心    

ap = AffinityPropagation(damping=0.5,max_iter=500,convergence_iter=30,
                         preference=p).fit(Xn)
cluster_centers_indices = ap.cluster_centers_indices_

for idx in cluster_centers_indices:
    print(Xn[idx])
View Code

 

 6.AP算法的優點

    1） 不需要制定最終聚類族的個數 

    2） 已有的數據點作為最終的聚類中心，而不是新生成一個族中心。 

    3）模型對數據的初始值不敏感。 

    4）對初始相似度矩陣數據的對稱性沒有要求。 

    5）.相比與k-centers聚類方法，其結果的平方差誤差較小。

 

7.AP算法的不足

    1）AP算法需要事先計算每對數據對象之間的相似度，如果數據對象太多的話，內存放不下，若存在數據庫，頻繁訪問數據庫也需要時間。

    2）AP算法的時間復雜度較高，一次叠代大概O(N3)

    3）聚類的好壞受到參考度和阻尼系數的影響。

 


Tags: 對角線 可能性 親和力 黃色 能力

文章來源：

ads
ads
相關文章
ads

相關文章

ad