贝叶斯优化实战（四）（2）https://developer.aliyun.com/article/1516485
A.5 第六章：使用信息论与基于熵的策略
本章中有两个练习：
- 第一个练习涵盖了二分搜索的变体，其中在做出决策时可以考虑先验信息。
- 第二个练习将引导我们通过在前几章中看到的超参数调整问题中实现最大值熵搜索（MES）的过程。
A.5.1 练习 1：将先验知识纳入熵搜索
这个练习，实现在 CH06/02 - Exercise 1.ipynb 中，展示了在找到信息论最优决策时使用不同先验的一个实例，最终将帮助我们进一步欣赏熵搜索作为一种通用决策在不确定性下的过程的优雅和灵活性：
- 证明 Pr(X = 1) + Pr(X = 2) + … + Pr(X = 10) = 1。
我们可以通过简单地将概率相加来实现这一点：
1 / 2 + 1 / 4 + … + 1 / 2⁹ + 1 / 2⁹ = 1 / 2 + 1 / 4 + … + 1 / 2⁸ + 1 / 2⁸ = … = 1 / 2 + 1 / 2 = 1。 - 计算这个先验分布的熵。
请记住熵的公式为 –Σ*[i]*p[i] log p[i]。我们可以编写一个计算此和的 Python 函数：
def compute_entropy(first, last): entropy = 0 for i in range(first, last + 1): p = marginal_probability(i, first, last) ❶ entropy += -p * np.log2(p) ❷ return entropy
- ❶ 获取当前概率。
❷ 对项求和。
此函数将first和last作为参数，它们对应于 X 可能的最小和最大值（起始值为 1 和 10），分别。然后我们遍历first和last之间的数字，并累加 –p[i] log p[i] 项。这里，marginal_probability()是一个计算 Pr(X = n) 的辅助函数，我们实现如下：
def marginal_probability(floor, first, last): if floor == last: ❶ return 2 ** -(last - first) ❶ return 2 ** -(floor - first + 1)
- ❶ 当底层是可能的最高层时的边缘情况
运行compute_entropy(1,10)将给出 1.99609375. 这是 X 先验分布的熵。 - 鉴于在 1 到 10 之间定义的先验分布，从二楼掉落时手机会损坏的概率是多少？从第五层呢？第一层呢？
手机从二楼掉落时损坏的概率恰好是Pr(X = 1)，即 0.5。
手机从第五层掉落损坏的概率是X ≤ 4 的概率，即Pr(X = 1) + Pr(X = 2) + Pr(X = 3) + Pr(X = 4) = 15/16 = 0.9375。
这两个计算可以作为一个函数来实现：
def cumulative_density(floor, first, last): return sum( ❶ [ ❶ marginal_probability(i, first, last) ❶ for i in range(first, floor) ❶ ] ❶ ) ❶
- ❶ 对于小于阈值的X的概率求和
由于我们的先验知识规定，如果从一楼掉落，手机不会损坏，这个概率为 0。 - 计算在我们在第五层进行试验的两种情况下的虚拟后验分布的熵。
使用我们实现的compute_entropy()函数，我们可以计算两种情况下的熵。如果手机损坏，我们将first设为1，last设为4；否则，我们将first设为5，last设为10：
>>> compute_entropy(1, 4) Output: 1.75 >>> compute_entropy(5, 10) Output: 1.9375
- 鉴于先验分布，计算在第五层进行试验后的预期后验熵。
对于这个预期后验熵计算，我们已经进行了必要的计算。首先，手机从第五层掉落损坏的概率是 0.9375，在这种情况下，后验熵是 1.75。其次，手机从第五层掉落不损坏的概率是 1 - 0.9375 = 0.0625，在这种情况下，后验熵是 1.9375。
对两种情况取平均值得到 (0.9375) 1.75 + (0.0625) 1.9375 = 1.76171875。这是你在第五层进行试验后的预期后验熵。 - 计算其他楼层的这个预期后验熵。
我们可以实现一个函数来进行我们刚刚讨论的计算：
def compute_expected_posterior_entropy(floor, first, last): break_probability = cumulative_density ➥(floor, first, last) ❶ return ( ❷ break_probability * compute_entropy ❷ ➥(first, floor - 1) ❷ + (1 - break_probability) * compute_entropy ❷ ➥(floor, last) ❷ ) ❷
- ❶ 从给定楼层手机损坏的概率
❷ 对两种情况取平均值
使用这个函数，我们可以绘制出在 1 到 10 之间数字的预期后验熵。
这个图表显示在图 A.11 中，告诉我们我们第一次试验的信息论最佳地点是二楼，因为 2 给出了最低的预期后验熵（因此，不确定性最低）。
图 A.11 预期后验熵作为试验地点的函数
我们看到这与决策二进制搜索建议的不同，5。这是我们使用先验分布对X进行编码的领域知识的直接影响：因为X = 2 的可能性很高（50%），如果手机损坏，直接尝试这个数字可能会更好，因为这样我们可能会立即找到答案。
有趣的是，从一楼掉下手机不会减少熵。这是因为我们确定手机不会从这一层摔坏，所以在进行这次试验之后，我们对世界的认知不会改变。
A.5.2 练习 2：超参数调优的 BayesOpt
此练习在 CH06/ 03- Exercise 2.ipynb 笔记本中实现，将 BayesOpt 应用于模拟超参数调整任务中支持向量机模型的精度曲面的目标函数：
- 在 CH04/ 03- Exercise 2.ipynb 中重新创建 BayesOpt 循环，包括实现重复实验的外部循环。
- 运行 MES 策略。
考虑到我们的目标函数具有两个维度，我们应将 MES 使用的 Sobol 序列的大小设置为 2,000：
num_candidates = 2000 num_repeats = 10 incumbents = torch.zeros((num_repeats, num_queries)) for trial in range(num_repeats): ... ❶ for i in tqdm(range(num_queries)): ... ❷ sobol = torch.quasirandom.SobolEngine(1, scramble=True) candidate_x = sobol.draw(num_candidates) candidate_x = bounds[0] + (bounds[1] - bounds[0]) * candidate_x policy = botorch.acquisition.max_value_entropy_search.qMaxValueEntropy( model, candidate_x ) ... ❸
- ❶ 随机生成初始训练数据
❷ 记录现任值并重新训练模型
❸ 寻找最大化收获值的点，查询目标函数，并更新训练数据
图 A.12 显示了 MES 的优化效果。我们可以看到，该策略在迄今为止学到的所有 BayesOpt 策略中都具有竞争力。
图 A.12 展示了 MES 的优化效果。该策略在所显示的四个策略中表现最佳。
A.6 第七章：通过批量优化最大化吞吐量
本章包含两个练习：
- 第一个部分涵盖了在批量设置下实现 TS 的方法。
- 第二个部分介绍了如何在一个四维气动结构优化问题上运行 BayesOpt 策略。
A.6.1 练习 1：通过重新采样将 TS 扩展到批量设置
请记住，在我们在第 5.3 节学到的顺序设置中，TS 是从当前 GP 关于目标函数的信念中抽取一个样本，并查询最大化该样本的数据点。在批次设置中，我们只需重复此过程多次，从而多次采样 GP 并最大化样本，以组装所需大小的批次查询。此练习的代码可以在 CH07/ 02- Exercise 1.ipynb 笔记本中找到：
- 重新创建 CH05 / 01- BayesOpt loop.ipynb 笔记本中的批次 BayesOpt 循环。
- 根据第 5.3 节介绍的方法，使用 Sobol 采样器实现 TS。
我们按照以下方式实施该策略，其中我们使用 2,000 元素的 Sobol 序列，并指定样本数作为批次大小：
num_candidates = 2000 ❶ ... ❷ for i in tqdm(range(num_iters)): ... ❸ sobol = torch.quasirandom.SobolEngine ➥(2, scramble=True) ❹ candidate_x = sobol.draw(num_candidates) ❹ candidate_x = (bounds[1] - bounds[0]) * ❹ ➥candidate_x + bounds[0] ❹ ts = botorch.generation.MaxPosteriorSampling(model, ➥replacement=False) next_x = ts(candidate_x, num_samples=batch_size) ❺ ... ❻
- ❶ 指定 Sobol 序列的长度
❷ 随机选择初始训练数据
❸ 重新训练 GP
❹ 初始化 Sobol 序列
❺ 从 GP 中绘制多个样本
❻ 查询目标函数并更新训练数据 - 运行此 TS 策略对超参数调整目标函数进行观察其性能。
运行批次 TS 后，我们可以绘制该策略对我们所学到的其他策略的进展，如图 A.13 所示。在仅进行第一批查询之后，TS 就能取得显著进展。
图 A.13 批处理 TS 在超参数调整示例中的进展。只查询了第一批后，该策略就取得了显著进展。
A.6.2 练习 2: 优化飞机设计
这个练习提供了一个目标函数，用于模拟基准测试飞机设计的过程。代码提供在 CH07/04 - Exercise 2.ipynb 笔记本中。完成以下步骤：
- 实现模拟性能基准测试的目标函数。
目标函数的代码已经提供，所以我们只需将其复制粘贴到程序中：
def flight_utility(X): X_copy = X.detach().clone() X_copy[:, [2, 3]] = 1 - X_copy[:, [2, 3]] X_copy = X_copy * 10 - 5 return -0.005 * (X_copy**4 - 16 * X_copy**2 + 5 * X_copy).sum(dim=-1) + 3
- 使用一个常数均值函数和 Matérn 2.5 核实现一个 GP 模型，输出范围由一个
gpytorch.kernels.ScaleKernel对象实现。
这个 GP 模型的类实现和之前大部分相同，只需要在 ARD 核中指定正确的维度数：
class GPModel( gpytorch.models.ExactGP, botorch.models.gpytorch.GPyTorchModel, botorch.models.model.FantasizeMixin ): _num_outputs = 1 def __init__(self, train_x, train_y, likelihood): super().__init__(train_x, train_y, likelihood) self.mean_module = gpytorch.means.ConstantMean() self.covar_module = gpytorch.kernels.ScaleKernel( gpytorch.kernels.MaternKernel( ❶ nu=2.5, ❶ ard_num_dims=4 ❶ ) ❶ ) def forward(self, x): ...
- ❶ 四维的 Matérn 2.5 核
- 实现一个辅助函数，用于在给定的训练数据集上训练 GP。
我们可以直接从本章其他笔记本（即 02 - Exercise 1.ipynb）中复制相同的辅助函数fit_gp_model()，因为我们不需要在此辅助函数中进行修改。 - 定义优化问题的设置。
我们先定义搜索空间的边界：
lb = 0 ub = 1 bounds = torch.tensor([[lb] * 4, [ub] * 4], dtype=torch.float)
- 然后我们指定可以进行的查询数量、批次大小和每个策略要重复的实验次数：
num_experiments = 5 num_queries = 100 batch_size = 5 num_iters = num_queries // batch_size
- 运行本章学习的每个批处理 BayesOpt 策略在先前实现的目标函数上。
我们首先使用这段代码实现优化循环和重复每个策略实验的外部循环。具体来说，对于每个单独的实验，我们随机在搜索空间内取一个数据点，然后运行每个 BayesOpt 策略，直到我们用完查询次数。下一步我们将看到每个策略是如何定义的：
incumbents = torch.zeros((num_experiments, num_iters)) pbar = tqdm(total=num_experiments * num_iters) for exp in range(num_experiments): torch.manual_seed(exp) ❶ train_x = bounds[0] + (bounds[1] - ❶ ➥bounds[0]) * torch.rand(1, 4) ❶ train_y = flight_utility(train_x) ❶ for i in range(num_iters): incumbents[exp, i] = train_y.max() ❷ model, likelihood = fit_gp_model ➥(train_x, train_y) ❷ ... ❸ next_y = flight_utility(next_x) ❹ train_x = torch.cat([train_x, next_x]) ❹ train_y = torch.cat([train_y, next_y]) ❹ pbar.update()
- ❶ 随机初始化训练数据
❷ 跟踪优化进展并更新预测模型
❸ 定义策略并查找下一个要查询的批次
❹ 查询策略推荐的点并更新训练数据
对于 PoI 策略，我们使用以下代码：
policy = botorch.acquisition.monte_carlo.qProbabilityOfImprovement( model, best_f=train_y.max() )
- 对于 EI 策略，我们使用以下代码：
policy = botorch.acquisition.monte_carlo.qExpectedImprovement( model, best_f=train_y.max() )
- 对于 UCB 策略，我们使用以下代码：
policy = botorch.acquisition.monte_carlo.qUpperConfidenceBound( model, beta=2 )
- 然后可以使用以下代码对这三个策略进行优化：
next_x, acq_val = botorch.optim.optimize_acqf( policy, bounds=bounds, q=batch_size, num_restarts=100, raw_samples=200, )
- 否则，对于 TS 或 MES，我们需要先定义 Sobol 序列：
sobol = torch.quasirandom.SobolEngine(4, scramble=True) ❶ candidate_x = sobol.draw(5000) candidate_x = (bounds[1] - bounds[0]) * candidate_x + bounds[0]
- ❶ 指定维度数为 4
对于 TS 策略，我们使用以下代码：
ts = botorch.generation.MaxPosteriorSampling(model, replacement=False) next_x = ts(candidate_x, num_samples=batch_size)
- 对于 MES，我们使用以下代码来实现循环优化，其中使用了辅助函数
optimize_acqf_cyclic()。请注意，我们指定循环优化的最大迭代次数为 5：
policy = botorch.acquisition.max_value_entropy_search.qMaxValueEntropy( model, candidate_x ) next_x, acq_val = botorch.optim.optimize_acqf_cyclic( policy, bounds=bounds, q=batch_size, num_restarts=40, raw_samples=100, cyclic_options={"maxiter": 5} ❶ )
- ❶ 指定循环优化的最大迭代次数
- 绘制我们运行的 BayesOpt 策略的优化进展并观察其性能。
图 A.14 显示了我们实现的策略得到的优化结果。我们看到大多数策略是可比的，除了 TS；批量 PoI 稍微领先一点。
图 A.14 各种 BayesOpt 策略在飞机设计优化示例中的进展。大多数策略是可比的，除了 TS。
贝叶斯优化实战（四）（4）
