使用贝叶斯优化工具实践XGBoost回归模型调参

0. 关于调参

0.1. 超参数

在机器学习的上下文中，超参数（hyper parameters）是在开始学习过程之前设置值的参数，而不是通过训练得到的参数数据。通常情况下，需要对超参数进行优化，给学习机选择一组最优超参数，以提高学习的性能和效果。

例如深度学习中的学习率、批次大小（batch_size）、优化器（optimizer）等，XGBoost算法中的最大的树深度（max_depth）、子样本比例（subsample）、子节点最小权重和（min_child_weight）等。

0.2. 调参方法

超参数设置通常调参方法如下：

• 随机搜索 （RandomSearch），搜索超参数的值不是固定，是在一定范围内随机的值；
• 网格搜索（GridSearch ），给定超参和取值范围，遍历所有组合得到最优参数。首先你要给定一个先验的取值，不能取得太多，否则组合太多，耗时太长，可以启发式的尝试；
• 贝叶斯优化（Bayesian optimization），采用高斯过程迭代式的寻找最优参数，每次迭代都是在上一次迭代基础上拟合高斯函数上，寻找比上一次迭代更优的参数。

本文主要是分享贝叶斯优化调参方法，其他略。

1. 贝叶斯优化理论

贝叶斯优化是一种逼近思想，当计算非常复杂、迭代次数较高时能起到很好的效果，多用于超参数确定。

贝叶斯优化算法主要包含两个核心部分——概率代理模型（probabilistic surrogate model）和采集函数
(acquisition function)。

• 概率代理模型包含先验概率模型和观测模型：先验概率模$p(f)$；观测模型描述观测数据生成的机
  制，即似然分布$p(D|f)$。更新概率代理模型意味着根据公式 (1)得到包含更多数据信息的后验概率分布
  $p(f|D_i)$。
• 采集函数是根据后验概率分布（高斯过程）构造的，通过最大化采集函数来选择下一个最有 “潜 力”的评估点。同时 ，有效的采集函数能够保证选择的评估点序列使得总损失(loss)最小．损失有 时表示为regret：
  $$r_n=|y^{\ast }-y_n|$$
  或者累计表示为：$$R_i=\sum _{i=1}^n$$其中，$y^{\ast }$表示当前最优解。

掌握贝叶斯优化调参，须要从三个部分讲起：

• 贝叶斯定理
• 高斯过程，用以拟合优化目标函数
• 贝叶斯优化，包括了“开采”和“勘探”，用以花最少的代价找到最优值

1.1. 贝叶斯定理

贝叶斯优化名称的由来是因为参数优化中使用了贝叶斯定理。

$$p(f|D_i)=\frac{p(D_i|f)p(f)}{p(D_i)}(1)$$

其中，$f$表示未知目标函数（或者标识参数模型中的参数）：

