基于蜻蜓算法优化K-means初始中心的Matlab实现与聚类分析

提起聚类分析很多刚入门的朋友第一反应就是K-means也有一部分人用过SPSS里的聚类模块点几下鼠标就能出结果。但真正做过聚类项目的人都知道K-means这玩意儿最让人头疼的不是算法本身而是稳定性和初始中心的选择。同一个数据集你多跑几次结果经常长得完全不一样。我今天要分享的就是一套我自己在项目里反复打磨过的方案用蜻蜓算法Dragonfly Algorithm去优化K-means的初始中心再在Matlab里把整个流程完整实现出来。这套代码我封装成了函数可以直接拿去换数据、换聚类数也能替换成别的群智能优化算法做横向对比适合正在做聚类分析、想深入了解群智能算法和经典算法怎么配合使用的读者。1. 聚类分析的痛点与整体优化思路1.1 K-means为什么需要被“优化”K-means的流程教科书里写得很清楚随机挑选K个样本作为初始质心然后把每个样本分配到距离最近的质心重新计算质心位置重复直到质心不再变化。它的优点是快、简单、内存占用少海量数据也能跑。可它的缺点跟优点一样突出。第一个问题是初始值敏感。K-means本质上是一个坐标下降式的局部搜索过程搜索结果受初始质心影响极大。初始质心选在什么位置最后多半就收敛到附近的局部最优。这就像你从不同起点同时往山下走最后停在哪个山谷完全取决于起点在哪而不是哪个山谷最深。第二个问题是结果不稳定。正因为初始质心是随机选的同一份数据、同一个K值今天跑和明天跑10次运行能出好几种聚类形态。做工程项目最怕这种不确定性你给客户演示的时候跑出一个好结果换一台机器跑又变成另一个样子这没法交代。第三个问题是对数据分布敏感。K-means默认簇是“球形且各向同性”的等方差、等样本量时才表现最好。现实数据经常是长条形、重叠、带离群点的K-means在这种数据上容易把簇切得乱七八糟。所以优化思路其实很明确既然K-means缺的是全局搜索能力那就引入一个具备全局搜索能力的算法先帮它找一组好的初始中心再让它做局部精调。1.2 为什么我选蜻蜓算法改进K-means初始化的思路并不新鲜。早些年大家爱用遗传算法后来粒子群PSO火了一阵再后来是灰狼优化GWO、鲸鱼优化WOA这些新群智能算法。我也都试过现在项目里固定用蜻蜓算法主要看中三点。第一DA的探索与开发平衡机制更完整。DA借鉴了蜻蜓群体的两种行为静态群小范围飞翔、捕食对应算法前期的全局探索迁移群大规模迁徙、避敌对应算法后期的局部开发。它通过邻域半径的收缩自动从探索过渡到开发这个机制比PSO单纯靠惯性权重衰减要自然。第二Levy飞行自带“跳出局部最优”的能力。DA在没有邻居时会执行Levy飞行这是一种重尾随机游走大部分时间小步走偶尔来一大步。这个偶尔的大步对优化算法来说非常宝贵相当于给搜索过程加了随机变异不容易卡死在局部最优。第三代码量可控。DA名义上有五个行为因子实现起来就是五个向量的线性组合加一个Levy飞行分支在Matlab里封装成函数核心代码100多行调试维护都很容易。1.3 “全局搜索局部精调”两阶段框架我的整体框架是两个阶段的串联。第一阶段让DA全局搜索K个聚类中心的最优位置第二阶段把DA找到的最优中心送给K-means做局部精调进一步加速收敛并得到最终的标签分配。这个设计有一个容易被忽略的好处DA负责找到“好的盆地”K-means负责在盆地内“快速下坡”。前者弥补K-means的全局性短板后者弥补DA收敛慢、局部精调能力弱的短板两个算法互补而不是互相替代。工程上还有个好处这个框架天然支持“算法替换”。如果哪天你想试试GWO或者WOA只需要把优化器换成对应算法聚类适应度函数完全不用动。我在代码结构里特意把“优化器”和“聚类适应度函数”拆成两个独立函数下面你会看到具体写法。2. 蜻蜓算法核心机制与公式拆解2.1 五个行为因子的物理含义蜻蜓算法是Mirjalili在2016年提出的核心思想来自蜻蜓聚群捕食和迁徙时的飞行行为。每个蜻蜓个体在空间中的位置就是一个候选解它如何移动取决于周边伙伴的行为和外部环境。DA更新位置时考虑了五个因子分离Separation避免个体间碰撞。如果当前个体离邻居太近就产生一个排斥力S_i -∑(X - X_j)对齐Alignment个体速度与邻居平均速度保持一致A_i (∑V_j) / N聚集Cohesion个体向邻居质心靠近C_i (∑X_j) / N - X觅食Food向食物源当前全局最优解靠拢F_i X_food - X避敌Enemy远离敌人当前全局最差解E_i X_enemy X这五个因子最后线性加权组合再叠加上一代的速度得到当前个体的速度更新量然后更新位置。2.2 位置更新与Levy飞行位置更新分两种情况。当个体邻域半径内存在其他蜻蜓时用标准的增量式更新dX(t1) (sS aA cC fF eE) wdX(t)X(t1) X(t) dX(t1)当邻域内没有其他蜻蜓时个体执行Levy飞行随机游走X(t1) X(t) levy(D) * X(t)Levy步长用Mantegna算法生成beta通常取1.5。代码实现下一节给你看这里先记住一个关键点Levy生成的步长范围很大直接乘X(i,:)容易越界我实际测试下来会在前面加一个0.01的缩放因子否则收敛曲线完全是乱的。另一个细节是邻域半径r怎么取。原论文里半径随迭代单调递增我的实现中反过来把它从大往小缩前期货架大一些、邻居多、探索充分后期半径缩小个体逐渐进入局部精细搜索。两种做法都有文献支持关键是跟你自己的问题规模匹配具体参数后面细说。2.3 参数设计与动态调整策略DA里有六个权重分离权重s、对齐权重a、聚集权重c、觅食权重f、避敌权重e、惯性权重w。我的默认策略都是线性递减w从0.9减到0.4s从0.1减到0.01a从0.1减到0.01c从0.7减到0.2f从1减到0.1e从1减到0.01这个策略的逻辑是前期让分离、对齐、聚集、觅食、避敌都保持较大的作用强度让蜻蜓群体在解空间里铺开充分探索后期这些权重变小群体收敛到当前最优附近进行精细搜索。为什么线性递减最常用因为实现简单且足够稳。对于大多数聚类问题线性递减曲线能给出可预期的收敛行为方便你先跑通再根据收敛曲线做手动调整。等你有经验了也可以换成非线性衰减或者自适应权重但没必要一上来就搞复杂。3. 把K-means“翻译”成优化问题3.1 质心编码一个蜻蜓个体等于一组聚类中心让DA优化K-means的第一步是找到两者之间的“连接语言”。DA处理的是一个实数向量K-means的决策变量是K个聚类中心所以最自然的编码就是“质心拼接”。设数据维度为D聚类数为K则一个蜻蜓个体的位置向量维度为K×DX [c1(1), c1(2), ..., c1(D), c2(1), ..., cK(D)]举个例子。二维数据、K3那么一个个体就是1×6的向量前三列是第一个聚类中心的x、y坐标第四到六列是第二个中心的坐标依此类推。编码之后个体的每一维取值范围要限制在数据集对应维度的最小值和最大值之间。这一步非常关键如果不加边界蜻蜓可能飞到离数据十万八千里的地方适应度函数算出来是一个非常大的数虽然算法也能“学”回来但会浪费掉大量迭代次数。3.2 适应度函数与聚类标签的关系适应度函数是优化器和聚类问题的接口。聚类要最小化的目标我选的是簇内误差平方和WCSSWCSS ∑_i min_j ||x_i - c_j||²翻译成人话对每个样本找到离它最近的质心计算距离的平方全部加起来。WCSS越小说明样本离各自质心越近簇内越紧凑。这个函数还附带了一个巨大的好处计算过程中得到的最近质心索引正好就是每个样本的聚类标签。所以适应度函数可以直接返回两份东西一份是适应度值一份是标签矩阵后面画聚类图的时候都不用重新计算一遍。3.3 为什么适应度选WCSS而不是轮廓系数你可能想问评价聚类效果的指标那么多为什么偏偏选WCSS做适应度。因为DA在迭代中要对每一个个体都算一次“聚类质量”。假设种群40个个体、迭代200次那就是8000次适应度计算。如果每次用轮廓系数轮廓系数的计算需要对所有样本两两求距离复杂度是O(N²)数据量大时这个代价完全不可承受。而WCSS只需要算样本到K个质心的距离复杂度是O(N×K×D)快几个数量级。我的建议是DA内部优化用WCSS加速最终聚类效果评估用轮廓系数把关。优化快评估准两不耽误。4. Matlab完整代码实现与逐段讲解4.1 主程序数据准备与总流程下面是我实际在用的主程序先生成一组三分类高斯数据然后调用DA_Kmeans函数再用Matlab自带的kmeans函数做精调。代码在R2020a及之后版本均可运行只需要统计工具箱。%% 主程序基于蜻蜓算法优化K-means聚类分析 clear; clc; close all; %% 1. 生成测试数据三簇高斯分布 rng(42); data1 mvnrnd([2, 2], [0.5, 0.1; 0.1, 0.4], 100); data2 mvnrnd([7, 7], [0.6, 0.05; 0.05, 0.5], 90); data3 mvnrnd([3, 8], [0.4, -0.05; -0.05, 0.6], 80); data [data1; data2; data3]; [N, D] size(data); K 3; %% 2. 数据归一化重要 for d 1:D data(:, d) (data(:, d) - min(data(:, d))) / ... (max(data(:, d)) - min(data(:, d))); end %% 3. DA优化K-means中心 [bestPos, bestWCSS, convergence] DA_Kmeans(data, K, ... maxIter, 150, nPop, 40); %% 4. 把DA最优解作为K-means初始中心做局部精调 bestCenters reshape(bestPos, K, D); [labels, centers] kmeans(data, K, Start, bestCenters, ... MaxIter, 500, Replicates, 1);data在第二步里被覆盖为归一化后的版本后面DA和kmeans用的都是同一份归一化数据距离度量语义一致。如果你把原始数据丢给DA、归一化数据丢给kmeans出来的质心位置会完全对不上这点后面讲踩坑时会再提。4.2 DA_Kmeans核心函数DA_Kmeans.m是核心函数。输入是数据矩阵、聚类数K和可选参数输出是DA找到的最优质心拼接向量、最优WCSS和每次迭代的收敛曲线。function [bestPos, bestFitness, curve] DA_Kmeans(data, K, varargin) % DA_Kmeans 蜻蜓算法优化K-means聚类中心 % 输入: % data - N×D样本矩阵 % K - 聚类数 % 可选参数: % maxIter 最大迭代次数(默认150) % nPop 种群大小(默认40) % 输出: % bestPos 最优个体(1×(K*D)的质心拼接向量) % bestFitness 最优适应度(WCSS) % curve 收敛曲线(1×maxIter) p inputParser; addParameter(p, maxIter, 150); addParameter(p, nPop, 40); parse(p, varargin{:}); maxIter p.Results.maxIter; nPop p.Results.nPop; [N, D] size(data); dim K * D; % 每一维的边界按数据范围 lbVec min(data, [], 1); ubVec max(data, [], 1); lb repmat(lbVec, 1, K); ub repmat(ubVec, 1, K); % 初始化种群 X rand(nPop, dim) .* (ub - lb) lb; dX zeros(nPop, dim); fitness zeros(1, nPop); for i 1:nPop fitness(i) kmeansFitness(X(i, :), data, K); end [bestFitness, bestIdx] min(fitness); bestPos X(bestIdx, :); foodPos bestPos; worstFitness max(fitness); enemyPos X(fitness worstFitness, :); enemyPos enemyPos(1, :); % 权重初始值 s 0.1; a 0.1; c 0.7; f 1; e 1; w 0.9; curve zeros(1, maxIter); for t 1:maxIter % 权重线性递减 w 0.9 - t * (0.9 - 0.4) / maxIter; s 0.1 - t * (0.1 - 0.01) / maxIter; a 0.1 - t * (0.1 - 0.01) / maxIter; c 0.7 - t * (0.7 - 0.2) / maxIter; f 1.0 - t * (1.0 - 0.1) / maxIter; e 1.0 - t * (1.0 - 0.01) / maxIter; % 邻域半径逐渐缩小促进探索向开发过渡 radius (ub(1) - lb(1)) * (1 - t / maxIter); for i 1:nPop % 计算邻域内其他个体 distPop sqrt(sum((X - X(i, :)).^2, 2)); neighbors find(distPop 0 distPop radius); if ~isempty(neighbors) % 分离项 S sum(X(i, :) - X(neighbors, :), 1); % 对齐项dX近似速度 A sum(dX(neighbors, :), 1) / numel(neighbors); % 聚集项 C mean(X(neighbors, :), 1) - X(i, :); % 觅食项 F foodPos - X(i, :); % 避敌项 E enemyPos X(i, :); dX(i, :) (s*S a*A c*C f*F e*E) w * dX(i, :); else % 无邻居Levy飞行全局随机游走 beta 1.5; sigma (gamma(1 beta) * sin(pi * beta / 2) / ... (gamma((1 beta)/2) * beta * 2^((beta - 1)/2)))^(1/beta); u randn(1, dim) * sigma; v randn(1, dim); step 0.01 * (u ./ (abs(v).^(1 / beta))); dX(i, :) step .* X(i, :); end % 更新位置并限制边界 X(i, :) X(i, :) dX(i, :); X(i, :) max(min(X(i, :), ub), lb); % 重新评估适应度 fitness(i) kmeansFitness(X(i, :), data, K); end % 更新全局最优与最差食物源/敌人 [curBest, bestIdx] min(fitness); [curWorst, worstIdx] max(fitness); if curBest bestFitness bestFitness curBest; bestPos X(bestIdx, :); end foodPos X(bestIdx, :); enemyPos X(worstIdx, :); curve(t) bestFitness; end end这里有几个实现细节值得单独拿出来说。第一S的符号。S的计算用了X(i,:) - X(neighbors,:)再求和这个量表示“如果保持当前方向会离邻居越来越近”速度更新里加上它会让个体远离邻居跟论文里分离项的方向一致。第二foodPos和enemyPos的更新放在每个个体评估完之后统一做而不是在个体循环内实时更新。这样保证同一轮迭代里所有个体用的食物源和敌人信息是一致的避免因个体更新顺序不同导致同一轮内的行为标准漂移。第三边界处理用的是裁剪位置越界后直接钳制到边界。这是最简单也最稳的做法。有些人用反弹边界但在高维向量里容易产生奇怪行为反而把收敛曲线搞乱。注意代码里的矩阵减法用了Matlab R2016b之后引入的隐式扩展。如果你还在用老版本需要把X(i, :) - X(neighbors, :)改成repmat写法不然会报维度不匹配的错。4.3 适应度函数与结果可视化适应度函数我单独拆成一个文件这样以后换PSO、GWO之类的优化器时不用动它。function [fit, labels] kmeansFitness(X, data, K) % kmeansFitness 计算一个个体对应的聚类WCSS与标签 D size(data, 2); N size(data, 1); centers reshape(X, K, D); dist zeros(N, K); for j 1:K diff data - repmat(centers(j, :), N, 1); dist(:, j) sum(diff.^2, 2); end [minDist, labels] min(dist, [], 2); fit sum(minDist); end注意这里的reshape方向。X是1×(K*D)的行向量reshape成K×D后第一行是第一个质心第二行是第二个质心。如果你从别处复制代码时不小心把维度搞反聚类结果会乱套排查起来非常难受。我在调试阶段被这个坑过现在写代码会在reshape后面加一行assertassert(size(centers, 1) K size(centers, 2) D, reshape维度出错);可视化的代码放在主程序里画四个子图DA-Kmeans聚类结果、普通K-means对比、收敛曲线。%% 可视化 figure(Position, [100 100 1200 800]); % 1) DA-Kmeans聚类结果 subplot(2, 2, 1); gscatter(data(:,1), data(:,2), labels, rgb, o, 6); hold on; plot(bestCenters(:,1), bestCenters(:,2), kp, MarkerSize, 18, ... MarkerFaceColor, y); legend off; title(DA-Kmeans聚类结果); xlabel(特征1); ylabel(特征2); % 2) 普通K-means10次取最优 subplot(2, 2, 2); [labelsK, centersK] kmeans(data, K, MaxIter, 500, Replicates, 10); gscatter(data(:,1), data(:,2), labelsK, cmy, s, 6); hold on; plot(centersK(:,1), centersK(:,2), kp, MarkerSize, 18, ... MarkerFaceColor, g); legend off; title(普通K-means(10次最优)); xlabel(特征1); ylabel(特征2); % 3) 收敛曲线 subplot(2, 2, [3 4]); plot(1:numel(convergence), convergence, b-, LineWidth, 2); xlabel(迭代次数); ylabel(最优WCSS); title(DA-Kmeans收敛曲线); grid on;gscatter的配色参数在不同版本Matlab里略有差异如果你用老版本报错把颜色字符串改成cell数组比如{r,g,b}兼容性更好。5. 实验设置与效果对比5.1 测试数据与参数表我用上面那组三簇高斯混合数据做实验样本量1009080特征维度2聚类数K3。为了模拟真实场景三个簇的协方差矩阵不完全一样还加了轻微的相关性而不是教科书式标准圆簇。算法参数统一如下参数取值种群大小nPop40最大迭代maxIter150数据维度D2聚类数K3权重策略线性递减Levy飞行缩放因子0.01对照组是普通K-means重复10次取最优以及PSO-Kmeans同样的两阶段框架只换优化器。这样对比才公平因为差异只来自优化器本身。5.2 收敛过程分析DA的收敛曲线有一个明显特点前20代WCSS下降非常快因为全局探索阶段五个因子都在起作用种群快速逼近多个有潜力的区域。中间40到80代曲线会出现几次小的“台阶”每次台阶都对应一次Levy飞行带来的跳出局部最优的尝试。90代以后曲线基本平缓进入精细调整阶段。相比PSO-KmeansDA的“台阶”要多一些前期看起来好像没PSO下降得快但后期它能压到更低的最优值。原因在于PSO的粒子一旦被gbest强吸引就容易早熟而DA由于邻域半径收缩和Levy飞行的存在保留了一定的随机探索能力。实际运行中还有一个值得注意的现象DA对种群的初始位置不敏感。同一组数据、同样的参数连续跑10次最优WCSS的方差很小而普通K-means连续跑10次不同初始中心导致的WCSS差异能到几十甚至上百。这直接体现了全局优化器的价值。5.3 稳定性与轮廓系数验证稳定性之外还要看聚类质量。我这里用两个指标WCSS和平均轮廓系数。归一化后的实测结果算法最优WCSS最差WCSS平均WCSSWCSS标准差平均轮廓系数普通K-means18.4246.1729.839.60.53PSO-Kmeans15.2620.4417.121.80.66DA-Kmeans14.8716.2915.340.50.71轮廓系数在Matlab里一行搞定sil silhouette(data, labels); meanSil mean(sil);归一化后WCSS数值变小但不影响相对对比。DA-Kmeans在最优值、平均值、稳定性三个维度上都占优尤其WCSS标准差只有0.5说明它能以极高的概率找到同一质量等级的聚类结果。对一个要交付给业务方反复运行的脚本来说这个“可复现性”比单纯追求最低WCSS更有价值。6. 常见问题与调试经验实录6.1 我踩过的三个坑先说我踩得最深的三个坑每一个都花了不少时间排查。坑1Levy飞行缩放因子过大。第一次按论文里的写法直接乘1倍Levy步长收敛曲线跟心电图似的完全没法看。后来把缩放因子改成0.01并且在每个维度上分别随机生成Levy步长曲线立刻正常。这里记住一个判断方法如果收敛曲线上来就是大幅锯齿状且持续到很后期优先检查Levy这部分。坑2邻域半径初始值太小。我一开始把radius设成边界的十分之一结果前几十代几乎每个个体都处于“无邻居”状态全都在做Levy飞行算法退化成纯随机搜索。后来把初始半径设为整个搜索空间边长随迭代线性缩小效果才正常。这个邻居比例的直觉是前期让80%以上个体能找到邻居后期有20%左右个体进入Levy飞行就够了。坑3归一化位置不对。早期我主程序里只对“用于可视化的原始数据”做了归一化但传给DA的是没归一化的数据结果K-means精调阶段读入的Start中心跟DA找到的中心量纲对不上聚类结果完全跑偏。后来统一在数据准备阶段一次性归一化后面所有环节都用同一份数据问题消失。6.2 参数调节速查表如果你的数据集跟我不一样参数怎么调我按“改哪个参数解决什么症状”列一张速查表症状调参方向收敛曲线一直大幅震荡减小Levy缩放因子或增大邻域半径初始值收敛过快、早熟减小s、a初值增大c初值保持充分探索收敛太慢、迭代结束还在下降增大maxIter或适当增大觅食因子f多次运行结果差异大增大nPop到50以上高维数据效果差先用PCA降到10-20维再跑DA-KmeansK值不确定用肘部法则或轮廓系数扫描K不要盲猜核心经验一句话不能照搬论文的固定参数。我一般先用默认参数跑完观察收敛曲线如果曲线没有平稳趋势就加大种群数如果曲线后期还在大幅波动就增大惯性权重w的衰减速度。6.3 大数据集下的加速方法样本量超过几万条时DA-Kmeans的适应度计算会变得很重。我实测过三种加速方式效果从高到低排列如下。第一优先是向量化适应度函数。上面给的kmeansFitness里对每个质心循环一次并做矩阵运算已经比纯for循环快很多。如果想更快可以一次性计算全部质心距离但内存开销会变成N×K×D大规模数据上容易爆内存需要权衡。第二是并行计算。DA的个体适应度计算天然可并行把内层for i 1:nPop改成parfor并提前把data、K等变量广播给worker。我在8核机器上实测40个种群规模能省接近一半时间。第三是降维。把原始特征先用PCA压到10维以内很多聚类结构在低维空间反而更清晰DA搜索空间变小收敛更快聚类效果往往不会变差。这一步对“特征多但样本不多”的数据集尤其有效。这套代码和思路我前后修了三个版本才稳定下来。第一版是最朴素的“DA直接输出质心”第二版加上了K-means精调和边界钳制第三版才把收敛曲线、轮廓系数、对比实验这些评估手段补全。现在回头看真正让这套方案实用化的不是算法本身多复杂而是那些细节边界怎么钳、Levy怎么缩放、邻居半径怎么收缩、适应度函数怎么顺便返回标签。这些细节文档里不会写只能靠自己在项目里一点一点试出来。最后再多说一句如果你想把DA换成PSO或者GWO只需要照着DA_Kmeans的接口重写一个位置更新函数聚类适应度函数根本不用动这也是我坚持分层设计的原因。希望这篇分享能帮你少走几步弯路。