Gradient Boosting Decision Tree

---

借鉴博客 GBDT原理与Sklearn源码分析-回归篇 （★★★★★）

---

1 GBDT概述

  梯度提升树属于Boosting集成学习算法的一种，其思想不同于随机森林、Bagging的并行化、投票的流程，GBDT模型所输出的结果是由其包含的若干棵决策树累加而得到的，每一棵子决策树都是实现对先前决策树组预测残差的拟合，是对先前模型的结果的一种“修正”。梯度提升树既可以用于回归问题（此时被称为CART回归树），也可以被用于解决分类问题（GBDT分类树）。本文主要介绍GBDT回归的原理！

  在GBDT的迭代中，假设我们前一轮迭代得到的强学习器是$f_{m-1}(x)$, 损失函数是$L(y,f_{m-1}(x))$， 我们本轮迭代的目标是找到一个CART回归树模型的弱学习器$h_m(x)$，让本轮的损失$L(y,f_m(x))=L(y,f_{m-1}(x)+h_m(x))$最小。也就是说，本轮迭代找到决策树，要让样本的损失尽量变得更小。

  GBDT的思想可以用一个通俗的例子解释，假如有个人30岁，我们首先用20岁去拟合，发现损失有10岁，这时我们用6岁去拟合剩下的损失，发现差距还有4岁，第三轮我们用3岁拟合剩下的差距，差距就只有一岁了。如果我们的迭代轮数还没有完，可以继续迭代下面，每一轮迭代，拟合的岁数误差都会减小。

  从上面的例子看这个思想还是蛮简单的，但是有个问题是这个损失的拟合不好度量，损失函数各种各样，怎么找到一种通用的拟合方法呢？

2 GBDT的负梯度拟合

  在上一节中，我们介绍了GBDT的基本思路，但是没有解决损失函数拟合方法的问题。针对这个问题，大牛Freidman提出了用损失函数的负梯度来拟合本轮损失的近似值，进而拟合一个CART回归树。第$m$轮的第$i$个样本的损失函数的负梯度表示为
$$r_{im}=-{\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\right]}_{f(x)=f_{m-1}(x)}$$

  利用 $(x_i,r_{im}),(i=1,2,..N)$，我们可以拟合一颗CART回归树，得到了第 $m$ 颗回归树，其对应的叶节点区域 $R_{jm},j=1,2,...,J_m$。其中$J_m$ 为叶子节点的个数。

  针对每一个叶子节点里的样本，我们求出使损失函数最小，也就是拟合叶子节点最好的的输出值${\gamma }_{jm}$如下:
$${\gamma }_{jm}=\underset{\gamma }{\mathrm{argmin}}\sum _{x_i\in R_{jm}}L(y_i,f_{m-1}(x_i)+\gamma )$$

  这样我们就得到了本轮的决策树拟合函数如下：
$$h_m(x)=\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$

  从而本轮最终得到的强学习器的表达式如下：
$$f_m(x)=f_{m-1}(x)+\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$

  通过损失函数的负梯度来拟合，我们找到了一种通用的拟合损失误差的办法，这样无轮是分类问题还是回归问题，我们通过其损失函数的负梯度的拟合，就可以用GBDT来解决我们的分类回归问题。区别仅仅在于损失函数不同导致的负梯度不同而已。

3 GBDT回归算法

  好了，有了上面的思路，下面我们总结下GBDT的回归算法。
  输入：训练集样本{ $(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$}， 最大迭代次数$M$, 损失函数$L$。
  输出：强学习器 $\hat{f}(x)$

  1) 初始化弱学习器
$$f_0(x)=\underset{\gamma }{\mathrm{argmin}}\sum _{i=1}^{N}L(y_i,\gamma )$$

  2) 对迭代轮数 $m=1,2,...M$ 有：
    a) 对样本 $i=1,2，...N$，计算负梯度
