GBDT算法原理及其代码实现

前言

GBDT算法，全名为梯度提升树，是以CART决策树为弱学习器并基于Boosting思想的一种集成学习算法。其与传统的Boosting算法有较大差别，但是也是一种弱学习器之间相互关联，只能串行训练的机器学习算法。GBDT在机器学习算法中具有举足轻重的地位，在深度学习流行之前，在kaggle，天池等竞赛中常常占据榜首，并且在工业界也被经常使用。本文将介绍GBDT的算法原理及其sklearn代码实现。

注：本文基于决策树与集成学习，对此不了解的可以移步本站的决策树与集成学习的博客，决策树上篇，决策树下篇，集成学习简介。另外，本文也会谈及到AdaBoost算法，想要了解AdaBoost算法可以移步本站的AdaBoost算法的博客，AdaBoost算法。

GBDT算法思想

提升思想

GBDT基于Boosting，但不完全是Boosting，相比于AdaBoost这类完全遵循Boosting框架的算法，GBDT对传统的Boosting框架做了一定的改动，原Boosting算法框架图如下

相比于传统的Boosting采用的更新训练集权重，GBDT采用的方法是直接修改训练集，每一次训练使用的数据集为上一次训练出的学习器的拟合残差。也就是说，对于之前的学习器所没有拟合好的样本，在这一轮中会被新的学习器继续缩小差距。这就是所谓的提升思想。

举例来说，假设真实值为100，第一个学习器预测值为80，那残差就是20，第二个学习器的任务就是拟合这个残差，也就是尽量靠近20，假设第二个学习器预测值为15，则残差为5，送到第三个学习器继续学习，最终将所有学习器的预测值相加即可。

综上，GBDT的流程可以总结为下图

负梯度拟合

虽然上文提到GBDT是采用不断拟合残差的训练方法，也就是说$r_{ki}$代表第$k$轮训练时使用的第$i$个样本为$(x_i,y_i - f(x_i))$。但是对于一般的损失函数而言，计算每一步优化往往不是那么简单，因此GBDT使用了损失函数的负梯度来近似残差，具体定义为

$$r_{ki}=-\Big[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\Big]_{f(x)=f_{k-1}(x)}$$

也就是说第$k$轮训练时使用的第$i$个样本为$(x_i,r_{ki})$，$L(y,f(x))$代表损失函数，$f(x) = f_{k - 1}(x)$代表前面$k - 1$轮训练出的学习器。当损失函数为平方误差时，负梯度与残差相等（严格来说是相差常数倍）

$$-\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)} = -\frac{\partial (y_i - f(x_i))^2}{\partial f(x_i)} = -2(y_i - f(x_i))(-1) = 2(y_i - f(x_i))$$

补充

关于负梯度，上文我们是直接给出，这里给出负梯度的构造过程。首先，我们来回顾一下泰勒展开公式，泰勒展开是对于一个函数的近似，其展开阶数越高，拟合程度越高，并且其在给定的$x = x_0$附近的拟合效果是最好的。这里我们回顾一下一阶的泰勒展开

$$f(x) \approx f(x_0)+f^{\prime}(x_0)(x-x_0)$$

这里我们不关注余项。将损失函数$L(y,f(x))$在$f(x) = f_{k-1}(x)$处进行一阶泰勒展开

$$L(y,f(x)) \approx L(y,f_{k-1}(x)) + \Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)} \left( f(x) - f_{k-1}(x) \right)$$

令$f(x) = f_k(x)$，且$f_k(x) = f_{k-1}(x) + h_k(x)$，得

$$L(y,f_k(x)) \approx L(y,f_{k-1}(x)) + \Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)} \left( f_k(x) - f_{k-1}(x) \right) \\ = L(y,f_{k-1}(x)) + \Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)} h_k(x)$$

移项可得

$$L(y,f_{k}(x)) - L(y,f_{k-1}(x)) \approx \Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)} h_k(x)$$

现在我们的目的是求解$h_k(x)$使得$L(y,f_k(x))$最小，由于$L(y,f_{k-1}(x))$是一个定值，等同于使得$L(y,f_{k}(x)) - L(y,f_{k-1}(x))$最小，故$\Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)} h_k(x)$可以作为目标函数的近似，$\Big[\frac{\partial L(y,f(x))}{\partial f(x)}\Big]_{f(x)=f_{k-1}(x)}$为$h_k(x)$的梯度，要使得$h_k(x)$取到最小，那么可以往梯度的负方向优化，其思想类似于梯度下降，这也就构造出了负梯度。

