AP(Affinity Propagation)亲和力传播算法

原文链接 https://www.cnblogs.com/lc1217/p/6908031.html

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

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] = preturn simi

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 Rdef 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_relse: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_rprint("max_r:"+str(np.max(R)))#print(np.min(R))return R

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 += 0A[i][k] = max3##带入阻尼系数更新A值A[i][k] = (1-lam)*A[i][k] +lam*old_aelse :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] = 0else :A[i][k] = R[k][k] + max4##带入阻尼系数更新A值A[i][k] = (1-lam)*A[i][k] +lam*old_aprint("max_a:"+str(np.max(A)))#print(np.min(A))return A

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 += 1curr_iter += 1print(curr_iter)if curr_iter >= max_iter or curr_comp > max_comp :breakreturn class_cen

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)

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] = preturn 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 Rdef 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_relse: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_rprint("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 += 0A[i][k] = max3##带入阻尼系数更新A值A[i][k] = (1-lam)*A[i][k] +lam*old_aelse :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] = 0else :A[i][k] = R[k][k] + max4##带入阻尼系数更新A值A[i][k] = (1-lam)*A[i][k] +lam*old_aprint("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 += 1curr_iter += 1print(curr_iter)if curr_iter >= max_iter or curr_comp > max_comp :breakreturn class_cenif __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()

迭代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值)取不同值时的聚类中心的数目在代码中注明了。

View Code

6.AP算法的优点

1）不需要制定最终聚类族的个数

2）已有的数据点作为最终的聚类中心，而不是新生成一个族中心。

3）模型对数据的初始值不敏感。

4）对初始相似度矩阵数据的对称性没有要求。

5）.相比与k-centers聚类方法，其结果的平方差误差较小。

7.AP算法的不足

1）AP算法需要事先计算每对数据对象之间的相似度，如果数据对象太多的话，内存放不下，若存在数据库，频繁访问数据库也需要时间。

2）AP算法的时间复杂度较高，一次迭代大概O(N3)

3）聚类的好坏受到参考度和阻尼系数的影响。

AP(Affinity Propagation)亲和力传播算法相关推荐

Affinity propagation 近邻传播算法
终日乾乾,与时偕行. - 周易⋅\cdot乾⋅\cdot文言近邻传播算法是一种基于代表点的聚类方法,它会同时考虑所有数据点都是潜在的代表点,通过结点之间的信息传递,最后得到高质量的聚类.这个信息的传 ...
聚类方法学习（三）亲和力传播算法AP
AP简介亲和力传播算法无需指定聚类的结果,使用的是AffinityPropagation,详情请参考:https://www.cnblogs.com/lc1217/p/6908031.html 算法 ...
近邻传播 Affinity Propagation(AP) 聚类算法原理及实现
1.概述机器学习中 AP(Affinity Propagation)通常被称为近邻传播算法或者密切度传播或类同传播算法,由 Frey 与Dueck于2007年在Science首次提出.AP算法的基本 ...
常见聚类算法及使用--Affinity Propagation（AP）
文章目录前言一些准备亲和力传播(Affinity Propagation) 相似度矩阵(Similarity Matrix) 吸引度矩阵(Responsibility Matrix) 归属度矩阵 ...
（一）反向传播算法理解（Back-Propagation）
本文参考deeplearningbook.org一书第六章 6.5 Back-Propagation and Other Diﬀerentiation Algorithms 反向传播算法分为两篇来讲 ...
AP近邻传播聚类算法(Affinity Propagation)
Affinity Propagation (AP) 聚类是2007年在Science杂志上提出的一种新的聚类算法.它根据N个数据点之间的相似度进行聚类,这些相似度可以是对称的,即两个数据点互相之间的相 ...
Affinity Propagation (AP)近邻传播聚类
近邻传播聚类:根据 N 个数据点之间的相似度聚类,相似度可以是对称的,即两个数据点互相之间的相似度一样(如欧氏距离):也可以是不对称的,即两个数据点互相之间的相似度不等.这些相似度组成 N×N 的相似 ...
Affinity Propagation聚类算法详解
欢迎关注"生信修炼手册"! Affinity Propagation简称AP, 称之为近邻传播算法, 是一种基于图论的聚类算法.将所有样本点看做是一个网络中的节点,图示如下在样本 ...
聚类算法实践——PCCA、SOM、Affinity Propagation
转载自:http://www.itongji.cn/article/0R52E22013.html 这篇日志是这个系列里算法部分的最后一篇,关注的是几个相对另类一点的聚类算法:PCCA.SOM和Aff ...

AP(Affinity Propagation)亲和力传播算法

AP(Affinity Propagation)亲和力传播算法相关推荐

最新文章

热门文章