从CART到Xgboost - 原理及代码实践

CART (Classification and Regression Truee)

CART,即分类回归树。顾名思义，CART既可以用于分类任务，又可以用于回归任务.决策树的建模思想非常简单，将特征空间划分位互不相交的子空间，也就是将样本进行分组；然后位每一组样本分配一个预测值。在预测阶段，对一个新来的样本，位样本找到分组，然后给出预测值。

• 分类任务：预测目标是离散值，分类树采用Gini值、entropy、information gain作为节点分裂的依据。采用叶子节点里概率最大的类别作为当前节点的预测类别
• 回归任务： 预测目标是连续值，采用均方误差MSE作为节点分类的依据。则采用叶子节点的均值作为预测值

分类树和回归树的生成步骤主要有三步：

1. 特征选择
2. 树的生成
3. 树的剪枝

分类树

样本集合$D$的基尼指数 (值越大，不确定性越大):

$$Gini(D)=1-\sum _{k=1}^K(\frac{|C_k|}{|D|}{)}^2$$

其中, $C_k$是属于第$k$类的子集.
假定现在将样本集合$D$按照最优切分特征$x_j$进行划分，最优切分特征值为$s$，划分为$x_j>s$和$x_j<s$两部分$D_1$, $D_2$。计算划分之后的基尼系数:
$$Gini(D,x_j)=\frac{|D_1|}{|D|}Gini(D_1)+\frac{D_2}{D}Gini(D_2)$$
对于每一个结点，遍历所有的特征和特征值，找到最小的$Gini(D,A)$. 然后分治递归知道满足终止条件 （例如误差小于某阈值、样本个数小于某阈值、没有更多特征等）

回归树

回归树的预测值:
$$f(x)=\sum _{m=1}^{M}I(x\in R_m)$$
即输入空间划分为$M$个单元，$R_1,...,R_M$, 单元$R_M$对应的输出为$c_m$，$I$是指示函数. 输出$c_m$取$I(x\in R_m)$输出的均值 ：
$$c_m=ave(y_i|x_i\in R_m)$$

假定现在将样本集合$D$按照最优切分特征$x_j$进行划分，最优切分特征值为$s$, 输入空间被划分为两部分:

R a ( j , s ) = { x ∣ x j ≤ s } R^a(j, s) = \{x|x^{j} \leq s \} Ra(j,s)={x∣xj≤s} 和 R b ( j , s ) = { x ∣ x j > s } R^b(j, s) = \{x|x^{j} > s \} Rb(j,s)={x∣xj>s}

每个被划分的空间的输出是$c^a$和$c^b$。因此我们的优化目标是:
$$\underset{j,s}{\min }[\sum _{x_i\in R^a(j,s)}(y_i-c^a{)}^2+\sum _{x_i\in R_b(j,s)}(y_i-c^b{)}^2]$$

遍历所有的特征$x_j$和特征切分点$s$。找到满足上式最优的$x_j,s$。之后分治递归直到满足终止条件 (例如误差小于某阈值、样本个数小于某阈值、没有更多特征等)

树的剪枝

首先在生成完整的决策树之后，获取该决策树的所有子树序列。对于任意的一个子树都可以用如下方式计算其代价函数：
$$J_{\lambda }(T)=J(T)+\lambda |T|$$

式中的$T$就是目标子树，$J(T)$是对训练数据的预测误差（平方误差或者基尼系数）， $|T|$为子树的叶结点个数， $\lambda$ 则是正则项系数。利用代价函数进行剪枝的具体过程如下：

1. 首先对完整的生成树$T_0$进行剪枝，对其中任意一个内部结点$t$都有：
   • 以$t$为根节点的子树(不剪枝)的代价函数为$J_{\lambda }(T_t)=J(T_t)+\lambda |T_t|$
   • 以$t$为单节点的子树(不剪枝)的代价函数为$J_{\lambda }(t)=J(T)+\lambda$
