多核LSSVM与PSO参数优化：MATLAB回归预测实战指南

简介：MATLAB实现的多核最小二乘支持向量机（LSSVM）优化代码包，面向机器学习和数据挖掘方向的研究者与学生，用于处理复杂非线性数据的分类与回归建模。代码引入粒子群优化（PSO）对多核LSSVM的正则化参数和核函数参数进行自动寻优，从而提升模型泛化能力与预测精度。压缩包共7个文件，以4个mat数据文件和3个m脚本文件组成，整体仅10KB；mat文件存放训练/测试样本及对应标签，m文件则涵盖PSO参数优化、多核LSSVM模型构建以及自适应改进函数等核心模块。目前已有877人学习下载。资源在原单核LSSVM基础上扩展到多核版本，并结合PSO自动调参，既适合入门者理解多核学习与智能优化算法结合的原理，也可为进阶研究提供可运行的轻量级代码参考，便于进行模型对比与二次开发，是一份实用价值较高的算法学习资料。

1. 多核LSSVM与PSO-LSSVM到底在解决什么问题

做回归预测建模时，经常遇到一个很尴尬的情况：LSSVM的单一RBF核在局部区间拟合得很好，换一个区间就出现系统性偏差；把核换成多项式，整体趋势对了，局部细节又丢掉。与其来回换核函数，不如把多个核组合起来。多核LSSVM让核矩阵变成多个基核的加权和，但代价是待调参数从两三个变成五六个，网格搜索直接爆炸。PSO-LSSVM的做法就是用粒子群算法去搜这一组参数，包括LSSVM的正则化参数γ、基核宽度σ和各核权重w。下面这套代码可以直接在MATLAB里运行：不依赖matlab优化工具箱，手写PSO主循环，用3折交叉验证做适应度评估，最后输出预测结果和最优核权重。适合需要在回归、软测量这类任务里快速搭一套可解释调参模型的人。

2. 多核LSSVM的建模原理与待优化参数

2.1 从SVM到LSSVM：二次规划换成线性方程组

LSSVM（Least Squares Support Vector Machine）和标准SVM的核心差别在于误差项的写法。标准SVM用ε-不敏感带加松弛变量，求解的是带不等式约束的二次规划；LSSVM直接把误差写成等式约束：

min J(w, e) = 1/2 * ||w||^2 + γ/2 * Σ e_i^2 s.t. y_i = w^T * φ(x_i) + b + e_i, i = 1,...,n

每个等式约束配一个拉格朗日乘子α_i，KKT条件整理后得到一个 (n+1) 维线性方程组：

[ 0 1^T ] * [ b ] [ 0 ] [ 1 Ω + γ^{-1} I] [ α ] = [ y ]

其中 Ω(i,j) = K(x_i, x_j)，1 是 n 维全 1 向量。这个方程组直接用 MATLAB 的反斜杠\求解即可，不需要任何优化迭代。这里要注意 γ 的物理含义：γ 出现在核矩阵对角线的偏移项1/γ上，γ 越大，对角线偏移越小，模型对训练样本的拟合越激进；γ 越小，模型越平滑。这也是后面对 γ 做对数编码的原因——它实际起作用的是数量级。

LSSVM 没有稀疏解，几乎每个训练样本的 α_i 都不为零，这是它和标准 SVM 的最大差异。对于回归预测这类不追求支持向量稀疏性的任务，这个代价完全可接受，换来的是训练速度快一个量级。预测公式为：

y_new = Σ α_i * K(x_i, x_new) + b

在 MATLAB 里，这等价于测试核矩阵乘上 α 再加 b。理解了这条主线，多核 LSSVM 就只剩一个改动：把 K 从单一核换成多个核的加权和。

2.2 多核LSSVM的核组合机理

多核 LSSVM 的做法是把核矩阵写成若干基核的加权和：

K(i,j) = Σ w_k * K_k(x_i, x_j), w_k ≥ 0, Σ w_k = 1

