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

Kriging这个名字听起来很有学术感，但实际它的定位特别朴素：用少量仿真样本训练一个低成本近似模型，去替代那些动不动就跑几小时的高精度仿真。我做结构优化和参数标定时，最头疼的不是优化算法选型，而是“一次仿真三小时起步，优化算法跑两百次”的成本账。后来把Kriging代理模型嵌进流程，真实仿真调用次数砍到原来的五分之一，总耗时从几周压缩到一天以内。这篇文章就聊聊如何在MATLAB里搭建Kriging代理模型，从手写核心代码到内置fitrgp落地，再到用它做贝叶斯优化，最后分享几个实践中绕不开的坑。适合正在做仿真优化、代理模型研究或设计空间探索的工程师和研究生参考。

1. 为什么选Kriging当代理模型：从“仿真太贵”到“插值推理”

1.1 代理模型怎么就成了“仿真替代品”

仿真成本高是优化问题里最常见的痛点。以CFD和有限元分析为例，一次计算可能花费几小时甚至整夜，而参数优化、容差分析、可靠性评估这些工作天然需要“反复问模型要答案”。如果每次都调用原仿真模型，计算资源根本扛不住。

代理模型的思路很简单：先用少量样本点（比如拉丁超立方采样生成的三十组参数组合）调用真实仿真，拿到一批输入输出数据，然后训练一个计算成本极低的数学模型，用它来预测任意新参数组合下的输出。这个模型不需要完美还原仿真的物理细节，只要趋势和关键响应足够准确，就能支撑后续的优化搜索。

Kriging在这里的特殊优势在于它不只是给一个预测值，还顺带给出预测方差——也就是“这个位置的预测有多大把握”。这个不确定性信息非常宝贵，后面讲EGO采集函数时会看到它如何被用来平衡探索与利用。简单说，Kriging天然是给“代价高昂的黑箱函数”设计的代理模型。

1.2 Kriging和其他常见代理模型的对标

做代理模型有很多选择，响应面、径向基、支持向量回归都能干。我整理了一张常用代理模型的对比表，方便你判断什么时候该选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) * y

R是训练样本之间的相关矩阵，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当成超参数一起优化，这意味着得到的模型不会严格穿过训练点，而是留了一点噪声平滑。如果你面对的是确定性仿真数据，希望代理模型精确插值，可以手动固定一个很小的Sigma：

gprMdl = 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*n

sigma2是过程方差，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价值的场景是把它嵌入到优化流程里。经典的EGO（Efficient 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这种封装好的工具。这样即使遇到封装函数抛出的看不懂的报错，也能从底层逻辑判断问题出在数据上、相关函数上还是超参数优化上。