$$r_{im}=-{\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\right]}_{f(x)=f_{m-1}(x)}$$

    b) 利用 $(x_i,r_{im}),(i=1,2,..N)$，我们可以拟合一颗CART回归树，得到了第 $m$ 颗回归树，其对应的叶节点区域 $R_{jm},j=1,2,...,J_m$。其中 $J_m$ 为叶子节点的个数。

    c) 对叶子区域 $j=1,2,..J_m$，计算最佳拟合值
$${\gamma }_{jm}=\underset{\gamma }{\mathrm{argmin}}\sum _{x_i\in R_{jm}}L(y_i,f_{m-1}(x_i)+\gamma )$$
    d) 更新强学习器
$$f_m(x)=f_{m-1}(x)+\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$

  3) 得到强学习器f(x)的表达式
$$\hat{f}(x)=f_M(x)=f_0(x)+\sum _{m=1}^M\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$

4 Demo

• 对于公式 $$r_{im}=-{\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\right]}_{f(x)=f_{m-1}(x)}$$
  当 loss function $L(y_i,f(x_i))=\frac{1}{2}\ast (y_i-f(x_i){)}^2$，也即是 Least-square 时，$\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}=y_i-f(x_i)$，代入当前模型$f(x)=f_{m-1}(x)$，则有 $$\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}=y_i-f_{m-1}(x_i)$$所以我们能看到，当损失函数选用 Least-square 时，每一次拟合的值就是（真实值-当前模型的值）。
  【其它 loss function 的形式请参考 GBDT原理与Sklearn源码分析-回归篇 （★★★★★）】

• 对于初始化弱分类器 $f_0(x)$
  $$f_0(x)=\underset{\gamma }{\mathrm{argmin}}\sum _{i=1}^{N}L(y_i,\gamma )$$
  为什么需要初始化？很简单，因为每次在计算负梯度值时需要用到前一个模型 $f_{m-1}(x_i)$ 预测的值。对于我们训练的第一个模型 $m=1$ 而言需要有 $f_0(x_i)$ 的存在。那么 $f_0(x)$ 初始化为多少？这个取决于loss function的选择，当 loss function 选择 MSE 时： $$f_0(x)=\bar{y}$$
  【其它 loss function 的形式请参考 GBDT原理与Sklearn源码分析-回归篇 （★★★★★）】

• 对叶子区域 $j=1,2,..J_m$，计算最佳拟合值
  $${\gamma }_{jm}=\underset{\gamma }{\mathrm{argmin}}\sum _{x_i\in R_{jm}}L(y_i,f_{m-1}(x_i)+\gamma )$$
  那么叶子节点的取值为多少？也就是这颗树到底输出多少？ 在Friedman的论文中有这部分的推导。这里简单总结一下，叶子节点的取值和所选择的loss function有关。对于不同的Loss function，叶子节点的值也不一样。比如选择 MSE 作为 loss function 时：
  $${\gamma }_{jm}=av{e}_{x_i\in R_{jm}}\left(y_i-f_{m-1}\left(x_i\right)\right)$$
  【其它 loss function 的形式请参考 GBDT原理与Sklearn源码分析-回归篇 （★★★★★）】

掌握了以上技巧，我们下面来试试一个小例子

$x_i$	1	2	3	4	5	6	7	8	9	10
$y_i$	5.56	5.7	5.91	6.4	6.8	7.05	8.9	8.7	9	9.05

1）选择 MSE 做为建树的分裂准则
2）选择 MSE 作为误差函数
3）树的深度设置为1

• 初始化弱学习器 $$f_0(x)=\bar{y}=7.307$$

• 拟合第一颗树（$m=1$）$$r_{im}=-{\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\right]}_{f(x)=f_{m-1}(x)}=y_i-f_{m-1}(x_i)=y_i-f_0(x_i)$$

• 利用 $(x_i,r_{im}),(i=1,2,..N)$，我们可以拟合一颗CART回归树。

$x_i$	1	2	3	4	5	6	7	8	9	10
$r_{i1}$	-1.747	-1.607	-1.397	-0.907	-0.507	-0.257	1.593	1.393	1.693	1.743

