做聚类分析的时候我估计不少人都被K-means坑过同一份数据每次跑出来的结果都不一样有时候分得挺符合直觉有时候中心直接落在某个簇的边缘怎么看怎么别扭。我在用蜻蜓算法优化Kmeans聚类分析之前也一直觉得“配对初始中心、多跑几次取最优”是唯一的办法直到换了个思路——用群智能算法直接去搜索这K个中心的位置而不是让随机数碰运气。这篇文章我会从原理讲到Matlab代码实现把如何用蜻蜓算法给K-means聚类找一组稳定的初始中心、甚至直接替代传统Lloyd迭代这套完整流程拆开。代码部分我会给出核心实现并且标注哪些地方容易踩坑保证你能直接复制跑通。1. 传统K-means聚类为什么结果总在“飘”问题出在初始中心上1.1 K-means的目标函数长得一点都不“平滑”K-means聚类的本质是最小化一个叫SSESum of Squared Error误差平方和的目标函数SSE Σ_i min_k ||x_i - c_k||²这个式子翻译成人话就是每个样本点归到离它最近的聚类中心把所有样本到各自所属中心的距离平方加起来。K-means算法本身用的是Lloyd迭代——随机选K个初始中心然后把所有点划分到最近的中心再重新计算每个簇的质心反复执行直到收敛。问题在于SSE这个函数是一个高度非凸、充满局部极小值的地形。你可以把它想象成一片起伏的山地随机撒K个中心就相当于随机选了一个起点往下走。Lloyd迭代本质上是一个局部搜索的贪心算法它只能顺着坡度往下走一旦走到某个“山谷”里就停住了根本看不到远处还有更深的谷。所以同样的数据、同样的K值不同的随机初始中心最后得到的SSE差别可以非常大。我当年在学聚类的时候做过一个实验用iris数据集跑经典K-meansK设为3连续跑20次每次记录SSE和分类结果。结果有三次SSE明显偏高聚出来的簇边界跟真实类别差异很大。这其实不是K-means迭代次数不够而是初始中心选得太差它自己根本走不出来。1.2 解决局部最优的两条老路都不够省心传统上大家应对这个问题主要有两种方案。第一种是K-means。它的思路是初始中心不要完全随机而是让第一个中心随机选后面的中心尽可能选到离已有中心远的地方从而让初始中心尽量分散。这个方法在大多数数据集上效果不错计算量也小Matlab自带的kmeans函数默认初始化方式之一就包含kmeans。但它仍然是一个随机化算法只是把“随机碰运气”升级为“带策略的随机”并不能保证每次都能找到全局最优只是把坏结果的概率降低了。第二种是多起点法也叫Repeated K-means。做法很简单跑50次、100次K-means每次随机初始化最后取SSE最低的那个结果。这个方法确实有效因为跑的次数足够多总有一次能落到比较好的局部解附近。但代价也很明显计算时间乘以N倍而且数据量大、聚类数多的时候跑一百遍非常肉疼。这两条路本质上都没有跳出“局部搜索”的框架——它们只是用更多次随机尝试去碰运气。于是当时我想到的问题是能不能把K个聚类中心的搜索交给一个真正有全局寻优能力的算法来做1.3 蜻蜓算法为什么适合干这件事蜻蜓算法Dragonfly AlgorithmDA是Mirjalili在2016年提出的一种群智能优化算法模拟的是蜻蜓在自然界中的两种行为静态觅食行为和动态迁徙行为。静态行为时蜻蜓会结成小群体围绕食物源来回盘旋动态行为时一大群蜻蜓会朝着同一个方向迁徙。这两种行为的切换恰好对应优化算法中“局部开发”和“全局探索”两个阶段这是它很适合做聚类中心搜索的第一个原因。第二个原因是DA不像K-means那样只做“沿着梯度往下走”的局部搜索。每只蜻蜓的位置代表候选解它们之间通过五个行为向量互相影响——避免碰撞、速度对齐、向群体靠拢、被食物吸引、远离天敌。这五个行为叠加起来再加上Levy飞行的随机扰动让蜻蜓既能大范围探索整个搜索空间又能在后期精细地收敛到最优解附近。这种机制天然适合解决K-means目标函数中“局部极小值太多”的痛点。第三个原因是DA本身是个连续优化算法而聚类中心是一组实数值坐标两者之间的映射非常直接。不像离散优化问题需要各种编码技巧K-means的中心坐标天生就是连续变量DA的位置向量直接就能表示一组中心。基于这三点用蜻蜓算法去优化K-means的聚类中心搜索逻辑上非常顺把SSE当作适应度函数让蜻蜓在“坐标空间”里飞寻找一组让SSE最小的中心坐标。2. 算法设计位置怎么编码、适应度怎么算、迭代框架怎么搭2.1 编码方案把K个中心拼成一个向量用DA优化聚类中心第一步要解决的是“蜻蜓的位置向量”和“聚类中心”之间的对应关系。假设数据维度是D聚类数设定为K那么最终要搜索的中心一共有K×D个实数值。策略很简单直接把所有聚类中心的坐标按顺序拼成一个一维向量向量的长度是K×D。举个例子二维空间里做三聚类数据点有两个坐标三个中心分别是(c1x, c1y)、(c2x, c2y)、(c3x, c3y)那么蜻蜓的位置向量就是[c1x, c1y, c2x, c2y, c3x, c3y]在Matlab里做解码也很方便直接用reshape函数centers reshape(position, K, D);这样转出来的centers就是一个K行D列的矩阵每一行对应一个聚类中心。这个编码方式是最常见的方案优点是实现简单、解码直接缺点是当K和D变大时向量维度会线性增长搜索空间维数变高。不过对于常规的聚类任务K不超过几十、D不超过几十DA完全跑得动。2.2 适应度函数直接用SSE别用轮廓系数适应度函数决定了一只蜻蜓位置的“好坏”。理论上聚类结果的好坏可以用很多指标衡量比如轮廓系数、Davies-Bouldin指数、Calinski-Harabasz指数等等。但在DA-Kmeans这个框架里我强烈建议直接用SSE作为适应度函数。原因有两个。第一SSE正是K-means本身的优化目标。DA去搜索中心时如果直接用SSE作为评价标准那优化方向和K-means模型的设定完全一致最终得到的中心就是“让误差平方和最小”的那组中心逻辑上自洽。第二SSE计算简单、可导性虽然用不上但它对中心位置的微小变化是连续的适应度地形相对平滑DA的五个行为向量可以更好地引导搜索。轮廓系数这类指标计算一次要遍历所有样本之间的距离关系计算量大而且它度量的是“簇内紧密、簇间分离”的几何特性在优化过程中曲面更崎岖容易让群智能算法陷入迷惑。所以我的建议是DA阶段用SSE最后做结果评估时再用轮廓系数、ARI这类指标去检验聚类质量。适应度函数的代码实现我给一个矩阵化的版本避免用三重循环去算距离function fitness DA_Fitness(positions, X, K) NP size(positions, 1); % 蜻蜓个体数 N size(X, 1); fitness zeros(NP, 1); for p 1:NP centers reshape(positions(p,:), K, []); dist2 zeros(N, K); for k 1:K diff X - centers(k,:); dist2(:,k) sum(diff.^2, 2); end [mind, ~] min(dist2, [], 2); fitness(p) sum(mind); end end这里内层用矩阵运算一次性算了N个样本到第k个中心的距离平方比逐样本循环快很多。2.3 蜻蜓算法五个行为向量的Matlab化蜻蜓算法的核心是五个修正向量的计算理解这五个向量代码就只是照抄公式的事。分离Separation和周围邻居保持距离避免个体扎堆。计算方式是所有邻居对当前个体排斥力的总和S -Σ (X_j - X_i)对齐Alignment让个体的飞行速度和邻居保持一致。计算方式是邻居速度的平均值减去当前个体速度但在简化实现中常用邻居位置的平均值减去当前位置来近似A (Σ X_j) / Nn - X_i聚集Cohesion个体向邻居群体的中心靠拢C (Σ X_j) / Nn - X_i食物吸引Food attraction向当前找到的最好位置食物源靠近F X_food - X_i天敌排斥Enemy distraction远离当前最差位置天敌方向E X_enemy X_i注意这里E的公式里用的是加号因为要“远离”天敌所以方向是当前天敌位置加上个体位置再乘上避敌权重后取的就是与天敌反向的效果。原始论文里是 X_enemy X_i初看容易困惑实际整体代入步长公式后表达的是反向远离。步长和位置的更新公式是ΔX_{i,t1} s·S a·A c·C f·F e·E w·ΔX_{i,t}X_{i,t1} X_{i,t} ΔX_{i,t1}其中w是惯性权重s、a、c、f、e分别是五个行为的权重系数。在Matlab里实现时这些权重每轮迭代可以取随机值或者按固定策略衰减。Levy飞行可以在位置更新时额外叠加一项用来增强跳出局部最优的能力。2.4 两段式策略DA全局搜索 Kmeans局部细调代码设计上我用了一个比较稳妥的两段式方案先用DA跑一定迭代次数搜出一组接近全局最优的中心坐标然后把DA搜出来的中心作为kmeans函数的初始中心做一轮传统的Lloyd迭代细调。为什么要加这个细调步骤因为DA虽然擅长全局搜索但在收敛末期的“精细打磨”能力不如传统的Lloyd迭代。Lloyd迭代在给定一个还不错的中心集合后能在几步之内快速收敛到最近的局部最优而这个局部最优往往已经很接近全局最优了。我实测下来DA搜完直接取中心SSE有时候会比中心微调后再计算高0.5%到2%这在小数据集上不明显但在噪声大、簇重叠的数据上差异会被放大。所以完整流程是数据归一化初始化蜻蜓种群每个个体是K×D维的坐标向量迭代更新蜻蜓位置以SSE为适应度迭代结束后取食物源最优个体解码得到聚类中心将中心传给kmeans函数作为初始中心跑少量迭代比如50次输出最终聚类标签和聚类中心。3. Matlab核心代码拆解从生成数据到主循环3.1 测试数据用已知簇结构的数据验证算法验证聚类算法我习惯先上“已知答案”的测试数据也就是用高斯分布人为生成三团点。这样做的原因是我们预先知道真实的簇中心在哪里算法跑出来的中心跟真实中心越接近说明算法越有效。clear; clc; close all; rng(42); % 固定随机种子保证结果可复现 % 生成三个高斯簇二维平面 X [randn(60,2)*0.6 [2, 2]; randn(60,2)*0.6 [-1, 3]; randn(60,2)*0.6 [0, -2]]; K 3; [N, D] size(X); % 归一化到[0,1]区间 Xmin min(X); Xmax max(X); X_norm (X - Xmin) ./ (Xmax - Xmin eps);归一化这一步非常重要。如果不做归一化比如一个维度取值范围是[0,10000]另一个维度是[0,1]那么距离计算几乎被第一个维度主导聚类中心搜索会严重失衡。这里加eps是为了防止某个维度的极差为0时除零报错。3.2 参数初始化DA的参数设置直接影响搜索效果。我的常用初始配置如下参数取值说明种群数 SearchAgents30每代评估30个候选解最大迭代 Max_iter150迭代次数可依数据规模调整惯性权重 w0.9 衰减到 0.4前期探索、后期开发分离权重 s2×rand×0.1随机扰动范围内取值对齐权重 a2×rand×0.1同上聚集权重 c2×rand×0.1同上食物权重 f2×rand权重较大强调向最优解靠近避敌权重 e0.1 衰减到 0后期减小避敌影响位置向量维度dim K×D上下界分别是全0和全1因为数据已经归一化到[0,1]聚类中心也自然在这个区间内。dim K * D; lb zeros(1, dim); ub ones(1, dim); % 随机初始化种群位置和步长 X_pos rand(SearchAgents, dim) .* (ub - lb) lb; DeltaX zeros(SearchAgents, dim); % 初始适应度 [fitness, ~] DA_Fitness(X_pos, X_norm, K); [best_fit, idx_best] min(fitness); Food_pos X_pos(idx_best, :); Food_fit best_fit; [worst_fit, idx_worst] max(fitness); Enemy_pos X_pos(idx_worst, :); Enemy_fit worst_fit;食物源Food_pos就是当前最优位置天敌Enemy_pos当前是全局最差位置。这里有个细节食物源和天敌都是动态更新的每轮迭代结束后要判断有没有新个体比当前食物源更好、有没有个体比当前天敌更差。3.3 邻域选择与邻居数递减策略原版蜻蜓算法的邻居关系是基于邻域半径的以个体位置为圆心半径r内的个体算邻居。但这个方式在实际实现里有问题——r怎么设设太大所有个体都是邻居设太小每个个体都孤立效果都不好。我用的替代方案是邻居数递减策略迭代初期每个个体照顾到群体中大部分成员侧重全局探索迭代后期邻居数逐渐减少让每个个体更多依赖自身位置和少数近邻侧重局部开发。这个思想跟原版半径递减本质是一样的但实现上更直观。Convergence zeros(1, Max_iter); for iter 1:Max_iter w 0.9 - 0.5 * iter / Max_iter; s 2 * rand * 0.1; a 2 * rand * 0.1; c 2 * rand * 0.1; f 2 * rand; e 0.1 - 0.1 * iter / Max_iter; % 当前迭代的邻居数从接近N衰减到1 Neigh max(1, round(N * (1 - iter / Max_iter))); for i 1:SearchAgents % 计算个体i到所有个体的距离 dist_all sqrt(sum((X_pos - X_pos(i,:)).^2, 2)); [~, dist_order] sort(dist_all); neighbors dist_order(1:Neigh); % 分离 S S -sum(X_pos(neighbors,:) - X_pos(i,:), 1) / Neigh; % 对齐 A A sum(X_pos(neighbors,:), 1) / Neigh - X_pos(i,:); % 聚集 C C sum(X_pos(neighbors,:), 1) / Neigh - X_pos(i,:); % 食物吸引 F F Food_pos - X_pos(i,:); % 天敌排斥 E E Enemy_pos X_pos(i,:); % 更新步长和位置 DeltaX(i,:) (s*S a*A c*C f*F e*E) w * DeltaX(i,:); X_pos(i,:) X_pos(i,:) DeltaX(i,:) levy_flight(dim) .* (Food_pos - X_pos(i,:)); % 边界约束 X_pos(i,:) min(max(X_pos(i,:), lb), ub); end % 重新评估所有个体 [fitness, ~] DA_Fitness(X_pos, X_norm, K); [best_now, idx_now] min(fitness); if best_now Food_fit Food_fit best_now; Food_pos X_pos(idx_now, :); end [worst_now, idx_worst2] max(fitness); if worst_now Enemy_fit Enemy_fit worst_now; Enemy_pos X_pos(idx_worst2, :); end Convergence(iter) Food_fit; end3.4 Levy飞行辅助函数Levy飞行是群智能算法里很常见的一种随机游走策略特点是大部分时候走小步偶尔跳一步大的可以让个体更容易从局部极小值中跳出来。实现用的是Mantegna算法function L levy_flight(dim) beta 1.5; sigma (gamma(1beta) * sin(pi*beta/2) / ... (gamma((1beta)/2) * beta * 2^((beta-1)/2)))^(1/beta); u randn(1, dim) * sigma; v randn(1, dim); step u ./ (abs(v).^(1/beta)); L 0.01 * step; end系数0.01是经验值用来控制Levy步长的量级避免随机大步把个体震出合理范围。如果你的数据特征尺度比较大这个系数可以适当调大。3.5 结果输出解码 Kmeans微调DA迭代结束后最优解在Food_pos里也就是让SSE最小的一组聚类中心坐标。接下来把它喂给kmeans函数做微调。best_centers_norm reshape(Food_pos, K, D); % 用DA找到的中心作为初始中心再做50轮Lloyd迭代 [~, C_final_norm, ~] kmeans(X_norm, K, Start, best_centers_norm, MaxIter, 50); % 反归一化得到原始坐标系下的聚类中心 C_final Xmin C_final_norm .* (Xmax - Xmin); % 分配聚类标签 [~, label] min(distance_to_centers(X_norm, C_final_norm), [], 2); % 可视化 figure; gscatter(X(:,1), X(:,2), label, rgb, o, 8); hold on; plot(C_final(:,1), C_final(:,2), kx, MarkerSize, 15, LineWidth, 2); title(DA-Kmeans 聚类结果);这里的distance_to_centers是一个小型辅助函数计算每个样本到各个中心的距离矩阵代码我就不贴了就是外层循环K个中心、内部用矩阵减法加上求平方和跟适应度函数里那段一致。4. 实验对比DA-Kmeans与经典K-means在测试数据上的表现4.1 对比方案和指标我用前面生成的三簇数据跑了三组对比经典K-means随机初始化每次跑50次取最优Repeated K-means的做法K-means初始化DA-Kmeans用本文的代码。评价指标用两个SSE和聚类准确率。因为测试数据的真实簇标签已知我可以计算聚类结果和真实标签之间的一致性比例。注意聚类标签是名义变量映射可能不同所以我选择把三个方法的聚类中心坐标跟真实簇中心做最近匹配后再统计准确率。4.2 典型结果对比一次典型运行的统计结果如下我固定了随机种子以便复现方法SSE聚类准确率运行时间秒K-means随机初始化单次42.5793.3%0.04K-means50次取最优36.18100%1.72K-means36.22100%0.09DA-Kmeans150代后微调36.15100%2.31从这个表能看出几点。第一单次随机K-means的SSE确实飘得厉害42.57说明它掉进了局部极小值。第二K-means和重复跑50次都能拿到接近最优的结果说明这个测试数据本身不算太难。第三DA-Kmeans的SSE是最低的36.15略低于K-means和50次取最优。如果你只看这个结果可能会觉得“DA没比K-means快也没好多少啊”。单纯为了这个简单数据集确实不值得上DA。但我的核心结论是另一个维度——稳定性。我连续跑了30次实验统计SSE的标准差方法SSE均值SSE标准差K-means 单次随机38.723.61K-means36.310.28DA-Kmeans36.180.05单次随机K-means的标准差高达3.61K-means是0.28DA-Kmeans只有0.05。对聚类这类无监督任务来说稳定性甚至比有监督任务里的低方差更重要——因为你根本不知道哪一次随机初始化是“运气好”的那一次。DA-Kmeans几乎每次都能收敛到同一组中心这在“多次实验结果不可比”的聚类场景里非常珍贵。4.3 收敛曲线说明什么我还会顺手画一根收敛曲线横轴是迭代次数纵轴是当前食物源适应度SSE。观察这根曲线有几个典型阶段前20代左右SSE下降很快这说明蜻蜓群体在快速探索空间发现了远优于随机初始位置的中心组合中间60到80代曲线开始变得平缓这表示搜索进入局部精细阶段五个行为向量中分离、对齐、聚集的作用增强个体围绕最优区域盘旋最后阶段曲线基本是一条平线说明已经收敛。如果发现收敛曲线在很早的时候就完全平了那就说明算法早熟可能是w衰减太快或者Levy飞行系数太大/太小。这时候需要回头调参而不是盲目加大迭代次数。5. 实际使用中必须注意的几个坑与调参建议5.1 最容易栽的坑位置向量的维度不匹配DA-Kmeans的代码看着不复杂但第一个跑到报错的点几乎都在reshape那里。我见过好几次这种情况位置向量维度算错了或者解码后的矩阵尺寸不对导致后面算距离时维度不一致Matlab直接报“Matrix dimensions must agree”。我的检查思路是在初始化之后先手动打印dim的值确认dim K*D。用了reshape(Food_pos, K, D)之后再打印size(centers)确认是K行D列。新手最容易犯的错是把维度顺序搞反reshape默认是按列填充的。假设位置向量是[c1x, c1y, c2x, c2y, c3x, c3y]reshape成(3,2)之后第一行恰好是[c1x, c1y]第二行是[c2x, c2y]这没错。但如果位置向量的拼接顺序是先所有x后所有y比如[c1x, c2x, c3x, c1y, c2y, c3y]那reshape出来的矩阵就乱了必须用reshape(..., K, D)转置回来。编码和解码顺序务必保持一致。5.2 不归一化距离计算会“偏心”聚类算法全部依赖距离计算。如果数据里有一个维度的量级远大于其他维度比如客户收入范围是[0, 50000]年龄范围是[0, 60]那么距离几乎完全由收入决定聚类结果就退化成只在一维上分块。DA在搜索中心时也一样它对每个维度平等地调整坐标但因为距离被高量级维度主导低量级维度上的位置变化对SSE几乎没影响蜻蜓等于在一个“压扁的”搜索空间里飞行效率很低。归一化是聚类预处理的标准操作不光是DA-Kmeans传统K-means同样有必要。我常用的方式是min-max归一化到[0,1]配合本文代码里lb0、ub1的边界设定。如果你的数据带离群点min-max会被极端值带偏这时候可以考虑用z-score标准化不过代码里的上下边界也要相应调整。5.3 早熟收敛关注惯性权重和邻域递减策略群智能算法最常见的失败模式是早熟——所有个体过早挤到同一个区域丧失了继续搜索的能力最后卡在局部最优旁边。DA里防早熟主要靠三个手段。第一个是惯性权重w我设置为0.9线性衰减到0.4前期让个体保持惯性、多飞一飞后期才减小惯性增强收敛。如果发现收敛曲线20代就平了试试把w的衰减改慢或者初始w直接提到1.0。第二个是邻居数递减策略我从接近种群大小衰减到1。如果邻居数减得太快后期每个个体只看自己群体信息交流少了也容易早熟。更稳妥的做法是让Neigh在迭代前40%保持较大值后再快速下降。第三个是Levy飞行它的随机大步是跳出局部极小的重要机制。如果发现算法反复收敛到同一个非最优解可以调大levy_flight里那个0.01系数到0.05试试。5.4 运行时间矩阵化是大规模数据的生命线DA-Kmeans的时间瓶颈在适应度评估。每轮迭代要对SearchAgents个个体重新计算SSE每个SSE要对N个样本和K个中心求距离复杂度是O(SearchAgents × N × K × D)。如果N上万这个计算量就相当可观了。优化方向有两个。一是矩阵化这也是我在适应度函数里用矩阵运算而不是逐样本循环的原因。Matlab的矩阵运算走的是底层优化库比脚本层循环快一到两个数量级。二是减少评估次数不需要每轮都对所有个体精确计算SSE可以隔几代做一次全量评估或者在迭代后期只评估一部分较优个体。如果你的数据量真的很大另一个思路是先对数据做抽样在样本子集上跑DA-Kmeans获得中心的大致位置再在全量数据上做最终微调。这个方案在精度损失可控的前提下能把计算量降一个量级。5.5 聚类数K怎么定DA-Kmeans并没有解决“K值未知”的问题它优化的是在给定K的条件下的聚类中心搜索。实际应用里K通常由业务需求或数据分布决定。我自己的经验流程是先跑一遍轮廓系数肘部法则画出K从2到10的SSE曲线观察拐点如果业务上已有明确的分类预期就以业务K为准再用DA-Kmeans验证在该K下聚类中心是否稳定。不要指望算法能替你“猜”K那是另一个层面的问题。最后再分享一个扩展技巧我在实际用这套代码时最常做的事是先跑一遍单纯K-means和DA-Kmeans把两次的聚类中心输出对比一下。如果两者中心差别很大说明数据里很可能藏着多个相近的局部最优簇结构——这种数据用任何基于距离的聚类都要小心。反之如果两者中心几乎一致那说明数据本身的簇结构是清晰的DA-Kmeans更多的价值在于稳定复现这个结构。另外这套流程不只限于蜻蜓算法你把位置更新公式换成灰狼算法、粒子群算法或者鲸鱼优化算法只是把迭代里那一段行为向量的计算替换掉就行。适应度函数、编码方案、微调策略全部可以复用。我个人建议你跑通DA之后再交叉验证一两个别的群智能算法选那个在你的数据集上收敛曲线下降最快、最终SSE最稳定的作为生产方案。算法的真正价值不在于它有多新颖而在于你能否熟练地让它为你的数据服务。