python 实现 AP近邻传播聚类算法(Affinity Propagation)

Affinity Propagation (AP) 聚类是2007年在Science杂志上提出的一种新的聚类算法。它根据N个数据点之间的相似度进行聚类,这些相似度可以是对称的,即两个数据点互相之间的相似度一样(如欧氏距离);也可以是不对称的,即两个数据点互相之间的相似度不等。这些相似度组成N×N的相似度矩阵S(其中N为有N个数据点)。

AP算法不需要事先指定聚类数目,相反它将所有的数据点都作为潜在的聚类中心,称之为 exemplar。以S矩阵的对角线上的数值s (k, k)作为k点能否成为聚类中心的评判标准,这意味着该值越大,这个点成为聚类中心的可能性也就越大,这个值又称作参考度p ( preference) 。聚类的数量受到参考度p的影响,如果认为每个数据点都有可能作为聚类中心,那么p就应取相同的值。如果取输入的相似度的均值作为p的值,得到聚类数量是中等的。如果取最小值,得到类数较少的聚类。

AP算法中传递两种类型的消息,(responsiility)和(availability) 。r(i,k)表示从点i发送到候选聚类中心k的数值消息,反映k点是否适合作为i点的聚类中心。a(i,k)则从候选聚类中心k发送到i的数值消息,反映i点是否选择k作为其聚类中心。r (i, k)与a (i, k)越强,则k点作为聚类中心的可能性就越大,并且i点隶属于以k点为聚类中心的聚类的可能性也越大。AP算法通过迭代过程不断更新每一个点的吸引度和归属度值,直到产生m个高质量的exemplar,同时将其余的数据点分配到相应的聚类中。

 

在这里介绍几个文中常出现的名词：

exemplar：指的是聚类中心。

similarity：数据点i和点j的相似度记为S(i，j)。是指点j作为点i的聚类中心的相似度。

preference：数据点i的参考度称为P(i)或S(i,i)。是指点i作为聚类中心的参考度。一般取S相似度值的中值。

Responsibility：R(i,k)用来描述点k适合作为数据点i的聚类中心的程度。

Availability：A(i,k)用来描述点i选择点k作为其聚类中心的适合程度。

Damping factor：阻尼系数，主要是起收敛作用的。

AP聚类算法是将每个数据看成图中的一个节点，迭代的过程即是在图中通过传播信 息来找到聚类集合。本文计算两个数据点的相似度采用距离的负数，也就是说距离越近，相似度越大。相似矩阵S中i到j的相似度就是刚刚所说的距离的负数。但是主对角线上的那些数表示的是某个点和自身的相似度，但是这里我们不能直接用0来表示。根据算法要求，主对角线上的值s(k,k)一般称为偏向参数，一般 情况下对所有k，s(k,k)都相等，取非主对角线上的所有数的中位数。这个值很重要，他的大小与最后得到的类的数目有关，一般而言这个数越大，得到的类的数目就越多。

这里为什么要设定一个偏向参数而不直接用0来算呢，估计是因为AP聚类算法是要 用图论的一些东西来理解的，它把所有的点都看成一个图中的节点，通过节点之间的信息传递来达到聚类的效果。具体比较复杂，形象一点说就是我告诉你我和这些人是死党，如果你认为你也是我死党的话，那你就加入我们这一堆人里面来吧！

有一些详细的原理上的东西就不说了，直接说计算过程吧。。聚类就是个不断迭代的过程，迭代的过程主要更新两个矩阵，代表(Responsibility)矩阵R = [r(i,k)]_N×N和适选(Availabilities)矩阵A=[a(i,k)]_N×N。这两个矩阵才初始化为0，N是所有样本的数目。r(i,k)表示第k个样本适合作为第i个样本的类代表点的代表程度，a(i,k)表示第i个样本选择第k个样本作为类代表样本的适合程度。迭代更新公式如下：

每次更新后就可以确定当前样本i的代表样本(exemplar)点k，k就是使{a(i,k)+r(i,k)}取得最大值的那个k，如果i=k的话，那么说明样本i就是自己这个cluster的类代表点，如果不是，那么说明i属于k所属的那个cluster。

当然，迭代停止的条件就是所有的样本的所属都不在变化为止，或者迭代了n次都还没有变化（n的值可以自己取）。

说起来还有一种判断点属于属于哪一类的方法，就是找出所有决策矩阵主对角线元素{a(k,k)+r(k,k)}大于0的所有点，这些点全部都是类代表点，之后在决定其余的点属于这里面的一类。这两种方法的结果我没比较过诶，不知是不是一样的。

另外还有一点就是AP聚类算法迭代过程很容易产生震荡，所以一般每次迭代都加上一个阻尼系数λ：

r_new(i,k) = λ*r_old(i,k) + (1-λ)*r(i,k)

a_new(i,k) = λ*a_old(i,k) + (1-λ)*a(i,k)

[html]  view plain  copy