这里简单介绍一下决策树建树的过程：

决策树学习最关键的步骤就是选择最优划分属性，一般而言，随着划分不过程不断的进行，我们希望决策树的分支节点所包含的样本尽可能属于同一类别（方差小）。通常，我们会选择一个准则来评价划分的质量，比如回归树中经常使用的 MSE（这种方法属于启发式的）

对于连续值，我们可以穷尽每个值 $v$，把每个值 $v$ 作为一个分裂点（$\le v$ 和 $>v$），然后计算两个分支的$MS{E}_{left}$、$MS{E}_{right}$，选择最小的$MS{E}_{sum}=MS{E}_{left}+MS{E}_{right}$ 的分裂点 $v$
选不同 $x_i$ （$x_1-x_9$）做为分裂点的结果如下：

[1.7470098765432096, 1.5140984375000002, 1.2069133786848074, 
0.99636875000000047, 0.78226400000000029, 0.32765763888888899, 
1.1579092970521545, 1.4673937500000001, 1.748733333333333] 

可以得到当选择 6 作为分裂点时 $MS{E}_{sum}=0.3276$ 最小

• 得到了第 $m$ 颗回归树，其对应的叶节点区域 $R_{jm},j=1,2,...,J_m$。其中 $J_m$ 为叶子节点的个数
  落地为：得到了第 $1$ 颗回归树，其对应的叶节点区域 $R_{j1},j=1,2$。
  $R_{11}$ 为 $x_1-x_6$
  $R_{21}$ 为 $x_7-x_{10}$
  $${\gamma }_{11}=\frac{\overline{r_{11}}+\overline{r_{21}}+\overline{r_{31}}+\overline{r_{41}}+\overline{r_{51}}+\overline{r_{61}}}{6}=-1.0703$$
  $${\gamma }_{21}=\frac{\overline{r_{71}}+\overline{r_{81}}+\overline{r_{91}}}{6}=1.6055$$

• 更新强学习器
  $$f_m(x)=f_{m-1}(x)+\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$
  落地为，更新 $f_1(x_i)$
  $$f_1(x_i)=f_0(x_i)+\sum _{j=1}^2{\gamma }_{j1}I(x_i\in R_{j1})$$
  例如更新 $x_1$ 的预测值，$x_1$ 落在 $R_{11}$，对应 ${\gamma }_{11}$，$f_1(x_1)=f_0(x_1)+{\gamma }_{11}=7.307-1.0703=6.2367$
  在GBDT里，我们通常不会直接把上一个轮的预测值 $f_{m-1}(x)$ 直接加上 ${\sum }_{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$，而是会在 ${\sum }_{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$ 上乘上一个学习率，可以理解，因为如果每次完全加上（学习率为1）本轮模型的预测值容易导致过拟合。所以通常在GBDT中的做法（也叫Shrinkage)是：$$f_m(x)=f_{m-1}(x)+\eta \sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$
  $\eta$ 为学习率，所以，当$\eta =0.1$时，上面的计算结果变为
  $$f_1(x_1)=f_0(x_1)+0.1\ast {\gamma }_{11}=7.307-0.10703=7.19997$$
  
  所有样本的更新结果如下，也即是强学习 $f_1(x)$：

[7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 
7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 
7.4675500000000001, 7.4675500000000001, 7.4675500000000001, 
7.4675500000000001] 

至此一轮迭代（第一个颗树拟合）完成，下面开始第二轮迭代（第二颗树拟合）。

---

• 强学习 $f_1(x)$

[7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 
7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 
7.4675500000000001, 7.4675500000000001, 7.4675500000000001, 
7.4675500000000001] 

• 拟合第二颗树（$m=2$）$$r_{im}=-{\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\right]}_{f(x)=f_{m-1}(x)}=y_i-f_{m-1}(x_i)=y_i-f_1(x_i)$$

