智能优化算法优化核极限学习机预测:原理与Matlab实现
做预测类课题的同学十有八九都干过这么一件事换一个智能优化算法去调一个机器学习模型的参数然后发一篇“某某算法优化某某模型预测”的论文。哈里斯鹰算法、鲸鱼算法、粒子群算法、蝴蝶算法再加上核极限学习机这几个词往一起一放基本就是一个标准的“优化KELM预测”组合拳。这篇博文就把这套东西掰开揉碎讲一遍重点说清楚四类元启发式算法在优化核极限学习机时各自怎么运作、Matlab代码怎么写、实际跑起来会遇到哪些坑以及怎么设计一次公平的对比实验。内容面向正在做回归预测、软测量、时间序列预测的研究生和工程师也适合想快速搭建一套可复现预测流水线的朋友。1. 先把KELM搞清楚为什么偏偏优化C和γ1.1 ELM到KELM一个直观的演化极限学习机的核心思路其实非常暴力。传统BP神经网络靠梯度下降一层层反传误差慢且容易陷入局部最优。ELM直接把隐层输入权重和偏置随机生成、保持固定然后只求解输出层权重整个训练过程变成一个最小二乘问题一步到位。假设训练样本是XN行d列随机生成输入权重W和偏置b之后隐层输出矩阵H就是N行L列L为隐层节点数那么输出权重β可以直接算出来beta pinv(H) * Y;如果不加正则化这就是最小二乘解。但ELM有一个尴尬之处隐层节点数L怎么选选少了欠拟合选多了过拟合L本身又变成一个需要调的参数。KELM的思路更聪明——既然H是随机映射那干脆不要H了直接引入核函数用核矩阵Omega替代H*H避开了显式隐层映射也就不用纠结隐层节点数。KELM的训练就是在求解Omega kernel_matrix(X, X, gamma); Beta (Omega eye(N) / C) \ Y;这里Omega是N×N的核矩阵C是正则化系数gamma是核参数。预测时对一个新的测试样本x先计算它和所有训练样本的核向量Ktest kernel_matrix(x, X, gamma)再乘上Beta即可。1.2 核参数与正则化系数的实际含义KELM的两个关键参数C和gamma一个决定模型复杂度惩罚的强度一个决定核函数的“瘦高”程度。C是正则化系数。C越小对输出权重β的惩罚越重模型越平滑防止过拟合C越大越信任训练数据本身的拟合但过大时模型容易跟着训练集的噪声走。实际操作里C从10的-3次方到10的3次方这个范围内变化效果差异会非常明显。gamma是RBF核的参数表达式为exp(-gamma * ||xi - xj||²)。gamma越大核函数越“尖”只有距离很近的样本之间才会有较大核值模型容易被训练样本“钉死”泛化能力下降gamma越小核矩阵越接近常数矩阵模型趋于线性拟合能力不足。这两者的最佳组合在不同数据集上完全不一样手调费时费力这就给了智能优化算法一个发挥空间。1.3 为什么这是个适合优化算法解决的连续问题你可能会有疑问两个参数而已用网格搜索不行吗理论上是可行的10×10的网格就是100次KELM训练数据集不大时几十秒就能跑完。但实际做课题时数据集动辄上千样本一次KELM训练就要计算核矩阵并解一个N×N的线性方程组复杂度接近O(N³)。100次全量训练可能就要十几分钟甚至更久。网格搜索还需要预先确定范围和步长范围给窄了最优值可能在范围外范围给宽了步长又不够精细。更重要的一点在论文里单纯用网格搜索调参写进“方法对比”一节的说服力远不如使用群智能优化算法。这也是这类标题泛滥的根本原因——学术发表场景有对比需求。从工程角度来看C和gamma构成的二维连续搜索空间是相对光滑的元启发式算法往往几十次评估就能逼近一个很不错的区域比网格搜索高效得多。2. 四种算法横向拆解谁收敛快、谁容易早熟2.1 PSO最快找到好解但容易早熟粒子群算法是这里最经典的一个。它的思想很朴素一群粒子在解空间里飞行每个粒子记住自己找到过的最好位置pbest同时整个种群共享全局最好位置gbest。每次迭代粒子根据三个因素更新速度自身惯性、朝个体最优的倾向、朝全局最优的倾向。v w * v c1 * rand(size(x)) .* (pbest - x) c2 * rand(size(x)) .* (gbest - x); x x v;w是惯性权重常见的做法是从0.9线性递减到0.4。c1和c2是加速因子通常取1.5到2之间。PSO的优势是简单、收敛快几轮迭代就能把适应度压下来。缺点是太容易早熟种群在前期一旦被某个较优个体吸引所有粒子会迅速靠拢失去探索能力。优化KELM参数时PSO往往表现不差但容易停在局部最优附近反复跑多次结果波动较大。2.2 WOA包围与螺旋的切换鲸鱼算法模拟座头鲸的气泡网捕食策略核心机制只有三种包围猎物、随机搜索、螺旋气泡网。每次迭代根据随机数p和系数A的值决定走哪条更新路线。a 2 - 2 * t / maxgen; A 2 * a * rand - a; C 2 * rand; if p 0.5 if abs(A) 1 % 包围猎物向当前最优个体靠拢 else % 随机搜索跳出局部区域 end else % 螺旋更新模拟气泡网攻击 endWOA胜在结构简单参数少搜索行为比较随机全局探索能力强于PSO。在KELM参数优化中WOA的表现通常比较平稳不会出现PSO那种非常快的收敛但最终精度往往比PSO好一点。它的弱点是后期收敛速度偏慢需要足够的迭代次数才能充分逼近最优解。2.3 HHO从探索到开发的动态转移哈里斯鹰算法是这里面机制最复杂的一个。它引入了“逃逸能量”E这个概念用E的大小动态决定当前是探索还是开发E0 2 * rand - 1; E 2 * E0 * (1 - t / maxgen);|E|大于等于1时鹰群处于探索阶段随机在广阔区域寻找猎物|E|小于1时转向开发阶段再根据随机数r和|E|的大小细分四种追击策略软包围、硬包围、渐进式快速俯冲。这里的参数多策略多意味着HHO对不同问题的适应能力更强。我在实际使用中的感受是HHO在迭代前期探索得很充分不容易掉进局部最优后期开发阶段又有精细的局部搜索收敛精度往往在三四个算法里名列前茅。代价是代码量更大、单次迭代计算更重同时有些策略分支对边界处理比较敏感需要仔细写。2.4 BOA香味强度引导的全局局部切换蝴蝶优化算法是这里面相对小众的一个。它的灵感来自蝴蝶利用香味定位食物和配偶。每只蝴蝶的位置代表一个候选解适应度转换成“香味强度”f c * I^a其中I是刺激强度即适应度值a和c是常量。每次迭代以概率p决定是全局搜索还是局部搜索r2 rand * rand; if rand p % 全局搜索向全局最优靠近 x_new x r2 * (gbest - x) * f_i; else % 局部搜索在随机两只蝴蝶之间游走 x_new x r2 * (x_j - x_k) * f_i; endBOA的局部搜索机制很特别它没有固定地朝最优靠拢而是在种群内部随机抽取两个个体做差分牵引这样种群多样性保持得不错。实际跑下来BOA在KELM参数优化上的发挥不太稳定有时能超过WOA有时还不如PSO和它的切换概率p以及参数a、c的设置关系很大。2.5 为什么不直接用贝叶斯优化很多人会问Matlab自带的bayesopt不就能调C和gamma吗答案是能但场景不同。贝叶斯优化适合评估次数极少、每次评估很昂贵的场景比如深度网络调参。它需要维护一个代理模型还要选择采集函数代码抽象程度高且老版本Matlab上没有内置支持。智能优化算法实现的逻辑更直观、可控性强迭代过程可以清晰打印和绘制收敛曲线。另外做算法对比实验时把四种算法放在同一套框架下跑只需要替换更新公式写起来非常顺手。这也是这类方案在科研场景里长盛不衰的现实原因。3. Matlab实现关键设计参数编码、适应度与代码骨架3.1 参数编码与搜索边界第一步要想清楚一个问题优化算法的搜索空间是什么。C和gamma的取值范围动辄跨越几个数量级如果直接在原始尺度上搜索比如C从0.001变到1000那么小数值区域几乎不会被充分探索。我强烈建议使用对数编码。% 每个个体的两个维度log10(C), log10(gamma) lb [-3, -3]; % 对应 C0.001, gamma0.001 ub [3, 3]; % 对应 C1000, gamma1000个体位置x是二维向量适应度函数内部用C 10^x(1)、gamma 10^x(2)还原成真实参数。这样做的好处是搜索空间的尺度均匀-3到3相当于C从千分之一到一千倍每一段数量级的范围都被同等对待。实际项目里如果发现最优解总是贴着边界再把边界扩展到[-5, 5]也不迟。3.2 适应度函数与验证集划分适应度函数是整个优化框架的核心。每次算法生成一组候选参数就调用一次KELM训练并返回预测误差作为适应度值优化算法朝误差减小的方向搜索。最容易犯的错误是拿测试集做适应度计算。这会直接导致“测试集泄漏”最终指标看着漂亮一放到新数据上就露馅。正确做法是把原始数据拆成三份trainData % 训练集占60%-70% valData % 验证集占15%-20%用于优化过程中评价参数 testData % 测试集占15%-20%只在所有优化结束后评估一次训练集用来训练KELM验证集用来算适应度测试集只在最终模型确定后才碰。如果数据量很小也可以用K折交叉验证的均值作为适应度代价是每次评估要训练K次KELM运行时间会线性增长。适应度函数代码如下function rmse kelm_objfun(param, Xtrain, Ytrain, Xval, Yval) C 10^param(1); gamma 10^param(2); model kelm_train(Xtrain, Ytrain, C, gamma); Yp kelm_predict(model, Xval); rmse sqrt(mean((Yval - Yp).^2)); end3.3 核矩阵的向量化实现KELM每次训练都要算核矩阵这是整个流程最耗时的部分。新手通常会写双重for循环逐个算核值遇到上千样本时速度惨不忍睹。用向量化写法能快一两个数量级function K kernel_matrix(A, B, gamma) % A: na x d, B: nb x d返回 na x nb 的RBF核矩阵 na size(A,1); nb size(B,1); A2 sum(A.^2, 2); B2 sum(B.^2, 2); D2 max(repmat(A2, 1, nb) repmat(B2., na, 1) - 2 * (A * B.), 0); K exp(-gamma * D2); end这里利用了欧氏距离的展开公式||a-b||² ||a||² ||b||² - 2ab。最后套一层max(...,0)是防止浮点误差产生负的小值导致开方或exp出现NaN。3.4 目标函数统一封装与算法主循环骨架为了保证四种算法对比公平我会把所有算法共用的部分封装成统一接口。无论PSO、WOA、HHO还是BOA主循环结构都是同一个模板% 初始化种群位置x随机生成 % 初始化个体最优 pbest全局最优 gbest for t 1:maxgen for i 1:N % 算法各自的位置更新公式 % 越界处理x max(min(x, ub), lb); % 计算新位置的适应度 % 更新该个体的pbest end % 更新全局gbest % 记录本次迭代的最优适应度 end不同算法只替换中间“位置更新公式”那一块其余完全一样。这样跑出来的结果差异就可以归因于算法本身而不是实现细节不同。我常用的优化设置种群数N30迭代次数maxgen50。这个规模下四种算法都要做30×501500次KELM训练评估。加上数据量不算特别大的话整体运行时间在几分钟到十几分钟完全可接受。如果追求更稳的结果可以把maxgen提到80到100。4. 完整流程与一次实测对比记录4.1 数据预处理与训练流程以某工程软测量场景为例数据集有1200个样本输入特征10维输出1维。我按照0.7比0.15比0.15拆分训练集、验证集、测试集。数据预处理这块输入特征做Z-score标准化输出做min-max归一化都只用训练集的统计量来变换验证集和测试集复用同一组均值和标准差。输出归一化有一个容易被忽视的好处适应度RMSE的值会稳定在一个较小的量级四种算法的收敛曲线更容易对比。最终评估测试集时把预测值反归一化回原始量纲再计算RMSE等指标这样报出来的数字对读者来说更有物理意义。完整流程按顺序走一遍数据加载与划分输入标准化、输出归一化设置优化算法参数种群30、迭代50、边界[-3,3]用验证集RMSE作为适应度跑四种优化算法用每种算法得到的最优C和gamma重新用训练集验证集训练最终KELM在测试集上完成预测并计算指标。4.2 训练与预测核心代码KELM训练和预测部分很短但需要注意几个细节点function model kelm_train(X, Y, C, gamma) n size(X, 1); Omega kernel_matrix(X, X, gamma); Beta (Omega eye(n) / C) \ Y; model.Beta Beta; model.Xtrain X; model.gamma gamma; end function Yp kelm_predict(model, Xtest) Ktest kernel_matrix(Xtest, model.Xtrain, model.gamma); Yp Ktest * model.Beta; end这里用反斜杠运算符求解线性方程组而不要写成inv(Omega eye(n)/C) * Y数值稳定性会差很多。C非常小时Omega eye(n)/C的对角线数值很大会造成病态矩阵反斜杠配合矩阵分解能更好地处理这种情况。还有一个细节Beta的大小是N×1单输出情况预测时是KtestM×N乘BetaN×1完全避开了显式计算逆矩阵。4.3 一次实测结果记录在同一份数据、同一种子、同一评估标准下跑完四种算法后得到的测试集结果如下。这里展示的数值是一次有代表性的运行结果具体数值会随数据集不同而变化但量级和排序趋势有参考意义。算法组合 RMSE MAE R² KELM(默认参数) 21.64 16.32 0.9134 PSO-KELM 8.57 6.12 0.9821 WOA-KELM 7.93 5.66 0.9846 HHO-KELM 7.14 5.08 0.9873 BOA-KELM 8.23 5.91 0.9838默认参数C1、gamma1的KELM在标准化后的输入上表现很一般RMSE高达21.64。经过优化后RMSE都大幅下降。HHO这次拿到了最好的结果WOA紧随其后PSO和BOA基本打平。这个排序不是我提前预设好的而是和算法本身的搜索特性有关HHO的探索阶段足够充分WOA的螺旋机制保证了局部搜索精细度PSO虽然收敛快但在这个解空间中提前落入了某个局部区域。值得注意的是如果我把某个算法的种群数或迭代次数刻意调大它的结果就会变成最好——这种情况在论文汇报里要特别小心解释清楚自己在用相同计算预算做对比否则数据不具有说服力。4.4 收敛曲线解读收敛曲线记录的是每次迭代后全局最优适应度值的变化。四种算法画在一起时有几点很值得观察。PSO的曲线在前期下降非常快通常第5到第10代就把误差压到很低之后曲线趋于平缓甚至完全不再变化。这说明种群已经聚集早熟特征明显。WOA的曲线下降曲线比较平滑前期略慢但一直有持续下降的趋势。HHO的曲线会呈现“阶梯式”下降前期探索阶段适应度波动较大进入开发阶段后突然连续下降几次然后稳定在一个较低水平。BOA曲线通常比较跳跃因为局部搜索阶段引入的随机差分扰动较大适应度会有反复但全局最优值还是能缓慢下降。看自己的运行结果时如果某条曲线长时间没有变化不代表已经找到最优解更可能是种群多样性丧失了。这时候可以检查种群位置是否高度集中或者通过增大种群数来改善。5. 常见问题与调试排查实录5.1 优化后精度反而更差这是最打击人的情况。明明跑了半天优化最后测试集RMSE比默认参数还高。我排查这类问题时第一步就是看适应度函数有没有写错尤其是C和gamma的还原过程。比如预期是10^param(1)结果代码里直接用了param(1)这和初始位置的范围[-3,3]完全不匹配导致KELM的参数变成负数或非法的值KELM训练出来的模型自然一塌糊涂。另外验证集划分不合理也会埋雷。如果验证集只有几十个样本且与训练集分布差异大适应度评估本身就不可靠优化过程就会随机乱跳。最小验证集样本量我建议不低于100个。5.2 收敛曲线不动或震荡收敛曲线完全不动通常有两种原因一是位置更新后全部粒子越界越界处理把它们全部拉回同一个边界点种群实际上变成了同一个点二是适应度函数内部报错导致返回的适应度是一个固定的大值算法误以为“所有位置都一样差”。震荡不断、无法收敛常见于BOA中a和c取值不合适香味强度的尺度差异过大导致更新步长要么太大要么太小。我的建议是对于KELM参数优化这种很简单的问题BOA里的f_i可以直接设成1只保留随机差分项稳定得多。5.3 搜索边界与数据尺度不匹配如果你做完标准化以后最优gamma总是落在10的上边界附近说明当前边界设置不够。反之如果最优gamma集中在0.001附近说明核函数太“尖”了可以检查数据是否做了标准化。核函数依赖距离计算不同特征的尺度差异极大时gamma的最优值会发生明显的偏移。所以数据标准化做没做是排在第一位的检查项。我在一个项目里遇到过类似情况有一列输入特征是温度数值在300到600之间另一列是流量占比在0到1之间。不做标准化的KELMgamma不管怎么调都不理想标准化之后同样的优化流程跑一遍结果立刻好起来了。5.4 KELM训练慢与内存爆炸核矩阵是N×N的稠密矩阵内存占用为8N²字节。N2000时约32MB还能承受N10000时就要800MB开始紧张。对大数据集我建议先用训练集的一部分子样本训练KELM并优化参数等参数确定后再用全部数据训练最终模型。这是近似处理但工程上非常实用。时间开销方面Omega的向量化计算已经是比较快的状态但每一次适应度评估都要做一次矩阵分解。1500次评估乘上N1200时的分解耗时整体运行时间可能在10到20分钟左右。如果嫌慢可以适度降低种群数和迭代次数到20×30或者把适应度改成每5次迭代才重新计算一次牺牲一点精度换速度。5.5 对比实验不公平对比实验最容易翻车的是没有使用相同的计算预算。比如PSO跑了30代HHO跑了100代这根本不是公平对比。正确的做法是固定种群数N和迭代次数maxgen完全一致始位置和随机种子保持一致适应度函数完全一致只在位置更新公式上有所区别。这样才能说明“在相同计算预算下算法A比算法B表现更好”。另一个容易忽视的地方是随机种子。不同算法用不同的随机数种子会引入额外方差。建议在初始化时统一设置rng(1)之类的固定种子保证每次运行结果可复现。论文写作时为了说服力可以把每种算法独立跑30次用平均值和标准差报结果。6. 这套框架还能怎么扩展6.1 把框架迁移到其他模型“优化算法KELM”这套组合只是冰山一角。同样的封装方式把kelm_train和kelm_predict替换成LSSVM、SVR、GRNN甚至简单的BP网络就能搭出“PSO优化SVR预测”“HHO优化LSSVM分类”之类的组合。我见过很多更高阶的做法用HHO同时优化KELM的C、gamma再加上核函数类型选择把离散变量也编码进位置向量里进化到一个范围就取整对应一种核函数。这种连续离散混合编码的方案在一两个项目里确实带来了额外提升但也增加了搜索复杂度并不总是值回票价。6.2 换核、换特征、换数据划分KELM不一定非要用RBF核。多项式核、Sigmoid核在不同数据特征下表现不一。最省事的扩展是让优化算法在几个候选核函数之间也做选择或者分别跑四组实验对比不同核函数下的最优结果。特征选择也可以并进来把位置向量加长到“特征维度2”每一位表示是否选择该特征C和gamma照旧一步完成“特征选择参数优化”的联合搜索。这种方案在特征维度较高时效果很明显但要注意搜索维度的增加会让优化难度明显上升。6.3 一个小技巧坐标轮换搜索最后分享一个我在实际项目中很爱用的小技巧。如果你不想上优化算法或者只是想做一次快速验证可以先固定gamma1只优化C找到合适的C后固定它再优化gamma。这种坐标轮换搜索对二维问题非常有效通常迭代三轮就能得到一个不错的参数组合。它的逻辑在于C和gamma的交互作用没有想象中那么强至少在KELM里这两个参数对模型的影响相对独立。很多同学写代码时默认要同时优化两个参数其实分步做你不仅更容易理解每个参数的作用而且调起来也更快。我在实际使用中还有一个体会不要迷信某一种算法永远最优。每次换数据集、换特征组合最优算法都可能发生变化。真正重要的是把评估流程和代码框架搭好让所有算法在同一条起跑线上公平竞争。谁的适应度最低就用谁的结果这才是最务实的做法。