极端提升树XGBoost及其实践

XGBoost stands for eXtreme Gradient Boosting.
XGBoost is an implementation of gradient boosted decision trees designed for speed and performance.

Speed: 并行计算每个样本的一、二次梯度，CART分裂时，也可并行计算不同分裂点，从中选择最优分裂点。

Performance: 修剪技术、正则项、考虑二次梯度。

目标函数：训练损失 + 正则项

目标函数由两部分组成：训练损失$L(\theta )$和正则项$\Omega (\theta )$。
$$\text{obj}(\theta )=L(\theta )+\Omega (\theta )$$
训练损失用于测量模型在训练集的预测性能，训练损失可以是 mean squared error:
$$L(\theta )=\sum _i(y_i-{\hat{y}}_i{)}^2$$
或者是 logistic loss:
$$L(\theta )=\sum _i[y_i\ln (1+e^{-{\hat{y}}_i})+(1-y_i)\ln (1+e^{{\hat{y}}_i})]$$
正则项控制模型复杂度，以避免overfitting，考虑下列图片中对同一数据的不同拟合：

The model marked in red visually seems a reasonable fit to you (bias-variance tradeoff).
Understanding the process in a formalized way also helps us to understand the objective that we are learning and the reason behind the heuristics such as pruning and smoothing

集成决策树

以下是使用CART决策树预测是否喜欢电脑游戏的例子：

家族成员被分配到不同叶子节点，CART中不同叶节点对应不同分数，使用分数可以使用统一的优化方法（后续章节介绍）。通常情况下，单棵树在实践中效果不好，可利用集成思想，即建立多棵树，并把这些树的结果组合起来作为最终输出：

以上集成两棵树，同一样本在两棵树中对应叶节点的分数值相加作为最终输出。从例中可看到，两棵树互补，从数学角度可解释为
$${\hat{y}}_i=\sum _{k=1}^K{f}_k(x_i),f_k\in \mathcal{F}$$

其中$K$是树的数量，$f$是函数空间$\mathcal{F}$中的函数，$\mathcal{F}$是所有可能CARTs的集合。

目标函数的优化由下列公式给出：
$$\text{obj}(\theta )=\sum _i^n\ell (y_i,{\hat{y}}_i)+\sum _{k=1}^K\Omega (f_k)$$

XGBoost 训练

模型参数是什么？ 从目标函数可以看出，模型参数是多棵叶节点输出分数值的决策树。

树优化不能使用传统的欧式空间中的优化方法，因此我们采用累加策略：迭代学习每一棵子树用于修正当前损失（残差学习），并把所学子树组合到总体。若样本$(x_i,y_i)$在时间步$t$的预测为${\hat{y}}_i$，则模型学习过程（逐步逼近）为
$$\begin{array}{rl}{\hat{y}}_i^{(0)}& =0\\ {\hat{y}}_i^{(1)}& =f_1(x_i)={\hat{y}}_i^{(0)}+f_1(x_i)\\ {\hat{y}}_i^{(2)}& =f_1(x_i)+f_2(x_i)={\hat{y}}_i^{(1)}+f_2(x_i)\\ & \dots \\ {\hat{y}}_i^{(t)}& =\sum _{k=1}^t{f}_k(x_i)={\hat{y}}_i^{(t-1)}+f_t(x_i)\end{array}$$
子树如何学习？ 新的子树可降低总体损失，$t$时刻，我们定义目标函数为：
$$\begin{array}{rl}{\text{obj}}^{(t)}& =\sum _{i=1}^n\ell (y_i,{\hat{y}}_i^{(t)})+\sum _{i=1}^t\Omega (f_i)\\ & =\sum _{i=1}^n\ell (y_i,{\hat{y}}_i^{(t-1)}+f_t(x_i))+\Omega (f_t)+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\end{array}$$

我们将$f_t(x)$作为增量，使用泰勒公式将损失函数$\ell$在${\hat{y}}^{(t-1)}$处二次展开得
$$\ell (y,{\hat{y}}^{(t-1)}+f_t(x))=\ell (y,{\hat{y}}^{(t-1)})+g_t{f}_t(x)+\frac{1}{2}{h}_t{f}_t^2(x)$$
式中一阶导$g_t=\partial \ell /\partial {\hat{y}}^{(t-1)}$，二阶导$h_t=\partial {\ell }^2/{\partial }^2{\hat{y}}^{(t-1)}$，$\ell (y,{\hat{y}}^{(t-1)})$为常数。

