MATLAB中Kriging代理模型实战：从手写代码到贝叶斯优化

Kriging这个名字听起来很有学术感但实际它的定位特别朴素用少量仿真样本训练一个低成本近似模型去替代那些动不动就跑几小时的高精度仿真。我做结构优化和参数标定时最头疼的不是优化算法选型而是“一次仿真三小时起步优化算法跑两百次”的成本账。后来把Kriging代理模型嵌进流程真实仿真调用次数砍到原来的五分之一总耗时从几周压缩到一天以内。这篇文章就聊聊如何在MATLAB里搭建Kriging代理模型从手写核心代码到内置fitrgp落地再到用它做贝叶斯优化最后分享几个实践中绕不开的坑。适合正在做仿真优化、代理模型研究或设计空间探索的工程师和研究生参考。1. 为什么选Kriging当代理模型从“仿真太贵”到“插值推理”1.1 代理模型怎么就成了“仿真替代品”仿真成本高是优化问题里最常见的痛点。以CFD和有限元分析为例一次计算可能花费几小时甚至整夜而参数优化、容差分析、可靠性评估这些工作天然需要“反复问模型要答案”。如果每次都调用原仿真模型计算资源根本扛不住。代理模型的思路很简单先用少量样本点比如拉丁超立方采样生成的三十组参数组合调用真实仿真拿到一批输入输出数据然后训练一个计算成本极低的数学模型用它来预测任意新参数组合下的输出。这个模型不需要完美还原仿真的物理细节只要趋势和关键响应足够准确就能支撑后续的优化搜索。Kriging在这里的特殊优势在于它不只是给一个预测值还顺带给出预测方差——也就是“这个位置的预测有多大把握”。这个不确定性信息非常宝贵后面讲EGO采集函数时会看到它如何被用来平衡探索与利用。简单说Kriging天然是给“代价高昂的黑箱函数”设计的代理模型。1.2 Kriging和其他常见代理模型的对标做代理模型有很多选择响应面、径向基、支持向量回归都能干。我整理了一张常用代理模型的对比表方便你判断什么时候该选Kriging。代理模型类型预测均值不确定性估计小样本表现调参难度典型场景多项式响应面有无一般低低阶趋势分析径向基RBF有无较好低快速插值拟合支持向量回归SVR有无标准版一般中高维监督回归Kriging/高斯过程有有好中高仿真代理、贝叶斯优化Kriging最打动我的三点一是插值特性训练样本点上的预测值会精确回到仿真结果这对无噪声确定性仿真非常友好二是它能给出局部置信区间让后续优化算法知道哪里“还没探明白”三是它对样本量的容忍度比神经网络高很多三五十个样本就能搭起一个能用的模型而深度学习方法动辄需要上千数据点。当然Kriging也有短板。训练过程涉及协方差矩阵求逆和超参数优化样本量超过几千以后计算开销明显上升预测外推能力弱几乎不能指望它预测训练样本范围之外的行为对输入特征尺度敏感不同维度取值范围差异过大时相关长度参数很容易优化失败。这些短板不是劝退理由但一定要在设计流程时提前规避后面我会给出具体处理办法。2. 手写Kriging核心相关矩阵、克里金方程组与预测函数2.1 数学模型只用记三行很多教程把Kriging讲得很玄剥开来看核心其实就三件事设定相关函数、估计常数均值、解线性方程组。Kriging假设响应函数的形式为y(x) mu Z(x)其中mu是全局常数均值Z(x)是一个零均值的高斯过程。任意两点之间的协方差由相关函数决定最常用的是平方指数形式k(xi, xj) exp(-sum(theta .* (xi - xj).^2))theta就是相关长度参数它的每一维对应一个输入变量。theta越大两点相关性随距离衰减越快模型越“敏感”theta越小相关性衰减越慢模型越平滑。给定训练数据X和y之后普通Kriging的预测均值写成y_hat(x) mu r * R^(-1) * (y - 1 * mu) mu (1 * R^(-1) * 1)^(-1) * 1 * R^(-1) * yR是训练样本之间的相关矩阵r是新点x与训练样本之间的相关向量。预测方差为s^2(x) sigma^2 * [1 - r*R^(-1)*r (1 - 1*R^(-1)*r)^2 / (1*R^(-1)*1)]其中sigma^2用最大似然估计得到sigma^2 (y - 1*mu) * R^(-1) * (y - 1*mu) / n这套公式看着不复杂实际编程时最容易出错的地方是求解和数值稳定性的处理。直接写inv(R)是最糟糕的做法样本稍多或样本点稍近就会让矩阵接近奇异求逆结果直接飞掉。2.2 MATLAB核心实现普通Kriging下面的代码是我在MATLAB里手写Kriging时常用的版本去掉了花哨功能保留最核心的预测能力。为了代码可读性我没有做向量化加速但如果训模样本上千建议把相关矩阵改成bsxfun或者pdist2写法。function [yhat, s2, mu, sigma2] kriging_predict(x, Xtr, ytr, theta, nugget) % 普通Kriging预测x为单个新样本点Xtr为训练输入ytr为训练输出 % theta为相关长度向量nugget为数值稳定的小量 n size(Xtr, 1); R corr_mtx(Xtr, Xtr, theta) nugget * eye(n); % Cholesky分解R L * L L chol(R, lower); one ones(n, 1); a L \ one; % a L^(-1) * 1 b L \ ytr; % b L^(-1) * y mu (a * b) / (a * a); % 常数均值估计 res ytr - mu * one; c L \ res; lambda L \ c; % lambda R^(-1) * (y - mu*one) r corr_mtx(x, Xtr, theta); % 1 x n 相关向量 yhat mu r * lambda; sigma2 res * (L \ (L \ res)) / n; p L \ r; rRinvr p * p; oneRinvo a * a; oneRinvr a * p; s2 sigma2 * (1 - rRinvr (1 - oneRinvr)^2 / oneRinvo); s2 max(s2, 0); % 由于数值误差可能轻微为负 end function R corr_mtx(X1, X2, theta) % 平方指数相关函数 n1 size(X1, 1); n2 size(X2, 1); R zeros(n1, n2); for i 1:n1 for j 1:n2 d X1(i,:) - X2(j,:); R(i,j) exp(-sum(theta .* d.^2)); end end end用一段简单数据测一下rng(2025); Xtr lhsdesign(25, 2); ytr sin(3*Xtr(:,1)) .* cos(2*Xtr(:,2)) 0.1*Xtr(:,1).*Xtr(:,2); theta0 [0.1, 0.1]; [yhat, s2] kriging_predict([0.3, 0.7], Xtr, ytr, theta0, 1e-8);训练样本点的预测值会精确回到原仿真值这正是Kriging插值性的体现。如果你只需要预测均值不关心方差代码里最后几行方差计算可以划掉但做优化迭代时建议保留方差它是贝叶斯优化和加点策略的核心。2.3 为什么用Cholesky而不是inv我在初学Kriging时直接写了lambda inv(R) * (y - mu)结果在样本点分布稍密时预测曲线疯狂震荡检查半天才发现是矩阵求逆引入了不可接受的数值误差。相关矩阵R理论上是对称正定的但实际计算中由于浮点误差和样本点距离过近条件数可能爆炸到10的12次方以上。Cholesky分解先把R拆成L * L然后再用两次回代求解。这样做速度快、数值稳定而且不需要在代码里出现inv函数。另一个重要补充是nugget参数也就是在R的对角线上加一个很小的量比如1e-8到1e-6。它既能让矩阵更健康又相当于给模型加了一点点噪声容忍度。对于无噪声仿真数据nugget应尽可能小否则会破坏精确插值特性对于有噪声实验数据nugget本身就是一个需要优化的超参数。3. 用fitrgp快速落地参数配置、精度验证与实践边界3.1 fitrgp快速拟合代码如果你不想反复调试手写代码MATLAB自带的fitrgp就是高斯过程回归的工业级实现。fitrgp底层就是Kriging的高斯过程视角还内置了超参数优化、交叉验证、预测区间输出等功能。对于多数工程项目直接用它比手写版本稳得多。一个典型的最小代码示例rng(2025); Xtr lhsdesign(30, 2); ytr sin(3*Xtr(:,1)) .* cos(2*Xtr(:,2)) 0.05*randn(30,1); gprMdl fitrgp(Xtr, ytr, ... KernelFunction, ardsquaredexponential, ... Standardize, true, ... HyperparameterOptimizationOptions, struct( ... AcquisitionFunctionName, expected-improvement, ... MaxObjectiveEvaluations, 30)); Xeval lhsdesign(100, 2); [ypred, ysd] predict(gprMdl, Xeval);predict函数的第二个输出ysd就是预测标准差对应手写公式里的sqrt(s2)。如果想验证模型精度可以留一部分样本做测试或者使用fitrgp内置的交叉验证gprMdlCV fitrgp(Xtr, ytr, KFold, 5, Standardize, true); rmseCV kfoldLoss(gprMdlCV, Mode, average);kfoldLoss返回的是均方误差开方后就是RMSE。这个数值可以作为代理模型与真实仿真的平均偏差参考。3.2 关键参数取舍fitrgp的默认配置能跑通大多数问题但要得到可靠结果有几个参数建议手动过一遍。核函数选择上ardsquaredexponential是最常用的起点。ARD表示每个输入维度独立估计相关长度适合特征尺度不同、重要程度不同的场景。如果数据响应比较粗糙matern32或matern52会比平方指数更好Matérn核能控住样本点附近的光滑程度。Standardize必须设为true这相当于把输入输出缩放到统一尺度。别小看这一步自编Kriging时没做归一化导致收敛困难的例子我见过太多次了。FitMethod和PredictMethod默认值是exact适合样本量几百以内的场景。如果样本量上万exact计算会非常吃力可以切换为sd子集数据点方案但预测不确定性会打折。一个容易忽略的坑是噪声Sigma。fitrgp默认把Sigma当成超参数一起优化这意味着得到的模型不会严格穿过训练点而是留了一点噪声平滑。如果你面对的是确定性仿真数据希望代理模型精确插值可以手动固定一个很小的SigmagprMdl fitrgp(Xtr, ytr, ... KernelFunction, squaredexponential, ... Sigma, 1e-6, ... Standardize, true);这样做的好处是训练点预测误差基本为零坏处是如果仿真本身有数值噪声模型容易被带偏。判断依据很简单你的数据源头是实验测量还是纯仿真模型。4. 超参数优化实战初始值、边界、nugget与收敛判断4.1 负对数似然理解优化目标手写Kriging时最核心的工作是估计相关长度theta和噪声参数。fitrgp内部已经做了这步但如果你用自己写的预测函数就需要手动优化超参数。最常用的目标函数是负对数边际似然NLL。它衡量的是“在给定超参数下当前训练数据出现的概率有多大”。概率越大NLL越小超参数越好。省略常数项后可以写成nll n/2 * log(sigma2) sum(log(diag(L))) 0.5*nsigma2是过程方差L是相关矩阵的Cholesky下三角因子。用这种profile形式的好处是sigma2可以解析消掉只需要优化theta向量简化了问题。下面是一段可运行的目标函数代码function nll kriging_nll(logtheta, Xtr, ytr, nugget) theta exp(logtheta); % 在对数空间优化保证正数 n size(Xtr, 1); R corr_mtx_fast(Xtr, Xtr, theta) nugget * eye(n); L chol(R, lower); one ones(n, 1); a L \ one; b L \ ytr; mu (a * b) / (a * a); res ytr - mu * one; c L \ res; sigma2 res * (L \ (L \ res)) / n; nll n/2 * log(sigma2) sum(log(diag(L))) 0.5*n; end通过对数变换约束theta始终为正比直接在原空间加边界约束更稳。优化时可以用fminsearch也可以用fmincon加边界。4.2 调参流程与初始值经验超参数优化最怕两件事初始值离谱导致陷入局部最优以及theta冲到边界导致模型失效。我的实际操作流程基本固定为四步。第一步是归一化输入。所有输入特征缩放到[0, 1]区间这一步能让theta初始值有统一量纲。我常用的是Xnor (Xtr - min(Xtr)) ./ (max(Xtr) - min(Xtr));第二步是设定theta初始值。归一化之后theta0取0.1乘全一向量基本不会出大错。太大会让相关矩阵对角占优模型退化成“只认识样本点”太小会让相关矩阵接近全一矩阵模型退化成多项式回归。第三步是设定优化边界。我把theta的搜索范围放在0.001到100之间这是在对数空间下很宽的区间。如果优化结果落在边界上说明数据本身或特征选择可能有问题不是单纯调参能解决的。第四步是训练后验证。把优化得到的theta放回训练集看看交叉验证RMSE是否合理同时检查相关矩阵条件数cond(R)如果条件数超过1e10基本可以判断样本点存在近重复或过于密集的情况这时候nugget应往上调整或者对输入做去重。4.3 一个实际调试例子有一次数值标定问题输入是两个材料参数输出是一个响应指标。我用了二十五个样本点初始theta设[0.1, 0.1]fmincon优化后得到theta约为[0.37, 0.15]。看起来第二维相关长度更小代表第二个参数在较大距离上仍有较强相关性模型对第二个参数的变化更敏感。验证集RMSE是0.034相对于输出幅值0.4已经足够支撑后续优化。另一个项目里我偷懒没做归一化第一维范围是0到500第二维范围是0到1。theta的优化结果反复震荡NLL曲线锯齿状最后发现超参数把大部分注意力放在了第一维的尺度上第二维几乎被忽略。归一化之后同样的问题一次收敛。现在我做代理模型前会把归一化当成强制步骤而不是可选优化技巧。5. 不止拟合Kriging驱动的加点策略与贝叶斯优化5.1 从“预测”到“采集”EGO的思路拟合代理模型只是第一步真正发挥Kriging价值的场景是把它嵌入到优化流程里。经典的EGOEfficient Global Optimization思路是先做一批初始样本训练Kriging然后通过最大化采集函数挑选下一个最有价值的样本点去跑真实仿真再把结果加入训练集重新拟合Kriging如此循环。采集函数的代表是期望改进量EI。它把Kriging的预测均值和预测方差揉在一起形成一个关于“新点能比当前最优解好多少”的期望值。EI大意味着两种可能要么预测均值显著优于当前最优点要么预测方差很大代表这块区域还没探索明白。EGO会自然地在“探索未知区域”和“开发已知低点”之间做平衡。这种下一点选择逻辑非常聪明。它避免了一次性铺满整个设计空间的浪费也不需要人工指定每个迭代步在哪里采样。对高成本仿真来说通常二三十轮EGO迭代就能找到接近全局最优的解而直接跑遗传算法可能要消耗几百次仿真。这也解释了为什么Kriging几乎成了贝叶斯优化的默认代理模型。5.2 bayesopt一把梭如果你不想自己实现EI公式和加点循环MATLAB的bayesopt函数就是现成的EGO工业化实现。它内部使用高斯过程代理模型并提供不同的采集函数直接用起来省心很多。下面是用bayesopt做两变量黑箱函数最小化的完整示例f (x) expensiveSimulation(x.Ca, x.T); results bayesopt(f, ... [optimizableVariable(Ca, [0.5, 2.5]), ... optimizableVariable(T, [300, 500])], ... AcquisitionFunctionName, expected-improvement, ... MaxObjectiveEvaluations, 30, ... IsObjectiveDeterministic, true, ... Verbose, 1);这里的expensiveSimulation可以是调用真实仿真的函数。MaxObjectiveEvaluations设置最多调用真实仿真多少次相对于直接在优化器里跑几百次这已经是相当“省钱”的预算了。运行结束后用results.XAtMinObjective查看找到的最优参数组合用results.MinObjective查看最优目标值。还可以画一下results的曲线看看每轮迭代目标值的变化轨迹能直观感受到Kriging代理模型驱动优化时收敛有多快。需要注意一点bayesopt对目标函数是随机噪声还是确定性仿真有不同的推荐设置。如果是仿真结果没有随机波动把IsObjectiveDeterministic设为true会更贴合Kriging的插值假设如果目标是实验测量或有随机噪声保持默认false更稳妥。6. 实测中绕不开的坑与我的操作习惯6.1 病态矩阵是最隐蔽的罪魁祸首手写Kriging过程中我踩过最大的坑就是相关矩阵病态。症状非常典型预测结果在样本点附近剧烈振荡NLL怎么优化都不收敛或者优化出来的theta完全贴到下边界。病态的来源通常是样本点过于密集或者存在近似重复的样本。比如用LHS采样后两个点可能在高维空间的某几个维度上距离几乎为零相关矩阵的某两行就会近似线性相关。另一种情况是nugget设得太大相关矩阵退化成对角占优矩阵预测结果完全忽略邻近样本的影响。处理流程我建议按顺序排查先统计样本是否有重复或近重复有就删掉只保留其一然后检查数据是否归一化接着看条件数cond(R)超过1e8就调大nugget到1e-6量级最后再看优化结果是否贴边界。这四步走完九成病态问题都能解决。6.2 数据质量与样本数量代理模型的精度上限由数据决定。Kriging不是魔法如果初始样本没有覆盖设计空间的边缘它很难凭空预测边缘区域的极端行为。我见过不少人训练完模型只看验证集RMSE漂亮就直接拿去做优化结果最优解落在训练样本覆盖范围之外代理模型预测的是一个“外推值”完全失真。我的经验是做代理模型前先想清楚三个问题样本点是否在设计空间内均匀铺开、边界附近是否有样本点、目标响应的关键非线性区域是否有足够的样本密度。如果发现边界样本稀疏先补几个边界点再说不要急着训练。样本量方面从低维问题看n取输入维数的5到10倍通常可以建立基本可靠的Kriging模型也就是两三个变量用三十个左右样本起步。变量数超过十个以后建议先做敏感性分析筛掉一部分参数别指望Kriging在一个二十维问题上还能用几十个样本翻出浪花。另外要记住Kriging本质上是一种插值方法内插可信外推只能作为粗估绝不能作为设计决策依据。每次新增样本点前我都会看一眼训练样本落在设计空间的哪些位置确保每次仿真调用都花在真正“信息量大”的区域而不是重复验证模型已经熟知的地方。这个部分没有太多高深理论全是实操中反复踩出来的经验。如果你也在用MATLAB做代理模型相关项目建议先把手写版本调通理解每一步在算什么再切换到fitrgp和bayesopt这种封装好的工具。这样即使遇到封装函数抛出的看不懂的报错也能从底层逻辑判断问题出在数据上、相关函数上还是超参数优化上。