2. 随着$\lambda$不断增大，一定有$J_{\lambda }(T_t)=J_{\lambda }(t)$, 此时$\lambda =\frac{J(t)-J(T_t)}{|T_t|-1}$. 两种情况下代价函数值相等，而$t$比$T_t$节点要少，复杂度更低，相应的泛化能力越强，所以此时对$T_t$进行剪枝
3. 由上，对完全生成树$T_0$ 的每一个中间结点都计算出这个$\lambda$ 值，选取该值最小的结点进行剪枝，得到一个子树$T_1$ ，并将该 值设为 ${\lambda }_1$ 。对该子树递归的使用上述步骤，直至剪枝至只剩根结点。此时，就有了子树序列 { T 0 , T 1 , . . . , T n } \{ T_0, T_1,...,T_n \} {T0​,T1​,...,Tn​} 以及参数序列 { λ 0 , λ 1 , . . . , λ n } \{ \lambda_0, \lambda_1,...,\lambda_n \} {λ0​,λ1​,...,λn​} 。此时，利用独立的一个验证用数据集，对序列中的所有子树进行测试，计算其平方误差或者基尼系数。取误差最小的子树为最优决策树，对应的参数值也就是最优参数值。至此CART的剪枝完成

XGboost (Extreme Gradient Boosting)

原文地址：https://arxiv.org/pdf/1603.02754.pdf

XGBoost是基于决策树的集成机器学习算法, 以梯度提升（Gradient Boost）为框架。
XGBoost是由由GBDT发展而来，同样是利用加法模型与前向分步算法实现学习的优化过程，但与GBDT是有区别的。主要区别包括以下几点：

• 目标函数：XGBoost的损失函数添加了正则化项，使用正则用以控制模型的复杂度，正则项里包含了树的叶子节点个数、每个叶子节点权重（叶结点的socre值）的平方和。
• 优化方法：GBDT在优化时只使用了一阶导数信息，XGBoost在优化时使用了一、二介导数信息。
• 缺失值处理：XBGoost对缺失值进行了处理，通过学习模型自动选择最优的缺失值默认切分方向。
• 防止过拟合: XGBoost除了增加了正则项来防止过拟合,还支持行列采样的方式来防止过拟合
• 并行计算：Block存储可以并行计算等

XGboost数学原理

Xgboost使用 CART 作为基学习器，也可以使用线性分类器作为基学习器。以CART作为基学习器时，其决策规则和决策时时一样的，但CART的每个节点都具有一个权重，也就是叶节点的得分或者时预测值。

Regularized Learning Objective

给定数据集 D = { ( x i , y i ) } , ( ∣ D ∣ = n , x i ∈ R m , y i ∈ R ) D=\{(x_i, y_i)\}, (|D|=n, x_i \in R^m, y_i \in R) D={(xi​,yi​)},(∣D∣=n,xi​∈Rm,yi​∈R)，有$n$个样本，每个样本$x_i$有$m$个特征。集成树模型使用$K$个加法函数来预测输出:

$$\hat{y}=\sum _{k=1}^K{f}_k(x_i),f_k\in \mathcal{F}$$,

F = { f ( x ) = w q ( x ) } , ( q : R m − > T , w ∈ R T ) \mathcal{F} = \{ f(x) = w_{q(x)} \}, (q: R^m->T, w \in R^T) F={f(x)=wq(x)​},(q:Rm−>T,w∈RT)是回归树的空间，$q$函数将每一个样本映射到对应的叶子节点，$T$是叶子节点数量。$w$是$T$维向量，用$w_i$代表第$i$个叶子节点。对于一个给定的样本，将每一棵树的预测值叠加起来作为最终的预测值。

为了学习每一颗树的参数，定义如下目标函数
$$L=\sum _{i}l(y_i,\widehat{y_i})+\sum _k\Omega (f_k)$$

$$\Omega (f)=\gamma T+\frac{1}{2}\lambda ||w|{|}^2$$