此时可以得到近似目标函数
$${\text{obj}}^{(t)}\approx \sum _{i=1}^n[g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}_t^2(x_i)]+\Omega (f_t)+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}$$
最小化目标函数求解第$t$步子树，第$t$步损失函数仅包含之前损失对样本上的一、二阶导，XGBoost支持自定义损失函数。

Taylor Expression
$$f(x+\Delta x)=f(x)+f^{\prime }(x)\Delta x+\frac{1}{2}{f}^{\prime \prime }(x)(\Delta x{)}^2$$

令$q(x)$表示样本$x$所属叶节点，$w$表示所有叶节点对应的分数向量，$T$表示叶节点数量，则决策树$f_t(x)$表示为
f t ( x ) = w q ( x ) w ∈ R T , q : R d → { 1 , 2 , ⋯ , T } . \color{darkcyan} f_t(x)=w_{q(x)} \quad w∈R^T,q:R_d→\{1,2,⋯,T\}. ft​(x)=wq(x)​w∈RT,q:Rd​→{1,2,⋯,T}.
XGBoost的复杂度定义为
$$\Omega (f)=\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2$$
定义树复杂度的方式有很多，但是这种方式的定义在实践中表现很好。

现在我们将损失函数中对样本损失求和转化为对叶节点损失求和，我们这样做的原因是因为同一叶节点中的样本分数相同，且样本数一定小于节点数，这样可以简化计算量或支持并行计算，因此，目标函数在时间$t$的损失为
$$\begin{array}{rl}{\text{obj}}^{(t)}& \approx \sum _{i=1}^n[g_i{w}_{q(x_i)}+\frac{1}{2}{h}_i{w}_{q(x_i)}^2]+\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2\\ & =\sum _{j=1}^T\left[w_j{\sum }_{i\in I_j}{g}_i+\frac{1}{2}{w}_j^2\left(\lambda +{\sum }_{i\in I_j}{h}_i\right)\right]+\gamma T\end{array}$$
式中 I j = { i ∣ q ( x i ) = j } \color{darkcyan}I_j = \{i|q(x_i)=j\} Ij​={i∣q(xi​)=j}是节点$j$中的样本集，令$G_j={\sum }_{i\in I_j}{g}_i,H_j={\sum }_{i\in I_j}{h}_i$，目标函数简化为
$${\text{obj}}^{(t)}=\sum _{j=1}^T[G_j{w}_j+\frac{1}{2}(H_j+\lambda )w_j^2]+\gamma T$$
其中$w_j$相互独立，在特定树结构$q(x)$中，最优叶节点输出值和目标函数为
$$\begin{array}{rl}{w}_j^{\ast }& =-\frac{G_j}{H_j+\lambda }\\ {\text{obj}}^{\ast }& =-\frac{1}{2}\sum _{j=1}^T\frac{G_j^2}{H_j+\lambda }+\gamma T\end{array}$$
剩下最后一个问题：怎么学习子树的结构$q(x)$？

学习到的子树应能极小化目标函数，可枚举所有可能的子树，并选择使损失最小的子树，但计算量太大。因此，我们通过贪婪方式学习子树：每个节点的分裂均能极小化目标函数，直到达到分类终止条件，如达最大深度、或分裂之后损失不在下降。

我们首先可获得每个样本的$g_i$和$h_i$，在特定的树结构$q(x)$下，可以很容易计算出每个节点的$G_j$和$H_j$。根据单节点二分裂前后对应的最优目标函数$\text{obj}$，可以计算出分裂后降低的损失（增益） 为
$$Gain=\frac{1}{2}\left[\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 }\right]-\gamma$$
公式中四项的解释：划分后左子节点的损失、划分后右子节点的损失、原始节点的损失、叶节点上的正则项系数。
重要特性：额外增加的$\gamma$用于控制复杂度，若增益小于$\gamma$，则停止分裂，也称为决策树 “pruning (修剪)技术”。