• 利用 $(x_i,r_{im}),(i=1,2,..N)$，我们可以拟合一颗CART回归树。

$x_i$	1	2	3	4	5	6	7	8	9	10
$r_{i2}$	-1.63996667	-1.49996667	-1.28996667	-0.79996667	-0.39996667	-0.14996667	1.43245	1.23245	1.53245	1.58245

 [-1.639966666666667, -1.4999666666666664, -1.2899666666666665, 
 -0.79996666666666627, -0.3999666666666668, -0.1499666666666668, 
 1.4324500000000002, 1.2324499999999992, 1.5324499999999999, 
 1.5824500000000006] 

• 决策树的构建
  选不同 $x_i$ （$x_1-x_9$）做为分裂点的结果如下：

1.4289876354595334, 1.2145779184027778, 0.94102838038548786, 
0.7759147700617286, 0.63424046777777798, 0.32765763888888905, 
0.9921468829365081, 1.2236044283854166, 1.451869445987654

可以得到当选择 6 作为分裂点时 $MS{E}_{sum}$ 最小

• 得到了第 $m$ 颗回归树，其对应的叶节点区域 $R_{jm},j=1,2,...,J_m$。其中 $J_m$ 为叶子节点的个数
  落地为：得到了第 $2$ 颗回归树，其对应的叶节点区域 $R_{j2},j=1,2$。
  $R_{12}$ 为 $x_1-x_6$
  $R_{22}$ 为 $x_7-x_{10}$
  $${\gamma }_{12}=\frac{\overline{r_{12}}+\overline{r_{22}}+\overline{r_{32}}+\overline{r_{42}}+\overline{r_{52}}+\overline{r_{62}}}{6}=-0.9633$$
  $${\gamma }_{22}=\frac{\overline{r_{72}}+\overline{r_{82}}+\overline{r_{92}}}{6}=1.44495$$

• 更新强学习器
  $$f_m(x)=f_{m-1}(x)+\sum _{j=1}^{J_m}{\gamma }_{jm}I(x\in R_{jm})$$
  落地为，更新 $f_2(x_i)$
  $$f_2(x_i)=f_1(x_i)+\sum _{j=1}^2{\gamma }_{j2}I(x_i\in R_{j2})$$
  例如更新 $x_1$ 的预测值，$x_1$ 落在 $R_{12}$
  $$f_2(x_1)=f_1(x_1)+0.1\ast {\gamma }_{12}=7.19996+0.1\ast (-0.9633)=7.10363$$
  
  所有样本的更新结果如下，也即是强学习 $f_2(x)$：

[7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 
7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 
7.6120450000000002, 7.6120450000000002, 7.6120450000000002, 
7.6120450000000002] 

当只有两颗树的时候，$f_2(x)$即为预测的结果。
如果多颗树的话，拟合残差的过程如下：

---

计算上列过程的代码如下

import numpy as np
# 用 MSE 算最优分裂点
def MSE(residual):
    mse = []
    for i in range(1,len(residual)):
        mse.append(np.var(residual[:i])+ np.var(residual[i:]))
    print('mse:',mse,'\n')
    split = mse.index(min(mse))+1
    print('split:',split)
    r1 = np.mean(residual[:mse.index(min(mse))+1])
    r2 = np.mean(residual[mse.index(min(mse))+1:])
    return r1,r2,split

两棵树，深度为1的计算过程

y = [5.56,5.7,5.91,6.4,6.8,7.05,8.9,8.7,9,9.05]
print('y:',y,'\n')
F0 = np.average(y)
print('F0:',F0,'\n')
residual = y - F0
print("residual:",residual,'\n')
r11,r21,split = MSE(residual)
print('r11:',r11)
print('r21:',r21,'\n')

# F1 
F1 = []
for i in range(0,10):
    if i<split:
        F1.append(F0+0.1*r11) # learning rate 0.1
    else:
        F1.append(F0+0.1*r21)
print('F1:',F1,'\n')

