Matlab实现贝叶斯全局优化：高斯过程与采集函数实战-编程实验室

简介：本资源是一套面向本硕博阶段科研与工程学习者的贝叶斯全局优化实践材料，聚焦高斯过程建模与采集函数（如UCB、EI）驱动的黑箱函数优化问题，适用于超参数调优、实验设计、仿真系统寻优等典型应用场景。压缩包共11个文件，含9个核心MATLAB函数（实现高斯过程拟合、边际似然计算、采集函数求解及参数最大化等关键模块）、1个说明文本和1个全程操作录屏AVI视频，整体仅173KB，轻量易部署。已有2025人下载学习，配套视频详细演示Runme.m主入口运行流程、路径设置要点及各子模块协同逻辑，避免常见运行错误；所有代码经MATLAB 2021a及以上版本实测，结构清晰、注释完备，便于理解贝叶斯优化的数学原理与工程实现细节。

1. 项目概述：当贝叶斯遇上全局优化

在工程优化、机器学习调参、甚至是实验设计里，我们常常会碰到一个让人头疼的问题：目标函数“黑盒化”。你只知道输入一些参数，然后得到一个输出结果（比如性能指标、成本、误差），但这个函数内部具体长什么样、有多少个坑（局部最优）、哪里是高峰（全局最优），你一无所知。更麻烦的是，每次评估这个“黑盒”函数都可能代价高昂——可能是跑一次仿真需要几小时，做一次物理实验要花几天，或者调用一次大型模型消耗大量计算资源。这时候，传统的优化方法，比如梯度下降（你得知道梯度）或者网格搜索（计算量爆炸），就显得力不从心了。

Bayesian全局优化（Bayesian Optimization, BO）就是为了解决这类问题而生的“聪明”策略。它的核心思想非常直观：与其盲目尝试，不如利用已有的观测数据，建立一个关于目标函数的概率模型，来描述我们对这个未知函数的“信念”。然后，基于这个模型，我们精心挑选下一个最有“希望”的评估点——这个点既要可能具有很高的函数值（ exploitation，利用已知的好区域），又要能帮助我们探索不确定性高的区域（ exploration，探索未知），以期更快地找到全局最优解。这个挑选下一个点的准则，就是所谓的“采集函数”。

在这个项目中，我们聚焦于使用高斯过程作为这个概率模型，在Matlab环境中实现一套完整的Bayesian全局优化仿真框架。高斯过程特别适合BO，因为它不仅能给出未知点的预测值（均值），还能给出预测的不确定性（方差），这正好为平衡“利用”与“探索”提供了量化依据。通过这个项目，你不仅能理解BO的理论之美，更能亲手实现它，并看到它如何一步步“学会”寻找复杂黑盒函数的最优点。

2. 核心原理：高斯过程与贝叶斯优化的默契配合

要理解这个项目，我们需要拆解两个核心部件：高斯过程回归模型，以及基于该模型的贝叶斯优化循环。

2.1 高斯过程：为未知函数建模

你可以把高斯过程想象成一个“函数的概率分布”。我们平常说的正态分布是描述一个随机变量取值的分布，而高斯过程则是描述一个函数 $f(x)$ 在所有输入点 $x$ 上取值的联合分布。它由均值函数 $m(x)$ 和协方差函数（核函数）$k(x, x’)$ 完全确定。

均值函数 $m(x)$：通常我们假设一个简单的常数均值（比如0），因为GP的灵活性主要来自核函数，数据标准化后这个假设通常是合理的。
核函数 $k(x, x’)$：这是GP的灵魂。它定义了函数在不同点 $x$ 和 $x’$ 之间的相似性。最常用的是平方指数核： $k(x, x’) = \sigma_f^2 \exp\left(-\frac{1}{2l^2} |x - x’|^2\right)$ 这里，$\sigma_f^2$ 是信号方差，控制函数整体的波动幅度；$l$ 是长度尺度，控制函数变化的“平滑度”。$l$ 越大，函数变化越缓慢，认为相距较远的点也更相关。