$l$是衡量预测值$\widehat{y_i}$和真实值$y$的误差。第二项$\Omega$则表示对CART模型复杂度的惩罚，对于复杂度的衡量，考虑了每棵树的叶子结点数量以及叶子结点得分的平方和。

Gradient Boosting

上述的目标函数并不能使用传统的优化算法(例如Gradient descent). 我们使用贪心算法，累加的去优化每一颗树。举个例子，${\widehat{y_i}}^{(t)}$代表第$i$个样本在第$t$次迭代（第$t$颗树）的预测值, 那么再加入$f_t$之后要优化如下目标函数:

$$L^{(t)}=\sum _{i=1}^{n}l(y_i,{\widehat{y_i}}^{(t-1)}+f_t(x_i))+\Omega (f_t)$$

表示在第$t$步，贪心的添加$f_t$使得模型得目标函数减小最多。使用二阶泰勒展开对上式进行近似, 二阶泰勒展开:

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

对于第$t$次迭代, $l(y_i,{\widehat{y_i}}^{(t-1)}+f_t(x_i))$中, $y_i,{\widehat{y_i}}^{(t-1)}$都是确定值，如果把${\widehat{y_i}}^{(t-1)}$看成$x$， $f_t(x)$看成$\Delta x$， 那么对目标函数得到如下近似:

$$L^{(t)}=\sum _{i=1}^{n}l(y_i,{\widehat{y_i}}^{(t-1)})+g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}^2(x_i)+\Omega (f_t)$$

其中$g_i={\partial }_{{\widehat{y_i}}^{(t-1)}}l(y_i,{\widehat{y_i}}^{(t-1)})$, $h_i={\partial }_{{\widehat{y_i}}^{(t-1)}}^{2}l(y_i,{\widehat{y_i}}^{(t-1)})$.
继续对上述目标函数进行推导。由于$l(y_i,{\widehat{y_i}}^{(t-1)})$和第$t$步的优化没有关系，所以可以去除。 设 I j = { i ∣ q ( x i ) = j } I_j=\{i |q(x_i)=j\} Ij​={i∣q(xi​)=j} 为第$j$个叶子结点的样本集合，故我们的目标函数可以写成:
$$\begin{gathered} L^{(t)}=\sum _{i=1}^n[g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}^2(x_i)]+\Omega (f_t) \\ =\sum _{i=1}^n[g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}^2(x_i)]+\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2 \\ =\sum _{j=1}^T[\sum _{i\in I_j}{g}_i{w}_j+\frac{1}{2}{h}_i{w}_j^2]+\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2 \\ =\sum _{j=1}^T[\sum _{i\in I_j}{g}_i{w}_j+\frac{1}{2}(\sum _{i\in I_j}{h}_i+\lambda )w_j^2]+\gamma T \end{gathered}$$

遍历所有的样本后求每个样本的损失函数，但样本最终会落在叶子节点上，所以我们也可以遍历叶子节点，然后获取叶子节点上的样本集合，最后在求损失函数。即我们之前样本的集合，现在都改写成叶子结点的集合，由于一个叶子结点有多个样本存在，因此才有了${\sum }_{i\in I_j}{g}_i$ 和 ${\sum }_{i\in I_j}{h}_i$ 这两项， $w_j$ 为第 $j$个叶子节点取值.

为了简化表达式，定义$G_j={\sum }_{i\in I_j}{g}_i$, $H_j={\sum }_{i\in I_j}{h}_i$ , 得到简化的表达式如下：

$$L^{(t)}=\sum _{j=1}^T[G_j{w}_j+\frac{1}{2}(H_i+\lambda )w_j^2]+\gamma T$$

需要注意的是$G_j$和$H_j$是前$t-1$步得到的结果，其值可以视为常数, 只有最后一棵树的叶子结点$w_j$不确定，将目标函数$L^{(t)}$对$w_j$求导，可以求得$w_j$：
$$w_j^{\ast }=-\frac{G_j}{H_j+\lambda }$$

