机器学习之十：Xgboost

1、算法简介
    Xgboost 的全称是eXtreme Gradient Boosting，它是Gradient Boosting Machine的一个c++实现，作者为正在华盛顿大学研究机器学习的大牛陈天奇 。xgboost最大的特点在于，它能够自动利用CPU的多线程进行并行，同时在算法上加以改进提高了精度。首先回顾一下上一篇中的提升树模型目标函数：

$$Obj(\Theta)=\sum_{i=1}^nL(y_i,\hat y_i)+\sum_{k=1}^K\Omega (f_k)$$
其中 $n$为样本个数，$i$表示第 $i$个样本，$L(y_i,\hat y_i)$表示第 $i$个样本的预测误差，$K$表示有 $K$棵树。
    （1）第一部分是训练误差，也就是大家相对比较熟悉的如平方误差，而第二部分是每棵树的复杂度的和。因为现在我们的参数可以认为是在一个函数空间里面，我们不能采用传统的如SGD之类的算法来学习我们的模型，因此我们会采用一种叫做Boosting的 方式，每一次保留原来的模型不变，加入一个新的函数$f$到我们的模型中。
$$\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) \\ ... ... \\ \hat y_i^{(t)}=\sum_{k=1}^tf_k(x_i)=\hat y_i^{(t-1)}+f_t(x_i)$$
    其中 $x_i$表示第 $i$个样本点，$\hat y_i^{(t)}$表示第 $t$轮的预测模型，$f_t(x_i)$表示第 $t$轮加入的函数，我们如何选择每一轮加入什么$f$呢？答案是非常直接的，选取一个f来使得我们的目标函数尽量最大地降低，因此改写我们的目标函数：
$$Obj^{(t)}=\sum_{i=1}^nL(y_i,\hat y_i^{(t)})+\sum_{i=1}^t\Omega (f_i) \\= \sum_{i=1}^nL(y_i,\hat y_i^{(t-1)}+f_t(x_i))+\Omega (f_t) +const$$
注： $\sum_{i=1}^{t-1}\Omega (f_i)$在目标函数中变成了常数项 $const$， $n$为样本个数，$t$为树的棵数，也表示总共 $t$轮预测。
    当$L$是平方误差的时候，这时我们的目标可以写成下面的二次函数：
$$Obj^{(t)}=\sum_{i=1}^n(y_i-\hat y_i^{(t-1)}-f_t(x_i))^2+\Omega (f_t) +const \\ =\sum_{i=1}^n(y_i-\hat y_i^{(t-1)}-f_t(x_i))^2+\Omega (f_t) +const \\ =\sum_{i=1}^n[2(\hat y_i^{(t-1)}-y_i)f_t(x_i)+f_t(x_i))^2]+\Omega (f_t) +const$$
注：( $\hat y_i^{(t-1)}-y_i$)一般叫做残差，括号的平方展开后，常数项加入了最后的 $const$中， $n$为样本个数。
    更一般的，对于不是平方误差的情况，我们会采用如下的泰勒展开近似来定义一个近似的目标函数，方便我们进行这一步的计算。