GP的预测过程：假设我们已经观测到一组数据 $\mathcal{D} = {X, y}$，其中 $X$ 是输入点矩阵，$y$ 是对应的观测值（可能含噪声）。对于一个新的输入点 $x_$，GP可以给出其函数值 $f_$ 的后验分布，这是一个正态分布： $p(f_* | X, y, x_) = \mathcal{N}(\mu_, \sigma_*^2)$

其中， $\mu_* = k(x_, X)[K(X, X) + \sigma_n^2 I]^{-1} y$ $\sigma_^2 = k(x_, x_) - k(x_, X)[K(X, X) + \sigma_n^2 I]^{-1} k(X, x_)$

这里，$K(X, X)$ 是所有观测点之间的协方差矩阵，$\sigma_n^2$ 是观测噪声的方差。$\mu_$ 就是我们对 $f(x_)$ 的最佳预测，而 $\sigma_*^2$ 则量化了这个预测的不确定性。这个“预测均值+预测方差”的二元输出，是后续优化决策的基础。

注意：核函数的选择和超参数 $(\sigma_f, l, \sigma_n)$ 的设定至关重要。在实际代码中，我们通常通过最大化边缘似然来优化这些超参数，让GP模型最好地拟合现有数据。

2.2 贝叶斯优化循环：序贯决策的艺术

有了GP这个强大的代理模型，BO就可以展开一个优雅的迭代循环：

初始化：在搜索空间内，随机选择或通过拉丁超立方采样选取少量初始点，评估目标函数，得到初始观测数据集 $\mathcal{D}$。
拟合GP模型：用当前数据集 $\mathcal{D}$ 拟合（训练）高斯过程模型，即优化其超参数。
最大化采集函数：利用训练好的GP模型，在整个搜索空间上计算采集函数 $a(x)$ 的值。然后，选择使 $a(x)$ 最大的点作为下一个评估点 $x_{next}$。 $x_{next} = \arg\max_{x \in \mathcal{X}} a(x)$
评估目标函数：在 $x_{next}$ 处运行昂贵的“黑盒”函数（在仿真中就是计算目标函数值），得到观测值 $y_{next}$。
更新数据集：将新的数据对 $(x_{next}, y_{next})$ 加入到观测集 $\mathcal{D}$ 中。
循环判断：如果未达到最大迭代次数或精度要求，则回到步骤2；否则，终止循环，并输出历史观测中找到的最佳点作为全局最优解的估计。

这个循环的核心在于第3步的采集函数。它负责平衡“利用”和“探索”。最常见的采集函数有：

期望改进：衡量新点比当前最优值改进的期望值。
上置信界：直接使用预测均值加上一个系数乘上预测标准差。这个系数控制了探索的积极性。
概率改进：衡量新点比当前最优值改进的概率。

在Matlab实现中，我们需要为这个循环的每一步编写清晰的函数或脚本模块。

3. Matlab仿真实现全流程拆解

下面，我们以一个经典的一维测试函数——Forrester函数为例，来详细拆解在Matlab中实现BO仿真的每一步。这个函数表达式为 $f(x) = (6x-2)^2\sin(12x-4)$，定义域 $x \in [0, 1]$。它有一个全局最小值和一个局部最小值，非常适合演示BO的寻优能力。

3.1 环境准备与目标函数定义

首先，我们需要一个干净的Matlab工作环境。建议为这个项目单独创建一个文件夹。

% 清空环境 clear; close all; clc; % 添加当前文件夹到路径，方便调用自定义函数 addpath(genpath(pwd));

接下来，定义我们的目标函数。在真实场景中，这里应该替换成你的“黑盒”仿真或实验接口。

function y = forrester_function(x) % Forrester 函数，常用于贝叶斯优化测试 % 输入 x 可以是标量或向量 y = (6*x - 2).^2 .* sin(12*x - 4); end

为了可视化，我们先画出这个函数的真实形状。

% 定义搜索空间 bounds = [0, 1]; % 生成密集的点用于绘图 x_plot = linspace(bounds(1), bounds(2), 1000)'; y_plot = forrester_function(x_plot); figure(‘Position‘， [100, 100, 800, 400]); plot(x_plot, y_plot, ‘b-‘, ‘LineWidth‘, 2); xlabel(‘输入 x‘); ylabel(‘输出 f(x)‘); title(‘Forrester 函数 (真实形状)‘); grid on; hold on;