residual2 = list(np.array(y) - np.array(F1))
print("residual2:",residual2,'\n')
r12,r22,split = MSE(residual2)
print('r12:',r12)
print('r22:',r22,'\n')

# F2
F2 = []
for i in range(0,10):
    if i<split:
        F2.append(F1[i]+0.1*r12) # learning rate 0.1
    else:
        F2.append(F1[i]+0.1*r22)
print('F2:',F2,'\n')

output

y: [5.56, 5.7, 5.91, 6.4, 6.8, 7.05, 8.9, 8.7, 9, 9.05] 

F0: 7.307 

residual: [-1.747 -1.607 -1.397 -0.907 -0.507 -0.257  1.593  1.393  1.693  1.743] 

mse: [1.7470098765432096, 1.5140984375000002, 1.2069133786848074, 0.99636875000000047, 0.78226400000000029, 0.32765763888888899, 1.1579092970521545, 1.4673937500000001, 1.748733333333333] 

split: 6
r11: -1.07033333333
r21: 1.6055 

F1: [7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 7.1999666666666666, 7.4675500000000001, 7.4675500000000001, 7.4675500000000001, 7.4675500000000001] 

residual2: [-1.639966666666667, -1.4999666666666664, -1.2899666666666665, -0.79996666666666627, -0.3999666666666668, -0.1499666666666668, 1.4324500000000002, 1.2324499999999992, 1.5324499999999999, 1.5824500000000006] 

mse: [1.4289876354595334, 1.2145779184027778, 0.94102838038548786, 0.7759147700617286, 0.63424046777777798, 0.32765763888888905, 0.9921468829365081, 1.2236044283854166, 1.451869445987654] 

split: 6
r12: -0.9633
r22: 1.44495 

F2: [7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 7.1036366666666666, 7.6120450000000002, 7.6120450000000002, 7.6120450000000002, 7.6120450000000002] 

【Scikit-learn源码分析请参考 GBDT原理与Sklearn源码分析-回归篇 （★★★★★）】

5 优缺点

1）优点：

• 可以灵活处理各种类型的数据，包括连续值和离散值。

• 在相对少的调参时间情况下，预测的准备率也可以比较高。这个是相对SVM来说的。

• 使用一些健壮的损失函数，对异常值的鲁棒性非常强。比如 Huber损失函数和Quantile损失函数。

2）缺点：

• 由于弱学习器之间存在依赖关系，难以并行训练数据。不过可以通过自采样的SGBT来达到部分并行。

【附录】GBDT和XGBoost的区别

1） GBDT是机器学习算法，XGBoost是该算法的工程实现
2） XGBoost 加入了正则项
3） GBDT只用了 cost function 的一阶导信息，XGBoost 对 cost fucntion 进行了泰勒展开，可同时使用一阶和二阶导
4） 传统 GBDT 用 CART作为基分类器，XGBoost 支持多种基分类器，比如线性分类器
5） 传统的 GBDT每次迭代使用所有的数据，XGBoost 采用了与随机森林相似的策略，支持对数据进行采样
6） 传统的 GBDT 没有设置对缺损值的处理，XGBoost 能自动的学习出缺失值的处理策略

参考

【1】 自我代码提升之梯度提升树
【2】干货｜从零开始学习Gradient Boosting算法
【3】 梯度提升树(GBDT)原理小结
【4】 30分钟学会用scikit-learn的基本回归方法（线性、决策树、SVM、KNN）和集成方法（随机森林，Adaboost和GBRT）
【5】A Kaggle Master Explains Gradient Boosting
【6】结合Scikit-learn介绍几种常用的特征选择方法
【7】scikit-learn 梯度提升树(GBDT)调参小结
【8】scikit-learn (sklearn) 官方文档中文版
【9】scikit-learn (sklearn) 官方文档
【10】GBDT算法原理以及实例理解（★★★★）
【11】GBDT原理与Sklearn源码分析-回归篇 （★★★★★）
【12】GBDT源码剖析
【13】sklearn.ensemble.GradientBoostingRegressor