AP近邻传播聚类算法：原理、Matlab实现与调参实战

近邻传播聚类算法AP算法Affinity Propagation是我这几年做无监督聚类时最常用的一招。它最大的好处就写在标题里不需要预先确定聚类数目也不需要手动指定聚类中心跑完一轮迭代点和类全给你分好了。之前带过一个做基因表达谱聚类的学生数据没有标签导师只丢一句你自己看着分他拿K-means画肘部图画了半天也很难说清K到底取3还是6换AP算法十几行代码跑完聚类中心和簇数自己冒出来反而省了最多的时间。这篇文章我就把AP算法的原理、Matlab实现、调参手感、踩过的坑一次讲清楚给想在Matlab里直接跑通的人一份能抄作业的完整代码。1. 不用预设聚类数之后K-means留下的坑谁来填1.1 K-means的K值困境为什么这么折磨人传统聚类里K-means的问题是出了名的K值不给定就根本没法迭代给错了又全盘皆错。很多人依赖肘部法则画SSE曲线找拐点但真实数据的拐点往往不锐利误差平方和从K2到K10一路平滑下降你说拐在哪不是搞数据的人可能觉得这只是参数选取的小问题但业务上K值往往就是最核心的先验知识——你不知道用户群体到底分几类不知道图像里有多少个物体不知道人脸库里有多少个身份。K-means还有第二个隐含痛点它算出来的聚类中心是虚拟均值点不是真实样本。你聚类酒店客户中心点是个谁都不像的合成客户聚类行人轨迹中心轨迹也不存在于任何真实轨迹库中。很多业务场景其实希望聚类中心是真实存在的样本这样可以直接拿去当代表、当模板、当锚点。1.2 AP算法是怎么绕开这两个坑的AP算法全称Affinity Propagation2007年发表在Science上思路完全换了个方向它不预设中心数量而是让每个数据点当选民在迭代中互相发送消息评价你适不适合当我的聚类中心最终自己涌现出一组真实样本作为聚类中心论文里叫exemplar也叫代表点。这相当于把我要分几类这个宏观问题转化成了每个点该选谁当代表的微观投票问题宏观结果由微观互动自然涌现。这是AP和K-means最本质的区别前者先假设簇数再求解簇后者先求解簇再统计簇数。1.3 AP算法必须先认识的三个核心概念要跑通AP必须先分清三个东西相似度矩阵SS(i,j)表示点j适合作为点i聚类中心的程度。通常取负欧氏距离平方值越大表示j越适合当i的代表。参考度preference矩阵对角线S(k,k)代表点k对自己当中心的意愿。preference越大越容易出现更多聚类中心反之聚类中心变少。消息传递的两个量责任度R和可用度A下文详细讲。这三者决定AP的全过程。你不用理解太深也能跑代码但如果看不懂R和A的更新公式遇到不收敛、中心数不合理时就会一头雾水。下面用最直白的方式拆解消息传递机制。2. 责任度与可用度AP算法内部的拉票和投票2.1 责任度R(i,k)点i告诉点k我看好你责任度R(i,k)描述的是数据点i在考虑了其他候选中心之后有多大的把握觉得点k适合当自己的中心。注意这里面有个竞争逻辑——点i不是单纯看k对自己好不好而是看k对我的吸引力减去除k之外的所有候选中心里对我吸引力最强的那一个R(i,k) S(i,k) - max_{k≠k} [A(i,k) S(i,k,i等效)]这个式子的直译是你k作为我的中心比我可能选的其他任何中心强多少如果其他候选中心的总吸引力都比你高R就会偏负只有在k确实是局部最优时R才可能为正。你可以把它理解成点i向点k喊话在所有候选人里我目前最挺你。2.2 可用度A(i,k)点k反过来跟点i说你可以站过来光有人喊人选不够候选中心还要获得群体的背书这就是可用度A(i,k)当i≠k时A(i,k) min{0, R(k,k) Σ_{i≠i,k} max(0, R(i,k))}当ik时A(k,k) Σ_{i≠k} max(0, R(i,k))可用度衡量的是除了点i自己之外还有多少其他点也支持k当中心。如果一大片点都对k表达了正的责任度k这时候就有底气站出来说你们都过来吧。公式里有个R(k,k)这是k的自评——它认为自己适不适合当中心后面那一大串是别人对k的支持累计。两者相加足够高k才配得到正的可用度。直白说就是我不但觉得自己行而且还有很多人愿意选我那我才是真的中心。2.3 阻尼因子到底在干什么AP迭代很容易在两个状态之间反复横跳就像一群人争论谁当组长今天推张三、明天推李四。为了压制这种震荡算法引入阻尼因子damping常规做法是把新旧值加权混合R_new (1 - damping) * R_raw damping * R_oldA_new (1 - damping) * A_raw damping * A_olddamping取值在0到1之间一般取0.9震荡严重时可以取0.95甚至0.99。它本质上是给R和A更新加了惯性让消息变化更平滑。很多人一开始跑AP发现不收敛八成就是阻尼因子设得太低或者干脆没做阻尼更新。2.4 聚类中心揭晓的时刻RA才是最终裁判迭代收敛后把责任度和可用度相加E R A。每个数据点i最终归属到使E(i,k)最大的那个点k而真正的聚类中心是那些满足E(k,k) 0的点。这个判据很有意思它是自评加上他评。一个点只有当自己愿意当中心同时别人也认可自己当中心时才会成为最终的中心。所以AP选出来的中心不是虚构均值而是真实存在的样本点这对很多业务场景非常宝贵。3. Matlab代码实现从数据生成到AP主循环3.1 演示数据不要只用两个圆球很多教程喜欢用两个高斯分布当聚类数据两个圆球谁都能分。这里我构造三个有区别的簇两个高斯簇加一个拉长的条状簇。拉长的条状簇对K-means不太友好但AP基于样本间的相似度只要相似度矩阵设计合理往往能分出更有解释力的簇。数据规模选520个点方便观察聚类中心和计算速度。3.2 相似度矩阵的构造方式相似度矩阵我直接取负欧氏距离平方S(i,j) -||x_i - x_j||^2距离越小相似度越大负值越接近0。在Matlab里一行搞定S -pdist2(X, X).^2;注意这里需要Statistics and Machine Learning Toolbox。如果没有这个工具箱可以用下面这段替代n size(X,1); D zeros(n); for i 1:n D(i,:) sum((X - X(i,:)).^2, 2); end S -D;3.3 AP主函数完整可运行的Matlab代码下面这段是我反复用的一套实现结构清晰适合学习和改造。代码里对角元preference有三种模式取中位数、取最小值、直接传数值。function [labels, exemplars, R, A] apcluster(S, pref_mode, damping, maxit, convit, tol) % apcluster 近邻传播聚类 % 输入: % S : n x n 相似度矩阵S(i,j)表示点j适合作为点i中心的程度 % pref_mode: median或min或数值参考度设置 % damping : 阻尼因子默认0.9 % maxit : 最大迭代次数默认1000 % convit : 连续稳定多少次后停止默认100 % tol : 收敛容差默认1e-5 % 输出: % labels : n x 1每个点归属的聚类中心索引 % exemplars : 聚类中心索引列表 % R, A : 责任度与可用度矩阵 if nargin 6, tol 1e-5; end if nargin 5, convit 100; end if nargin 4, maxit 1000; end if nargin 3, damping 0.9; end if nargin 2, pref_mode median; end n size(S, 1); % 设置参考度对角元 if ischar(pref_mode) switch pref_mode case median p median(S(:)); case min p min(S(:)); case max p max(S(:)); otherwise error(未知的pref_mode); end else p pref_mode; end S(1:n1:end) p; R zeros(n); % 责任度矩阵 A zeros(n); % 可用度矩阵 stable 0; for iter 1:maxit % ---------- 更新责任度 R ---------- AS A S; % 候选中心吸引力 支持度 相似度 [max1, idxmax] max(AS, [], 2); % 每行最大吸引力 max2 zeros(n, 1); for i 1:n row AS(i, :); row(idxmax(i)) -Inf; max2(i) max(row); % 每行次大吸引力 end Rnew S - max1; % 非最大位置用最大值 lin sub2ind([n n], (1:n), idxmax); Rnew(lin) S(lin) - max2; % 最大位置用次大值 % ---------- 更新可用度 A ---------- Rp max(Rnew, 0); % 只保留正责任度 colsum sum(Rp, 1); % 每个候选中心获得的总支持 Rdiag diag(Rnew); Anew zeros(n); for k 1:n col Rdiag(k) colsum(k) - Rp(:, k) - Rp(k, k); col min(col, 0); % 非对角项有上界0 col(k) colsum(k) - Rp(k, k); % 对角项直接取总支持 Anew(:, k) col; end % ---------- 阻尼更新 ---------- Rnew (1 - damping) * Rnew damping * R; Anew (1 - damping) * Anew damping * A; % ---------- 收敛检测 ---------- diffR max(abs(Rnew(:) - R(:))); diffA max(abs(Anew(:) - A(:))); R Rnew; A Anew; if max(diffR, diffA) tol stable stable 1; if stable convit break; end else stable 0; end end if iter maxit warning(达到最大迭代次数可能未完全收敛建议增大maxit或调整damping); end % ---------- 提取聚类中心与标签 ---------- E R A; [~, labels] max(E, [], 2); exemplarFlag diag(E) 0; exemplars find(exemplarFlag); end3.4 这段代码为什么这样写先说责任度更新那一步。很多人第一次写AP会用三重循环对每个i、每个k都算一次去掉k之后的最大值N稍微大一点就慢得没法看。这里用了一个经典技巧每行先找最大值max1和它的位置idxmax再找次大值max2。对非最大值位置竞争者是max1对最大值位置竞争者变成max2。这样一次就解决一整行的R更新把O(N^3)降到O(N^2)。这个技巧是AP代码高性能运行的关键我建议直接背下来。可用度更新里那段循环看似绕其实逻辑是这样每个候选中心k先统计全列正责任度之和colsum(k)然后减去自己的责任度贡献再减去点i自己的贡献剩下的就是别人对k的支持。加上k的自评R(k,k)得到候选中心的整体可用性同时通过min(0,.)把非对角项压成有界值避免数值无限增长。对角项单独处理因为k自己的可用度只取决于别人对它的支持不应该被自己限制在0以下。3.5 调用示例与可视化下面这段demo脚本可以直接复制运行看看AP在你的机器上落地效果% demo_ap.m rng(42); % 簇1高斯中心(2,2) X1 randn(200,2) * 0.6 [2 2]; % 簇2高斯中心(-2,-2) X2 randn(200,2) * 0.6 [-2 -2]; % 簇3拉长条状 t linspace(-3, 3, 120); X3 [t*0.4 0.2, t*0.06 1.2] randn(120,2)*0.05; X [X1; X2; X3]; % 相似度矩阵 S -pdist2(X, X).^2; % 跑AP [labels, exemplars] apcluster(S, median, 0.9, 2000, 50); % 可视化 figure; gscatter(X(:,1), X(:,2), labels, [], .); hold on; plot(X(exemplars,1), X(exemplars,2), ro, MarkerSize, 14, LineWidth, 3); hold off; title(AP聚类结果红色圆点为自动发现的中心); xlabel(x1); ylabel(x2); disp([自动确定的聚类中心个数: num2str(length(exemplars))]);我这边用Matlab R2021a实测preference取中位数时这个数据会稳定收敛到3个簇三个中心分别落在三个簇内部拉长条状簇也能被单独识别出来不会像K-means那样经常把它拦腰切两半。4. preference和damping两个旋钮的调参手感4.1 preference取中位数是最合理的起点preference设多少直接决定聚类数。我推荐先用median(S(:))也就是所有相似度值的中位数。这个取值通常会让AP呈现数据中最自然的簇结构不多不少。取太小会倾向于合并成一个大簇取太大会倾向于把每个点都当成孤立中心。实际操作里我会先跑一遍中位数看聚类数如果想要更多细节就往大调如果想要更少的粗粒度簇就往小调。这个旋钮比K-means的K值直观得多因为它是连续量变化是平滑的。4.2 damping的取值决定了收敛稳定性damping低于0.7时我的经验是很容易出现震荡。有一次跑一个1000点的数据damping0.8连续100次迭代都不收敛R和A一直在两个状态之间跳。后来把damping加到0.95收敛得很稳定。这不是说你永远都要0.95而是如果发现有震荡倾向优先去调damping而不是去动preference。对绝大多数数据集0.9到0.95之间是甜点区。4.3 用一个扫描实验理解preference和聚类数的关系可以写个循环从相似度矩阵的第10百分位扫到第90百分位观察聚类中心数目的变化pvals prctile(S(:), 10:10:90); kcounts zeros(size(pvals)); for i 1:length(pvals) [~, ex] apcluster(S, pvals(i), 0.9, 1000, 50); kcounts(i) length(ex); end plot(10:10:90, kcounts, -o); xlabel(preference百分位); ylabel(聚类中心个数);跑完这张图你会发现一个规律preference从小到大聚类数大体是从少到多单调递增的。这给了你一个清晰的操作思路想要几类就去对应的preference区间里找。K-means要反复试探K值AP则是在连续区间里做选择每个数都有对应的连续位置。4.4 配合轮廓系数自动选preference如果不想人工盯着选可以结合轮廓系数做自动搜索。对扫描出的每个preference调用Matlab的silhouette函数计算平均轮廓值选最大的一组bestP pvals(1); bestScore -inf; for i 1:length(pvals) [labels, ~] apcluster(S, pvals(i), 0.9, 1000, 50); if length(unique(labels)) 1 s silhouette(X, labels); score mean(s); if score bestScore bestScore score; bestP pvals(i); end end end注意silhouette函数的性能在大样本下会比较慢建议N不超过3000时使用。轮廓值衡量的是簇内紧致和簇间分离的平衡用它选preference比纯靠肉眼看图靠谱。5. 我在实际项目中踩过的坑5.1 相似度矩阵的方向性搞反聚类结果直接乱套AP的相似度矩阵S(i,j)表示点j适合作为点i中心的程度这个方向很多人上手时会犯迷糊。如果你把S定义成S(i,j) -||x_i - x_j||^2那它就是对称矩阵方向上怎么理解都不会有影响因为距离是天然对称的。我建议统一用对称距离型相似度这样即使方向理解错了结果也不会出错。除非你用的是不对称相似度比如有向图上的传播强度那时候方向就极其重要更新公式里的i和k不能调换。我在一个网页推荐相似度项目里用过非对称相似度一开始R更新写反了结果聚类中心全落到边上排查了很久才发现问题是谁适合当谁的中心方向搞反了。5.2 内存爆炸N5000时S矩阵有多占地方AP的空间复杂度是O(N^2)这一点必须提前心里有数。double类型的相似度矩阵样本数NS矩阵内存占用10008 MB200032 MB5000200 MB10000800 MB200003.2 GBN到10000时光S矩阵就占800MB再加R和A单机直接吃满。超过5000个点我基本不会直接跑原始AP而是先采样一部分算中心再用最近邻把剩余点映射到簇上。如果业务上要求全量样本考虑用Fuzzy AP或者分块近邻传播。还有一个小技巧相似度矩阵用single类型内存直接减半在距离精度要求不高时影响可以忽略。5.3 迭代不收敛震荡、重复样本和preference取值震荡有几个常见来源。第一是damping太低解决方法是调高到0.95以上。第二是数据里有大量重复样本或者几乎相同的样本它们之间相似度太高导致R和A在几个点之间反复拉扯可以先做去重或者给坐标加微小噪声。第三是preference取到极端值系统容易处于临界状态比如某个点左右摇摆一会儿当中心一会儿不当中心。这种情况可以把convit设大一点比如连续稳定300次再停同时把tol放宽到1e-4反而更容易收敛。5.4 高维特征一定要先归一化再算相似度AP对相似度尺度非常敏感。如果你直接拿原始特征的欧氏距离平方当相似度量纲大的特征会完全主导相似度矩阵。比如特征图里一个维度是灰度值0到255另一个维度是纹理占比0到1后者基本被淹没。正确做法是先做标准化比如z-score或者min-max归一化再算相似度。当然如果你用余弦相似度做语义向量聚类一般不需要归一化因为余弦本身已经刻画方向关系这在文本句向量聚类里面尤其好用。5.5 大数据量怎么办先采样再映射当N很大时我用的方案是先对数据做随机采样得到几千个点跑AP得到中心然后对每个中心用最近邻规则把剩余样本映射到最近的聚类中心。这个策略在几万个点的项目上也跑得动。如果连采样后的几千个点都嫌慢还可以先做PCA降维到低维空间再用AP速度会明显提升聚类结果在高维数据上通常也不会差太多。6. 和K-means、DBSCAN、层次聚类放在一起比6.1 四类聚类方法核心差异方法是否自动定K中心是什么主要参数复杂度稳定性K-means否需指定K虚拟均值点K、初始中心O(NKI)依赖初值多次运行结果波动AP是自动涌现真实样本点preference、dampingO(N^2T)稳定结果可复现DBSCAN是基于密度连通无显式中心eps、minPtsO(N^2)依赖邻域参数对密度变化敏感层次聚类否需切树状图各层簇的代表链接准则O(N^2)~O(N^3)稳定可复现这个表一眼就能看出AP的独特定位。K-means胜在极限速度但K是个大麻烦DBSCAN不需要K但eps这个参数对密度变化非常敏感一组参数在稠密区域效果好到稀疏区域就全被标成离群点层次聚类稳定但算距离矩阵加不断合并大样本下速度吃亏。AP夹在中间不需要K中心是真实样本结果稳定还能通过preference连续调节粗细粒度。6.2 什么场景我推荐AP如果业务上满足下面任意一条我首先想到AP聚类中心必须是真实样本比如推荐系统里要点名这个用户是这群人的代表人脸聚类里要用真实人脸图当样板。不知道分几类且没有可靠的先验标签参考。比如新业务用户行为探索老板说你先看看分几拨人。需要结果可复现不同人跑出来一样。K-means换个种子结果就变AP在相同参数下是确定性的。数据规模在几千这个量级再大还能用采样方案兜底。文本句向量聚类、图像区域分割、人脸聚类、社交网络群体发现这些场景AP都有很成熟的应用案例。6.3 什么场景我劝你别用AP没有银弹AP也有明显边界样本超过几万内存和耗时都扛不住优先mini-batch K-means或BIRCH这类流式方案。数据形状非常复杂比如螺旋形、环形嵌套AP基于点到点相似度的思路很难正确聚出这些结构这时候DBSCAN或谱聚类更合适。数据维度极高又对欧氏距离不敏感切到余弦相似度再看AP偶尔表现不稳定需要先降维。所以我的经验是选聚类算法先回答两个问题中心需不需要是真实样本以及数据规模有多大。答案指向真实样本、中等规模时AP几乎是最省心的一档。我自己现在遇到类似你看着分的无标签数据第一反应就是先跑一次AP拿基线结果再根据业务需求做微调。跑得通就用跑不通再转K-means或者DBSCAN这样做基本没有翻车的案例。文章里的代码我放在一起就是一套可直接复用的工具箱下一次遇到聚类问题直接改相似度矩阵和preference很快就能看到结果。