Python實現聚類算法AP

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：R_t+1(i,k)表示新的R(i,k)，R_t(i,k)表示舊的R(i,k)，也許這樣說更容易理解。其中λ是阻尼係數，取值[0.5,1)，用於算法的收斂

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

             （公式二）

          sklearn官網的公式是：

                 （公式三）

            我試了這兩種公式以後，發現仍是公式一的聚類效果最好。一樣的數據都採起S的中值做爲參考度，我本身寫的算法聚類中心是5個，sklearn提供的算法聚類中心是十三個，可是若是把參考度設置爲p=-50，則我本身寫的算法聚類中心不少，sklearn提供的聚類算法產生標準的3個聚類中心(由於數據是圍繞三個中心點產生的)，目前還不清楚這個p=-50是怎麼獲得的。

    2）歸屬度迭代公式

      

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

4.算法流程

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

def init_sample():
    """
    第一步：生成測試數據
        1.生成實際中心爲centers的測試樣本300個，
        2.Xn是包含150個(x,y)點的二維數組
        3.labels_true爲其對應的真是類別標籤
    """
    # 生成的測試數據的中心點
    centers = [[1, 1], [-1, -1], [1, -1]]
    # 生成數據
    X, label_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5, random_state=0)
    return X, label_true

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

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

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

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

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

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

完整代碼

# -*- coding: utf-8 -*-

"""
@Datetime: 2019/3/31
@Author: Zhang Yafei
"""
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets.samples_generator import make_blobs


def init_sample():
    """
    第一步：生成測試數據
        1.生成實際中心爲centers的測試樣本300個，
        2.Xn是包含150個(x,y)點的二維數組
        3.labels_true爲其對應的真是類別標籤
    """
    # 生成的測試數據的中心點
    centers = [[1, 1], [-1, -1], [1, -1]]
    # 生成數據
    X, label_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5, random_state=0)
    return X, label_true


class AP(object):
    """ AP聚類 """

    def __init__(self):
        self.Xn = None
        self.Xn_len = None
        self.R = None
        self.A = None
        self.simi_matrix = None
        self.class_cen = None

    def fit(self, data):
        self.Xn = data
        self.Xn_len = len(data)
        # 初始化R、A矩陣
        self.R = np.zeros((self.Xn_len, self.Xn_len))
        self.A = np.zeros((self.Xn_len, self.Xn_len))
        # 計算類似度
        self.cal_simi()
        # 輸出聚類中心
        self.class_cen = self.cal_cls_center()

    def cal_simi(self):
        """
        計算類似度矩陣
            這個數據集的類似度矩陣，最終是二維數組
            每一個數字與全部數字的類似度列表，即矩陣中的一行
            採用負的歐式距離計算類似度
        :return:
        """
        simi = [[-np.sqrt((m[0] - n[0]) ** 2 + (m[1] - n[1]) ** 2) for n in self.Xn] for m in self.Xn]

        # 設置參考度，即對角線的值，通常爲最小值或者中值
        # p = np.min(simi)   ##11箇中心
        # p = np.max(simi)  ##14箇中心

        p = np.median(simi)  ##5箇中心
        for i in range(self.Xn_len):
            simi[i][i] = p

        self.simi_matrix = simi

    def iter_update_R(self, old_r=0, lam=0.5):
        """
        計算吸引度矩陣，即R
               公式1：r(n+1) =s(n)-(s(n)+a(n))-->簡化寫法，具體參見上圖公式
               公式2：r(n+1)=(1-λ)*r(n+1)+λ*r(n)
        迭代更新R矩陣
        :param old_r: 更新前的某個r值
        :param lam: 阻尼係數，用於算法收斂
        :return:
        """
        # 此循環更新R矩陣
        for i in range(self.Xn_len):
            for k in range(self.Xn_len):
                old_r = self.R[i][k]
                if i != k:
                    max1 = self.A[i][0] + self.R[i][0]  ##注意初始值的設置
                    for j in range(self.Xn_len):
                        if j != k:
                            if self.A[i][j] + self.R[i][j] > max1:
                                max1 = self.A[i][j] + self.R[i][j]
                    ##更新後的R[i][k]值
                    self.R[i][k] = self.simi_matrix[i][k] - max1
                    ##帶入阻尼係數從新更新
                    self.R[i][k] = (1 - lam) * self.R[i][k] + lam * old_r
                else:
                    max2 = self.simi_matrix[i][0]  ##注意初始值的設置
                    for j in range(self.Xn_len):
                        if j != k:
                            if self.simi_matrix[i][j] > max2:
                                max2 = self.simi_matrix[i][j]
                    ##更新後的R[i][k]值
                    self.R[i][k] = self.simi_matrix[i][k] - max2
                    ##帶入阻尼係數從新更新
                    self.R[i][k] = (1 - lam) * self.R[i][k] + lam * old_r
        print("max_r:" + str(np.max(self.R)))

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

                    # 上圖公式中的min部分
                    if self.R[k][k] + max4 > 0:
                        self.A[i][k] = 0
                    else:
                        self.A[i][k] = self.R[k][k] + max4

                    # 帶入阻尼係數更新A值
                    self.A[i][k] = (1 - lam) * self.A[i][k] + lam * old_a
        print("max_a:" + str(np.max(self.A)))

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

    def c_list(self):
        # 根據聚類中心劃分數據
        c_list = []
        for m in self.Xn:
            temp = []
            for j in self.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)
        print(c_list)
        return c_list


def plot(class_cen, X, c_list):
    # 畫圖
    colors = ['red', 'blue', 'black', 'green', 'yellow']
    plt.figure(figsize=(8, 6))
    plt.xlim([-3, 3])
    plt.ylim([-3, 3])
    for i in range(len(X)):
        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.savefig('AP 聚類.png')
    plt.show()


if __name__ == '__main__':
    # 初始化數據
    Xn, labels_true = init_sample()
    ap = AP()
    ap.fit(data=Xn)
    class_cen = ap.class_cen
    # for i in class_cen:
    #    print(str(i)+":"+str(Xn[i]))
    c_list = ap.c_list()
    plot(class_cen=class_cen, X=Xn, c_list=c_list)

AP.py

效果圖

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值)取不一樣值時的聚類中心的數目在代碼中註明了。

# -*- coding: utf-8 -*-

"""
@Datetime: 2019/3/31
@Author: Zhang Yafei
"""

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


def init_sample():
    """
    第一步：生成測試數據
        1.生成實際中心爲centers的測試樣本300個，
        2.Xn是包含150個(x,y)點的二維數組
        3.labels_true爲其對應的真是類別標籤
    """
    # 生成的測試數據的中心點
    centers = [[1, 1], [-1, -1], [1, -1]]
    # 生成數據
    X, label_true = make_blobs(n_samples=150, centers=centers, cluster_std=0.5, random_state=0)
    return X, label_true


def simi_matrix(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)
    return simi


if __name__ == '__main__':
    Xn, label_true = init_sample()
    simi_matrix = simi_matrix(Xn)

    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_
    print(ap.labels_)
    for idx in cluster_centers_indices:
        print(Xn[idx]

sklearn_AP.py

 6.AP算法的優勢

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

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

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

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

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

7.AP算法的不足

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

    2）AP算法的時間複雜度較高，一次迭代大概O(N³)

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