3.2 高斯过程回归模型的实现

Matlab的统计和机器学习工具箱提供了fitrgp函数来拟合高斯过程回归模型，这大大简化了我们的工作。但为了理解原理，我们也可以手动实现核心部分。这里我们先使用工具箱函数。

步骤1：生成初始观测数据我们通常使用空间填充设计，如拉丁超立方采样，来获取有代表性的初始点。

% 设置随机种子保证可重复性 rng(42); % 初始点数量 n_init = 5; % 使用拉丁超立方采样生成初始点 X_init = lhsdesign(n_init, 1); % 生成[0,1]区间的样本 X_init = bounds(1) + (bounds(2)-bounds(1)) * X_init; % 映射到实际边界 % 评估目标函数（模拟昂贵操作） y_init = forrester_function(X_init); % 在图上标出初始点 plot(X_init, y_init, ‘ro‘, ‘MarkerSize‘, 10, ‘MarkerFaceColor‘, ‘r‘); legend(‘真实函数‘, ‘初始观测点‘);

步骤2：定义并拟合GP模型我们将使用平方指数核函数。

% 使用 fitrgp 拟合高斯过程模型 % ‘Basis‘, ‘constant‘: 使用常数基函数（均值函数） % ‘KernelFunction‘, ‘squaredexponential‘: 使用平方指数核 % ‘FitMethod‘, ‘exact‘: 使用精确推断（适用于数据量不大的情况） % ‘PredictMethod‘, ‘exact‘: 使用精确预测 % ‘Standardize‘, 1: 标准化输出，有助于数值稳定 gp_model = fitrgp(X_init, y_init, ... ‘Basis‘, ‘constant‘, ... ‘KernelFunction‘, ‘squaredexponential‘, ... ‘FitMethod‘, ‘exact‘, ... ‘PredictMethod‘, ‘exact‘, ... ‘Standardize‘, 1); disp(‘GP模型拟合完成。‘); disp([‘长度尺度参数 l: ‘, num2str(gp_model.KernelInformation.KernelParameters(1))]); disp([‘信号标准差 σ_f: ‘, num2str(sqrt(gp_model.KernelInformation.KernelParameters(2)))]); disp([‘噪声标准差 σ_n: ‘, num2str(gp_model.Sigma)]);

步骤3：利用GP进行预测和可视化拟合好模型后，我们可以用它来预测整个搜索空间上的函数行为。

% 在绘图点上进行预测 [ypred, ysd, yint] = predict(gp_model, x_plot); % 可视化GP的预测 figure(‘Position‘， [100, 100, 1000, 400]); subplot(1,2,1); plot(x_plot, y_plot, ‘b-‘, ‘LineWidth‘, 1.5, ‘DisplayName‘, ‘真实函数‘); hold on; plot(X_init, y_init, ‘ro‘, ‘MarkerSize‘, 10, ‘MarkerFaceColor‘, ‘r‘, ‘DisplayName‘, ‘观测点‘); plot(x_plot, ypred, ‘k--‘, ‘LineWidth‘, 2, ‘DisplayName‘, ‘GP预测均值‘); % 绘制置信区间（通常取均值 ± 2*标准差） fill([x_plot; flipud(x_plot)], [yint(:,1); flipud(yint(:,2))], ... ‘k‘, ‘FaceAlpha‘, 0.2, ‘EdgeColor‘, ‘none‘, ‘DisplayName‘, ‘95% 置信区间‘); xlabel(‘x‘); ylabel(‘f(x)‘); title(‘高斯过程回归拟合结果 (初始阶段)‘); legend(‘Location‘, ‘best‘); grid on; % 绘制预测标准差（不确定性） subplot(1,2,2); plot(x_plot, ysd, ‘m-‘, ‘LineWidth‘, 2); xlabel(‘x‘); ylabel(‘预测标准差‘); title(‘GP预测的不确定性‘); grid on;

