简介MATLAB实现GMDH自组织网络模型多输入回归预测的完整源码与数据包面向需要完成多输入回归预测的开发者数据包含7个特征输入与1个输出变量可直接用于教学实验、算法对比或工程预研。压缩包共11个文件以6个m源码文件为核心覆盖网络构建、神经元生成、结果绘图等关键环节docx说明文档介绍了使用方式xlsx提供可直接运行的多输入样例数据另含2张效果示意图整体仅232KB便于快速下载与部署。运行环境需为MATLAB2018b及以上源码若出现乱码多为版本编码差异用记事本打开复制到自己的文件中即可正常使用。目前已有583人学习下载适合具备一定MATLAB基础、希望快速掌握GMDH建模流程的读者。通过替换自身数据即可快速验证模型效果同时可深入理解GMDH自组织网络的结构演化与神经元生成逻辑为后续算法改进或工程应用提供可运行的基础代码。1. 多输入回归预测为什么绕不开 GMDH 自组织网络样本量只有两三百个、输入却有七八个维度的回归任务第一反应往往是直接上 MLP 或随机森林。训练到一半才发现验证集 RMSE 比线性回归还难看。GMDH 自组织网络模型在这种小样本、多输入、非线性交互的场景里反而更稳它不拍脑袋去拟合一整张高阶多项式表而是让每一层自动把两个变量组合成一个二次多项式神经元再按验证误差挑出有用的神经元层层生成直到误差不再下降。这个思路最早来自 Ivakhnenko今天用 MATLAB 复现并不复杂核心训练加预测函数可以控制在百行以内特别适合做多输入回归预测、特征筛选和快速基线模型。下面从建模原理开始再落到参数、代码与常见坑。2. GMDH 自组织网络的建模原理与数学表达2.1 从 Kolmogorov-Gabor 多项式到二阶神经元GMDH 的数学起点是 Kolmogorov-Gabor 多项式y a0 Σ ai xi Σ aij xi xj Σ aijk xi xj xk ...这个式子理论上可以描述任意连续映射但 m 个输入一旦超过 10交叉项的数目会迅速膨胀直接求解既不现实也容易过拟合。GMDH 的做法是把目标拆成大量小的二阶多项式每个神经元只接收两个输入输出一个中间变量z w0 w1 ui w2 uj w3 ui^2 w4 uj^2 w5 ui uj其中 ui、uj 可以是原始输入也可以是上一层神经元的输出。每个这样的二阶多项式就是一个神经元它的六个系数用最小二乘在训练集上估计。如果输入之间存在乘法或平方关系这类神经元能比纯线性模型更早把它们抓出来同时因为每层只做二阶组合不会出现三阶、四阶多项式在边界处剧烈振荡的问题。在 MATLAB 里这个设计矩阵可以写成% poly2_design.m function D poly2_design(u, v) % u, v 是两个输入变量的列向量长度等于样本数 % 列顺序常数项, u, v, u^2, v^2, u*v D [ones(size(u)), u, v, u.^2, v.^2, u.*v]; end拿到 D 之后系数就是b D \ yMATLAB 会走 QR/最小二乘路径。如果特征列之间严重共线建议把D \ y改成pinv(D) * y避免数值警告导致系数不稳定。2.2 自组织选择训练集拟合、验证集筛选GMDH 里的自组织主要体现在不对最终网络结构做先验假设。每一层都会产生一批候选神经元再由外部准则决定哪些进入下一层常见流程是把样本分成拟合集和验证集当前层有 p 个输入变量枚举 C(p, 2) 个变量对每个变量对生成一个二阶多项式神经元用拟合集估计六个系数用验证集计算该神经元的输出 RMSE把 C(p, 2) 个候选按 RMSE 排序保留前 K 个作为下一层输入。这个“拟合在训练集上、筛选在验证集上”的分离是 GMDH 和普通多项式回归最大的区别。普通多项式回归只做一个全局最小二乘而 GMDH 会把每一层的网络规模、深度、神经元连接方式全部交给数据去选择。外部准则的选择会直接影响模型稳定性准则计算方式适用场景RMSEsqrt(mean((y - yhat).^2))最常用直接对应回归误差修正 AIC/BIC在 RMSE 基础上加入参数惩罚样本少、特征多时防止过拟合最小偏差准则比较不同样本划分下的模型输出差数据噪声大时需要稳定性实际项目中我一般先用验证集 RMSE 做筛选等模型层数确定后再用 AIC/BIC 比较最终候选这样计算量小结果也容易解释。2.3 层数停止与组合爆炸控制GMDH 不会无限制叠层。每层输入变量数是上一层的 K所以网络宽度通常被限制在一个可控范围内但第一层的组合数仍可能很大。假设输入 p10C(10,2)45 对p50C(50,2)1225 对p200C(200,2)19900 对。如果每个候选都做一次六列最小二乘计算量会快速上涨。一种常见做法是在每层随机抽一部分变量对参与竞争避免全枚举。代码里可以这样限流p size(Xcur, 2); pairs nchoosek(1:p, 2); max_pairs 500; if size(pairs, 1) max_pairs pairs pairs(randperm(size(pairs, 1), max_pairs), :); end这里不是简单降精度GMDH 的每一层都有冗余候选随机抽样相当于对所有可能交互项做了一次子采样配合后面保留前 K 个节点的筛选仍然能稳定找到有效特征组合。只有 p 很小的时候才建议全枚举。层数停止用提前终止实现如果当前层最小验证 RMSE 比上一层没有明显下降就回退到上一层模型。这个逻辑和神经网络里的 early stopping 本质上是一回事只不过 GMDH 的每一层就是一次显式的结构增长。3. 用 MATLAB 搭建 GMDH 自组织网络核心函数3.1 数据标准化与训练/验证集划分多输入回归里特征量纲差异可能很大。GMDH 的二次多项式会同时出现 x、x^2、x1*x2 这样的项量纲不统一会让设计矩阵的条件数变差系数估计也容易失真。所以第一步做 z-score 标准化。标准化的均值和标准差只能在训练集上计算验证集和测试集必须复用同一组参数。如果验证集单独算均值会造成信息泄漏验证误差会偏乐观。mu_x mean(Xtr, 1); sd_x std(Xtr, 0, 1); Xtr_n (Xtr - mu_x) ./ sd_x; Xva_n (Xva - mu_x) ./ sd_x; mu_y mean(ytr); sd_y std(ytr); ytr_n (ytr - mu_y) ./ sd_y; yva_n (yva - mu_y) ./ sd_y;std(Xtr, 0, 1)里的第二个参数 0 表示除以 n-1和默认行为一致。输出 y 也做标准化这样验证 RMSE 是在无量纲尺度上比较不会受某个大数值输出列主导。3.2 GMDH 主循环组合、二次回归、保留 K 个节点下面是 GMDH 训练函数的核心骨架。它把每层生成的所有候选节点记为列计算验证 RMSE保留最好的 K 个并记录全局最优层和最优节点。为了控制文章篇幅这里省略了矩阵预分配之外的防御性检查但整体逻辑可直接运行。function [model, hist] gmdh_train(Xtr, ytr, Xva, yva, ... max_layers, keep, tol, max_pairs) % GMDH 自组织网络训练 % Xtr, ytr: 标准化后的训练输入输出 % Xva, yva: 标准化后的验证输入输出 % max_layers: 最大层数 % keep: 每层保留的节点数 % tol: 验证 RMSE 提升阈值小于该值则停止 % max_pairs: 每层最多尝试的变量对数 model.layers {}; model.best_layer 0; model.best_node 0; best_rmse inf; hist []; Xcur_tr Xtr; % 当前层训练输入 Xcur_va Xva; % 当前层验证输入 for layer 1:max_layers p size(Xcur_tr, 2); if p 2 break; end % 生成变量对超过上限则随机抽样 pairs nchoosek(1:p, 2); if size(pairs, 1) max_pairs pairs pairs(randperm(size(pairs, 1), max_pairs), :); end npairs size(pairs, 1); coefs cell(npairs, 1); ztr_all zeros(size(Xcur_tr, 1), npairs); zva_all zeros(size(Xcur_va, 1), npairs); rmse zeros(npairs, 1); for k 1:npairs i pairs(k, 1); j pairs(k, 2); Dtr poly2_design(Xcur_tr(:, i), Xcur_tr(:, j)); Dva poly2_design(Xcur_va(:, i), Xcur_va(:, j)); b Dtr \ ytr; % 训练集拟合 ztr Dtr * b; % 训练输出 zva Dva * b; % 验证输出 coefs{k} b; ztr_all(:, k) ztr; zva_all(:, k) zva; rmse(k) sqrt(mean((zva - yva).^2)); % 外部准则 end [sorted_rmse, idx] sort(rmse); layer_best_rmse sorted_rmse(1); layer_best_node idx(1); if layer_best_rmse best_rmse best_rmse layer_best_rmse; model.best_layer layer; model.best_node layer_best_node; end keep_idx idx(1:min(keep, npairs)); model.layers{layer} struct(... pairs, pairs, coefs, {coefs}, ... keep, keep_idx, rmse, rmse); % 提前停止提升小于 tol 就退出 if layer 1 (hist(end, 3) - layer_best_rmse) tol break; end Xcur_tr ztr_all(:, keep_idx); Xcur_va zva_all(:, keep_idx); hist [hist; layer, sorted_rmse(1), mean(sorted_rmse(1:min(keep, npairs)))]; end end这里的核心是keep_idx和Xcur_tr的更新。keep_idx保存当前层应该继续向下传的节点位置Xcur_tr则变成这些节点的输出。下一层枚举变量对时面对的是上一层的中间变量而不是原始特征这就是 GMDH 所谓“自组织网络结构”的体现。model.best_layer和model.best_node记录全局验证误差最小的那个节点。即使继续叠层导致验证误差反弹最终预测仍会回到这个最佳位置相当于隐式剪枝。参数建议范围参数作用推荐起点max_layers最大网络深度8 到 15keep每层保留节点数3 到 8tol验证 RMSE 提升阈值1e-4 到 1e-3max_pairs每层最大变量对数200 到 6003.3 前向预测与反标准化训练完成后预测要按同样的层结构逐层前向传播。如果还没达到best_layer需要把该层保留的 K 个节点全部算出来作为下一层输入到达best_layer后只输出最佳节点即可。function yhat gmdh_predict(model, X) % GMDH 前向预测X 必须是按训练集统计量标准化后的输入 Xcur X; for layer 1:model.best_layer lay model.layers{layer}; pairs lay.pairs; coefs lay.coefs; keep lay.keep; if layer model.best_layer out zeros(size(Xcur, 1), numel(keep)); for r 1:numel(keep) k keep(r); D poly2_design(Xcur(:, pairs(k, 1)), ... Xcur(:, pairs(k, 2))); out(:, r) D * coefs{k}; end Xcur out; else k model.best_node; D poly2_design(Xcur(:, pairs(k, 1)), ... Xcur(:, pairs(k, 2))); yhat D * coefs{k}; end end endgmdh_train和gmdh_predict都依赖poly2_design.m这个函数在第二章已经给出。预测完成后再把标准化后的输出还原回原始量纲yhat_original yhat * sd_y mu_y;这里的sd_y、mu_y必须是训练集上计算的那一组不能用验证集重新计算否则还原后的数值会带着偏差。4. 多输入回归预测样例数据准备、训练与误差对比4.1 用模拟数据跑通最小流程先构造一组包含线性和交叉项的数据用来确认 GMDH 能捕捉非线性关系。实际项目里把下面这段里的模拟生成换成readmatrix读 CSV 即可rng(7); n 400; p 6; X randn(n, p); % 真实关系包含线性项、交叉项和平方项 y 1.2 * X(:, 1) - 0.9 * X(:, 2) 0.7 * X(:, 3) .* X(:, 4) ... 0.2 * X(:, 5).^2 0.05 * randn(n, 1); cv cvpartition(n, HoldOut, 0.25); Xtr X(cv.training, :); ytr y(cv.training, :); Xva X(cv.test, :); yva y(cv.test, :);然后做标准化并训练mu_x mean(Xtr, 1); sd_x std(Xtr, 0, 1); Xtr_n (Xtr - mu_x) ./ sd_x; Xva_n (Xva - mu_x) ./ sd_x; mu_y mean(ytr); sd_y std(ytr); ytr_n (ytr - mu_y) ./ sd_y; yva_n (yva - mu_y) ./ sd_y; [model, hist] gmdh_train(Xtr_n, ytr_n, Xva_n, yva_n, ... 8, 5, 1e-4, 300); yhat gmdh_predict(model, Xva_n); yhat yhat * sd_y mu_y; rmse_gmdh sqrt(mean((yva - yhat).^2)); fprintf(GMDH 验证集 RMSE: %.4f\n, rmse_gmdh);这里yva用的是原始量纲因为预测已经反标准化。hist里记录了每层最小验证 RMSE可以画出来观察网络是否还有必要加深。如果拿到的是发布包里的 CSV常见格式是最后一列是输出前面所有列是特征读取方式为data readmatrix(多输入回归预测数据.csv); X data(:, 1:end-1); y data(:, end);readmatrix 在 R2020a 之后都可用处理带表头的 CSV 时可以加NumHeaderLines参数具体行数按文件实际情况调整。4.2 和线性回归、BP 神经网络对比线性回归作为基线能直接告诉我们数据里的非线性项对误差贡献有多大lm fitlm(Xtr_n, ytr_n); y_lm predict(lm, Xva_n) * sd_y mu_y; rmse_lm sqrt(mean((yva - y_lm).^2));如果数据里只有纯线性关系GMDH 的 RMSE 会和线性回归接近如果存在交叉项或平方项GMDH 通常会有明确优势。前面构造的模拟数据里GMDH 因为能选出 x3*x4 和 x5^2 两个交互结构RMSE 一般比线性回归低 30% 左右。想和神经网络比可以用 MATLAB 里的feedforwardnetnet feedforwardnet(10); net.trainFcn trainlm; net train(net, Xtr_n, ytr_n); y_nn net(Xva_n); y_nn y_nn * sd_y mu_y; rmse_nn sqrt(mean((yva - y_nn).^2));注意神经网络数据要转置成“特征数 × 样本数”的列排列。和 BP 相比GMDH 的明显优势是训练结果稳定不需要设置学习率、动量项或隐藏层节点数主要参数只有层数和每层保留节点数。误差对比可以直接放在一张表里模型验证 RMSE说明线性回归由rmse_lm输出只含线性结构GMDH由rmse_gmdh输出自动构造二阶多项式节点BP 神经网络由rmse_nn输出受随机初始化影响较大5. GMDH 自组织网络的调参、过拟合与模型边界5.1 四个直接影响结果的参数GMDH 需要重点调的不是学习率而是网络规模相关的参数。keep是最关键的一个它控制每层有多少候选神经元能进入下一层。keep设太大会让下一层变量数变大组合数指数增长设太小又会过早丢掉可能有用的中间变量。常见做法是先固定keep5跑通后再尝试keep3和keep8观察验证 RMSE 是否明显变化。max_layers不要一开始就设成很大。多数表格型回归问题里GMDH 的有效层数在 3 到 8 层之间超过 10 层后基本都在过拟合。tol控制早停灵敏度设成 1e-4 通常够用如果你发现训练曲线抖动可以把tol调到 1e-3让模型更早停下来。如果想自动搜索参数可以用 MATLAB 优化工具箱里的surrogateopt把keep和max_layers当作整数变量处理。不过 GMDH 训练本身已经很快很多时候直接做两层网格搜索反而更直观。5.2 用训练历史判断过拟合GMDH 不像神经网络那样能直接画 loss 曲线但hist里已经保存了每层最小验证 RMSE。最直接的诊断方法是把hist第二列画出来figure; plot(hist(:, 2), o-); xlabel(层数); ylabel(最小验证 RMSE);如果曲线先下降、后上升说明网络在某一层之后开始记忆验证集噪声。gmdh_train里已经通过best_layer自动回退到最优层所以最终预测用的不是最后一层而是历史上验证误差最小的那一层。这个设计能抵消一部分过拟合但不能完全替代验证集划分的合理性。更稳妥的做法是外层再做一次 k 折交叉验证用每次折的验证误差重新选keep和max_layers。GMDH 训练开销小交叉验证代价通常可以接受。5.3 哪些场景不适合 GMDHGMDH 并不适合所有多输入回归问题。输入特征来自图像像素、文本向量或高频序列时原始特征之间缺乏显式可组合的交互结构GMDH 的二次多项式候选会变得很笨重这时用卷积网络或循环网络更合理。另一个边界是特征数量特别大的场景。当 p 超过 100 时即便加max_pairs限流随机抽样也会带来不确定性这时需要先做一遍特征筛选或者先用 LASSO 压缩到二三十个变量再交给 GMDH。还有GMDH 本身不擅长处理缺失值和类别型变量类别特征必须先做 One-hot 或嵌入编码。6. 让 GMDH 在 MATLAB 里更稳特征解释、保存与再训练6.1 用扰动法看特征重要性GMDH 中间层保存的是变量组合关系直接观察第一层的最佳节点只能看到两个中间变量未必能映射回原始特征。更实用的做法是扰动法把验证集某一列打乱观察预测 RMSE 变化多少。变化越大说明该特征对模型输出的影响越大。y0 gmdh_predict(model, Xva_n); imp zeros(1, p); for f 1:p Xp Xva_n; Xp(:, f) Xp(randperm(n), f); yp gmdh_predict(model, Xp); imp(f) sqrt(mean((yp - y0).^2)); end bar(imp); xlabel(特征编号); ylabel(扰动后 RMSE 变化);由于每次扰动都走完整的前向网络这个结果比查看单一系数更可靠。注意特征被完全打乱后如果模型本来就很少使用它RMSE 变化会接近 0如果变化很大说明该特征参与了多个层级的组合是模型的主要依赖变量。6.2 保存模型时必须连同标准化参数一起存实际项目里模型训练完通常要落盘到.mat文件。GMDH 的预测流程依赖训练时生成的层次结构也依赖mu_x、sd_x、mu_y、sd_y这四组统计量。只保存model而丢掉标准化参数是新手最容易踩的坑。save(gmdh_model.mat, model, mu_x, sd_x, mu_y, sd_y);部署预测时S load(gmdh_model.mat); x_new (new_row - S.mu_x) ./ S.sd_x; y_new gmdh_predict(S.model, x_new) * S.sd_y S.mu_y;这样保存的模型和参数绑定在一起换一台机器或换一个 MATLAB 版本加载后都能直接预测。GMDH 的模型体积很小中间层存储的无非是变量对索引和六个系数比同等精度神经网络更适合放进嵌入式环境或实时评分服务。本文还有配套的精品资源点击获取 GMDH自组织网络：MATLAB实现多输入回归预测的完整指南
简介MATLAB实现GMDH自组织网络模型多输入回归预测的完整源码与数据包面向需要完成多输入回归预测的开发者数据包含7个特征输入与1个输出变量可直接用于教学实验、算法对比或工程预研。压缩包共11个文件以6个m源码文件为核心覆盖网络构建、神经元生成、结果绘图等关键环节docx说明文档介绍了使用方式xlsx提供可直接运行的多输入样例数据另含2张效果示意图整体仅232KB便于快速下载与部署。运行环境需为MATLAB2018b及以上源码若出现乱码多为版本编码差异用记事本打开复制到自己的文件中即可正常使用。目前已有583人学习下载适合具备一定MATLAB基础、希望快速掌握GMDH建模流程的读者。通过替换自身数据即可快速验证模型效果同时可深入理解GMDH自组织网络的结构演化与神经元生成逻辑为后续算法改进或工程应用提供可运行的基础代码。1. 多输入回归预测为什么绕不开 GMDH 自组织网络样本量只有两三百个、输入却有七八个维度的回归任务第一反应往往是直接上 MLP 或随机森林。训练到一半才发现验证集 RMSE 比线性回归还难看。GMDH 自组织网络模型在这种小样本、多输入、非线性交互的场景里反而更稳它不拍脑袋去拟合一整张高阶多项式表而是让每一层自动把两个变量组合成一个二次多项式神经元再按验证误差挑出有用的神经元层层生成直到误差不再下降。这个思路最早来自 Ivakhnenko今天用 MATLAB 复现并不复杂核心训练加预测函数可以控制在百行以内特别适合做多输入回归预测、特征筛选和快速基线模型。下面从建模原理开始再落到参数、代码与常见坑。2. GMDH 自组织网络的建模原理与数学表达2.1 从 Kolmogorov-Gabor 多项式到二阶神经元GMDH 的数学起点是 Kolmogorov-Gabor 多项式y a0 Σ ai xi Σ aij xi xj Σ aijk xi xj xk ...这个式子理论上可以描述任意连续映射但 m 个输入一旦超过 10交叉项的数目会迅速膨胀直接求解既不现实也容易过拟合。GMDH 的做法是把目标拆成大量小的二阶多项式每个神经元只接收两个输入输出一个中间变量z w0 w1 ui w2 uj w3 ui^2 w4 uj^2 w5 ui uj其中 ui、uj 可以是原始输入也可以是上一层神经元的输出。每个这样的二阶多项式就是一个神经元它的六个系数用最小二乘在训练集上估计。如果输入之间存在乘法或平方关系这类神经元能比纯线性模型更早把它们抓出来同时因为每层只做二阶组合不会出现三阶、四阶多项式在边界处剧烈振荡的问题。在 MATLAB 里这个设计矩阵可以写成% poly2_design.m function D poly2_design(u, v) % u, v 是两个输入变量的列向量长度等于样本数 % 列顺序常数项, u, v, u^2, v^2, u*v D [ones(size(u)), u, v, u.^2, v.^2, u.*v]; end拿到 D 之后系数就是b D \ yMATLAB 会走 QR/最小二乘路径。如果特征列之间严重共线建议把D \ y改成pinv(D) * y避免数值警告导致系数不稳定。2.2 自组织选择训练集拟合、验证集筛选GMDH 里的自组织主要体现在不对最终网络结构做先验假设。每一层都会产生一批候选神经元再由外部准则决定哪些进入下一层常见流程是把样本分成拟合集和验证集当前层有 p 个输入变量枚举 C(p, 2) 个变量对每个变量对生成一个二阶多项式神经元用拟合集估计六个系数用验证集计算该神经元的输出 RMSE把 C(p, 2) 个候选按 RMSE 排序保留前 K 个作为下一层输入。这个“拟合在训练集上、筛选在验证集上”的分离是 GMDH 和普通多项式回归最大的区别。普通多项式回归只做一个全局最小二乘而 GMDH 会把每一层的网络规模、深度、神经元连接方式全部交给数据去选择。外部准则的选择会直接影响模型稳定性准则计算方式适用场景RMSEsqrt(mean((y - yhat).^2))最常用直接对应回归误差修正 AIC/BIC在 RMSE 基础上加入参数惩罚样本少、特征多时防止过拟合最小偏差准则比较不同样本划分下的模型输出差数据噪声大时需要稳定性实际项目中我一般先用验证集 RMSE 做筛选等模型层数确定后再用 AIC/BIC 比较最终候选这样计算量小结果也容易解释。2.3 层数停止与组合爆炸控制GMDH 不会无限制叠层。每层输入变量数是上一层的 K所以网络宽度通常被限制在一个可控范围内但第一层的组合数仍可能很大。假设输入 p10C(10,2)45 对p50C(50,2)1225 对p200C(200,2)19900 对。如果每个候选都做一次六列最小二乘计算量会快速上涨。一种常见做法是在每层随机抽一部分变量对参与竞争避免全枚举。代码里可以这样限流p size(Xcur, 2); pairs nchoosek(1:p, 2); max_pairs 500; if size(pairs, 1) max_pairs pairs pairs(randperm(size(pairs, 1), max_pairs), :); end这里不是简单降精度GMDH 的每一层都有冗余候选随机抽样相当于对所有可能交互项做了一次子采样配合后面保留前 K 个节点的筛选仍然能稳定找到有效特征组合。只有 p 很小的时候才建议全枚举。层数停止用提前终止实现如果当前层最小验证 RMSE 比上一层没有明显下降就回退到上一层模型。这个逻辑和神经网络里的 early stopping 本质上是一回事只不过 GMDH 的每一层就是一次显式的结构增长。3. 用 MATLAB 搭建 GMDH 自组织网络核心函数3.1 数据标准化与训练/验证集划分多输入回归里特征量纲差异可能很大。GMDH 的二次多项式会同时出现 x、x^2、x1*x2 这样的项量纲不统一会让设计矩阵的条件数变差系数估计也容易失真。所以第一步做 z-score 标准化。标准化的均值和标准差只能在训练集上计算验证集和测试集必须复用同一组参数。如果验证集单独算均值会造成信息泄漏验证误差会偏乐观。mu_x mean(Xtr, 1); sd_x std(Xtr, 0, 1); Xtr_n (Xtr - mu_x) ./ sd_x; Xva_n (Xva - mu_x) ./ sd_x; mu_y mean(ytr); sd_y std(ytr); ytr_n (ytr - mu_y) ./ sd_y; yva_n (yva - mu_y) ./ sd_y;std(Xtr, 0, 1)里的第二个参数 0 表示除以 n-1和默认行为一致。输出 y 也做标准化这样验证 RMSE 是在无量纲尺度上比较不会受某个大数值输出列主导。3.2 GMDH 主循环组合、二次回归、保留 K 个节点下面是 GMDH 训练函数的核心骨架。它把每层生成的所有候选节点记为列计算验证 RMSE保留最好的 K 个并记录全局最优层和最优节点。为了控制文章篇幅这里省略了矩阵预分配之外的防御性检查但整体逻辑可直接运行。function [model, hist] gmdh_train(Xtr, ytr, Xva, yva, ... max_layers, keep, tol, max_pairs) % GMDH 自组织网络训练 % Xtr, ytr: 标准化后的训练输入输出 % Xva, yva: 标准化后的验证输入输出 % max_layers: 最大层数 % keep: 每层保留的节点数 % tol: 验证 RMSE 提升阈值小于该值则停止 % max_pairs: 每层最多尝试的变量对数 model.layers {}; model.best_layer 0; model.best_node 0; best_rmse inf; hist []; Xcur_tr Xtr; % 当前层训练输入 Xcur_va Xva; % 当前层验证输入 for layer 1:max_layers p size(Xcur_tr, 2); if p 2 break; end % 生成变量对超过上限则随机抽样 pairs nchoosek(1:p, 2); if size(pairs, 1) max_pairs pairs pairs(randperm(size(pairs, 1), max_pairs), :); end npairs size(pairs, 1); coefs cell(npairs, 1); ztr_all zeros(size(Xcur_tr, 1), npairs); zva_all zeros(size(Xcur_va, 1), npairs); rmse zeros(npairs, 1); for k 1:npairs i pairs(k, 1); j pairs(k, 2); Dtr poly2_design(Xcur_tr(:, i), Xcur_tr(:, j)); Dva poly2_design(Xcur_va(:, i), Xcur_va(:, j)); b Dtr \ ytr; % 训练集拟合 ztr Dtr * b; % 训练输出 zva Dva * b; % 验证输出 coefs{k} b; ztr_all(:, k) ztr; zva_all(:, k) zva; rmse(k) sqrt(mean((zva - yva).^2)); % 外部准则 end [sorted_rmse, idx] sort(rmse); layer_best_rmse sorted_rmse(1); layer_best_node idx(1); if layer_best_rmse best_rmse best_rmse layer_best_rmse; model.best_layer layer; model.best_node layer_best_node; end keep_idx idx(1:min(keep, npairs)); model.layers{layer} struct(... pairs, pairs, coefs, {coefs}, ... keep, keep_idx, rmse, rmse); % 提前停止提升小于 tol 就退出 if layer 1 (hist(end, 3) - layer_best_rmse) tol break; end Xcur_tr ztr_all(:, keep_idx); Xcur_va zva_all(:, keep_idx); hist [hist; layer, sorted_rmse(1), mean(sorted_rmse(1:min(keep, npairs)))]; end end这里的核心是keep_idx和Xcur_tr的更新。keep_idx保存当前层应该继续向下传的节点位置Xcur_tr则变成这些节点的输出。下一层枚举变量对时面对的是上一层的中间变量而不是原始特征这就是 GMDH 所谓“自组织网络结构”的体现。model.best_layer和model.best_node记录全局验证误差最小的那个节点。即使继续叠层导致验证误差反弹最终预测仍会回到这个最佳位置相当于隐式剪枝。参数建议范围参数作用推荐起点max_layers最大网络深度8 到 15keep每层保留节点数3 到 8tol验证 RMSE 提升阈值1e-4 到 1e-3max_pairs每层最大变量对数200 到 6003.3 前向预测与反标准化训练完成后预测要按同样的层结构逐层前向传播。如果还没达到best_layer需要把该层保留的 K 个节点全部算出来作为下一层输入到达best_layer后只输出最佳节点即可。function yhat gmdh_predict(model, X) % GMDH 前向预测X 必须是按训练集统计量标准化后的输入 Xcur X; for layer 1:model.best_layer lay model.layers{layer}; pairs lay.pairs; coefs lay.coefs; keep lay.keep; if layer model.best_layer out zeros(size(Xcur, 1), numel(keep)); for r 1:numel(keep) k keep(r); D poly2_design(Xcur(:, pairs(k, 1)), ... Xcur(:, pairs(k, 2))); out(:, r) D * coefs{k}; end Xcur out; else k model.best_node; D poly2_design(Xcur(:, pairs(k, 1)), ... Xcur(:, pairs(k, 2))); yhat D * coefs{k}; end end endgmdh_train和gmdh_predict都依赖poly2_design.m这个函数在第二章已经给出。预测完成后再把标准化后的输出还原回原始量纲yhat_original yhat * sd_y mu_y;这里的sd_y、mu_y必须是训练集上计算的那一组不能用验证集重新计算否则还原后的数值会带着偏差。4. 多输入回归预测样例数据准备、训练与误差对比4.1 用模拟数据跑通最小流程先构造一组包含线性和交叉项的数据用来确认 GMDH 能捕捉非线性关系。实际项目里把下面这段里的模拟生成换成readmatrix读 CSV 即可rng(7); n 400; p 6; X randn(n, p); % 真实关系包含线性项、交叉项和平方项 y 1.2 * X(:, 1) - 0.9 * X(:, 2) 0.7 * X(:, 3) .* X(:, 4) ... 0.2 * X(:, 5).^2 0.05 * randn(n, 1); cv cvpartition(n, HoldOut, 0.25); Xtr X(cv.training, :); ytr y(cv.training, :); Xva X(cv.test, :); yva y(cv.test, :);然后做标准化并训练mu_x mean(Xtr, 1); sd_x std(Xtr, 0, 1); Xtr_n (Xtr - mu_x) ./ sd_x; Xva_n (Xva - mu_x) ./ sd_x; mu_y mean(ytr); sd_y std(ytr); ytr_n (ytr - mu_y) ./ sd_y; yva_n (yva - mu_y) ./ sd_y; [model, hist] gmdh_train(Xtr_n, ytr_n, Xva_n, yva_n, ... 8, 5, 1e-4, 300); yhat gmdh_predict(model, Xva_n); yhat yhat * sd_y mu_y; rmse_gmdh sqrt(mean((yva - yhat).^2)); fprintf(GMDH 验证集 RMSE: %.4f\n, rmse_gmdh);这里yva用的是原始量纲因为预测已经反标准化。hist里记录了每层最小验证 RMSE可以画出来观察网络是否还有必要加深。如果拿到的是发布包里的 CSV常见格式是最后一列是输出前面所有列是特征读取方式为data readmatrix(多输入回归预测数据.csv); X data(:, 1:end-1); y data(:, end);readmatrix 在 R2020a 之后都可用处理带表头的 CSV 时可以加NumHeaderLines参数具体行数按文件实际情况调整。4.2 和线性回归、BP 神经网络对比线性回归作为基线能直接告诉我们数据里的非线性项对误差贡献有多大lm fitlm(Xtr_n, ytr_n); y_lm predict(lm, Xva_n) * sd_y mu_y; rmse_lm sqrt(mean((yva - y_lm).^2));如果数据里只有纯线性关系GMDH 的 RMSE 会和线性回归接近如果存在交叉项或平方项GMDH 通常会有明确优势。前面构造的模拟数据里GMDH 因为能选出 x3*x4 和 x5^2 两个交互结构RMSE 一般比线性回归低 30% 左右。想和神经网络比可以用 MATLAB 里的feedforwardnetnet feedforwardnet(10); net.trainFcn trainlm; net train(net, Xtr_n, ytr_n); y_nn net(Xva_n); y_nn y_nn * sd_y mu_y; rmse_nn sqrt(mean((yva - y_nn).^2));注意神经网络数据要转置成“特征数 × 样本数”的列排列。和 BP 相比GMDH 的明显优势是训练结果稳定不需要设置学习率、动量项或隐藏层节点数主要参数只有层数和每层保留节点数。误差对比可以直接放在一张表里模型验证 RMSE说明线性回归由rmse_lm输出只含线性结构GMDH由rmse_gmdh输出自动构造二阶多项式节点BP 神经网络由rmse_nn输出受随机初始化影响较大5. GMDH 自组织网络的调参、过拟合与模型边界5.1 四个直接影响结果的参数GMDH 需要重点调的不是学习率而是网络规模相关的参数。keep是最关键的一个它控制每层有多少候选神经元能进入下一层。keep设太大会让下一层变量数变大组合数指数增长设太小又会过早丢掉可能有用的中间变量。常见做法是先固定keep5跑通后再尝试keep3和keep8观察验证 RMSE 是否明显变化。max_layers不要一开始就设成很大。多数表格型回归问题里GMDH 的有效层数在 3 到 8 层之间超过 10 层后基本都在过拟合。tol控制早停灵敏度设成 1e-4 通常够用如果你发现训练曲线抖动可以把tol调到 1e-3让模型更早停下来。如果想自动搜索参数可以用 MATLAB 优化工具箱里的surrogateopt把keep和max_layers当作整数变量处理。不过 GMDH 训练本身已经很快很多时候直接做两层网格搜索反而更直观。5.2 用训练历史判断过拟合GMDH 不像神经网络那样能直接画 loss 曲线但hist里已经保存了每层最小验证 RMSE。最直接的诊断方法是把hist第二列画出来figure; plot(hist(:, 2), o-); xlabel(层数); ylabel(最小验证 RMSE);如果曲线先下降、后上升说明网络在某一层之后开始记忆验证集噪声。gmdh_train里已经通过best_layer自动回退到最优层所以最终预测用的不是最后一层而是历史上验证误差最小的那一层。这个设计能抵消一部分过拟合但不能完全替代验证集划分的合理性。更稳妥的做法是外层再做一次 k 折交叉验证用每次折的验证误差重新选keep和max_layers。GMDH 训练开销小交叉验证代价通常可以接受。5.3 哪些场景不适合 GMDHGMDH 并不适合所有多输入回归问题。输入特征来自图像像素、文本向量或高频序列时原始特征之间缺乏显式可组合的交互结构GMDH 的二次多项式候选会变得很笨重这时用卷积网络或循环网络更合理。另一个边界是特征数量特别大的场景。当 p 超过 100 时即便加max_pairs限流随机抽样也会带来不确定性这时需要先做一遍特征筛选或者先用 LASSO 压缩到二三十个变量再交给 GMDH。还有GMDH 本身不擅长处理缺失值和类别型变量类别特征必须先做 One-hot 或嵌入编码。6. 让 GMDH 在 MATLAB 里更稳特征解释、保存与再训练6.1 用扰动法看特征重要性GMDH 中间层保存的是变量组合关系直接观察第一层的最佳节点只能看到两个中间变量未必能映射回原始特征。更实用的做法是扰动法把验证集某一列打乱观察预测 RMSE 变化多少。变化越大说明该特征对模型输出的影响越大。y0 gmdh_predict(model, Xva_n); imp zeros(1, p); for f 1:p Xp Xva_n; Xp(:, f) Xp(randperm(n), f); yp gmdh_predict(model, Xp); imp(f) sqrt(mean((yp - y0).^2)); end bar(imp); xlabel(特征编号); ylabel(扰动后 RMSE 变化);由于每次扰动都走完整的前向网络这个结果比查看单一系数更可靠。注意特征被完全打乱后如果模型本来就很少使用它RMSE 变化会接近 0如果变化很大说明该特征参与了多个层级的组合是模型的主要依赖变量。6.2 保存模型时必须连同标准化参数一起存实际项目里模型训练完通常要落盘到.mat文件。GMDH 的预测流程依赖训练时生成的层次结构也依赖mu_x、sd_x、mu_y、sd_y这四组统计量。只保存model而丢掉标准化参数是新手最容易踩的坑。save(gmdh_model.mat, model, mu_x, sd_x, mu_y, sd_y);部署预测时S load(gmdh_model.mat); x_new (new_row - S.mu_x) ./ S.sd_x; y_new gmdh_predict(S.model, x_new) * S.sd_y S.mu_y;这样保存的模型和参数绑定在一起换一台机器或换一个 MATLAB 版本加载后都能直接预测。GMDH 的模型体积很小中间层存储的无非是变量对索引和六个系数比同等精度神经网络更适合放进嵌入式环境或实时评分服务。本文还有配套的精品资源点击获取