麻雀搜索算法优化核极限学习机SSA-KELM的MATLAB小样本回归实现

我最早碰这个组合，是处理一个只有180多条样本的工业过程预测任务。传统BP神经网络在这种数据规模下确实不占优势，训练不稳定、超参数多、还容易过拟合。ELM倒是快，但随机生成输入权重这个设计让结果每次都有些差异，复现性很差。后来我把模型换成核极限学习机（KELM），再用麻雀搜索算法（SSA）自动搜索它的核参数和正则化系数，问题才算是干净利落地解决了。这篇文章就把这套MATLAB实现细节、代码逻辑和调参心得完整地整理出来，给正在做小样本回归预测的朋友做个参考。

1. 为什么我会选SSA-KELM解决小样本回归问题

1.1 ELM省去了反向传播，但留下了一个隐性缺陷

极限学习机（ELM）的核心思路是：输入层到隐含层的权重随机生成，隐含层到输出层的权重通过最小二乘一步求解。训练速度极快，不需要迭代优化，这是它最大的卖点。

但问题恰恰出在这个"随机"上。输入权重和偏置一旦随机生成，隐含层输出矩阵H就随之确定了，而H的质量完全取决于抽签运气。同一个数据集，跑十次ELM可能得到十个差异不小的结果。特征映射不稳定，意味着你没法判断模型性能到底是好还是坏，也无法保证下次部署时结果还在可接受范围内。

另一个隐性缺陷是隐含层神经元数量不好选。神经元太少，特征表达能力不够；神经元太多，隐含层输出矩阵容易出现共线性，最小二乘求解时数值稳定性变差。这个超参数在实际调参中非常折磨人，往往要配合交叉验证一起试。

所以ELM只是"快"，并不代表"稳"，在追求稳定复现的工程任务里，这个缺陷是很致命的。

1.2 KELM用核函数补上了随机性这个洞

KELM（核极限学习机）的思路就聪明得多：它不再显式构造隐含层映射函数，而是直接用核函数计算训练样本两两之间的相似度，得到一个N×N的核矩阵Ω，然后用这个核矩阵替代隐含层输出矩阵去做最小二乘回归。

这样做有三个直接好处：

• 消除了随机权重的不确定性，同样的数据和参数，训练结果完全可复现；
• 核矩阵包含了样本间的相似结构信息，对非线性关系的表达能力天然强于手工设置的有限维特征映射；
• 只需要调节两个参数——核参数γ和正则化系数C，调参维度大幅下降。

从数学上看，ELM的隐含层输出矩阵H是N×L的（L是隐含层神经元数），而KELM用核矩阵Ω = HHᵀ代替它，维度变成了N×N。这里有个很关键的推论：核矩阵的维度只跟样本数有关，跟隐含层神经元数量无关。而求逆计算是O(N³)的复杂度，所以KELM特别适合样本量几十到几千的小规模场景；样本量一旦到几万，核矩阵的内存占用和求逆耗时都会变得不可接受。

1.3 SSA把参数搜索变成自动化过程

KELM参数少，但手动调参仍然是个繁琐事。核参数γ决定了样本间相似度随距离衰减的快慢，C是正则化系数，两者都存在一个合适区间，且它们的取值会互相影响。

麻雀搜索算法（SSA）是2020年前后提出的一种群智能优化算法，模拟麻雀的觅食和反捕食行为。把种群划分为发现者、加入者和警戒者三种角色，发现者负责搜索优质食物区域，加入者跟随发现者获取食物，警戒者监视天敌并决定是否发出警报。

我对比过PSO、GA和SSA在KELM超参数搜索上的表现，PSO收敛速度不错但容易早熟，GA全局搜索能力强但收敛偏慢，SSA综合来看在中后期的收敛精度更好，而且MATLAB实现非常简单——一个函数就能写完三种角色的位置更新逻辑，不需要复杂的算子设计。

1.4 这套组合到底适合哪些场景

直接给结论：SSA-KELM最适合20到2000个样本之间的回归任务，尤其是特征维度适中、样本量不足以支撑深度网络训练、但又有一定非线性关系的场景。

一个很典型的情况是做仿真模型替代（代理模型）。比如用仿真软件算一次结果可能要几分钟，你只能算几百个样本点，这时用SSA-KELM拟合出一个近似模型替代仿真计算，误差和效率都比较理想。常见的热词里也提到过"适合小样本仿真数据预测的模型"，高斯过程回归确实也适合这类任务，但GPR的推断复杂度同样跟样本数相关，而且核函数的超参数需要通过极大似然估计迭代求解，在工程部署上不如KELM直接。KELM的训练就是一次矩阵运算，部署和维护成本低很多。