目标：
$$\ Obj^{(t)}=\sum_{i=1}^nL(y_i,\hat y_i^{(t-1)}+f_t(x_i))+\Omega (f_t) +const$$
泰勒公式：
$$f(x+\Delta x)\approx f(x)+f^{'}(x)\Delta x+\frac12\quad f^{''}(x)\Delta x^2$$
定义：
$$g_i=\frac {\partial L(y_i,\hat y^{(t-1)})}{\partial \hat y^{(t-1)}}\quad \,\,\,\, h_i=\frac {\partial^2 L(y_i,\hat y^{(t-1)})}{\partial \hat y^{(t-1)}}\quad$$
展开时把 $y_i$看做常数， $\hat y_i^{(t-1)}$看做 $x$，$f_t(x_i)$看做 $\Delta x$，泰勒展开后为：
$$Obj^{(t)}\approx \sum_{i=1}^n[L(y_i,\hat y_i^{(t-1)})+g_if_t(x_i)+\frac12\quad h_if_t^2(x_i)]+\Omega (f_t) +const$$
    当我们把常数项移除之后，我们会发现如下一个比较统一的目标函数。这一个目标函数有一个非常明显的特点，它只依赖于每个数据点的在误差函数上的一阶导数和二阶导数:
$$Obj^{(t)}\approx \sum_{i=1}^n[g_if_t(x_i)+\frac12\quad h_if_t^2(x_i)]+\Omega (f_t) +const$$
    （2）以上是目标函数中训练误差的部分，接下来定义树的复杂度。 对于模型树 $f$的定义做一下细化，把树拆分成结构函数$q$（输入样本 $x$输出叶子节点索引）和叶子权重部分$w$（输入叶子节点索引输出表示节点的数值这是回归树的东西，分类树对应的是类别），结构函数 $q$把输入映射到叶子的索引号上面去，而$w$给定了每个索引号对应的叶子数值是什么。
$$f_t(x_i)=w_{q(x)}，w\in R^T，q:R_d->\{1,2,3...,T\}$$
其中 $T$为叶节点个数，$w$是叶子数值向量， $q$是树的结构。
    当我们给定了如上定义之后，我们可以定义一棵树的复杂度如下。这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的$L2$模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。
$$\Omega (f_t) =\gamma T+\frac 12\quad \lambda \sum_{j=1}^Tw_j^2$$
其中 $T$表示叶子节点的个数，$w_j$表示节点的数值（这是回归树的东西，分类树对应的是类别）, $\gamma$和 $\lambda$是超参数，用来调整前后两项所占的比重。
    直观上看，目标要求预测误差尽量小，叶子节点尽量少，节点数值尽量不极端。这个怎么看，如果某个样本label数值为4，那么第一个回归树预测3，第二个预测为1；另外一组回归树，一个预测2，一个预测2，那么倾向后一种，为什么呢？前一种情况，第一棵树学的太多，太接近4，也就意味着有较大的过拟合的风险。
    接下来是最关键的一步，在这种新的定义下，我们可以把目标函数进行如下改写，其中 $I_j$被定义为每个叶子 $j$上面样本集合$I_j={i|q(x_i)=j}$，则有：
$$Obj^{(t)}\approx \sum_{i=1}^n[g_if_t(x_i)+\frac12\quad h_if_t^2(x_i)]+\Omega (f_t) \\ =\sum_{i=1}^n[g_iw_{q(x_i)}+\frac12\quad h_iw^2_{q(x_i)}]+\gamma T+\frac 12\quad \lambda \sum_{j=1}^Tw_j^2 \\ =\sum_{j=1}^T[(\sum_{i\in I_j}g_i)w_j+\frac12\quad(\sum_{i\in I_j}h_i+\lambda)w^2_j]+\gamma T$$
这一个目标包含了 $T$个相互独立的单变量二次函数。我们可以定义
$$G_j=\sum_{i\in I_j}g_i \,\,\,\, H_j=\sum_{i\in I_j}h_i$$
那么这个目标函数可以进一步改写成如下的形式，假设我们已经知道树的结构 $q$，我们可以通过这个目标函数来求解出最好的$w$，以及最好的 $w$对应的目标函数最大的增益:
$$Obj^{(t)}=\sum_{j=1}^T[G_jw_j+\frac12\quad (H_jh_i+\lambda)w^2_j]+\gamma T$$
    这里只涉及到了如何求一个一维二次函数的最小值的问题，这两个的结果对应如下，左边是最好的w，右边是这个w对应的目标函数的值：
$$w^*_j=-\frac{G_j}{H_j+\lambda}\quad \,\,\,\,Obj=-\frac12\quad \sum_{j=1}^T\frac{G_j}{H_j+\lambda}\quad+\gamma T$$
    （3）预测值的计算
    Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)。你可以认为这个就是类似吉尼系数一样更加一般的对于树结构进行打分的函数。所以我们的算法也很简单，我们不断地枚举不同树的结构，利用这个打分函数来寻找出一个最优结构的树，加入到我们的模型中，再重复这样的操作。不过枚举所有树结构这个操作不太可行，所以常用的方法是贪心法，每一次尝试去对已有的叶子加入一个分割。对于一个具体的分割方案，我们可以获得的增益可以由如下公式计算：
