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

時間  2019-11-06
標籤 機器 學習 python 實現 算法 欄目 Python 简体版
原文   原文鏈接

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)越大，表示兩個點距離越近，類似度也就越高。使用負的歐式距離，類似度是對稱的，若是採用其餘算法，類似度可能就不是對稱的。app

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

     3)  preference：指點i做爲聚類中心的參考度(不能爲0)，取值爲S對角線的值(圖1紅色標註部分)，此值越大，最爲聚類中心的可能性就越大。可是對角線的值爲0，因此須要從新設置對角線的值，既能夠根據實際狀況設置不一樣的值，也能夠設置成同一值。通常設置爲S類似度值的中值。(有的說設置成S的最小值產生的聚類最少，可是在下面的算法中設置成中值產生的聚類是最少的) ide

     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）聚類的好壞受到參考度和阻尼係數的影響。

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程