集成算法Boosting之 ---- XGBoost

1. XGBoost简介

---

XGBoost: eXtreme Gradient Boosting
创造者：陈天奇博士 http://homes.cs.washington.edu/~tqchen/
原论文地址：https://arxiv.org/pdf/1603.02754.pdf
项目地址：https://github.com/dmlc/xgboost
PPT地址：https://homes.cs.washington.edu/~tqchen/data/pdf/BoostedTree.pdf
官方文档：https://xgboost.readthedocs.io/en/latest/index.html

XGBoost是基于GradientBoosting算法的一个优化版本，既可用于分类也可以用于回归问题中。何为GradientBoosting? 就是将弱学习器组合成强学习器，由于新的弱学习器的加入，可以纠正前面已有弱学习器的残差，最终形成的模型准确率会比单独的学习器要高(boosting的含义)，并在添加新的弱学习器时使用了梯度下降法来最小化损失函数(gradient的含义)。

XGBoost的实力已被很多机器学习和数据挖掘的比赛所广泛认可，例如在Kaggle经常可以见到选手用它作为最终模型。可以说，很少有别的集成方法效果能超过XGBoost。

2. 算法原理推导

---

2.1 构造目标函数

给定一个有n个样例和m个特征的数据集$D=\{(x_i,y_i)\}(|D|=n,x_i\in R^m,y_i\in R)$，集成方法使用K个弱学习器相加来预测目标变量：

$$\hat{y_i}=\phi(x_i)=\sum_{k=1}^Kf_k(x_i), f_k\in F$$
其中 $F=\{f(x)=w_{q(x)}\}(q:R^m\rightarrow T, w\in R^T)$是CART树的集合，每一个弱分类器 $f_k$有着独立的树结构 $q$和叶子节点的权值$w$ (我们用 $w _i$来代表第 $i$个叶子节点的权值)。因此最后分类某个样本的依据是将它在每个弱分类器中被预测的各个结果的权值汇总在一起，来进行最后的判断。
为了训练得到弱分类器，定义了一个具有正则项的损失函数，我们通过最小化这个函数来得到分类器：
$${\cal L}(\phi)=\sum_il(\hat{y_i}, y_i)+\sum_k\Omega(f_k)$$
其中 $\Omega(f)=\gamma T+\frac{1}{2}\lambda||w||^2$。
上式中， $l$衡量了预测目标值$\hat{y_i}$和真实目标值 $y_i$之间的误差， $\Omega$对整个模型的复杂度进行惩罚，以防止过拟合。正则项 $\Omega$里的T是叶子节点的个数， $\gamma$是这一项的系数， $\lambda$是所有叶子节点的权值的l2正则之和的系数。因此，最小化该损失函数倾向于去选择一个比较简单且预测准确度高的模型。(注：当正则项为0时，就变成了一个传统的GBDT算法)

对数据集以及目标函数和模型架构有了初步的了解以后，我们来进一步深入了解其中的原理。
为了得到最终的强学习器，整个模型的训练是以additive的迭代方式进行的，即每次训练一个弱学习器，再把它们整合在一起。
整个过程可以用描述为：

$$初始化：\hat{y_i}^{(0)}=0\\ 第一轮：\hat{y_i}^{(1)}=f_1(x_i)=\hat{y_i}^{(0)}+f_1(x_i)\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 第二轮：\hat{y_i}^{(2)}=f_1(x_i)+f_2(x_i)=\hat{y_i}^{(1)}+f_2(x_i)\\ \quad\quad\ \ \ \ \ ......\\ \quad\ \ \ \ 第t轮：\hat{y_i}^{(t)}=\sum_{k=1}^tf_k(x_i)=\hat{y_y}^{(t-1)}+f_t(x_i)$$
对于第t轮迭代过程，需要确定的就是 $f_t(x_i)$。假设 $\hat{y_i}^{(t)}$是第 $i$个样本在第t轮的目标预测值，我们需要新加入一个弱学习器$f_t$来最小化以下函数：

