一文理解贝叶斯优化

在线运行和调试文中代码

一个小故事

💭第一天：出师不利

你是石油大亨贝爷，来到一片富饶的土地开采石油。虽然线报告诉你土地下有着丰富的石油，但是人生地不熟的你，并不知道应该在哪里打井，你决定先随意找片“风水宝地”打个井试试。

然而，你的运气不大好，这一钻下去并没有发现油田。“至少我知道这里没有油！”你总是很乐观，安慰自己到，“明天又是新的一天了！”

💭第二天：小试牛刀

吸取昨天的教训，你估摸着中部大概没什么油田。你投了个硬币，打算去西边打个井，看看那儿有没有油田。

运气不错，钻井里涌出黑色的石油原浆！这里确实有个小油田！“至少这一趟采油之旅不会无功而返了！”你舒了口气，“现在我知道西边有个小油田，中部是贫油区了！”

💭第三天：一场豪赌

你昨晚回家算了算账，这个小油田只够抵过这一趟来的成本，并赚不到什么钱，你不是很甘心，打算再打个井碰碰运气。“不过这一次，我应该在西边的小油田附近再打个井看看有没有其他小油田呢？还是到尚未探索的东边碰碰运气看看有没有大油田呢？”你嘀咕到，“不过，既然我已经在西边找到了一个小油田，即使再有油田的话，应该也不会特别大。那我不如去东边碰碰运气，看看有没有大油田呢！”

你下定决心，重整旗鼓，大步迈向神秘的东方。随着钻井打下去，浓厚的石油蓬勃而出！“风浪越大鱼越贵！”这一次你赌对了，兴奋地喊了起来！

💭第四天：粗中有细

因为昨天的成功，你的钻井队都跃跃欲试，想再接再厉、做大做强，挥动着最后一个钻头，问你再去哪里探个大油田！这时你冷静地思考了会儿，觉着昨天确实运气不错，应该是钻到一个大油田了，但是昨天过于兴奋，其实还不是百分百确定这个油田是不是真的很大呢！与其再去那些不确定的地方碰运气，不如把最后这次钻井机会还是放在昨天探到的“大油田”附近，确认一下这个油田是不是确实绵延不绝、货真价实。如果真有预想中的那么大，根据历史经验，这就是这个片区最大的油田了！

钻进队直呼你不仅神机妙算，更是粗中有细！在昨天的大油田附近又打了一眼油井！“有油呦！”随着钻进队传来的喜报，你知道今年“阿不福思名人榜——最成功石油大亨”的名号又非你莫属了！

💭第五天：贝爷’s传奇

你的大成功在石油界引起了轰动，果然阿不福思名人榜的记者闻讯赶来，采访你的成功经验。

你春风得意，娓娓道来：“

• 其实我初来乍到的时候，对这片土地一概不知，这片土地上的任何位置出现石油的可能性都是一样的。
• 第一天在中部打下油井，没有发现油田的时候，我至少知道中部出现油田的概率会比东部和西部更低，因此我接下来就不会继续探索中部了。
• 于是我决定第二天去西部碰碰运气，结果发现了一个小油田。这时候在整片土地上，我知道西部有小油田，中部没有油，只有东部对我来说是未知的了。
• 因为西部的小油田不大，而且因为我已经发现一个小油田，所以附近即使还有油田，也不会特别大。所以继续开发西部对我来说回报不高，我决定不如去探索一下风险与机遇并存的东部。所以第三天我带着钻进队去了神秘的东方，结果运气不错，确实被我赌对了，我们发现了一个大油田！
• 最后一天，我还剩一个钻头，还有最后一次钻井的机会。我觉得在其他地方发现比东部更大的油田概率不高，于是我决定在东部大油田附近的位置再打个井，好好开发一下东部大油田，确认是不是真的足够大。最后事实证明，这应该就是这个片区最大的油田了。

我觉得我最大的优势在于平衡了探索和开发，充分利用有限的钻井次数，既挑战风险探索了不确定性高的区域，又抓住机遇开发了回报高的油田。”

记者崇拜地听完了你的介绍，迅速整理了稿件——《贝爷’s传奇》～

💭第六天：一波拉踩

显然你的成功经验介绍引发了更大范围的关注！“那你怎么评价其他的石油大亨呢？”记者总喜欢搞事情，拉踩最容易博人眼球！

你清了清嗓：“

