先说实话这个标题看起来就是一个课程作业级别的项目但真正动手写过之后我发现AGNES这种最直观的聚类算法想在C里写得既正确又高效坑远比想象中多。网上能找到的示例大部分是Python调包稍微改改就能出结果一旦落到C距离矩阵怎么存、合并之后怎么更新、标签怎么维护每一处都是取舍。这篇文章把我自己从零实现AGNESAGglomerative NESting凝聚式层次聚类的完整过程、踩坑记录和优化思路整理出来代码可以直接编译运行适合要在生产系统里嵌入聚类能力、或者正在做相关课程设计的朋友参考。1. 先说清楚AGNES到底在干什么适用场景与算法全景1.1 从需求场景看层次聚类的存在价值很多人问过我一个问题聚类用k-means不就行了为什么还要研究层次聚类这个问题的答案取决于你要的是一组互不重叠的类别标签还是一份带有从细到粗演进关系的类别层级。举几个我实际遇到过的场景。给电商用户做分层运营时运营同学不只是想知道用户分成几类他们更希望看到哪两类人群可以先合并合并之后和再上一级的群体差异在哪这样才能设计梯度化的运营策略。再比如做基因表达谱分析样本之间的层级关系本身就有生物学意义相近的样本在进化树上更靠近这种关系用k-means直接输出扁平标签是体现不了的。还有一类场景是k未知且分布不均衡比如图像颜色量化、文档主题分组你没法一句话说出到底分几类最合适这时候一个完整的合并树可以让你在不同粒度上反复观察。AGNES作为层次聚类里最常用的自底向上算法思路朴素得有点不像机器学习算法一开始每个样本算一个簇然后反复找距离最近的两个簇并合并直到簇的数量达到你预设的目标或者合并距离超过某个阈值。整个过程记录每一次合并的先后顺序和距离就形成了一棵层次树俗称树状图dendrogram。这棵树本身就是信息你可以从中发现数据的自然结构、异常值、乃至噪声的分布状况。1.2 AGNES的执行流程拆解把AGNES的执行逻辑写成一个循环其实只有四步初始化把每个样本点单独看作一个簇。计算/查找根据当前所有簇两两之间的距离找到最近的那一对。合并把这一对簇合并成一个新簇。更新计算这个新簇到其余簇的距离替换掉原来的两行距离信息。重复步骤2到4直到剩余簇数等于目标k。步骤1和步骤2都有标准做法步骤3和步骤4是真正的实现重心也是我要花最多篇幅讲的地方。1.3 与k-means、DBSCAN的核心差异我整理了一张表方便你一眼看清楚AGNES在算法家族里的位置维度AGNES层次聚类k-meansDBSCAN是否需要预设k需要或者用树状图事后选k必须预设k不需要但要调邻域半径和密度阈值输出形式聚类标签 完整合并树扁平标签扁平标签 噪声点标记对簇形状的假设依赖链接规则单链接可处理长条簇偏向球形簇可处理任意形状适合密度不均稳定性每次结果确定若距离度量固定受初始中心影响受参数影响时间复杂度O(n^3)朴素实现O(nkiter)有索引时接近O(n log n)这张表不是说要你抛弃k-means而是想说明AGNES适合的是那些结果分布的形状没有先验假设、且你需要层级关系来辅助决策的场景。它的缺点也很明显——朴素实现是O(n^3)的样本量上到几千以后就必须考虑优化手段。这一点我会在第5章展开讲。2. 写实现之前必须确定的三件事链接规则、更新公式、标签方案2.1 簇间距离的三种链接规则算法第二步找到最近的两个簇这里的难点在于簇不是单个点两个簇之间距离最近并没有唯一标准。常见的有三种链接规则我在实现里把它们全部做成了可配置项单链接single linkage两个簇之间的最短距离等于各自包含的点之间所有距离的最小值。全链接complete linkage两个簇之间的最长距离等于各自包含的点之间所有距离的最大值。平均链接average linkage两个簇之间所有点对距离的平均值。它们的差异用一句话概括就是单链接只看最接近的那条边因此容易被离群点带偏但能发现细长的条形簇全链接只看最远的那条边倾向于生成紧凑、球形的簇对噪声更敏感平均链接在两者之间取平衡实践中最常用。这里有一个新手容易绕晕的地方单链接并不等价于当前距离矩阵中取最小值。距离矩阵每一行对应一个簇到其他簇的距离这个值本身就已经是簇间关系比如单链接下它记录的是min距离合并时需要把两个旧簇到第三个簇的min、max或avg算出来更新到矩阵里而不是重新去遍历原始样本点。2.2 平均链接的更新公式UPGMA与简单平均的差异如果你要在不同教程之间对比代码会发现平均链接有两种写法简单平均newDist (d[clusterA][C] d[clusterB][C]) / 2加权平均UPGMAunweighted pair group method with arithmetic meannewDist (sizeA * d[clusterA][C] sizeB * d[clusterB][C]) / (sizeA sizeB)很多演示代码为了让写法简单直接用简单平均。但严格来说AGNES里标准意义的平均链接是UPGMA它按两个簇包含的样本数量做加权。原因很直观假设簇A有10个点簇B只有1个点簇C有5个点。A到C的平均距离是2B到C的距离是10。如果不加权合并A、B之后到C的距离会算成6实际上合并后的簇里10个点和C挨得很近只有1个点离得远整体距离应该更接近2才对。UPGMA算出来是(10*21*10)/11≈2.73明显更合理。这个差异会影响合并次序进而影响整棵树的形状。我在第4章的踩坑部分会用一个具体例子说明这个问题。2.3 标签维护并查集的简化版还是直接改标签合并过程中最容易被忽略的是簇的标识怎么管理。系统里每个样本的簇归属是不断在变的一个常见做法是维护一个parent数组合并时把被合并簇的parent指向保留簇最终通过向上找根的方式获得簇标签。这其实就是并查集Union-Find的思想只是AGNES里面没有查找按秩合并那种复杂操作因为合并顺序完全由距离决定我们只是需要一个数组记录每个点的最终归属。另一个简单方案是每次都给所有点的label做一次重编号但这样每轮合并都要全量扫描没必要。我采用的方式是parent[i]初始等于i表示每个点自成一簇。合并簇j到簇i时只执行parent[j] i。打完标签后对每个点执行一次findRoot逻辑把链式关系收拢到最终簇根。这种方案天然能处理簇的簇这种层级结构而且合并历史也不丢。下面讲核心实现时你会看到它和距离矩阵更新配合得非常顺畅。3. 核心代码拆解合并循环里的每个关键分支3.1 数据结构和距离矩阵初始化整个类的骨架如下#include cmath #include limits #include vector #include iostream #include iomanip #include algorithm struct MergeRecord { int left; int right; double distance; int mergedSize; }; class AgnesClusterer { public: enum Linkage { SINGLE 0, COMPLETE 1, AVERAGE 2 // UPGMA 加权平均 }; AgnesClusterer(const std::vectorstd::vectordouble points, int targetK, int linkageType AVERAGE) : data(points), k(targetK), linkage(linkageType) { } std::vectorint run(); const std::vectorMergeRecord history() const { return mergeRecords; } private: std::vectorstd::vectordouble data; int k; int linkage; std::vectorint parent; std::vectorMergeRecord mergeRecords; double euclidean(const std::vectordouble a, const std::vectordouble b) const { double sum 0.0; for (size_t i 0; i a.size(); i) { double diff a[i] - b[i]; sum diff * diff; } return std::sqrt(sum); } int findRoot(int x) const { while (parent[x] ! x) { x parent[x]; } return x; } };我特意把run()拆到类外面实现因为它的逻辑较长放外面更好注释。距离矩阵用二维vector存虽然连续内存上可能不如一维数组高效但可读性好且对于教学实现完全够用。初始化阶段做的事情很简单计算所有样本之间的欧氏距离存入一个对称矩阵。这里有个小优化——对角线不用算我们只使用上三角和下三角。std::vectorint AgnesClusterer::run() { const int n static_castint(data.size()); if (n 0 || k 0 || k n) { std::vectorint empty; return empty; } // 1. 欧氏距离矩阵 std::vectorstd::vectordouble dist(n, std::vectordouble(n, 0.0)); for (int i 0; i n; i) { for (int j i 1; j n; j) { double d euclidean(data[i], data[j]); dist[i][j] d; dist[j][i] d; } } // 2. 初始化簇状态 parent.resize(n); for (int i 0; i n; i) parent[i] i; std::vectorint clusterSize(n, 1); std::vectorbool active(n, true); int clusterCount n; mergeRecords.clear(); // 3. 主循环合并直到剩余簇数等于 k while (clusterCount k) { // 4. 查找距离最近的两个活跃簇 double bestDist std::numeric_limitsdouble::max(); int bestI -1, bestJ -1; for (int i 0; i n; i) { if (!active[i]) continue; for (int j i 1; j n; j) { if (!active[j]) continue; if (dist[i][j] bestDist) { bestDist dist[i][j]; bestI i; bestJ j; } } } if (bestI 0) break; // 理论上不会出现防御性检查 // 5. 记录合并历史合并 j 到 i int si clusterSize[bestI]; int sj clusterSize[bestJ]; MergeRecord rec{bestI, bestJ, bestDist, si sj}; mergeRecords.push_back(rec); parent[bestJ] bestI; clusterSize[bestI] si sj; active[bestJ] false; // 6. 更新合并簇到其他簇的距离 for (int m 0; m n; m) { if (!active[m] || m bestI) continue; double newDist; if (linkage SINGLE) { newDist std::min(dist[bestI][m], dist[bestJ][m]); } else if (linkage COMPLETE) { newDist std::max(dist[bestI][m], dist[bestJ][m]); } else { // UPGMA 加权平均避免简单平均的偏差 newDist (si * dist[bestI][m] sj * dist[bestJ][m]) / (si sj); } dist[bestI][m] newDist; dist[m][bestI] newDist; } --clusterCount; } // 7. 展平标签并输出 std::vectorint labels(n); for (int i 0; i n; i) { labels[i] findRoot(i); } // 把根标签重编号为 0, 1, 2, ..., clusterCount-1方便使用 std::vectorint renumberDict; std::vectorint result(n); for (int i 0; i n; i) { int root labels[i]; auto it std::find(renumberDict.begin(), renumberDict.end(), root); int newLabel; if (it renumberDict.end()) { newLabel static_castint(renumberDict.size()); renumberDict.push_back(root); } else { newLabel static_castint(it - renumberDict.begin()); } result[i] newLabel; } return result; }3.2 主循环里最烧脑的更新逻辑上面代码第6步的更新逻辑是整个AGNES实现里最容易被写错的地方。注意几个细节第一更新距离时我们用si和sj保存的是合并前的簇大小而不是合并后的。因为UPGMA公式里si * dist[bestI][m] sj * dist[bestJ][m]需要用到合并前的统计量。如果你在合并后才读取clusterSize[bestI]这个值已经变成了sisj代入公式会算错。我在第一次实现时就在这里踩了坑后面会细讲。第二对于单链接和全链接距离更新只涉及dist[bestI][m]和dist[bestJ][m]两个值不需要簇大小。但有一个隐藏问题单链接的min操作如果遇到距离矩阵中保存的是非最近距离因为之前合并更新过这时候dist[bestI][m]和dist[bestJ][m]本身已经是簇间距离的表示直接取min是没问题的。这一点如果你用原始点对距离去求反而是错的。第三active[bestJ] false意味着bestJ这一行一列在后续查找中会被跳过但距离矩阵里那一行一列的数据还在。这对于我们这种合并后把新簇保留在bestI行、bestJ行作废的策略是必要的。如果你的实现选择创建一个全新的簇对象来代表合并结果那么距离矩阵需要动态增删行列代码复杂度会明显上升我实测下来反而不如上面的固定行号active标记方案好调试。3.3 findRoot为什么不能偷懒直接改label合并操作执行的是parent[bestJ] bestI这会让parent数组形成一棵棵树。假设有5个点合并顺序是(1,2)、(3,4)、(2,3)那么最终parent可能是2的父是21的父是24的父是33的父是2。这时你不能再简单地看parent[i]来判定簇归属而必须向上追溯到根。findRoot函数的循环实现就是从任意节点出发一直向上走直到parent[x] x为止。有朋友问过我为什么不直接维护一个label数组每次合并把所有label等于bestJ的元素改成bestI对于小数据量这样确实直观。但在合并次数多、层级深的情况下每次全量扫描label数组是O(n)整个算法就是O(n^3)里再乘一个n变成了O(n^4)。虽然对于几千个点来说也不至于完全跑不动但在可扩展性上远不如并查集思想。而且合并历史记录里保存的左右簇编号也天然依赖这个parent结构。3.4 合并历史记录的用途mergeRecords数组记录每一次合并的左右簇编号、合并距离和合并后的簇大小。这个数组可以说是整棵树的压缩存档用途有两个第一生成树状图。每个合并记录对应一个内部节点它的左右孩子分别是left和right簇最终代表的那棵子树distance就是高度。拿这个数组可以直接画出一棵标准的横向树状图。我在第5章的扩展部分会给出绘制思路。第二辅助选择k。你不需要重新运行算法只需要翻看合并历史找出合并距离突然跳变的那个位置就能知道合理的簇数。比如某次合并距离是2.1下一次突然变成8.7说明此次合并把两个相对紧凑的簇强行拼成了一个松散的大簇这时候就应当把k定在跳变之前。这个思路在实践里非常有用。4. 实测结果与踩坑复盘三种链接策略的差异4.1 模拟数据构造与结果对比为了验证实现我用两组高斯分布生成了一个二维测试集一组200个点中心在(2, 2)另一组100个点中心在(8, 8)同时往数据里撒了5个远离两簇的噪声点。每类数据的分布都经过标准化量纲一致。分别用SINGLE、COMPLETE、AVERAGE三种链接跑了一遍设定k2结合实际输出总结如下链接方式主体簇划分噪声点归属树状图特征单链接能正确分开两个高斯簇常被拉进某个簇会出现明显的链式结构全链接能正确分开噪声会形成单独簇或影响合并顺序簇间合并距离偏大平均链接正确且稳定噪声影响较小合并层次比较均衡这个结果和我经验里的预期一致。单链接在这组数据上虽然聚类结果正确但观察树状图可以发现噪声点与某个簇的距离和簇间距离接近导致合并次序不稳定。如果数据里噪声更多单链接很容易把两个簇靠桥接点连起来产生长条噪声道。这也是为什么很多系统里默认不用单链接。4.2 踩坑一简单平均链接写出来的结果和标准实现对不上这是我踩得最久的一个坑。最开始我图省事在AVERAGE分支里写的是newDist (dist[bestI][m] dist[bestJ][m]) / 2.0;跑出来的结果在大多数平缓数据上和UPGMA差别不大但只要两类簇的样本量差距明显就会出问题。我给一个极端测试A簇200个点B簇5个点C簇100个点让A和C的距离很近、B和C很远。简单平均结果在合并A、B之后新簇到C的距离被B拉高导致下一步合并对象完全是另一个簇而UPGMA不会出现这种少数点绑架多数点的问题。对比代码只有一行之差但你去和标准库结果核对聚类标签时差异会非常明显。建议所有实现AGNES的朋友直接用UPGMA公式不要贪图那点代码简洁。4.3 踩坑二距离矩阵更新顺序错了导致整个树错位另一个关键坑是更新顺序。官方流程是先合并再更新新簇到其他簇的距离。这个顺序在概念上没有问题但如果你在更新距离之前修改了clusterSize[bestI]就会像我在3.2节提到的那样si和sj已经不再是合并前的值。我当时是用一个调试打印才发现问题的某个簇的簇大小在两次迭代里从10跳到了21第二次迭代里它到其他簇的距离参考值全部偏大。解决的思路是在进入更新循环之前先把si和sj用局部变量保存下来所有更新公式一律引用局部变量而不是引用clusterSize。这样即便后面改了clusterSize[bestI]也不影响公式的准确性。这种先把状态快照拿出来的习惯在处理很多增量更新场景时都很管用。4.4 特征缩放的隐形影响第一个版本我直接拿原始特征跑欧氏距离结果发现聚类结果总是被数值最大的那个维度主导。比如一份数据里第一个特征是年龄20到60第二个特征是消费金额1000到10000那么欧氏距离基本由消费金额决定年龄特征形同虚设。这不是算法的问题是距离度量的度量衡问题。解决方法是聚类前对每个特征做标准化最常见的做法是z-score标准化(x - mean) / std。更简单粗暴一点也可以用min-max归一化把每个特征压缩到[0,1]。这步操作发生在距离矩阵计算之前属于数据预处理不属于聚类算法本身但忘了它AGNES的效果会大打折扣。5. 性能瓶颈与优化路线什么时候该放弃纯AGNES5.1 复杂度来源O(n^3)到底卡在哪朴素AGNES的三层复杂度分别来自初始化距离矩阵O(n^2)。主循环最多跑n-k轮每轮查找最近簇对是O(n^2)。每轮更新距离矩阵是O(n)。所以总复杂度是O(n^3)存储是O(n^2)。听起来很吓人但实际上n在1000以内完全没有问题我测过1500个点平均链接版本几秒内跑完。真正开始卡的是n上到5000以上这时候每轮查找都要比较几千万次主循环数千轮就非常吃力了。5.2 优先队列优化的真实体验一个直观的优化是用优先队列最小堆维护每个簇的最近邻居这样取最近簇对可以做到O(1)或接近O(log n)。听起来很美但实践中有个大麻烦每次合并之后涉及bestI和bestJ到其他簇的距离全变了堆里存着大量过期的pair。要处理过期条目常见方案是惰性删除——pop时发现堆顶pair对应的距离和当前距离矩阵不一致就丢掉重新pop。这个做出来代码复杂度上升一个台阶而且合并后更新堆时还需要遍历所有相关簇实际收益远没有想象中大。数据量在5000以内我建议乖乖用O(n^2)扫描配合编译器优化速度和正确性都好控制。5.3 更现实的优化方案如果你确实面对上万样本我建议不要死磕纯AGNES而是考虑下面几条路先抽样再用AGNES比如先用k-means或随机采样得到几千个代表点做层次聚类再把剩余点分配给最近的簇。这个思路本质上是粗聚类细分配实践中非常常见。使用BIRCH算法它专门为大样本设计通过CF树先做数据压缩再在叶节点上做层次聚类能处理百万级数据。分块策略把数据按维度分块用分布式框架并行计算距离矩阵但实现门槛很高小项目不划算。5.4 树状图的可视化思路最后给想看层次结构的朋友一个简单方案把mergeRecords输出成文本格式如下节点ID, 左孩子, 右孩子, 合并距离, 合并后规模比如300, 12, 87, 1.234, 23 301, 300, 45, 2.100, 56用任意支持树状图绘制的脚本读入这个文件就可以画出横向树状图。即使不画图你也可以在控制台里按合并距离排序输出观察合并距离跳跃点这就是选k的直接依据。我在项目里就是靠这份文本记录来确定最终分群数的比盯着扁平标签猜k要直观得多。5.5 现阶段实现的不足与扩展方向这套实现的不足也很明确它不能自动选k、不支持加权样本、只实现了欧氏距离。如果你想进一步扩展可以按这个优先级来增加距离函数抽象把欧氏距离换成余弦距离、曼哈顿距离、马氏距离。支持样本权重让部分关键样本对聚类结果有更大影响。在run()里增加一个最大合并距离参数当最近簇对距离超过阈值时提前停止这样连k都不用传。其中距离函数抽象对实际业务最有用因为我遇到过不少场景特征不是连续数值而是高维稀疏向量这时候欧氏距离效果很差需要换成余弦距离。改造方法也不复杂把euclidean函数改成函数指针或模板参数即可。最后提一下我个人的使用体会。写这个C版本的AGNES最初只是为了应付一个嵌入式相关的聚类需求系统里没有Python环境只能靠C搞定。写完跑通之后我反而对层次聚类有了更直观的理解它慢、它朴素但它的输出你可能绕了好几圈还是想要那一棵树。如果你也正打算在C里实现这个算法先把我踩过的这两个公式坑绕开——UPGMA别用简单平均簇大小快照记得先保存——剩下的事情按代码一步步走就能跑通。