GBDT回归树

从上文中我们可以注意到，GBDT是从回归问题的角度出发的，因此我们从回归问题开始讲起。

GBDT回归流程

根据上文中给出的GBDT思想，我们可以给出构建GBDT回归树的流程

输入：训练集$T = \{ (x_i, y_i), i = 1,2,\cdots,m \}$，训练次数$K$

输出：强学习器$f(x)$

(1) 使用训练集$T$训练初始弱学习器$f_0(x) = h_0(x) = \min_{c} \sum_{i = 1}^{m} L(y_i,c)$，其中$h_0(x)$为只有根节点的CART回归树

(2) 对于第$k$轮训练，计算负梯度$r_{ki}=-\Big[\frac{\partial L(y_i,f(x_i)))}{\partial f(x_i)}\Big]_{f(x)=f_{k-1}(x)}$，得到训练集$T_k = \{ (x_i, r_{ki}), i = 1,2,\cdots,m \}$

(3) 使用训练集$T_k$训练得到CART回归树$h_k(x)$，更新强学习器为$f_k(x) = f_{k - 1}(x) + h_k(x)$

(4) 经过$K$轮训练后得到强学习器$f(x) = f_K(x) = \sum_{k = 0}^{K} h_k(x)$

回归常用损失函数

对于回归问题，常用的损失函数有四种，分别为

平方误差（均方差）：最常用的损失函数

$$L(y,f(x)) = (y - f(x))^2$$

绝对误差：在$x$为0处不可导，因此没有解析解，但是对常用数值计算方法求近似解的机器学习算法来说，无伤大雅

$$L(y,f(x)) = |y - f(x)|$$

Huber损失：这是平方误差和绝对损失的折中产物，对于远离中心的异常点，采用绝对损失，而中心附近的点采用均方差，这个界限一般用分位数点度量。

$$L(y,f(x))=\begin{cases}\frac12(y-f(x))^2&|y-f(x)|\leq\delta\\\delta(|y-f(x)|-\frac\delta2)&|y-f(x)|>\delta&\end{cases}$$

分位数损失：即分位数回归的损失函数

$$L(y,f(x))=\sum_{y\geq f(x)}\theta|y-f(x)|+\sum_{y<f(x)}(1-\theta)|y-f(x)|$$

其中$\theta$为分位数。

注：一般使用前两种损失函数，后两种主要用于数据异常点较多的情况，用于提高模型的健壮性。

GBDT分类树

而对于分类问题来说，由于输出$y$是离散值，这使得我们无法描述残差。GBDT使用的做法为，将分类问题也看做回归问题对待，这类似于logistic回归使用线性回归模型来处理分类问题，弱学习器依旧使用CART回归树，使用损失函数来描述残差。

二分类问题

logistic回归回顾

与逻辑回归问题相同，我们先从二分类问题开始讲起，这里我们简单回顾一下逻辑回归，其思想为通过sigmoid函数，将线性回归的输出映射到$(0,1)$代表概率，定义二分类变量$y \in \{ 0,1 \}$，则逻辑回归模型可以写为

$$f(x) = \sigma(h(x)) = \frac{1}{1 + e^{-h(x)}}$$

其中$h(x)$为线性回归模型，$f(x)$可以理解为$x$属于类别1的概率，而$x$属于类别0的概率自然就是$1 - f(x)$，记为

$$p_1 = f(x) \\ p_0 = 1 - f(x)$$

上述两个概率值可以统一写为

$$p(x) = f(x)^{y} (1 - f(x))^{1 - y} = p_1^{y} p_0^{1 - y}$$

可以注意到，当$y = 1$时，$p(x) = p_1$，即判断为第一类的概率，反之同理。因此$p(x)$的意义为判断正确的概率。因此逻辑回归问题的目标函数为

$$\prod_{i=1}^m f(x_i)^{y_i}(1-f(x_i))^{(1-y_i)}$$

也就是令所有样本的判断正确概率最大，为了方便计算，我们将其取对数，并加上负号，得到损失函数为