$${\cal L}^{(t)}=\sum_{i=1}^nl(y_i, \hat{y_i}^{(t)})+\sum_{i=1}^{t}\Omega(f_i)\\ = \sum_{i=1}^nl(y_i, \hat{y_i}^{(t-1)}+f_t(x_i))+\Omega(f_t)+constant$$
通过优化第t轮迭代的目标函数来确定该轮的 $f_t(x_i)$，目标函数的变量是 $\hat{y_i}^{(t-1)}$。
求解上式的最小值仍旧比较复杂，这里我们考虑使用泰勒二次展开进行简化。

注：泰勒公式是一个用函数在某点的信息描述其附近取值的公式。如果函数足够平滑的话，在已知函数在某一点的各阶导数值的情况下，泰勒公式可以用这些导数值做系统构建一个多项式来近似函数在这一点的邻域中的值。若$f(x)$在包含$x_0$的闭区间$[a,b]$上具有n阶导数，且在开区间$(a,b)$上具有n+1阶导数，则对闭区间$[a,b]$上任意一点x，成立：

$$f(x)=\frac{f(x_0)}{0!}+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+Rn(x)$$
在实际工作中，泰勒公式需要截断，只取有限项。

将上述目标函数${\cal L}^{(t)}$用泰勒二次展开近似得到：

$${\cal L}^{(t)}\approx\sum_{i=1}^n[l(y_i, \hat{y_i}^{(t-1)})+g_i\cdot f_t(x_i)+\frac{1}{2}h_i{f_t}^2(x_i)]+\Omega(f_t)+constant$$
令 $g_i=\partial_{\hat{y}^{(t-1)}}l(y_i,\hat{y_i}^{(t-1)})$， $h_i=\partial_{\hat{y}^{(t-1)}}^{2}l(y_i,\hat{y_i}^{(t-1)})$。即 $g_i$是一阶导， $h_i$是二阶导，都是关于上一轮的损失函数的导数。
为了计算方便，我们可以把式中的常数项移掉：
$${\cal L}^{(t)}\approx\sum_{i=1}^n[g_i\cdot f_t(x_i)+\frac{1}{2}h_i{f_t}^2(x_i)]+\Omega(f_t)$$

2.2 定义树结构

现在用数学方法来定义CART树的结构，对于第t棵树，存在一个映射函数能够把一个样本映射到某个叶子节点，这个方法为$q(x)$，存在一组向量$w$能够表示叶子节点的权值，那么第t个分类器可以表示为：

$$f_t(x)=w_{q(x)}$$
其中$w\in R^T$，$q:R^d\rightarrow \{1,2,...,T\}$

定义正则化项为(也可以采用别的可行的定义)：

$$\Omega(f_t)=\gamma T+\frac{1}{2}\lambda\sum_{j=1}^{T}w_j^2$$
正则化项囊括了叶子节点数和节点权值的L2范数。

假设第j个叶子节点对应的样本集合为$I_j=\{i|q(x_i)=j\}$，因为所有的样本映射到了某个叶子上，所以目标函数可以从样本求和转化为叶子的求和：

$${\cal L}^{(t)}=\sum_{i=1}^n[g_iw_{q(x_i)}+\frac{1}{2}h_iw_{q(x_1)}^2]+\gamma T+\lambda \frac{1}{2}\sum_{j=1}^Tw_j^2\\ =\sum_{j=1}^T[(\sum_{i\in I_j}g_i)w_j+\frac{1}{2}(\sum_{i\in I_j}h_i)w_j^2]+\gamma T+\frac{1}{2}\lambda \sum_{j=1}^Tw_j^2\\ =\sum_{j=1}^T[(\sum_{i\in I_j}g_i)w_j+\frac{1}{2}(\sum_{i\in I_i}h_i+\lambda)w_j^2]+\gamma T$$
目标函数转化为求一元二次方程的和。
令 $G_j=\sum_{i\in I_j}g_i$， $H_j=\sum_{i\in I_j}h_i$, 则上式转化为：
$${\cal L}^{(t)}=\sum_{j=1}^T[G_j+\frac{1}{2}(H_j+\lambda)w_j^2]+\gamma T$$
该式在 $w_j=-\frac{G_j}{H_j+\lambda}$时取得最小值，目标函数 ${\cal L}^{(t)}=-\frac{1}{2}\sum_{j=1}^T\frac{G_j^2}{H_j+\lambda}+\gamma T$。这个目标函数的值可以用来作为一个评分标准衡量树结构 $q$的好坏，这个值越小越好。