2. MATLAB里不装额外工具箱也能跑通的KELM代码

2.1 核矩阵计算：从双层循环到矩阵运算

KELM中最核心的代码块就是核矩阵的计算。以RBF核为例，公式是：

K(x_i, x_j) = exp(-γ * ||x_i - x_j||²)

最朴素的写法是双层循环，一个样本一个样本地算距离。虽然直观，但在MATLAB里效率太低，完全违背了KELM"一次矩阵运算"的优势。

我用的是矩阵化写法：

function Omega = kernel_matrix(X, gamma) % X: 每一行是一个样本 n = size(X, 1); D2 = sum(X.^2, 2) * ones(1, n) + ones(n, 1) * sum(X.^2, 2)' - 2 * (X * X'); Omega = exp(-gamma * D2); end

这个写法的核心是用欧氏距离的展开式：||a-b||² = ||a||² + ||b||² - 2abᵀ。第一项是每行样本自身的平方和，第二项转置过来，第三项是样本间内积的负两倍。这样做的好处是只用矩阵乘法和广播加法，对几百个样本的规模运算效率非常高。

需要注意，D2里对角线元素理论上应该为零，但浮点运算可能产生微小的负值，导致exp(-gamma * D2)里出现0.9999之类的数值。一般影响不大，但如果核矩阵不是严格半正定的，后续求逆可能出现警告。严谨一点可以在D2后面加一行D2(D2 < 0) = 0;。

2.2 KELM训练的核心逻辑：正则化求逆

KELM的输出权重计算公式是：

β = (Ω + I/C)⁻¹ T

其中T是训练集的输出标签向量，C是正则化系数。这个公式跟岭回归非常像，本质就是在核空间里做带L2正则化的最小二乘回归。

理解这个公式只需要把握两点：

• 为什么加I/C：核矩阵Ω本身可能奇异，加入单位矩阵的缩放项以后，能保证求逆的数值稳定性，同时这个C也起到了控制模型复杂度的作用。C越大，正则化越弱，模型对训练数据的拟合越充分，但过拟合风险也越高；C越小，模型越平滑，泛化性通常更好但可能出现欠拟合。
• 为什么只需要矩阵运算：因为没有迭代过程，不需要梯度下降，一步直接解出解析解。

对应代码：

function [OutputWeight, Omega] = KELM_train(X, Y, C, gamma) Omega = kernel_matrix(X, gamma); n = length(Y); OutputWeight = (Omega + eye(n) / C) \ Y; end

这行代码里用\而不是inv()，也是一个值得说的细节。MATLAB的\运算符会根据矩阵结构自动选择求解算法，对于对称正定阵会采用Cholesky分解，计算效率和数值稳定性都优于直接求逆矩阵再相乘。在小样本场景下两者差异不明显，但遇到核矩阵接近奇异的情况，\会表现得稳健得多。

2.3 KELM预测函数与反归一化操作

训练完成后，对新样本的预测输出为：

f(x) = [K(x, x_1), ..., K(x, x_N)] · β

也就是计算新样本与每一个训练样本的核函数值，组成一行向量，再与输出权重β做内积。对应代码：

function y_pred = KELM_predict(X_train, X_test, OutputWeight, gamma) nt = size(X_test, 1); n = size(X_train, 1); Omega_test = zeros(nt, n); for i = 1:nt d2 = sum((X_test(i,:) - X_train).^2, 2); Omega_test(i, :) = exp(-gamma * d2)'; end y_pred = Omega_test * OutputWeight; end

这里的循环是绕不开的，因为新样本与训练集的距离矩阵结构不属于简单的向量化操作。不过测试样本量通常不会大到离谱，这个循环的耗时可以接受。

预测出来的是归一化后的值，别忘了反归一化：

