刚接触聚类那会儿最让我头疼的不是算法怎么写而是每次用K-means之前都得先回答一个问题你到底想分几类业务上经常跟我说“大概5到8类吧”可数据跑出来手肘图跟滑梯似的SSE拐点不明显K值怎么选都有点道理选完还得想办法解释为什么是6类而不是7类。后来翻论文看到近邻传播聚类算法行内常叫AP算法Affinity Propagation第一反应是“不用预设聚类数目天底下还有这种好事”。等真在Matlab里把代码写出来跑通发现它省掉的不仅仅是选K的步骤还顺手把怎么定义聚类中心、怎么处理噪声点这些事全换了个思路。这篇就当是一份AP算法的落地笔记把核心原理、Matlab实现和我在里面踩过的几个坑一起梳理清楚适合被K-means选K折磨过、或者刚接触聚类想换种思路的朋友。1. 选K选到头秃AP算法就是来改这个命的1.1 K-means的死穴聚类数目必须事先给定先说一个很现实的问题K-means本身是个相当简洁的算法几步就能写完选K个初始中心算距离分配样本更新中心位置重复到收敛。但“选K个初始中心”和“每个簇更新均值”这两步都隐含了一个前提——你必须先把K定下来。K是几出来的中心就是几个均值点K是5聚类结果就绝对不会出现第6个类别哪怕数据里真的藏着第6个深层的子群体。于是大家习惯用肘部法画出不同K值下的SSE找拐点用轮廓系数遍历K选平均轮廓最大的或者干脆业务拍板。肘部法的毛病是拐点经常不明显尤其是高维数据或者簇与簇有重叠时SSE曲线几乎单调下降所谓“拐点”纯属心理安慰。轮廓系数在数据紧致时好用但真实数据里噪声一多轮廓值普遍偏低选出来的K也未必符合业务直觉。就算K选对了K-means对初始中心非常敏感同一个K跑十次可能出十种差不离但不同的结果这导致你很难判断当前这个分法是真的结构还是初始化碰巧走运。说白了K-means把最难的问题扔给了使用者。1.2 AP算法把“定K”换成了“定尺度”AP算法的思路完全不一样它不要求你事先给聚类数目而是把所有样本点同时当作潜在的聚类中心。算法通过样本与样本之间的“消息传递”自己判断哪些点有资格当老大哪些点该跟着谁混。整个过程只有两个核心参数需要关注——preference和damping factor。你不需要回答“分几类”只需要回答“尺度多细”聚类数目会在迭代过程中自动浮现出来。这个特性在实际项目里很舒服。比如拿到一批用户行为数据你根本不知道分成几类合适业务也只能给个模糊描述。AP先帮你跑出一个“数据自己觉得合理的答案”你再结合业务去调preference的颗粒度整个过程有点像先让算法交一份参考答案再由人来批改而不是让人先猜答案再交给算法验证。工程上AP属于典型的O(N^2)复杂度算法适合几千个样本以内、希望聚类中心能用真实样本点解释的场景。它的聚类中心不是一个抽象的均值向量而是某个实际存在的数据点这在文本聚类、图像聚类、面波数据分群这类可解释性要求高的场景里特别讨喜。我用一张表快速对比一下K-means和AP在实际使用中的差异维度K-meansAP近邻传播聚类数目必须事先给定自动确定受preference影响聚类中心各簇均值向量虚拟点某个真实样本点初始化依赖对初始中心敏感基于消息传递对初始点不敏感可解释性中心是均值需要人工解读中心就是样本可直接查看原始记录适用规模大样本也很快几千样本以内比较舒服2. 消息传递干了什么吸引力、归属度和偏好值2.1 相似度矩阵AP眼里数据长什么样AP算法处理的第一步是把所有样本两两之间的“亲密程度”算出来形成相似度矩阵S。S(i,k)的含义是点k有多适合当点i的聚类中心。数值越大代表越合适数值越小代表越不合适。最经典的定义是负的欧氏距离平方S(i,k) -||xi - xk||²也就是说两个点越靠近这个值越大越接近0离得越远负得越厉害。如果数据是三维以上的特征向量直接算欧氏距离平方然后取负就行。如果你没有统计和机器学习工具箱pdist2用不了也可以用下面的向量化写法替代function S negativeSqDist(X) % X: n-by-d 数据矩阵 n size(X, 1); % ||xi - xj||^2 ||xi||^2 ||xj||^2 - 2*xi*xj sq sum(X.^2, 2); % 每个样本的模长平方 S 2 * (X * X) - sq - sq; S -S; % 取负变成相似度 end这里有个容易被忽略的细节矩阵对角线S(k,k)不能直接是0因为0意味着“点k自己当自己的中心”和“距离为0的点当中心”收益一样算法会分不清谁该当老大。AP的惯例是用preference偏好值填充对角元表示每个点“自荐当中心”的底气和意愿。这个值设得越大越多的点倾向于独立成簇设得越小大家就越愿意投奔别人。后面调参章节会细讲。2.2 R和A两个矩阵一场得体的互相试探AP算法的迭代过程核心是在两个矩阵之间反复交换信息吸引度矩阵RResponsibility常译作吸引度和归属度矩阵AAvailability常译作归属度。R(i,k)回答的问题是点i对点k说——“你够格当我的代表吗”它等于点i从k那里得到的相似度收益减去从其他所有候选中心里能拿到的最好收益除掉k自己。如果k是i看来最合适的中心R(i,k)就是正的如果点i觉得别人比k更合适R就变成负的。公式是R(i,k) S(i,k) - max_{k ≠ k} { A(i,k) S(i,k) }A(i,k)回答的则是另一个方向的表态点k对点i说——“大家都说你愿意支持我我收下你真的没问题吗”它等于点k自身的“自荐意愿”R(k,k)加上所有其他点对k的“支持票”总和再被限制到最高为0。也就是说哪怕大家一片热忱归属度也不会超过0一旦超过就说明k确实很适合当中心。公式是A(i,k) min(0, R(k,k) Σ_{i ≠ k} max(0, R(i,k)))对角元的归属度比较特殊A(k,k)是点k对自己作为中心的认可等于所有其他点对k的正向吸引度总和。你可以把两者的关系想象成一场选秀R是选手在大声喊“选我选我”A是观众在权衡“你那么出色我支持你”。一开始所有人的消息都是0谁也不服谁然后每一轮每个点根据收到的信息更新自己对其他人的看法。几轮之后有些点的净收益越来越高有些点慢慢认命选了别人最终整个网络达到一个稳定的“权力分配”。2.3 中心判定谁的对角元为正谁就是老大每轮迭代结束时算法会计算E R A代表每个点作为中心的综合实力。点i最终归属的中心就是那个让E(i,k)最大的k。如果某个点k满足“自己选自己”也就是E(k,k)大于0那它就被认定为聚类中心。实际代码里通常用下面两种方式之一来判定一种是用max(E, [], 2)找到每个点归属的中心索引然后检查哪些点的索引指向了自己另一种是直接看E对角元大于0的点。前者更直观后者更接近论文的数学表达两者在收敛后基本等价。迭代过程中还有一个必不可少的稳定器叫阻尼系数damping。R和A的更新如果不加限制消息可能来回震荡像两个人吵架越吵越激动永远停不下来。阻尼的数学形式是把上一轮的值和本轮的计算值按比例混合R_new (1 - damping) × R_calculated damping × R_oldA_new (1 - damping) × A_calculated damping × A_olddamping取值范围在0到1之间通常取0.5到0.9。值越大更新越保守越不容易震荡但收敛也越慢。3. Matlab手搓一个AP聚类器完整代码与跑通示范3.1 造一份能看出效果的数据集为了验证代码正确性我先造三簇二维高斯数据保证肉眼就能看出应该有三类。生成数据并计算相似度矩阵的代码如下%% 生成含三簇的模拟数据 rng(42); n1 80; n2 70; n3 90; X [randn(n1, 2) [2, 3]; randn(n2, 2) [-1, -2]; randn(n3, 2) [4, -1]]; % 相似度负欧氏距离平方 S -pdist2(X, X).^2; if isempty(which(pdist2)) S negativeSqDist(X); % 没有工具箱时的替代函数 end % 查看相似度矩阵量级后面设preference用 fprintf(相似度范围: %.2f ~ %.2f\n, min(S(:)), max(S(:)));跑完这段代码你会看到相似度矩阵的值大概落在-20到0之间。因为每簇内部点对之间距离小相似度更接近0而簇与簇之间距离大负得很厉害。这个量级信息在调preference时非常有用。3.2 AP主循环完整实现下面这个函数是我按论文标准流程实现的AP聚类不需要任何额外工具箱只用基础矩阵操作。整个逻辑包括填充对角元、初始化R和A、迭代更新、阻尼、收敛判断、输出结果。为了方便你直接抄作业我尽量把注释写全function [idx, centers, netsim] apcluster_demo(S, p, damping, maxits, convits) % 近邻传播聚类实现AP算法 % 输入 % S : n*n 相似度矩阵数值越大越适合互为中心 % p : preference标量或n维向量不传时取S中位数 % damping : 阻尼系数常用0.5~0.9 % maxits : 最大迭代轮数 % convits : 连续多少轮归属不变视为收敛 % 输出 % idx : n*1 簇标签从0开始编号 % centers : 聚类中心点的原始索引 % netsim : 收敛时所有样本到对应中心的相似度总和 n size(S, 1); if nargin 2 || isempty(p) p median(S(:)); % 默认偏置取中位数 end if nargin 3 || isempty(damping) damping 0.5; end if nargin 4 || isempty(maxits) maxits 1000; end if nargin 5 || isempty(convits) convits 100; end % 把preference写入对角元代表每个样本自荐当中心的意愿 if isscalar(p) S(1:n1:end) p; else S(1:n1:end) p(:); end R zeros(n, n); % 吸引度 Responsibility A zeros(n, n); % 归属度 Availability lastAssign zeros(n, 1); unchanged 0; for iter 1:maxits Rold R; Aold A; % ---- 更新吸引度 R ---- % R(i,k) S(i,k) - max_{k~k} { A(i,k) S(i,k) } for i 1:n item A(i,:) S(i,:); % 所有候选中心对i的综合收益 [bestVal, bestIdx] max(item); tmp item; tmp(bestIdx) -Inf; secondVal max(tmp); % 排除最优后剩下的最大收益 for k 1:n if k bestIdx R(i,k) S(i,k) - secondVal; else R(i,k) S(i,k) - bestVal; end end end % ---- 更新归属度 A ---- % A(i,k)min(0, R(k,k)sum_{i~k} max(0,R(i,k))) % A(k,k)sum_{i~k} max(0,R(i,k)) for k 1:n posSum sum(max(0, R(:,k))) - max(0, R(k,k)); A(k,k) posSum; for i 1:n if i ~ k A(i,k) min(0, R(k,k) posSum); end end end % ---- 阻尼 ---- if iter 1 R (1 - damping) * R damping * Rold; A (1 - damping) * A damping * Aold; end % ---- 收敛判断 ---- E R A; [~, assign] max(E, [], 2); if isequal(assign, lastAssign) unchanged unchanged 1; if unchanged convits break; end else unchanged 0; lastAssign assign; end end % ---- 整理输出 ---- E R A; [~, assign] max(E, [], 2); centers find(assign (1:n)); % 自己选自己的点就是聚类中心 [~, idx] ismember(assign, centers); idx idx - 1; % 转成0-based类别标签 netsim sum(S(sub2ind([n, n], (1:n), assign(:)))); end这段代码我实测下来跑三簇二维数据大概几十轮就收敛了速度很快。N在几千以内都还算能接受超过几千之后每次迭代的循环会明显变慢这点后面专门讲。3.3 跑通示例三簇数据聚类与可视化用我上面的测试数据调用函数看看结果%% 运行AP聚类 p median(S(:)); % preference 默认取中位数 [idx, centers, netsim] apcluster_demo(S, p, 0.7, 1000, 50); fprintf(自动聚类数: %d\n, numel(centers)); fprintf(聚类中心索引: %s\n, mat2str(centers(:))); fprintf(netsim: %.4f\n, netsim); %% 可视化 figure; gscatter(X(:,1), X(:,2), idx); hold on; plot(X(centers, 1), X(centers, 2), kp, MarkerSize, 14, LineWidth, 2); hold off; title(AP聚类结果黑色方块为自动选出的聚类中心);在我的环境里跑出来聚类数是3中心索引落在三个簇各自的密集区域。当你用真实数据跑时聚类数不一定会等于你内心的预期这是正常现象。AP是在优化“样本到中心的相似度总和”这个目标它没有“标签正确”的概念只有“数据里的强结构是什么”。所以别急着说算法错了先看看中心的样本在业务上意味着什么。4. 参数调优实战preference定尺度damping治振荡4.1 preference怎么选从默认值到“按需定制”preference是AP算法最关键的参数直接控制聚类数目的走向。默认值取S全矩阵的中位数理由是中位数对量级有很强的鲁棒性不会因为个别极端大或极端小的相似度拉动整体判断。你可以把中位数理解成“一个普通点自荐当中心的平均底气”比它高就更有底气比它低就更想抱团。调preference的经验规律很简单preference设置聚类趋势适用场景明显小于中位数甚至接近min聚类数减少、颗粒变粗只要大方向分块忽略局部细节取中位数附近默认适中第一次摸底、不确定时先用默认明显大于中位数甚至接近max聚类数增多、颗粒变细细粒度挖掘或存在大量局部子群有个很容易记混的点preference越大每个点越倾向于自立门户所以聚类数会变多preference越小大家更愿意投奔别人聚类数会变少。我第一次用的时候就反着调了怎么调都跟预期拧着来后来把对角元公式重新想了一遍才理顺。如果你想让某些特定样本更容易成为聚类中心可以把preference从一个标量改成向量。比如你事先知道某些点是可信的“种子”就把它们对应的p值设高一点其他点设低一点。这个技巧在异常检测和类别不均衡场景里很有用相当于给算法注入专家先验。4.2 damping、convits和迭代次数的权衡damping的取值直接决定算法是顺滑收敛还是来回震荡。经验上0.5能跑通大部分数据0.7到0.9是更稳妥的选择。如果运行过程中聚类中心数量一直在跳先把damping提到0.9试一下。升高阻尼的代价是收敛变慢迭代次数可能要翻倍但结果稳定比跑得快更重要。convits设多少也要看数据。设太小比如10轮算法可能在某个分配上来回微调时被误判为收敛结果还没稳定就停了设太大比如500轮又明显浪费时间。我的习惯是先用100轮做初筛看到结果再决定要不要加大。迭代次数maxits同理默认1000对付大部分场景够用。如果循环跑满maxits还没收敛函数不会报错但输出的聚类结果不可信入手点是检查damping和convits而不是硬调maxits。4.3 想要“恰好6类”用二分法搜索preference虽然AP卖点是不用预设聚类数目但很多时候业务会反问一句“你能不能给我分6类”这时候可以用二分搜索来反推preference。思路是preference变大聚类数变多变小聚类数变少那就在可行区间里不断二分找到让聚类数接近目标值的plo min(S(:)); hi max(S(:)); targetK 6; for t 1:30 mid (lo hi) / 2; [idx, ~] apcluster_demo(S, mid, 0.9, 1000, 100); k numel(unique(idx)); fprintf(p%.4f %d 类\n, mid, k); if k targetK break; elseif k targetK hi mid; % 类太多降低preference else lo mid; % 类太少提高preference end end注意AP的聚类数对preference不是严格单调的数据里有特殊结构时可能出现小幅反复。但这个二分法在实际项目里作为“粗调”非常好用跑几轮就能锁定数量级再手工微调就轻松了。5. 实际跑数据时踩过的坑与排查思路5.1 数据不标准化相似度矩阵直接带歪结果这个坑在我第一次跑真实多维数据时踩得很疼。AP用的是欧氏距离平方如果某个特征量纲范围很大比如年薪从3万到300万而年龄只有20到60距离平方几乎全被年薪主导年龄信息在相似度矩阵里约等于不存在。聚类结果出来每个簇的年龄分布完全是乱的业务根本看不下去。解决办法很直接进AP之前先把特征标准化X zscore(X); % 每列零均值单位方差如果你的数据是类别型特征或者数值型和类别型混着欧氏距离就不太合适了。可以考虑用Gower距离、余弦距离等替代方案再转换成相似度。核心原则是相似度矩阵必须真实反映你希望算法关注的结构而不是被某个量纲大的特征绑架。5.2 不收敛、振荡、中心数全变——三步排查链路如果运行结果不稳定先别急着换参数按照下面的顺序排查第一步检查相似度矩阵里有没有NaN或Inf。数据里有一个缺失值pdist2算出来的距离就可能出现NaNNaN在迭代里会传染最终聚类一团糟。运行前用any(isnan(S(:)))扫一遍有脏数据先清洗。第二步确认preference真的写进了对角线。很多人把S算完直接丢进AP函数忘了函数内部要设置S(1:n1:end) p。如果对角元一直是0算法会把“自己当中心”和“与距离为0的点同中心”等同起来中心判定会短路的。这也是为什么我在函数实现里特意把填充对角元放在迭代之前。第三步把damping拉高到0.9。振荡的典型表现是E对角元的符号在连续几轮里反复翻转也就是点一会儿自认是中心一会儿又投奔别人。高阻尼配合更长的convits能压住大部分振荡。如果0.9还不行试试0.95代价只是慢一点。5.3 数据量一大就内存爆炸AP的相似度矩阵是N乘N的稠密矩阵double类型每个元素占8字节。N2000时S矩阵大约32MBN5000就到200MBN10000直接800MB。这还没算R和A两个同样大小的矩阵三个加一起内存压力非常大。我在笔记本上跑5000个样本就已经能听到风扇狂转了。实际项目中如果数据超过几千条我常用的办法是先对数据进行一次下采样比如随机抽3000条跑AP得到聚类中心后再把全量数据按“离哪个中心最近”分配到对应簇。这样谱系结构来自AP大规模分配用简单的最近邻复杂度从O(N^2)降到O(N*K)。如果你非要在全量上做至少把矩阵数据类型改成single内存能省一半代价是相似度精度略降。另外一点如果数据本身有明确的相似度上限可以考虑稀疏矩阵。很多真实场景里距离较远的样本之间相似度非常低对最终聚类几乎没贡献。只保留每个样本最近的k个邻居的相似度把其他位置填0或-inf然后用稀疏矩阵重写迭代逻辑可以撑到更大规模。不过这就是AP走向工程化的话题了基础版本先把几千样本的坑摸清楚更实在。最后再分享一点个人体会我后来做聚类项目已经习惯先拿AP跑一轮再根据它的中心点往回看特征分布。这个过程不一定是最终交付的结果但它像一面镜子能照出数据里那些你没意识到的边界。它不会取代K-means更不会取代你对业务的理解但在“完全不知道有几类”的摸底阶段这种不需要预设聚类数目的算法确实帮我少熬了很多个选K的夜。如果你也被肘部图折磨过建议直接复制上面的代码跑一遍你自己的数据看看聚类中心落在哪些样本上再回头调preference你会对这份数据有新的理解。