2.3 选择分裂点

一般情况下，穷举所有可能的树结构$q$是不太现实的，在实际应用中，我们用一种贪心算法来建立树结构：从树顶层开始，对于每一个节点，尝试增加一个分支。假设 $I_L$和 $I_R$分别是分裂点的左右节点的实例，则 $I=I_L\bigcup I_R$，现在来算一下因为分裂产生的增益：

$$Gain=\frac{1}{2}[\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}]-\gamma$$
式中 $\frac{G_L^2}{H_L+\lambda}$表示左子树分数， $\frac{G_R^2}{H_R+\lambda}$表示右子树分数， $\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}$表示如果不分裂我们可以得到的原分数，最后一项 $\gamma$表示加入新叶子节点引入的复杂度。其实 $\gamma$相当于阈值，它是正则项里叶子节点数T的系数，所以XGBoost在优化目标函数的同时相当于做了预剪枝，另外 $\lambda$是正则项里叶子权值的L2模平方的系数，对叶子权值做了平滑，也起到了防止过拟合的作用。
那么如何找到最好的分裂点呢？我们肯定是希望找到信息增益最大的分裂点，最简单的方法就是穷举(Exact Greedy Algorithm):

每个样本根据属性值的大小进行排列（共有m个属性，排列m次），对每个样本计算它的 $g_i,h_i$，从左到右尝试设置分裂点并计算增益，找到每个属性增益最大值的分裂点，取所有属性分裂点的增益最大的属性作为最先分裂点，以此来构造树结构。
Exact Greedy Algorithm是非常强大且精确的，因为它枚举了所有可能的分裂点，但是往往会造成很大的计算量和效率低下，尤其是样本量很大的时候要把所有数据读入内存可能会有困难。因此陈天奇博士在他的论文里提到了另一种寻找分裂点的近似方法——Approximate Algorithm。

简单地说，就是根据特征k的分布来确定 $l$个候选切分点$S_k$，然后根据这些候选切分点把相应的样本进行划分放入对应的桶中，对每个桶的 $G,H$进行累加，最后在候选切分点集合上贪心查找（类似于Exact Greedy Algorithm）。关于如何根据特征的分布来确定候选切分点，陈天奇博士提到了“Weighted Quantile Sketch”加权分位图，详情可以见附录，在这里不做展开。
那么何时进行分裂点的选取呢？XGBoost有两种策略，全局策略(global)和局部策略(local)。全局策略意味着在学习每棵树之前，提出候选切分点，然后每次分裂都使用相同的候选切分点；局部策略是在每次分裂前，重新提出候选切分点。可以看出，全局策略会比局部策略更简单点，但是因为每次分裂后没有调整候选切分点，全局策略会使用更多的候选切分点，局部策略会更适用于较深的树。

从上图可以看出来，全局策略的切分点个数够多的时候，和Exact Greedy Algorithm算法性能相当，局部策略的切分点个数不需要那么多，因为每一次分裂都重新进行了选择。

2.4 稀疏值处理

在实际操作中，因为缺失值，大量的0值和one-hot编码，我们的特征值往往是稀疏的。XGBoost能够对稀疏值进行自动处理，主要是通过对缺失值自动学习出它该被划分的方向（左子树或右子树）。那么我们该如何判断划分的方向？论文中提到了一个Sparsity-aware Split Finding算法：

简而言之，这个算法跟之前提到的Exact Greedy Algorithm类似，但是只遍历非缺失值。可以分类似的两步骤：
1) 将特征k的所有缺失值都放在右子树，然后枚举切分点，计算出最大的增益；
2) 将特征k的所有缺失值都放在左子树，然后枚举切分点，计算出最大的增益。
这样在求最大增益的同时，也知道了缺失值的样本应该放在左边还是右边。使用该方法，相当于比传统方法多遍历了一次，但是它只在非缺失值上进行迭代，因此复杂度与非缺失值的样本呈线性关系。经过测验，发现Sparsity-aware分裂算法比传统的算法快将近50倍。