• 靠过往经验吃饭的任工，以前在炼丹炉里练出过火眼金睛。他在熟悉的地盘上确实能通过经验快速找到大油田，不过一旦到了不熟悉的地方就发挥不出他的学识了。
• 艺高人胆大闭着眼睛指沙盘的隋姐，有时候就运气爆棚能找到大油田，有时候运气就不大好，看天吃饭啦，发挥不稳定。
• 兢兢业业的王哥比较严谨，会把地盘划分成好多个等分，基本上只要划地够密，一个一个眼打下去就肯定能找到大油田啦。但是现在打井成本越来越高了，他们公司让他少打点井，成功率就下去了。
• 还有那个有着魔法的的提督，天天举着法典算来算去，在我们这种两眼一抹黑的问题上行不通的啦。

”

你双手在胸前交叉，抬头挺胸、气宇轩昂地总结道：“对于这种每次评估成本都很高，而且没法获得什么额外信息的这种黑盒问题，我就是最牛的！”

正经的贝叶斯优化介绍：一种不基于梯度的全局优化算法

以上的小故事其实蕴含了一种不基于梯度的全局优化算法——“贝叶斯优化”的核心思想，而贝叶斯优化被广泛应用于复杂黑盒函数的优化问题上。

接下来，就让我们一起来了解一下贝叶斯优化吧！

贝叶斯优化（Bayesian Optimization）是一种高效的全局优化方法，主要用于优化具有高度非线性、高维度和未知性质的目标函数。在机器学习和深度学习领域，贝叶斯优化经常用于调整模型的超参数，以提高模型的性能。贝叶斯优化的优势在于它可以通过利用已有的观测数据，来指导对下一个采样点的选择，从而加速优化过程。
——ChatGPT

1. 黑盒优化问题（Black Box Optimization）

通俗的讲，贝叶斯优化适用于最大化一个评估起来很昂贵的黑盒函数$f:\mathcal{X}\to \mathbb{R}$类型的问题，即求解

$$x^{\ast }=\arg \underset{x\in \mathcal{X}}{\max }f(x)$$

我们无法获得函数$f$的函数形式，唯一的办法是通过昂贵的代价获得$f$在一系列测试点上的评估值，并期望在少量评估次数后，确定接近最佳值的点$x^{\ast }$。

# 如下，我们可以定义一个非常非常复杂的黑盒函数
from typing import Dict
import numpy as np

def black_box_function(x):
    return - x * np.sin(x)
# 我们可以一起看看这个黑盒函数长什么样

import matplotlib.pyplot as plt

# 让我们看看black_box_function长什么样
x = np.linspace(0, 10, 100)
y = black_box_function(x)

plt.plot(x,y)

font = {'family': 'serif',
        'color':  'black',
        'weight': 'normal',
        'size': 15}

plt.xlabel("$x$", fontdict=font)
plt.ylabel("$y$", fontdict=font)
plt.xlim(0,10)

plt.show()

<Figure size 640x480 with 1 Axes>

2. 代理模型（surrogate model）

为了在少量评估中优化$f$，我们需要一种方法来推断我们在尚未评估的点$x$处$f$可能的表现。这个时候，我们可以使用一个模型，在一些已经评估过的点上，去拟合$f$，以获得对其他未评估的点的预测能力。

在贝叶斯优化中，这个模型一般被称为代理模型（surrogate model）。特别的，贝叶斯优化中要求代理模型能通过后验分布（posterior distribution），去量化其在点$x$处$f(x)$预测值的不确定度（uncertainty）。待会儿会讲到为什么需要去量化预测值的不确定度，不要着急！

为了实现这一目标，贝叶斯优化中最常用的代理模型一般是高斯过程（Gaussian Process)模型。高斯过程模型的好处是它在任何一个有限集合中的后延都可以由一系列多元正态分布给出。

# 这里，我们使用高斯过程模型来作为代理模型
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
kernel = Matern()
gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, normalize_y=True, n_restarts_optimizer=5)

n_initial_samples = 10
X_sample = np.random.uniform(0, 10, size=(n_initial_samples, 1))
Y_sample = black_box_function(X_sample)

gp.fit(X_sample, Y_sample)
# 我们可以看到高斯过程模型除了能给出预测值（MEAN）外，还能给出不确定度的估计（VAR）

GaussianProcessRegressor(alpha=1e-06, kernel=Matern(length_scale=1, nu=1.5),
                         n_restarts_optimizer=5, normalize_y=True)

