XGBoost原理理解

作者: 单调不减 | 来源:发表于2019-11-14 23:19 被阅读0次

XGBoost = Extreme Gradient Boosting，也就是说，XGBoost 把 Gradient Boosting 做到了极致，具体来说，这种极致体现在训练精度高、所需计算资源少。正如陈天奇大神在 Quora 的相关回答中说：

The name xgboost, though, actually refers to the engineering goal to push the limit of computations resources for boosted tree algorithms. Which is the reason why many people use xgboost.

下图是基于树的算法的发展路径图，可以看到 XGBoost 是 Gradient Boosting 的（终极）强化版（图自here）。

接下来我们按照原论文的组织顺序来大概记录一下对 XGBoost 的原理理解。

1、Model

首先模型可以表示如下：

$\hat{y}_{i}=\phi\left(\mathbf{x}_{i}\right)=\sum_{k=1}^{K} f_{k}\left(\mathbf{x}_{i}\right), \quad f_{k} \in \mathcal{F}$

这里各个 $f_k$ 表示回归树模型， $\mathcal{F}=\left\{f(\mathrm{x})=w_{q(\mathrm{x})}\right\}\left(q: \mathbb{R}^{m} \rightarrow T, w \in \mathbb{R}^{T}\right)$ 表示回归树空间。 $q$ 表示一个树结构对应的映射，这个映射是将一个样本映射到所划分的叶子结点的编号。 $w_{q(x)}$ 则表示样本所属结点的得分或权重。 $T$ 表示树的叶子结点数。

听起来很复杂的样子，实际上就是建立 $K$ 棵回归树，然后对一个输入 $X$ ，在 $K$ 颗树中都找到其所属叶子结点的得分，然后加起来就得到最终的预测 $\hat{y}$ 。原论文中有清晰的图解：

2、Loss Function & Quality Score

有了模型，下一步就是定义损失函数：

$\begin{aligned} & \mathcal{L}(\phi) =\sum_{i} l\left(\hat{y}_{i}, y_{i}\right)+\sum_{k} \Omega\left(f_{k}\right) \\ &\text { where }\ \ \ \Omega(f) =\gamma T+\frac{1}{2} \lambda\|w\|^{2} \\ \end{aligned}$

这里 $l$ 是一个可导的凸损失函数，衡量模型预测值和真实值之间的差异， $\Omega$ 是针对模型复杂度的罚项，可以视作正则项。这里的 $\Omega(f)$ 由两部分组成，第一项是限制叶子结点数，第二项则是限制各叶子结点的分数，目的都是防止过拟合。删去正则项则损失函数退化至普通 Gradient Boosting 的情形。

需要注意的是，损失函数 $\mathcal{L}(\phi)$ 的输入是一个函数，因此不能用像SGD这样的方法去找到各个 $f_k$ ，因为他们是树而不是数值向量。解决问题的方法就是加法训练(Additive Training)。定义 $\hat{y_{i}}^{(t)}$ 为第 $i$ 个样本在第 $t$ 轮迭代得到的模型上做出的预测。则假设当前我们已经进行了 $t-1$ 轮迭代，则我们希望得到一个 $f_t$ 使得下面的损失函数最小：

$\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)$

我们用 Taylor 公式对上式进行二次近似，之所以用二次而不是一次，是为了优化时使得损失更快下降。这里我们回顾一下泰勒公式：