$$L = -\sum_{i=1}^m y_i \ln(f(x_i))+(1-y_i) \ln(1-f(x_i))$$

该损失函数被称为对数似然损失函数。

GBDT二分类算法的流程

现在我们将逻辑回归的理论套用至回归树上，只需要将$h(x)$替换为CART回归树即可，损失函数依旧为

$$L(y,f(x)) = - y \ln(f(x)) - (1-y) \ln(1-f(x))$$

注：这里写的是单个样本的损失函数，上面写的是所有样本的损失函数，因为这里我们要用单个样本的损失函数来求负梯度，实际上两个损失函数是一样的。

有了损失函数之后，就是求负梯度

$$r_{ki}=-\Big[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}\Big] = \frac{1 - y_i}{f(x_i)} - \frac{y_i}{f(x_i)} = \frac{f(x_i) - y_i}{f(x_i)(1 - f(x_i))}$$

至此，我们可以给出GBDT二分类算法的流程，这与GBDT回归树相比仅仅是换了一个损失函数

输入：训练集$T = \{ (x_i, y_i), i = 1,2,\cdots,m \}$，训练次数$K$

输出：强学习器$f(x)$

(1) 使用训练集$T$训练初始弱学习器$f_0(x) = h_0(x) = \min_{c} \sum_{i = 1}^{m} L(y_i,c)$，其中$h_0(x)$为只有根节点的CART回归树

(2) 对于第$k$轮训练，计算负梯度$r_{ki}=-\Big[\frac{\partial L(y_i,f(x_i)))}{\partial f(x_i)}\Big]_{f(x)=f_{k-1}(x)} = \frac{f(x_i) - y_i}{f(x_i)(1 - f(x_i))}$，得到训练集$T_k = \{ (x_i, r_{ki}), i = 1,2,\cdots,m \}$

(3) 使用训练集$T_k$训练得到CART回归树$h_k(x)$，更新强学习器为$f_k(x) = f_{k - 1}(x) + h_k(x)$

(4) 经过$K$轮训练后得到强学习器$f(x) = f_K(x) = \sum_{k = 0}^{K} h_k(x)$

多分类问题

softmax回归回顾

而对于多分类问题，我们可以参照一般化的logistic回归，也就是softmax回归。其使用的是one-vs-rest的策略，也就是对于$R$个类别的问题，构建$R$个逻辑回归模型，但是每一个模型的输出意义不同，第$r$个模型的输出$f_r(x)$代表样本属于第$r$类的概率，反之$1 - f_r(x)$代表样本不属于第$r$类的概率。最后将其结合为一个向量$Y$，向量$Y$的第$r$个分量$Y^{(r)}$代表该样本属于第$r$类的概率，为了满足概率的要求，需要该向量做归一化，以保证向量和为1。理想情况下的$Y$应该是只有正确类别对应的分量为1，其余分量均为0，因此可以对训练集的标签进行独热编码，将每一个样本的标签转化为$1 \times R$的向量。

softmax回归所使用的softmax函数，就是完成了从线性回归到概率向量的转化，也就是将第$r$个线性回归模型的输出$h_r(x)$直接转化为向量$Y$的第$r$个分量$Y^{(r)}$。该函数定义为

$$\mathrm{softmax}(x_i) = \frac{x_i}{\sum_{j=1}^{R} x_j}$$

也就是说，在得到$R$个$h(x)$后，可以计算得向量$Y$的各分量为

$$Y^{(r)} = \mathrm{softmax}(h_r(x)) = \frac{h_r(x)}{\sum_{j=1}^{R} h_j(x)}$$

综上，softmax回归的流程为

接着是损失函数，我们记真实值的标签向量为$Y$，这是一个仅有一个分量为1，其余分量均为0的向量，记概率向量为$P$，其每一个分量都代表样本被分到该类的概率，并且向量和为1。很自然的，我们会希望两个向量的距离越近越好，也就是希望向量的模$||Y-P||_2$越小。但是在分类问题中，我们其实并不用追求正确的分类对应的分量取值为1，如一个归属与第$r$类的样本，只要$P^{(r)} > 1/R$，其实就可以判断正确，所以平方损失函数有些过于严苛了，我们可以使用交叉熵来衡量