可以看到，在观测点附近，GP的预测非常自信（不确定性低，置信区间窄）；而在远离观测点的区域，不确定性很高。这正是探索的潜在区域。

3.3 采集函数的选择与最大化

接下来是实现BO的核心——采集函数。我们以实现期望改进为例。EI的公式为： $EI(x) = (\mu(x) - f(x^+) - \xi)\Phi(Z) + \sigma(x)\phi(Z)$，如果 $\sigma(x) > 0$ 其中，$f(x^+)$ 是当前最佳观测值，$\Phi$ 和 $\phi$ 分别是标准正态分布的累积分布函数和概率密度函数，$Z = \frac{\mu(x) - f(x^+) - \xi}{\sigma(x)}$，$\xi$ 是一个小的正数，用于平衡探索。

function ei = expected_improvement(x, gp_model, y_best, xi) % 计算期望改进采集函数 % x: 待评估的点 (可以是矩阵，每行一个点) % gp_model: 拟合好的GP模型 % y_best: 当前最佳观测值 f(x+) % xi: 探索参数，默认为0.01 if nargin < 4 xi = 0.01; end [mu, sigma] = predict(gp_model, x); sigma = max(sigma, 1e-12); % 避免除零错误 z = (mu - y_best - xi) ./ sigma; ei = (mu - y_best - xi) .* normcdf(z) + sigma .* normpdf(z); % 对于 sigma 为0的点，EI设为0 ei(sigma <= 1e-12) = 0; end

现在，我们需要在搜索空间内找到使EI最大的点。对于一维问题，我们可以简单地在密集网格上计算；对于高维问题，则需要使用全局优化器（如fmincon配合多起点）。

% 找到当前最佳观测值 [y_best, idx_best] = min(y_init); x_best = X_init(idx_best); % 在密集网格上计算EI ei_values = expected_improvement(x_plot, gp_model, y_best, 0.01); % 找到使EI最大的点 [ei_max, idx_ei_max] = max(ei_values); x_next = x_plot(idx_ei_max); % 可视化EI函数和下一个建议点 figure; yyaxis left; plot(x_plot, ypred, ‘k--‘, ‘LineWidth‘, 1.5, ‘DisplayName‘, ‘GP预测‘); hold on; plot(X_init, y_init, ‘ro‘, ‘MarkerSize‘, 10, ‘MarkerFaceColor‘, ‘r‘, ‘DisplayName‘, ‘观测点‘); errorbar(x_plot, ypred, 2*ysd, ‘k:‘, ‘HandleVisibility‘, ‘off‘); % 置信区间 ylabel(‘f(x) / EI(x)‘); yyaxis right; plot(x_plot, ei_values, ‘g-‘, ‘LineWidth‘, 2, ‘DisplayName‘, ‘期望改进 (EI)‘); plot(x_next, ei_max, ‘g*‘, ‘MarkerSize‘, 15, ‘LineWidth‘, 2, ‘DisplayName‘, ‘下一个建议点‘); xlabel(‘x‘); ylabel(‘EI(x)‘); title(‘采集函数 (期望改进) 与下一个评估点‘); legend(‘Location‘, ‘best‘); grid on;

3.4 整合完整的贝叶斯优化循环

将以上步骤整合到一个循环中，就构成了完整的BO算法。