$f(x+\Delta x) \approx f(x) + f^{'}(x)\Delta x + \frac{1}{2}f^{''}(x)\Delta x^2$

把 $\mathcal{L}^{(t)}$ 中的 $\hat{y_{i}}^{(t-1)}$ 当作上式中的 $x$ ， $f_{t}\left(\mathbf{x}_{i}\right)$ 当作上式中的 $\Delta x$ ，则：

$\mathcal{L}^{(t)} \simeq \sum_{i=1}^{n}\left[l\left(y_{i}, \hat{y}^{(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)$ 为二阶导。值得注意的是，如果这里只做一次近似，那么当 $f_t(X)$ 沿 $-g$ 方向下降最快，因此可以用 $f_t(X)$ 拟合 $-g$ ，这实际上就是 Gradient Boosting。

由于 $l\left(y_{i}, \hat{y}^{(t-1)}\right)$ 和当前要优化的函数 $f_t$ 无关，因此可以作为常数项去掉，从而：

$\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\}$ 为第 $j$ 个叶结点包含的样本，则上式可以写作：

$\begin{aligned} \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{aligned}$

对于固定的树结构 $q(x)$ ，可以计算出各叶子结点的最优分数：

$w_{j}^{*}=-\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$

于是上式就可以作为衡量一个树结构 $q$ 好坏的分数，这个分数类似于决策树中的基尼系数或信息增益，可以用于判断结点如何分裂。

假设当前树中一个叶子结点的样本集为 $I$ ，分裂后得到的左孩子和右孩子的样本为 $I_L$ 和 $I_R$ ，则结点分裂后得到的树结构的最优损失为：

$\tilde{\mathcal{L}}^{(t)}(q_{split})=-\frac{1}{2} \sum_{j=1}^{T} \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} \right]+\gamma (T+1)$

上面两式相减，损失减少量为：

$\mathcal{L}_{s p l i t}=\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$

当然遍历所有可能的树结构是不可能的，我们只能采用贪心的算法，即从初始的根节点开始，每一步都采用当前的最优结点划分策略，直至达到终止条件。这部分的算例原论文也有给出：

而上面所述的贪心算法的正规化表达如下：

3、Additional Techniques

在防止过拟合的措施方面，除了上面 loss function 中正则项对树模型的限制，还可以使用 Shrinkage 和 Column Subsampling 的方法。

所谓 Shrinkage，就是在每一轮添加一棵树的时候加入一个收缩系数 $\epsilon$ ，类似于梯度下降算法中的学习率：

$\hat{y}_{i}^{(t)}=\hat{y}_{i}^{(t-1)}+\epsilon f_{t}\left(x_{i}\right)$

收缩系数 $\epsilon$ 通常设置为大约 0.1。这样做的原因很简单，就是减小单个树模型对结果的影响，并为未来的树对模型的改进留足空间，从而有效缓解过拟合。

所谓 Column Subsampling ，其实就是 Feature Subsampling，特征抽样。这个防止过拟合的方法就是 Random Forest 中所采用的，论文中提到，Column Subsampling 在防止过拟合方面的效果通常优于 Row Subsampling (Sample Subsampling)。

4、Approximate Algorithm

上面提到的贪心算法是确定性的，因为每次结点划分前，我们都要遍历各个特征以及各特征可能的划分点，所以计算消耗较大，尤其是数据分布很集中的情形会造成很多不必要的计算，且对于数据量大到无法全部存在内存中的情形，计算效率很低（因为每个结点判断如何划分的过程中都要进行内存的存取操作）。

为了加速计算，论文中提出了相应的近似算法：

概括一下就是，枚举所有特征，根据特征，比如第 $k$ 个特征的分位数来决定出 $l$ 个候选切分点 $S_k=\{s_{k1},s_{k2},⋯s_{kl}\}$ ，然后根据这些候选切分点把相应的样本映射到对应的桶中，对每个桶的 $G,H$ 进行累加。最后在候选切分点集合上贪心查找。

下面是一个三分位数的例子，图自知乎wepon。

论文中提到了该近似算法的两种形式，一种是 global，另一种是 local。

所谓 global，就是在建树之前预先将数据进行全局（global）分桶，所谓 local，就是在各个结点分裂时重新局部（local）分桶。

直觉上讲，global 需要设置更小的分位数间隔（这里用 $ϵ$ 表示，3分位的分位数间隔就是 $\frac{1}{3}$ ），产生更多的桶。而 local 可以设置较大的 $ϵ$ ，产生更少的桶。

论文给出了两种近似算法效果的比较：

可以看到，两种近似算法在设置合适的分位数间隔的情况下都可以达到与确定性算法相当的效果，这证明了近似算法的有效性。此外，论文中还提到，local 算法更适合树较深的场景。

5、Weighted Quantile Sketch

上面我们提到的近似算法中很重要的步骤就是对各个特征找到合适的分位数作为候选划分点，从所举的3分位数的例子中我们也可以看到，对于各个样本我们是一视同仁的，也就是说，我们在计算分位数的时候把各个样本的权重都视为1。论文中提出了一种更加合理的 Weighted Quantile Sketch，即带权的分位方案。

记 $D_k = {(x_{1k}, h_1),(x_{2k}, h_2)· · ·(x_{nk}, h_n)}$ 为各样本的第 $k$ 个特征的值以及其对应的二阶梯度构成的点对集合。

在此基础上，我们可以定义一个排序函数 $r_k:R→[0,1]$ ：

$r_{k}(z)=\frac{1}{\sum_{(x, h) \in \mathcal{D}_{k}} h} \sum_{(x, h) \in \mathcal{D}_{k}, x<z} h$

上式表示的就是该特征的值小于 $z$ 的样本权重和占总样本权重和的比例。我们发现，这里我们把二阶梯度的值 $h$ 作为样本的权重，为什么这样做接下来会提到。

于是我们就能用下面这个不等式来寻找分裂候选点 $\{s_{k1},s_{k2},s_{k3},⋯,s_{kl}\}$ ：

$\left|r_{k}\left(s_{k, j}\right)-r_{k}\left(s_{k, j+1}\right)\right|<\epsilon, \quad s_{k 1}=\min _{i} \mathbf{x}_{i k}, s_{k l}=\max _{i} \mathbf{x}_{i k}$

设定合适的 $ϵ$ ，并控制相邻两个候选分裂点相差不超过某个值 $ϵ$ ，则 $\frac{1}{ϵ}$ 的整数值就代表几分位，比如 $ϵ=\frac{1}{3}$ 就是三分位，即有 2 个候选分裂点。

下面我们来看一下为什么用二阶导 $h$ 作为样本的权重是合理的，回到之前我们定义的 loss function 并做一下变形：

$\begin{aligned} O b j^{(t)} & \approx \sum_{i=1}^{n}\left[g_{i} f_{t}\left(x_{i}\right)+\frac{1}{2} h_{i} f_{t}^{2}\left(x_{i}\right)\right]+\Omega\left(f_{t}\right) \\ &=\sum_{i=1}^{N} \frac{1}{2} h_{i}\left(2 \frac{g_{i}}{h_{i}} f_{t}\left(\mathbf{x}_{i}\right)+f_{t}^{2}\left(\mathbf{x}_{i}\right)\right)+\Omega\left(f_{t}\right) \\ &=\sum_{i=1}^{N} \frac{1}{2} h_{i}\left(\frac{g_{i}^{2}}{h_{i}^{2}}+2 \frac{g_{i}}{h_{i}} f_{t}\left(\mathbf{x}_{i}\right)+f_{t}^{2}\left(\mathbf{x}_{i}\right)\right)+\Omega\left(f_{t}\right)-\frac{g_{i}^{2}}{2 h_{i}} \\ &=\sum_{i=1}^{N} \frac{1}{2} h_{i}\left(f_{t}\left(\mathbf{x}_{i}\right)-\left(-\frac{g_{i}}{h_{i}}\right)\right)^{2}+\Omega\left(f_{t}\right)-\frac{g_{i}^{2}}{2 h_{i}} \\ &=\sum_{i=1}^{N} \frac{1}{2} h_{i}\left(f_{t}\left(\mathbf{x}_{i}\right)-\left(-\frac{g_{i}}{h_{i}}\right)\right)^{2}+\Omega\left(f_{t}\right)-\text {constant } \end{aligned}$

注意，这里我们说 $\frac{g_i^2}{2h_i^2}$ 是常数是因为 $g_i$ 和 $h_i$ 都是由上一轮迭代的结果得到的，与当前的 $f_t$ 无关。此外，这里原论文中把 $-\frac{g_i}{h_i}$ 写成了 $\frac{ g_i}{h_i}$ ，少了一个负号，这里应该是一个小小的错误，stackexchange上的相关回答也指出了这一点。

总之，变形后的 loss function 可以视作标签为 $-\frac{g_i}{h_i}$ ，权重为 $h_i$ 的平方损失。

那么为什么我们在选取候选分裂点时要考虑平方损失的系数呢？一种直观的理解是，越是权重大的样本我们越要谨慎对待，因为同样的预测偏差，权重大的样本会带来更大的损失，因此对权重大的样本我们应该细粒度划分，即尽可能考虑每个大权重样本应该划入左孩子还是右孩子结点，而对小权重的样本则不必这么谨慎，可以采取粗粒度划分。

一个实际的例子如下（图自知乎wepon）：

不难看出，随着样本权重的增加，我们划分的粒度越来越细。

6、 Sparsity-aware Split Finding

在实际问题中，数据稀疏性是十分常见的，可能的原因有：1）数据缺失；2）数据零值过多；3）特征工程的结果，比如One-hot编码。

在训练集出现大量数据稀疏的情况下，XGBoost 在建树时可能由于当前样例数据值缺失无法判断被分配到哪个分支，最直接的想法就是给定一个默认方向，当数据值缺失时直接将样例分配给默认方向的分支。至于默认方向是左还是右则需从数据集中学习得到。

论文中给出的算法如下：

上述算法的主要特点在于只需访问数据值无缺失的样例，算法的计算复杂度与无缺失数据集的大小呈线性关系，在数据稀疏的情况下能大大减少计算量。

可以看到，我们对每个特征 $k$ 都做两次遍历，分别按照无缺失数据集的升序和降序来遍历无缺失的值，从中寻找最优的分裂点。因为算法在计算 $G$ 和 $H$ 时既包含了缺失值的样例，又包括了无缺失值的样例，而遍历的时候只包含了无缺失值的样例，所以两次遍历的结果是不一样的。

我们把当前结点的数据集 $I$ 划分成在特征 $k$ 上无缺失值的 $I_k$ 和有缺失值的 $I_{k_{NA}}$ ，则上述算法考虑两种情形，第一种是将 $I_k$ 按升序排列，并遍历候选分裂点，得到 $I_L$ 和 $(I_{k_{NA}},I_R)$ 两个分支，然后计算这样分裂结点的得分，找到其中的最大得分作为默认方向为右的得分；第二种是将 $I_k$ 按降序排列，并遍历候选分裂点，得到 $(I_{k_{NA}},I_L)$ 和 $I_R$ 两个分支，然后计算这样分裂结点的得分，找到其中的最大得分作为默认方向为左的得分。最后比较这两个得分，选择得分大的方向作为默认方向。