$$H(Y,P) = - \sum_{j = 1}^{R} Y^{(j)} \ln(P^{(j)})$$

注意到向量$Y$仅有一个分量为1，因此对属于第$r$类的样本

$$- \sum_{j = 1}^{R} Y^{(j)} \ln(P^{(j)}) = -\ln(P^{(r)})$$

也就是说交叉熵仅关注正确类别对应的概率，这很符合我们的需求。

GBDT多分类算法的流程

GBDT多分类算法就是在softmax回归的基础上，将线性回归模型替换为CART回归树，也就是说对于有$R$个类别的分类问题，我们需要构建$R$颗CART回归树，也就是说，我们需要将GBDT二分类算法做的事情做$R$次。那就剩最后一个问题，就是确定损失函数与负梯度。接下来我们来看对于第$r$颗树的训练，损失函数使用softmax回归中使用的交叉熵损失函数

$$H(Y,P) = - \sum_{j = 1}^{R} Y^{(j)} \ln(P^{(j)})$$

由于$P^{(j)} = \mathrm{softmax}(f^{(j)}(x))$，故损失函数为

$$L(Y^{(j)},f^{(j)}(x)) = - \sum_{j = 1}^{R} Y^{(j)} \ln \left( \frac{e^{f^{(j)}(x)}}{\sum_{i=1}^{R} f^{(i)}(x)} \right) = - \sum_{j = 1}^{R} Y^{(j)} \left[ f^{(j)}(x) - \ln \left( \sum_{i=1}^{R} f^{(i)}(x) \right) \right]$$

由于$\sum_{j = 1}^{R} Y^{(j)} = 1$，故

$$L(Y^{(j)},f^{(j)}(x)) = \ln \left( \sum_{i=1}^{R} f^{(i)}(x) \right) - \sum_{j = 1}^{R} Y^{(j)} f^{(j)}(x)$$

最后是求解负梯度，这里需要注意的是，当我们求解第一个GBDT分类器$f^{(1)}(x)$时，我们是将损失函数对$f^{(1)}(x)$求偏导，也就是

$$-\Big[\frac{\partial L(Y^{(j)},f^{(j)}(x))}{\partial f^{(1)}(x)}\Big] = Y^{(1)} - \frac{e^{f^{(1)}(x)}}{\sum_{i=1}^{R} f^{(i)}(x)} = Y^{(1)} - \mathrm{softmax}(f^{(1)}(x))$$

因此，在我们训练第一个GBDT分类器时，第$k$轮训练时使用的第$i$个样本为$(x_i,r_{ki})$，其中

$$r_{ki} = -\Big[\frac{\partial L(Y^{(j)},f^{(j)}(x))}{\partial f^{(1)}(x)}\Big]_{f^{(1)}(x)=f^{(1)}_{k-1}(x)} = Y^{(1)} - \mathrm{softmax}(f^{(1)}(x))$$

也就是说，在我们训练第$r$个GBDT分类器时，第$k$轮训练时使用的第$i$个样本为$(x_i,r_{ki})$，其中

$$r_{ki} = -\Big[\frac{\partial L(Y^{(j)},f^{(j)}(x))}{\partial f^{(r)}(x)}\Big]_{f^{(r)}(x)=f^{(r)}_{k-1}(x)} = Y^{(r)} - \mathrm{softmax}(f^{(r)}(x))$$

至此，我们可以得到GBDT多分类算法的流程为

输入：训练集$T = \{ (x_i, y_i), i = 1,2,\cdots,m \}$，训练次数$K$，样本类别数$R$

输出：强学习器$F(x)$

(1) 对训练集$T$进行独热编码，得到$R$个训练集$T^{(r)},r = 1,2,\cdots,R$

(2) 对于第$r$个学习器的训练，使用训练集$T^{(r)}$训练初始弱学习器$f_0^{(r)}(x) = h_0(x) = \min_{c} \sum_{i = 1}^{m} L(y_i,c)$，其中$h_0(x)$为只有根节点的CART回归树

(3) 对于第$k$轮训练，计算负梯度$r_{ki} = -\Big[\frac{\partial L(Y^{(j)},f^{(j)}(x))}{\partial f^{(r)}(x)}\Big]_{f^{(r)}(x)=f^{(r)}_{k-1}(x)}$，得到训练集$T_k = \{ (x_i, r_{ki}), i = 1,2,\cdots,m \}$