$$Gain=\frac12\quad[\frac{G^2_L}{H_L+\lambda}\quad+\frac{G^2_R}{H_R+\lambda}\quad-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\quad]-\lambda$$
其中第一项为左子树分数，第二项为右子树分数，第三项为不分割我们可以拿到的分数。
通俗解释贪心策略：
    就是决策时刻按照当前目标最优化决定，说白了就是眼前利益最大化决定，“目光短浅”策略，他的优缺点细节大家自己去了解，经典背包问题等等。这里是怎么用贪心策略的呢，刚开始你有一群样本，放在第一个节点，这时候 $T$=1，$w$多少呢，不知道，是求出来的，这时候所有样本的预测值都是 $w$（这个地方自己好好理解，决策树的节点表示类别，回归树的节点表示预测值）,带入样本的label数值，此时函数就是一个关于$w$的二次函数求最小值，取最小值的点就是这个节点的预测值，最小的函数值为最小损失函数。暂停下，这里你发现了没，二次函数最优化！ 要是损失函数不是二次函数咋办，哦，泰勒展开式会否？，不是二次的想办法近似为二次。
    接着来，接下来要选个feature分裂成两个节点，变成一棵弱小的树苗，那么需要：1)确定分裂用的feature，how？最简单的是粗暴的枚举，选择loss function效果最好的那个；2)如何确立节点的 w <script type="math/tex" id="MathJax-Element-1349">w</script>以及最小的loss function，怎么做呢？对，二次函数的求最值（细节的会注意到，计算二次最值是不是有固定套路，导数=0的点）那么节奏是，选择一个feature分裂，计算loss function最小值，然后再选一个feature分裂，又得到一个loss function最小值…你枚举完，找一个效果最好的，把树给分裂，就得到了小树苗。
    在分裂的时候，你可以注意到，每次节点分裂，loss function被影响的只有这个节点的样本，因而每次分裂，计算分裂的增益（loss function的降低量）只需要关注打算分裂的那个节点的样本，就得到上述的公式。
    （4）Xgboost的优点：
Xgboost相对于普通GBDT的实现，可能具有以下的一些优势：
    1) 显式地将树模型的复杂度作为正则项加在优化目标;
    2) 公式推导里用到了二阶导数信息，而普通的GBDT只用到一阶;
    3) 允许使用column(feature) sampling来防止过拟合，借鉴了Random Forest的思想，sklearn里的gbm好像也有类似实现;
    4) 实现了一种分裂节点寻找的近似算法，用于加速和减小内存消耗;
    5) 节点分裂算法能自动利用特征的稀疏性;
    6) data事先排好序并以block的形式存储，利于并行计算;
    7) cache-aware, out-of-core computation，这个我不太懂;
    8) 支持分布式计算可以运行在MPI，YARN上，得益于底层支持容错的分布式通信框架rabit。
原文参考： https://blog.csdn.net/github_38414650/article/details/76061893
            http://www.52cs.org/?p=429&_360safeparam=88181703
2、Xgboost的实现
（1）安装Xgboost：
https://wang-shuo.github.io/2017/02/21/%E5%9C%A8Windows%E4%B8%8B%E5%AE%89%E8%A3%85XGBoost/
（2）测试用例

#coding=utf-8
import xgboost as xgb
import numpy as np

data = np.random.rand(5,10) # 5 entities, each contains 10 features
label = np.random.randint(2, size=5) # binary target
dtrain = xgb.DMatrix( data, label=label)
dtest = dtrain
param = {'bst:max_depth':2, 'bst:eta':1, 'silent':1, 'objective':'binary:logistic' }
param['nthread'] = 4
param['eval_metric'] = 'auc'
evallist  = [(dtest,'eval'), (dtrain,'train')]
num_round = 10
bst = xgb.train( param, dtrain, num_round, evallist )
print bst