1. # coding:utf-8  
2.   
3. from sklearn.datasets.samples_generator import make_blobs  
4. import matplotlib.pyplot as plt  
5. import random  
6. import numpy as np  
7.   
8. ##############################################################################  
9. # Generate sample data  
10. centers = [[1, 1], [-1, -1], [1, -1]]  
11. X, labels_true = make_blobs(n_samples=1000, centers=centers, cluster_std=0.4,  
12.                             random_state=0)  
13.   
14. ##############################################################################  
15.   
16.   
17. def euclideanDistance(X, Y):  
18.     """计算每个点与其他所有点之间的欧几里德距离"""  
19.     X = np.array(X)  
20.     Y = np.array(Y)  
21.     # print X  
22.     return np.sqrt(np.sum((X - Y) ** 2))  
23.   
24.   
25.   
26. def computeSimilarity(datalist):  
27.   
28.     num = len(datalist)  
29.   
30.     Similarity = []  
31.     for pointX in datalist:  
32.         dists = []  
33.         for pointY in datalist:  
34.             dist = euclideanDistance(pointX, pointY)  
35.             if dist == 0:  
36.                 dist = 1.5  
37.             dists.append(dist * -1)  
38.         Similarity.append(dists)  
39.   
40.     return Similarity  
41.   
42.   
43. def affinityPropagation(Similarity, lamda):  
44.   
45.     #初始化 吸引矩阵 和 归属 矩阵  
46.     Responsibility = np.zeros_like(Similarity, dtype=np.int)  
47.     Availability = np.zeros_like(Similarity, dtype=np.int)  
48.   
49.     num = len(Responsibility)  
50.   
51.     count = 0  
52.     while count < 10:  
53.         count += 1  
54.         # update 吸引矩阵  
55.   
56.         for Index in range(num):  
57.             # print len(Similarity[Index])  
58.             kSum = [s + a  for s, a in zip(Similarity[Index], Availability[Index])]  
59.             # print kSum  
60.             for Kendex in range(num):  
61.                 kfit = delete(kSum, Kendex)  
62.                 # print fit  
63.                 ResponsibilityNew = Similarity[Index][Kendex] - max(kfit)  
64.                 Responsibility[Index][Kendex] = lamda * Responsibility[Index][Kendex] + (1 - lamda) * ResponsibilityNew  
65.   
66.         # print "Responsibility", Responsibility  
67.   
68.   
69.         # update 归属矩阵  
70.   
71.         ResponsibilityT = Responsibility.T  
72.   
73.         # print ResponsibilityT, Responsibility  
74.   
75.         for Index in range(num):  
76.   
77.             iSum = [r for r in ResponsibilityT[Index]]  
78.   
79.             for Kendex in range(num):  
80.   
81.                 # print Kendex  
82.                 # print "ddddddddddddddddddddddddddd", ResponsibilityT[Kendex]  
83.                 #  
84.                 ifit = delete(iSum, Kendex)  
85.                 ifit = filter(isNonNegative, ifit)   #上面 iSum  已经全部大于0  会导致  delete 下标错误  
86.   
87.                 #   k == K 对角线的情况  
88.                 if Kendex == Index:  
89.                     AvailabilityNew  = sum(ifit)  
90.                 else:  
91.                     result = Responsibility[Kendex][Kendex] + sum(ifit)  
92.                     AvailabilityNew = result if result > 0 else 0  
93.                 Availability[Kendex][Index] = lamda * Availability[Kendex][Index] + (1 - lamda) * AvailabilityNew  
94.         print "###############################################"  
95.         print Responsibility  
96.         print Availability  
97.         print "###############################################"  
98.         return Responsibility + Availability  
99.   
100. def computeCluster(fitable, data):  
101.     clusters = {}  
102.     num = len(fitable)  
103.     for node in range(num):  
104.         fit = list(fitable[node])  
105.         key = fit.index(max(fit))  
106.         if not clusters.has_key(key):  
107.             clusters[key] = []  
108.         point = tuple(data[node])  
109.         clusters[key].append(point)  
110.   
111.     return clusters  
112. ##############################################################################  
113.   
114. """切片删除 返回新数组"""  
115. def delete(lt, index):  
116.     lt = lt[:index] + lt[index+1:]  
117.     return lt  
118.   
119. def isNonNegative(x):  
120.     return x >= 0  
121.   
122.   
123. ##############################################################################  
124.   
125. Similarity = computeSimilarity(X)  
126.   
127. Similarity = np.array(Similarity)  
128.   
129. print "Similarity", Similarity  
130.   
131. fitable = affinityPropagation(Similarity, 0.34)  
132.   
133. print fitable  
134.   
135. clusters = computeCluster(fitable, X)  
136.   
137. # print clusters  
138.   
139. ##############################################################################  
140. clusters = clusters.values()  
141.   
142. print len(clusters)  
143.   
144. ##############################################################################  
145. def plotClusters(clusters, title):  
146.     """ 画图 """  
147.     plt.figure(figsize=(8, 5), dpi=80)  
148.     axes = plt.subplot(111)  
149.     col=[]  
150.     r = lambda: random.randint(0,255)  
151.     for index in range(len(clusters)):  
152.         col.append(('#%02X%02X%02X' % (r(),r(),r())))  
153.     color = 0  
154.     for cluster in clusters:  
155.         cluster = np.array(cluster).T  
156.         axes.scatter(cluster[0],cluster[1], s=20, c = col[color])  
157.         color += 1  
158.     plt.title(title)  
159.     # plt.show()  
160. ##############################################################################  
161. plotClusters(clusters, "clusters by affinity propagation")  
162. plt.show()  
163.   
164. ##############################################################################