权重非负且和为 1 不是随便定的。核函数要满足 Mercer 条件才能保证对应的核矩阵半正定，而非负系数的线性组合能保持这一性质。如果你把某个权重设成负数，核矩阵可能出现负特征值，LSSVM 的线性方程组求解就不稳定了。

选择什么样的基核跟数据分布相关。我一般会固定三到五个差异足够大的基核，而不是堆十几个。差异大意味着它们的特征映射不在同一个方向上，组合起来才有意义。常见组合是 RBF 核 + 多项式核 + 线性核：RBF 负责局部细节，多项式负责全局趋势，线性核在特征维度不高时作为保底项。如果三个核本身高度相似，多核组合的收益就很有限。

多核真正带来的问题是参数维度上升。RBF 核有宽度 σ，多项式核有次数和偏移，每个核还有权重 w。再加上 LSSVM 本身的 γ，参数空间很容易到 5 维以上。网格搜索在这种维度下不可行，所以需要 PSO 这种不依赖梯度、对目标函数解析形式没有要求的实数优化算法。

2.3 PSO-LSSVM的待优化参数与选型

把整套 PSO-LSSVM 的待优化参数列出来，方便后面写代码时一一对应：

PSO 的选型理由很直接：目标函数是交叉验证 RMSE，它关于 γ 和 σ 不可导，而且存在大量局部极小。粒子群的全局搜索能力不依赖初始点，实现也就几十行代码。相比 matlab 自带的优化工具箱，手写 PSO 能把每一次迭代里核矩阵的计算复用做到极致——这一点在下一章会体现出来。

3. 多核LSSVM的MATLAB核心代码与核矩阵预计算

3.1 数据准备与基核矩阵预计算

先造一个带非线性、非平稳特性的回归数据集，方便验证多核的效果。这里用 150 个样本，输入是一维的，输出包含周期项、趋势项和噪声：

clear; clc; rng(42); N = 150; X = linspace(-4, 4, N)'; Y = 0.8 * sin(2*X) ./ (0.6 + abs(X)) + 0.25*X + 0.1*randn(N, 1); % 训练/测试 60/40 划分 idx_train = randperm(N, round(0.6*N)); idx_test = setdiff(1:N, idx_train); Xtr = X(idx_train); Ytr = Y(idx_train); Xte = X(idx_test); Yte = Y(idx_test);

数据生成里sin(2*X)./(0.6+abs(X))保证了 X=0 处不会出现除零，同时周期和趋势叠加的结构会让单核模型“偏科”：RBF 核能抓住局部周期，但趋势外推弱；多项式核趋势好，但局部拟合不细腻。这正是多核发挥作用的场景。

接下来预计算三个基核矩阵。这一步是整个方案里最关键的提速点：PSO 每次评估适应度时，只需要对三个固定好的矩阵做加权求和，不需要重新算核函数。

