
简介这份资源面向需要掌握多输入多输出回归建模与模型可解释性的机器学习学习者和工程人员提供一套基于MATLAB的完整实现方案。内容围绕BP神经网络展开覆盖回归预测、SHAP可解释分析以及新数据预测三大环节配套Excel格式的多输入多输出数据集可直接替换数据复现实验。压缩包共8个文件包含4个m脚本、3个xlsx数据表和1个txt说明文件整体约55KB脚本分别承担主流程回归、SHAP值计算、新样本预测与核心函数封装等职责数据表则提供训练与待预测样本。已有92人学习关注。读者可借此理解BP网络在多输出任务中的搭建与调参思路掌握SHAP方法对特征贡献的量化解释流程并学会将训练好的模型迁移到新数据上完成预测适合作为课程设计、论文实验或工程原型的参考模板。1. 从一张“黑箱”预测表说起BP神经网络多输入多输出回归到底在解决什么你手里有一批实验数据输入是 6 个工艺参数输出是 3 个性能指标领导要你“建个模型既能预测新样本又能说清楚哪个输入影响最大”。这时候单输出模型要训 3 次每次还得单独调参改一个输入维度就得全部重跑。BP 神经网络多输入多输出回归就是冲着这个场景来的一个网络同时吐出多个目标值共享隐含层特征训练一次搞定多目标。但纯 BP 的预测结果没人敢信因为它是黑箱——SHAP 可解释分析就是那把撬开黑箱的螺丝刀用博弈论里的 Shapley 值给每个输入特征分配贡献度让你能对着图说“第 3 个输入对第 2 个输出的影响占了 40%”。MATLAB 完整源码和数据意味着你不用从零搭轮子改改数据接口就能跑自己的项目。这套组合适合做实验数据回归、工艺参数优化、多指标预测的工程师尤其是样本量在几百到几千条、输入输出维度不超过 20 的场景。下面我把从数据组织到 SHAP 解释再到新数据预测的完整链路拆开讲中间踩过的坑一并奉上。2. 多输入多输出 BP 网络的数据组织与网络搭建2.1 输入输出矩阵怎么摆MATLAB 里的行列约定BP 网络在 MATLAB 里最容易被数据维度搞翻车。feedforwardnet或fitnet默认把每一列当作一个样本每一行当作一个特征。也就是说如果你有 500 个样本、6 个输入、3 个输出输入矩阵X应该是 6×500输出矩阵Y应该是 3×500。很多人从 Excel 读进来是 500×6直接丢进去训练结果网络把 500 个特征、6 个样本拿去学训练误差看着降了预测全是垃圾。这是血泪经验里排第一的坑。正确的数据组织方式如下% 假设 rawData 是 500×9 的矩阵前6列输入后3列输出 rawData readmatrix(data.xlsx); % 500×9 X rawData(:, 1:6); % 转置为 6×500每列一个样本 Y rawData(:, 7:9); % 转置为 3×500每列一个样本 % 检查维度 fprintf(输入维度: %d×%d\n, size(X,1), size(X,2)); fprintf(输出维度: %d×%d\n, size(Y,1), size(Y,2));这段代码的核心就两步读数据、转置。readmatrix是 MATLAB R2019a 之后推荐的表格读取函数比xlsread干净。转置之后size(X,1)是特征数size(X,2)是样本数后面所有操作都围绕这个约定。如果你用的是.mat文件直接load进来后检查变量名别假设它一定叫data。提示转置后一定用size打印确认别凭感觉。我见过有人转置了两次等于没转训练了半小时才发现。2.2 网络结构选型隐含层节点数不是越多越好多输入多输出 BP 网络的结构设计核心就三个决策几个隐含层、每层多少节点、用什么训练函数。对于输入输出维度都在 20 以内、样本量几百到几千的问题一个隐含层足够。理论上有万能逼近定理撑着两个隐含层只在函数复杂度极高时才需要而且更容易过拟合。隐含层节点数的经验公式有好几个我一般用这个起步nInput size(X, 1); % 输入维度比如 6 nOutput size(Y, 1); % 输出维度比如 3 nSample size(X, 2); % 样本数比如 500 % 经验公式sqrt(输入输出) 调节项 nHidden round(sqrt(nInput nOutput)) 5; % 约 8 % 或者用 2*输入1 起步 nHidden_alt 2 * nInput 1; % 13 % 搭建网络 net feedforwardnet(nHidden, trainlm); net.trainParam.epochs 1000; net.trainParam.goal 1e-5; net.trainParam.lr 0.01; net.trainParam.showWindow false; % 批量跑的时候关掉窗口feedforwardnet的第一个参数是隐含层节点数第二个是训练函数。trainlm是 Levenberg-Marquardt 算法收敛快适合中小规模网络但内存占用比trainscg高。如果样本超过一万条换trainscg更稳。epochs设 1000 是上限实际训练中如果验证集误差连续 6 次不降trainlm会自动早停这个默认参数是net.trainParam.max_fail 6。goal设 1e-5 是目标误差别设太小否则容易过拟合。节点数怎么定我的做法是从sqrt(nInputnOutput)5开始跑三次不同随机种子看验证集 MSE 的均值和方差。如果方差大说明节点数偏多减 2 到 3 个再试。如果均值高加 2 个。这个过程一般迭代 3 到 4 轮就能找到稳定区间。别用遗传算法或粒子群去优化节点数对于这个规模的问题手动试比自动搜索快。2.3 数据划分与归一化别让量纲差异毁了训练多输入场景下不同输入的量纲可能差几个数量级。比如温度是 200 到 800压力是 0.1 到 0.5如果不归一化梯度下降会被大量纲特征主导小量纲特征几乎不更新。MATLAB 的feedforwardnet默认在训练前自动做mapminmax归一化把数据映射到 [-1, 1]训练后再反归一化输出。但这个自动处理有个坑它是在train函数内部做的你拿到的net对象里归一化参数存在net.inputs{1}.processSettings里新数据预测时必须手动调用同样的归一化参数否则预测结果完全不对。% 手动划分训练/验证/测试集 net.divideParam.trainRatio 0.7; net.divideParam.valRatio 0.15; net.divideParam.testRatio 0.15; % 训练 [net, tr] train(net, X, Y); % 训练集预测 Y_train_pred net(X(:, tr.trainInd)); % 测试集预测 Y_test_pred net(X(:, tr.testInd)); % 计算测试集 MSE mse_test perform(net, Y(:, tr.testInd), Y_test_pred); fprintf(测试集 MSE: %.6f\n, mse_test);divideParam的三个比例加起来必须是 1。tr.trainInd、tr.valInd、tr.testInd是训练完成后tr结构体里的索引直接拿来切数据最可靠。perform函数自动处理了归一化和反归一化算出来的 MSE 是原始量纲下的。如果你自己手算mean((Y_pred - Y_true).^2)记得先反归一化否则数值对不上。注意train函数每次调用会重新随机划分数据想复现结果就在train之前设rng(42)固定种子。3. SHAP 可解释分析从黑箱里挖出特征贡献度3.1 SHAP 值在回归问题里的数学含义SHAP 的核心思想来自合作博弈论里的 Shapley 值把每个特征看作一个“玩家”模型预测值看作“总收益”每个特征分到的收益就是它对预测结果的贡献。对于回归问题SHAP 值满足三个性质可加性所有特征 SHAP 值之和等于预测值减去基线值、对称性两个贡献相同的特征 SHAP 值相同、一致性特征贡献变大时 SHAP 值不会减小。这些性质保证了归因的合理性比简单的特征重要性排序靠谱得多。对于 BP 网络这种非线性模型精确计算 SHAP 值需要遍历所有特征子集计算量是 2 的 n 次方。实际用的是 KernelSHAP 或 DeepSHAP 近似算法。KernelSHAP 把 SHAP 值计算转化为一个加权线性回归问题对每个样本采样若干特征子集用模型预测值拟合。DeepSHAP 则利用神经网络的反向传播把 SHAP 值分解到每一层效率更高但要求网络结构已知。在 MATLAB 里没有官方 SHAP 工具箱常见做法有两种一是调用 Python 的shap库通过 MATLAB 的 Python 接口传数据二是自己实现 KernelSHAP 的核心逻辑。我一般用第一种因为 Python 的shap库成熟稳定MATLAB 只负责训练网络和导出预测函数。3.2 用 MATLAB 训练网络并导出预测接口要让 Python 的 SHAP 库能调用 MATLAB 训练好的网络最干净的方式是把网络导出为可独立调用的函数。MATLAB 提供了genFunction函数可以把训练好的网络转成纯 MATLAB 代码不依赖神经网络工具箱。% 训练完成后导出网络为函数 genFunction(net, bpNetPredict, MatrixOnly, yes); % 测试导出的函数 Y_check bpNetPredict(X); fprintf(导出函数与网络预测最大差异: %.2e\n, max(abs(Y_check(:) - net(X)(:))));genFunction生成的bpNetPredict.m文件包含了网络的所有权重、偏置和归一化参数输入输出都是矩阵格式。MatrixOnly设为yes表示只接受矩阵输入不接受元胞数组这样在 Python 里调用更方便。导出后一定要用max(abs(...))验证一下差异应该在 1e-10 量级如果大了说明导出过程有问题。接下来在 Python 里通过matlab.engine调用这个函数import matlab.engine import numpy as np import shap # 启动 MATLAB 引擎 eng matlab.engine.start_matlab() eng.cd(rC:\your_project_path, nargout0) # 准备数据X_py 是 numpy 数组形状 (n_samples, n_features) X_py np.load(X_for_shap.npy) X_matlab matlab.double(X_py.tolist()) # 调用 MATLAB 预测函数 Y_pred eng.bpNetPredict(X_matlab) Y_pred np.array(Y_pred) # 用 KernelSHAP 解释 # 注意这里需要一个包装函数输入 numpy 返回 numpy def model_predict(X): X_m matlab.double(X.tolist()) Y eng.bpNetPredict(X_m) return np.array(Y).T # 转置为 (n_samples, n_outputs) # 对第一个输出做 SHAP 分析 explainer shap.KernelExplainer( lambda x: model_predict(x)[:, 0], # 只取第一个输出 shap.sample(X_py, 50) # 用 50 个背景样本 ) shap_values explainer.shap_values(X_py[:100], nsamples200)这段代码的关键点matlab.double把 numpy 数组转成 MATLAB 能识别的双精度矩阵model_predict包装函数负责在 Python 和 MATLAB 之间转换数据格式shap.KernelExplainer的第一个参数是预测函数第二个参数是背景数据集背景样本数一般取 50 到 100太少会导致 SHAP 值方差大太多计算慢。nsamples200是每个样本采样的特征子集数越大越精确但计算时间线性增长。3.3 SHAP 图怎么看从 summary plot 到 dependence plotSHAP 分析跑完后核心产出是三类图summary plot、dependence plot 和 force plot。summary plot 把每个特征的 SHAP 值分布画成蜂群图横轴是 SHAP 值纵轴是特征名颜色表示特征值高低。看这张图能快速判断哪些特征重要SHAP 绝对值大、影响方向是什么特征值高时 SHAP 正还是负。import matplotlib.pyplot as plt # Summary plot shap.summary_plot(shap_values, X_py[:100], feature_names[fX{i1} for i in range(6)]) plt.savefig(shap_summary.png, dpi300, bbox_inchestight) # Dependence plot看第 3 个特征对第 1 个输出的影响 shap.dependence_plot(2, shap_values, X_py[:100], feature_names[fX{i1} for i in range(6)]) plt.savefig(shap_dependence_X3.png, dpi300, bbox_inchestight)summary_plot的feature_names参数建议用有物理意义的名称比如[温度, 压力, 流速, ...]别用X1、X2否则图给领导看的时候还得解释。dependence_plot的第一个参数是特征索引从 0 开始。这张图能看出特征与 SHAP 值的关系是线性还是非线性如果散点呈现明显的曲线说明 BP 网络捕捉到了非线性效应这正是用神经网络而不是线性回归的理由。提示SHAP 值有正负正表示该特征把预测值推高负表示推低。summary plot 里如果某个特征的 SHAP 值集中在 0 附近说明这个特征对模型几乎没贡献可以考虑剔除后重新训练简化模型。4. 新数据预测从单条样本到批量推理的完整链路4.1 新数据预处理的三个必须对齐新数据预测翻车十有八九是预处理没对齐。训练时用的归一化参数、缺失值处理方式、异常值截断阈值在新数据上必须一模一样。MATLAB 的genFunction导出的函数已经包含了训练时的归一化参数所以只要新数据的原始量纲和训练数据一致直接调用就行。但如果你在训练前手动做过缺失值填充或异常值替换新数据也得走同样的流程。% 新数据一条样本6 个输入 newSample [350, 0.35, 12.5, 80, 2.1, 0.9]; % 直接调用导出的函数 prediction bpNetPredict(newSample); fprintf(预测输出: %.4f, %.4f, %.4f\n, prediction(1), prediction(2), prediction(3)); % 批量预测100 条新样本 newBatch rand(100, 6) .* [500, 0.5, 20, 100, 3, 1.5]; % 模拟新数据 predBatch bpNetPredict(newBatch); fprintf(批量预测维度: %d×%d\n, size(predBatch,1), size(predBatch,2));注意newSample是 1×6 的行向量转置后变成 6×1 的列向量符合网络输入要求。predBatch是 3×100每列一个样本的三个输出。如果新数据的某个特征超出了训练数据的范围BP 网络会外推但外推可靠性随超出程度增加而下降。我一般会检查新数据每个特征是否在训练数据的 [min, max] 范围内超出 20% 以上的样本标记出来人工复核。4.2 预测结果的置信区间估计BP 网络给出的是点预测没有置信区间。但在工程决策里光有点预测不够还需要知道预测的不确定性。常用做法是集成多个不同初始化的网络用预测值的均值和标准差作为置信区间的近似。% 训练 10 个不同初始化的网络 nEnsemble 10; Y_ensemble zeros(nOutput, size(X_new, 2), nEnsemble); for i 1:nEnsemble rng(i * 100); % 不同随机种子 net_i feedforwardnet(nHidden, trainlm); net_i.trainParam.showWindow false; net_i train(net_i, X, Y); genFunction(net_i, sprintf(bpNetPredict_%d, i), MatrixOnly, yes); Y_ensemble(:, :, i) feval(sprintf(bpNetPredict_%d, i), X_new); end % 计算均值和标准差 Y_mean mean(Y_ensemble, 3); Y_std std(Y_ensemble, 0, 3); % 95% 置信区间近似 Y_lower Y_mean - 1.96 * Y_std; Y_upper Y_mean 1.96 * Y_std; fprintf(第一个输出的 95%% 置信区间宽度均值: %.4f\n, mean(Y_upper(1,:) - Y_lower(1,:)));这段代码训练 10 个网络每个用不同随机种子预测结果取均值和标准差。std的第二个参数 0 表示按 N-1 归一化第三个参数 3 表示沿第三维集成维度计算。置信区间宽度反映了模型在这个样本上的不确定性宽度大的样本建议人工复核。这个方法的计算成本是单网络的 10 倍如果训练一个网络要 5 分钟集成就要 50 分钟适合离线批量预测不适合实时推理。4.3 把预测和 SHAP 解释串成一条流水线实际项目里新数据预测和 SHAP 解释往往需要一起交付。比如给一批新样本既要预测值又要知道每个样本的预测主要受哪个特征驱动。这时候可以把预测和 SHAP 分析串成一个脚本输入原始数据输出预测表加解释图。% 完整流水线新数据预测 SHAP 解释 function [Y_pred, shap_values] predictWithExplanation(X_new, model_path) % 加载导出的预测函数 addpath(model_path); % 预测 Y_pred bpNetPredict(X_new); % 导出新数据供 Python SHAP 使用 writematrix(X_new, X_new_for_shap.csv); writematrix(Y_pred, Y_pred_for_shap.csv); % 调用 Python 脚本做 SHAP 分析 system(python run_shap_analysis.py); % 读取 SHAP 结果 shap_values readmatrix(shap_values.csv); fprintf(预测完成SHAP 分析完成\n); end这个函数把 MATLAB 预测和 Python SHAP 分析串起来中间用 CSV 文件交换数据。writematrix和readmatrix是 MATLAB 里最稳定的 CSV 读写函数。system调用 Python 脚本时确保 Python 环境里装了shap、numpy、matlab.engine等依赖。如果 Python 脚本报错MATLAB 这边不会自动捕获建议在system调用后检查返回状态码。注意system调用 Python 时工作目录要和 Python 脚本里读写文件的路径一致否则会找不到文件。我一般用绝对路径省得排查路径问题。5. 避坑与排查多输入多输出 BP SHAP 的五个高频翻车点5.1 训练集 MSE 很低但测试集 MSE 爆炸现象训练完看tr.best_perf是 1e-6 量级但拿测试集一算 MSE 是 0.5差了五个数量级。原因过拟合。隐含层节点太多、训练轮数太多、样本量太少三者占一个就会这样。解决先减隐含层节点从sqrt(nInputnOutput)5减到sqrt(nInputnOutput)再把max_fail从 6 降到 4让早停更激进如果样本确实少用trainbr贝叶斯正则化替代trainlm它自带正则项抗过拟合能力强。5.2 SHAP 值全为正或全为负现象summary plot 里所有特征的 SHAP 值都在零线同一侧看起来每个特征都在推高或推低预测。原因背景数据集选得不对。KernelSHAP 的基线是背景数据集的平均预测值如果背景数据集和解释数据集分布差异大SHAP 值会整体偏移。解决背景数据集从训练集里随机采样别从测试集或新数据里采。样本数 50 到 100 之间太少方差大太多计算慢。另外检查model_predict函数返回的维度是否和shap_values期望的一致多输出时只取一个输出做解释。5.3 新数据预测结果全是 NaN现象bpNetPredict(newSample)返回 NaN。原因新数据里有 NaN 或 Inf。BP 网络的矩阵运算遇到 NaN 会传播到整个输出。解决预测前检查any(isnan(newSample))和any(isinf(newSample))有的话先填充或剔除。另外检查新数据的量纲是否和训练数据一致如果训练时输入是 0 到 1新数据是 0 到 100归一化后可能超出 [-1, 1] 范围但不会产生 NaN只会预测不准。5.4 MATLAB 和 Python 数据交换时维度对不上现象Python 里matlab.double(X.tolist())传过去后 MATLAB 报维度错误。原因tolist()把 numpy 数组转成嵌套列表matlab.double默认按行优先解释而 MATLAB 是列优先。如果 numpy 数组是 (n_samples, n_features)转过去 MATLAB 看到的是 (n_features, n_samples)正好转置了。解决在 Python 里转置一下matlab.double(X.T.tolist())或者在 MATLAB 里再转置一次。我一般约定 Python 端传转置后的数据MATLAB 端不再转减少混乱。5.5 genFunction 导出的函数预测结果和原网络不一致现象bpNetPredict(X)和net(X)的结果差很多。原因genFunction默认不包含归一化参数或者导出时MatrixOnly设成了no。解决导出时明确指定MatrixOnly, yes并且检查生成的.m文件里是否有mapminmax_apply和mapminmax_reverse的调用。如果没有说明归一化没导出需要手动在genFunction之前设置net.inputs{1}.processFcns和net.outputs{2}.processFcns确保包含mapminmax。导出后必须用max(abs(...))验证差异大于 1e-8 就要查。6. 进阶技巧用 SHAP 交互值定位特征协同效应单特征 SHAP 值只能告诉你每个特征独立贡献了多少但多输入场景下特征之间的交互效应往往才是关键。比如温度高时压力对输出的影响可能比温度低时大得多这种协同效应单特征 SHAP 图看不出来。SHAP 交互值SHAP interaction values能拆解出每对特征的联合贡献计算量是单特征 SHAP 的 n 倍n 是特征数6 个特征就是 6 倍还能接受。在 Python 的shap库里有两种方式算交互值。一是shap.TreeExplainer自带shap_interaction_values方法但只支持树模型。BP 网络得用KernelExplainer加shap_interaction_values函数计算更慢但通用。# 计算 SHAP 交互值只对前 20 个样本计算量大 shap_interaction explainer.shap_interaction_values(X_py[:20]) # shap_interaction 形状: (n_samples, n_features, n_features) # 对角线是单特征 SHAP 值非对角线是交互值 # 提取第 0 个样本的第 2 和第 4 个特征的交互值 interaction_2_4 shap_interaction[0, 2, 4] print(f特征3和特征5的交互 SHAP 值: {interaction_2_4:.4f}) # 画交互热力图 import seaborn as sns mean_interaction np.mean(np.abs(shap_interaction), axis0) sns.heatmap(mean_interaction, annotTrue, fmt.3f, xticklabels[fX{i1} for i in range(6)], yticklabels[fX{i1} for i in range(6)]) plt.title(SHAP Interaction Heatmap) plt.savefig(shap_interaction.png, dpi300, bbox_inchestight)shap_interaction_values返回一个三维数组第一维是样本第二维和第三维是特征对。对角线元素就是单特征 SHAP 值非对角线元素是交互值。热力图里颜色越深表示交互效应越强。如果发现某对特征的交互值很大说明这两个特征对输出的影响不是简单叠加而是有协同或拮抗。这时候可以在工艺上重点关注这两个参数的匹配关系而不是单独调一个。我一般会把这个热力图和工艺知识对照如果两个特征在物理上确实有关联比如温度和压力在热力学上耦合那 SHAP 交互值大是合理的说明模型学到了真实规律如果两个特征物理上无关但交互值大可能是数据里的伪相关需要检查采样过程是否有偏差。最后一个习惯每次跑完 SHAP 分析我都会把shap_values和原始数据一起存成.mat文件命名带上日期和模型版本。因为 SHAP 计算耗时下次想复现某张图不用重跑。这个习惯帮我省了至少几十个小时的重复计算。希望帮到你。本文还有配套的精品资源点击获取