xgboost原理介绍与代码实现

1.引言

自从2015年以来，各大算法比赛网站上频繁出现了一个神器——xgboost，xgboost是陈天奇提出来的一种梯度提升树的方法，据作者统计，在2015年kaggle上的29个比赛中，有17个冠军就是用了xgboost模型。因此，本文将重点介绍一下这个神器的原理。

• 论文地址：XGBoost: A Scalable Tree Boosting System

2. xgboost原理

2.1 正则化目标函数

给定一个数据集，假设样本量为$n$，特征数为$m$， D = { ( x i , y i ) } ( ∣ D ∣ = n , x i ∈ R m , y i ∈ R ) \mathcal{D}=\left\{\left(\mathbf{x}_{i}, y_{i}\right)\right\}\left(|\mathcal{D}|=n, \mathbf{x}_{i} \in \mathbb{R}^{m}, y_{i} \in \mathbb{R}\right) D={(xi​,yi​)}(∣D∣=n,xi​∈Rm,yi​∈R)，对于每一个样本，xgboost使用$K$棵树的预测值加和作为模型最后的预测值：
$${\hat{y}}_i=\varphi \left({\mathbf{x}}_i\right)=\sum _{k=1}^K{f}_k\left({\mathbf{x}}_i\right),f_k\in \mathcal{F}$$
其中，$\mathcal{F}=\left\{f(\mathbf{x})=w_{q(\mathbf{x})}\right\}\left(q:{\mathbb{R}}^m\to T,w\in {\mathbb{R}}^T\right)$，$q$表示树的结构或者决策规则，即将一个样本映射到对应的叶节点序号，说白了就是判断一个样本属于哪个叶节点。$T$表示叶节点的数量，$w$表示一个$T$维向量，即整棵决策树所有叶节点对应的值，对于第$i$个叶节点的预测值，可以表示为$w_i$。

xgboost使用正则化的目标函数如下：
$$\mathcal{L}(\varphi )=\sum _{i}l\left({\hat{y}}_i,y_i\right)+\sum _k\Omega \left(f_k\right)$$
其中，$\Omega (f)=\gamma T+\frac{1}{2}\lambda \Vert w{\Vert }^2$，这里$l$表示一个可微凸损失函数，用于计算预测值与真实值直接的差异，$\Omega$表示模型复杂度的惩罚，这里对每棵树的惩罚项主要包含两部分，一部分是树的结点树，另一部分是树每个叶节点预测值的L2值，这里主要是希望每棵树可以相对均匀地预测目标值，防止出现过拟合的情况。当惩罚项的参数为0时，则目标函数就变为传统梯度提升树的目标函数。

2.2 梯度提升树

如果直接对上面的目标函数进行最小化求解，很难对参数进行估计。由于xgboost是对每棵树的预测结果进行加总，因此，对于第$t$次迭代时第$i$个样本的预测值，可以记为${\hat{y}}_i^{(t)}$，此时，对于第$i$棵树的参数求解，只需要最小化如下的目标函数：
$${\mathcal{L}}^{(t)}=\sum _{i=1}^{n}l\left(y_i,{\hat{y}}_i^{(t-1)}+f_t\left({\mathbf{x}}_i\right)\right)+\Omega \left(f_t\right)$$
其中，${\hat{y}}_i^{(t-1)}$表示前面$t-1$棵树的总预测值，对于第$t$棵树来说，此时相当于一个常数，为了使得对不同的$l$函数都可以方便地进行求解，作者采用二次泰勒展开对目标函数进行了近似：
$${\mathcal{L}}^{(t)}\simeq \sum _{i=1}^n\left[l\left(y_i,{\hat{y}}_i^{(t-1)}\right)+g_i{f}_t\left({\mathbf{x}}_i\right)+\frac{1}{2}{h}_i{f}_t^2\left({\mathbf{x}}_i\right)\right]+\Omega \left(f_t\right)$$
其中，$g_i={\partial }_{\hat{y}}(t-1)l\left(y_i,{\hat{y}}^{(t-1)}\right)$，$h_i={\partial }_{{\hat{y}}^{(t-1)}}^{2}l\left(y_i,{\hat{y}}^{(t-1)}\right)$分别表示损失函数函数的一阶和二阶偏导数。去掉常数项后，目标函数可以简化为：
$${\tilde{\mathcal{L}}}^{(t)}=\sum _{i=1}^n\left[g_i{f}_t\left({\mathbf{x}}_i\right)+\frac{1}{2}{h}_i{f}_t^2\left({\mathbf{x}}_i\right)\right]+\Omega \left(f_t\right)$$
记 I j = { i ∣ q ( x i ) = j } I_{j}=\left\{i | q\left(\mathbf{x}_{i}\right)=j\right\} Ij​={i∣q(xi​)=j}为第$t$棵树第$j$个叶节点的样本索引集合，则目标函数可以进一步表示为：
$$\begin{array}{rl}{\tilde{\mathcal{L}}}^{(t)}& =\sum _{i=1}^n\left[g_i{f}_t\left({\mathbf{x}}_i\right)+\frac{1}{2}{h}_i{f}_t^2\left({\mathbf{x}}_i\right)\right]+\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2\\ & =\sum _{j=1}^T\left[\left(\sum _{i\in I_j}{g}_i\right)w_j+\frac{1}{2}\left(\sum _{i\in I_j}{h}_i+\lambda \right)w_j^2\right]+\gamma T\end{array}$$
当决策树的结构已经确定时，也就是决策树的分支规则已经确定时，此时，目标函数相当于$w_j$的二次函数，因此，对其求偏导易得最优解为：
$$w_j^{\ast }=-\frac{{\sum }_{i\in I_j}{g}_i}{{\sum }_{i\in I_j}{h}_i+\lambda }$$
此时，目标函数达到最优值：
$${\tilde{\mathcal{L}}}^{(t)}(q)=-\frac{1}{2}\sum _{j=1}^T\frac{{\left({\sum }_{i\in I_j}{g}_i\right)}^2}{{\sum }_{i\in I_j}{h}_i+\lambda }+\gamma T$$
上面这个目标函数的最优值表达式也作为节点分裂的评价指标，通常来说，要考虑所有可能的分支规则是不可能的，因此，一般会提供贪婪算法，即从一个节点开始，每次选择一个最优的变量作为分裂规则，将数据划分为左右两部分，然后递归下去，直到达到结束的条件。假设一个节点的数据被划分为左右两部分$I_L$、$I_R$，则根据上面目标最优值的计算公式，可以得到分裂后，目标函数的下降的值为：
$${\mathcal{L}}_{split}=\frac{1}{2}\left[\frac{{\left({\sum }_{i\in I_L}{g}_i\right)}^2}{{\sum }_{i\in I_L}{h}_i+\lambda }+\frac{{\left({\sum }_{i\in I_R}{g}_i\right)}^2}{{\sum }_{i\in I_R}{h}_i+\lambda }-\frac{{\left({\sum }_{i\in I}{g}_i\right)}^2}{{\sum }_{i\in I}{h}_i+\lambda }\right]-\gamma$$
通过该函数，可以在每次分支时选择可以使目标函数下降最大的变量和分裂节点作为最优的变量和分裂节点。

