机器学习：Python实现聚类算法（二）之AP算法

机器学习：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():
## ⽣成的测试数据的中⼼点
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
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值。

终⽌条件是聚类中⼼在⼀定程度上不再更新或者达到最⼤迭代次数
##进⾏聚类，不断迭代直到预设的迭代次数或者判断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)点的⼆维数组
bels_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）聚类的好坏受到参考度和阻尼系数的影响。