3. 系统设计

---

3.1 分块并行

建树过程中最耗时的就是数据排序，为了减少时耗，原文中提到了用Block结构存储数据。Block中的数据以稀疏格式CSC (compressed column) 进行存储，并且按照特征值进行排序 (不对缺失值排序)，bBock中的特征还需要存储指向样本的索引，这样才能根据特征的值来取梯度，一个block中可以存储一个或多个特征的值。

可以看出，只需要在建树前排序一次，后面节点分裂时可以直接根据索引得到梯度信息。
$\cdot$在ExactGreedy算法中，将整个数据集存放在一个Block中，复杂度从原来的的$O(Hd||x||_0logn)$降为$O(Hd||x||_0+||x||_0logn)$，其中$||x||_0$为训练集中非缺失值的个数。这样，Exact greedy算法就省去了每一步中的排序开销。
$\cdot$在Approximate算法中，使用多个Block，每个Block对应原来数据的子集，不同的Block可以在不同的机器上进行计算。这个方法对local策略很有效，因为local策略每次分支都重新生成候选切分点。
Block结构还有其他好处，数据按列存储，可以同时访问所有的列，很容易实现并行的寻找切分点算法。缺点是空间消耗大。

3.2 缓存优化

使用Block结构的一个缺点是通过索引获取梯度时，梯度的获取顺序是按照特征的大小，这将导致非连续的内存访问，可能使得CPU cache缓存命中率低，从而影响算法效率。
因此，Exact Greedy Algorithm使用缓存预取，即对每个线程分配一个连续的Buffer，读取梯度信息并存入Buffer中，以实现非连续到连续的转化，然后再统计梯度信息，这种方法在训练样本数打的时候很管用。在Approximate Algorithm中，对Block的大小进行了合理的设置，定义Block的大小为Block中最多的样本数。设置过大的Block容易导致命中率低，过小则容易导致并行化效率不高，经过试验，发现2^16比较合适。

3.3 Out-of-core 计算

当数据量太大不能全部放入主内存的时候，将数据划分为多个Block并存放在磁盘上。主要的设计思想是：
1) 计算的时候，使用独立的线程预先将Block放入主内存，因此可以在计算的同时读取磁盘；
2) Block压缩，按列进行压缩，读取的时候用另外的线程解压。对于行索引，只保存第一个索引值，然后用16位的整数保存与该Block第一个索引的差值；
3) Block Sharding，将数据划分到不同的硬盘上，提高磁盘吞吐率。

4. XGBoost优点以及与GBDT的区别

---

1) 传统的GBDT以CART作为弱分类器，XGBoost还支持线性分类器，这时它相当于带L1和L2正则化项的逻辑斯蒂回归（分类问题）或者线性回归（回归问题）；
2) GBDT在优化时只用到了一阶导数信息，XGBoost则对代价函数进行了二阶泰勒展开，用了一阶和二阶导数；
3) XGBoost在代价函数里加入了正则项，用于控制模型的复杂度；
4) 列抽样，XGBoost借鉴了随机森林的做法，支持列抽样，不仅能降低过拟合，还能减少计算；
5) 对缺失值的处理，XGBoost可以自动学习出它的分裂方向；
6) 支持并行，由于XGBoost在训练前，预先对数据进行了排序，然后保存为Block结构，后面的迭代中重复地是用这个结构，大大减小了计算量，这个Block结构也使得并行成为了可能，在进行节点分裂时，需要计算每个特征的增益，各个特征的增益计算就可以用多线程进行。

5. API及参数介绍

---

介绍Scikit-Learn的API接口。

默认参数：
xgboost.XGBRegressor(max_depth=3, learning_rate=0.1, n_estimators=100, silent=True, objective=’reg:linear’, booster=’gbtree’, n_jobs=1, nthread=None, gamma=0, min_child_weight=1, max_delta_step=0, subsample=1, colsample_bytree=1, colsample_bylevel=1, reg_alpha=0, reg_lambda=1, scale_pos_weight=1, base_score=0.5, random_state=0, seed=None, missing=None, **kwargs)

