scikit-learn 源码解读之Affinity Propagation聚类 2016-04-17 # scikit-learn 源码解读之Affinity Propagation聚类 [TOC] 顺着scikit-learn 文档往下读，上一篇[scikit-learn 源码解读之Kmeans](http://midday.me/article/f8d29baa83ae41ec8c9826401eb7685e)解读的是Kmeans 。Affinity Propagation 算法(下文简写AP算法)是2007年发表在Science 上的，可以在[这里](http://www.psi.toronto.edu/index.php?q=affinity%20propagation)找到作者提供的论文，和matlab，R的实现。资了很详细。本文从scikit-learn 的实现解读其中的原理，和实现细节。 ## 原理 AP算法相对于Kmeans ,的优势是不需要指定聚类数量，对初始值的不敏感。 AP 算法用样本间的相似矩阵作为输入，这个相似矩阵可以是对称或不对称的，scikit-learn 的实现中默认会采用欧式距离来计算相似矩阵，如果需要采用其他方式计算这个相似矩阵就必须提取计算好相似矩阵作为算法的输入。相似矩阵对角线上的元素叫preference 它描述的每个样本作为聚类中心的适合程度，通常是取相似矩阵的中值意思就是所有样本作为中心的可能性是相等的，同时preference 也会影响聚类数量的多少，越小聚类数就会相对较少，反之也成立。 AP算法也被翻译成“吸引子传播算法”，我想主要是因为算法本身主要考虑样本间的消息传递，主要有两种： 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个高质量的聚类中心同时将其余的数据点分配到相应的类中。 *如果要真正理解算法的原理我还离得很远，因为还不太能理解作者为何要量化r(i,k)和a(i,k)这其中有什么物理含义还是相对费解的。* 下面是a(i,k) 和r(i,k) 的定义，以及迭代过程. 看数学公式是不复杂的，复杂的是为什么要这样设计..........如果你知道请一定邮件给我（在首页点头像能找到我的邮箱） $$ r(i,k) = S(i,k)-max(a(i,\hat{k}),S(i, \hat{k}))) \qquad k\neq{\hat{k} } \quad (1)$$ $$a(i,k)=min(0,r(k,k)+\sum(max(0,\hat{i},k))) \qquad\hat{i}\neq{k,i} \qquad (2)$$ $$a(k,k)=max(0,r(\hat{i},k)) \qquad \hat{i}\neq{k} \qquad(3)$$ 上面就是如何量化的公式，其中S就是相似矩阵。通过下面的迭代来实现不断的更新，R,和A由于编辑器对公式的支持感觉还有Bug在这先贴个图片  通过判断迭代过程的结果是否一直没有发生变化，如连续15次迭代结果都没变化就退出，和大于最大迭代次数两个条件来结束迭代。最后选择R+A>0 的样本作为聚类中心。 很多地方是可以自行修改的，比如preference,和距离矩阵都可以给定权重。 ## 算法实现 ### 分析 看完Paper 要实现算法其实最大的难度是将上面的公式很好的转化成程序的语言，这个转化过程决定了你最终代码的样子，上面的公式定义的是两个样本之间的信息度量，实际上就是矩阵中一个元素值的定义，而且整个公式也可以用矩阵计算，前提是要理解清楚。在这种时候才是体现数学和编程完美结合的地方，如果数学不行，你写的代码绝对是不能看的。下面这段代码是我在看scikit-learn 源码之前写的，奇丑无比。 ``` def affinity_propagation(X, max_iter=300, lam=0.6): S = euclidean_distances(X,X,squared=True) print S n_sample = S.shape[0] A_t1 = np.zeros((n_sample, n_sample)) A_t2 = np.zeros((n_sample, n_sample)) A_t3 = np.zeros((n_sample, n_sample)) R_t1 = np.zeros((n_sample, n_sample)) R_t2 = np.zeros((n_sample, n_sample)) R_t3 = np.zeros((n_sample, n_sample)) for step in range(max_iter): for i in range(n_sample): for j in range(n_sample): sum_temp = [] for k in range(n_sample): if k==j: continue sum_temp.append(A[i,k]+S[i,k]) R_t2[i][j] = S[i][j]-max(sum_temp) R_t3 = (1-lam)*R_t2 + lam*R_t1 R_t1 = R_t2 for i in range(n_sample) : for j in range(n_sample): if i==j: sum_j = 0 for t in range(n_sample): if t==j: continue sum_j += max(0, R_t3[t,j]) A_t2[i][j] = sum_j continue sum_j = 0 for t in range(n_sample): if t==i or t==j: continue sum_j += max(0,R_t3[t,j]) A_t2[i][j] = min(0, sum_j+R_t3[j][j]) ``` 基本实现了迭代的功能，但是看着惨不忍睹。所以有人问数学和编程有多大关系，现在我想说是的是，你的数学能力会决定你写什么层次的程序，有些片面，如果有人能很好的把公式转化成程序语言，实现其实是很简单的，看scikit-learn 源码整个实现总共326行，还有大部分注释。 ### 源码解读 #### 用例 ``` from sklearn.cluster import AffinityPropagation from sklearn.datasets.samples_generator import make_blobs centers = [[1, 1], [-1, -1], [1, -1]] X, labels_true = make_blobs(n_samples=300, centers=centers, cluster_std=0.5, random_state=0) af = AffinityPropagation(preference=-50).fit(X) cluster_centers_indices = af.cluster_centers_indices_ labels = af.labels_ n_clusters_ = len(cluster_centers_indices) ``` 第二行利用scikit-learn 提供的工具产生了测试数据，是有300个样本以3个中心分布的随机数 接口提供了AffinityPropagation 类，fit 方法用于训练数据，这里采用了默认的距离计算方式，即欧式距离，如果要采用其他的计算方式需要自行算好S,然后设置参数affinity=precomputed就可以。 #### 代码结构 在scikit-learn 项目cluster 下的affinity_propagation.py文件中， ``` class AffinityPropagation(BaseEstimator, ClusterMixin) ``` 接口类一如既往的继承了两个基础的Class， 在fit 方法里会把计算距离矩阵赋值给**affinity_matrix_**变量, 最终的聚类算法在**affinity_propagation**中实现，所已只需要看这个function 就能搞清楚所有实现的细节 #### 核心注释 先贴上源代码后面会详细解读，为了排版注释和代码是分开的。 ``` def affinity_propagation(S, preference=None, convergence_iter=15, max_iter=200, damping=0.5, copy=True, verbose=False, return_n_iter=False): S = as_float_array(S, copy=copy) n_samples = S.shape[0] if S.shape[0] != S.shape[1]: raise ValueError("S must be a square array (shape=%s)" % repr(S.shape)) if preference is None: preference = np.median(S) if damping < 0.5 or damping >= 1: raise ValueError('damping must be >= 0.5 and < 1') random_state = np.random.RandomState(0) # Place preference on the diagonal of S S.flat[::(n_samples + 1)] = preference A = np.zeros((n_samples, n_samples)) R = np.zeros((n_samples, n_samples)) # Initialize messages # Intermediate results tmp = np.zeros((n_samples, n_samples)) # Remove degeneracies S += ((np.finfo(np.double).eps * S + np.finfo(np.double).tiny * 100) * random_state.randn(n_samples, n_samples)) # Execute parallel affinity propagation updates e = np.zeros((n_samples, convergence_iter)) ind = np.arange(n_samples) for it in range(max_iter): # tmp = A + S; compute responsibilities np.add(A, S, tmp) # I = A+S 每行最大值所在的列下标 I = np.argmax(tmp, axis=1) #Y=A+S 每行最大值组成的向量 Y = tmp[ind, I] # np.max(A + S, axis=1) # 把tmp 中最大值赋值为无穷小 tmp[ind, I] = -np.inf # Y2 = A+S 每行中第二大的元素 Y2 = np.max(tmp, axis=1) # tmp = Rnew， subtract 计算矩阵相减 np.subtract(S, Y[:, None], tmp) # 计算最大值部分 tmp[ind, I] = S[ind, I] - Y2 # Damping tmp *= 1 - damping R *= damping R += tmp # tmp = Rp np.maximum(R, 0, tmp) # 对角线上元素 tmp.flat[::n_samples + 1] = R.flat[::n_samples + 1] # tmp = -A tmp -= np.sum(tmp, axis=0) dA = np.diag(tmp).copy() tmp.clip(0, np.inf, tmp) tmp.flat[::n_samples + 1] = dA # Damping tmp *= 1 - damping A *= damping A -= tmp # Check for convergence E = (np.diag(A) + np.diag(R)) > 0 e[:, it % convergence_iter] = E K = np.sum(E, axis=0) if it >= convergence_iter: se = np.sum(e, axis=1) unconverged = (np.sum((se == convergence_iter) + (se == 0)) != n_samples) if (not unconverged and (K > 0)) or (it == max_iter): if verbose: print("Converged after %d iterations." % it) break else: if verbose: print("Did not converge") I = np.where(np.diag(A + R) > 0)[0] K = I.size # Identify exemplars if K > 0: c = np.argmax(S[:, I], axis=1) c[I] = np.arange(K) # Identify clusters # Refine the final set of exemplars and clusters and return results for k in range(K): ii = np.where(c == k)[0] j = np.argmax(np.sum(S[ii[:, np.newaxis], ii], axis=0)) I[k] = ii[j] c = np.argmax(S[:, I], axis=1) c[I] = np.arange(K) labels = I[c] # Reduce labels to a sorted, gapless, list cluster_centers_indices = np.unique(labels) labels = np.searchsorted(cluster_centers_indices, labels) else: labels = np.empty((n_samples, 1)) cluster_centers_indices = None labels.fill(np.nan) if return_n_iter: return cluster_centers_indices, labels, it + 1 else: return cluster_centers_indices, labels ``` 参数： >preference ：即距离矩阵对角线上的元素 convergencd_iter: 如果连续迭代convergence_iter 次结果无变化就停止 max_iter ： 最大迭代次数 damping: 公式中的$$\lambda$$ 应该在0.5和1直接取值 return_n_iter ：是否返回迭代次数 copy, verbose, 是scikit-learn 常见参数不在这里解释了 1. 10，11行判断和计算preference 的值，如果为空会使用中值 2. 18行把preference 赋值给S的对角线元素，这里就有一个numpy 的技巧了，第一次见到绝对搞不清楚在干什么，所以一定要很好的理解numpy 的切片是看懂scikit-learn 的基础 3. 20至34行做了初始化，A.R 都初始为零矩阵 4. 34 到54是在迭代R ， >这部分的注释在上面的代码中，由于发现公式中有个不等式条件的限制，公式里面max里面的表达式是计算A+S 矩阵第i行的最大值，其实只有每行最大值的那个元素才会受到等式的限制。(请好好理解这段话，这是理解上面代码计算过程的关键所在) 5. 57 到69是对A的迭代 >这里也要中分理解公式，a 的公式里面最后一部分是在求sum(max(R第k行每个元素与0)) 在加上r(k,k),考虑到不等式的限制 ，考虑到不等式的限制就是第k行的所有元素和-R(i,k) 代码中的实现是计算-A ，逻辑是相同的。 6. 72行至86是在判断算法是否已经收敛 这个判断的思想真的很巧妙,把连续convergence_iter次计算结构存入e，换句话解释收敛就是连续convergence_iter次迭代过程中，是聚类中心的点没有变化，自然不是的也没有变化，再看下面这行 ``` (np.sum((se == convergence_iter) + (se == 0)) != n_samples) ``` 很容易理解了，如果有变化就是会小于n_samples ，不会有大于n_samples的情况出现。 后面剩下的就相对容易了，自己读吧。 ## 总结 AP 算法我看在实际应用中也不会有多少用途，因为对大数据量很难做，可解释性也很弱。