import pandas as pd
X_test = np.random.uniform(0, 10, size=(10, 1))
Y_test = black_box_function(X_test)

mean, var = gp.predict(X_test, return_std=True)

df = pd.DataFrame({"LABEL": Y_test.reshape(-1), "MEAN": mean.reshape(-1), "VAR": var.reshape(-1)})
df

	LABEL	MEAN	VAR
0	-1.474424	-1.475075	0.264146
1	2.393722	-2.818221	1.859164
2	-0.552635	-3.443158	1.357517
3	4.798706	4.060867	1.170749
4	3.827768	-2.583765	2.060149
5	3.564018	3.582482	0.147019
6	2.824129	1.853092	1.773595
7	-1.805934	-1.809172	0.155350
8	-0.966256	-0.970078	0.222544
9	3.927909	2.827580	1.740191

3. 采集函数（Acquisition Functions）

这里就是贝叶斯优化的核心了！

对于一般的机器学习模型拟合优化问题而言，我们的优化步骤一般是这样的：

1. 在给定的数据集上，训练一个机器学习模型；
2. 将训练好的机器学习模型，在要预测的数据上进行推理，获得预测值。

如果这时你有额外的预算来标记更多数据，那你还会有以下几步：

1. 在给定的数据集上，训练一个机器学习模型；
2. 将训练好的机器学习模型，在要预测的数据上进行推理，获得预测值；
3. 额外的挑选一些样本，进行数据标注；
4. 将新标注的数据添加到训练集中，重复123步，直到模型收敛或预算用尽。

这里将引入一个额外的问题，也就是**“如何选择额外的样本点，来帮助模型预测获得更大的提升”**。

• 一般来说，对模型预测不确定度高的点进行观测，能有效提升模型的表现——这可以通过高斯过程模型去获得每个预测点的不确定度，我们标记为$\sigma$；
• 额外的，我们需要注意我们的优化目标是找到接近最佳值的点$x^{\ast }$，所以我们还需要额外考虑当前模型在每一个点的预测值的大小**，这也可以通过高斯过程模型的均值求得，我们标记为$\mu$；

因此，在这个找最佳值的问题中，我们既希望去观察均值大的点来实现我们找到接近最佳值点的目标（贝叶斯优化中把这个过程叫做exploritation——开发），我们同时又需要去挑选那些方差大的点进行额外的观测，来整体提升模型表现，避免漏掉发现全局最优解的机会（贝叶斯优化中把这个过程叫做exploration——探索）。那么究竟什么时候开发什么时候探索，并且开发和探索各占多少比例呢？

这就是**采集函数（Acquisition Functions）**要发挥作用的地方了！

采集函数是一个权衡探索和开发的函数，一般以预测值的均值$\mu$和方差$\sigma$作为入参，去获得每个预测点是否值得额外评估的综合得分，以最为经典的GP-UCB（Gaussian Process-Upper Confidence Bound）为例，它的函数形式为：

$$x=\arg \underset{x\in \mathbb{R}}{\max }\alpha (x)=\arg \underset{x\in \mathbb{R}}{\max }(\mu (x)+{\beta }^{1/2}\sigma (x{)}^2)$$

很容易发现，$\alpha (x)$其实就是均值$\mu$和方差$\sigma$的加权和，起到了权衡探索和开发的作用。同时，我们可以通过对超参数$\beta$的调整，来控制当前的采样应该更注重探索还是开发。

## GP-UCB的定义
def gaussian_process_upper_confidence_bound(x, beta=1):
    x = x.reshape(-1,1)
    mean, sigma = gp.predict(x, return_std=True)
    return mean.reshape(-1) + beta**1/2 * sigma.reshape(-1)**2
## 我们可以一起看看刚才的预测值经过采集函数得到的打分

df["GP-UCB"] = gaussian_process_upper_confidence_bound(X_test)
df

	LABEL	MEAN	VAR	GP-UCB
0	-1.474424	-1.475075	0.264146	-1.440188
1	2.393722	-2.818221	1.859164	-1.089975
2	-0.552635	-3.443158	1.357517	-2.521731
3	4.798706	4.060867	1.170749	4.746194
4	3.827768	-2.583765	2.060149	-0.461658
5	3.564018	3.582482	0.147019	3.593289
6	2.824129	1.853092	1.773595	3.425912
7	-1.805934	-1.809172	0.155350	-1.797105
8	-0.966256	-0.970078	0.222544	-0.945315
9	3.927909	2.827580	1.740191	4.341712

