贝叶斯优化实战（四）（3）https://developer.aliyun.com/article/1516486
A.7 第八章：通过受限优化满足额外约束
本章有两个练习：
- 第一个验证我们从 BoTorch 对受限 EI 策略的实现得到的结果是否与常规 EI 分数和可行性概率的乘积相同。
- 第二部分向我们展示了如何在一个四维气动结构优化问题上运行受限 BayesOpt。
A.7.1 练习 1：受限 EI 的手动计算
受限 EI 策略的获取分数是 EI 分数和可行性概率的乘积。虽然 BoTorch 的 ConstrainedExpectedImprovement 类提供了受限 EI 分数的实现，但实际上我们可以手动执行计算。在这个练习中，我们探索这种手动计算，并将我们的结果与 ConstrainedExpectedImprovement 类的结果进行验证。此练习的解决方案在 CH08/02 - Exercise 1.ipynb 笔记本中，并可以解释如下：
- 重新创建 CH08/01 - Constrained optimization.ipynb 中使用的受限 BayesOpt 问题，包括目标函数、成本函数、GP 实现以及在一些训练数据上训练 GP 的辅助函数
fit_gp_model()。 - 创建一个 PyTorch 张量，它是在 -5 到 5 之间的密集网格。这个张量将作为我们的测试集。我们使用
torch.linspace()来创建一个密集网格：
lb = -5 ub = 5 xs = torch.linspace(lb, ub, 201)
- 通过从我们的搜索空间中随机抽样三个数据点（在 -5 到 5 之间）来创建一个玩具训练数据集，并在这些点上评估目标和成本函数。我们使用
torch.rand()在 0 和 1 之间随机采样，然后将样本缩放到我们的搜索空间：
n = 3 torch.manual_seed(0) ❶ train_x = bounds[0] + (bounds[1] - bounds[0]) * torch.rand(n) ❷ train_utility = objective(train_x) train_cost = cost(train_x)
- ❶ 为了可重现性而固定种子
❷ 在 0 和 1 之间进行采样，然后将样本缩放到我们的搜索空间 - 使用辅助函数
fit_gp_model()在目标函数数据和成本函数数据上训练一个 GP。
utility_model, utility_likelihood = fit_gp_model( ❶ train_x.unsqueeze(-1), train_utility ❶ ) ❶ cost_model, cost_likelihood = fit_gp_model( ❷ train_x.unsqueeze(-1), train_cost ❷ ) ❷
- ❶ 在目标函数数据上训练一个 GP
❷ 在成本函数的数据上训练一个 GP - 使用在成本函数数据上训练的 GP 来计算测试集中每个点的可行性概率。
我们首先计算成本 GP 在我们的测试集上的预测分布：
with torch.no_grad(): cost_pred_dist = cost_likelihood(cost_model(xs)) cost_pred_mean = cost_pred_dist.mean cost_pred_lower, cost_pred_upper = \ cost_pred_dist.confidence_region()
- 然后，我们初始化一个正态分布对象，其均值和标准差对应于
cost_pred_dist的均值和标准差：
normal = torch.distributions.Normal(cost_pred_mean, cost_pred_dist.stddev)
- 最后，我们在这个对象上调用
cdf()方法来计算可行性的概率。这个方法所取的参数是我们成本约束的上限，即 0：
feasible_prob = normal.cdf(torch.zeros(1))
- 初始化常规 EI 策略，其中
model参数是在目标函数数据上训练的 GP，而best_f参数是当前的可行入围者。
我们使用train_utility[train_cost <= 0].max()计算当前的可行入围者：
ei = botorch.acquisition.analytic.ExpectedImprovement( model=utility_model, best_f=train_utility[train_cost <= 0].max(), )
- 然后，通过在
xs[:, None, None]上调用 EI 策略对象来计算 EI 分数，该对象是为确保其形状适当而重新塑造的测试密集网格：
with torch.no_grad(): ei_score = ei(xs[:, None, None])
- 初始化受限 EI 策略，并为测试集中的每个点计算受限 EI 分数：
constrained_ei = botorch.acquisition.analytic.ConstrainedExpectedImprovement( model=botorch.models.model_list_gp_regression.ModelListGP( utility_model, cost_model ), best_f=train_utility[train_cost <= 0].max(), objective_index=0, constraints={1: [None, 0]} )
- 我们还使用重塑后的测试集计算受限 EI 分数：
with torch.no_grad(): constrained_ei_score = constrained_ei(xs[:, None, None])
- 计算 EI 分数和可行性概率的乘积，并验证这种手动计算是否与 BoTorch 的实现结果相同。 运行断言以确保所有相应的术语匹配：
assert torch.isclose( ei_score * feasible_prob, constrained_ei_score, atol=1e-3 ).all()
- 在图中绘制 EI 分数和受限 EI 分数，并直观地验证前者始终大于或等于后者。 证明这是正确的。
我们如下绘制迄今为止计算的分数：
plt.plot(xs, ei_score, label="EI") plt.plot(xs, constrained_ei_score, label="BoTorch constrained EI")
- 此代码生成图 A.15，显示 EI 分数确实始终至少等于受限 EI 分数。
图 A.15 EI 的获取分数（实线）和受限 EI（虚线）。 前者始终大于或等于后者。
我们可以通过注意到受限 EI 分数等于正常 EI 分数乘以可行性概率来数学证明这一点。 可行性概率始终最大为 1，因此 EI 分数始终大于或等于受限 EI 分数。
A.7.2 练习 2：飞机设计的受限优化
在这个练习中，我们使用第七章练习 2 中的飞机效用目标函数来解决一个受限制的优化问题。 这个过程允许我们在一个高维问题上运行受限 BayesOpt，在这个问题中，不明显的是可行性最优解在哪里。 这个练习的解决方案包含在 CH08/03 - Exercise 2.ipynb 笔记本中：
- 重新创建在 CH07/04 - Exercise 2.ipynb 笔记本中使用的 BayesOpt 问题，包括名为
flight_utility()的飞机效用目标函数、我们搜索空间的边界（四维单位超立方体）、GP 实现以及在一些训练数据上训练 GP 的辅助函数fit_gp_model()。 - 实现以下成本函数，模拟制造由四维输入指定的飞机设计的成本：
def flight_cost(X): X = X * 20 - 10 part1 = (X[..., 0] - 1) ** 2 i = X.new(range(2, 5)) part2 = torch.sum(i * (2.0 * X[..., 1:] ** 2 - X[..., :-1]) ** 2, ➥dim=-1) return -(part1 + part2) / 100_000 + 2
- 我们的目标是在遵循成本小于或等于 0 的约束条件的情况下最大化目标函数
flight_utility()。
为此，我们将 BayesOpt 策略每个实验的查询次数设置为 50，并指定每个策略需要运行 10 次重复实验：
num_queries = 50 num_repeats = 10
- 如果找不到可行解，则默认值量化优化进度应设置为-2。
default_value = -2 feasible_incumbents = torch.ones((num_repeats, num_queries)) * default_value
- 在此问题上运行受限 EI 策略以及常规 EI 策略；可视化并比较它们的平均进展（以及误差线）。
我们以与第八章相同的方式实现了受限 EI 策略，其中我们将best_f设置为当前可行的最优解（如果找到了可行解）或默认值-2（如果没有找到可行解）。我们的模型列表包含目标 GP，其索引为 0，和成本 GP，其索引为 1：
if (train_cost <= 0).any(): ❶ best_f = train_utility[train_cost <= 0].max() ❶ else: ❶ best_f = torch.tensor(default_value) ❶ policy = botorch.acquisition.analytic.ConstrainedExpectedImprovement( model=botorch.models.model_list_gp_regression.ModelListGP( ❷ utility_model, cost_model ❷ ), ❷ best_f=best_f, objective_index=0, ❸ constraints={1: [None, 0]} ❹ )
- ❶ 找到当前最优解的适当值
❷ GP 模型列表
❸ 目标模型的索引
❹ 约束模型的索引和下限和上限
我们如下实现常规 EI 策略：
policy = botorch.acquisition.analytic.ExpectedImprovement( model=utility_model, best_f=train_utility.max(), )
- 图 A.16 显示了我们实现的两个先前策略获得的优化结果。我们看到，受限 EI 通过考虑我们施加的成本约束完全支配了常规 EI 策略。
图 A.16 各种贝叶斯优化策略在受限飞机设计优化示例中所取得的进展。与常规 EI 相比，受限变体平均找到更可行的解决方案。
A.8 第九章：用多信度优化平衡效用和成本
本章有两个练习：
- 练习 1 介绍了跨多个实验测量和可视化优化策略平均性能的过程。
- 练习 2 将我们所知的优化策略应用于具有三个可查询函数的二维问题。
A.8.1 练习 1：在多信度优化中可视化平均性能
在本练习中，我们多次运行优化循环，并学习如何取得平均性能以获得更全面的比较：
- 复制 CH09/03 - 测量性能.ipynb 笔记本中的问题设置和多信度优化循环，并添加另一个变量，表示我们要运行的实验次数（默认为 10 次）。
- 为了方便重复实验，在优化循环代码中添加一个外部循环。这应该是一个具有 10 次迭代的
for循环，每次生成一个不同的随机观察值：
num_repeats = 10 ❶ for trial in range(num_repeats): torch.manual_seed(trial) ❷ train_x = bounds[0] + (bounds[1] - bounds[0]) ❷ ➥* torch.rand(1, 1) ❷ train_x = torch.cat( ❷ [train_x, torch.ones_like(train_x) ❷ ➥* fidelities[0]], dim=1 ❷ ) ❷ train_y = evaluate_all_functions(train_x) ❷ current_budget = 0 ❸ while current_budget < budget_limit: ❸ ... ❸
- ❶ 重复实验 10 次
❷ 生成特定于当前迭代的随机初始训练集
❸ 内循环，直到我们耗尽预算 - 将变量
recommendations和spent_budget的每个变量都设为一个列表的列表，其中每个内部列表跟踪单个实验的优化性能。我们将上一步中嵌套循环的代码添加如下所示：
num_repeats = 10 recommendations = [] ❶ spent_budget = [] ❶ for trial in range(num_repeats): torch.manual_seed(trial) train_x = bounds[0] + (bounds[1] - bounds[0]) * torch.rand(1, 1) train_x = torch.cat( [train_x, torch.ones_like(train_x) * fidelities[0]], dim=1 ) train_y = evaluate_all_functions(train_x) current_budget = 0 recommendations.append([]) ❷ spent_budget.append([]) ❷ while current_budget < budget_limit: ... rec_x = get_final_recommendation(model) ❸ recommendations[-1].append ❸ ➥(evaluate_all_functions(rec_x).item()) ❸ spent_budget[-1].append(current_budget) ❸ ...
- ❶ 每个变量都是一个（当前为空的）列表的列表。
❷ 向每个列表的列表添加一个空列表，以供下一个实验使用
❸ 将优化进度统计信息添加到每个变量的最新列表中 - 在我们的优化问题上运行多保真度 MES 策略及其单保真度版本。
- 我们首先制作正则网格和当前空白的插值推荐值，稍后我们将填写其中：
xs = np.arange(budget_limit) interp_incumbents = np.empty((num_repeats, budget_limit))
- 然后我们遍历
recommendations中的每个列表（在我们的代码中重命名为incumbents）和spend_budget，计算线性插值，然后填写interp_incumbents中的值：
for i, (tmp_incumbents, tmp_budget) in enumerate( zip(incumbents, spent_budget) ): interp_incumbents[i, :] = np.interp(xs, tmp_budget, tmp_incumbents)
- 使用线性插值值绘制我们运行的两个策略的平均性能和误差条，并比较它们的性能。比较可视化在图 A.17 中展示，我们可以看到多保真度 MES 策略大大优于其单保真度竞争对手。
图 A.17 显示了在 10 次实验中单一保真度和多保真度 MES 策略在 Forrester 函数上的平均优化进展。多保真度策略大大优于单保真度策略。 - 绘制线性插值曲线，代表各个运行的优化进展，以及平均性能和误差条。比较可视化在图 A.18 中呈现。事实上，我们每次运行的优化进展，由最大后验平均推荐值来衡量，并不是单调递增的。
图 A.18 线性插值曲线代表了在 10 次实验中各个运行的优化进展。我们每次运行的优化进展，由最大后验平均推荐值来衡量，并不是单调递增的。
贝叶斯优化实战（四）（5）
