优化算法学习——贝叶斯优化(Bayesian Optimization)

暮雨٩(๑˃̵ᴗ˂̵๑)۶终将落下 发布于 2024-11-04 1819 次阅读

概述

贝叶斯优化算法是机器学习领域一种经典的优化算法，本文从一个实际的问题出发$$ y = x^2 + z^2 $$ 使用python代码， 一步步实现BO算法。

BO算法基本流程

一个经典的Bayesian Optimization流程如下。

graph LR A[初始化] --> B[选择代理模型] B --> C[定义目标函数] C --> D[选择初始样本点] D --> E[代理模型拟合] E --> F[选择采集函数] F --> G[优化采集函数] G --> H[评估新的样本点] H --> I{终止条件} I -->|是| J[输出最优解] I -->|否| E

高斯过程

我们知道，高斯过程是用来替代昂贵的目标函数的，也就是流程图里的代理模型。比如说，我们的机器性能很差，计算$$ y = x^2 + z^2 $$要花费1个小时（当然，这里只是做一个假设，因为实际工程上的目标函数十分复杂，计算一次要非常久），所以我们就要使用一个其它的万能函数来代替目标函数，这个过程就是高斯过程回归。

class GaussianProcess:
    def __init__(self, kernel_variance=1.0, noise_variance=1e-10):
        self.kernel_variance = kernel_variance
        self.noise_variance = noise_variance
        self.X_train = None
        self.y_train = None
        self.alpha = None
        self.L = None

    def kernel(self, X1, X2):
        sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
        return self.kernel_variance * np.exp(-0.5 * sqdist)

    def fit(self, X_train, y_train):
        self.X_train = X_train
        self.y_train = y_train
        K = self.kernel(X_train, X_train) + self.noise_variance * np.eye(len(X_train))
        self.L = np.linalg.cholesky(K)
        self.alpha = np.linalg.solve(self.L.T, np.linalg.solve(self.L, y_train))

    def predict(self, X_s):
        K_s = self.kernel(self.X_train, X_s)
        K_ss = self.kernel(X_s, X_s) + 1e-8 * np.eye(len(X_s))
        v = np.linalg.solve(self.L, K_s)
        mu_s = K_s.T.dot(self.alpha)
        cov_s = K_ss - v.T.dot(v)
        return mu_s, np.diag(cov_s)

我们看上面这个高斯过程，init函数用来初始化所有的变量，高斯过程涉及到三个比较重要的函数

核函数

核函数（也称为协方差函数）用于衡量输入数据点之间的相似性，这里我们使用RBF核函数$$k(x_i, x_j) = \sigma^2 \exp\left(-\frac{|x_i - x_j|^2}{2\ell^2}\right)$$

拟合函数 fit

这一步骤也可以叫做回归过程，就像最小二乘法求一次函数系数也叫做回归一样，我们要使用高斯函数去预测未知点的均值$\mu_s $ 和$ \Sigma_s$，我们也必须要求出高斯过程中的一些“回归系数”，这里我们需要按照下面的式子计算未知点的值

\begin{equation}
\label{eq:mu}
\mu_s = \mathbf{K}_s^T \mathbf{K}^{-1} \mathbf{y}_{\text{train}}
\end{equation}
\begin{equation}
\label{eq:sigma}
\mathbf{\Sigma}_s = \mathbf{K}_{ss} - \mathbf{K}_s^T \mathbf{K}^{-1} \mathbf{K}_s
\end{equation}

对于上述公式中的 $ \mathbf{K}^{-1}\mathbf{y}_{\text{train}} $ ,我们使用 $ \alpha $ 表示, $ \alpha $可以从方程 $ \eqref{eq:alpha} $ 解出。

\begin{equation}
\label{eq:alpha}
\mathbf{K\alpha} = \mathbf{y}_{train}
\end{equation}