% 训练集核矩阵 K_rbf = exp(-pdist2(Xtr, Xtr).^2 / (2*1.0^2)); % RBF核，宽度固定为1.0 K_poly = (Xtr * Xtr' + 1).^2; % 多项式核，次数2，偏移1 K_lin = Xtr * Xtr'; % 线性核 % 测试集对训练集的核矩阵，用于最终预测 K_rbf_te = exp(-pdist2(Xte, Xtr).^2 / (2*1.0^2)); K_poly_te = (Xte * Xtr' + 1).^2; K_lin_te = Xte * Xtr'; K_cell = {K_rbf, K_poly, K_lin}; K_cell_te = {K_rbf_te, K_poly_te, K_lin_te};

这里把 RBF 宽度固定为 1.0 是有意的。第一版先让 PSO 只搜 γ 和核权重，核宽度后面再放开。pdist2 是 MATLAB 基础函数，不需要任何工具箱。三个核矩阵的尺寸都是 90×90（训练集约 90 个样本），测试核矩阵是 60×90，内存开销完全可以忽略。

3.2 LSSVM训练与预测的MATLAB函数

核心求解器写在两个局部函数里。MATLAB R2016b 及以后版本允许把函数直接放在脚本末尾，不需要单独建文件：

function [alpha, b] = lssvm_train(K, y, gamma) n = size(K, 1); A = [0, ones(1, n); ones(n, 1), K + eye(n)/gamma]; sol = A \ [0; y(:)]; b = sol(1); alpha = sol(2:end); end function yp = lssvm_predict(K_pred, alpha, b) yp = K_pred * alpha + b; end

A \ [0; y]就是求解那个 (n+1) 维线性方程组。矩阵 A 的第一行第一列是 0，对应约束 Σα_i = 0；右下角块K + eye(n)/gamma就是前面公式里的Ω + γ^{-1}I。gamma 越大，eye(n)/gamma越小，核矩阵越接近原始 K，模型越容易过拟合。预测函数里K_pred是测试样本与训练样本的核矩阵，尺寸 m×n，乘上 α 再加上偏置 b 得到 m×1 预测值。

注意y(:)先把标签拉成列向量，避免行向量输入导致广播维度出错。这是 MATLAB 里很常见的坑——如果你从表格或元胞数组里取数据，y 可能是行向量。

3.3 核权重的softmax解码，避免边界约束

直接优化 w1、w2、w3 有个麻烦：PSO 是无约束优化，粒子的位置更新后很可能跑出 [0,1] 范围，导致权重和不为 1。虽然可以用“钳制+归一化”处理，但边界上的梯度信息很别扭。

我一般用 softmax 把无约束的实数 t 映射成有效权重：

function w = softmax_weights(t) z = [t(:)', 0]; % 第3个核的t固定为0作为基准 ez = exp(z - max(z)); % 减max防溢出 w = ez / sum(ez); end

这个函数接收一个 2 维向量 t，输出 3 维权重，且恒满足 w_k ≥ 0、Σw_k = 1。PSO 搜索的其实是 t1、t2，而不是权重本身。理由是 softmax 的输入落在实数轴任意位置都能映射到合法权重，粒子在任何区域的解都可用。

3.4 待搜索参数与范围设置

这一版 PSO 搜索 3 个维度，每个维度的意义和边界如下：

t 的范围给 [-10, 10] 看似很宽，但 softmax 输出对输入差值不敏感，±10 已经能把权重压缩到接近 0 或接近 1，不需要精确搜索边界。log10(γ) 范围对应 γ ∈ [0.001, 1000]，基本覆盖了回归任务里“欠平滑到过拟合”的完整区间。

4. PSO寻优多核LSSVM的MATLAB主程序与交叉验证

4.1 粒子编码与PSO参数范围

粒子向量为pos = [log10(γ), t1, t2]，维度 dim = 3。PSO 的固定参数设置如下：

nP = 30; nIter = 30; dim = 3; lb = [-3, -10, -10]; ub = [3, 10, 10]; pos = rand(nP, dim) .* (ub - lb) + lb; % 均匀初始化 vel = zeros(nP, dim); % 初始速度为0

惯性权重在这里需要随迭代线性衰减，从 0.9 降到 0.4。前期大惯性权重保证全局探索，后期小权重收敛到局部精细搜索。c1 和 c2 都取 1.5，这是 PSO 经验参数里最常用的默认值，对 3~6 维问题表现稳定。

在进行 PSO 迭代之前，要先固定交叉验证的折划分。否则每次调用适应度函数都随机重划分，同一个粒子两次评估的 RMSE 会因数据划分不同而抖动，PSO 会误把噪声当作适应度变化，收敛明显变差。

4.2 适应度函数：3折交叉验证的RMSE

我用 3 折交叉验证计算 RMSE 作为适应度。折数太少评估噪声大，折数太多每次要训练 5 次 LSSVM，30 个粒子迭代 30 代就是 4500 次训练，没必要。3 折是效率与稳定的平衡点。

fold = mod(randperm(length(Ytr)), 3) + 1; % 固定划分，保证可复现 function rmse_cv = fitness_A(x, Ytr, K_cell, fold) gamma = 10^x(1); w = softmax_weights(x(2:3)); % 核矩阵加权组合 K = zeros(size(K_cell{1})); for k = 1:3 K = K + w(k) * K_cell{k}; end sse = 0; for f = 1:3 ite = (fold == f); % 测试折 itr = ~ite; % 训练折 [alpha, b] = lssvm_train(K(itr, itr), Ytr(itr), gamma); yp = K(ite, itr) * alpha + b; sse = sse + sum((Ytr(ite) - yp).^2); end rmse_cv = sqrt(sse / length(Ytr)); end

核矩阵组合只做三次矩阵数乘再加起来，共用size(K_cell{1})的零矩阵。三个基核预计算后，适应度评估里最重的运算是那个 60×60 线性方程组的求解，速度很快。归一化 RMSE 用所有样本的总误差除以总样本数，而不是把每折的 RMSE 做算术平均，这样折大小不一致时评估更稳。

提示：如果适应度出现 NaN，先检查 gamma 是否解析为 Inf 或 0。log10 编码一般不会，但粒子速度过大时位置可能溢出，边界钳制要加上。

4.3 PSO主循环与边界钳制

完整的 PSO 主循环如下。每代先评估所有粒子的适应度，更新个体最优 pbest 和全局最优 gbest，再用标准速度更新公式推进粒子：

pbest = pos; pbest_fit = inf(nP, 1); gbest = pos(1, :); gbest_fit = inf; for it = 1:nIter w_i = 0.9 - (0.9 - 0.4) * (it - 1) / (nIter - 1); % 惯性权重线性衰减 for i = 1:nP fit = fitness_A(pos(i, :), Ytr, K_cell, fold); if fit < pbest_fit(i) pbest_fit(i) = fit; pbest(i, :) = pos(i, :); end if fit < gbest_fit gbest_fit = fit; gbest = pos(i, :); end end r1 = rand(nP, dim); r2 = rand(nP, dim); vel = w_i * vel + 1.5 * r1 .* (pbest - pos) + 1.5 * r2 .* (gbest - pos); pos = pos + vel; pos = min(max(pos, lb), ub); % 边界钳制 fprintf('iter %d, best RMSE_cv = %.5f\n', it, gbest_fit); end

这里速度更新里1.5 * r1 .* (pbest - pos)是认知项，1.5 * r2 .* (gbest - pos)是社会项。两者分别把粒子拉向自身历史最优和全局最优。注意速度vel没有显式上限，所以位置钳制到[lb, ub]是必要的，否则粒子可能直接飞到 log10(γ) = 20 这种无意义区域，softmax 输入过大时 exp 溢出。

迭代 30 代后，gbest 里保存的就是最优的[log10(γ), t1, t2]。可以在迭代过程中把gbest_fit画出来观察收敛情况：正常情况下前 5 代 RMSE 快速下降，后 20 代缓慢波动，波动幅度逐步收窄。

4.4 用gbest重建多核LSSVM并对比单核模型

PSO 得到的最优参数还需要重建模型，在留出的测试集上做一次最终评估。同时用相同 γ 跑一个纯 RBF 核的 LSSVM 作为 baseline，判断多核到底带来多少收益：

gamma_best = 10^gbest(1); w_best = softmax_weights(gbest(2:3)); % 重建最优多核矩阵 K_best = zeros(size(K_cell{1})); K_best_te = zeros(size(K_cell_te{1})); for k = 1:3 K_best = K_best + w_best(k) * K_cell{k}; K_best_te = K_best_te + w_best(k) * K_cell_te{k}; end % 多核模型 [alpha, b] = lssvm_train(K_best, Ytr, gamma_best); Yp = lssvm_predict(K_best_te, alpha, b); rmse_mk = sqrt(mean((Yte - Yp).^2)); % 单RBF核baseline [alpha_rbf, b_rbf] = lssvm_train(K_cell{1}, Ytr, gamma_best); Yp_rbf = lssvm_predict(K_cell_te{1}, alpha_rbf, b_rbf); rmse_rbf = sqrt(mean((Yte - Yp_rbf).^2)); fprintf('多核LSSVM test RMSE = %.5f, weights = [%.3f %.3f %.3f]\n', ... rmse_mk, w_best); fprintf('单RBF核LSSVM test RMSE = %.5f\n', rmse_rbf);

这里有个细节很容易被忽略：baseline 用的 γ 和 PSO 搜出来的最优 γ 相同。这样对比的是“多核权重”带来的增益，而不是“调参带来的增益”。如果你发现单 RBF 核的 RMSE 反而更低，先不要急着否定多核，可能是三个基核的宽度设置不合理，RBF 核宽度 1.0 在这个数据上恰好太合适了。这种场景下真正有价值的做法是把 σ 也放进粒子搜索，下一章展开。

5. PSO-LSSVM的验证与进阶：核宽度搜索与稳健求解

5.1 先做两组对比，验证核权重起了作用

多核权重容易给人“参数越多越好”的错觉。要验证权重真正有用，建议做两组对照：第一组是“等权重 + 最优 γ”，第二组是“PSO 搜出的权重 + 最优 γ”。如果两者测试 RMSE 差距很小，说明三个基核本身冗余，多核收益有限；如果差距明显，说明 PSO 确实找到了更合理的核组合方向。我在实际项目中更关注权重的分布形态：某个核权重被压到 0.05 以下，基本可以判断这个核贡献不足，下一步应该换别的基核，而不是继续加权重。

5.2 核矩阵数值问题与稳健求解

LSSVM 的线性方程组在 gamma 偏大时容易病态，尤其当样本存在重复或高度相关特征时。A\b内部会做列主元消去，大多数情况没问题，但如果求解结果出现符号抖动或预测值异常大，可以把核矩阵显式对称化：

H = K + eye(n) / gamma; H = (H + H') / 2; % 消除数值不对称 sol = [0, ones(1,n); ones(n,1), H] \ [0; y(:)];

还可以观察核矩阵的条件数。条件数超过 1e12 时，建议先对输入特征做标准化，再考虑调小 gamma。PSO 里对 gamma 做对数编码本身就能缓解这个问题，因为 log10 范围限制了粒子不可能直接取到 1e8 这种危险值。

5.3 让PSO连同核宽度一起搜索

固定核宽度只是第一版方案。当核宽度本身需要调优时，把粒子扩展到 5 维：[log10(γ), log10(σ1), log10(σ2), t1, t2]。σ1 是 RBF 核宽度，σ2 是多项式核的偏移系数。这时核矩阵每次都要按 σ 重新计算，之前预计算的方案失效，需要用缓存避免重复计算：

cache = containers.Map('KeyType', 'char', 'ValueType', 'any'); key = sprintf('s%.4f', sigma); if cache.isKey(key) K_rbf = cache(key); else K_rbf = exp(-pdist2(Xtr, Xtr).^2 / (2*sigma^2)); cache(key) = K_rbf; end

核矩阵缓存让同一个 σ 在多次适应度评估中只算一次。粒子数 30、迭代 30 代时，缓存命中率通常能到 60% 以上。σ 的搜索范围我一般设为 log10 域 [-1, 1]，对应 σ ∈ [0.1, 10]，覆盖从“几乎线性”到“高度局部”的完整过渡；这个范围再放宽后，你在测试集上看到的 RMSE 下降会比单纯调核权重更大一些。

本文还有配套的精品资源，点击获取