集成学习之Boosting —— XGBoost

集成学习之Boosting —— AdaBoost

集成学习之Boosting —— Gradient Boosting

集成学习之Boosting —— XGBoost

Gradient Boosting 能够看作是一个整体的算法框架，起始于Friedman 的论文 [Greedy Function Approximation: A Gradient Boosting Machine] 。XGBoost (eXtreme Gradient Boosting) 是于2015年提出的一个新的 Gradient Boosting 实现，由华盛顿大学的 陈天奇 等人开发，在速度和精度上都有显著提高，于是近年来在 Kaggle 等各大数据科学比赛中都获得了普遍应用。本文主要对其原理进行阐述，并将其与传统的 GBDT 进行比较。

大致来看，XGBoost 在原理方面的改进主要就是在损失函数上做文章。一是在原损失函数的基础上添加了正则化项产生了新的目标函数，这相似于对每棵树进行了剪枝并限制了叶结点上的分数来防止过拟合。二是对目标函数进行二阶泰勒展开，以相似牛顿法的方式来进行优化（事实上早在 [Friedman, J., Hastie, T. and Tibshirani, R., 1999] 中就已有相似方案，即利用二阶导信息来最小化目标函数，陈天奇在论文中也提到了这一点）。

在上一篇文章中，了解到 Gradient Boosting 的思想能够理解为经过函数空间的梯度降低来最小化目标函数 \(L(f) = \sum\limits_{i=1}^NL(y_i,f_m(x_i))\)，其只使用了一阶导信息，而 XGBoost 引入二阶导的一大好处是能够推导出一种新的增益计算方法，事实证实采用新的增益计算方法在优化目标函数上更加有效，精确度上也赛过传统的 GBDT。因此下面也从目标函数的优化入手进行推导。

XGBoost的目标函数

与 GBDT 同样，XGBoost 一样采用加法模型，设基学习器为\(f(x)\)，预测值为\(\hat{y}\)：
\[ \hat{y} = \sum\limits_{m=1}^M f_m(x)\;, \quad f_m \in \mathcal{F} \tag{1.1} \]
在第m步，前m-1个基学习器是固定的，于是目标函数为
\[ \mathcal{L}^{(m)} = \sum\limits_{i=1}^N L(y_i, \;\hat{y_i}^{(m-1)} + f_m(x_i)) + \Omega(f_m) \tag{1.2} \]
在传统的 GBDT 中，是没有正则化项 \(\Omega\) 的，而在XGBoost中采用的是
\[ \Omega(f) = \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \]
其中 \(J\) 为叶结点数目，\(b_j\) 为各个叶结点上的值。该项限制了树的复杂度，这样至关于使叶结点的数目变小，同时限制叶结点上的分数，由于一般分数越大学得越快，就越容易过拟合。