(4) 使用训练集$T_k$训练得到CART回归树$h_k(x)$，更新强学习器为$f_k^{(r)}(x) = f_{k - 1}^{(r)}(x) + h_k(x)$

(5) 对于第$r$个学习器的训练，经过$K$轮训练后得到学习器$f^{(r)}(x) = f_K^{(r)}(x) = \sum_{k = 0}^{K} h_k(x)$

(6) 经过$R$次GBDT学习器的构建，将$R$个GBDT学习器结合为一个向量$f(x) = (f^{(1)}(x),f^{(2)}(x),\cdots,f^{(R)}(x))$

(7) 最终输出为$f(x)$中取值最大的分量对应的类别，即$F(x) = max_r f^{(r)}(x)$

也可以总结为下图

sklearn代码实现

最后我们来介绍一下sklearn代码实现，运行环境为jupyter notebook，需要用到的库为

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

from sklearn.ensemble import AdaBoostClassifier  # AdaBoost分类
from sklearn.ensemble import AdaBoostRegressor  # AdaBoost回归
from sklearn.tree import DecisionTreeClassifier  # 决策树分类
from sklearn.tree import DecisionTreeRegressor  # 决策树回归

from sklearn.datasets import load_wine  # 红酒数据集
from sklearn.datasets import load_boston  # 波士顿房价数据集

from sklearn.model_selection import train_test_split  # 切分数据

from scipy import stats

import warnings
warnings.filterwarnings("ignore")

使用红酒数据集与波士顿房价数据集作为示例数据

# 读取数据
wine = load_wine()
boston = load_boston()

参数列表

在sklearn中，实现GBDT的函数为GradientBoostingClassifier()与GradientBoostingRegressor()，其常用参数列表为

参数	说明
n_estimators	训练次数，也就是弱学习器数量，在任何集成学习算法中都是重要的参数。
learning_rate	学习率，默认为1，较大的学习率可以完成模型的快速收敛，较小的学习率可以使得模型的学习更加的精细。learning_rate和n_estimators之间存在权衡关系，一般较小的学习率要配合较大的训练次数。
loss	损失函数。对于分类问题，可填{‘deviance’, ‘exponential’}，默认为’deviance’，’deviance’指对数似然损失函数，‘exponential’指指数损失函数，也就是AdaBoost中所使用的。对于回归问题，可填{‘ls’, ‘lad’, ‘huber’, ‘quantile’}，默认为‘ls’，其中‘ls’指平方损失函数，‘lad’指绝对损失函数，‘huber’指Huber损失函数，‘quantile’指分位数损失函数。

由于GBDT使用决策树作为其弱学习器，与随机森林一样，其继承了决策树的参数，常用的有

max_depth：最大深度

min_samples_leaf & min_samples_split：节点或子节点需要包含的最少训练样本

min_impurity_decrease：限制信息增益

max_features：限制特征数

random_state：随机数种子，用于随机选择特征，保证每次训练得到的结果一致。

下面我们已经以红酒数据集与波士顿房价数据集作为示例数据集，来跑分类与回归的demo

分类

# 提取特征与标签
data = pd.DataFrame(wine.data)
data.columns = wine.feature_names
target = pd.DataFrame(wine.target)
target.columns = ['class']
# 划分训练集与测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(data, target, test_size=0.2)

# GBDT
clf = GradientBoostingClassifier(n_estimators=200
                                ,learning_rate=0.5
                                ,max_depth=4
                                ,random_state=1008)
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)

作为以决策树为弱学习器的集成学习方法，GBDT自然也继承了决策树的一些重要属性，如特征重要性

clf.feature_importances_  # 特征重要性

回归

# 提取特征与标签
data = pd.DataFrame(boston.data)
data.columns = boston.feature_names
target = pd.DataFrame(boston.target)
target.columns = ['MEDV']
# 划分训练集与测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(data, target, test_size=0.2)

# GBDT
clf = GradientBoostingRegressor(n_estimators=200
                               ,learning_rate=0.5
                               ,max_depth=3
                               ,min_samples_split=40
                               ,min_samples_leaf=40
                               ,random_state=1008)
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)