对应到代码里，我们可以看到除了 $ \mathbf{\alpha} $ ，我们还将 $ \mathbf{K} $ 分解成为了一个下三角矩阵 $ \mathbf{L} $ ,之所以这样做是因为矩阵 $\mathbf{K}$ 是拟合过程中所有点的协方差，衡量输入点的相似度，该矩阵中的某些数值可能会导致直接求逆出现数值不稳定，所以我们将其分解为下三角矩阵。同时由于该矩阵是正定对称矩阵，我们在代码中使用Cholesky分解可以更快得到结果。

    def fit(self, X_train, y_train):
        self.X_train = X_train
        self.y_train = y_train
        K = self.kernel(X_train, X_train) + self.noise_variance * np.eye(len(X_train))
        self.L = np.linalg.cholesky(K)
        self.alpha = np.linalg.solve(self.L.T, np.linalg.solve(self.L, y_train))

预测函数 predict

预测函数predict的作用是当输入一个不知道y值的x，我们根据高斯过程中已知的量和方程 $ \eqref{eq:mu} \eqref{eq:sigma} $ 去求出对应的预测值。这里的 $ \mu_s $ 和 $ \Sigma_s $ 分别代表预测出来的平均值和方差

Bayesian Optimization

有了这个高斯过程之后我们来理解BO过程，贝叶斯优化算法，说白了也就是生成一大堆数据点，然后找到数据点里值最小的点，但是一般我们目标函数计算开销很大，所以我们使用高斯模型计算这些点。
这里有一个问题，我们的高斯模型计算出来的值是平均值和方差，所以要有一个acquisition_function利用平均值和方差计算预测的值，当我们求最小值时，我们可以直接用平均值减去方差。
高斯过程给出的预测的值并不是实际值，我们拿到acquisition_function给出的预测值后，找到最小的点，然后求出该点的实际值，把该点加入GP过程的fit函数，对代理模型进行更新，重复这个过程直到达到最大迭代次数，或者值的精度达到了要求。

本文用的python代码

import numpy as np

def objective_function(x, z):
    return x**2 + z**2

class GaussianProcess:
    def __init__(self, kernel_variance=1.0, noise_variance=1e-10):
        self.kernel_variance = kernel_variance
        self.noise_variance = noise_variance
        self.X_train = None
        self.y_train = None
        self.alpha = None
        self.L = None

    def kernel(self, X1, X2):
        sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
        return self.kernel_variance * np.exp(-0.5 * sqdist)

    def fit(self, X_train, y_train):
        self.X_train = X_train
        self.y_train = y_train
        K = self.kernel(X_train, X_train) + self.noise_variance * np.eye(len(X_train))
        self.L = np.linalg.cholesky(K)
        self.alpha = np.linalg.solve(self.L.T, np.linalg.solve(self.L, y_train))

    def predict(self, X_s):
        K_s = self.kernel(self.X_train, X_s)
        K_ss = self.kernel(X_s, X_s) + 1e-8 * np.eye(len(X_s))
        v = np.linalg.solve(self.L, K_s)
        mu_s = K_s.T.dot(self.alpha)
        cov_s = K_ss - v.T.dot(v)
        return mu_s, np.diag(cov_s)

def acquisition_function(mu, sigma, kappa=2.0):
    return mu - kappa * sigma

def bayesian_optimization(n_iter=500, bounds=(-10, 10)):
    # 初始化
    X_train = np.random.uniform(bounds[0], bounds[1], (5, 2))
    y_train = np.array([objective_function(x[0], x[1]) for x in X_train])

    gp = GaussianProcess()
    gp.fit(X_train, y_train)

    for _ in range(n_iter):
        # 预测
        X_sample = np.random.uniform(bounds[0], bounds[1], (100, 2))
        mu, sigma = gp.predict(X_sample)
        # 计算采集函数
        acquisition_values = acquisition_function(mu, sigma)
        # 选择下一个评估点
        next_sample = X_sample[np.argmin(acquisition_values)]
        # 评估目标函数
        y_next = objective_function(next_sample[0], next_sample[1])
        # 更新数据集
        X_train = np.vstack((X_train, next_sample))
        y_train = np.append(y_train, y_next)
        # 重新拟合高斯过程
        gp.fit(X_train, y_train)
    return X_train[np.argmin(y_train)], np.min(y_train)

# 运行贝叶斯优化
best_point, best_value = bayesian_optimization()
print(f"Optimal value of x and z: {best_point}")
print(f"Minimum value of objective function: {best_value}")

程序输出结果如下：

Optimal value of x and z: [ 0.00952015 -0.0112046 ]
Minimum value of objective function: 0.0002161764355009679