蜻蜓算法优化Kmeans：Matlab聚类稳定性提升实践

做聚类分析的项目多了，你会慢慢体会到一件事：Kmeans这个算法，入门门槛确实低，但调起来相当心累。同样是iris数据集、同一个K值，换一次初始质心，聚类结果就可能完全两样，甚至把本该分开的两个簇硬生生合并成一个。最近我在Matlab里重做一批聚类实验，把蜻蜓算法（DA）和Kmeans结合起来做聚类分析，效果比我想象中扎实不少。这篇就把DA-Kmeans的核心原理、Matlab实现步骤、参数调法以及我踩过的坑一次性写清楚，给同样在折腾聚类方案的朋友一份可以直接参考的作业。

先说结论：蜻蜓算法优化的不是聚类过程本身，而是Kmeans的两个命门——初始质心选择和局部最优逃逸。它用群体智能的方式在解空间里搜索一组更合理的质心，再把这组质心交给Kmeans去精炼。这样组合下来，既能保留Kmeans的低开销特性，又能显著提高聚类结果的稳定性和质量。下面我按"为什么需要它—算法原理—代码实现—实验对比—调参避坑"这个顺序展开。

1. Kmeans的宿命：为什么聚类这件事总栽在初值手里

1.1 Kmeans的本质与SSE目标

先说清楚Kmeans到底在干什么。给定n个样本，每个样本有d维特征，Kmeans的任务是把这些样本分到K个簇里，使得每个样本到它所属簇中心的距离平方和最小。这个目标通常用SSE（Sum of Squared Errors）表示：

SSE = Σ Σ ||x_i - c_j||^2

其中第一个求和是对全部K个簇，第二个求和是对簇内所有样本；c_j是第j个簇的质心。

所以Kmeans本质是一个连续优化问题：需要找到K个质心位置，让所有样本到最近质心的总距离平方最小。Kmeans的标准做法是交替迭代：先根据当前质心给每个样本分配最近的簇，再根据簇内样本均值更新质心，不断重复直到质心不再变化。

这个过程本身有一个很好的性质：每一步都在客观上降低SSE，所以算法必然收敛。但"收敛"不等于"收敛到全局最优"，它更像一个盲人爬坡，只要脚下是下坡路就继续走，最后停在哪座山头完全取决于起点在哪里。

1.2 局部最优从哪来

Kmeans最大的问题可以概括为一句话：它是对初值极其敏感的贪心算法。随机初始化的质心如果落在某个数据密集区的角落里，后续迭代很可能被锁定在一个局部凹陷处，那里虽然是某个区域的局部低点，但远不是全局最佳划分。

我做过一个简单的统计实验：在同样的iris数据上，把Kmeans跑20次，每次都用不同的随机种子初始化质心，结果SSE从最小的97左右到最大的124都有，波动幅度接近30%。这个数字直接说明了问题——如果你只跑一次Kmeans就下结论，结果的好坏有一半取决于运气。

这也是为什么大家会尝试很多改进方案：Kmeans++、二分聚类、多起点随机重启。这些方法的本质都是"想办法让起点更好一点"或"多试几次看哪个结果最稳"。蜻蜓算法走的也是这条思路，但它不是盲目试，而是用群体智能在解空间里有方向地搜索，搜索效率比随机重启高得多。

2. 蜻蜓算法的仿生逻辑与更新公式拆解

2.1 五个行为因子的数学表达

蜻蜓算法（Dragonfly Algorithm，简称DA）是Seyedali Mirjalili在2016年提出的一种群体智能优化算法，灵感来自蜻蜓在自然界中的觅食和迁徙行为。蜻蜓在飞行过程中每个个体会表现出五种基本行为：分离、对齐、凝聚、觅食和躲避天敌。

这五个行为在算法里对应五个不同的调整向量，分别控制个体下一时刻的移动方向：

• 分离（Separation）：避免自己和邻居靠得太近。计算方式是当前蜻蜓与所有邻居之间向量差的负方向累加。直观理解就是"别挤我"。
• 对齐（Alignment）：让自己的飞行方向和邻居的平均方向保持一致。这是群体迁徙时形成整齐队形的关键，对应"跟大家一个方向飞"。
• 凝聚（Cohesion）：向邻居群体的重心靠拢。对应"别离队伍太远"。
• 觅食（Food Attraction）：向整个种群当前找到的最优位置靠近。这个最优位置在聚类优化里就是当前SSE最小的质心组合。
• 躲避天敌（Enemy Avoidance）：远离种群当前最差的位置，避免坠入已知的劣质解。在最优化语境下，天敌代表了"已知的坏解"。