贝叶斯优化的完整流程

综合上面的介绍，我们可以得到一个贝叶斯优化的完整流程如下：

1. 我们首先选择一个代理模型来对真实函数$f$建模并定义其先验；
2. 给定一组观察值（函数评估），请使用贝叶斯规则获取后验；
3. 使用采集函数$\alpha (x)$，基于后验的均值和方差来确定下一个采样点；
4. 将新采样的数据添加到观测值集中，然后执行步骤2，直到收敛或用尽预算。

也就是说贝叶斯优化实际上是：

• 由于目标函数无法直接优化，那我们找一个代理函数来优化；
• 为了找到当前目标函数的合适的代理函数，赋予了代理函数概率分布；
• 然后根据当前已有的先验知识，使用采集函数来评估如何选择新的采样点；
• 获得新的采样点之后更新先验知识，然后继续迭代找到最合适的替代函数，最终求得全局最优解。

import matplotlib.pyplot as plt 

def plot_gp(gp, X_sample, Y_sample, iteration):
    x = np.linspace(0, 10, 100).reshape(-1, 1)
    y = black_box_function(x)
    
    mu, std = gp.predict(x, return_std=True)
    upper_confidence_bound = mu.reshape(-1) + 2 * std.reshape(-1)
    lower_confidence_bound = mu.reshape(-1) - 2 * std.reshape(-1)
    gpucb = gaussian_process_upper_confidence_bound(x)
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
    
    ax1.plot(x, y, label="Objective Function")
    ax1.scatter(X_sample, Y_sample, color='red', label="Sampled Points")
    ax1.plot(x, mu, label="GP Mean")
    ax1.fill_between(
        x.flatten(), 
        lower_confidence_bound.flatten(), 
        upper_confidence_bound.flatten(), 
        alpha=0.3, label="2 Std Dev")
    ax1.scatter(X_sample[-1], Y_sample[-1], color='green', label="New Sampled Point")
    ax1.legend()
    ax1.set_title(f"Iteration {iteration}", fontdict=font)
    ax1.set_xlabel("$x$", fontdict=font)
    ax1.set_ylabel("$y$", fontdict=font)
    ax1.set_xlim(0,10)
    
    ax2.plot(x, gpucb, label="Acquisition Function (UCB)")
    ax2.set_title(f"Iteration {iteration}")
    ax2.set_title(f"Iteration {iteration}", fontdict=font)
    ax2.set_xlabel("$x$", fontdict=font)
    ax2.set_ylabel("$y$", fontdict=font)
    ax2.set_xlim(0,10)
    ax2.legend()
    
    plt.show()

## 接下来我们就可以定义这个优化流程了！
from scipy.optimize import minimize

def bayesian_optimization(n_iterations=20, n_init_samples=5):
    X_sample = np.random.uniform(0, 10, size=(n_init_samples, 1))
    Y_sample = black_box_function(X_sample)
    for i in range(n_iterations):
        gp.fit(X_sample, Y_sample)
        plot_gp(gp, X_sample, Y_sample, i+1)
        next_x = minimize(gaussian_process_upper_confidence_bound, x0=np.random.uniform(0, 10), bounds=[(0, 10)]).x
        next_y = black_box_function(next_x)
        X_sample = np.vstack([X_sample, next_x])
        Y_sample = np.vstack([Y_sample, next_y])
    return X_sample, Y_sample

## 让我们来看看贝叶斯优化的表现！
X_opt, Y_opt = bayesian_optimization(10, 2)
min_index = np.argmin(Y_opt)
min_x = X_opt[min_index]
min_y = Y_opt[min_index]
print("Minimum value found at x =", min_x, "with y =", min_y)

<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>



<Figure size 1200x300 with 2 Axes>

🔖以上就是贝叶斯优化一个简单的介绍

当你遇到的场景有以下特点时，你可以试试贝叶斯优化：

• 你面向的优化目标的评估成本特别高；
• 评估函数$f$没有解析表达式，无法获取梯度等额外信息；
• 目标是找到函数$f$在自变量$x$上的全局最优解。

例如，神经网络的超参数（模型大小、学习率等）搜索就是一个典型的贝叶斯优化的适用场景！快来试试吧！