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：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():
"""
第一步：生成測試數據
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(N3)
3）聚類的好壞受到參考度和阻尼系數的影響。
Python實現聚類算法AP