% 贝叶斯优化主循环参数设置 max_iter = 20; % 最大迭代次数（不含初始点） X_obs = X_init; % 观测点集合 y_obs = y_init; % 观测值集合 y_best_history = zeros(max_iter+1, 1); % 记录每次迭代后的历史最佳值 y_best_history(1) = y_best; x_best_history = zeros(max_iter+1, 1); x_best_history(1) = x_best; % 为动画或过程记录做准备 figure(‘Position‘， [100, 100, 1200, 500]); for iter = 1:max_iter % 1. 用所有现有数据重新拟合GP模型 gp_model = fitrgp(X_obs, y_obs, ... ‘Basis‘, ‘constant‘, ... ‘KernelFunction‘, ‘squaredexponential‘, ... ‘FitMethod‘, ‘exact‘, ... ‘PredictMethod‘, ‘exact‘, ... ‘Standardize‘, 1); % 2. 计算当前最佳观测值 [current_best, idx] = min(y_obs); current_best_x = X_obs(idx); % 3. 在整个搜索空间上计算采集函数（这里用EI） % 对于高维问题，这一步需要用优化器，这里演示仍用网格搜索 ei_vals = expected_improvement(x_plot, gp_model, current_best, 0.01); [~, idx_next] = max(ei_vals); x_next = x_plot(idx_next); % 4. 评估目标函数（模拟昂贵评估） y_next = forrester_function(x_next); % 5. 更新观测数据集 X_obs = [X_obs; x_next]; y_obs = [y_obs; y_next]; % 6. 更新历史最佳记录 if y_next < current_best y_best_history(iter+1) = y_next; x_best_history(iter+1) = x_next; else y_best_history(iter+1) = current_best; x_best_history(iter+1) = current_best_x; end % --- 可视化当前迭代状态 --- clf; % 子图1：函数真实值、GP预测、观测点、下一个点 subplot(2, 3, [1, 2, 4, 5]); [ypred, ysd] = predict(gp_model, x_plot); plot(x_plot, y_plot, ‘b-‘, ‘LineWidth‘, 1, ‘DisplayName‘, ‘真实函数‘); hold on; plot(X_obs(1:end-1), y_obs(1:end-1), ‘ko‘, ‘MarkerSize‘, 8, ‘MarkerFaceColor‘, ‘k‘, ‘DisplayName‘, ‘历史观测点‘); plot(X_obs(end), y_obs(end), ‘ro‘, ‘MarkerSize‘, 10, ‘MarkerFaceColor‘, ‘r‘, ‘DisplayName‘, ‘新评估点‘); plot(x_plot, ypred, ‘k--‘, ‘LineWidth‘, 1.5, ‘DisplayName‘, ‘GP预测均值‘); fill([x_plot; flipud(x_plot)], [ypred-2*ysd; flipud(ypred+2*ysd)], ... ‘k‘, ‘FaceAlpha‘, 0.1, ‘EdgeColor‘, ‘none‘, ‘DisplayName‘, ‘±2σ 区间‘); xlabel(‘x‘); ylabel(‘f(x)‘); title(sprintf(‘贝叶斯优化迭代 %d / %d‘, iter, max_iter)); legend(‘Location‘, ‘best‘); grid on; xlim(bounds); % 子图2：采集函数EI subplot(2,3,3); plot(x_plot, ei_vals, ‘g-‘, ‘LineWidth‘, 2); hold on; plot(x_next, ei_vals(idx_next), ‘g*‘, ‘MarkerSize‘, 15, ‘LineWidth‘, 2); xlabel(‘x‘); ylabel(‘EI(x)‘); title(‘期望改进采集函数‘); grid on; xlim(bounds); % 子图3：历史最佳值收敛情况 subplot(2,3,6); plot(0:iter, y_best_history(1:iter+1), ‘b-o‘, ‘LineWidth‘, 1.5, ‘MarkerFaceColor‘, ‘b‘); xlabel(‘迭代次数‘); ylabel(‘最佳观测值‘); title(‘最优值收敛历史‘); grid on; drawnow; pause(0.5); % 暂停一下以便观察动画效果 end % 输出最终结果 [final_best_y, final_idx] = min(y_obs); final_best_x = X_obs(final_idx); fprintf(‘\n======= 贝叶斯优化完成 =======\n‘); fprintf(‘总评估次数: %d\n‘, length(y_obs)); fprintf(‘找到的最优点 x*: %.6f\n‘, final_best_x); fprintf(‘对应的最优值 f(x*): %.6f\n‘, final_best_y); % 对于Forrester函数，理论最小值在 x≈0.757 处，f≈ -6.020 fprintf(‘理论最小值 f(x) ≈ -6.020 (x≈0.757)\n‘);

运行这段代码，你将看到一个动态的优化过程。GP模型随着新数据的加入而不断更新，采集函数引导着评估点从高不确定性区域（探索）逐渐聚焦到可能的最优区域（利用）。最终，算法应该能以远少于网格搜索的次数，逼近函数的全局最优点。