y_pred_norm = KELM_predict(X_train_norm, X_test_norm, beta, gamma); y_pred = mapminmax('reverse', y_pred_norm', Y_ps);

这里有个常踩的坑：mapminmax在训练时的输入要求是每列一个样本，但KELM的数据矩阵约定是每行一个样本。所以传入mapminmax之前要转置，反归一化之后也记得转置回来。我在第一次写代码时忽略了行列方向，结果预测曲线形状对，但数值全部偏离一个量级，排查了半小时才发现是维度方向问题。

2.4 数据归一化必须注意的两个细节

第一个细节是归一化的范围统一。有人习惯把输入和输出各自做归一化，这没问题，但一定要记录归一化参数（mapminmax返回的X_ps和Y_ps），测试集的归一化必须用训练集算出的参数来变换，不能对测试集单独做归一化。否则预测值在反归一化时会对不上。

第二个细节是如果输出是单列，直接把Y映射到[0,1]区间；如果输出是多列，mapminmax按行处理，需要转置。我的写法是统一用：

[X_norm, X_ps] = mapminmax(X', 0, 1); [Y_norm, Y_ps] = mapminmax(Y', 0, 1);

这样把每一行特征当作一个样本传入，符合mapminmax的习惯。

3. 麻雀搜索算法在MATLAB中的代码级拆解

3.1 种群初始化与适应度函数设计

SSA先初始化一个种群，每个个体代表一组KELM超参数(C, γ)。假如种群数N=30，每个个体就是一个长度为2的向量：

N = 30; % 种群规模 T = 100; % 最大迭代次数 dim = 2; % 待优化的参数维度: [C, gamma] lb = [0.01, 0.01]; % 参数下界 ub = [100, 100]; % 参数上界 X = repmat(lb, N, 1) + rand(N, dim) .* repmat((ub - lb), N, 1);

每个个体的适应度，就是把这组参数代进KELM训练并预测之后的回归误差。一般用均方根误差（RMSE）作为适应度，越小越好：

function fitness = objfun(params, X_train, Y_train, X_val, Y_val) C = params(1); gamma = params(2); [beta, ~] = KELM_train(X_train, Y_train, C, gamma); y_val_pred = KELM_predict(X_train, X_val, beta, gamma); fitness = sqrt(mean((Y_val - y_val_pred).^2)); end

这里还有个取舍：适应度用训练集误差还是验证集误差？我在第5节会专门展开。简单说，用训练集误差速度快、收敛容易，但有说过拟合的风险；用验证集误差更接近真实泛化能力，但代价是SSA每次迭代都要多算一批样本的预测。两者都有人用，关键是一旦选定就不要中途切换。

3.2 三种角色的位置更新逻辑

麻雀算法把种群按适应度排序后分成三类角色。我的实现里，每一步迭代的伪代码逻辑如下：

发现者（适应度排名前20%的个体）负责全局探索，占据食物资源更好的位置，更新策略是：

R2 = rand; for i = 1:floor(N * PD) % PD = 0.2 if R2 < ST % ST = 0.8，说明当前安全 X(i, j) = X(i, j) * exp(-i / (alpha * T)); else % 发现危险，飞向安全区域 X(i, j) = X(i, j) + Q; % Q为正态分布随机数 end end

加入者（剩余80%）跟踪发现者觅食。排名靠前的加入者跟随最优个体，排名靠后的加入者则飞向更远的地方重新搜索：

for i = (floor(N * PD) + 1):N if i > N/2 X(i, :) = Q * exp((X_worst - X(i, :)) / i^2); else X(i, :) = X_best + rand(1, dim) .* abs(X(i, :) - X_best); end end

警戒者（随机选20%个体）负责反捕食。如果这个警戒者本身适应度较差，它会跳到全局最优附近重新搜索；如果它已经是适应度较好的个体，则会就地做一个小扰动，避免陷入局部最优：

for i = 1:floor(N * SD) idx = randi(N); if fitness(idx) > fitness_best X(idx, :) = X_best + beta * abs(X(idx, :) - X_best); else X(idx, :) = X(idx, :) + K_delta * abs(X(idx, :) - X_worst) / ... (fitness(idx) - fitness_worst + eps); end end

理解这套逻辑，可以借助一个比喻。麻雀群里总有几个胆大的"侦察兵"负责探索新地盘，如果发现安全就招呼大家过来饱餐一顿；普通麻雀不主动探索，只跟着前面的人走；但始终有一小撮"哨兵"注意四周，一旦发现危险，立刻逃到安全地点，同时整个族群会因此迅速调整觅食方向。三种角色互相制衡，既有全局探索又有局部开发，所以它的收敛效果在低维连续优化问题上表现得相当均衡。

3.3 边界处理与早熟收敛问题

麻雀个体在飞行迭代中很容易飞出参数边界，比如C变成负数、gamma跑到上千。常用的处理方式有三种：

• 直接设为边界值
• 在边界内随机重新初始化
• 用边界和随机系数把越界的值拉回来

我试过前两种，体验是：直接设边界值会让大量个体堆积在边界上，优化过程变得很钝；随机重新初始化保留了多样性，但会导致收敛曲线偶尔出现跳跃。折中方案是在每次更新后加一个边界检查：

X = max(X, lb); X = min(X, ub);

同时在越界比例超过50%的那一代，随机重置一部分个体到搜索空间内部。这个做法对C和gamma这种取值范围宽、上下界相差几个数量级的参数很实用，因为它能避免种群多样性快速丢失。

3.4 SSA主循环的完整结构

把前面几段串起来，就是整个SSA优化器的骨架：

%% SSA参数设置 N = 30; T = 100; dim = 2; lb = [1e-2, 1e-3]; ub = [1e2, 1e2]; PD = 0.2; SD = 0.2; ST = 0.8; %% 初始化 X = repmat(lb, N, 1) + rand(N, dim) .* repmat((ub - lb), N, 1); fitness = zeros(N, 1); for i = 1:N fitness(i) = objfun(X(i, :), X_train_norm, Y_train_norm, X_train_norm, Y_train_norm); end [best_fit, best_idx] = min(fitness); best_pos = X(best_idx, :); %% 主迭代 for iter = 1:T [sorted_fit, sort_idx] = sort(fitness); X_sorted = X(sort_idx,:); X_best = X_sorted(1, :); X_worst = X_sorted(end, :); fitness_best = sorted_fit(1); fitness_worst = sorted_fit(end); % 发现者更新 ... % 加入者更新 ... % 警戒者更新 ... % 边界处理 ... % 重新计算适应度并更新全局最优 ... end

整个循环的核心就三句话：排序找最优、按角色更新位置、重新算适应度。复杂度很低。真正花费时间在适应度评估，因为每评估一个个体都要训练一次KELM。假设样本数500，每次训练KELM需要计算500×500核矩阵和一次500×500矩阵求逆，MATLAB里大约在10毫秒级别。种群30、迭代100次，总耗时在30秒上下，完全可接受。

4. 实战：在公开数据集上跑通SSA-KELM全流程

4.1 数据准备与实验设置

我用UCI的Energy Efficiency数据集做了个验证实验，取前200条样本模拟小样本场景。特征有8个，预测输出是加热负荷。具体实验设置如下：

• 样本量：200
• 训练集170个，测试集30个（按随机划分，固定随机种子保证可复现）
• 归一化方式：min-max映射到[0,1]
• SSA种群数N=30，最大迭代T=100，搜索范围C∈[0.01,100]，γ∈[0.001,100]
• 适应度函数：训练集RMSE

主脚本中最容易出错的地方，是SSA每次迭代前要把当前最优个体记录下来，而不是用上一轮的最优值更新其他个体。我在初版代码里就是在警戒者更新时误用了上一轮的X_best，导致每轮最优个体不继续迭代更新，收敛曲线出现明显的平台期。

4.2 收敛曲线与优化结果

跑完100次迭代后，SSA的收敛曲线呈现很典型的快速下降—平稳徘徊形态。前20代RMSE从初始的1.8左右迅速降到0.6附近，40代以后下降速度明显减缓，大约在70代左右找到最优解。最优参数落在C≈38.5，γ≈2.7，对应的训练集RMSE约为0.43。

这里有个值得注意的现象：因为C和γ的搜索范围跨了四个数量级，SSA在前期探索阶段会把大量个体推向边界，导致前20代最优适应度变化比较剧烈。等到种群发现边界附近并不是最优区域后，位置分布才逐渐往内部收缩。所以如果收敛曲线前段出现锯齿状跳动，恰恰说明种群多样性还保持得不错，不必慌张。

最优参数下，测试集RMSE约为0.51，R²约为0.93。考虑到这个数据集的真实规律本身存在一定噪声，这个效果对这个样本量级别来说已经比较理想。

4.3 与原始ELM和KELM的精度对比

为了验证SSA到底有没有贡献，我做了对照实验。ELM使用50个隐含层神经元，KELM手工调参，SSA-KELM自动调参。评价指标用测试集RMSE、MAE、R²、MAPE。

SSA-KELM在四项指标上全部最优。这个结果也说明了手动调参的局限性：手动调参时，绝大多数人会倾向于把γ取得偏大，因为核宽度越窄，模型在训练集上的拟合能力越强，看起来效果越好，但测试集误差反而失控。SSA搜索出来的γ取在2到3之间，核宽度适中，相当于自动找到了一个泛化性能更优的平衡点。

4.4 结果可视化与误差分析代码

图是论文和报告里最直观的部分，推荐画三张图：

• 拟合对比图：真实值vs预测值，测试集样本按顺序连线
• 收敛曲线图：迭代次数vs最优适应度
• 误差分布图：预测值-真实值的分布直方图

对应的核心绘图代码：

figure; plot(Y_test, 'ko-', 'LineWidth', 1.2); hold on; plot(y_pred, 'b*-', 'LineWidth', 1.2); legend('真实值', 'SSA-KELM预测值'); xlabel('样本序号'); ylabel('测试集输出'); figure; plot(Convergence_curve, 'r-', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('最优适应度(RMSE)'); figure; histogram(Y_test - y_pred, 15, 'FaceColor', [0.5, 0.8, 0.9]); xlabel('预测误差'); ylabel('频数');

注意第四列误差直方图的中心应当接近零，偏差过大的话说明模型在某一段区域存在系统性偏差，这时可以回头检查训练集和测试集的归一化是否使用了不同参数。

5. 跑了50次之后，关于这套代码的调参与避坑心得

5.1 参数搜索范围应该怎么划

C和γ的搜索范围直接决定优化结果的走向。我跑过几组不同边界的对比实验，感受非常明显：

• C的范围太宽（比如0.01到10000），SSA前期大量个体在无效区域浪费时间，收敛很慢；
• γ的范围太宽（比如0.0001到1000），核宽度跨度太大，极端值会让核矩阵要么几乎退化为单位阵，要么全为1，适应度都很差；
• 比较稳妥的做法是先把C设到[0.01,100]、γ设到[0.001,100]，根据第一次搜索结果将最优参数附近的范围缩小后再重复一轮搜索。

还有一个好用的小技巧：对C和γ这类跨度大的参数，搜索更新时按对数坐标更合理。也就是说初始化时先在对数域采样，更新位置也在对数域进行，最后用10的幂次还原。这样做的好处是让算法在量级变化上步进更均匀，不会让C从50跳到70但γ从0.1跳到50这种不平衡的步长出现。

5.2 适应度用训练集误差还是交叉验证误差

我一开始图省事，适应度直接用训练集RMSE，收敛非常快，20代就能落到很低值。但把最优参数拿去做测试时，发现测试误差波动比较大，有时候R²只有0.85。后来做了对比实验，改用5折交叉验证误差作为适应度，虽然每次迭代的计算量大了将近5倍，但最终模型的测试集R²稳定在0.93以上。

结论是：如果只是快速验证流程，训练集RMSE可以用；如果要追求泛化性能和稳定性，建议切换到K折交叉验证误差。K折的实现不需要自己写循环，用MATLAB的cvpartition就能生成分折索引。

换适应度函数之后SSA的最优参数也会变，这是正常现象。用训练集RMSE求出的C往往偏大，模型偏向过拟合；用交叉验证误差求出的C会小一些，γ也会选得更温和。

5.3 常见报错与对应解决方案

这几类报错我在调试过程中遇到得最多，列出来供参考：

• 矩阵接近奇异的警告：Matrix is close to singular or badly scaled。原因通常是C取值过大，或者核矩阵里某些列高度相关。解决方案是把C的上界调低一些，或者给核矩阵加一个微小的对角扰动+ 1e-8 * eye(n)。
• mapminmax维度不匹配：The number of rows in Y must match...。几乎都是因为忘了传转置矩阵。统一约定用mapminmax(X', 0, 1)即可避开。
• SSA种群个体大量重合：迭代后期适应度曲线不再下降的同时，种群中不同个体的C和γ差异很小。这是早熟收敛的信号。可按5.3节的方法重置一部分个体。
• 预测结果全是一个常数：比如所有预测值都是0.5。这种症状通常是测试集的归一化复用参数时出了问题，测试数据被归一化到了错误的分布上。检查X_ps是否来自训练集。

5.4 后续可以考虑的扩展方向

这套SSA-KELM框架的最大价值在于它的插件化设计。可以做的扩展很多，比如把KELM的RBF核换成多项式核或线性核；把SSA换成其他群智能优化算法做对比；也可以把单输出扩展到多输出KELM。另外，SSA的适应度函数也可以改造成包含时间复杂度的多目标形式，在预测精度和训练效率之间做权衡。

我在使用中发现，SSA-KELM在Matlab里跑通的完整链路其实不到200行代码，很多时间其实花在数据清洗和结果分析上。如果把这个框架封装成一个函数，输入训练X、训练Y、测试X和参数范围，直接返回最优参数和预测结果，后面遇到新数据集就能一键复用。这也是我建议想用这套方法的人最值得先做的一件事。