每个行为都有对应的权重系数：s、a、c、f、e。蜻蜓的最终移动步长是这五个向量乘上各自权重后相加，再加上上一次移动步长的惯性项。位置更新公式可以写成：

deltaX(t+1) = s·S_i + a·A_i + c·C_i + f·F_i + e·E_i + w·deltaX(t)

X(t+1) = X(t) + deltaX(t+1)

这里w是惯性权重，用于平衡全局搜索和局部搜索。w大时蜻蜓倾向于沿原来的方向飞行，全局探索能力强；w小时蜻蜓受当前邻居和食物影响更大，局部开采能力更强。实际实现中通常让w从0.9线性衰减到0.2，前期多做全局搜索，后期专注精细收敛。

2.2 位置更新与适应度对接

用DA去优化Kmeans，最关键的一步是搞清楚"蜻蜓的位置"和"Kmeans的解"之间怎么映射。

在标准DA里，每个蜻蜓个体的位置是一个D维向量，D就是待优化问题的变量个数。对于Kmeans来说，一个完整解就是K个质心坐标拼接起来，每个质心维度等于特征维度d，所以单个蜻蜓的位置向量长度就是K×d。

举个例子：iris数据集有4个特征，设定K=3，那么每个蜻蜓个体的编码长度就是12。前4个数是第一个簇的质心，中间4个数是第二个簇的质心，最后4个数是第三个簇的质心。解码时只需要对这个向量做一个reshape(3, 4)，就能得到完整的质心矩阵。

适应度函数直接采用Kmeans的目标：计算当前质心划分下所有样本到最近质心的SSE。SSE越小，说明这组质心越好，这个蜻蜓个体就越接近食物位置。

这种编码方式的好处是与算法天然解耦：DA只负责在质心空间里搜索"更小的SSE"，完全不关心聚类过程的内部细节；Kmeans只负责在DA给出的质心基础上做一次标准迭代精炼。两端各做各的，替换起来也方便——你想把Kmeans换成模糊C均值，只需要改适应度函数里的目标计算就可以了。

3. 把蜻蜓塞进Kmeans：编码方式与Matlab实现步骤

3.1 个体编码与维度计算

编码这一步是很多新手翻车的地方，我必须单独强调一次。

如果你的样本X是n行d列的矩阵，K是聚类数，那么每个蜻蜓个体不是K个质心的集合，而是这K个质心按行拼接后展成的一维行向量。也就是说，个体维度dim = K * d。

实用中我会用两个辅助变量：

• ub：每个维度取值的上界，通常设为样本数据各特征的最大值向量，复制K份拼成一个1×dim的行向量；
• lb：每个维度取值的下界，设为样本各特征的最小值向量，同样复制K份。

这样初始化时直接用X = lb + rand(pop, dim) .* (ub - lb)就能生成一个初始种群，把每个蜻蜓个体的初始位置限制在样本范围内，避免一开始就跑到数据范围外去。

3.2 主循环核心代码

先看适应度函数。输入一个蜻蜓个体的位置向量x，把它解码成质心矩阵，再计算SSE：

function sse = calSSE(x, data, K, d) centers = reshape(x, K, d); n = size(data, 1); distMat = zeros(n, K); for i = 1:n distMat(i, :) = sum((repmat(data(i, :), K, 1) - centers).^2, 2); end [minDist, ~] = min(distMat, [], 2); sse = sum(minDist.^2); end

注意这里minDist直接取的是每个样本到最近质心的欧氏距离，平方之后再累加，得到的就是标准SSE。

主循环我按原始DA的思路实现了一个简化但完整的版本，使用全局邻居假设，即所有蜻蜓互相视为邻居：