4. 关键参数调优与实操心得

一个鲁棒的BO实现，离不开对几个关键参数的深入理解和调优。这些参数直接影响着优化的效率和效果。

4.1 核函数超参数的优化

在之前的代码中，我们使用了fitrgp的默认超参数优化。但理解其原理很重要。通常，我们通过最大化对数边缘似然来优化超参数 $\theta = (l, \sigma_f, \sigma_n)$。

$\log p(y|X, \theta) = -\frac{1}{2}y^T(K+\sigma_n^2I)^{-1}y - \frac{1}{2}\log|K+\sigma_n^2I| - \frac{n}{2}\log 2\pi$

fitrgp内部就是通过最小化负对数边缘似然来实现的。手动实现这个优化可以帮助我们设置边界，防止过拟合或欠拟合。

% 手动设置优化选项示例 gp_model_manual = fitrgp(X_obs, y_obs, ... ‘Basis‘, ‘constant‘, ... ‘KernelFunction‘, ‘squaredexponential‘, ... ‘FitMethod‘, ‘exact‘, ... ‘PredictMethod‘, ‘exact‘, ... ‘Standardize‘, 1, ... ‘Optimizer‘, ‘quasinewton‘, ... % 使用拟牛顿法优化 ‘OptimizeHyperparameters‘, ‘auto‘, ... % 自动优化所有超参数 ‘HyperparameterOptimizationOptions‘, ... struct(‘AcquisitionFunctionName‘, ‘expected-improvement-per-second-plus‘, ... ‘MaxObjectiveEvaluations‘, 50)); % 对超参数进行贝叶斯优化！

实操心得：对于复杂或高维问题，GP超参数的优化本身可能陷入局部最优。一个技巧是使用多起点初始化来优化超参数。另一个常见问题是数值不稳定，特别是当观测点非常接近时，协方差矩阵 $K$ 可能接近奇异。添加一个小的“抖动”到矩阵对角线（如K = K + 1e-6 * eye(n)）可以保证其正定性。Matlab的fitrgp通常能很好地处理数值问题。

4.2 采集函数参数的选择

以EI函数中的 $\xi$ 参数为例。它控制着探索的倾向性：

$\xi = 0$：倾向于在已有好点附近做局部改进（偏向利用）。
$\xi > 0$：更愿意探索预测不确定性高但均值不一定最好的区域（偏向探索）。

一个动态调整的策略是在优化初期使用较大的 $\xi$（鼓励探索），随着迭代进行逐渐减小 $\xi$（聚焦利用）。例如：

function xi = get_dynamic_xi(iter, max_iter) % 动态调整EI的xi参数 xi_start = 0.1; xi_end = 0.001; xi = xi_start * (xi_end/xi_start)^(iter/max_iter); end

对于UCB采集函数 $a(x) = \mu(x) + \beta \sigma(x)$，参数 $\beta$ 扮演着类似的角色。理论上有一些渐进最优的 $\beta$ 设置，但在实践中，将其设置为一个常数（如2.0或5.0）或根据迭代次数衰减，通常效果也不错。

4.3 处理高维与复杂约束问题

当输入维度增加时，BO面临“维数灾难”。搜索空间呈指数增长，GP模型在高维下的拟合和优化采集函数都变得异常困难。

降维与特征选择：如果可能，利用领域知识进行降维。
使用更适合高维的核函数：如ARD核，它为每个输入维度分配独立的长度尺度，可以自动进行特征重要性排序。
```
gp_model_ard = fitrgp(X, y, ‘KernelFunction‘, ‘ardsquaredexponential‘);
```
优化采集函数的策略：在高维空间全局最大化采集函数本身就是一个难题。可以使用多起点局部优化、随机采样结合局部优化、甚至基于梯度的优化（如果采集函数可导）。

对于带约束的优化问题，例如 $g(x) \leq 0$，我们需要将约束信息融入BO框架。一种常见方法是约束概率：用一个独立的GP模型来建模约束函数 $g(x)$，然后计算一个点可行的概率 $p(g(x) \leq 0)$。最终的采集函数可以设计为 $a(x) = EI(x) \times p(\text{feasible})$，即只考虑可行区域内的期望改进。

