AP聚类算法：消息传递机制与Python实战解析

简介AP聚类亲和传播是一种无需预设簇数量的无中心聚类算法这份资源适合想要理解并落地AP算法的Python开发者与数据分析学习者尤其适合在处理K-means难以确定K值的场景时参考。压缩包内共有1个文件为纯Python脚本体积约1KB代码精简便于直接阅读和复用。目前已有781人学习下载。脚本覆盖了AP算法的完整实现链路输入数据转换、欧氏距离矩阵构建、职责与可用性消息的迭代更新、收敛判断以及最终的簇中心分配和标签输出同时涉及numpy与scipy的典型用法可作为课程作业或项目原型的轻量参考。通过研读这份实现读者可以直观对照AP算法公式与代码步骤理解消息传递机制如何自动确定簇数并在此基础上根据自身数据调整相似度度量或停止条件。1. 从不确定簇数到自动发现为什么我改用AP聚类做数据分群时最常遇到的尴尬是客户画像还没看清先要回答到底分几类。K-means依赖肘部法则DBSCAN对密度差异敏感而当我们面对电商用户行为、基因表达谱这类结构未知的数据时预设类别数本身就是一个强假设。AP聚类Affinity Propagation的思路完全不同——它把每个样本点都当作潜在聚类中心通过样本间的消息传递自动收敛出簇的数量和成员关系。我最初接触apcluster.py是为了处理一份没有先验标签的用户行为数据当时用K-means反复调K值调到怀疑人生换成AP聚类后一次就跑出稳定的分群结果这种体验让我重新审视了这套2007年提出的算法——它并不新但在特征维度不高、样本量在千到万级别的场景下依然很能打。这篇博文会从消息传递的数学直觉讲起直接拆解apcluster.py的实现路径然后给出可直接运行的Python代码、参数调节经验和踩过的坑。无论是正在做聚类项目选型的初学者还是想评估AP聚类边界的老手这篇都能提供实际可用的参考。我和几个同事在实际项目中已经把这套代码用于用户分群和图像分割的预处理下面聊的每一条都是跑过真实数据后的结论。2. AP算法的消息传递机制责任度R与可用度A的迭代博弈2.1 为什么传统聚类在簇数不确定时失效K-means的目标函数是组内平方误差和WCSS优化过程依赖预定的K值层次聚类虽然不需要预设K但树状图的截断位置依然是人为主观决定的。DBSCAN用密度连通性定义簇但对eps和min_samples两个参数高度敏感遇到密度不均匀的数据几乎必然顾此失彼。这些方法的共同问题是簇结构假设先行数据拟合在后。AP算法换了个玩法——它把每个数据点都当作潜在的聚类中心exemplar通过点与点之间的信息传播让簇中心自己冒出来。这意味着簇的数量、簇的形状、中心的选取全部由数据本身决定不需要任何先验输入。核心创新在于两条消息的交替更新责任度responsibility和可用度availability这两个矩阵的博弈过程决定了最终的聚类结果。2.2 责任度R我有多适合代表你责任度矩阵R中的元素R(i,k)表示点k作为点i的聚类中心的合适程度。注意这个方向性它不是对称的。更新公式为# 责任度更新R(i,k) S(i,k) - max_{k ! k} [ A(i,k) S(i,k) ] def update_responsibility(S, A): S: 相似度矩阵 (n_samples, n_samples) A: 可用度矩阵 (n_samples, n_samples)初始为全零 返回更新后的责任度矩阵R n S.shape[0] R np.zeros_like(S) for i in range(n): # 对每个点i找到除当前k外最大的 A S 值 for k in range(n): AS A[i] S[i] AS[k] -np.inf # 排除当前k R[i, k] S[i, k] - np.max(AS) return R这段代码的逻辑是点i对点k的信任程度等于两者的原始相似度S(i,k)减去其他人对点k的信任总和的最大值。直观理解就是——虽然i觉得k不错但如果k已经在忙着代表别人那k对i的价值就要打个折扣。AS[k] -np.inf这个操作是为了在比较时排除当前候选中心k自身这是容易出错的细节很多人第一次写AP算法都在这里翻了车。2.3 可用度A别人是否认可你当中心可用度矩阵A中的元素A(i,k)表示点i选择k作为聚类中心的合适程度它综合了其他点对k的支持度。更新公式分两部分# 可用度更新 # 对角元A(k,k) sum_{i ! k} max(0, R(i,k)) # 非对角元A(i,k) min(0, R(k,k) sum_{i ! i,k} max(0, R(i,k))) def update_availability(R): R: 责任度矩阵 返回更新后的可用度矩阵A n R.shape[0] A np.zeros_like(R) # 先计算每个候选中心k的支持者总和 for k in range(n): # 正责任度贡献 positive_sum np.maximum(R[:, k], 0).sum() # 对角元除去k自己 A[k, k] positive_sum - max(R[k, k], 0) # 非对角元 for i in range(n): if i k: continue # 除去i和k之外的正责任度之和 others_sum positive_sum - max(R[i, k], 0) - max(R[k, k], 0) A[i, k] min(0, R[k, k] others_sum) return A可用度的含义很直观如果一堆点都认为k是优秀的中心即R(i,k)为正那么k作为中心的声望就高。非对角元的min(0, ...)限制确保可用度不会虚高——即使别人都支持kk对某个具体点i的可用度最高也只能到0相当于k可以带你玩但不欠你什么。2.4 迭代收敛两个矩阵的角力如何自然停止R和A的更新交替进行R更新时把A当作背景A更新时把R当作输入。这个循环一直持续到满足以下任一条件连续iter次迭代中簇中心不再变化达到预设的最大迭代次数R和A的变化量小于设定阈值def ap_cluster_step(S, R, A, damping0.5): 单步迭代交替更新R和A并应用阻尼系数防止震荡 R_new update_responsibility(S, A) A_new update_availability(R_new) # 阻尼当前值与上一步值的加权平均防止数值震荡 R (1 - damping) * R_new damping * R A (1 - damping) * A_new damping * A return R, A阻尼系数damping factor是AP算法收敛稳定性的关键默认取0.5实际使用中我通常调到0.9——因为更新过程中R和A的绝对值容易剧烈波动阻尼相当于给消息传递加了惯性。关于收敛判据我一般同时监控簇中心的变化和R矩阵对角元R(k,k)的取值后者直接决定哪些点最终成为exemplarR(k,k) 0意味着点k对自己有信心才能成为簇中心。3. apcluster.py的工程化实现从相似度矩阵到分群结果3.1 相似度矩阵的构建欧氏距离之外的选择AP算法的输入不是原始特征矩阵而是相似度矩阵S其中S(i,k)表示点i与点k之间的相似度数值越大表示两者越接近。最自然的做法是用负欧氏距离但实际项目中要分场景讨论。import numpy as np from scipy.spatial.distance import cdist def compute_similarity(X, metricneg_euclidean): 将特征矩阵转换为相似度矩阵 参数说明: X: (n_samples, n_features) 输入特征 metric: 相似度度量方式 if metric neg_euclidean: # 欧氏距离取负保证数值越大越相似 dist cdist(X, X, metriceuclidean) S -dist elif metric cosine: # 余弦相似度直接作为相似度 normed X / np.linalg.norm(X, axis1, keepdimsTrue) S np.dot(normed, normed.T) return S选型建议特征维度低且各向同性时用负欧氏距离文本向量或高维稀疏特征用余弦相似度。我在处理用户行为数据时发现归一化后使用欧氏距离与使用余弦相似度的分群结果差异很大——前者倾向于按行为强度的L2范数分后者纯按方向分业务上往往余弦相似度更合理。数据量超过5000条时S矩阵的存储就是一关5000×5000的float64矩阵占200MB这个成本要在建模前就心里有数。3.2 偏好值preference唯一需要人工设定的关键参数preference是S矩阵对角线上的取值表示每个点作为聚类中心的先验倾向。论文里的默认做法是取S非对角元素的中位数但这里有个反直觉的坑preference越大簇数越多preference越小簇数越少。中位数只是一个相对保守的起点实际业务场景中需要对preference做网格搜索。def fit_ap(S, preferenceNone, max_iter200, damping0.9, convergence_iter15): 完整AP聚类拟合流程 参数: S: 相似度矩阵 preference: 偏好值None时取中位数 max_iter: 最大迭代次数 damping: 阻尼系数 convergence_iter: 连续多少轮簇中心不变视为收敛 n S.shape[0] if preference is None: # 取非对角元的中位数作为默认偏好值 off_diag S[~np.eye(n, dtypebool)] preference np.median(off_diag) # 初始化 np.fill_diagonal(S, preference) R np.zeros_like(S) A np.zeros_like(S) # 迭代直到收敛 for step in range(max_iter): R_prev R.copy() R update_responsibility(S, A) A update_availability(R) R (1 - damping) * R damping * R_prev # 提取当前簇中心 exemplars np.where(np.diag(R) np.diag(A) 0)[0] if step convergence_iter: # 检查连续convergence_iter轮中心是否一致 break return exemplars, R, A关于preference的调参逻辑当业务目标偏向粗粒度分群时把preference调小让更少的点具备当中心的资格偏向细粒度时调大。这里注意一个边界——preference取S矩阵的最小值时所有点都会倾向于自己做中心最终每个样本点单独成簇取最大值时全部归为一类。fit_ap函数中的np.diag(R) np.diag(A) 0是提取簇中心的核心判据两个矩阵对角元之和大于0说明该点对自己成为中心的信心足够。3.3 簇分配一键输出labels当簇中心确定后每个非中心点分配到哪个簇取决于它和哪个中心的吸引度最高这个值用RA来衡量def assign_labels(S, exemplars, R, A): 为所有样本分配簇标签 返回: labels: 长度n的数组每个位置是对应的簇中心索引 n S.shape[0] labels np.zeros(n, dtypeint) # 每个点到每个簇中心的吸引力 R A for i in range(n): # 候选中心集合中找最大吸引力 scores R[i, exemplars] A[i, exemplars] best_idx np.argmax(scores) labels[i] exemplars[best_idx] return labels注意labels存的是簇中心的原始索引不是篡改成0到K-1的连续编号。这种保留原始索引的做法对后续分析更友好——簇中心本身是一个真实的样本点直接索引它对应的原始记录即可还原出代表样本的特征。实际分析时我通常把每个簇的中心点特征单独拉出来对照原始业务表做交叉验证。4. 实战用apcluster.py做客户分群并评估效果4.1 数据准备和AP聚类完整执行流程下面用一套模拟数据走完整流程。假设我们有1000条用户特征每条包含消费频次、客单价、活跃天数三个维度希望对用户自动分群不预设类别数。import numpy as np import pandas as pd from sklearn.preprocessing import StandardScaler # 生成模拟数据 np.random.seed(42) n 1000 features np.c_[np.random.gamma(2, 2, n), # 消费频次 np.random.uniform(20, 500, n), # 客单价 np.random.beta(2, 5, n) * 30] # 活跃天数 scaler StandardScaler() X_scaled scaler.fit_transform(features) # 计算相似度矩阵 S compute_similarity(X_scaled, metricneg_euclidean) # 进行AP聚类 exemplars, R, A fit_ap(S, preferenceNone, max_iter300, damping0.9) labels assign_labels(S, exemplars, R, A) # 统计结果 unique_exemplars np.unique(labels) print(f自动发现的簇数量: {len(unique_exemplars)}) for ex in unique_exemplars: cluster_size np.sum(labels ex) print(f簇中心: {ex}, 成员数量: {cluster_size})这里所有代码都遵循先标准化再计算相似度的流程否则量纲差异会主导距离计算分群结果毫无业务解释性。StandardScaler归一化后3个特征的量纲一致欧氏距离才有意义。4.2 评估AP聚类结果轮廓系数为什么只做参考AP聚类没有真实的标签做对照时评估主要依赖轮廓系数Silhouette Coefficient和簇内/簇间距离比对from sklearn.metrics import silhouette_score # 计算轮廓系数 sil silhouette_score(X_scaled, labels) print(f轮廓系数: {sil:.3f}) # 检查是否有异常簇成员数过少的簇需要关注 cluster_sizes np.bincount(labels, minlengthlen(np.unique(labels))) small_clusters [i for i, s in enumerate(cluster_sizes) if s 5] if small_clusters: print(f警告: 存在过小簇(成员5): {small_clusters})注意np.bincount要求labels是连续的0到K-1编号所以需要先用pd.factorize转换。轮廓系数在AP聚类中的参考价值有限——AP算法本身并不以欧氏距离的紧致性为目标函数它更看重消息传递的一致性所以轮廓系数低不等于结果差。我的经验是重点看簇规模分布是否合理以及簇中心是否有业务含义。簇中心是真实样本点直接看它们的原始特征值就能知道这簇人的典型画像这是K-means做不到的——K-means的中心是虚拟的平均点。预防过小簇的方法是在fit_ap后过滤掉成员数少于3的簇并把对应样本分给次近的簇中心。4.3 网格搜索preference的实操方法既然preference直接决定簇数量实际项目中我会在[p50, p75, p90]的分位数附近做网格搜索。注意术语对齐这里的分位数是S矩阵中非对角元素的分位数前文提到的中位数对应np.median。S值本身是负距离所以数值越大表示相似度越高preference取高分位点时会有更多点具备当中心的资格簇数随之增多。def grid_search_preference(S, p_vals[0.1, 0.3, 0.5, 0.7, 0.9]): 对不同preference做网格搜索返回簇数和评估指标 off_diag S[~np.eye(S.shape[0], dtypebool)] results [] for p in p_vals: pref np.quantile(off_diag, p) exemplars, _, _ fit_ap(S, preferencepref, max_iter300, damping0.9) n_clusters len(exemplars) results.append((p, pref, n_clusters)) return results网格搜索的价值在于建立参数-簇数-业务解释的映射关系。比如某次项目中preference取第50分位数得到5个簇簇间分离度好第70分位数得到9个簇出现两个成员极少的碎片簇。这时候宁可回到5个簇的配置也不要为了追求细粒度而牺牲稳定性。5. 收敛判据、性能优化和与sklearn实现的边界5.1 震荡问题的定位与阻尼调参策略AP算法在迭代后期最常见的异常表现是簇中心在两个状态之间来回跳这背后的原因是R与A的更新步长过大或者消息更新存在周期性振荡。关键现象是convergence_iter设得再大中心也无法稳定收敛。定位方法是在每次迭代时监控np.diag(R) np.diag(A) 0的参数变化def check_convergence_trend(R, A, history[]): 追踪簇中心变化趋势便于判断是否震荡 centers set(np.where(np.diag(R) np.diag(A) 0)[0]) history.append(centers) if len(history) 10: recent history[-10:] unique_center_sets len(set(map(tuple, map(sorted, recent)))) # 若最近10轮出现了超过3种不同的中心组合判定为震荡 return unique_center_sets 3 return False震荡的标准解法是调高阻尼系数到0.95。我在实践中还发现将preference的绝对值调大即让S矩阵的对角元素更极端有助于削弱震荡——极端偏好值会让消息更新更快分出胜负减少来回摆动的空间。如果调完阻尼依然震荡另一个思路是改用分块APFrey等人提出的策略不过工程上我一般建议先减少样本量到2000以内AP的收敛稳定性在规模增大时会显著下降。5.2 大规模数据的瓶颈相似度矩阵的内存占用AP的空间复杂度是O(n²)存储相似度矩阵加R、A两个矩阵总内存需求是3×n²×8字节。以10000个样本为例单精度也要2.4GB。工程上的常用解法有两个方向def mini_batch_ap(S_full, batch_size500, max_iter100): 小批量AP聚类策略 先用少量样本确定簇中心再分配其余样本 n S_full.shape[0] # 取出一个子集进行AP计算 idx_sample np.random.choice(n, batch_size, replaceFalse) S_sample S_full[np.ix_(idx_sample, idx_sample)] exemplars, R, A fit_ap(S_sample, max_itermax_iter) sample_labels assign_labels(S_sample, exemplars, R, A) # 对剩余样本计算到各簇中心的相似度分配到最近的簇 remaining np.setdiff1d(np.arange(n), idx_sample) center_indices idx_sample[exemplars] S_rem S_full[np.ix_(remaining, center_indices)] rem_labels np.argmax(S_rem, axis1) return remaining, center_indices, rem_labels小批量的代价是可能丢掉小簇。另外还可以用scipy的稀疏矩阵存储稀疏相似度比如k近邻截断但AP算法的消息更新在非全连接图上效果会打折扣——截断后一些潜在的簇中心关系被切断分群质量通常不如全集。所以我会根据数据规模分档处理5000样本以内跑全量5000到20000用小批量或多轮抽样取稳定中心超过20000直接换Mini-Batch K-means或谱聚类。5.3 和sklearn AffinityPropagation的实际差异前面全部是自实现的版本和sklearn.cluster.AffinityPropagation的差距要注意sklearn的实现在版本之间有多次更新且convergence_iter用的是簇集合的稳定次数而不是联合概率变化所以相同参数和自实现的结果可能有差异。需要明确的一点是版本差异影响API参数名和默认行为但不影响上述调参方向。在1.2.x以上版本中sklearn默认的AffinityPropagation对damping默认依然是0.5不建议用默认值调到0.9更稳。from sklearn.cluster import AffinityPropagation # sklearn接口 af AffinityPropagation( damping0.9, max_iter300, convergence_iter15, affinityeuclidean # 内置计算负欧氏距离 ) labels_sk af.fit_predict(X_scaled)自实现版本的优势在于可以直接打印每一轮的R和A矩阵监控数值走向方便教学更容易嵌入自定义相似度计算逻辑比如带业务权重的相似度。sklearn版本胜在C加速用稀疏图矩阵时效率高不少而且少写很多底层代码。如果是项目落地我建议用sklearn做生产版本把apcluster.py保留为教学和理解用途——这也是这个文件最有价值的场景。5.4 一个容易误用的细节preference的取值方向再强调一次preference与簇数的关系这是AP算法最容易误用的参数方向。S矩阵中对角元的初始值代表每个点认为自己适合当中心的程度这个值越大竞争中心的门槛越低所以簇数多越小只有相似度最高的少数点能冒头所以簇数少。有人用负无穷做preference发现所有点都成了中心这是因为每个点和自己相似度最高在无竞争者的情况下全被选上——这在数据有重复或近重复点时是个真实陷阱用参数网格搜索比凭直觉设定靠谱得多。最后提供一条性能优化的实用建议。当发现自实现函数在大数据上迭代很慢时优先检查两段代码一是update_responsibility中的双循环是否可以用向量化替代二是每次迭代拷贝R矩阵的开销是否值得用原地更新节省def update_responsibility_vectorized(S, A): 向量化责任度更新避免显式双循环 n S.shape[0] AS A S # 先找到全局最大值再处理对角线排除 max_indices np.argmax(AS, axis1) max_vals AS[np.arange(n), max_indices] # 对角线元素单独处理 max_vals max_vals 1e-10 # 避免除零 R S - max_vals[:, None] # 恢复对角线公式中排除k自身后的最大值 for i in range(n): AS_i AS[i].copy() AS_i[i] -np.inf R[i, i] S[i, i] - np.max(AS_i) return R向量化版本大约有2-3倍加速如果数据超过2000样本建议直接走sklearn的C加速实现把自实现版本当作理解算法的脚手架。整体来看自实现版本在学术验证和教研场景价值高生产环境做分群时我更倾向于用sklearn配合自研的相似度矩阵传入affinityprecomputed来兼顾两头的优势。本文还有配套的精品资源点击获取