贝叶斯优化（Bayesian Optimization）在Python中的实现：高斯过程与采集函数
大家好，今天我们来聊聊贝叶斯优化，这是一种非常强大的全局优化方法，尤其适用于目标函数计算代价昂贵，且没有显式表达式的情况。我们将深入探讨贝叶斯优化的核心组成部分：高斯过程和采集函数，并通过Python代码演示如何实现它们。
1. 贝叶斯优化简介
想象一下，你要调整一个机器学习模型的超参数，例如学习率、正则化系数等。每次评估一组超参数的性能都需要训练模型并在验证集上进行测试，这个过程可能非常耗时。传统的网格搜索或随机搜索效率较低，因为它们没有利用之前的评估结果来指导下一步搜索。
贝叶斯优化正是为了解决这类问题而生的。它通过构建目标函数的概率模型（通常是高斯过程），并利用采集函数来决定下一个要评估的点，从而在尽可能少的迭代次数内找到全局最优解。
核心思想：
- 代理模型（Surrogate Model）： 使用一个易于计算的概率模型（如高斯过程）来近似目标函数。
- 采集函数（Acquisition Function）： 根据代理模型，选择下一个最有希望改进目标函数的点。采集函数平衡了探索（exploration）和利用（exploitation）。
2. 高斯过程（Gaussian Process）
高斯过程是贝叶斯优化的核心。它是一个随机过程，其有限个随机变量的联合分布都服从多元高斯分布。换句话说，对于任何一组输入点，高斯过程都会给出一个关于这些点对应函数值的联合高斯分布。
2.1 高斯过程的数学定义：
一个高斯过程由均值函数 m(x) 和协方差函数 k(x, x') 唯一确定，记作：
f(x) ~ GP(m(x), k(x, x'))
其中：
m(x)是输入x处的均值，通常设为0。k(x, x')是输入x和x'之间的协方差，描述了函数在不同点之间的相关性。协方差函数也称为核函数。
2.2 常用的核函数：
-
RBF (Radial Basis Function) 核 (也称为 Gaussian 核):
k(x, x') = σ^2 * exp(-||x - x'||^2 / (2 * l^2))其中：
σ^2是信号方差，控制函数值的幅度。l是长度尺度，控制函数的光滑程度。较小的l意味着函数变化更快。
-
Matérn 核:
Matérn 核是一类更一般的核函数，RBF 核是 Matérn 核在
ν -> ∞时的极限情况。Matérn 核的形式如下：k(x, x') = σ^2 * (2^(1-ν) / Γ(ν)) * (√(2ν) * ||x - x'|| / l)^ν * K_ν(√(2ν) * ||x - x'|| / l)其中：
ν是一个正数，控制函数的光滑程度。ν = 5/2和ν = 3/2是常用的特殊情况。Γ(ν)是 Gamma 函数。K_ν是第二类修正贝塞尔函数。
-
线性核:
k(x, x') = x^T * x'
2.3 高斯过程的预测：
假设我们已经观测到一些数据点 X = [x_1, x_2, ..., x_n] 及其对应的函数值 y = [f(x_1), f(x_2), ..., f(x_n)]。现在，我们想预测新点 x_* 的函数值 f(x_*)。根据高斯过程的定义，[y, f(x_*)] 服从联合高斯分布：
[y, f(x_*)] ~ N([m(X), m(x_*)] , [[K(X, X), K(X, x_*)], [K(x_*, X), K(x_*, x_*)]])
其中：
K(X, X)是X中所有点对之间的协方差矩阵。K(X, x_*)是X中的点与x_*之间的协方差向量。K(x_*, x_*)是x_*与自身之间的协方差。
根据高斯过程的性质，f(x_*) 的后验分布也是高斯分布：
f(x_*) | X, y, x_* ~ N(μ_*, σ_*^2)
其中：
μ_* = m(x_*) + K(x_*, X) * K(X, X)^-1 * (y - m(X))σ_*^2 = K(x_*, x_*) - K(x_*, X) * K(X, X)^-1 * K(X, x_*)
μ_* 是 f(x_*) 的后验均值，可以作为预测值。σ_*^2 是 f(x_*) 的后验方差，表示预测的不确定性。
2.4 Python 代码实现高斯过程：
import numpy as np
from scipy.linalg import solve
from scipy.stats import norm
class GaussianProcess:
def __init__(self, kernel, noise_level=1e-8):
self.kernel = kernel
self.noise_level = noise_level
self.X = None
self.y = None
self.K = None
self.L = None # Cholesky decomposition
def fit(self, X, y):
"""
根据观测数据拟合高斯过程。
"""
self.X = X
self.y = y
self.K = self.kernel(X, X) + np.eye(len(X)) * self.noise_level
self.L = np.linalg.cholesky(self.K) # Cholesky 分解
def predict(self, X_star, return_std=False):
"""
预测新点的函数值和标准差。
"""
K_star = self.kernel(self.X, X_star)
# 计算后验均值
alpha = solve(self.L.T, solve(self.L, self.y)) #solve(A,b) 返回方程 Ax=b 的解
mu_star = K_star.T @ alpha
# 计算后验方差
v = solve(self.L, K_star)
var_star = self.kernel(X_star, X_star) - v.T @ v
std_star = np.sqrt(np.diag(var_star))
if return_std:
return mu_star, std_star
else:
return mu_star
def rbf_kernel(X1, X2, l=1.0, sigma_f=1.0):
"""
RBF (Gaussian) 核函数。
"""
sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)
# 示例用法
if __name__ == '__main__':
# 训练数据
X = np.array([-4, -3, -2, -1, 1, 2]).reshape(-1, 1)
y = np.sin(X)
# 测试数据
X_star = np.linspace(-5, 5, 100).reshape(-1, 1)
# 创建高斯过程对象
gp = GaussianProcess(kernel=rbf_kernel)
# 拟合高斯过程
gp.fit(X, y)
# 预测
mu_star, std_star = gp.predict(X_star, return_std=True)
# 打印预测结果
print("预测均值:", mu_star[:5])
print("预测标准差:", std_star[:5])
# 可视化（需要安装 matplotlib）
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.plot(X, y, 'rx', markersize=10, label='Observations')
plt.plot(X_star, mu_star, 'b-', label='Mean Prediction')
plt.fill_between(X_star.flatten(), mu_star - 1.96 * std_star, mu_star + 1.96 * std_star, alpha=0.2, color='b', label='95% Confidence Interval')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Gaussian Process Regression')
plt.legend()
plt.grid(True)
plt.show()这个代码首先定义了一个 GaussianProcess 类，包含了 fit 和 predict 方法。fit 方法根据观测数据计算协方差矩阵并进行 Cholesky 分解，predict 方法利用后验分布公式计算预测均值和标准差。 同时定义了RBF核函数，作为高斯过程的协方差函数。最后，代码演示了如何使用 GaussianProcess 类进行回归预测，并用 matplotlib 可视化结果。
3. 采集函数（Acquisition Function）
采集函数是贝叶斯优化的另一个关键组成部分。它利用高斯过程的后验分布，来衡量在某个点进行评估的“价值”。采集函数平衡了探索（exploration）和利用（exploitation）：
- 利用（Exploitation）： 选择预测均值高的点，希望找到更好的解。
- 探索（Exploration）： 选择预测方差大的点，希望了解更多关于目标函数的信息。
3.1 常用的采集函数：
-
概率改进 (Probability of Improvement, PI):
PI(x) = P(f(x) > f(x^+))其中，
f(x^+)是当前找到的最佳函数值。PI 衡量的是在x处评估函数值超过当前最佳值的概率。 -
期望改进 (Expected Improvement, EI):
EI(x) = E[max(0, f(x) - f(x^+))]EI 衡量的是在
x处评估函数值超过当前最佳值的期望值。EI 比 PI 更关注改进的幅度。 -
置信上限 (Upper Confidence Bound, UCB):
UCB(x) = μ(x) + κ * σ(x)其中，
μ(x)是预测均值，σ(x)是预测标准差，κ是一个调节参数，控制探索和利用的平衡。较大的κ鼓励更多的探索。
3.2 Python 代码实现采集函数：
from scipy.stats import norm
def expected_improvement(X, gp, xi=0.01):
"""
计算期望改进 (Expected Improvement)。
"""
mu, sigma = gp.predict(X, return_std=True)
mu_sample_opt = np.max(gp.y)
with np.errstate(divide='warn'):
imp = mu - mu_sample_opt - xi
Z = imp / sigma
ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
ei[sigma == 0.0] = 0.0
return ei
def upper_confidence_bound(X, gp, kappa=1.96):
"""
计算置信上限 (Upper Confidence Bound)。
"""
mu, sigma = gp.predict(X, return_std=True)
return mu + kappa * sigma
# 示例用法
if __name__ == '__main__':
# 创建高斯过程对象 (使用之前的代码)
X = np.array([-4, -3, -2, -1, 1, 2]).reshape(-1, 1)
y = np.sin(X)
gp = GaussianProcess(kernel=rbf_kernel)
gp.fit(X, y)
# 测试数据
X_star = np.linspace(-5, 5, 100).reshape(-1, 1)
# 计算期望改进
ei = expected_improvement(X_star, gp)
# 计算置信上限
ucb = upper_confidence_bound(X_star, gp)
# 打印结果
print("期望改进:", ei[:5])
print("置信上限:", ucb[:5])
# 可视化（需要安装 matplotlib）
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.plot(X_star, ei, 'g-', label='Expected Improvement')
plt.plot(X_star, ucb, 'm-', label='Upper Confidence Bound')
plt.xlabel('x')
plt.ylabel('Acquisition Function Value')
plt.title('Acquisition Functions')
plt.legend()
plt.grid(True)
plt.show()这段代码实现了 expected_improvement 和 upper_confidence_bound 两个采集函数。代码中，首先计算高斯过程的预测均值和标准差，然后根据公式计算采集函数值。最后，代码演示了如何使用这两个采集函数，并用 matplotlib 可视化结果。
4. 贝叶斯优化的完整流程
现在，我们已经了解了高斯过程和采集函数，可以将它们组合起来，实现一个完整的贝叶斯优化流程：
- 初始化： 随机选择一些点进行评估，得到初始数据
X和y。 - 拟合高斯过程： 使用
X和y拟合高斯过程模型。 - 选择下一个评估点： 使用采集函数选择下一个最有希望改进目标函数的点
x_*。通常可以使用优化算法（如 L-BFGS-B）来最大化采集函数。 - 评估目标函数： 在
x_*处评估目标函数，得到函数值f(x_*)。 - 更新数据： 将
x_*和f(x_*)添加到X和y中。 - 重复步骤 2-5： 重复上述步骤，直到达到最大迭代次数或满足其他停止条件。
4.1 Python 代码实现贝叶斯优化：
from scipy.optimize import minimize
def bayesian_optimization(func, bounds, n_iter=10, gp_params=None, acq_func='ei'):
"""
贝叶斯优化主函数。
Args:
func: 目标函数，接受一个 numpy 数组作为输入，返回一个标量。
bounds: 搜索空间的边界，一个列表，每个元素是一个包含最小值和最大值的元组。
n_iter: 迭代次数。
gp_params: 高斯过程的参数，一个字典。
acq_func: 采集函数，'ei' (Expected Improvement) 或 'ucb' (Upper Confidence Bound)。
Returns:
(X, y): 找到的最佳输入和对应的函数值。
"""
# 默认高斯过程参数
if gp_params is None:
gp_params = {'kernel': rbf_kernel, 'noise_level': 1e-8}
# 初始化
X = np.array([[np.random.uniform(b[0], b[1]) for b in bounds] for _ in range(5)]) # 随机初始化5个点
y = np.array([func(x) for x in X])
gp = GaussianProcess(**gp_params)
gp.fit(X, y)
# 优化循环
for i in range(n_iter):
# 定义采集函数
if acq_func == 'ei':
def acquisition(x):
return -expected_improvement(x.reshape(1, -1), gp) # 负号是因为要最大化采集函数
elif acq_func == 'ucb':
def acquisition(x):
return -upper_confidence_bound(x.reshape(1, -1), gp)
else:
raise ValueError("Invalid acquisition function: {}".format(acq_func))
# 使用优化算法最大化采集函数
x0 = np.array([np.random.uniform(b[0], b[1]) for b in bounds]) # 随机初始化
result = minimize(acquisition, x0=x0, bounds=bounds, method='L-BFGS-B') # 约束优化算法
x_next = result.x
# 评估目标函数
y_next = func(x_next)
# 更新数据
X = np.vstack((X, x_next))
y = np.append(y, y_next)
# 重新拟合高斯过程
gp.fit(X, y)
print(f"Iteration {i+1}: x={x_next}, y={y_next}")
# 返回结果
best_index = np.argmax(y)
return X[best_index], y[best_index]
# 示例目标函数
def objective_function(x):
"""
一个简单的目标函数。
"""
return - (x[0]**2 + (x[1] - 1)**2) # 负号是因为贝叶斯优化是寻找最大值
# 示例用法
if __name__ == '__main__':
# 搜索空间边界
bounds = [(-5, 5), (-5, 5)]
# 运行贝叶斯优化
best_x, best_y = bayesian_optimization(objective_function, bounds, n_iter=20, acq_func='ei')
# 打印结果
print("Best x:", best_x)
print("Best y:", best_y)这段代码实现了一个通用的 bayesian_optimization 函数。它接受目标函数、搜索空间边界、迭代次数、高斯过程参数和采集函数作为输入。函数使用 L-BFGS-B 算法来最大化采集函数，并返回找到的最佳输入和对应的函数值。示例代码定义了一个简单的二维目标函数，并演示了如何使用 bayesian_optimization 函数来寻找其最大值。
5. 贝叶斯优化的优点与缺点
优点：
- 效率高： 特别适用于目标函数计算代价昂贵的情况。
- 全局优化： 能够找到全局最优解，而不仅仅是局部最优解。
- 无需梯度信息： 不需要目标函数的梯度信息，适用于黑盒优化。
- 自适应性： 能够根据之前的评估结果自适应地调整搜索策略。
缺点：
- 计算复杂度高： 高斯过程的计算复杂度随着数据点的增加而增加，可能成为瓶颈。
- 参数敏感： 高斯过程的核函数和采集函数的参数对优化结果有很大影响。
- 维度限制： 在高维空间中，高斯过程的性能可能会下降。
6. 总结与思考
今天，我们深入探讨了贝叶斯优化的原理和实现，包括高斯过程和采集函数。通过 Python 代码，我们了解了如何构建高斯过程模型，如何选择合适的采集函数，以及如何实现一个完整的贝叶斯优化流程。
贝叶斯优化适用于目标函数计算代价昂贵且无显式表达式的场景，通过构建代理模型和采集函数，在探索和利用之间取得平衡，最终找到全局最优解。
掌握贝叶斯优化的核心概念和实现细节，能够帮助我们更好地解决实际问题，例如超参数优化、实验设计等。
更多IT精英技术系列讲座，到智猿学院