5. 性能评估、对比与结果分析

为了证明BO的有效性，我们需要将其与基准方法进行对比。

5.1 对比方法：随机搜索与网格搜索

我们实现两种朴素的全局优化方法作为基准。

% 随机搜索 function [x_best_rs, y_best_rs, history_rs] = random_search(obj_func, bounds, n_eval) dim = size(bounds, 1); x_best_rs = []; y_best_rs = inf; history_rs.best_val = zeros(n_eval, 1); history_rs.all_x = zeros(n_eval, dim); history_rs.all_y = zeros(n_eval, 1); for i = 1:n_eval x = bounds(:,1) + (bounds(:,2)-bounds(:,1)) .* rand(dim, 1); y = obj_func(x); history_rs.all_x(i, :) = x‘; history_rs.all_y(i) = y; if y < y_best_rs y_best_rs = y; x_best_rs = x; end history_rs.best_val(i) = y_best_rs; end end % 网格搜索（以一维为例） function [x_best_gs, y_best_gs, history_gs] = grid_search(obj_func, bounds, n_grid) x_grid = linspace(bounds(1), bounds(2), n_grid)‘; y_grid = obj_func(x_grid); [y_best_gs, idx] = min(y_grid); x_best_gs = x_grid(idx); % 网格搜索的历史就是按顺序评估所有点 history_gs.best_val = cummin(y_grid); % 累积最小值 end

5.2 运行对比实验

我们使用相同的总评估次数（例如30次，含5个初始点）来比较BO、随机搜索和网格搜索。

% 定义实验参数 total_eval = 30; n_init = 5; bounds = [0; 1]; obj_func = @forrester_function; % 运行贝叶斯优化 (使用我们之前实现的循环，但控制总评估次数) % ... (这里嵌入之前的主循环代码，将 max_iter 设置为 total_eval - n_init) % 假设运行后得到 history_bo % 运行随机搜索 [x_best_rs, y_best_rs, history_rs] = random_search(obj_func, bounds‘, total_eval); % 运行网格搜索 n_grid = total_eval; % 网格点数等于总评估次数 [x_best_gs, y_best_gs, history_gs] = grid_search(obj_func, bounds‘, n_grid); % 绘制收敛曲线对比图 figure(‘Position‘， [100, 100, 900, 500]); plot(1:total_eval, history_bo.best_val, ‘b-o‘, ‘LineWidth‘, 2, ‘MarkerFaceColor‘, ‘b‘, ‘DisplayName‘, ‘贝叶斯优化‘); hold on; plot(1:total_eval, history_rs.best_val, ‘r-s‘, ‘LineWidth‘, 1.5, ‘MarkerFaceColor‘, ‘r‘, ‘DisplayName‘, ‘随机搜索‘); plot(1:total_eval, history_gs.best_val, ‘g-^‘, ‘LineWidth‘, 1.5, ‘MarkerFaceColor‘, ‘g‘, ‘DisplayName‘, ‘网格搜索‘); xlabel(‘函数评估次数‘); ylabel(‘当前找到的最佳值‘); title(‘优化方法收敛曲线对比 (Forrester函数)‘); legend(‘Location‘, ‘best‘); grid on; set(gca, ‘YScale‘, ‘log‘); % 使用对数坐标更清晰地显示差距

5.3 结果分析与解读

运行对比实验后，你通常会看到：

贝叶斯优化：曲线初期下降很快，能在很少的评估次数内找到接近最优的解。这是因为它智能地利用了历史信息来指导搜索。
随机搜索：曲线呈阶梯状下降，下降速度慢于BO，但最终随着评估次数增加也可能找到不错解。其性能方差较大。
网格搜索：曲线呈单调下降（因为累积最小值），但下降速度取决于网格的密度。在总评估次数固定时，其最终解的质量通常不如BO，甚至可能不如好的随机搜索，因为它缺乏对函数形状的适应性。

关键指标：