接下来将 \(\Omega\) 代入\((1.2)\) 式并对目标函数在 \(\hat{y}_i^{(m-1)}\) 处进行二阶泰勒展开：
\[ \mathcal{L}^{(m)} \simeq \sum\limits_{i=1}^N \left[ L(y_i, \hat{y}_i^{(m-1)}) + g_if_m(x_i) + \frac12 h_i f_m^2(x_i)\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \tag{1.3} \]
其中 \(g_i = \frac{\partial L(y_i, \,\hat{y}_i^{m-1})}{\partial \, \hat{y}^{(m-1)}}\;\;,\;\; h_i = \frac{\partial^2 L(y_i, \,\hat{y}_i^{(m-1)})}{(\partial\, \hat{y}^{(m-1)})^2}\)

在上一篇 Gradient Boosting 中提到决策树将特征空间划分为各个独立区域，每一个样本只属于其中一个区域，于是单颗决策树可表示为 \(f(x) = \sum\limits_{j=1}^J b_j I(x \in R_j)\)，若是将全部样本加起来，则可表示为 \(\sum\limits_{i=1}^N f(x_i) = \sum\limits_{j=1}^J \sum\limits _{x \in R_j} b_j\) ，代入 \((1.3)\) 式并将常数项 \(L(y_i, \hat{y}_i^{(m-1)})\) 移去：
\[ \begin{align*} \mathcal{L}^{(m)} &= \sum\limits_{i=1}^N \left[ g_if_m(x_i) + \frac12 h_i f_m^2(x_i)\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \\ &= \sum\limits_{j=1}^J \left[ \sum\limits_{x_i \in R_j}g_ib_j + \frac12 \sum\limits_{x_i \in R_j} h_i b_j^2\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \\ &= \sum\limits_{j=1}^J \left[ \sum\limits_{x_i \in R_j}g_ib_j + \frac12 (\sum\limits_{x_i \in R_j} h_i + \lambda) \,b_j^2 \right] + \gamma \, J \\ & = \sum\limits_{j=1}^J(G_j b_j + \frac12 (H_j + \lambda) \,b_j^2) + \gamma\, J \tag{1.4} \end{align*} \]

其中 \(G_j = \sum\limits_{x_i \in R_j} g_i \;\; , \;\; H_j = \sum\limits_{x_i \in R_j} h_i\)

XGBoost的增益计算和树分裂方法

通过一系列推导后，如今咱们的目标变为最小化 \((1.4)\) 式，并求得相应的树结构 \(\{R_j\}^J_1\) 和叶结点上的值 \(\{b_j\}^J_1\) 。

精确地划分 \(\{R_j\}^J_1\) 是一个 NP hard 问题，现实中常使用贪心法，遍历每一个特征的每一个取值，计算分裂先后的增益，并选择增益最大的进行分裂。


对于具体增益的衡量标准，在几种决策树算法中，ID3 采用了信息增益：
\[ I(D,A) = H(D) - H(D|A) = H(D) - \sum\limits_{v=1}^V\frac{|D^v|}{|D|} H(D^v) \]
其中 V 表示特征 A 有 V 个可能的取值，\(D^v\) 表示第 v 个取值上的样本数量。

C4.5 采用了信息增益比：
\[ I_R(D,A) = \frac{I(D, A)}{H_A(D)} \]
其中 \(H_A(D) = -\sum\limits_{v=1}^V \frac{|D^v|}{|D|} log_2 \frac{|D^v|}{|D|}\) 。

CART 分类树采用了基尼系数：
\[ Gini(D,A) = \sum\limits_{v=1}^V \frac{|D^v|}{|D|} Gini(D^v) = \sum\limits_{v=1}^V \frac{|D^v|}{|D|} \left(1 - \sum\limits_{k=1}^K \left\lbrace\frac{|C_k|}{|D|}\right\rbrace^2\right) \]
其中 K 为类别个数，\(|C_k|\) 为 \(D\) 中属于第k类的样本数量。

CART 回归树采用了均方偏差：
\[ \mathop{min}_{A,s}\Bigg[\mathop{min}_{c_1}\sum\limits_{x_i \in R_1(A,s)}(y_i - c_1)^2 + \mathop{min}_{c_2}\sum\limits_{x_i \in R_2(A,s)}(y_i - c_2)^2\Bigg] \]
其中 \(c_1\)为特征 A 的切分点 s 划分出来的 \(R_1\) 区域的样本输出均值，\(c_1 = mean(y_i | x_i \in R_1(A,s))\)，\(c_2\)为 \(R_2\) 区域的样本输出均值，\(c_2 = mean(y_i | x_i \in R_2(A,s))\) 。

而XGBoost则提出了一种新的增益计算方法。

若是已经肯定了树的结构 \(\{R_j\}^J_1\) ，则直接对 \((1.4)\) 式求导，最优值为：
\[ b^*_j = - \frac{G_j}{H_j + \lambda} \]
再代入\((1.4)\)式获得此时最小的为目标函数为：
\[ \mathcal{L^{(m)}} = -\frac12 \sum\limits_{j=1}^J \frac{G_j^2}{H_j+ \lambda} + \gamma\, J \tag{1.5} \]
式 \((1.5)\) 能够认为是 XGBoost 的打分函数，该值越小，说明树的结构越好，下图示例了该式的计算方法：

该式的优势是除了能做为树分裂的衡量标准外，还能使 XGboost 适用于各类不一样的损失函数，因此 XGBoost 包中支持自定义损失函数，但前提是一阶和二阶可导。

从另外一角度看， \((1.5)\) 式就相似于传统决策树中的不纯度指标，在决策树中咱们但愿一次分裂后两边子树的不纯度越低越好，对应到XGBoost中则是但愿 \((1.5)\)式 越小越好，即 \(\frac{G_j^2}{H_j+ \lambda}\) 越大越好，这样分裂先后的增益定义为：
\[ Gain = \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 \tag{1.6} \]
\(\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} $ 为分裂前的分数，\(Gain\) 越大说明越值得分裂。固然 \((1.6)\) 式中 \(Gain\) 可能会变成负的，这个时候分裂后的目标函数不会减小，但这并不意味着不会分裂 。事实上 XGBoost 采用的是后剪枝的策略，建每棵树的时候会一直分裂到指定的最大深度(max_depth)，而后递归地从叶结点向上进行剪枝，对以前每一个分裂进行考察，若是该分裂以后的 \(Gain \leqslant 0\)，则咔咔掉。 \(\gamma\) 是一个超参数，具备双重含义，一个是在 \((1.3)\) 式中对叶结点数目进行控制的参数；另外一个是 \((1.6)\) 式中分裂先后 \(Gain\) 增大的阈值，固然两者的目的是同样的，即限制树的规模防止过拟合。

接下来考察决策树创建的过程。若是是使用贪心法，就是遍历一个叶结点上的全部候选特征和取值，分别计算 \(Gain\) ，选择 \(Gain\) 最大的候选特征和取值进行分裂，以下树分裂算法流程 (注意这是单个叶结点的分裂流程)：

有了上述单个叶结点上的分裂流程后，咱们能够总结下整个 XGBoost 的学习算法：

1. 初始化 \(f_0(x)\)
2. for m=1 to M :
   (a) 计算损失函数在每一个训练样本点的一阶导数 \(g_i = \frac{\partial L(y_i, \,\hat{y}_i^{m-1})}{\partial \, \hat{y}^{(m-1)}}\) 和二阶导数 \(h_i = \frac{\partial^2 L(y_i, \,\hat{y}_i^{(m-1)})}{(\partial\, \hat{y}^{(m-1)})^2}\) , $ i = 1,2 \cdots N$
   (b) 递归地运用树分裂算法生成一颗决策树 \(f_m(x)\)
   (c) 把新生成的决策树添加到模型中， $\hat{y_i}^{(m)} = \hat{y_i}^{(m-1)} + f_m(x) $

若是把上述 XGBoost 的学习算法和上一篇中传统 GBDT 的学习算法做比较，XGBoost 的主要优点是在损失函数中加入正则化项后使得学习出来的树更加简单，防止过拟合，但除此之外并不能体现出 XGBoost 的速度优点。XGBoost 之因此快的一大缘由是在工程上实现了 Column Block 方法，使得并行训练成为了可能。

Column Block for Parallel Learning

对于决策树来讲，对连续值特征进行划分一般比较困难，由于连续值特征每每取值众多。一般的作法是先按特征取值对样本排序，再按顺序计算增益选择划分点。若每次分裂都要排一次序，那是很是耗时的，因此 XGBoost 的作法是在训练以前，预先按特征取值对样本进行了排序，而后保存为block结构，采用CSC格式存储，每一列(一个特征列)均升序存放，这样，一次读入数据并排好序后，之后都可重复使用，大大减少计算量。

因为已经预先排过序，因此在树分裂算法中，经过一次线性扫描 (linear scan) 就能找出一个特征的最优分裂点，以下图所示，同时能够看到缺失值并无保存：

陈天奇论文里有关键的一句话：Collecting statistics for each columns can be parallelized 。因为已经预先保存为block 结构，因此在对叶结点进行分裂时，每一个特征的增益计算就能够开多线程进行，训练速度也由此提高了不少。并且这种 block 结构也支持列抽样，只要每次从全部 block 特征中选择一个子集做为候选分裂特征就能够了，据个人使用经验，列抽样大部分时候都比行抽样的效果好。

近似分裂算法

当数据量很是大难以被所有加载进内存时或者在分布式环境下时，上文的树分裂算法将再也不适用，于是做者提出了近似分裂算法，不只解决了这个问题，也同时能提高训练速度，具体流程见下图：

近似分裂算法其实就是对连续性特征进行离散化，对于某个特征 \(k\)，算法首先根据特征分布找到 \(l\) 个分位点的候选集合 \(S_k = \{s_{k1}, s_{k2}, ... ,s_{kl} \}\) ，而后根据集合 \(S_k\) 中的分位点将相应样本划分到桶(bucket)中。在遍历该特征时，只需遍历各个分位点，对每一个桶内的样本统计值 \(g\)、\(h\) 进行累加统计，寻找最佳分裂点进行分裂。该算法可分为 global 近似和 local 近似，global 就是在新生成一棵树以前就对各个特征计算分位点并划分样本，以后在每次分裂过程当中都重复利用这些分位点进行划分，而 local 就是每一次结点分裂时都要从新计算分位点。

总结

最后总结一下 XGBoost 与传统 GBDT 的不一样之处：

• 传统 GBDT 在优化时只用到一阶导数信息，XGBoost 则对目标函数进行了二阶泰勒展开，同时用到了一阶和二阶导数。另外 XGBoost 工具支持自定义损失函数，只要函数可一阶和二阶求导。

• XGBoost 在损失函数中加入了正则化项，用于控制模型的复杂度，防止过拟合，从而提升模型的泛化能力。

• 传统 GBDT 采用的是均方偏差做为内部分裂的增益计算指标（由于用的都是回归树），而 XGBoost 使用的是通过优化推导后的式子，即式 \((1.6)\) 。

• XGBoost 借鉴了随机森林的作法，支持列抽样，不只能下降过拟合，还能减小计算量，这也是 XGBoost 异于传统 GBDT 的一个特性。

• XGBoost 添加了对稀疏数据的支持，在计算分裂增益时不会考虑带有缺失值的样本，这样就减小了时间开销。在分裂点肯定了以后，将带有缺失值的样本分别放在左子树和右子树，比较二者分裂增益，选择增益较大的那一边做为默认分裂方向。

• 并行化处理：因为 Boosting 自己的特性，没法像随机森林那样树与树之间的并行化。XGBoost 的并行主要体如今特征粒度上，在对结点进行分裂时，因为已预先对特征排序并保存为block 结构，每一个特征的增益计算就能够开多线程进行，极大提高了训练速度。

• 传统 GBDT 在损失再也不减小时会中止分裂，这是一种预剪枝的贪心策略，容易欠拟合。XGBoost采用的是后剪枝的策略，先分裂到指定的最大深度 (max_depth) 再进行剪枝。并且和通常的后剪枝不一样， XGBoost 的后剪枝是不须要验证集的。 不过我并不以为这是“纯粹”的后剪枝，由于通常仍是要预先限制最大深度的呵呵。

说了这么多 XGBoost 的优势，其固然也有不完美之处，由于要在训练以前先对每一个特征进行预排序并将结果存储起来，对于空间消耗较大。另外虽然相比传统的 GBDT 速度是快了不少，但和后来的 LightGBM 比起来仍是慢了很多，不知之后还会不会出现更加快的 Boosting 实现。

---

Reference

1. Tianqi Chen and Carlos Guestrin. XGBoost: A Scalable Tree Boosting System
2. Tianqi Chen. Introduction to Boosted Trees
3. https://zhuanlan.zhihu.com/p/34534004
4. https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/
5. https://stackoverflow.com/questions/52672116/what-is-xgboost-pruning-step-doing

附录： 泰勒公式

XGBoost 推导的关键一步是二阶泰勒展开，这里做一下拓展。泰勒公式的简单解释就是用多项式函数去逼近原函数，由于用多项式函数每每求解更加容易，而多项式中各项的系数则为原函数在某一点的n阶导数值除以n阶乘。这里力荐 3Blue1Brown 关于泰勒公式的视频 微积分的本质 - 泰勒级数 ，讲得很是形象 。

已知函数 \(f(x)\) 在 \(x=x_0\) 处n阶可导，那么：
\[ \begin{aligned} f(x) &= f(x_0) + f'(x_0)(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 + \cdots + \frac{f^{(n)}(x_0)}{n!}(x-x_0)^n \\ & = \sum\limits_{n=0}^{\infty}\frac{f^{(n)}x_0}{n!}(x - x_0)^n \end{aligned} \]
例如，\(x_0 = 0, \;\; f(x) = e^x\)时，\(e^x = \sum\limits_{n=0}^{\infty}\frac{x^n}{n!} = 1+x+\frac{x^2}{2!} + \frac{x^3}{3!} + \cdots\)

在机器学习中泰勒公式的一个典型应用就是牛顿法。在牛顿法中，将损失函数 \(L(\theta)\) 在 \(\theta^{k}\) 处进行二阶泰勒展开(考虑 \(\theta\) 为标量)：
\[ L(\theta) =L(\theta^k) + L'(\theta^k)(\theta - \theta^k) + \frac{L''(\theta^k)}{2}(\theta - \theta^k)^2 \]
将一阶导和二阶导分别记为 \(g_k\) 和 \(h_k\)，为使上式最小化，令
\[ \frac{\partial\left[g(\theta - \theta^k) + \frac{h}{2}(\theta - \theta^k)^2\right]}{\partial \,\theta} = 0, \qquad 求得\; \theta-\theta^k = -\frac{g_k}{h_k}\;, \]
则参数更新公式为 \(\theta^{k+1} = \theta^k - \frac{g_k}{h_k} \;\)，推广到向量形式则为 \(\boldsymbol{\theta}^{k+1} = \boldsymbol{\theta}^k - \mathbf{H}^{-1}_k\mathbf{g}_k\;\)，\(\mathbf{H}_k\) 为海森矩阵(Hessian matrix)，\(\mathbf{g}_k\) 为梯度向量：
\[ \nabla_\theta L = \mathbf{g}_k = \begin {bmatrix} \frac{\partial L}{\partial \theta_1}\\ \frac{\partial L}{\partial\theta_2}\\ \vdots \\ \frac{\partial L}{\partial\theta_n} \end {bmatrix} \;\;,\;\; \nabla_{\theta}^2 L = \mathbf{H}_k = \begin {bmatrix} \frac{\partial^2 L}{\partial\, \theta_1^2} & \frac{\partial^2 L}{\partial\,\theta_1 \partial\,\theta_2} & \cdots & \frac{\partial^2 L}{\partial\,\theta_1 \partial\,\theta_n} \\ \frac{\partial^2 L}{\partial\,\theta_2 \partial\,\theta_1} & \frac{\partial^2 L}{\partial\, \theta_2^2} & \cdots & \frac{\partial^2 L}{\partial\,\theta_2 \partial\,\theta_n} \\ \vdots & \vdots & \ddots \\ \frac{\partial^2 L}{\partial\,\theta_n \partial\,\theta_1} & \frac{\partial^2 L}{\partial\,\theta_n \partial\,\theta_2} & \cdots & \frac{\partial^2 L}{\partial\, \theta_n^2} \end {bmatrix} \]

/