MATLAB手写SVR：从二次规划到核函数的最小可运行实现

简介面向机器学习初学者和需要实现回归预测的开发者这份PDF围绕支持向量回归展开系统比较多元线性回归、BP神经网络与决策向量机的原理与目标函数差异并重点讲解BP神经网络与决策向量机在优化思路和学习效率上的区别。文档给出了支持向量回归在Matlab中的完整实现从数据导入、主元分析到参数网格寻优再到模型训练与预测的代码示例同时针对径向基核、多项式核、线性核三种核函数分别展示优化前后的对比结果给出了训练集预测值与实际值的绘图比较及均方误差评估并整理了两部分便于迁移使用的代码。内容还对比了BP神经网络与决策向量机在学习效率上的差异指出支持向量机采用数学方法和优化技术区别于神经网络的学习机制文档末尾对不同核函数下的实验结果进行归纳辅助读者在实际项目中合理选择核函数。这份压缩包内含一个PDF文件整体大小约632KB内容精炼、结构清晰特别适合用于算法入门、课程设计或项目选型参考。资料目前已有478人学习下载对想在Matlab环境下快速完成支持向量回归建模、比较核函数效果并用于课程设计或科研实践的读者很有帮助无论是初学还是实战都能从中受益。1. MATLAB 解决 SVR为什么最小可运行代码比调包更有用MATLAB 自带的fitrsvm一句话就能训练支持向量回归可很多人搜「matlab解决svr代码」时真正缺的并不是这个封装函数而是能看懂、能改核、能复现原理的完整实现。SVR 的难点集中在拉格朗日乘子、KKT 条件和 epsilon 不敏感带这三件事上文档里通常一句话带过可一旦落成矩阵、落到quadprog的调用参数上才会发现有大量细节没对上。我把 epsilon-SVR 的最小实现拆开讲一遍先立住二次规划的对偶写法再给出能直接运行并出图的 MATLAB 脚本最后把参数搜索和正确性验证讲透。适合要复现论文、做课程设计以及想把手写 SVR 融入自己优化流程的工程师新手照着步骤能跑通熟手也能从这里找到构造核矩阵和验证收敛的边界细节。2. SVR 原理与 MATLAB 求解路径从 epsilon 不敏感损失到二次规划2.1 原始问题里三个关键项目标、不敏感带、松弛变量先把原问题写成代码注释里常见的形式后面每一行 MATLAB 代码都能对应到这里的某一项min 0.5 * ||w||^2 C * sum(xi_i xi*_i) s.t. y_i - (w * phi(x_i) b) eps xi_i (w * phi(x_i) b) - y_i eps xi*_i xi_i 0, xi*_i 0第一项0.5 * ||w||^2控制模型复杂度防止权重向量范数过大第二项是惩罚项松弛变量xi_i和xi*_i表示样本被允许越过不敏感带的程度。约束条件说的是同一件事预测值w * phi(x) b与真实值y_i的误差超过eps的部分才计入损失误差落在[-eps, eps]内的样本不产生任何损失这就是 epsilon 不敏感损失。C是惩罚系数C越大模型越不肯让样本落到带外也就越容易过拟合。这里有个初学者容易绕进去的点两个方向的违反是不对称出现的一个样本最多只会出现在上带外部或下带外部但两个松弛变量同时大于零的情况在理论推导里是被允许的实际求解时 KKT 条件会自动避免这种冗余。所以代码里不必显式限制xi_i * xi*_i 0二次规划解出来自然满足。2.2 对偶问题与 KKT 条件alpha 减 alpha* 决定稀疏性直接求解原始问题涉及映射phi(x)在特征空间里显式展开代价很高。SVR 的标准做法是转对偶问题让核函数K(x_i, x_j) phi(x_i) * phi(x_j)直接进入计算。对偶形式下决策函数写成f(x) sum_i (alpha_i - alpha*_i) * K(x_i, x) b其中alpha_i、alpha*_i是两组拉格朗日乘子都在[0, C]内且满足sum(alpha_i) sum(alpha*_i)。预测时真正起作用的是beta_i alpha_i - alpha*_i绝大部分beta_i是零只有落在不敏感带边界上或带外的样本对应的beta_i才非零这些样本就是支持向量。这也解释了为什么 SVR 的解是稀疏的噪声点多数落在带内系数为零不参与预测。KKT 条件里有一个对写代码特别重要的结论若某个样本满足0 alpha_i C且alpha*_i 0那么它恰好落在上边界即y_i - f(x_i) eps若0 alpha*_i C且alpha_i 0则落在下边界f(x_i) - y_i eps。这两类样本称为自由支持向量第 3 章求偏置b必须靠它们这也是手写 SVR 时最容易写错的位置。2.3 fitrsvm 与自写二次规划的取舍先给一个直接可用的对比表方便按场景决定要不要放弃调包。对比点fitrsvm自写 quadprog自定义核函数支持函数句柄但调参、调试不透明核函数就是你自己写的矩阵完全可控获取 alpha 乘子需要额外接口且乘子形态是内部约定alpha、alpha* 直接是二次规划解向量增量训练不支持重新训练代价高训练过程拆成核矩阵加 QP可自由组合对求解器的控制黑盒无法查看目标值、对偶间隙所有中间量可见方便排错代码量几行一个完整脚本约几十行fitrsvm对于几万样本、高斯核能覆盖大多数常规回归任务但如果你要自定义核函数、要拿乘子做特征筛选、要把 SVR 嵌进自己的优化循环调包反而难用。还有一个现实限制自写二次规划路径依赖quadprog它属于 MATLAB 优化工具箱安装 MATLAB 时如果没勾选这个组件调用会直接报Undefined function。这是我建议先确认环境的原因和“matlab安装”这一步是同一个问题装完再跑脚本能省很多时间。2.4 构造二次规划系数H、f、约束怎么排进代码对偶问题整理成 MATLABquadprog的标准形式后求解变量为z [alpha; alpha_star]维度是2n。核心组装代码如下这一段是第 3 章完整脚本的地基% 变量排列z(1:n) 为 alphaz(n1:2n) 为 alpha_star n size(K, 1); % 二次项0.5 * z * H * z来自 0.5*(alpha-alpha*)*K*(alpha-alpha*) H [K, -K; -K, K]; % 2n x 2n 对称矩阵 % 一次项f * z来自 -(y*(alpha-alpha*) - eps*sum(alphaalpha*)) f [epsilon - y; epsilon y]; % 2n x 1 % 边界约束0 alpha, alpha* C lb zeros(2 * n, 1); ub C * ones(2 * n, 1); % 等式约束sum(alpha) - sum(alpha*) 0 Aeq [ones(1, n), -ones(1, n)]; beq 0;H的块结构不是随便拼的展开(alpha - alpha*) * K * (alpha - alpha*)会得到四项其中交叉项正好对应右上和左下的-K。f里epsilon - y对应alpha的线性项epsilon y对应alpha*的线性项这两项符号很容易写反写反的后果是预测曲线整体被拉向错误方向。等式约束来自对偶问题对b求偏导后得到的条件缺了它b会无法识别。遇到quadprog报“矩阵必须为正定”的错先对核矩阵做K K 1e-8 * eye(n)的修正重复样本会让K奇异这是最常见原因。3. 用 MATLAB 跑通 SVR 最小实现核矩阵、训练与预测3.1 高斯核矩阵的向量化写法高斯核的公式是K(x_i, x_j) exp(-gamma * ||x_i - x_j||^2)。直接写双重循环在 MATLAB 里效率太低工程上一般用平方距离展开式function K rbf_kernel(X1, X2, gamma) % X1: n1xd 矩阵X2: n2xd 矩阵 % 返回 K: n1xn2K(i,j) exp(-gamma * ||X1(i,:) - X2(j,:)||^2) n1 size(X1, 1); n2 size(X2, 1); x1_sq sum(X1.^2, 2); % n1x1每行向量的平方和 x2_sq sum(X2.^2, 2); % n2x1 % ||a-b||^2 ||a||^2 ||b||^2 - 2*a*b dist2 x1_sq x2_sq - 2 * (X1 * X2); % n1xn2 K exp(-gamma * dist2); end这段代码里x1_sq x2_sq利用 MATLAB 的广播机制生成距离平方矩阵X1 * X2替代循环算点积。dist2是n1 x n2的稠密矩阵所以当训练样本超过一万时内存会快速增长此时需要改成按块计算但在手写 SVR 的场景里样本量通常几千以内这个写法是最稳的。gamma是核宽度参数传入前要确定它控制相似度随距离衰减的速度后面第 4 章会专门讲怎么扫值。3.2 完整训练脚本 svr_demo.m造数据、解 QP、求偏置下面给一个可以直接复制运行的完整脚本数据用带噪声的sin曲线目标是把训练、求解、求b一步走通% svr_demo.m % 手写 epsilon-SVR基于高斯核和 quadprog % 生成带噪声的 sin 数据 rng(42); x (0:0.1:4); y sin(x) 0.15 * randn(size(x)); % 超参数初始化 gamma 1.0; % 高斯核宽度 C 10; % 惩罚系数 epsilon 0.05; % 不敏感带宽度 % 1. 训练核矩阵加 jitter 防止数值奇异 K rbf_kernel(x, x, gamma); K K 1e-8 * eye(length(x)); % 2. 组装二次规划 n length(y); H [K, -K; -K, K]; f [epsilon - y; epsilon y]; lb zeros(2 * n, 1); ub C * ones(2 * n, 1); Aeq [ones(1, n), -ones(1, n)]; beq 0; % 3. 求解二次规划 opts optimoptions(quadprog, Display, off, ... Algorithm, interior-point-convex); z quadprog(H, f, [], [], Aeq, beq, lb, ub, [], opts); % 4. 拆出两组乘子 alpha z(1:n); alpha_star z(n1:end); beta alpha - alpha_star; % 5. 用自由支持向量求偏置 b free1 find(alpha 1e-6 alpha C - 1e-6 alpha_star 1e-6); free2 find(alpha_star 1e-6 alpha_star C - 1e-6 alpha 1e-6); b_candidates []; if ~isempty(free1) % 上边界y_i - f(x_i) eps b_candidates [b_candidates; ... y(free1) - epsilon - K(free1, :) * beta]; end if ~isempty(free2) % 下边界f(x_i) - y_i eps b_candidates [b_candidates; ... y(free2) epsilon - K(free2, :) * beta]; end b mean(b_candidates); % 6. 训练阶段预测与展示 K_all rbf_kernel(x, x, gamma); yhat K_all * beta b; figure; plot(x, y, ko, MarkerSize, 4); hold on; plot(x, yhat, r-, LineWidth, 1.5); xlabel(x); ylabel(y); legend(观测数据, SVR预测, Location, northwest); title(sprintf(手写 SVR: C%g, gamma%g, eps%g, C, gamma, epsilon));脚本里最关键的是第 5 步求b。它不是简单把所有支持向量的误差平均而必须区分上边界自由支持向量和下边界自由支持向量上边界的关系是y_i - f(x_i) eps下边界是f(x_i) - y_i eps。两套式子分别算出候选b再取均值能抵消单侧边界上的数值偏差。如果free1和free2都为空说明epsilon选得太大所有样本都在带内此时b无法唯一确定实际中需要调小epsilon。quadprog里的interior-point-convex是较新 MATLAB 默认提供的凸 QP 算法对H对称半正定的情况支持良好。3.3 预测函数封装新样本点怎么过核矩阵训练完成后对新输入X_new的预测公式是f(X_new) K(X_new, X_train) * beta b。封装成函数避免每次预测重新拼装整个训练过程function yhat svr_predict(X_new, X_train, beta, b, gamma) % X_new: 预测点矩阵X_train: 训练样本矩阵 % beta: 训练得到的 alpha - alpha_starb: 偏置 K rbf_kernel(X_new, X_train, gamma); yhat K * beta b; end调用时只要把训练脚本里的beta、b、gamma传进来比如对x_test (0:0.01:4)预测就会得到平滑曲线。这里要注意预测阶段构造的核矩阵维度是n_test x n_train行是新点、列是训练样本方向不能反。如果反了K * beta的维度对不上MATLAB 会直接报矩阵维度错误。支持向量的可视化也很简单找出abs(beta) 1e-6的下标在scatter里用不同颜色标出即可它们应当集中在曲线转折和噪声较大的区域。4. SVR 参数调优C、gamma、epsilon 的搜索范围与 K 折交叉验证4.1 三个超参数分别控制拟合的哪些行为手写 SVR 的调试热点和调包一样最终都落在C、gamma、epsilon这三个值上。先把各自的行为边界说清楚。参数控制对象取值太小取值太大常用搜索范围C对带外样本的惩罚强度模型欠拟合预测曲线过于平缓过拟合几乎每个点都被逼近2^-5 到 2^15按 2 的幂次扫gamma高斯核的作用半径核函数过于平滑所有点相似度趋同核衰减过快预测曲线剧烈抖动1 / (d * var(X)) 附近的 0.1 到 10 倍epsilon不敏感带宽度支持向量占比升高模型复杂带过宽预测过于平滑甚至退化为常数y 标准差的 5% 到 20%gamma的初始值一般取1 / (特征维度 * 训练集方差)这是一个能保证核矩阵不整体趋近 0 或 1 的经验起点比随手填 0.01 或 100 靠谱得多。epsilon则要结合目标变量y的尺度来定如果y在 0 到 1 之间0.05 是合理的起点如果y是千量级同样取 0.05 会让所有样本都落到带外支持向量占比接近 100%模型彻底失去稀疏性。这三者的调节顺序也有讲究先固定epsilon到一个合理值再扫C和gamma最后回到epsilon细分能显著减少组合次数。4.2 把训练过程封成函数做网格搜索网格搜索的代码不能重复粘贴训练脚本否则嵌套循环里维护局部变量会非常痛苦。先封装一个训练函数返回模型结构体function model svr_train_quadprog(x, y, gamma, C, epsilon) % 返回 model: 含 x、beta、b、gamma、support_idx n length(y); K rbf_kernel(x, x, gamma) 1e-8 * eye(n); H [K, -K; -K, K]; f [epsilon - y; epsilon y]; lb zeros(2*n, 1); ub C * ones(2*n, 1); Aeq [ones(1,n), -ones(1,n)]; beq 0; opts optimoptions(quadprog, Display, off, ... Algorithm, interior-point-convex); z quadprog(H, f, [], [], Aeq, beq, lb, ub, [], opts); alpha z(1:n); alpha_star z(n1:end); beta alpha - alpha_star; free1 find(alpha 1e-6 alpha C - 1e-6 alpha_star 1e-6); free2 find(alpha_star 1e-6 alpha_star C - 1e-6 alpha 1e-6); b_candidates []; if ~isempty(free1) b_candidates [b_candidates; y(free1) - epsilon - K(free1,:) * beta]; end if ~isempty(free2) b_candidates [b_candidates; y(free2) epsilon - K(free2,:) * beta]; end b mean(b_candidates); model.x x; model.beta beta; model.b b; model.gamma gamma; model.support_idx find(abs(beta) 1e-6); end然后写 K 折交叉验证循环。这里用 5 折对每一组参数组合计算验证集 MSE% 网格搜索主脚本 rng(42); x (0:0.1:4); y sin(x) 0.15 * randn(size(x)); % 数据标准化避免 y 量纲影响 epsilon y_mean mean(y); y_std std(y); y_norm (y - y_mean) / y_std; epsilon_base 0.1 * y_std; gammas [0.1, 0.5, 1, 2, 5]; Cs [1, 10, 50, 100]; epsilons [0.5, 1.0, 2.0] * epsilon_base; Kfold 5; n length(x); idx crossvalind(Kfold, n, Kfold); % 需要 Bioinformatics Toolbox best_mse inf; best_params []; for g gammas for Ci Cs for ep epsilons mse_sum 0; for k 1:Kfold test_mask (idx k); train_mask ~test_mask; model svr_train_quadprog(... x(train_mask), y_norm(train_mask), g, Ci, ep); yhat svr_predict(... x(test_mask), model.x, model.beta, model.b, g); yhat yhat * y_std y_mean; mse_sum mse_sum mean((y(test_mask) - yhat).^2); end mse mse_sum / Kfold; if mse best_mse best_mse mse; best_params [g, Ci, ep]; end end end end fprintf(最优参数: gamma%.3f, C%.3f, epsilon%.4f, MSE%.4f\n, ... best_params(1), best_params(2), best_params(3), best_mse);这里把y标准化为均值为 0、标准差为 1 的序列再训练epsilon的候选值直接用标准化之后的尺度设定能保证不敏感带宽度与数据量纲无关。crossvalind来自 Bioinformatics Toolbox如果没装可以用randperm手动分组第 3 章那种 100 行以内的数据完全够用。网格搜索看起来暴力但在几千样本、几百组参数组合的场景下每次求解一个 2n 维 QP 的开销并不高反而是最稳的调参方式。4.3 调参顺序与常见误配置数据缩放和 epsilon 先验网格搜索里最容易翻车的不是循环写错而是数据没缩放。翻 csdn svr讲解时经常看到有人把C取 100、gamma取 0.1然后直接套在自己的数据上结果没收敛或完全欠拟合。原因通常是x的量纲和y的量纲跨度太大核矩阵里dist2的值被某个特征主导gamma的缩放完全失效。手写 SVR 时我会先把x的每一列标准化到均值为 0、方差为 1y也做同样的处理训练完成后再把预测结果反标准化回去。这样gamma的初值可以固定在 0.1 到 10 之间C也可以放心从 1 扫到 1000。还有一个顺序问题不要一开始就同时扫三个参数。先固定epsilon为0.1 * y_std只扫C和gamma找到较优区间后再在0.05 * y_std到0.2 * y_std之间细扫epsilon。原因在于epsilon对支持向量数量的影响是跳跃性的它会直接改变 QP 解的非零乘子分布与C的交互也最复杂把它的搜索放到最后能避免网格组合爆炸。5. 验证 SVR 代码正确性的三个技巧干净数据、不敏感带与对偶间隙5.1 用干净数据和无噪声回归做回归测试参数调完不等于代码正确。第一步验证是用无噪声数据设置C1e6、epsilon1e-6对y sin(x)训练预测曲线应该几乎完全穿过每个训练点。如果此时曲线仍出现明显偏移问题通常在b的符号或f的组装上。第二步验证看支持向量比例正常调优后支持向量占比应在 20% 到 60% 之间如果接近 100%说明epsilon过小或数据噪声过大如果低于 10%多半是epsilon太宽模型已经退化到近乎线性。5.2 epsilon 不敏感带可视化顺带和 BP 网络对比把不敏感带画出来比看任何数值指标都直观% 用训练好的 model 画预测曲线和不敏感带 x_plot (0:0.02:4); yhat svr_predict(x_plot, model.x, model.beta, model.b, model.gamma); y_std std(y); % 若 y 被标准化过这里要映射回原尺度 figure; plot(x, y, ko, MarkerSize, 4); hold on; plot(x_plot, yhat, b-, LineWidth, 1.5); plot(x_plot, yhat epsilon, k--, LineWidth, 1); plot(x_plot, yhat - epsilon, k--, LineWidth, 1); % 标出支持向量 sv_x model.x(model.support_idx); sv_y y(model.support_idx); plot(sv_x, sv_y, ro, MarkerFaceColor, r, MarkerSize, 5); legend(数据, SVR预测, 不敏感带上界, 不敏感带下界, 支持向量);多数训练点应当落在带内支持向量集中在曲线转折处或带边界附近。带内点太少说明模型在追噪声带外点比例过高说明epsilon偏小。画出这条带之后如果再叠一条用 BP 神经网络拟合的曲线能看到一个很典型的差异SVR 的解显著稀疏预测曲线在低密度区域更平稳而 BP 网络在同样数据量下更容易出现局部抖动。这不是说谁绝对好而是验证 SVR 是否保持了它应有的平滑性。5.3 对偶间隙比看损失更可靠的收敛检查quadprog返回后光看fval不足以判断解的质量。一个可靠的做法是计算原始目标函数值和对偶目标值之间的间隙间隙越小说明当前解越接近最优。代码里可以直接利用返回的fval% 训练完成后追加以下检查 K_train rbf_kernel(x, x, gamma); yhat_train K_train * beta b; slack max(abs(y - yhat_train) - epsilon, 0); % 每个样本的松弛量 primal 0.5 * beta * K_train * beta C * sum(slack); dual -fval; % quadprog 最小化的目标等于负对偶目标 rel_gap abs(primal - dual) / (abs(primal) 1e-12); fprintf(相对对偶间隙: %.6f\n, rel_gap);rel_gap小于1e-3通常说明收敛正常如果偏大优先检查核矩阵是否加了足够的 jitter其次检查free1和free2是否为空导致b的估计偏差。把这五段检查代码附到 SVR 脚本末尾每次改动参数后自动跑一遍调试成本会明显下降。本文还有配套的精品资源点击获取