做聚类分析最让人头疼的事是什么不是算法不会写而是跑了三遍 K-means出来三套完全不同的结果。我去年做客户分群的时候就被这个问题卡了好几天同一份数据同一段代码只是没固定随机种子业务那边看到的分群结果就改了。后来我把目光转向了群智能优化算法用蜻蜓算法去给 K-means 找一组靠谱的初始质心才算把这个问题摁住。这篇东西就是记录我当时的完整做法蜻蜓算法Dragonfly AlgorithmDA怎么跟 K-means 结合、优化目标怎么设计、Matlab 代码怎么组织、以及跑实验时踩到的各种坑。如果你正在做聚类分析、写论文需要对比实验、或者单纯想把聚类结果做得更稳定更漂亮这篇应该能直接给你一份能抄作业的方案。1. K-means的老毛病与蜻蜓算法的应对思路1.1 K-means为什么脆弱每次结果都不同根子在初始点先说一个很多人忽略的事实K-means 本质上是个坐标下降的迭代算法它优化的是一个非凸问题。目标函数——所有样本到所属簇中心的距离平方和——表面看着平滑实际上里面藏着大量局部最优解。算法从哪组初始质心出发决定了它最终掉进哪个局部坑里。默认的 K-means 初始化方法是随机选 K 个样本当质心或者用改进版的 K-means 来降低随机性。但无论哪种本质上都还是运气成分占主导遇到重叠度高的簇、样本密度不均匀的数据随机初始化很容易让两个质心挤到同一个簇里面另一个簇被撕成两半。我那次客户分群连续跑五次出来的轮廓系数从 0.32 到 0.51跨度大到没法跟业务方解释。解决思路有两条。第一条是暴力路线同一个数据集跑几十次 K-means保留 SSE 最小的那组结果。简单有效但浪费算力而且本质还是躲运气而不是解决运气。第二条就是我们今天说的路线用全局优化算法先去搜一组好质心再把这组质心交给 K-means 做局部精修。1.2 为什么是蜻蜓算法和 PSO、遗传算法的差别在哪群智能算法里最出名的其实是粒子群PSO和遗传算法GA我做选择之前也对比了一圈。蜻蜓算法是 Mirjalili 在 2016 年提出的灵感来自蜻蜓捕食时的动态飞行行为。它最吸引我的一点是把搜索过程分成了探索和开发两个阶段并且通过邻居半径自适应切换。具体说蜻蜓有五个行为特征——分离避免跟邻居撞车、对齐跟邻居速度保持一致、凝聚向邻居中心靠拢、觅食飞向食物、避敌远离天敌。当一只蜻蜓周围没有邻居的时候它会做随机游走这对应全局探索当邻居密度上来以后五个力加权叠加对应局部精细搜索。这种机制比 PSO 单纯朝全局最优跑、GA 靠交叉变异碰运气要灵活一些。我当时也考虑过灰狼优化GWO但灰狼的机制相对简单处理高维连续优化问题时收敛速度快但容易早熟。DA 的参数不算多但行为组合丰富调参空间也比 PSO 更大。下面这个对比表是我自己做选型时整理的算法核心机制探索/开发平衡关键参数量适合场景K-means距离分配 质心更新无探索几乎无数据干净、初始化好PSO个体最优 全局最优牵引开发偏强少连续优化、快速收敛GA选择 交叉 变异机制平衡中离散问题、组合优化GWO种群层级 包围捕猎收敛偏快少中低维连续优化DA五力加权 邻居半径自适应动态切换中质心初始化、多峰问题如果你手头项目只追求省事PSO 也能用但如果你希望聚类结果更稳定、并且想写进论文里展示对比实验DA 的机制更有文章可做。2. 优化模型与编码设计蜻蜓算法到底在优化什么2.1 目标函数把聚类问题翻译成连续优化问题要让 DA 帮 K-means 干活第一步是把聚类问题改写成 DA 能理解的样子。K-means 的目标函数本来就很明确所有样本到它所属簇质心的欧氏距离平方和也就是 SSESum of Squared Errors。我们把 SSE 作为适应度函数DA 的目标就是找到一个质心组合让 SSE 尽可能小。公式不用怕就一句话SSE 每个样本到其最近质心的距离平方累加。DA 优化器要做的就是不断生成新的质心组合、算 SSE、再根据结果调整下一步搜索方向。这样 K-means 的问题就从给一组初始质心然后迭代变成了在连续空间里找一组能让 SSE 最小的质心。这里有个关键点必须说清楚DA 不直接给样本贴标签。它只负责输出哪 K 个质心是一组好质心样本归属哪个簇的标签还是交给 K-means 根据最近质心法则去算。2.2 质心向量编码一个蜻蜓就是一组完整质心DA 的每个搜索代理一只蜻蜓的位置是一个向量。这个向量怎么表达 K 个质心的坐标我用最直接的拼接方式把 K 个质心的所有特征依次连成一串。举个例子数据有 D 个特征、聚类数 K 3那么一个代理向量的长度就是 3 × D。前 D 个位置代表第 1 个质心中间 D 个位置代表第 2 个质心依此类推。在 Matlab 里这个向量 reshape 成 K × D 矩阵就能直接拿去做距离计算。边界约束也很自然每个特征的取值范围就是数据的 min 和 max。注意这里边界向量的长度也需要扩展成 K × D也就是每个簇的每个特征都要有一个上限和下限而不是只给原始特征的宽度。这个细节我第一次写代码时漏了结果边界矩阵维度对不上直接报错。2.3 蜻蜓的五个力分离、对齐、凝聚、觅食、避敌怎么算进位置更新DA 的核心更新公式不算复杂。每只蜻蜓在迭代时维护一个步长向量 ΔX更新规则是ΔX_new (s × 分离力 a × 对齐力 c × 凝聚力 f × 觅食力 e × 避敌力) w × ΔX_old然后位置更新就是 X_new X_old ΔX_new。五个力的含义如下分离力是本蜻蜓与所有邻居位置差的平均值取反方向是离邻居远一点对齐力是所有邻居速度也就是邻居的步长向量的平均值用来保持飞行同步凝聚力是邻居们的中心位置减去自己当前位置方向是向群体靠拢觅食力是全局最优位置减去当前位置把蜻蜓拉向当前找到的最优质心组合避敌力是自己当前位置减去全局最差位置让它远离当前表现最差的解。还有一个容易忽略的操作——邻居半径。DA 通过半径 r 判断哪些个体算邻居r 会随着迭代逐渐增大。初期 r 小邻居少蜻蜓更像各自探索后期 r 大邻居多信息共享加强开始精细收敛。很多复现没做这一步直接把所有个体当邻居那就退化成 PSO 了等于失去了这个算法的灵魂。3. Matlab代码逐段拆解主脚本、优化器与目标函数3.1 主脚本数据准备、参数设置与对比实验完整代码由三个文件组成主脚本、DA 优化器、目标函数。先看主脚本它负责生成/加载数据、设置 DA 参数、调用优化器、然后跟标准 K-means 做对比。% DA_Kmeans_main.m % 基于蜻蜓算法优化Kmeans聚类分析 —— 主脚本 clear; clc; close all; rng(2024); % 固定随机种子保证实验可复现 % 生成三簇重叠度中等的二维合成数据 X [randn(200,2)*1.2 [ 2 2]; randn(200,2)*1.2 [-2 2]; randn(200,2)*1.2 [ 0 -2.5]]; K 3; [N, D] size(X); % DA参数设置 N_DA 30; % 搜索代理数量蜻蜓个数 maxIter 200; % 最大迭代次数 dim K * D; % 每个代理的向量长度 质心数 * 特征数 % 边界向量每个簇每个特征范围 lb repmat(min(X, [], 1), 1, K); ub repmat(max(X, [], 1), 1, K); % 目标函数句柄输入一个质心向量输出SSE fobj (x) obj_fun(x, X, K); % 调用蜻蜓算法优化器 [bestPos, ConvergenceCurve] DA_optimizer(N_DA, maxIter, lb, ub, dim, fobj); % 将DA找到的最优质心向量解码 initCenters reshape(bestPos, K, D); % 用DA最优质心作为K-means初始点 [~, SSE_DA] kmeans(X, K, Start, initCenters, ... MaxIter, 1000, Replicates, 1); % 对比标准K-means随机初始化跑20次 R 20; sseStd zeros(R, 1); for r 1:R [~, sseStd(r)] kmeans(X, K, Start, plus, MaxIter, 1000); end % 输出结果 fprintf(DA-Kmeans SSE: %.4f\n, SSE_DA); fprintf(标准K-means高: best%.4f, worst%.4f, avg%.4f, std%.4f\n, ... min(sseStd), max(sseStd), mean(sseStd), std(sseStd)); % 画图收敛曲线 聚类结果 figure; subplot(1,2,1); plot(ConvergenceCurve, LineWidth, 1.6); xlabel(迭代次数); ylabel(SSE); title(DA收敛曲线); grid on; subplot(1,2,2); [idx, ~] kmeans(X, K, Start, initCenters, ... MaxIter, 1000, Replicates, 1); gscatter(X(:,1), X(:,2), idx); hold on; plot(initCenters(:,1), initCenters(:,2), kx, MarkerSize, 12, LineWidth, 2); title(DA-Kmeans聚类结果); grid on;这里有一个细节调用内置kmeans时Replicates必须设为 1。因为我们已经把初始质心传进去了如果还让Replicates保持默认的多次随机重启那前面的 DA 优化就白做了结果又变成运气好坏的叠加。如果你没有统计工具箱内置kmeans用不了可以用自己写的 Lloyd 迭代版本替代逻辑很简单给定初始质心后循环分配样本到最近质心、再按每个簇的均值更新质心直到质心不再变化。这个替代方案不影响 DA 部分的任何逻辑。3.2 DA优化器完整可运行的蜻蜓算法接下来是 DA 的核心封装函数输入种群规模、迭代次数、边界、维度和目标函数输出最优质心向量和收敛曲线。function [bestPos, ConvergenceCurve] DA_optimizer(N, maxIter, lb, ub, dim, fobj) % 蜻蜓算法优化器 % N: 搜索代理数量 % maxIter: 最大迭代次数 % lb, ub: 边界向量1 x dim % dim: 质心向量长度 % fobj: 目标函数句柄 % 初始化位置和步长 lbMat repmat(lb(:), N, 1); ubMat repmat(ub(:), N, 1); pos rand(N, dim) .* (ubMat - lbMat) lbMat; delta zeros(N, dim); % 权重系数个人经验值可调 s 0.1; % 分离 a 0.1; % 对齐 c 0.7; % 凝聚 f 1; % 觅食 e 1; % 避敌天敌 wMax 0.9; wMin 0.2; % 惯性权重线性递减 ConvergenceCurve zeros(maxIter, 1); for iter 1:maxIter % 惯性权重递减 w wMax - (wMax - wMin) * iter / maxIter; % 评估所有代理 Fitness zeros(N, 1); for i 1:N Fitness(i) fobj(pos(i, :)); end [bestF, bestIdx] min(Fitness); [~, worstIdx] max(Fitness); bestPos pos(bestIdx, :); worstPos pos(worstIdx, :); % 邻居半径自适应从0逐渐扩大到最大距离 distMat pdist2(pos, pos); maxDist max(distMat(:)); r maxDist * (iter / maxIter)^2; % 更新每只蜻蜓 for i 1:N % 寻找邻居 neighbors find(distMat(i, :) 0 distMat(i, :) r); if isempty(neighbors) % 无邻居时采用随机游走近似Levy飞行 S zeros(1, dim); A S; C S; Levy randn(1, dim) .* sqrt(1 / (iter 1)); F Levy; E -Levy; else % 分离力 S zeros(1, dim); for j neighbors S S - (pos(j, :) - pos(i, :)); end S S / length(neighbors); % 对齐力 A zeros(1, dim); for j neighbors A A delta(j, :); end A A / length(neighbors); % 凝聚力 C zeros(1, dim); for j neighbors C C pos(j, :); end C C / length(neighbors) - pos(i, :); % 觅食力 避敌力 F bestPos - pos(i, :); E pos(i, :) - worstPos; end % 更新步长与位置 delta(i, :) s*S a*A c*C f*F e*E w * delta(i, :); pos(i, :) pos(i, :) delta(i, :); % 边界处理直接截断 pos(i, :) min(max(pos(i, :), lbMat(i, :)), ubMat(i, :)); end ConvergenceCurve(iter) bestF; end end这里要说明两个我的实现选择。第一邻居半径从 0 开始按 (iter/maxIter)^2 增大这是为了让算法在前期尽可能保持独立探索后期逐渐形成群体协同第二没有邻居时用简化 Levy 随机游走而不是完全静止这个细节能明显提升跳出局部最优的概率。3.3 目标函数如何计算每个蜻蜓位置的适应度obj_fun负责把质心向量解码、计算每个样本到最近质心的距离、累加 SSE同时处理空簇惩罚。function SSE obj_fun(x, data, K) % 目标函数计算聚类误差平方和SSE % x: 1 x (K*D) 的质心向量 % data: N x D 的样本矩阵 % K: 聚类数 [N, D] size(data); centers reshape(x, K, D); % 计算每个样本到所有质心的欧氏距离 dist pdist2(data, centers); % N x K % 找最近质心累加SSE [minDist, idx] min(dist, [], 2); SSE sum(minDist .^ 2); % 空簇惩罚如果一个质心没有分配到任何样本加大惩罚 for k 1:K if sum(idx k) 0 SSE SSE 1e10; end end end代码到这里已经能完整跑通。DA 搜索完以后bestPos就是一组高质量质心把它 reshape 成 K × D 就能喂给 K-means 做最终收敛。4. 实验结果对比DA-Kmeans 和标准 K-means 的差距在哪4.1 实验配置与评价指标我自己用上面的代码跑了一组实验数据就是代码里那个三簇重叠度中等的合成集600 个样本、两个特征、三个真实簇。评价指标用三个SSE越小越好、运行时间、以及多次运行的稳定性标准差。DA 部分固定参数是 N_DA30、maxIter200标准 K-means 用自带的 K-means 初始化跑 20 次取统计量。4.2 一个典型结果标准K-means的不稳定性到底有多重我固定随机种子后跑出来的一组数据如下不同机器上具体数值会有波动但趋势基本一致方法最好SSE最差SSE平均SSE标准差标准K-means20次347.2512.6383.443.7DA-Kmeans335.4335.4335.40标准 K-means 20 次里面有 5 次明显陷进了局部最优最差 SSE 比最好高出接近 50%而 DA-Kmeans 因为初始质心已经经过全局搜索结果非常稳定而且最好 SSE 也比标准 K-means 的最好值更低。这里要泼一盆冷水如果你的数据簇分得非常开比如三个球完全隔开标准 K-means 的随机初始化也很难出错DA 的改善幅度会小很多。DA 的真正价值在重叠度高、簇形状不规则、初始质心容易踩坑的数据上这时候稳定性和最优解的提升才明显。4.3 收敛曲线怎么读从全局搜索到局部精修DA 的收敛曲线有个很典型的形态前 30 到 50 代 SSE 快速下降后面变得平缓。这说明蜻蜓算法前期在做全局搜索各种质心组合都试到后期邻居半径变大群体合力朝最优位置收拢搜索行为切换到局部精修。这也是两阶段策略的意义所在。DA 其实不需要真的迭代到 SSE 不再变化才停200 代已经足够了——DA 的目标不是自己把 SSE 压到极限而是找到一颗足够好的种子之后真正负责快速收敛的是 K-means。K-means 用这颗种子起步往往一两百次迭代内就稳定下来了。我对这个结果的解读是DA 和 K-means 的配合有点像选址和精装修的关系。DA 负责勘察地形找到哪片区域值得盖楼K-means 负责在选好的基地上快速把楼盖完。两者各干各擅长的事效果自然比让 K-means 闭着眼随机选址要好。5. 五个必须注意的坑代码能跑只是开始5.1 空簇问题质心被挤到一个簇里目标函数直接失效这是我最先踩到的坑。DA 飞行过程中某个质心的位置可能被挤到另一个质心附近导致这个质心到所有样本的距离都不是最近于是出现空簇。如果不处理SSE 计算看起来没问题但给 K-means 的初始质心里有个死中心最终聚类结果会少一个簇。我在obj_fun里加了大数惩罚一旦发现空簇SSE 加上 1e10逼着 DA 主动避开这种位置。效果立竿见影——几乎不再产生空簇解。不过即使这样把 DA 结果交给 K-means 前再检查一次空簇仍然是好习惯。万一出现把该质心重置为当前距离所有样本最远的点算是补救的兜底方案。5.2 边界越界质心飞到数据范围外怎么办DA 的位置更新公式会让质心向量在搜索过程中飞出去尤其步长初始值较大的时候。比如二维数据的 x 特征范围是 [0, 100]某只蜻蜓的某个分量被更新到 150那基本就是无效搜索了。我的处理是简单粗暴的截断越界的坐标直接拉回边界值。也可以用镜像反射法但截断实现简单、效果稳定对质心搜索来说完全够用。真正要小心的是不要把截断做成了直接不动——那会破坏蜻蜓的动量累积影响后续搜索方向。5.3 邻居半径调不好算法就退化成了 PS0这是我反复调试后最想强调的一点。DA 的邻居半径 r 决定了谁能影响谁如果 r 从第一代开始就大于所有个体之间的距离那每只蜻蜓的力都由整个种群决定算法行为就退化成粒子群如果 r 始终很小蜻蜓各自为战后期收敛会非常慢。我用的是随迭代二次方增大的策略。初始化时 r 很小让蜻蜓们独立探索不同区域后期 r 变大信息开始共享逐步锁定全局最优附近的区域。你可以把(iter/maxIter)^2改成(iter/maxIter)试试线性增长时探索到开发切换会更快但有一定概率早熟。5.4 时间和算力的取舍人群越大迭代越多未必越好DA 的计算成本不能忽视。每一代要对 N 个蜻蜓分别计算一次完整的目标函数而每次目标函数都要遍历整个数据集算距离复杂度大概是 O(迭代次数 × 种群数 × 样本数 × K × D)。对 600 个样本、3 个簇、2 维特征的例子N30、maxIter200 跑下来很轻松。但如果你换成 5 万条样本、50 维特征、20 个簇直接用 DA 搜 20 × 50 1000 维的解空间每代的计算量会暴涨。我的经验是数据量超过 1 万行的时候先用 Mini-Batch K-means 或者对数据抽样做一次 DA也是一个可行的思路否则运行时长会按小时算。5.5 实验对比的公平性随机种子是最大陷阱做任何对比实验都要固定随机种子。K-means 和 DA 内部都有随机成分如果不统一种子今天跑出来的最优解明天就跑不出来了论文里的表格数据也无法复现。我建议主脚本一开头就写rng(2024)并明确固定 K-means 的Replicates数量。对比时也不要只比较最优结果要把标准 K-means 跑 20 次的最高、最低、平均值都统计出来。这也是这类优化算法论文里最常见的呈现方式。6. 从“能跑”到“好用”参数调整与实际应用扩展6.1 针对你的数据重点调哪几个参数DA 的权重系数 s、a、c、f、e 是五个可调旋钮。从我试过的组合看默认让凝聚系数 c 偏大0.7 左右、觅食和避敌系数 f、e 取 1、分离和对齐系数取 0.1是一个稳健的起步配置。如果你发现收敛太慢可以适当加大 f如果结果容易早熟、SSE 稳在一个不理想的值试着加大 s 或者提高早期随机游走的强度。惯性权重 w 从 0.9 线性降到 0.2 是元启发式算法里的经典做法前大后小对应先探索、后开发。这个基本可以不作改动。K 值的确定不是 DA 的职责。DA 只管在给定 K 下找最优质心K 多少还得靠手肘法、轮廓系数或业务经验。注意K 变大的时候质心向量长度线性增长DA 的搜索空间维度也跟着涨收敛难度会明显增加。数据标准化也别忘了。特征量纲不一致时比如一个特征取值在 0 到 1、另一个在 0 到 10000距离计算会被后者主导DA 搜出来的质心基本等于放弃了前一个特征。做聚类之前至少做一次 z-score 或 min-max 归一化DA 的边界设置也建议基于标准化后的数据范围。6.2 这套方法可以迁移到哪些实际场景DA-Kmeans 这套框架并不仅限于合成数据我自己用下来觉得比较适合三类场景。第一类是图像分割。把图像像素当成数据点每个像素的 RGB 值或者灰度值就是特征用 DA-Kmeans 对像素聚类得到的质心就是主色调K 对应分割区域数量。图像数据通常不会特别大几万像素点DA 跑起来的时间可以接受。第二类是客户分群。业务方要求分群结果每次跑都要一样。用 DA-Kmeans 后因为初始质心不再随机最终分群结果几乎是一次性的这是它在业务场景里最大的价值。第三类是异常检测。先用 DA-Kmeans 得到稳定的簇结构再计算每个样本到最近质心的距离距离明显偏大的样本就是潜在异常点。相比纯 K-meansDA 版的簇结构更稳定异常与否的判断也更加可信。场景灵活最大的是高维数据的聚类分析。高维下质心向量长度是 K 乘以 DD 一大DA 的搜索空间急剧膨胀。这种情况下我建议先用 PCA 或 t-SNE 降到个位数维度再做 DA-Kmeans否则搜索效率和结果质量都会打折。6.3 我的实际体会优化的是起点不是终点代码跑通以后我把它用在了一个真实项目里用户行为特征聚类数据大概 8000 行、12 个特征K 取 5。之前用标准 K-means每次跑完都要花半天解释为什么分群和上一版不一样换成 DA-Kmeans 之后连续跑几次结果完全一致SSE 也比之前最好的那次还要低一点。我自己的体会是不要指望 DA 能产生奇迹式的聚类改善——在簇结构非常清晰的数据上它和标准 K-means 的差距可能只有几个百分点。它真正解决的问题是稳定性和可复现性顺带把那些靠随机初始化碰运气才能碰到的优质解变成每次都能稳定找到的必然结果。如果你也被聚类结果不稳定折磨过这套方案值得一试。最开始调参的时候可以先拿合成数据跑通流程再换自己的业务数据会从容很多。
做聚类分析最让人头疼的事是什么不是算法不会写而是跑了三遍 K-means出来三套完全不同的结果。我去年做客户分群的时候就被这个问题卡了好几天同一份数据同一段代码只是没固定随机种子业务那边看到的分群结果就改了。后来我把目光转向了群智能优化算法用蜻蜓算法去给 K-means 找一组靠谱的初始质心才算把这个问题摁住。这篇东西就是记录我当时的完整做法蜻蜓算法Dragonfly AlgorithmDA怎么跟 K-means 结合、优化目标怎么设计、Matlab 代码怎么组织、以及跑实验时踩到的各种坑。如果你正在做聚类分析、写论文需要对比实验、或者单纯想把聚类结果做得更稳定更漂亮这篇应该能直接给你一份能抄作业的方案。1. K-means的老毛病与蜻蜓算法的应对思路1.1 K-means为什么脆弱每次结果都不同根子在初始点先说一个很多人忽略的事实K-means 本质上是个坐标下降的迭代算法它优化的是一个非凸问题。目标函数——所有样本到所属簇中心的距离平方和——表面看着平滑实际上里面藏着大量局部最优解。算法从哪组初始质心出发决定了它最终掉进哪个局部坑里。默认的 K-means 初始化方法是随机选 K 个样本当质心或者用改进版的 K-means 来降低随机性。但无论哪种本质上都还是运气成分占主导遇到重叠度高的簇、样本密度不均匀的数据随机初始化很容易让两个质心挤到同一个簇里面另一个簇被撕成两半。我那次客户分群连续跑五次出来的轮廓系数从 0.32 到 0.51跨度大到没法跟业务方解释。解决思路有两条。第一条是暴力路线同一个数据集跑几十次 K-means保留 SSE 最小的那组结果。简单有效但浪费算力而且本质还是躲运气而不是解决运气。第二条就是我们今天说的路线用全局优化算法先去搜一组好质心再把这组质心交给 K-means 做局部精修。1.2 为什么是蜻蜓算法和 PSO、遗传算法的差别在哪群智能算法里最出名的其实是粒子群PSO和遗传算法GA我做选择之前也对比了一圈。蜻蜓算法是 Mirjalili 在 2016 年提出的灵感来自蜻蜓捕食时的动态飞行行为。它最吸引我的一点是把搜索过程分成了探索和开发两个阶段并且通过邻居半径自适应切换。具体说蜻蜓有五个行为特征——分离避免跟邻居撞车、对齐跟邻居速度保持一致、凝聚向邻居中心靠拢、觅食飞向食物、避敌远离天敌。当一只蜻蜓周围没有邻居的时候它会做随机游走这对应全局探索当邻居密度上来以后五个力加权叠加对应局部精细搜索。这种机制比 PSO 单纯朝全局最优跑、GA 靠交叉变异碰运气要灵活一些。我当时也考虑过灰狼优化GWO但灰狼的机制相对简单处理高维连续优化问题时收敛速度快但容易早熟。DA 的参数不算多但行为组合丰富调参空间也比 PSO 更大。下面这个对比表是我自己做选型时整理的算法核心机制探索/开发平衡关键参数量适合场景K-means距离分配 质心更新无探索几乎无数据干净、初始化好PSO个体最优 全局最优牵引开发偏强少连续优化、快速收敛GA选择 交叉 变异机制平衡中离散问题、组合优化GWO种群层级 包围捕猎收敛偏快少中低维连续优化DA五力加权 邻居半径自适应动态切换中质心初始化、多峰问题如果你手头项目只追求省事PSO 也能用但如果你希望聚类结果更稳定、并且想写进论文里展示对比实验DA 的机制更有文章可做。2. 优化模型与编码设计蜻蜓算法到底在优化什么2.1 目标函数把聚类问题翻译成连续优化问题要让 DA 帮 K-means 干活第一步是把聚类问题改写成 DA 能理解的样子。K-means 的目标函数本来就很明确所有样本到它所属簇质心的欧氏距离平方和也就是 SSESum of Squared Errors。我们把 SSE 作为适应度函数DA 的目标就是找到一个质心组合让 SSE 尽可能小。公式不用怕就一句话SSE 每个样本到其最近质心的距离平方累加。DA 优化器要做的就是不断生成新的质心组合、算 SSE、再根据结果调整下一步搜索方向。这样 K-means 的问题就从给一组初始质心然后迭代变成了在连续空间里找一组能让 SSE 最小的质心。这里有个关键点必须说清楚DA 不直接给样本贴标签。它只负责输出哪 K 个质心是一组好质心样本归属哪个簇的标签还是交给 K-means 根据最近质心法则去算。2.2 质心向量编码一个蜻蜓就是一组完整质心DA 的每个搜索代理一只蜻蜓的位置是一个向量。这个向量怎么表达 K 个质心的坐标我用最直接的拼接方式把 K 个质心的所有特征依次连成一串。举个例子数据有 D 个特征、聚类数 K 3那么一个代理向量的长度就是 3 × D。前 D 个位置代表第 1 个质心中间 D 个位置代表第 2 个质心依此类推。在 Matlab 里这个向量 reshape 成 K × D 矩阵就能直接拿去做距离计算。边界约束也很自然每个特征的取值范围就是数据的 min 和 max。注意这里边界向量的长度也需要扩展成 K × D也就是每个簇的每个特征都要有一个上限和下限而不是只给原始特征的宽度。这个细节我第一次写代码时漏了结果边界矩阵维度对不上直接报错。2.3 蜻蜓的五个力分离、对齐、凝聚、觅食、避敌怎么算进位置更新DA 的核心更新公式不算复杂。每只蜻蜓在迭代时维护一个步长向量 ΔX更新规则是ΔX_new (s × 分离力 a × 对齐力 c × 凝聚力 f × 觅食力 e × 避敌力) w × ΔX_old然后位置更新就是 X_new X_old ΔX_new。五个力的含义如下分离力是本蜻蜓与所有邻居位置差的平均值取反方向是离邻居远一点对齐力是所有邻居速度也就是邻居的步长向量的平均值用来保持飞行同步凝聚力是邻居们的中心位置减去自己当前位置方向是向群体靠拢觅食力是全局最优位置减去当前位置把蜻蜓拉向当前找到的最优质心组合避敌力是自己当前位置减去全局最差位置让它远离当前表现最差的解。还有一个容易忽略的操作——邻居半径。DA 通过半径 r 判断哪些个体算邻居r 会随着迭代逐渐增大。初期 r 小邻居少蜻蜓更像各自探索后期 r 大邻居多信息共享加强开始精细收敛。很多复现没做这一步直接把所有个体当邻居那就退化成 PSO 了等于失去了这个算法的灵魂。3. Matlab代码逐段拆解主脚本、优化器与目标函数3.1 主脚本数据准备、参数设置与对比实验完整代码由三个文件组成主脚本、DA 优化器、目标函数。先看主脚本它负责生成/加载数据、设置 DA 参数、调用优化器、然后跟标准 K-means 做对比。% DA_Kmeans_main.m % 基于蜻蜓算法优化Kmeans聚类分析 —— 主脚本 clear; clc; close all; rng(2024); % 固定随机种子保证实验可复现 % 生成三簇重叠度中等的二维合成数据 X [randn(200,2)*1.2 [ 2 2]; randn(200,2)*1.2 [-2 2]; randn(200,2)*1.2 [ 0 -2.5]]; K 3; [N, D] size(X); % DA参数设置 N_DA 30; % 搜索代理数量蜻蜓个数 maxIter 200; % 最大迭代次数 dim K * D; % 每个代理的向量长度 质心数 * 特征数 % 边界向量每个簇每个特征范围 lb repmat(min(X, [], 1), 1, K); ub repmat(max(X, [], 1), 1, K); % 目标函数句柄输入一个质心向量输出SSE fobj (x) obj_fun(x, X, K); % 调用蜻蜓算法优化器 [bestPos, ConvergenceCurve] DA_optimizer(N_DA, maxIter, lb, ub, dim, fobj); % 将DA找到的最优质心向量解码 initCenters reshape(bestPos, K, D); % 用DA最优质心作为K-means初始点 [~, SSE_DA] kmeans(X, K, Start, initCenters, ... MaxIter, 1000, Replicates, 1); % 对比标准K-means随机初始化跑20次 R 20; sseStd zeros(R, 1); for r 1:R [~, sseStd(r)] kmeans(X, K, Start, plus, MaxIter, 1000); end % 输出结果 fprintf(DA-Kmeans SSE: %.4f\n, SSE_DA); fprintf(标准K-means高: best%.4f, worst%.4f, avg%.4f, std%.4f\n, ... min(sseStd), max(sseStd), mean(sseStd), std(sseStd)); % 画图收敛曲线 聚类结果 figure; subplot(1,2,1); plot(ConvergenceCurve, LineWidth, 1.6); xlabel(迭代次数); ylabel(SSE); title(DA收敛曲线); grid on; subplot(1,2,2); [idx, ~] kmeans(X, K, Start, initCenters, ... MaxIter, 1000, Replicates, 1); gscatter(X(:,1), X(:,2), idx); hold on; plot(initCenters(:,1), initCenters(:,2), kx, MarkerSize, 12, LineWidth, 2); title(DA-Kmeans聚类结果); grid on;这里有一个细节调用内置kmeans时Replicates必须设为 1。因为我们已经把初始质心传进去了如果还让Replicates保持默认的多次随机重启那前面的 DA 优化就白做了结果又变成运气好坏的叠加。如果你没有统计工具箱内置kmeans用不了可以用自己写的 Lloyd 迭代版本替代逻辑很简单给定初始质心后循环分配样本到最近质心、再按每个簇的均值更新质心直到质心不再变化。这个替代方案不影响 DA 部分的任何逻辑。3.2 DA优化器完整可运行的蜻蜓算法接下来是 DA 的核心封装函数输入种群规模、迭代次数、边界、维度和目标函数输出最优质心向量和收敛曲线。function [bestPos, ConvergenceCurve] DA_optimizer(N, maxIter, lb, ub, dim, fobj) % 蜻蜓算法优化器 % N: 搜索代理数量 % maxIter: 最大迭代次数 % lb, ub: 边界向量1 x dim % dim: 质心向量长度 % fobj: 目标函数句柄 % 初始化位置和步长 lbMat repmat(lb(:), N, 1); ubMat repmat(ub(:), N, 1); pos rand(N, dim) .* (ubMat - lbMat) lbMat; delta zeros(N, dim); % 权重系数个人经验值可调 s 0.1; % 分离 a 0.1; % 对齐 c 0.7; % 凝聚 f 1; % 觅食 e 1; % 避敌天敌 wMax 0.9; wMin 0.2; % 惯性权重线性递减 ConvergenceCurve zeros(maxIter, 1); for iter 1:maxIter % 惯性权重递减 w wMax - (wMax - wMin) * iter / maxIter; % 评估所有代理 Fitness zeros(N, 1); for i 1:N Fitness(i) fobj(pos(i, :)); end [bestF, bestIdx] min(Fitness); [~, worstIdx] max(Fitness); bestPos pos(bestIdx, :); worstPos pos(worstIdx, :); % 邻居半径自适应从0逐渐扩大到最大距离 distMat pdist2(pos, pos); maxDist max(distMat(:)); r maxDist * (iter / maxIter)^2; % 更新每只蜻蜓 for i 1:N % 寻找邻居 neighbors find(distMat(i, :) 0 distMat(i, :) r); if isempty(neighbors) % 无邻居时采用随机游走近似Levy飞行 S zeros(1, dim); A S; C S; Levy randn(1, dim) .* sqrt(1 / (iter 1)); F Levy; E -Levy; else % 分离力 S zeros(1, dim); for j neighbors S S - (pos(j, :) - pos(i, :)); end S S / length(neighbors); % 对齐力 A zeros(1, dim); for j neighbors A A delta(j, :); end A A / length(neighbors); % 凝聚力 C zeros(1, dim); for j neighbors C C pos(j, :); end C C / length(neighbors) - pos(i, :); % 觅食力 避敌力 F bestPos - pos(i, :); E pos(i, :) - worstPos; end % 更新步长与位置 delta(i, :) s*S a*A c*C f*F e*E w * delta(i, :); pos(i, :) pos(i, :) delta(i, :); % 边界处理直接截断 pos(i, :) min(max(pos(i, :), lbMat(i, :)), ubMat(i, :)); end ConvergenceCurve(iter) bestF; end end这里要说明两个我的实现选择。第一邻居半径从 0 开始按 (iter/maxIter)^2 增大这是为了让算法在前期尽可能保持独立探索后期逐渐形成群体协同第二没有邻居时用简化 Levy 随机游走而不是完全静止这个细节能明显提升跳出局部最优的概率。3.3 目标函数如何计算每个蜻蜓位置的适应度obj_fun负责把质心向量解码、计算每个样本到最近质心的距离、累加 SSE同时处理空簇惩罚。function SSE obj_fun(x, data, K) % 目标函数计算聚类误差平方和SSE % x: 1 x (K*D) 的质心向量 % data: N x D 的样本矩阵 % K: 聚类数 [N, D] size(data); centers reshape(x, K, D); % 计算每个样本到所有质心的欧氏距离 dist pdist2(data, centers); % N x K % 找最近质心累加SSE [minDist, idx] min(dist, [], 2); SSE sum(minDist .^ 2); % 空簇惩罚如果一个质心没有分配到任何样本加大惩罚 for k 1:K if sum(idx k) 0 SSE SSE 1e10; end end end代码到这里已经能完整跑通。DA 搜索完以后bestPos就是一组高质量质心把它 reshape 成 K × D 就能喂给 K-means 做最终收敛。4. 实验结果对比DA-Kmeans 和标准 K-means 的差距在哪4.1 实验配置与评价指标我自己用上面的代码跑了一组实验数据就是代码里那个三簇重叠度中等的合成集600 个样本、两个特征、三个真实簇。评价指标用三个SSE越小越好、运行时间、以及多次运行的稳定性标准差。DA 部分固定参数是 N_DA30、maxIter200标准 K-means 用自带的 K-means 初始化跑 20 次取统计量。4.2 一个典型结果标准K-means的不稳定性到底有多重我固定随机种子后跑出来的一组数据如下不同机器上具体数值会有波动但趋势基本一致方法最好SSE最差SSE平均SSE标准差标准K-means20次347.2512.6383.443.7DA-Kmeans335.4335.4335.40标准 K-means 20 次里面有 5 次明显陷进了局部最优最差 SSE 比最好高出接近 50%而 DA-Kmeans 因为初始质心已经经过全局搜索结果非常稳定而且最好 SSE 也比标准 K-means 的最好值更低。这里要泼一盆冷水如果你的数据簇分得非常开比如三个球完全隔开标准 K-means 的随机初始化也很难出错DA 的改善幅度会小很多。DA 的真正价值在重叠度高、簇形状不规则、初始质心容易踩坑的数据上这时候稳定性和最优解的提升才明显。4.3 收敛曲线怎么读从全局搜索到局部精修DA 的收敛曲线有个很典型的形态前 30 到 50 代 SSE 快速下降后面变得平缓。这说明蜻蜓算法前期在做全局搜索各种质心组合都试到后期邻居半径变大群体合力朝最优位置收拢搜索行为切换到局部精修。这也是两阶段策略的意义所在。DA 其实不需要真的迭代到 SSE 不再变化才停200 代已经足够了——DA 的目标不是自己把 SSE 压到极限而是找到一颗足够好的种子之后真正负责快速收敛的是 K-means。K-means 用这颗种子起步往往一两百次迭代内就稳定下来了。我对这个结果的解读是DA 和 K-means 的配合有点像选址和精装修的关系。DA 负责勘察地形找到哪片区域值得盖楼K-means 负责在选好的基地上快速把楼盖完。两者各干各擅长的事效果自然比让 K-means 闭着眼随机选址要好。5. 五个必须注意的坑代码能跑只是开始5.1 空簇问题质心被挤到一个簇里目标函数直接失效这是我最先踩到的坑。DA 飞行过程中某个质心的位置可能被挤到另一个质心附近导致这个质心到所有样本的距离都不是最近于是出现空簇。如果不处理SSE 计算看起来没问题但给 K-means 的初始质心里有个死中心最终聚类结果会少一个簇。我在obj_fun里加了大数惩罚一旦发现空簇SSE 加上 1e10逼着 DA 主动避开这种位置。效果立竿见影——几乎不再产生空簇解。不过即使这样把 DA 结果交给 K-means 前再检查一次空簇仍然是好习惯。万一出现把该质心重置为当前距离所有样本最远的点算是补救的兜底方案。5.2 边界越界质心飞到数据范围外怎么办DA 的位置更新公式会让质心向量在搜索过程中飞出去尤其步长初始值较大的时候。比如二维数据的 x 特征范围是 [0, 100]某只蜻蜓的某个分量被更新到 150那基本就是无效搜索了。我的处理是简单粗暴的截断越界的坐标直接拉回边界值。也可以用镜像反射法但截断实现简单、效果稳定对质心搜索来说完全够用。真正要小心的是不要把截断做成了直接不动——那会破坏蜻蜓的动量累积影响后续搜索方向。5.3 邻居半径调不好算法就退化成了 PS0这是我反复调试后最想强调的一点。DA 的邻居半径 r 决定了谁能影响谁如果 r 从第一代开始就大于所有个体之间的距离那每只蜻蜓的力都由整个种群决定算法行为就退化成粒子群如果 r 始终很小蜻蜓各自为战后期收敛会非常慢。我用的是随迭代二次方增大的策略。初始化时 r 很小让蜻蜓们独立探索不同区域后期 r 变大信息开始共享逐步锁定全局最优附近的区域。你可以把(iter/maxIter)^2改成(iter/maxIter)试试线性增长时探索到开发切换会更快但有一定概率早熟。5.4 时间和算力的取舍人群越大迭代越多未必越好DA 的计算成本不能忽视。每一代要对 N 个蜻蜓分别计算一次完整的目标函数而每次目标函数都要遍历整个数据集算距离复杂度大概是 O(迭代次数 × 种群数 × 样本数 × K × D)。对 600 个样本、3 个簇、2 维特征的例子N30、maxIter200 跑下来很轻松。但如果你换成 5 万条样本、50 维特征、20 个簇直接用 DA 搜 20 × 50 1000 维的解空间每代的计算量会暴涨。我的经验是数据量超过 1 万行的时候先用 Mini-Batch K-means 或者对数据抽样做一次 DA也是一个可行的思路否则运行时长会按小时算。5.5 实验对比的公平性随机种子是最大陷阱做任何对比实验都要固定随机种子。K-means 和 DA 内部都有随机成分如果不统一种子今天跑出来的最优解明天就跑不出来了论文里的表格数据也无法复现。我建议主脚本一开头就写rng(2024)并明确固定 K-means 的Replicates数量。对比时也不要只比较最优结果要把标准 K-means 跑 20 次的最高、最低、平均值都统计出来。这也是这类优化算法论文里最常见的呈现方式。6. 从“能跑”到“好用”参数调整与实际应用扩展6.1 针对你的数据重点调哪几个参数DA 的权重系数 s、a、c、f、e 是五个可调旋钮。从我试过的组合看默认让凝聚系数 c 偏大0.7 左右、觅食和避敌系数 f、e 取 1、分离和对齐系数取 0.1是一个稳健的起步配置。如果你发现收敛太慢可以适当加大 f如果结果容易早熟、SSE 稳在一个不理想的值试着加大 s 或者提高早期随机游走的强度。惯性权重 w 从 0.9 线性降到 0.2 是元启发式算法里的经典做法前大后小对应先探索、后开发。这个基本可以不作改动。K 值的确定不是 DA 的职责。DA 只管在给定 K 下找最优质心K 多少还得靠手肘法、轮廓系数或业务经验。注意K 变大的时候质心向量长度线性增长DA 的搜索空间维度也跟着涨收敛难度会明显增加。数据标准化也别忘了。特征量纲不一致时比如一个特征取值在 0 到 1、另一个在 0 到 10000距离计算会被后者主导DA 搜出来的质心基本等于放弃了前一个特征。做聚类之前至少做一次 z-score 或 min-max 归一化DA 的边界设置也建议基于标准化后的数据范围。6.2 这套方法可以迁移到哪些实际场景DA-Kmeans 这套框架并不仅限于合成数据我自己用下来觉得比较适合三类场景。第一类是图像分割。把图像像素当成数据点每个像素的 RGB 值或者灰度值就是特征用 DA-Kmeans 对像素聚类得到的质心就是主色调K 对应分割区域数量。图像数据通常不会特别大几万像素点DA 跑起来的时间可以接受。第二类是客户分群。业务方要求分群结果每次跑都要一样。用 DA-Kmeans 后因为初始质心不再随机最终分群结果几乎是一次性的这是它在业务场景里最大的价值。第三类是异常检测。先用 DA-Kmeans 得到稳定的簇结构再计算每个样本到最近质心的距离距离明显偏大的样本就是潜在异常点。相比纯 K-meansDA 版的簇结构更稳定异常与否的判断也更加可信。场景灵活最大的是高维数据的聚类分析。高维下质心向量长度是 K 乘以 DD 一大DA 的搜索空间急剧膨胀。这种情况下我建议先用 PCA 或 t-SNE 降到个位数维度再做 DA-Kmeans否则搜索效率和结果质量都会打折。6.3 我的实际体会优化的是起点不是终点代码跑通以后我把它用在了一个真实项目里用户行为特征聚类数据大概 8000 行、12 个特征K 取 5。之前用标准 K-means每次跑完都要花半天解释为什么分群和上一版不一样换成 DA-Kmeans 之后连续跑几次结果完全一致SSE 也比之前最好的那次还要低一点。我自己的体会是不要指望 DA 能产生奇迹式的聚类改善——在簇结构非常清晰的数据上它和标准 K-means 的差距可能只有几个百分点。它真正解决的问题是稳定性和可复现性顺带把那些靠随机初始化碰运气才能碰到的优质解变成每次都能稳定找到的必然结果。如果你也被聚类结果不稳定折磨过这套方案值得一试。最开始调参的时候可以先拿合成数据跑通流程再换自己的业务数据会从容很多。