如何选择分裂特征和分裂点？对于实值特征，将所有样本按此特征值顺序排列，从左到右扫描，可有效地找到最优分裂点：

对于N个样本M个实值特征（特征取值不重复）的训练集，需遍历M*(N-1)次，此外当数据无法全载入内存时，也不能使用这种方法。

由于模型整体是在不断迭代极小化损失，我们在学习子树时，也不一定需要找到最优子树，因此分裂点可采用等分位点方法，对于大数据集可显著降低计算量，如果子树较多，效果并不一定差。

缩减和随机特征采样

这里介绍两种额外的避免过拟合的技术。

第一种技术是Friedman提出的 “shrinkage (缩减)”，shrinkage是将所学新树乘以一个参数$\eta$之后再加到总体模型，shrinkage降低了单棵树对总体的影响，给出来子树学习提供了更多的空间。

第二种技术是 “column (feature) subsampling (随机特征)”，这种技术常用于随机森林（子树是在随机样本、随机特征中学习），column sub-sampling 比 “row sub-sampling (随机样本)” 更能防止过拟合。

XGBoost 糖尿病预测

# First XGBoost model for Pima Indians dataset
import numpy as np
import pandas as pd
from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from xgboost import XGBClassifier

# 数据集下载地址：https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv
df = pd.read_csv('pima-indians-diabetes.data.txt', header=None)
df.columns = ['preg_times', 'glucose', 'blood_pressure', 'skin_thick', 'insulin', 'BMI',
              'DPF', 'age', 'outcome']

print(df.describe())
print(df.shape)
print(df.head())

# glucose、blood_pressure、skin_thick、insulin、BMI 0值为缺失值，填充为均值
missing_columns = ['glucose', 'blood_pressure', 'skin_thick', 'insulin', 'BMI']
df[missing_columns] = df[missing_columns].replace(0, np.nan)
for column in missing_columns:
    df.loc[df[column].isna(), column] = df[column].mean()

print('missing filled.')
print(df.describe())

# load data
# dataset = loadtxt('pima-indians-diabetes.data.txt', delimiter=",")
dataset = df.values

# split data into X and y
X = dataset[:, 0:8]
y = dataset[:, 8]

# split data into train and test sets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=7, shuffle=True)

model = XGBClassifier(
    eta=0.08,  # 别名learning_rate，子树shrinkage
    gamma=0.001,  # 别名min_split_loss，分割损失小于改值时停止分割
    max_depth=3,  # 最大深度
    min_child_weight=2,  # 子节点样本权重和低于该值时停止分割
    max_delta_step=0,
    subsample=0.8,  # 每次训练时样本采样率
    sampling_method='uniform',  # 训练样本采样方法
    colsample_bytree=0.7,  # 每棵树的特征采样比
    colsample_bylevel=0.8,  # 每层的特征采样比
    colsample_bynode=0.7,  # 每个节点的特征采样比，对防止过拟合非常重要
    reg_alpha=0,  # L1正则化系数
    reg_lambda=1,  # L2正则化系数
    tree_method='auto',  # 树构造算法
    base_score=0.5,
    booster='gbtree',
    missing=None,
    n_estimators=80,
    n_jobs=4,
    objective='binary:logistic',
    random_state=2,
    seed=6,
    scale_pos_weight=1.0,
    silent=True,
    verbosity=2
)
model.fit(X_train, y_train)

# cv_result = cross_val_score(model, X, y, cv=KFold(
#     n_splits=10, random_state=22, shuffle=True))
# print('Cross validation accuracy: %.2f%%' % (cv_result.mean() * 100))

# train predictions
acc_train = accuracy_score(y_train, model.predict(X_train))
print("Train accuracy: %.2f%%" % (acc_train * 100.0))

# evaluate predictions
acc_test = accuracy_score(y_test, model.predict(X_test))
print("Test accuracy: %.2f%%" % (acc_test * 100.0))

预测结果：

Train accuracy: 84.20%
Test accuracy: 83.12%

Reference

1. XGBoost Tutorials
2. XGBoost: A Scalable Tree Boosting System