最终解质量：在相同评估预算下，BO找到的y_best通常最小（对于最小化问题）。
收敛速度：BO的曲线最早“触底”，说明其采样效率最高。
稳健性：可以多次运行随机搜索和BO（随机初始点不同），观察它们性能的方差。BO通常更稳定。

注意事项：BO的优势在评估代价高昂时最为明显。如果函数评估本身很快（比如毫秒级），那么简单的网格或随机搜索可能更简单直接。BO的计算开销主要来自GP模型的拟合和采集函数的优化，当观测数据超过几千个点时，精确GP的计算复杂度 $O(n^3)$ 会成为瓶颈，此时需要考虑稀疏GP等近似方法。

6. 常见问题排查与实战技巧

在实际实现和应用BO时，你可能会遇到以下典型问题。

6.1 数值不稳定与矩阵病态

问题：在拟合GP时，Matlab报错，提示协方差矩阵接近奇异或不是正定矩阵。原因：观测点中有两个点过于接近，导致协方差矩阵的行/列几乎线性相关。解决方案：

在调用fitrgp时，可以尝试增加‘Sigma‘参数的初始值或下界，这相当于强制增加一个观测噪声。
```
gp_model = fitrgp(X, y, ‘Sigma‘, 1e-4, ...);
```
或者在核函数中添加一个小的“白噪声”项（fitrgp会自动处理噪声参数Sigma）。
检查数据中是否有重复或极度接近的点，必要时进行去重或微扰。

6.2 优化陷入局部最优

问题：BO运行多次，但似乎总是收敛到同一个局部最优点，无法找到全局最优。原因：

初始点太少或位置不好，导致GP模型对整个空间的认知有偏。
采集函数的探索性不足（如EI的 $\xi$ 太小）。
搜索空间定义可能不正确，全局最优根本不在定义的范围内。解决方案：
增加初始点数量，并使用拉丁超立方采样确保初始点能较好地覆盖整个空间。
调整采集函数：在初期使用探索性更强的采集函数（如UCB，并设置较大的 $\beta$），或使用动态衰减的 $\xi$。
多次独立运行：从不同的随机初始点集开始运行BO，取多次运行中的最好结果。
如果问题维度不高，可以尝试在最大化采集函数这一步使用全局优化器（如patternsearch,ga）而非局部优化器，确保能找到采集函数的全局极大点。

6.3 高维问题性能下降

问题：当输入维度超过10维时，BO的优化效率急剧下降，甚至不如随机搜索。原因：这就是“维数灾难”。高维空间中，任何点的邻域都是近乎空白的，GP难以建立有效的相关性模型。解决方案：

使用ARD核：让模型自动学习每个维度的重要性。
主动降维：使用主成分分析或领域知识，将问题投影到更低维的空间进行优化。
考虑可加性结构：如果目标函数可以近似为多个低维函数的和，可以使用可加GP模型。
采用更高效的采集函数优化策略，如随机采样配合局部搜索。

6.4 与实际问题对接

问题：我的目标函数是一个需要运行1小时的仿真脚本，如何集成？解决方案：

封装函数接口：将你的仿真脚本封装成一个Matlab函数，该函数以参数向量为输入，输出一个标量性能指标（需要最小化或最大化）。确保该函数能处理可能的错误（如仿真不收敛），并返回一个惩罚值（如一个很大的数）。
异步与并行：标准的BO是串行的。但你可以修改采集函数，一次建议多个评估点（如通过q-EI），然后利用计算集群并行运行这些仿真。Matlab的并行计算工具箱可以帮助你实现这一点。
设置超时与容错：在BO主循环中，对目标函数的调用进行try-catch包装，并设置超时限制。如果评估失败，可以将其视为一个“糟糕”的值（例如，对于最小化问题，返回一个非常大的数），并继续优化流程。记录失败案例以供分析。

一个实用的调试技巧：在将BO应用到昂贵的真实仿真之前，先用一个快速代理模型（比如基于少量数据训练的简单神经网络或多项式模型）来测试你的整个BO流程。确保代码逻辑正确，参数设置合理，能够在这个廉价代理上找到最优解。这可以节省大量时间和计算资源。

本文还有配套的精品资源，点击获取