将$w_j^{\ast }$带入目标函数$L^{(t)}$，可以得到:

$$L^{(t)}=-\frac{1}{2}\sum _{j=1}^T\frac{G_j^2}{H_j+\lambda }+\gamma T$$

只需要将每个叶子结点的所有样本的一阶导数$g_j$和二阶导数$h_j$求和得到$G_j$和$H_j$, 然后用上述公式就可以算出目标函数值(越小越好)

Shrinkage and Column Subsampling

Shrinkage指的是对于新加的每一颗树，都乘上一个因子$\eta$. 类似于梯度下降中的学习率衰减。Shrinkage减小了每一颗新加树的影响从而提高模型得整体预测性能。
Column subsampling （feature subsampling）特征采样可以用于防止过拟合，在随机森林(random forest)中也有使用

Split finding algorithm

在决策树的生长过程中，一个非常关键的问题是如何找到叶子的节点的最优切分点，Xgboost支持两种分裂节点的方法——贪心算法(extrace greedy algorithm)和近似算法(approximate algorithm)。 首先要计算出每个特征分裂的收益. 假定在某一个结点上分裂，分裂前的目标函数可以写为：
$$L^{(1)}=-\frac{1}{2}\frac{(G_L+G_R{)}^2}{(H_L+H_R)+\lambda }+\gamma$$
分裂后的目标函数可以写为：
$$L^{(2)}=-\frac{1}{2}[\frac{G_L^2}{H_L+\lambda }+\frac{G_R^2}{H_R+\lambda }]+2\gamma$$
分裂后的收益为(目标函数越小越好)：

