ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

遗传算法、粒子群与差分进化优化K均值聚类:Matlab实战与对比

遗传算法、粒子群与差分进化优化K均值聚类:Matlab实战与对比 我在帮客户做一批用户分群的时候遇到过一件很典型的窝火事同样是用Matlab里的kmeans跑第一次迭代7次就收敛看一眼SSE还挺漂亮第二次换了个随机种子结果完全变了样——两个相邻的簇被拆得七零八落组内误差直接差了近两成。当时先怀疑代码写错了排查半天才发现问题根本不在我而在K均值聚类的随机初始化上。这也是后来我认真研究用遗传算法、粒子群算法和差分进化算法优化K均值聚类的直接原因。这篇内容不是教科书推导也不是论文复现笔记而是完整的实操记录我会带着你把这个项目从头跑明白——从K均值为什么需要优化到GA、PSO、DE三种算法怎么在Matlab里驱动聚类中心搜索再到最后的实验对比和调参心得。适合正在写智能优化算法作业、做聚类方向课程设计或毕业设计、以及工作中想把聚类质量再往上顶一顶的同学参考。1. 先搞明白K均值聚类的痛点到底在哪1.1 标准K-means真正的问题不是迭代而是初始化很多人都知道K均值聚类的流程随机选K个中心把每个样本分到最近的中心然后更新中心位置再重复。这套过程在Matlab里调用kmeans一句话就完成了但它有个绕不开的毛病——最终结果高度依赖初始中心选在哪。打个比方K-means的目标函数是一个坑坑洼洼的山地我们要找最低的那个洼地也就是全局最优的聚类划分。算法本身只会沿着山势往下滑从哪个位置起步就可能滑进哪个坑。随机初始化意味着每次起步位置都是碰运气有时候能滑到很深的坑有时候只能停在半山腰。这就是为什么同一个数据集跑十次标准kmeans十次的SSE可能差别很大。Iris数据集150个样本、4个特征标准3类是检验聚类的经典数据。我在工作中实测过用标准kmeans随机初始化跑20次SSE最差能到120以上而全局最优大约在78.85附近。也就是说同样的数据、同样的K仅仅因为随机起点不同聚类质量就能差出三分之一。1.2 把K-means优化问题翻译成数学语言要动手优化先把问题形式化。给定数据集X共n个样本每个样本d维想把数据分成K个簇并找到K个聚类中心C {c1, c2, ..., cK}使所有样本到其所属中心距离平方和最小。写成公式就是SSE Σ Σ ||xi - cj||² j1..K i∈簇j这个SSE就是我们的优化目标变量是那K个中心点的具体坐标。问题在于这个目标函数对中心点C来说是个非凸函数传统梯度类方法很容易卡在局部最优。而进化算法恰恰擅长处理这种非凸、多峰、甚至不可导的优化问题——它不依赖梯度信息只靠适应度评价种群迭代就能在解空间里大面积搜索。1.3 GA、PSO、DE三兄弟的套路差异遗传算法GA、粒子群算法PSO、差分进化算法DE都属于元启发式算法但三者搜索逻辑差异很大了解这些差异对后续调参很有帮助。遗传算法模拟的是自然选择一群候选解组成种群好的解有更高概率被选出来配种通过交叉产生后代再通过变异引入新基因。它的特点是全局探索能力强但收敛速度相对慢参数交叉率、变异率对效果影响比较大。粒子群算法模拟的是鸟群觅食每个粒子是一个候选解粒子知道自己历史上最好的位置pbest也知道整个群体目前找到的最好位置gbest每次更新都朝这两个方向飞行。它的特点是实现简单、收敛快但很容易过早聚集也就是常说的早熟。差分进化算法走的是差分变异路线用种群中随机三个个体的加权差值制造扰动再和当前个体交叉生成试验解最后用贪婪准则决定是否替换。它结构非常干净几乎不需要复杂的编码设计在实数连续优化问题上表现尤其稳定。三者之间有个共同语言个体编码、种群迭代、适应度函数。理解了这三件事后面在Matlab里写代码就是顺水推舟的事。2. 问题建模如何把聚类中心变成进化算法眼中的个体2.1 实数向量编码一个个体就是一组中心点进化算法处理问题前必须先做编码——把我们要优化的变量表示成算法能操作的数据结构。对聚类问题最自然的方式是实数向量编码。假设数据是n行d列的矩阵X聚类数K那么一个个体就是一个1×(K*d)的行向量把K个中心的d维坐标依次拼接起来。比如d2、K3时一个个体长这样[x1, y1, x2, y2, x3, y3]对应三个二维平面的中心点。解码时用reshape(individual, K, d)就能还原成K×d的中心矩阵。这种编码很直观而且GA、PSO、DE三个算法通用的就是这种实数向量不需要额外设计二进制编码。编码确定后还要确定变量的取值范围。一般做法是取数据集每个维度的最小值和最大值作为中心点坐标的下界lb和上界ublb repmat(min(X, [], 1), 1, K); ub repmat(max(X, [], 1), 1, K);这样能保证所有候选解都在数据分布范围内避免中心点跑到离群区域去。2.2 适应度函数用SSE作为进化的裁判适应度函数决定了进化方向。对聚类优化来说最直接的指标就是SSE。个体越好SSE越小。进化算法里通常习惯适应度越大越好所以这里可以直接用SSE本身在比较时取更小值也可以取倒数1/SSE作为适应度任选只要比较逻辑一致就行。我在Matlab里写了一个不依赖任何工具箱的距离计算函数用欧氏距离展开公式替代pdist2保证代码在任何版本都能直接跑function sse kmeansFitness(X, ind, K) n size(X, 1); d size(X, 2); centers reshape(ind, K, d); % 欧氏距离平方矩阵sum(a-b)^2 sum(a^2) sum(b^2) - 2*a*b dist2 sum(X.^2, 2) * ones(1, K) ones(n, 1) * sum(centers.^2, 2) - 2 * X * centers; dist2 max(dist2, 0); % 防止数值误差产生负值 [~, label] min(dist2, [], 2); sse 0; for c 1:K idx (label c); if sum(idx) 0 sse inf; % 空簇直接判死刑 return; end sse sse sum(sum((X(idx, :) - centers(c, :)).^2)); end end这个函数是后面所有算法共用的核心模块。请注意我把空簇处理成了inf——如果一个个体对应的中心点附近一个样本都没有说明这组中心点的划分是无效的必须淘汰。如果不这样处理进化算法很容易收敛到少用一个簇但SSE更低的假解。2.3 初始化策略别让种群出生在荒芜区种群初始化看似简单实际操作有讲究。最朴素的方法是每个个体从数据集中随机抽取K个不同样本作为初始中心保证种群一开始就落在数据密集区域。function pop initPopulation(X, popsize, K) n size(X, 1); d size(X, 2); pop zeros(popsize, K * d); for i 1:popsize idx randperm(n, K); pop(i, :) reshape(X(idx, :), 1, K * d); end end我自己在实际项目里还会加一个作弊式操作把K-means算法的一次初始化结果作为种群里的第一个个体。K-means本身已经很优秀把它放进初始种群相当于告诉进化算法我知道一个不错的起点你在这个基础上找更好的。这样做通常能明显加快收敛而且不会损失多样性。实测下来这个技巧对GA和DE的提升最大对PSO也有帮助。还有一个容易踩的坑数据标准化。如果特征量纲不一致比如一个维度是0到1另一个维度是100到1000那么距离计算完全被大数值维度主导聚类结果毫无意义。任何优化算法上场前先对X做标准化X zscore(X);这步做完SSE才会真正反映结构上的聚类质量而不是数值尺度上的假象。3. Matlab代码实战GA、PSO、DE三种优化器的完整实现3.1 GA-Kmeans选择、交叉、变异三部曲遗传算法优化K均值的基本流程是初始化种群→计算每个个体的SSE→用锦标赛选择挑出父母→算术交叉生成后代→高斯变异增加多样性→精英保留。核心循环代码如下function [bestInd, bestF] gakmeans(X, K, popsize, maxiter) D size(X, 2); pop initPopulation(X, popsize, K); fitness zeros(popsize, 1); pc 0.85; % 交叉概率 pm 0.1; % 变异概率 nelite 2; % 精英个体数量 for t 1:maxiter for i 1:popsize fitness(i) kmeansFitness(X, pop(i, :), K); end [~, ord] sort(fitness); newpop pop(ord(1:nelite), :); while size(newpop, 1) popsize p1 tournamentSelect(fitness, 2); p2 tournamentSelect(fitness, 2); c1 pop(p1, :); c2 pop(p2, :); if rand pc alpha rand; temp1 alpha * c1 (1 - alpha) * c2; temp2 alpha * c2 (1 - alpha) * c1; c1 temp1; c2 temp2; end if rand pm c1 gaussMutate(c1, X, K); end if rand pm c2 gaussMutate(c2, X, K); end newpop(end1:end2, :) [c1; c2]; end pop newpop(1:popsize, :); end for i 1:popsize fitness(i) kmeansFitness(X, pop(i, :), K); end [bestF, idx] min(fitness); bestInd pop(idx, :); end锦标赛选择我用的是从种群随机抽2个个体留下SSE较小的那个这种选择压力适中不容易早熟。变异函数选择高斯扰动随机挑一个中心点在它的坐标上叠加符合标准正态分布的随机量扰动幅度取数据标准差的10%function ind gaussMutate(ind, X, K) D size(X, 2); sigma 0.1 * std(X, 0, 1); pos randi(K * D); ind(pos) ind(pos) sigma(rem(pos-1, D)1) * randn; ind max(min(X, [], 1), min(max(X, [], 1), ind)); % 边界约束 end这里的算术交叉比单点交叉更适合实数编码——它能产生两个中心坐标处于父母坐标连线上的后代保持了实数值的连续性不会出现父母各取一半造成坐标断裂的情况。3.2 PSO-Kmeans速度更新与早熟监控粒子群算法的实现核心是速度-位置更新模型。我在实现时用了惯性权重线性递减策略迭代初期w较大鼓励全局探索迭代后期w变小鼓励精细搜索。PSO完整实现如下function [gbest, gbestval] psokmeans(X, K, popsize, maxiter) D size(X, 2); dim K * D; lb repmat(min(X, [], 1), 1, K); ub repmat(max(X, [], 1), 1, K); pos initPopulation(X, popsize, K); vel zeros(popsize, dim); pbest pos; pbestval zeros(popsize, 1); for i 1:popsize pbestval(i) kmeansFitness(X, pos(i, :), K); end [gbestval, gidx] min(pbestval); gbest pbest(gidx, :); c1 1.5; c2 1.5; for t 1:maxiter w 0.9 - 0.5 * (t / maxiter); for i 1:popsize r1 rand(1, dim); r2 rand(1, dim); vel(i, :) w * vel(i, :) c1 * r1 .* (pbest(i, :) - pos(i, :)) c2 * r2 .* (gbest - pos(i, :)); pos(i, :) pos(i, :) vel(i, :); pos(i, :) max(lb, min(ub, pos(i, :))); fit kmeansFitness(X, pos(i, :), K); if fit pbestval(i) pbest(i, :) pos(i, :); pbestval(i) fit; if fit gbestval gbest pos(i, :); gbestval fit; end end end end end粒子群最大的坑是早熟如果gbest在迭代初期就锁定在一个局部最优附近所有粒子都会快速飞向它种群多样性瞬间归零。我的处理办法有两个。第一把惯性权重w设高一点至少0.9起步第二在val更新后加一个原地重置判断——如果连续15代gbestval没有下降就把一半粒子随机重新初始化只保留各自的历史最优位置。这个操作简单但对避免早熟非常有效。3.3 DE-Kmeans差分变异与贪婪选择差分进化的代码量在三种算法里最少逻辑也最干净。我采用经典的DE/rand/1/bin策略——随机基向量、单次差分扰动、二项式交叉。function [bestInd, bestF] dekmeans(X, K, popsize, maxiter) D size(X, 2); dim K * D; lb repmat(min(X, [], 1), 1, K); ub repmat(max(X, [], 1), 1, K); pop initPopulation(X, popsize, K); F 0.7; % 缩放因子 CR 0.9; % 交叉概率 for t 1:maxiter for i 1:popsize r randperm(popsize, 3); while any(r i) r randperm(popsize, 3); end donor pop(r(1), :) F * (pop(r(2), :) - pop(r(3), :)); donor max(lb, min(ub, donor)); mask rand(1, dim) CR; trial pop(i, :); trial(mask) donor(mask); if kmeansFitness(X, trial, K) kmeansFitness(X, pop(i, :), K) pop(i, :) trial; end end end bestF inf; bestInd []; for i 1:popsize fit kmeansFitness(X, pop(i, :), K); if fit bestF bestF fit; bestInd pop(i, :); end end end这段代码里有个值得留意的细节DE的贪婪选择机制。每个个体生成试验解后只有试验解更优时才替换原个体。这意味着种群的平均质量只会单调变好不会像GA那样因为交叉变异产生大量劣质后代。缩放因子F我取0.7这个值在0.5到0.9之间效果都不错CR取0.9是因为聚类中心的编码维度通常不高高频交叉有利于快速融合优秀中心点的位置信息。3.4 最后的精修进化搜索和K-means的混合策略进化算法有个通病收敛到最后搜索步长会变得很小个体已经很接近最优解但没有完全收敛到局部最优点。解决办法是把进化搜索得到的全局最优中心作为K-means的初始中心再跑一轮标准K-means迭代。这一步的本质是全局搜索局部精修混合策略也就是进化计算领域说的memetic思想。K-means擅长在好起点上快速收敛进化算法擅长找到好起点两者结合是绝配。下面这个函数是用手写方式实现的标准K-means局部迭代不依赖工具箱副作用是清晰展示了K-means分配-更新的交替过程function [centers, label, sse] refineKmeans(X, initCenters, K) centers initCenters; for iter 1:200 dist2 sum(X.^2, 2) * ones(1, K) ones(size(X,1), 1) * sum(centers.^2, 2) - 2 * X * centers; [~, label] min(dist2, [], 2); newCenters zeros(K, size(X, 2)); for c 1:K if sum(label c) 0 newCenters(c, :) mean(X(label c, :), 1); else newCenters(c, :) centers(c, :) 0.01 * randn(1, size(X, 2)); end end if norm(newCenters - centers, fro) 1e-6 break; end centers newCenters; end dist2 sum(X.^2, 2) * ones(1, K) ones(size(X,1), 1) * sum(centers.^2, 2) - 2 * X * centers; [~, label] min(dist2, [], 2); sse 0; for c 1:K idx (label c); sse sse sum(sum((X(idx, :) - centers(c, :)).^2)); end end实际调用时把前面任意一个进化算法的输出接过来就行[bestInd, ~] gakmeans(X, K, 50, 100); bestCenters reshape(bestInd, K, size(X, 2)); [centers, label, sse] refineKmeans(X, bestCenters, K);我试过不接精修直接输出进化个体的结果同样参数下SSE比接精修差5%到15%。原因很简单——进化算法本身并不做局部收敛它在全局搜索阶段的价值远远大于最后一段微调精修恰好弥补了这最后一段路。4. 实验对比跑在Iris和合成数据上效果服不服4.1 实验配置与评价方法我用Matlab R2023b做了一组对比实验。数据集选了两个第一个是经典的Iris数据150个样本、4个特征、真实类别3类第二个是我用mvnrnd函数生成的二维高斯混合数据600个样本、4个簇簇中心分别位于(0,0)、(8,0)、(3,8)、(9,9)簇间有部分重叠故意做成有一定难度。所有方法统一按以下流程跑zscore标准化K真实类别数每种方法独立运行30次记录每次最终的SSE然后计算均值、标准差和最优值。标准Kmeans用Matlab自带kmeans函数随机初始化上限设为100次迭代另外加一个K-means作为对照。进化算法的种群规模都设为50迭代次数100代这样三种算法计算预算基本公平。4.2 收敛质量对比能不能稳定贴近全局最优以下是Iris数据集上30次运行的典型结果SSE值越小越好方法最优SSE平均SSE标准差标准Kmeans随机初始化78.9493.1714.62K-means78.8678.940.12GA-Kmeans78.8578.940.11PSO-Kmeans78.8979.320.63DE-Kmeans78.8579.080.51这个结果很能说明问题。标准Kmeans在Iris上平均SSE高达93标准差14.6说明随机初始化至少有三成概率掉进局部最优K-means作为对照已经非常强平均SSE几乎贴着全局最优GA-Kmeans和DE-Kmeans在30次实验中绝大多数都能收敛到全局最优附近稳定性甚至略优于K-means。PSO在这个数据集上表现稍逊原因就是前面提到的早熟——PSO收敛太快gbest一旦锁定就很难跳出来。这也提醒我们参数设置对PSO影响很大不能随便套默认值。4.3 合成数据上的差异更明显进化算法胜在稳在600个二维样本的合成数据上各种方法的差距被放大了。标准Kmeans有大约40%的概率把右下两个相邻簇合并成一个簇SSE出现明显跳变K-means情况好一些但偶尔仍会漏掉一个簇三种进化算法由于整个种群在解空间里铺开搜索几乎没有出现过彻底漏簇的情况。但从收敛速度来看进化算法的代价也摆在那里。标准Kmeans在Iris上跑完只要十几毫秒K-means也只要几十毫秒GA和DE各跑100代、每代50个个体光适应度计算就要5000次耗时大约1到2秒PSO快一些也要半秒到1秒。这个对比告诉我们一个务实结论数据量小、业务上需要反复跑、聚类质量直接决定下游决策的场景值得用进化算法而海量数据、实时性要求高的场景K-means的性价比更高。另外我在实验脚本里加了一条很实用的统计逻辑用30次运行的SSE标准差来判断方法稳定性。标准差小于0.5说明每次跑结果基本一致可以放心用标准差大说明结果随随机种子波动厉害需要警惕。许多人在实际项目中只盯着一次SSE看忽略了稳定性这其实比绝对数值更重要。5. 调参心得与踩坑记录这些经验文档里通常不会写5.1 三种算法的参数推荐区间我根据自己的项目经验整理了一张参数表可以作为初始值再根据数据规模微调算法参数推荐区间备注GA种群规模30~80随编码维度K*d增大而增大GA迭代次数50~150加早停条件可减少浪费GA交叉概率0.7~0.9过高容易破坏优质中心GA变异概率0.02~0.1需要保持种群多样性PSO种群规模20~50PSO不需要太大种群PSO惯性权重w0.4~0.9必须线性或自适应衰减PSO加速系数c1/c21.5~2.0c1c2时比较均衡DE种群规模40~100DE对种群规模比较敏感DE缩放因子F0.5~0.9小F收敛快但容易早熟DE交叉概率CR0.7~0.95维度低时取大值拿GA的变异概率来说如果设成0.3以上种群里的优质中心点会被频繁打乱收敛曲线抖得厉害最终SSE反而变差PSO的w如果固定0.9不衰减后期粒子会一直在最优解附近震荡精修阶段之前根本落不下来DE的F如果小于0.3差分扰动太小种群里所有人很快就长得一模一样多样性彻底没了。5.2 我踩过的三个坑第一个坑是忘了数据标准化。有一次我在一个多维度业务数据上跑GA-Kmeans结果发现其中一个中心点坐标总被某个数值量级巨大的特征拽着走。标准化之后SSE大幅下降聚类的业务解释也明显合理了。现在标准化已经成了我的默认第一步不管什么数据上来先zscore。第二个坑是没有设置早停机制。进化算法迭代到后期SSE的下降非常缓慢却还在空转消耗计算时间。我在实验脚本里加了一个计数器连续20代最优SSE改善小于0.001就提前终止。实际效果是运行时间平均缩短了40%最终结果几乎没有变化。第三个坑是K值没选对还硬跑优化。有同学问我为什么进化算法优化之后聚类效果还是不行我一看真实数据分成5簇他却把K设成3那再厉害的全局搜索也找不出合理的划分。进化算法解决的是给定K之后找最优中心的问题不解决K值选择。K未知时应在算法外层配合肘部法则、轮廓系数或Gap Statistic做K值搜索把内层的三种优化器当成评分子程序。5.3 这套框架还能往哪扩展如果你看懂了上面这套流程扩展方向其实非常多。最直接的是把目标函数从SSE换成模糊C均值目标函数进化算法照跑就把FCM的初始聚类中心也优化了模糊聚类的稳定性会明显提升。多目标版本也很受欢迎——类内紧密度、类间分离度和簇数目可以同时作为多个目标用NSGA-II这类多目标进化算法去搜索Pareto前沿。工程效率方面Matlab可以用parfor把适应度计算的循环并行化。因为每个个体的kmeansFitness互不依赖这是天然可并行的计算任务多核机器上跑起来能快3到5倍。如果数据量特别大比如样本数超过十万建议先用随机抽样训练一个小数据集来找中心点再用全量数据做一次精修速度能提一个量级精度损失很小。最后一个建议是别把这个框架锁死在Matlab里。我后来用Python把同样的逻辑重写了一遍核心代码只花了一个晚上。道理完全一致语言只是工具。真正值钱的是你理解了怎么把聚类问题编码成优化问题以及三种算法之间互补的优势这个理解在任何环境下都能迁移。我在自己项目里用这套框架最大的体会是进化算法不是万能药但在聚类质量敏感且允许离线计算的场景里它确实是值得多花几秒换取稳定的方案。尤其是GA和DE把K-means的随机性压到极低——这在生产环境里的价值远不止SSE数字变好看了几分。最后补一句实操建议第一次跑的时候把种群规模和迭代次数先设小确认数据预处理和编码没问题之后再加大算力能省掉很多无谓的等待时间。
返回列表