在训练过程中，xgboost还集成了另外两种技巧，一种是收缩（Shrinkage），即对每棵树的学习到的分数采用一个类似学习率的参数(\eta)进行放缩；另一种是变量抽样，借鉴于随机森林的技巧，即在每棵树学习时，对变量进行抽样，然后采用抽样后的变量作为候选分支变量，在一定程度上防止模型过拟合。

2.3 分裂查找算法

2.3.1 精确贪婪算法

精确贪婪算法其实跟CART的方法类似，就是穷举每个变量的每种可能取值作为分裂规则，将样本划分为左右两部分，然后计算分裂后目标函数的下降值，下降最大的分裂节点就作为此时的分裂节点。

2.3.2 近似算法

贪婪算法由于需要穷举每个变量的所有可能分裂节点，因此，当数据量比较大，可能会导致内存溢出，因此，作者又提出了两种近似算法，算法的大致思想是先对每个变量按照分位数进行分箱，然后计算每个分箱的$G$和$H$，在节点分支时再直接计算目标函数的下降值。

根据分箱的情况，作者划分了两种模式，即全局模式和局部模式，全局模式是在一棵树构建之前，就进行分箱，然后在每个节点分支时用统一一种分箱规则；局部模式是在每个节点分支时，重新对变量进行分箱。对于全局模式，减少了分箱的计算次数，运行速度快，但是一般需要分箱更细才能达到比较好的效果，对于局部模式，则一般只需要较少的分箱就可以达到比较好的效果，比较适合更深的树。

2.4 分位数的计算（Weighted Quantile Sketch）

上面的近似算法中，我们提到xgboost是采用分位数的形式确定候选分支节点，那么，分位数具体是怎么计算呢？作者采用的是上面目标函数中的二阶偏导值作为计算的指标，即对于每个变量，首先计算如下的配对数据集：
D k = { ( x 1 k , h 1 ) , ( x 2 k , h 2 ) ⋯ ( x n k , h n ) } \mathcal{D}_{k}=\left\{\left(x_{1 k}, h_{1}\right),\left(x_{2 k}, h_{2}\right) \cdots\left(x_{n k}, h_{n}\right)\right\} Dk​={(x1k​,h1​),(x2k​,h2​)⋯(xnk​,hn​)}
然后，定义如下的概率计算公式：
$$r_k(z)=\frac{1}{{\sum }_{(x,h)\in {\mathcal{D}}_k}h}\sum _{(x,h)\in {\mathcal{D}}_k,x<z}$$
其中，$z\in [0,+\infty )$，$r_k(z)$表示$x<z$的样本中，其二阶偏导$h$的和所占所有样本的比重。