• $D_i=(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$表示已观测集合，$x_i$表示决策向量，$y_i=f(x_i)+{\epsilon }_i$表示观测值误差；
• $p(D_i|f)$表示$y$的似然分布；
• $p(f)$表示$f$先验概率分布，对未知目标函数的假设；
• $p(D_i)$表示边际化$f$的边际似然分布或者“证据”，由于该边际似然存在概率密度函数的乘积和积 分，通 常难以得到明确的解析式，该边际似然在贝叶斯优化中主要用于优化超参数(hyper-parameter)；
• $p(f|D_i)$表示$f$的后验概率分布，后验概率分布描述通过己观测数据集对先验进行修正后未知目标函数的置信度。

贝叶斯定理是关于随机事件A和B的条件概率（或边缘概率）的一则定理。其中P(A|B)是在B发生的情况下A发生的可能性。

1.2. 高斯过程（Gaussian Process）

高斯过程（GP）全称是高斯分布随机过程，是一个可以被用来表示函数的分布情况的模型。当前，机器学习的常见做法是把函数参数化，然后用产生的参数建模来规避分布表示（如线性回归的权重）。但GP不同，它直接对函数建模生成非参数模型。由此产生的一个突出优势就是它不仅能模拟任何黑盒函数，还能模拟不确定性。这种对不确定性的量化是十分重要的，如当我们被允许请求更多数据时，依靠高斯过程，我们能探索最不可能实现高效训练的数据区域。这也是贝叶斯优化背后的主要思想。

那么高斯分布是如何计算的呢？

假设$f\sim GP(\mu ,K)$，对于$GP$高斯过程，其中$\mu$为均值，$K$为协方差核。所以预测也是服从正态分布的，即有$p(y|x,D)=N(\hat{y},{\hat{\sigma }}^2)$。

当随机变量的维度上升到有限$n$的时候，称之为高维高斯分布，$p(x)=N(\mu ,{\sum }_{n\times n})$。

对于一个$n$维的高斯分布而言，决定它的分布是两个参数：

• $n$维均值向量${\mu }_n$，表示$n$维高斯分布中各个维度随机变量的期望；
• $n\times n$的协方差矩阵${\sum }_{n\times n}$，表示高维分布中，各个维度自身的方差，以及不同维度间的协方差。

核函数(协方差函数)

核函数是一个高斯过程的核心，核函数决定了一个高斯过程的性质。核函数在高斯过程中起生成一个协方差矩阵(相关系数矩阵)来衡量任意两个点之间的“距离”。不同的核函数有不同的衡量方法，得到的高斯过程的性质也不一样。最常用的一个核函数为高斯核函数，也成为径向基函数 RBF。其基本形式如下。

$$K(x_i,x_j)={\sigma }^{2}exp(-\frac{||x_i-x_j|{|}_2^2}{2l^2})$$

其中$\sigma$和$l$是高斯核的超参数。

本文使用的BayesianOptimization的源代码中，高斯过程使用了sklearn.gaussian_process，如下截部分取代码所示。

from sklearn.gaussian_process.kernels import Matern
from sklearn.gaussian_process import GaussianProcessRegressor
......
        self._random_state = ensure_rng(random_state)

        # Data structure containing the function to be optimized, the bounds of
        # its domain, and a record of the evaluations we have done so far
        self._space = TargetSpace(f, pbounds, random_state)

        # queue
        self._queue = Queue()

        # Internal GP regressor
        self._gp = GaussianProcessRegressor(
            kernel=Matern(nu=2.5),
            alpha=1e-6,
            normalize_y=True,
            n_restarts_optimizer=5,
            random_state=self._random_state,
        )

        self._verbose = verbose
        self._bounds_transformer = bounds_transformer
        if self._bounds_transformer:
            self._bounds_transformer.initialize(self._space)

        super(BayesianOptimization, self).__init__(events=DEFAULT_EVENTS)

上图是一张高斯分布拟合函数的示意图，可以看到，它只需要九个点，就可以大致拟合出整个函数形状（图片来自：https://github.com/fmfn/BayesianOptimization）

1.3. 贝叶斯优化

参数优化基本思想是基于数据使用贝叶斯定理估计目标函数的后验分布，然后再根据分布选择下一个采样的超参数组合。它充分利用了前一个采样点的信息，其优化的工作方式是通过对目标函数形状的学习，并找到使结果向全局最大提升的参数。

假设一组超参数组合是$X=x_1,x_2,...x_n$($x_n$表示某一个超参数的值)，不同超参数会得到不同效果，贝叶斯优化假设超参数与最后模型需要优化的损失函数存在一个函数关系。

高斯过程用贝叶斯优化中对目标函数建模，得到其后验分布。

而目前机器学习其实是一个可见输入与输出的黑盒子(black box)，所以很难确直接定存在什么样的函数关系，所以我们需要将注意力转移到一个我们可以解决的函数上去，下面开始正式介绍贝叶斯优化。

贝叶斯优化的大体思路如下：

假设我们有一个函数$f:x\to \mathbb{R}$,我们需要在$x\subseteq X$内找到：
$$x^{\ast }=\underset{x\in X}{\mathrm{argmin}}f(x)(2)$$

注意上面的$x$表示的是超参数。以XGBoost分类任务为例，$x$可以是eta（学习率），max_depth（最大的树深度）等超参数的设置。而为了避免全文符号太多，所以将输入数据隐去了。

当$f$是凸函数且定义域$X$也是凸的时候，我们可以通过已被广泛研究的凸优化来处理，但是$f$并不一定是凸的，而且在机器学习中$f$通常是expensive black-box function，即计算一次需要花费大量资源。那么贝叶斯优化是如何处理这一问题的呢？

贝叶斯优化把搜索的模型空间假设为高斯分布，利用高斯过程，按迭代的方式每次计算得到比当前最优参数期望提升的新的最优参数。

例如详细算法，Sequential model-based optimization (SMBO) 是贝叶斯优化的最简形式，其算法思路如下：

Input（输入）：

• $f$是模型（所谓黑盒子），即输入一组超参数，得到一个输出值；
• $X$是超参数搜索空间等；
• $D$是表示一个由若干对数据组成的数据集，每一对数组表示为$(x,y)$，$x$是一组超参数,$y$表示该组超参数对应的结果（如精确率等）；
• $S$是参数选择算法Acquisition Function；
• $M$是对数据集$D$进行拟合得到的模型，可以用来假设的模型有随机森林、Tree Parzen Estimators、高斯模型等。本文实践是高斯拟合函数。

$InitSamples(f,x)\to D$

这一步骤就是初始化获取数据集$D=(x_1,y_1),...,(x_n,y_n)$，其中$y_i=f(x_i)$，这些都是已知的。

循环迭代参数$T$次

因为每次选出参数$x$后都需要计算$f(x)$，而正如前面介绍的每计算一次函数$f$，都会消耗大量资源，所以一般需要固定选参次数(或者是函数评估次数)。

$p(y|x,D)\leftarrow FITMODEL(M,D)$

首先，预先假设了模型$M$服从高斯分布，且已知了数据集$D$，所以可以通过计算得出具体的模型具体函数表示。

通过高斯过程建模之后，我们尝试抽样进行样本计算，而贝叶斯优化很容易在局部最优解上不断采样，这就涉及到了开发和探索之间的权衡。

• 开发 (exploitation)： 根据后验分布，在最可能出现全局最优解的区域进行采样, 开发高意味着均值高
• 探索 (exploration): 在还未取样的区域获取采样点， 探索高意味着方差高
  而如何高效的采样，即开发和探索，我们需要用到 Acquisition Function, 它是用来寻找下一个 x 的函数。

一般形式的Acquisition Funtion是关于x的函数，映射到实数空间R，表示改点的目标函数值能够比当前最优值大多少的概率，目前主要有以下几种主流的Acquisition Function：

• POI
• Expected Improvement
• Confidence bound criteria

2. XBoost回归算法使用贝叶斯优化调参

2.1. 高斯过程全局优化的Python实现

目前可以做贝叶斯优化的第三方python包很多，常用的有BayesianOptimization、bayesopt、skopt等等。本文使用BayesianOptimization为例，利用XGBoost的回归模型进行趋势预测。

安装第三方包：bayesian-optimization

• pip install bayesian-optimization

包的依赖：

• Numpy
• Scipy
• Scikit-learn

2.2. XGBoost参数调优实践案例

我们需要开始的就是实例化一个BayesianOptimization对象，指定要优化的函数$f$，代码中的 _xgb_evaluate以及它的参数及其相应的边界。这是一种受约束的优化技术，因此必须为每个参数指定可以探测的最小值和最大值，才能使其工作，例如XGBoost训练参数的边界参数：

 {'max_depth': (5, 11), 
   'gamma': (0, 1),
   'colsample_bytree': (0.5, 0.9),
   'min_child_weight': (5, 300)}

核心示例代码如下：

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
import xgboost as xgb
from xgboost import plot_importance
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score
from bayes_opt import BayesianOptimization

class ChurnDayModel(object):
    def __init__(self, df, model_filename):
        self.df = df
        self.model_filename = model_filename # 定义持久化模型文件路径及文件名称
           
    def evaluate(self):
        dtrain = xgb.DMatrix(self.x_train, label=self.y_train)
        
        def _xgb_evaluate(max_depth, gamma, colsample_bytree, min_child_weight):
            params = {'eval_metric': 'rmse',
                      'max_depth': int(max_depth),
                      'subsample': 0.8,
                      'eta': 0.3,
                      'gamma': gamma,
                      'colsample_bytree': colsample_bytree,
                      'min_child_weight': min_child_weight}
            # Used around 1000 boosting rounds in the full model
            cv_result = xgb.cv(params, dtrain, num_boost_round=100, nfold=3)    # nfold 分三组交叉验证

            # Bayesian optimization only knows how to maximize, not minimize, so return the negative RMSE
            return -1.0 * cv_result['test-rmse-mean'].iloc[-1]
                
        xgb_bo = BayesianOptimization(_xgb_evaluate, {'max_depth': (5, 11), 
                                             'gamma': (0, 1),
                                             'colsample_bytree': (0.5, 0.9),
                                             'min_child_weight': (5, 300)})
        # Use the expected improvement acquisition function to handle negative numbers
        # Optimally needs quite a few more initiation points and number of iterations
        xgb_bo.maximize(init_points=4, n_iter=25, acq='ei')
        
        print(xgb_bo.max)
        res = xgb_bo.max
        params_max = res['params']

        return params_max

    def set_train_data(self,flag_col):
        X = self.df.drop([flag_col],axis=1)
        Y = self.df[[flag_col]]
        #修改列名为中文
        X = X.rename(columns=chinesNames())
        # 归一化数据
        X = self._feature_StandardScaler(X, False) # True
        
        self.x_train,self.x_test, self.y_train, self.y_test = train_test_split(X,Y,test_size=0.3)

Model_Filename = 'XGB_Model/XGboostChurnDays.model'
df = pd.read_csv('card_churndays.csv')
# 数据特征
Columns_Name = ['carduser_id', 'balance', 'accsumamount', 'tradecount',......]

cols = Columns_Name.copy()
cols.remove('carduser_id')
cols.remove('churncount')

df0 = df[cols]
CM = ChurnDayModel(df0,Model_Filename)
CM.set_train_data('churndays')
params = CM.evaluate()

重点代码解析：

   maximize( init_points=5,
             n_iter=25,
             acq='ucb',

参数：

• init_points：初始点，int，可选（默认值=5），探索开始探索之前的迭代次数
• n_iter：最大迭代次数，int，可选（默认值=25），方法试图找到最大值的迭代次数
• acq：acq:{‘ucb’，‘ei’，‘poi’}，使用的采集方法。
  – “ucb”代表置信上限法
  –“ei”是预期的改进方法
  – *“poi”是改进的概率标准。

2.3. 优化输出最大值

输出最优结果符合模型的特点，例如最大的树深度，将对模型准确率有利，所以最优结果，其最大的树深度将与设置最大的树深度边界值基本一致。

3. 总结

对于缺乏经验的模型开发者，贝叶斯优化将为开发者提供一系列可用的较为优化的超参，指导开发者及早训练出符合要求的模型。对于有经验的模型开发者，贝叶斯优化将为开发者提供方便、自动化的调优，降低开发者的工作量，及早训练出符合要求的模型。

贝叶斯优化参数方法，给出的最优值，实际上未必是最优值。同时，比较耗费资源，需要比较长的时间进行模型训练。贝叶斯优化基于训练集数据，仅是使训练集上的模型最优，容易产生过拟合的情况，需要开发者注意降低过拟合参数边界设置。

在实现中，为了获得更好的训练效果，我们往往要做更多调整计算。你也许已经注意到了，GP包含两个非常重要的参数σ和l，如果你在之前采集样本的时候尝试改变过它们，那你会发现图像在垂直和水平方向上的神奇变化。例如，如果我们期望更大范围的输出，我们就需要相应地放大参数σ。事实上，和所有会用到核函数的方法一样，如果有需要，我们甚至可以完全改变核函数。

参考：
[1]. fmfn. https://github.com/fmfn/BayesianOptimization
[2]. BRIAN YUHAS. Bayesian Optimization with XGBoost . Kaggle. 2018.08
[3]. 杨睿. 强大而精致的机器学习调参方法：贝叶斯优化 . 博客园. 2018.07
[4]. wedo实验君. Python 机器学习：超参数调优. 微信公众号. Python中文社区, 2021.01
[5]. marsggbo. 贝叶斯优化(Bayesian Optimization)深入理解. 知乎. 2020.05
[6]. 你好世界炼丹师. 贝叶斯优化(Bayesian Optimization)只需要看这一篇就够了,算法到python实现, 知乎, 2020.07
[7]. SIGAI. 理解贝叶斯优化, 知乎. 2020.06
[8]. vector_xfy. 深度学习超参数介绍及调参. CSDN博客, 2019.07
[9]. 崔佳旭,杨 博. 贝叶斯优化方法和应用综述. 软件学报, 2018,29(10)
[10]. Alex Bridgland. 图文详解高斯过程（一）——含代码, 知乎. 公众号. 论智, 2017.12
[11]. 石溪. 如何通俗易懂地介绍 Gaussian Process？, 知乎, 2021.06
[12]. 王桂波. 高斯过程 Gaussian Processes 原理、可视化及代码实现. 知乎, 2021.08