function [bestCenters, bestSSE, history] = DA_Kmeans(data, K, pop, maxIter) [n, d] = size(data); dim = K * d; lbAll = repmat(min(data), 1, K); ubAll = repmat(max(data), 1, K); % 初始化种群和速度 X = lbAll + rand(pop, dim) .* (ubAll - lbAll); deltaX = zeros(pop, dim); fitness = zeros(pop, 1); w0 = 0.9; wEnd = 0.2; s = 0.1; a = 0.4; c = 0.7; f = 0.8; e = 0.2; history = zeros(maxIter, 1); for iter = 1:maxIter w = w0 - (w0 - wEnd) * iter / maxIter; for i = 1:pop fitness(i) = calSSE(X(i, :), data, K, d); end [bestSSE, bestIdx] = min(fitness); [~, worstIdx] = max(fitness); food = X(bestIdx, :); enemy = X(worstIdx, :); neighborR = max(ubAll - lbAll); % 全局邻域半径 for i = 1:pop % 分离项：与所有邻居反向 S = sum(repmat(X(i,:), pop, 1) - X, 1) / pop; S = -S / (norm(S) + eps); % 对齐项：邻居平均速度与自身速度差（简化取种群平均位置） avgPos = mean(X, 1); A = (avgPos - X(i,:)) / (norm(avgPos - X(i,:)) + eps); % 凝聚项：朝邻居重心 C = (mean(X, 1) - X(i,:)) / (norm(mean(X,1) - X(i,:)) + eps); % 觅食项 F = (food - X(i,:)) / (norm(food - X(i,:)) + eps); % 躲避天敌项 E = (X(i,:) - enemy) / (norm(X(i,:) - enemy) + eps); % 更新步长和位置 deltaX(i, :) = s*S + a*A + c*C + f*F + e*E + w*deltaX(i, :); X(i, :) = X(i, :) + deltaX(i, :); % 边界反弹处理 X(i, :) = max(min(X(i, :), ubAll), lbAll); end history(iter) = bestSSE; end [bestSSE, bestIdx] = min(fitness); bestCenters = reshape(X(bestIdx, :), K, d); end

这段代码有几个地方我要额外解释。

首先是"分离项"的计算。标准DA需要知道每个邻居的精确位置，我是用所有个体与当前个体的差向量累加后取平均，再取反方向归一化。这是一种全局邻居假设的简化，在不追求论文级复现的前提下完全够用，而且省掉了计算邻域矩阵的开销。

其次，我把对齐项和凝聚项都做了归一化。这一步很关键：原始DA公式里这些向量长度不同，直接加权相加时如果某个项模长特别大，会盖过其他项。归一化之后，每个行为因子只保留方向信息，权重s、a、c、f、e才能真正起到调节作用。

最后是边界反弹。这个处理容易被忽略，但如果不做，蜻蜓的质心坐标会随着迭代漂移到数据范围之外，导致适应度计算出现空簇或异常值。max(min(...), lbAll)这种写法等效于把越界分量拉回到边界上，简单有效。

3.3 参数初始化的建议

初始化时我一般把种群规模pop设在25到40之间，迭代次数maxIter设在100到200之间。这个量级在iris这种小数据集上只需要两三秒就跑完，但已经足够让DA找到接近全局最优的质心组合。

如果你处理的数据量比较大，建议先用小种群快跑一轮，观察适应度下降曲线是否在迭代后期仍明显下降。如果还在下降，说明迭代次数不够；如果中后期基本平坦，那继续增大迭代次数意义不大，应该增加种群规模扩大搜索广度。

4. 实验对比：DA-Kmeans和传统Kmeans差多少

4.1 测试数据集与评价指标

实验我选了三个经典数据集：iris（150×4，3类）、wine（178×13，3类）以及人工生成的四个高斯簇混合数据（400×2，4类）。这三个数据的共同点是已知真实簇数，适合用来检验聚类算法的稳定性和准确性。

评估指标除SSE之外，我同时记录了两个维度：一是50次独立运行里最优结果的SSE分布，用来衡量"能不能找到好解"；二是最优解对应的准确率（与真实标签对比的匹配度），用来衡量"好解是否真的好用"。

这里说句体外话：SSE低不代表聚类结果一定符合你的预期，因为它本质上是在衡量紧致性，没有考虑簇与簇之间的分离度。但在已知簇数和数据分布的情况下，SSE足够作为对比基准。

4.2 收敛性与稳定性对比

传统Kmeans的对照组设置是随机初始化跑50次；DA-Kmeans的对照组是种群30、迭代120、每个参数相同条件下跑50次。结果如下表：

最直观的差异在标准差：传统Kmeans的SSE在多次运行之间剧烈跳动，而DA-Kmeans的SSE几乎贴在同一个数值附近。这说明DA通过群体搜索，把"随机初始化碰运气"变成了一种可重复的稳定过程。

稳定性的提升在实操中意义很大。尤其是在业务场景里跑聚类脚本，你不可能每次跑完都人工检查一遍簇划分是否合理。如果算法每个批次输出的结果都差不多，后续的报表、标注、数据清洗才能放心地建立在聚类结果之上。

