集成学习之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 的学习算法：
- 初始化 \(f_0(x)\)
- 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
- Tianqi Chen and Carlos Guestrin. XGBoost: A Scalable Tree Boosting System
- Tianqi Chen. Introduction to Boosted Trees
- https://zhuanlan.zhihu.com/p/34534004
- https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/
- https://*.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}
\]
/