智能优化算法优化核极限学习机预测：原理与Matlab实现

做预测类课题的同学，十有八九都干过这么一件事：换一个智能优化算法，去调一个机器学习模型的参数，然后发一篇“某某算法优化某某模型预测”的论文。哈里斯鹰算法、鲸鱼算法、粒子群算法、蝴蝶算法，再加上核极限学习机，这几个词往一起一放，基本就是一个标准的“优化+KELM预测”组合拳。这篇博文就把这套东西掰开揉碎讲一遍，重点说清楚四类元启发式算法在优化核极限学习机时各自怎么运作、Matlab代码怎么写、实际跑起来会遇到哪些坑，以及怎么设计一次公平的对比实验。内容面向正在做回归预测、软测量、时间序列预测的研究生和工程师，也适合想快速搭建一套可复现预测流水线的朋友。

1. 先把KELM搞清楚：为什么偏偏优化C和γ

1.1 ELM到KELM：一个直观的演化

极限学习机的核心思路其实非常暴力。传统BP神经网络靠梯度下降一层层反传误差，慢且容易陷入局部最优。ELM直接把隐层输入权重和偏置随机生成、保持固定，然后只求解输出层权重，整个训练过程变成一个最小二乘问题，一步到位。

假设训练样本是X（N行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 % 螺旋更新，模拟气泡网攻击 end

WOA胜在结构简单，参数少，搜索行为比较随机，全局探索能力强于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; end

BOA的局部搜索机制很特别，它没有固定地朝最优靠拢，而是在种群内部随机抽取两个个体做差分牵引，这样种群多样性保持得不错。实际跑下来，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]; % 对应 C=0.001, gamma=0.001 ub = [3, 3]; % 对应 C=1000, gamma=1000

个体位置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)); end

3.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

不同算法只替换中间“位置更新公式”那一块，其余完全一样。这样跑出来的结果差异就可以归因于算法本身，而不是实现细节不同。

我常用的优化设置：种群数N=30，迭代次数maxgen=50。这个规模下，四种算法都要做30×50=1500次KELM训练评估。加上数据量不算特别大的话，整体运行时间在几分钟到十几分钟，完全可接受。如果追求更稳的结果，可以把maxgen提到80到100。

4. 完整流程与一次实测对比记录

4.1 数据预处理与训练流程

以某工程软测量场景为例，数据集有1200个样本，输入特征10维，输出1维。我按照0.7比0.15比0.15拆分训练集、验证集、测试集。数据预处理这块，输入特征做Z-score标准化，输出做min-max归一化，都只用训练集的统计量来变换，验证集和测试集复用同一组均值和标准差。

输出归一化有一个容易被忽视的好处：适应度RMSE的值会稳定在一个较小的量级，四种算法的收敛曲线更容易对比。最终评估测试集时，把预测值反归一化回原始量纲，再计算RMSE等指标，这样报出来的数字对读者来说更有物理意义。

完整流程按顺序走一遍：

1. 数据加载与划分；
2. 输入标准化、输出归一化；
3. 设置优化算法参数（种群30、迭代50、边界[-3,3]）；
4. 用验证集RMSE作为适应度跑四种优化算法；
5. 用每种算法得到的最优C和gamma重新用训练集+验证集训练最终KELM；
6. 在测试集上完成预测并计算指标。

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（单输出情况），预测时是Ktest（M×N）乘Beta（N×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

默认参数C=1、gamma=1的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之间。不做标准化的KELM，gamma不管怎么调都不理想，标准化之后同样的优化流程跑一遍，结果立刻好起来了。

5.4 KELM训练慢与内存爆炸

核矩阵是N×N的稠密矩阵，内存占用为8N²字节。N=2000时约32MB，还能承受；N=10000时就要800MB，开始紧张。对大数据集，我建议先用训练集的一部分子样本训练KELM并优化参数，等参数确定后再用全部数据训练最终模型。这是近似处理，但工程上非常实用。

时间开销方面，Omega的向量化计算已经是比较快的状态，但每一次适应度评估都要做一次矩阵分解。1500次评估乘上N=1200时的分解耗时，整体运行时间可能在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 一个小技巧：坐标轮换搜索

最后分享一个我在实际项目中很爱用的小技巧。如果你不想上优化算法，或者只是想做一次快速验证，可以先固定gamma=1，只优化C；找到合适的C后固定它，再优化gamma。这种坐标轮换搜索对二维问题非常有效，通常迭代三轮就能得到一个不错的参数组合。它的逻辑在于C和gamma的交互作用没有想象中那么强，至少在KELM里，这两个参数对模型的影响相对独立。很多同学写代码时默认要同时优化两个参数，其实分步做，你不仅更容易理解每个参数的作用，而且调起来也更快。

我在实际使用中还有一个体会：不要迷信某一种算法永远最优。每次换数据集、换特征组合，最优算法都可能发生变化。真正重要的是把评估流程和代码框架搭好，让所有算法在同一条起跑线上公平竞争。谁的适应度最低，就用谁的结果，这才是最务实的做法。