$$Gain=L^{(1)}-L^{(2)}=\frac{1}{2}[\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$$
该特征的分裂收益可以作为feature importance的重要依据。

Exact Greedy Algorithm

• 从深度为0的树开始，对每个叶节点枚举所有的可用特征；
• 针对每个特征，把属于该节点的训练样本根据该特征值进行升序排列，通过线性扫描的方式来决定该特征的最佳分裂点，并记录该特征的分裂收益；
• 选择收益最大的特征作为分裂特征，用该特征的最佳分裂点作为分裂位置，在该节点上分裂出左右两个新的叶节点，并为每个新节点关联对应的样本集
• 回到第 1 步，递归执行到满足特定条件为止

Approximate Algorithm

贪婪算法可以的到最优解，但当数据量太大时则无法读入内存进行计算，近似算法主要针对贪婪算法这一缺点给出了近似最优解。对于每个特征，只考察分位点可以减少计算复杂度。该算法会首先根据特征分布的分位数提出候选划分点，然后将连续型特征映射到由这些候选点划分的桶中，然后聚合统计信息找到所有区间的最佳分裂点。在提出候选切分点时有两种策略：

• Global：学习每棵树前就提出候选切分点，并在每次分裂时都采用这种分割；
• Local：每次分裂前将重新提出候选切分点。

直观上来看，Local 策略需要更多的计算步骤，而 Global 策略因为节点没有划分所以需要更多的候选点。

Sparsity-aware Split Finding

在现实问题中，输入是稀疏数据是很常见的。

XGBoost 在构建树的节点过程中为每个节点增加了一个缺省方向，当样本相应的特征值缺失时，可以被归类到缺省方向上，最优的缺省方向可以从数据中学到。至于如何学到缺省值的分支，其实很简单，分别枚举特征缺省的样本归为左右分支后的增益，选择增益最大的枚举项即为最优缺省方向。

Xgboost实践

xgboost文档：https://xgboost.readthedocs.io/en/latest/python/index.html

先来看一个多分类的例子

import copy
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn import metrics
import xgboost as xgb

def load_data():
    iris = load_iris()
    X = iris.data
    y = iris.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
    return X_train, X_test, y_train, y_test
    
def basic_modeling(X_train, X_test, y_train, y_test):
    params = {
    'booster': 'gbtree',					# 使用gradient boosting tree
    'objective': 'multi:softmax',			# 分类目标，这里是多分类
    'num_class': 3,							# 分类类别数
    'gamma': 0.1,                           # 最小loss减小量
    'max_depth': 3,							# 树的最大深度
    'monotone_constraints': (-1,-1,0,-1)，  # 单调性约束
    }
    dtrain = xgb.DMatrix(X_train, y_train)
    dtest = xgb.DMatrix(X_test)
    model = xgb.train(params, dtrain, 10)
    y_pred = model.predict(dtest)
    score = metrics.accuracy_score(y_test, y_pred)
    return model, score
  
if __name__ == "main":
    import matplotlib.pyplot as plt
    X_train, X_test, y_train, y_test = load_data()
    model, score = basic_modeling(X_train, X_test, y_train, y_test)
    print(score)

monotonic constraints

参考文档：https://xgboost.readthedocs.io/en/latest/tutorials/monotonic.html

在建立xgboost模型的过程中，如果对模型有很强的先验，即对模型的某个feature和response之前的关系有很强的先验，可以对模型添加约束， 从而提高模型得整体预测性能。
一种常见的约束类型是，某个特定的feature和response具有单调性。

$$f(x_1,x_2,...,x,...,x_n)\le f(x_1,x_2,...,x^{\prime },...,x_n)$$

whenever $x\le x^{\prime }$ is an increase constraints. 或者

$$f(x_1,x_2,...,x,...,x_n)\ge f(x_1,x_2,...,x^{\prime },...,x_n)$$

whenever $x\le x^{\prime }$ is an decrease constraints

来看一个例子。人为的构造数据集，$x_1$和$x_2$均服从正态分布，对训练集中的$y_{train}$， 添加噪音。 但是测试集合中$y_{test}$没有噪音。给定两个模型得参数，一个带有monotonic constraints, 一个没有montonic constraints, 发现给定单调约束的预测性能远远高于没有给定单调约束的模型

import math
import numpy as np
import pandas as pd
from sklearn import metrics
from sklearn.model_selection import train_test_split

# generate data
x1 = np.random.normal(0, 1, 1000)
x2 = np.random.normal(0, 1, 1000)

X = np.array([x1, x2]).T

X_train, X_test = train_test_split(X, test_size=0.25)

# add noise in train data
y_train = X_train[:,0] - X_train[:,1]  + np.random.normal(0, 2, X_train.shape[0])
# no noise in test data
y_test = X_test[:,0] - X_test[:,1]

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.scatter(X_train[:,0], y_train)
ax2.scatter(X_train[:,1], y_train)

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.scatter(X_test[:,0], y_test)
ax2.scatter(X_test[:,1], y_test)

import xgboost as xgb

params = {
'booster': 'gbtree',
'objective': 'reg:squarederror',
'gamma': 0.1,
'tree_method': "exact"
}

dtrain = xgb.DMatrix(X_train, y_train)
dtest = xgb.DMatrix(X_test, y_test)

model = xgb.train(params, dtrain, 10)
metrics.r2_score(y_test, model.predict(dtest)) # 0.629608701675472

params = {
'booster': 'gbtree',
'objective': 'reg:squarederror',
'gamma': 0.1,
'tree_method': "exact"
'monotone_constraints': (1, -1),
}
model = xgb.train(params, dtrain, 10)
metrics.r2_score(y_test, model.predict(dtest)) #     0.9429734575760605

只需要在params中加入"monotone_constraints"参数， 1代表单调递增，-1代表单调递减，0代表没有单调限制

• (1,0): An increasing constraint on the first predictor and no constraint on the second.
• (0,-1): No constraint on the first predictor and a decreasing constraint on the second.

feature importance

关于 feature importance可以参考 https://blog.csdn.net/John_xyz/article/details/85165064