parameters	定义
max_depth (*int*)	最大深度
learning_rate (*float*)	学习率
n_estimators (*int*)	弱学习器的数量
silent (*boolean*)	是否打印程序信息
objective (*string or callable*)	自定义的目标函数
booster (*string*)	弱分类器种类，gbtree, gblinear, dart
n_jobs (*int*)	线程数
gamma(*float*)	$\gamma$系数
min_child_weight (*int*)	叶子上的最小样本数
max_delta_step (*int*)	每棵树权重改变的最大步长
subsample (*float*)	每棵树随机采样的比例
colsample_bytree (*float*)	每棵树随机采样的列数的占比（每一列是一个特征）
colsample_bylevel (*float*)	用来控制树的每一级的每一次分裂，对列数的采样的占比
reg_alpha(*float*)	L1正则化项
reg_lambda (*float*)	L2正则化项
scale_pos_weight (*float*)	在各类别样本十分不均衡时，将参数设定为正值，可以使算法更快收敛
base_score	初始预测分数作为全局偏差
random_state (*int*)	随机种子
missing (*float*)	缺失值数据的设定

6. 附录

---

Weighted Quantile Sketch
章节2.3中提到了Approximate Algorithm中切分点$S_k=\{S_{k1}, S_{k2}, ...,S_{kl}\}$的选取。可以采用分位数，也可以直接构造梯度统计的近似直方图等。传统的分位数就是把概率分布划分为连续的区间，每个区间的概率相同，例如四分位数(把所有数值有小到大排列并分成四等份)。然而XGBoost采用的不是简单的分位数方法，而是对分位数进行加权(使用二阶梯度$h$)。

我们对特征$k$构造数据集：$D_k=(x_{1k}, h_1), (x_{2k}, h_2),...,(x_{nk}, h_n)$，其中$x_{ik}$表示样本i的特征$k$的值，$h_i$是对应的二阶梯度。定义一个rank function：

$$r_k(z)=\frac{1}{\sum_{(x,h)\in D_k}h}\sum_{(x,h)\in{D_k},x<z}h$$
该式表达了第 $k$个特征小于$z$的样本比例，类似于传统的分位数，不过这里按照二阶梯度进行统计。
候选切分点 $S_k=\{S_{k1}, S_{k2}, ...,S_{kl}\}$要求：
$$|r_k(s_{k,j})-r_k(s_{k, j+1})|<\varepsilon,\ \ \ \ \ \ S_{k1}=\min_ix_{ik}, \ \ S_{kl}=\max_ix_{ik}$$
即让相邻两个候选切分点相差不超过某个值 $\varepsilon$，因此总共会得到 $1/\varepsilon$个切分点。

为什么要用二阶梯度加权？
原因可以通过对泰勒二阶展开后的目标函数进行转化配方得出：

$${\cal L}^{(t)}\approx\sum_{i=1}^n[g_i\cdot f_t(x_i)+\frac{1}{2}h_i{f_t}^2(x_i)]+\Omega(f_t)\\ =\sum_{i=1}^N \frac{1}{2}h_i[2\frac{g_i}{h_i}f_t(x_i)+f_t^2(x_i)]+\Omega(f_t)\\ =\sum_{i=1}^N \frac{1}{2}h_i[\frac{g_i^2}{h_i^2}+2\frac{g_i}{h_i}f_t(x_i)+f_t^2(x_i)]+\Omega(f_t)\\ =\sum_{i=1}^N \frac{1}{2}h_i[f_t(x_i)-(-\frac{g_i}{h_i})]^2+\Omega(f_t)$$
(注：因为 $\frac{g_i^2}{h_i^2}是常量，所以可以加入到第三行$)
从上述推导中可以看出目标函数可以转化为权重为 $h_i$的平方损失，因此用 $h_i$加权。

7. 参考资料

---

$\cdot$ Chen, Tianqi, and Carlos Guestrin. “Xgboost: A scalable tree boosting system.” Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining. ACM, 2016.
$\cdot$ xgboost算法原理与实战
$\cdot$ 集成学习（三）XGBoost
$\cdot$ 知乎：机器学习算法中 GBDT 和 XGBOOST 的区别有哪些？
$\cdot$ XGBoost参数调优完全指南（附Python代码）