4.3 收敛性曲线怎么看

DA-Kmeans的收敛曲线有一个典型特征：前二三十代适应度快速下降，中后期趋于平缓。这是因为惯性权重w从0.9线性衰减到0.2，前期蜻蜓飞行范围大，搜索激进，很快就能发现一片优质区域；后期步长变小，群体围绕最优解做精细搜索。

如果收敛曲线在中后期仍然频繁出现明显起伏，通常有几种原因：步长因子s或躲避天敌因子e设得太大，导致个体反复跨过最优区域；或者种群多样性不足，所有个体挤在某个局部最优附近无法跳出。

对比传统Kmeans的"收敛"，差别在于Kmeans的收敛是机械式的——迭代到局部低点后stop；而DA的收敛是有逃逸能力的——即使接近当前最优，仍有可能通过群体协作和随机步长跳出陷阱。这恰好补上了Kmeans最弱的一环。

5. DA-Kmeans的调参与避坑记录

5.1 五个权重因子分别该怎么调

在DA算法里，s、a、c、f、e这五个权重直接决定了个体飞行时受哪种行为影响更大。我给出一个常规的调参起点和调整方向：

我自己的经验是：先固定f和c在一个中等偏高的水平，把s和a调小一点，让群体有足够的探索自由；等收敛曲线稳定后，再逐步增大f和c提高开采强度。这个顺序比五个参数一起乱试要高效得多。

5.2 容易翻车的几个场景

第一个坑是质心跑到数据范围之外。如果你像我最初那样忽略边界处理，运行20代后整个种群的位置向量可能包含很多远离样本云的坐标，导致某些簇一个样本都分不到。解决办法就是在每次更新后强制把所有个体拉回到lb和ub之间。

第二个坑是不管聚类数K。DA-Kmeans可以优化质心位置，但它不能替你决定K应该取几。实际操作中必须先用肘部法则、轮廓系数或业务先验确定K的值，再把K喂给算法。如果K选得不对，DA使再大劲也只是在错误的约束下找最优解。

第三个坑是空簇惩罚缺失。DA搜索质心的方式是连续的，因此在接近收敛时可能出现两个质心几乎重合的情况，这会让其中一个簇的样本数变成0。如果不处理，SSE的计算结果会给出一个"很漂亮"但实际上无意义的极小值。我的做法是在适应度函数里加一个软惩罚：如果某个簇的样本数小于1，就把SSE加上一个固定的大数（比如样本方差的100倍），迫使群体避开这种退化解。

第四个坑是计算复杂度。每一代迭代中，每个蜻蜓个体都要计算一次全样本到K个质心的距离矩阵，复杂度是O(pop × maxIter × n × K × d)。当n到万级、d到几十维时，运行时间会变得相当可观。这种情况下我建议先用小规模数据调参，再放大数据量；或者把SSE计算向量化，比如用pdist2替代循环计算距离矩阵，速度提升非常明显。

5.3 混合策略：让Kmeans自己收个尾

最后分享一个我觉得性价比最高的做法：不要用DA完全替代Kmeans，而是让DA负责"找好起点"，Kmeans负责"精细收敛"。

具体操作是：先用DA跑完，拿到一组最优质心，把它作为Kmeans的Start参数传入，再调用一次标准Kmeans迭代到收敛。这样DA阶段负责大范围搜索，Kmeans阶段负责在优秀起点附近做精确下降，两个算法各干各擅长的事。

options = statset('MaxIter', 200); [bestCenters, ~] = DA_Kmeans(data, K, 30, 120); [idx, centers] = kmeans(data, K, 'Start', bestCenters, 'MaxIter', 200, 'Options', options);

这种混合策略的收敛速度比纯DA快，稳定性比纯Kmeans高，而且代码改动量很小。我在实际项目中基本都采用这个方案，效果比单纯堆迭代次数好很多。

写在最后

如果你正在写论文的对比实验，或者只是想把聚类结果做得更稳，DA-Kmeans搭配一个"DA搜起点+Kmeans收尾"的流程，绝对值得一试。调试时记住一个简单的小技巧：把每一代食物位置的SSE画成曲线，如果曲线明显震荡，优先调小分离步长s和躲避天敌系数e；如果曲线过于平稳但最终SSE不够小，说明群体多样性不足，试试增大分离系数或者把每个蜻蜓个体按随机子集计算邻居。我在实际运行中靠这个方法节省了大量盲调参数的时间，希望它对你也有用。