那么，xgboost近似算法的目的是获得$l$个候选分裂节点 { s k 1 , s k 2 , ⋯ s k l } \left\{s_{k 1}, s_{k 2}, \cdots s_{k l}\right\} {sk1​,sk2​,⋯skl​}，那么，只需要：
$$\left|r_k\left(s_{k,j}\right)-r_k\left(s_{k,j+1}\right)\right|<ϵ,s_{k1}=\underset{i}{\min }{\mathbf{x}}_{ik},s_{kl}=\underset{i}{\max }{\mathbf{x}}_{ik}$$
其中，$ϵ=1/l$，这样，通过控制$ϵ$即可确定分裂的节点。

2.5 稀疏感知分裂算法（Sparsity-aware Split Finding）

在实际的建模过程中，我们往往还会遇到数据集稀疏的情况，对于稀疏数据集，会导致计算的时间变得更慢，一般导致数据集出现稀疏的原因主要有：①缺失值，②离散型变量转one-hot，③统计值太多0值。为了使得xgboost也适应与稀疏数据的场景，作者在每棵树的每个节点增加了一个默认的方向，如下图所示。当一个数据出现缺失值时，则在分支时，会将缺失值的样本统一分配到默认方向的节点。但是，这里有个问题，在每次分支时，是有左右两个分支的，那么，对于缺失值，应该是统一分到左分支还是右分支呢？xgboost是这样操作的，对于所有样本，首先将缺失值全部移到右边，然后计算分支后每个分裂节点的信息增益，然后接着又遍历一次，将缺失值全部移到左边，重新计算一次信息增益，最后比较哪种情景下信息增益更大则将缺失值归到那一边，其算法如下：

3. xgboost的代码实现

以上就是xgboost原理的大致介绍，作者为了提高模型的训练速度，其实还引入了缓存机制和核外计算等技巧，并开发了xgboost的包，支持python、R、c++、spark等语言，所以xgboost的实现直接使用作者提供的包即可。

• xgboost安装包

以下是笔者在训练过程整理的代码，加入了sklearn的网格搜索来确定xgboost的最优参数，感兴趣的朋友可以看下。

import pickle
import xgboost as xgb
from sklearn.model_selection import GridSearchCV


class XgbRegressor(object):
    def __init__(self, params):
        self.params = params

    def grid_search_cv(self, train_x=None, train_y=None, cv=3, n_jobs=8):
        """采用网格搜索获取最优参数

        :param train_x: 训练集特征变量
        :param train_y: 训练集目标变量
        :param cv:
        :return:
        """
        # 采用网格搜索获取最优参数
        regre = xgb.XGBRegressor()
        gsearchCV = GridSearchCV(regre, param_grid=self.params, scoring='r2', cv=cv, verbose=1, n_jobs=n_jobs)
        gsearchCV.fit(train_x, train_y)

        # 打印最优值和最优参数值
        print("best_score:")
        print("\t", gsearchCV.best_score_)
        print("\t", gsearchCV.best_params_)

        # 将最优参数值进行保存
        self.params = gsearchCV.best_params_

    def train(self, train_x=None, train_y=None, valid_x=None, valid_y=None,
              early_stopping_rounds=10, eval_metric='rmse',
              model_path='./saves/xgboost.model'):
        """根据最优参数训练模型

        :param train_x: 训练集特征变量
        :param train_y: 训练集目标变量
        :param valid_x: 验证集特征变量
        :param valid_y: 验证集目标变量
        :param early_stopping_rounds: 提前结束的迭代次数
        :param eval_metric: 评估指标
        :param model_path: 模型的保存路径
        :return:
        """
        # 初始化模型
        regre = xgb.XGBRegressor()
        for param in self.params.keys():
            setattr(regre, param, self.params[param])

        # 训练模型
        regre.fit(train_x, train_y, eval_metric=eval_metric, verbose=True, eval_set=[(valid_x, valid_y)],
                  early_stopping_rounds=early_stopping_rounds)

        # 保存模型
        pickle.dump(regre, open(model_path, "wb"))
        # regre.get_booster().save_model(model_path)
        return regre

    def load_model(self, save_path):
        """加载模型

        :param save_path: 模型的保存路径
        :return:
        """
        regre = pickle.load(open(save_path, "rb"))
        # regre = xgb.Booster()
        # regre.load_model(save_path)
        return regre

    def predict(self, bst, data):
        """预测

        :param bst:
        :param data:
        :return:
        """
        pred = bst.predict(data, ntree_limit=bst.best_ntree_limit)
        return pred

    def plot_importance(self, bst):
        xgb.plot_importance(bst)

    def plot_tree(self, bst):
        xgb.plot_tree(bst, num_trees=2)

4. 总结

最后，总结一下xgboost模型的优缺点：

• xgboost采用泰勒展开去近似目标函数的求解，在计算上速度更快，并且加入了正则项，使得模型不易于过拟合
• xgboost对大数据集、稀疏数据集的场景都可以很好地兼容，不需要太多地数据预处理过程
• xgboost由于是串行的方式，因此在计算上没法完全并行，模型的超参也比较多，在参数选优上会耗费比较长的时间。