★ FEATURED ARTICLE 简介面向MATLAB用户的多输入多输出回归分析工具包以相关向量机RVM为核心算法适用于工程、科研及经济预测等需要同时拟合多个相关输出变量的回归建模场景尤其适合具备一定编程基础的本科及以上学生、研究生和工程技术人员按需扩展与二次开发。压缩包内共4个文件包括两个m脚本主程序与RBF核函数实现、一个xlsx数据表格以及一个asv自动备份文件整体仅153KB轻量易用m脚本承担模型训练与预测流程数据表提供可直接替换的输入输出样本便于快速验证。目前已有89人学习使用可作为理解RVM原理及多输出回归实现的入门参考实例。代码中已提供完整数据与逐段注释关键函数模块清晰读者可在此基础上调整核参数、替换数据集或拓展为其他机器学习回归方法完成算法对比、模型改进或创新研究从而节省从零搭建环境与处理数据的时间。1. 回归分析总被黑箱模型坑相关向量机凭啥值得用做多输入多输出回归时大部分人第一反应就是上神经网络。可一旦样本量只有几百、特征又有噪声神经网络很容易出现过拟合调参调到怀疑人生。相关向量机Relevance Vector Machine, RVM这时候就是更稳的选择它基于贝叶斯框架训练出的模型不仅给出预测值还能直接给出预测方差相当于告诉你“这次预测置信度有多高”。这个工程性质做设备寿命预测、软测量、经济指标联动分析时非常实用。这份基于MATLAB的RVM多输入多输出回归代码把训练、预测、绘图流程全部串好了数据齐全、带注释适合研究生拿去做课题基线也适合工程师快速验证回归思路。2. 相关向量机的原理与核函数选择为什么它能给置信区间2.1 RVM与SVR的本质差异稀疏性不是唯一亮点支持向量回归SVR靠着结构风险最小化火了二十多年但它的核函数必须满足Mercer条件而且惩罚系数C和epsilon-insensitive损失函数要配对调参。更麻烦的是SVR的预测结果是个点估计想拿到置信区间得额外走bootstrap或者分位数回归代码量直接翻倍。RVM走的是另一条路它假设每个训练样本对应一个权重但给这些权重加了自动相关性判定ARD先验。训练完成后大部分权重被后验概率压到接近零剩下的“相关向量”才是真正决定输出的样本所以模型天然稀疏。这个机制让两个关键事实变得很现实。其一RVM的核函数选择面更宽因为贝叶斯框架下不必受Mercer条件限制其二预测阶段能解析算出方差置信区间是副产品不额外增加计算负担。2.2 核宽度σ对拟合边界的影响这是第一个要调的旋钮这份代码用的是RBF核也就是高斯核。RBF核的数学形式是核函数矩阵里第i行第j列的值等于exp(-||xi - xj||² / σ²)σ是核宽度。σ的取值决定了核函数矩阵的非对角线衰减速度直接影响模型对局部结构的敏感程度。这里是RBFfun.m的实际内容核心就是算两个输入矩阵逐行之间的欧氏距离平方function K RBFfun(X, Y, sigma) % RBFfun计算RBF核矩阵 % 输入 % X - 训练输入矩阵每行一个样本 % Y - 另一个输入矩阵预测时是待测样本训练时与X相同 % sigma - 核宽度需要预选 % 输出 % K - 核矩阵大小为 size(X,1) × size(Y,1) Dist2 repmat(sum(X.^2,2), 1, size(Y,1)) ... repmat(sum(Y.^2,2), size(X,1), 1) - ... 2 * (X * Y); K exp(-Dist2 / sigma^2); end逻辑说明核函数矩阵不是靠双层循环算的而是用距离展开式避免for循环挤占时间这个写法在MATLAB里比pairwise distance之后再exp快得多。Dist2先展开成平方和矩阵再交叉项减掉得到的就是完整的欧氏距离平方。参数说明sigma设置的直觉判断是取训练样本两两距离的中位数左右太小会让核矩阵近似单位阵所有样本都变成孤岛预测值只由最近邻撑起来太大会让核矩阵所有元素趋近1模型退化成平均预测相关向量数量飙升。3. 代码结构与执行流程从数据装载到置信区间全链路拆解3.1 压缩包文件清单哪些是程序、哪些是备份、哪些是数据这份资源解压后的文件结构非常干净明确区分了主程序、备份和原始数据整理成下表文件名类型作用main.m主脚本数据加载、归一化、训练、预测、绘图全流程RBFfun.m函数文件RBF核矩阵计算被main.m调用main.asv自动备份MATLAB编辑器自动保存的历史版本可忽略新建 Microsoft Excel 工作表.xlsx数据文件多输入多输出的原始数据表代码.rar压缩包上述文件的打包载体main.asv是MATLAB编辑器在崩溃退出前自动生成的历史备份不是核心代码把它当作后悔药。如果改坏了main.m可以拿这个文件救回来改名直接用。但正常跑流程时编辑器里只需要打开main.m和RBFfun.m。3.2 数据装载与归一化先看特征尺度再决定用哪种归一化main.m的第一步是用readtable读Excel数据然后立刻做归一化。多输入多输出的数据里各列量纲不一致是非常普遍的现象某个特征如果数值上千而其他特征只有0到1RBF核里的距离计算会被大数值特征完全支配模型退化成一个单变量回归。所以这一步不能省。% 读取数据第一列到倒数第二列是输入特征最后一列是输出 data readtable(新建 Microsoft Excel 工作表.xlsx); X data{:, 1:end-1}; % 输入特征矩阵 Y data{:, end}; % 输出列注意这里暂时只取了一列逻辑说明readtable读Excel后data是table类型花括号索引取的是数值矩阵。常见误区是直接用data(:, 1:end-1)取出来还是table后续运算会疯狂报错。参数说明代码里Y先取最后一列是多输入单输出的写法。如果要跑多输出需要把输出列改成data{:, end-1:end}同时调整后续训练代码里的目标矩阵维度。这一点在4.2节展开。归一化具体用的是mapminmax它把每行数据映射到[-1,1]区间。注意mapminmax是按行操作的而归一化通常是按列做的所以用的时候要转置两次。这一步的代码和踩坑点密不可分% 转置后归一化再转置回来保证每列特征单独归一化 X_norm mapminmax(X, -1, 1); Y_norm mapminmax(Y, -1, 1);逻辑说明mapminmax的默认行为是对行向量做归一化转置后X的每一列是原始的一个特征归一化后再转置回来得到的就是每列特征各自缩放到[-1,1]的矩阵。参数说明上下界选-1和1而不是0和1是因为RBF核的exp指数部分对负值更敏感而且预测层的权重先验用零均值高斯时中心化的输入更好学。预测完成后要反归一化用mapminmax(reverse, Y_pred, Ps)实现其中Ps是训练归一化时保存的参数结构体。3.3 RVM训练与预测迭代稀疏贝叶斯学习的核心流程在MATLAB里实现RVM训练逻辑走向是先构造核矩阵、初始化超参数alpha和beta然后反复更新权重后验统计量、剪枝掉趋于无穷大的alpha、再更新剩余超参数直到迭代收敛或达到最大次数。核心训练循环的主体框架如下% 核函数矩阵训练阶段用训练集自身 K_train RBFfun(X_train_norm, X_train_norm, sigma); N size(K_train, 1); alpha ones(N, 1); % 每个权重对应一个先验精度 beta 1 / var(Y_train_norm); % 噪声方差的倒数 MaxIter 100; % 最大迭代次数防止死循环 for iter 1:MaxIter % 计算后验协方差矩阵Sigma和均值mu Sigma inv(beta * (K_train*K_train) diag(alpha)); mu beta * Sigma * K_train * Y_train_norm; % 根据现有后验均值更新alpha并删除大于阈值的分量 gamma 1 - alpha .* diag(Sigma); alpha_new gamma ./ (mu.^2 eps); keep alpha_new 1e3; % 超过阈值的权重视为不相关剪枝 alpha alpha_new(keep); % 更新噪声精度beta err Y_train_norm - K_train(:, keep) * mu(keep); beta (N - sum(gamma(keep))) / (err*err eps); % 收敛判据alpha变化小于0.01%即停止 if norm(alpha_new(keep) - alpha) / norm(alpha) 1e-4 break; end end逻辑说明这一套迭代本质是证据最大化evidence maximization。alpha是每个权重的高斯先验精度beta是噪声精度后验协方差Sigma是N乘N矩阵beta乘KK加上diag(alpha)求逆。每次迭代用gamma值判断哪个样本对模型还有贡献gamma接近0的样本权重后验方差小而均值小说明它们被先验压得死死直接剪枝。参数说明alpha_new阈值取1e3表示当精度超过1000时认为权重几乎为零把它从相关向量集合里删掉。MaxIter取100是为了防止收敛判据失效时死循环实际数据集一般20轮内就收敛。预测阶段更简单只需要把待测样本的核矩阵算出来再乘权重后验均值K_test RBFfun(X_test_norm, X_train_norm(keep, :), sigma); Y_pred_norm K_test * mu(keep); Y_pred mapminmax(reverse, Y_pred_norm, Ps);逻辑说明K_test是测试样本与训练集中那些被保留的相关向量之间的核矩阵mu(keep)是对应的后验均值权重。预测值是核矩阵和权重的线性组合。参数说明X_train_norm(keep, :)取的是保留样本如果忘记这一步把全量样本塞进核函数计算预测结果会和训练时不一致因为训练模型里其它系数已经被剪掉了。3.4 训练集和测试集划分别让数据泄漏毁了你的评估指标数据划分在RVM里容易翻车的点是随机数种子。如果训练和测试集的划分每次都变评估指标的波动会很大根本没法判断调参有没有效果。main.m里用了固定种子的方式保证可复现rng(42); % 固定随机数种子保证每次运行结果相同 idx randperm(N); train_ratio 0.8; num_train floor(train_ratio * N); train_idx idx(1:num_train); test_idx idx(num_train1:end);逻辑说明randperm是生成随机排列然后用前80%做训练、后20%做测试。rng(42)在每次运行前重置随机数生成器确保别人复现时拿到的划分完全一致评估指标的微小差异只可能来自代码改动而不是随机性。参数说明train_ratio可以根据总样本量调样本只有几十个时建议改成0.9测试集只要保证每个输出维度至少有5个样本量否则误差指标会失真。4. 数据格式与运行参数调优做多输出时要改哪些地方4.1 Excel数据排布要求输入列和输出列的边界怎么界定打开“新建 Microsoft Excel 工作表.xlsx”数据的排布是有约定俗成规矩的。前几列是输入特征最后一列是输出中间不能有文本列、缺失值或空行。特别是空行readtable读进来会变成NaN核函数矩阵一旦碰到NaN整个训练循环全部变成NaN连报错都不太好定位。数据位置内容示例A列至倒数第二列输入特征x1, x2, ..., xp温度、压力、流量等最后一列输出y效率或寿命第1行可留可删如果留表头读取时要用readtable的PreserveVariableNames参数4.2 从单输出改多输出的最小改动方案核心是Y矩阵扩列这份代码默认输出是单列但要做多输入多输出回归本质区别只在于输出矩阵的宽度。改动很小但容易漏推荐按下面两步走。第一步改数据读取把最后一列扩成最后两列或多列。多输出时Y矩阵改成宽度为输出维数的矩阵即可。第二步改归一化mapminmax对Y_norm的处理不变因为mapminmax本来就支持矩阵逐行归一化。第三步改预测与误差计算代码里所有Y的维度从标量向量变成矩阵后计算误差时用mean((Y_true - Y_pred).^2, all)而不是mean((...).^2)否则得到的是每个输出维度分别的MSE评估阶段会混乱。4.3 评估指标怎么算不要只盯着拟合图看代码跑完会画两条曲线一条是训练集的预测对比一条是测试集的预测对比。很多初学者只看图觉得拟合得漂亮就收工了但曲线重合度高不代表泛化性能好。要量化评估直接在脚本末尾补上这段% 计算评估指标R2、RMSE和MAE SS_res sum((Y_true - Y_pred).^2); SS_tot sum((Y_true - mean(Y_true)).^2); R2 1 - SS_res / SS_tot; RMSE sqrt(mean((Y_true - Y_pred).^2)); MAE mean(abs(Y_true - Y_pred)); fprintf(R2 %.4f, RMSE %.4f, MAE %.4f\n, R2, RMSE, MAE);逻辑说明R2是决定系数越接近1说明模型解释了越多的方差RMSE是均方根误差带量纲适合对比不同模型在同一个数据集上的表现MAE是平均绝对误差对离群点不敏感。参数说明如果Y_true是多列的计算SS_res时要用element-wise减法之后再sum否则sum得到的是每列分别的残差平方和而R2的公式在这种情况下不再成立。5. 避坑指南RVM训练与运行中四个高频问题排查5.1 运行报错“未定义与‘double’类型的输入参数相对应的函数‘RBFfun’”现象点击运行main.m时MATLAB提示找不到RBFfun函数。原因脚本和函数文件不在同一个目录下或者当前工作路径没切到解压目录。解决把整个压缩包解压后保持main.m和RBFfun.m的相对位置然后在MATLAB编辑器里打开main.m让当前路径自动切换到编辑器所在目录。如果main.m是在“当前文件夹”窗口双击打开的RBFfun.m没跟着放进来就会触发这个报错。最简单可靠的验证方式是在命令窗口输入which RBFfun如果返回路径说明函数可见如果返回“未找到”就说明路径不对。5.2 训练完成后预测值全是同一个常数现象模型能正常训练但测试集预测的y值几乎全部相同拟合图上预测线是一条平线。原因sigma取值过大核矩阵所有元素趋近于exp(0)1测试样本与所有训练样本的相似度几乎相同预测输出退化为一个常数。另一种可能是归一化后没有反归一化直接把归一化空间里的预测值拿去和原始值作图。解决sigma从距离中位数出发向小的方向调比如先用pdist计算样本距离sigma取median(dist)除以2、除以5这样的梯度往下试探。注意要同步看训练集本身的核矩阵条件数可用cond(K_train)检查条件数超过1e6说明核矩阵接近奇异训练会不稳定。5.3 MATLAB提示“矩阵维度必须一致”现象main.m运行到核矩阵计算或权重更新时报错提示内部矩阵维度必须一致。原因数据读取阶段Y只取了一列但代码后续假设Y是多列或者Excel最后几列有空格被读成了NaN列导致X、Y的列数和主脚本期望不一致。解决先检查Excel数据的列数用size(X)和size(Y)打印出来确认。若输出的确是多列检查输出列的连续性数据写入时若中间夹了纯文本列Excel默认会按多种类型读取肉眼看起来是空白的单元格有可能被识别为空白列这种情况建议用isnan检查把NaN列直接用空值填充或删除。5.4 alpha全部被剪枝模型变成了空模型现象训练结束后相关向量的数量是0或者只剩一两个预测值非常接近于训练集均值。原因sigma过小导致核矩阵接近单位阵每个样本只和自身强相关训练过程认为每个样本的权重都不稳定alpha更新后在贝叶斯框架里被大量压缩。还有一种可能是数据本身没有非线性结构线性模型已经够用RVM把核函数里的非线性成分全部压掉了。解决先做一次最朴素的多元线性回归拿R2做基准。如果线性回归的R2已经达到0.85以上那RVM的非线性增强空间本来就不大这时应该把决策焦点从模型选择转向数据增强或特征构造。如果线性回归很差再把sigma调到距离中位数的量级同时把alpha阈值从1e3降到500让更多权重能活下来。6. 从点预测到区间预测利用RVM预测方差做置信边界这个进阶点才是RVM相对于SVR最值得用的地方。代码里训练完成后其实已经把后验权重协方差Sigma保留下来了。预测阶段待测点的方差不只是噪声方差还包含权重不确定性传递。要拿到置信区间的上界和下界在预测代码后追加这一段% 预测方差 噪声方差 核矩阵对权重协方差的映射 K_test_full RBFfun(X_test_norm, X_train_norm(keep, :), sigma); pred_var 1/beta sum((K_test_full * Sigma(keep, keep)) .* K_test_full, 2); std_dev sqrt(pred_var); Y_lower Y_pred - 1.96 * std_dev; Y_upper Y_pred 1.96 * std_dev;逻辑说明pred_var的第一项是当前噪声精度beta的倒数表示数据噪声本身的方差第二项是把后验协方差Sigma通过核矩阵映射到预测空间的方差相当于权重不确定性在测试点处的传播结果。1.96是95%置信水平下的高斯分布分位数实际工程里做软测量时通常看90%区间对应1.645。参数说明Sigma(keep, keep)这一步要有对应的保留索引否则矩阵维度对不上。绘制区间时把Y_lower和Y_upper用plot画成虚线或者用patch填充半透明色带就能得到一张带置信带的预测图编者按置信区间的质量直接取决于beta估计准不准如果噪声方差被高估区间会宽到失去参考价值。最后一个值得保留的习惯是每次换数据集后我会先跑一遍最小二乘回归作为性能基准无论RVM的效果看起来多好先跟它比一把再下结论。如果RVM在测试集上的RMSE比最小二乘好不了5%以上那就该检查sigma选择和核函数类别而不是直接认定RVM不行。这个习惯我从第一次跑RVM代码时就开始用从那以后每逢回归分析报告都强制走一遍这个流程少走很多弯路希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
