菜菜sklearn——XGBoost(2)

3 XGBoost的智慧

class xgboost.XGBRegressor (kwargs，max_depth=3, learning_rate=0.1, n_estimators=100, silent=True,objective=‘reg:linear’, booster=‘gbtree’, n_jobs=1, nthread=None, gamma=0, min_child_weight=1,max_delta_step=0, subsample=1, colsample_bytree=1, colsample_bylevel=1, reg_alpha=0, reg_lambda=1,scale_pos_weight=1, base_score=0.5, random_state=0, seed=None, missing=None, importance_type=‘gain’)

3.1 选择弱评估器：重要参数booster

梯度提升算法中不只有梯度提升树，XGB作为梯度提升算法的进化，自然也不只有树模型一种弱评估器。在XGB中，除了树模型，我们还可以选用线性模型，比如线性回归，来进行集成。虽然主流的XGB依然是树模型，但我们也可以使用其他的模型。基于XGB的这种性质，我们有参数“booster"来控制我们究竟使用怎样的弱评估器。
两个参数都默认为"gbtree"，如果不想使用树模型，则可以自行调整。当XGB使用线性模型的时候，它的许多数学过程就与使用普通的Boosting集成非常相似，因此我们在讲解中会重点来讲解使用更多，也更加核心的基于树模型的XGBoost。简单看看现有的数据集上，是什么样的弱评估器表现更好：

for booster in ["gbtree","gblinear","dart"]:
    reg = XGBR(n_estimators=180
               ,learning_rate=0.1
               ,random_state=420
               ,booster=booster).fit(Xtrain,Ytrain)
    print(booster)
    print(reg.score(Xtest,Ytest))
#gbtree
#0.9231068620728082
#gblinear
#0.6286510307485139
#dart
#0.923106843149575

3.2 XGB的目标函数：重要参数objective

梯度提升算法中都存在着损失函数。不同于逻辑回归和SVM等算法中固定的损失函数写法，集成算法中的损失函数是可选的，要选用什么损失函数取决于我们希望解决什么问题，以及希望使用怎样的模型。比如说，如果我们的目标是进行回归预测，那我们可以选择调节后的均方误差RMSE作为我们的损失函数。如果我们是进行分类预测，那我们可以选择错误率error或者对数损失log_loss。只要我们选出的函数是一个可微的，能够代表某种损失的函数，它就可
以是我们XGB中的损失函数。

在众多机器学习算法中，损失函数的核心是衡量模型的泛化能力，即模型在未知数据上的预测的准确与否，我们训练模型的核心目标也是希望模型能够预测准确。在XGB中，预测准确自然是非常重要的因素，但我们之前提到过，XGB的是实现了模型表现和运算速度的平衡的算法。普通的损失函数，比如错误率，均方误差等，都只能够衡量模型的表现，无法衡量模型的运算速度。回忆一下，我们曾在许多模型中使用空间复杂度和时间复杂度来衡量模型的运算效
率。XGB因此引入了模型复杂度来衡量算法的运算效率。因此XGB的目标函数被写作：传统损失函数 + 模型复杂度。
$$Obj=\sum _{i=1}^{m}l\left(y_i,{\hat{y}}_i\right)+\sum _{k=1}^K\Omega \left(f_k\right)$$
其中i代表数据集中的第i个样本，m表示导入第k棵树的数据总量，K代表建立的所有树(n_estimators)，当只建立了t棵树的时候，式子应当为${\sum }_{k=1}^t\Omega \left(f_k\right)$。第一项代表传统的损失函数，衡量真实标签$y_i$与预测值${\hat{y}}_i$之间的差异，通常是RMSE，调节后的均方误差。第二项代表模型的复杂度，使用树模型的某种变换$\Omega$表示，这个变化代表了一个从树的结构来衡量树模型的复杂度的式子，可以有多种定义。注意，我们的第二项中没有特征矩阵$x_i$的介入。我们在迭代每一棵树的过程中，都最小化Obj来力求获取最优的$\hat{y}$，因此我们同时最小化了模型的错误率和模型的复杂度，这种设计目标函数的方法不得不说实在是非常巧妙和聪明。(正则化项（1）防止过拟合（2）降低模型复杂度，提高模型运算速度)
其实并非如此，我们的第一项传统损失函数也是与已经建好的所有树相关的，相关在这里：
$${\hat{y}}_i^{(t)}=\sum _k^t{f}_k\left(x_i\right)=\sum _k^{t-1}{f}_k\left(x_i\right)+f_t\left(x_i\right)$$
我们的${\hat{y}}_i$中已经包含了所有树的迭代结果，因此整个目标函数都与K棵树相关。

我们还可以从另一个角度去理解我们的目标函数。回忆一下我们曾经在随机森林中讲解过的方差-偏差困境。在机器学习中，我们用来衡量模型在未知数据上的准确率的指标，叫做泛化误差（Genelization error）。一个集成模型(f)在未知数据集(D)上的泛化误差$E(f;D)$，由方差(var)，偏差(bais)和噪声(ε)共同决定，而泛化误差越小，模型就越理想。从下面的图可以看出来，方差和偏差是此消彼长的，并且模型的复杂度越高，方差越大，偏差越小。
方差可以被简单地解释为模型在不同数据集上表现出来地稳定性，而偏差是模型预测的准确度。那方差-偏差困境就可以对应到我们的Obj中了：
$$Obj=\sum _{i=1}^{m}l\left(y_i,{\hat{y}}_i\right)+\sum _{k=1}^K\Omega \left(f_k\right)$$
第一项是衡量我们的偏差，模型越不准确，第一项就会越大。第二项是衡量我们的方差，模型越复杂，模型的学习就会越具体，到不同数据集上的表现就会差异巨大，方差就会越大。所以我们求解Obj的最小值，其实是在求解方差与偏差的平衡点，以求模型的泛化误差最小，运行速度最快。我们知道树模型和树的集成模型都是学习天才，是天生过拟合的模型，因此大多数树模型最初都会出现在图像的右上方，我们必须通过剪枝来控制模型不要过拟合。现在XGBoost的损失函数中自带限制方差变大的部分，也就是说XGBoost会比其他的树模型更加聪明，不会轻易落到图像的右上方。可见，这个模型在设计的时候的确是考虑了方方面面，无怪XGBoost会如此强大了。

在应用中，我们使用参数“objective"来确定我们目标函数的第一部分中的$l\left(y_i,{\hat{y}}_i\right)$，也就是衡量损失的部分。
常用的选择有：

我们还可以选择自定义损失函数。比如说，我们可以选择输入平方损失$l\left(y_j,{\hat{y}}_j\right)={\left(y_j-{\hat{y}}_j\right)}^2$，此时我们的XGBoost其实就是算法梯度提升机器（gradient boosted machine）。在xgboost中，我们被允许自定义损失函数，但通常我们还是使用类已经为我们设置好的损失函数。我们的回归类中本来使用的就是reg:linear，因此在这里无需做任何调整。注意：分类型的目标函数导入回归类中会直接报错。现在来试试看xgb自身的调用方式。
由于xgb中所有的参数都需要自己的输入，并且objective参数的默认值是二分类，因此我们必须手动调节。试试看在其他参数相同的情况下，我们xgboost库本身和sklearn比起来，效果如何：

#默认reg:linear
reg = XGBR(n_estimators=180,random_state=420).fit(Xtrain,Ytrain)
reg.score(Xtest, Ytest)
#0.9231068620728082
MSE(Ytest,reg.predict(Xtest))
#7.155205217161047

#xgb实现法
import xgboost as xgb
#使用类DMatrix读取数据
dtrain = xgb.DMatrix(Xtrain,Ytrain) #特征矩阵和标签都进行一个传入
dtest = xgb.DMatrix(Xtest,Ytest)
#非常遗憾无法打开来查看，所以通常都是先读到pandas里面查看之后再放到DMatrix中
dtrain
#<xgboost.core.DMatrix at 0x2770de3bdd8>
pd.DataFrame(dtrain)
#报错
import pandas as pd
pd.DataFrame(Xtrain)

一般先把数据读到pandas里面，查看一下数据，进行探索、预处理，然后用sklearn分成测试集、训练集，再放到DMatrix

#写明参数
param = {'silent':True #默认为False，通常要手动把它关闭掉
         ,'objective':'reg:linear'
         ,"eta":0.1}
num_round = 180 #n_estimators
#类train，可以直接导入的参数是训练数据，树的数量，其他参数都需要通过params来导入
bst = xgb.train(param, dtrain, num_round)
#接口predict
preds = bst.predict(dtest)
preds
#array([ 6.4613175, 22.123888 , 30.755163 , 13.424351 ,  8.378565 ,
#       23.608477 , 14.2151165, 16.026499 , 15.498961 , 14.10649  ,
#       24.030867 , 34.36362  , 21.461111 , 28.839497 , 19.568035 ,
#
#       41.85804  , 17.44544  , 22.907183 , 21.02398  , 29.799366 ,
#       20.219465 , 12.404763 , 45.750965 , 25.56757  , 22.000706 ,
#       14.194921 , 27.102774 ], dtype=float32)
from sklearn.metrics import r2_score
r2_score(Ytest,preds)
#0.9260984298390122
MSE(Ytest,preds)
#6.87682821415069

看得出来，无论是从$R^2$还是从MSE的角度来看，都是xgb库本身表现更优秀，这也许是由于底层的代码是由不同团队创造的缘故。随着样本量的逐渐上升，sklearnAPI中调用的结果与xgboost中直接训练的结果会比较相似，如果希望的话可以分别训练，然后选取泛化误差较小的库。如果可以的话，建议脱离sklearnAPI直接调用xgboost库，因为xgboost库本身的调参要方便许多。

3.3 求解XGB的目标函数

有了如下目标函数，我们来看看如何求解它。时刻记得我们求解目标函数的目的：为了求得在第t次迭代中最优的树$f_t$。在逻辑回归和支持向量机中，我们通常先将目标函数转化成一种容易求解的方式（比如对偶），然后使用梯度下降或者SMO之类的数学方法来执行我们的最优化过程。之前我们使用了逻辑回归的迭代过程来帮助大家理解在梯度提升树中树是如何迭代的，那我们是否可以使用逻辑回归的参数求解方式来求解XGB的目标函数呢？

很遗憾，在XGB中我们无法使用梯度下降，原因是XGB的损失函数没有需要求解的参数。我们在传统梯度下降中迭代的是参数，而我们在XGB中迭代的是树，树$f_k$不是数字组成的向量，并且其结构不受到特征矩阵$x$取值大小的直接影响，尽管这个迭代过程可以被类比到梯度下降上，但真实的求解过程却是完全不同。

在求解XGB的目标函数的过程中，我们考虑的是如何能够将目标函数转化成更简单的，与树的结构直接相关的写法，以此来建立树的结构与模型的效果（包括泛化能力与运行速度）之间的直接联系。也因为这种联系的存在，XGB的目标函数又被称为“结构分数”。
$$\begin{array}{rl}{\hat{y}}_i^{(t)}=\sum _k^t{f}_k\left(x_i\right)& =\sum _k^{t-1}{f}_k\left(x_i\right)+f_t\left(x_i\right)\\ & ={\hat{y}}_i^{(t-1)}+f_t\left(x_i\right)\end{array}$$
首先，我们先来进行第一步转换。

其中$g_i$和$g_i$分别是在损失函数$l\left(y_i^t,{\hat{y}}_i^{(t-1)}\right)$上对${\hat{y}}_i^{(t-1)}$所求的一阶导数和二阶导数，他们被统称为每个样本的梯度统计量（gradient statisticts）。在许多算法的解法推导中，我们求解导数都是为了让一阶导数等于0来求解极值，而现在我们求解导数只是为了配合泰勒展开中的形式，仅仅是简化公式的目的罢了。所以GBDT和XGB的区别之中，GBDT求一阶导数，XGB求二阶导数，这两个过程根本是不可类比的。XGB在求解极值为目标的求导中也是求解一阶导数，后面会为大家详解。

$$f(x)=\frac{f(c)}{0!}+\frac{f^{\prime }(c)}{1!}(x-c)+\frac{f^{\prime \prime }(c)}{2!}(x-c{)}^2+\frac{f^{\prime \prime \prime }(c)}{3!}(x-c{)}^3+\dots$$
其中$f^{\prime }(c)$表示$f(x)$上对x求导后，令x的值等于c所取得的值。其中有假设：c与x非常接近，(x-c)非常接近0，于是我们可以将式子改写成：
$$f(x+x-c)\approx \frac{f(c)}{0!}+\frac{f^{\prime }(c)}{1!}(x-c)+\frac{f^{\prime \prime }(c)}{2!}(x-c{)}^2+\frac{f^{\prime \prime \prime }(c)}{3!}(x-c{)}^3+\dots$$
只 取 前 两 项 ， 取 约 等 于 ：
$$\approx \frac{f(c)}{0!}+\frac{f^{\prime }(c)}{1!}(x-c)+\frac{f^{\prime \prime }(c)}{2!}(x-c{)}^2$$$$\approx f(c)+f^{\prime }(c)(x-c)+\frac{f^{\prime \prime }(c)}{2}(x-c{)}^2$$
令 $x_1=x,x_2=x-c:$
$$f\left(x_1+x_2\right)\approx f\left(x_1\right)+f^{\prime }\left(x_1\right)\ast x_2+\frac{f^{\prime \prime }\left(x_1\right)}{2}\ast x_2^2$$$$F\left({\hat{y}}_i^{(t-1)}+f_t\left(x_i\right)\right)\approx F\left({\hat{y}}_i^{(t-1)}\right)+g_i\ast f_t\left(x_i\right)+\ast \frac{h_i}{2}{\left(f_t\left(x_i\right)\right)}^2$$
如刚才所说，x-c需要很小，与x相比起来越小越好，那在我们的式子中，需要很小的这部分就是$f_t(x_i)$。这其实很好理解，对于一个集成算法来说，每次增加的一棵树对模型的影响其实非常小，尤其是当我们有许多树的时候——比如，n_estimators=500的情况，$f_t(x_i)$与x相比总是非常小的，因此这个条件可以被满足，泰勒展开可以被使用。如此，我们的目标函数可以被顺利转化成：
$$Obj=\sum _{i=1}^m\left[f_t\left(x_i\right)g_i+\frac{1}{2}{\left(f_t\left(x_i\right)\right)}^2{h}_i\right)]+\Omega \left(f_t\right)$$
这个式子中，$g_i$和$h_i$只与传统损失函数相关，核心的部分是我们需要决定的树$f_t$。接下来，我们就来研究一下我们的$f_t$。

3.4 参数化决策树$f_k(x)$：参数alpha，lambda

class xgboost.XGBRegressor (max_depth=3, learning_rate=0.1, n_estimators=100, silent=True,objective=‘reg:linear’, booster=‘gbtree’, n_jobs=1, nthread=None, gamma=0, min_child_weight=1,max_delta_step=0, subsample=1, colsample_bytree=1, colsample_bylevel=1, reg_alpha=0, reg_lambda=1,scale_pos_weight=1, base_score=0.5, random_state=0, seed=None, missing=None, importance_type=‘gain’,*kwargs)

在参数化决策树之前，我们先来简单复习一下回归树的原理。对于决策树而言，每个被放入模型的任意样本 最终一个都会落到一个叶子节点上。对于回归树，通常来说每个叶子节点上的预测值是这个叶子节点上所有样本的标签的均值。但值得注意的是，XGB作为普通回归树的改进算法，在$\hat{y}$上却有所不同。

对于XGB来说，每个叶子节点上会有一个预测分数（prediction score），也被称为叶子权重。这个叶子权重就是所有在这个叶子节点上的样本在这一棵树上的回归取值，用$f_k(x)$或者w来表示。
当有多棵树的时候，集成模型的回归结果就是所有树的预测分数之和，假设这个集成模型中总共有 棵决策树，则整个模型在这个样本i上给出的预测结果为：
$${\hat{y}}_i^{(k)}=\sum _k^K{f}_k\left(x_i\right)$$

基于这个理解，我们来考虑每一棵树。对每一棵树，它都有自己独特的结构，这个结构即是指叶子节点的数量，树的深度，叶子的位置等等所形成的一个可以定义唯一模型的树结构。在这个结构中，我们使用$q(x_i)$表示样本$x_i$所在的叶子节点，并且使用$w_{q(x_i)}$来表示这个样本落到第t棵树上的第$q(x_i)$个叶子节点中所获得的分数，于是有：
$$f_t\left(x_i\right)=w_{q\left(x_i\right)}$$
这是对于每一个样本而言的叶子权重，然而在一个叶子节点上的所有样本所对应的叶子权重是相同的。设一棵树上总共包含了T个叶子节点，其中每个叶子节点的索引为j，则这个叶子节点上的样本权重是$w_j$。依据这个，我们定义模型的复杂度$\Omega (f)$为（注意这不是唯一可能的定义，我们当然还可以使用其他的定义，只要满足叶子越多/深度越大，复杂度越大的理论，我们可以自己决定我们的$\Omega (f)$要是一个怎样的式子）：
$$\Omega (f)=\gamma T+\text{ 正则项 (Regularization) }$$
如果使用L2正则项 ：
$$=\gamma T+\frac{1}{2}\lambda \Vert w{\Vert }^2$$$$=\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2$$
如果使用L1正则项 ：
$$=\gamma T+\frac{1}{2}\alpha |w|$$$$=\gamma T+\frac{1}{2}\alpha \sum _{j=1}^T\left|w_j\right|$$
还可以两个一起使用 ：
$$=\gamma T+\frac{1}{2}\alpha \sum _{j=1}^T\left|w_j\right|+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2$$
这个结构中有两部分内容，一部分是控制树结构的$\gamma T$，另一部分则是我们的正则项。叶子数量T可以代表整个树结构，这是因为在XGBoost中所有的树都是CART树（二叉树），所以我们可以根据叶子的数量T判断出树的深度，而$\gamma$是我们自定的控制叶子数量的参数。

至于第二部分正则项，类比一下我们岭回归和Lasso的结构，参数$\alpha$和$\lambda$的作用其实非常容易理解，他们都是控制正则化强度的参数，我们可以二选一使用，也可以一起使用加大正则化的力度。当$\alpha$和$\lambda$都为0的时候，目标函数就是普通的梯度提升树的目标函数。
来看正则化系数分别对应的参数：

根据我们以往的经验，我们往往认为两种正则化达到的效果是相似的，只不过细节不同。比如在逻辑回归当中，两种正则化都会压缩$\theta$参数的大小，只不过L1正则化会让$\theta$为0，而L2正则化不会。在XGB中也是如此。当 $\alpha$和$\lambda$越大，惩罚越重，正则项所占的比例就越大，在尽全力最小化目标函数的最优化方向下，叶子节点数量就会被压制，模型的复杂度就越来越低，所以对于天生过拟合的XGB来说，正则化可以一定程度上提升模型效果。

对于两种正则化如何选择的问题，从XGB的默认参数来看，我们优先选择的是L2正则化。当然，如果想尝试L1也不是不可。两种正则项还可以交互，因此这两个参数的使用其实比较复杂。在实际应用中，正则化参数往往不是我们调参的最优选择，如果真的希望控制模型复杂度，我们会调整$\gamma$而不是调整这两个正则化参数，因此大家不必过于在意这两个参数最终如何影响了我们的模型效果。对于树模型来说，还是剪枝参数地位更高更优先。大家只需要理解这两个参数从数学层面上如何影响我们的模型就足够了。如果我们希望调整$\alpha$和$\lambda$，我们往往会使用网格搜索来帮助我们。在这里，为大家贴出网格搜索的代码和结果供大家分析和学习。

#使用网格搜索来查找最佳的参数组合
from sklearn.model_selection import GridSearchCV
param = {"reg_alpha":np.arange(0,5,0.05),"reg_lambda":np.arange(0,2,0.05)}
gscv = GridSearchCV(reg,param_grid = param,scoring = "neg_mean_squared_error",cv=cv)
#======【TIME WARNING：10~20 mins】======#
time0=time()
gscv.fit(Xtrain,Ytrain)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
gscv.best_params_
gscv.best_score_
preds = gscv.predict(Xtest)
from sklearn.metrics import r2_score，mean_squared_error as MSE
r2_score(Ytest,preds)
MSE(Ytest,preds) #网格搜索的结果有什么样的含义呢？为什么会出现这样的结果？你相信网格搜索得出的结果吗？试着用数学和你对XGB的理解来解释一下吧。

3.5 寻找最佳树结构：求解w与T

在上一节中，我们定义了树和树的复杂度的表达式，树我们使用叶子节点上的预测分数来表达，而树的复杂度则是叶子数目加上正则项：
$$f_t\left(x_i\right)=w_{q\left(x_i\right)},\Omega \left(f_t\right)=\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2$$
假设现在第t棵树的结构已经被确定为q，我们可以将树的结构带入我们的损失函数，来继续转化我们的目标函数。转化目标函数的目的是：建立树的结构（叶子节点的数量）与目标函数的大小之间的直接联系，以求出在第t次迭代中需要求解的最优的树$f_t$。注意，我们假设我们使用的是L2正则化（这也是参数lambda和alpha的默认设置，lambda为1，alpha为0），因此接下来的推导也会根据L2正则化来进行。
橙色框中的转化是如何实现的？

于是我们可以有：
$$\sum _{i=1}^m{w}_{q\left(x_i\right)}\ast g_i=w_{q\left(x_1\right)}\ast g_1+w_{q\left(x_2\right)}\ast g_2+w_{q\left(x_3\right)}\ast g_3$$$$=w_1\left(g_1+g_2\right)+w_2\ast g_3$$$$=\sum _{j=1}^T\left(w_j\sum _{i\in I_j}{g}_i\right)$$
如此就实现了这个转化。

对于最终的式子，我们定义：
$$G_j=\sum _{i\in I_j}{g}_i,H_j=\sum _{i\in I_j}{h}_i$$
于是可以有：
$$Ob{j}^{(t)}=\sum _{j=1}^T\left[w_j{G}_j+\frac{1}{2}{w}_j^2\left(H_j+\lambda \right)\right]+\gamma T$$$$F^{\ast }\left(w_j\right)=w_j{G}_j+\frac{1}{2}{w}_j^2\left(H_j+\lambda \right)$$
其中每个j取值下都是一个以$w_j$为自变量的二次函数$F\ast$，我们的目标是追求让Obj最小，只要单独的每一个叶子j取值下的二次函数都最小，那他们的加和必然也会最小。于是，我们在$F\ast$上对$w_j$求导，让一阶导数等于0以求极值，可得：
$$\frac{\partial F^{\ast }\left(w_j\right)}{\partial w_j}=G_j+w_j\left(H_j+\lambda \right)$$$$0=G_j+w_j\left(H_j+\lambda \right)$$$$w_j=-\frac{G_j}{H_j+\lambda }$$
把这个公式带入目标函数，则有：
$$Ob{j}^{(t)}=\sum _{j=1}^T\left[-\frac{G_j}{H_j+\lambda }\ast G_j+\frac{1}{2}{\left(-\frac{G_j}{H_j+\lambda }\right)}^2\left(H_j+\lambda \right)\right]+\gamma T$$$$=\sum _{j=1}^T\left[-\frac{G_j^2}{H_j+\lambda }+\frac{1}{2}\ast \frac{G_j^2}{H_j+\lambda }\right]+\gamma T$$$$=-\frac{1}{2}\sum _{j=1}^T\frac{G_j^2}{H_j+\lambda }+\gamma T$$
到了这里，大家可能已经注意到了，比起最初的损失函数 + 复杂度的样子，我们的目标函数已经发生了巨大变化。我们的样本量i已经被归结到了每个叶子当中去，我们的目标函数是基于每个叶子节点，也就是树的结构来计算。所以，我们的目标函数又叫做“结构分数”（structure score），分数越低，树整体的结构越好。如此，我们就建立了树的结构（叶子）和模型效果的直接联系。

更具体一些，我们来看一个例子：
$$Obj=-\left(\frac{g_1^2}{h_1+\lambda }+\frac{g_4^2}{h_4+\lambda }+\frac{{\left(g_2+g_3+g_5\right)}^2}{h_2+h_3+h_5+\lambda }\right)+3\gamma$$
所以在XGB的运行过程中，我们会根据Obj的表达式直接探索最好的树结构，也就是说找寻最佳的树。从式子中可以看出，$\lambda$和$\gamma$是我们设定好的超参数，$G_j$和$H_j$是由损失函数和这个特定结构下树的预测结果${\hat{y}}_i^{(t-1)}$共同决定，而T只由我们的树结构决定。则我们通过最小化Obj所求解出的其实是T，叶子的数量。所以本质也就是求解树的结构了。

在这个算式下，我们可以有一种思路，那就是枚举所有可能的树结构1，然后一个个计算我们的Obj，待我们选定了最佳的树结构（最佳的T）之后，我们使用这种树结构下计算出来的$G_j$和$H_j$就可以求解出每个叶子上的权重$w_j$ ，如此就找到我们的最佳树结构，完成了这次迭代。

• 用w求解w？

一个大家可能会感到困惑的点是，$G_j$和$H_j$的本质其实是损失函数上的一阶导数$g_i$和二阶导数$h_i$之和，而一阶和二阶导数本质是：
$$g_i=\frac{\partial l\left(y_i,{\hat{y}}_i^{(t-1)}\right)}{\partial {\hat{y}}_i^{(t-1)}},h_i=\frac{{\partial }^{2}l\left(y_i,{\hat{y}}_i^{(t-1)}\right)}{\partial {\left({\hat{y}}_i^{(t-1)}\right)}^2}$$
$y_i$是已知的标签，但我们有预测值的求解公式：
$$\begin{array}{rl}{\widehat{y_i}}^{(k)}& =\sum _k^K{f}_k\left(x_i\right)\\ & =\sum _k^K{w}_{q\left(x_i\right)}\end{array}$$
这其实是说，$G_j$和$H_j$的计算中带有w，那先确定最佳的T，再求出$G_j$和$H_j$，结合$\lambda$求出叶子权重$w_i$的思路不觉得有些问题么？仿佛在带入w求解w？对于有这样困惑的大家，请注意我们的${\hat{y}}_i^{(t-1)}$与我们现在要求解的$w_j$其实不是在同一棵树上的。别忘记我们是在一直迭代的，我们现在求解的$w_j$是第t棵树上的结果，而${\hat{y}}_i^{(t-1)}$是前面的t-1棵树的累积w ，是在前面所有的迭代中已经求解出来的已知的部分。

• 求解第一棵树时，没有“前面已经迭代完毕的部分”，怎么办？

那第二个问题又来了：那我们如何求解第一棵树在样本i下的预测值${\hat{y}}_i^{(1)}$呢？在建立第一棵树时，并不存在任何前面的迭代已经计算出来的信息，但根据公式我们需要使用如下式子来求解$f_1\left(x_i\right)$，并且我们在求解过程中还需要对前面所有树的结果进行求导。
$${\hat{y}}_i^{(1)}={\hat{y}}_i^{(0)}+f_1\left(x_i\right)$$
这时候，我们假设${\hat{y}}_i^{(0)}=0$来解决我们的问题，事实是，由于前面没有已经测出来的树的结果，整个集成算法的结果现在也的确为0。

• 第0棵树的预测值假设为0，那求解第一棵树的$g_i$和$h_i$的过程是在对0求导？

这个问题其实很简单。在进行求导时，所有的求导过程都和之前推导的过程相一致，之所以能够这么做，是因为我们其实不是在对0求导，而是对一个变量${\hat{y}}_i^{(t-1)}$求导。只是除了求导之外，我们还需要在求导后的结果中带入这个变量此时此刻的取值，而这个取值在求解第一棵树时刚好等于0而已。更具体地，可以看看下面，对0求导，和对变量求导后，变量取值为0的区别：
$$\text{ 对常数 }0\text{ 进行求导: }\frac{\partial \left(x^2+x\right)}{\partial (0)}=0$$$$\text{ 对变量 }x\text{ 进行求导，但变量 }x\text{ 等于 }0:\frac{\partial \left(x^2+x\right)}{\partial (x)}=2x+1=2\ast 0+1=1$$
这些细节都理解了之后，相信大家对于先求解Obj的最小值来求解树结构T，然后以此为基础求解出$w_j$的过程已经没有什么问题了。回忆一下我们刚才说的，为了找出最小的Obj，我们需要枚举所有可能的树结构，这似乎又回到了我们最初的困惑——我们之所以要使用迭代和最优化的方式，就是因为我们不希望进行枚举，这样即浪费资源又浪费时间。那我们怎么办呢？来看下一节：贪婪算法寻找最佳结构。

3.6 寻找最佳分枝：结构分数之差

贪婪算法指的是控制局部最优来达到全局最优的算法，决策树算法本身就是一种使用贪婪算法的方法。XGB作为树的集成模型，自然也想到采用这样的方法来进行计算，所以我们认为，如果每片叶子都是最优，则整体生成的树结构就是最优，如此就可以避免去枚举所有可能的树结构。

回忆一下决策树中我们是如何进行计算：我们使用基尼系数或信息熵来衡量分枝之后叶子节点的不纯度，分枝前的信息熵与分治后的信息熵之差叫做信息增益，信息增益最大的特征上的分枝就被我们选中，当信息增益低于某个阈值时，就让树停止生长。在XGB中，我们使用的方式是类似的：我们首先使用目标函数来衡量树的结构的优劣，然后让树从深度0开始生长，每进行一次分枝，我们就计算目标函数减少了多少，当目标函数的降低低于我们设定的某个阈值时，就让树停止生长。

来个具体的例子，在这张图中，我们有中间节点“是否是男性”，这个中间节点下面有两个叶子节点，分别是样本弟弟和妹妹。我们来看看这个分枝点上，树的结构分数之差如何表示。

对于中间节点这一个叶子节点而言，我们的T=1 ，则这个节点上的结构分数为：
I = { 1 , 4 } G = g 1 + g 4 H = h 1 + h 4  Score  middle  = − 1 2 G 2 H + λ + γ \begin{aligned} I &=\{1,4\} \\ G &=g_{1}+g_{4} \\ H &=h_{1}+h_{4} \\ \text { Score }_{\text {middle }} &=-\frac{1}{2} \frac{G^{2}}{H+\lambda}+\gamma \end{aligned} IGH Score middle ​​={1,4}=g1​+g4​=h1​+h4​=−21​H+λG2​+γ​
对于弟弟和妹妹节点而言，则有：
$$\begin{array}{rl}{\text{ Score }}_{sis}& =-\frac{1}{2}\frac{g_4^2}{h_4+\lambda }+\gamma \\ {\text{ Score }}_{bro}& =-\frac{1}{2}\frac{g_1^2}{h_1+\lambda }+\gamma \end{array}$$
则分枝后的结构分数之差为：
$$\begin{array}{rl}\text{ Gain }& ={\text{ Score }}_{\text{sis }}+{\text{ Score }}_{\text{bro }}-{\text{ Score }}_{\text{middle }}\\ & =-\frac{1}{2}\frac{g_4^2}{h_4+\lambda }+\gamma -\frac{1}{2}\frac{g_1^2}{h_1+\lambda }+\gamma -\left(-\frac{1}{2}\frac{G^2}{H+\lambda }+\gamma \right)\\ & =-\frac{1}{2}\frac{g_4^2}{h_4+\lambda }+\gamma -\frac{1}{2}\frac{g_1^2}{h_1+\lambda }+\gamma +\frac{1}{2}\frac{G^2}{H+\lambda }-\gamma \\ & =-\frac{1}{2}\left[\frac{g_4^2}{h_4+\lambda }+\frac{g_1^2}{h_1+\lambda }-\frac{G^2}{H+\lambda }\right]+\gamma \\ & =-\frac{1}{2}\left[\frac{g_4^2}{h_4+\lambda }+\frac{g_1^2}{h_1+\lambda }-\frac{{\left(g_1+g_4\right)}^2}{\left(h_1+h_4\right)+\lambda }\right]+\gamma \end{array}$$
将 负 号 去 除 ：
$$=\frac{1}{2}\left[\frac{G_L^2}{H_L+\lambda }+\frac{G_R^2}{H_R+\lambda }-\frac{{\left(G_L+G_R\right)}^2}{H_L+H_R+\lambda }\right]-\gamma$$
CART树全部是二叉树，因此这个式子是可以推广的。从这个式子我们可以总结出，其实分枝后的结构分数之差为：
$$\text{ Gain }=\frac{1}{2}\left[\frac{G_L^2}{H_L+\lambda }+\frac{G_R^2}{H_R+\lambda }-\frac{{\left(G_L+G_R\right)}^2}{H_L+H_R+\lambda }\right]-\gamma$$
其中$G_L$和$H_L$从左节点（弟弟节点）上计算得出，$G_R$和$H_R$从有节点（妹妹节点）上计算得出，而$G_L+G_R$和$H_L+H_R$从中间节点上计算得出。对于任意分枝，我们都可以这样来进行计算。在现实中，我们会对所有特征的所有分枝点进行如上计算，然后选出让目标函数下降最快的节点来进行分枝。对每一棵树的每一层，我们都进行这样的计算，比起原始的梯度下降，实践证明这样的求解最佳树结构的方法运算更快，并且在大型数据下也能够表现不错。至此，我们作为XGBoost的使用者，已经将需要理解的XGB的原理理解完毕了。

3.7 让树停止生长：重要参数gamma

在之前所有的推导过程中，我们都没有提到$\gamma$这个变量。从目标函数和结构分数之差Gain的式子中来看，$\gamma$是我们每增加一片叶子就会被剪去的惩罚项。增加的叶子越多，结构分数之差Gain会被惩罚越重，所以$\gamma$又被称之为是“复杂性控制”（complexity control），所以$\gamma$是我们用来防止过拟合的重要参数。实践证明，$\gamma$是对梯度提升树影响最大的参数之一，其效果丝毫不逊色于n_estimators和防止过拟合的神器max_depth。同时，$\gamma$还是我们让树停止生长的重要参数。

在逻辑回归中，我们使用参数tol来设定阈值，并规定如果梯度下降时损失函数减小量小于tol下降就会停止。在XGB中，我们规定，只要结构分数之差$Gain$是大于0的，即只要目标函数还能够继续减小，我们就允许树继续进行分枝。也就是说，我们对于目标函数减小量的要求是：
$$\frac{1}{2}\left[\frac{G_L^2}{H_L+\lambda }+\frac{G_R^2}{H_R+\lambda }-\frac{{\left(G_L+G_R\right)}^2}{H_L+H_R+\lambda }\right]-\gamma >0$$$$\frac{1}{2}\left[\frac{G_L^2}{H_L+\lambda }+\frac{G_R^2}{H_R+\lambda }-\frac{{\left(G_L+G_R\right)}^2}{H_L+H_R+\lambda }\right]>\gamma$$
如此，我们可以直接通过设定$\gamma$的大小来让XGB中的树停止生长。$\gamma$因此被定义为，在树的叶节点上进行进一步分枝所需的最小目标函数减少量，在决策树和随机森林中也有类似的参数（min_split_loss，min_samples_split）。$\gamma$设定越大，算法就越保守，树的叶子数量就越少，模型的复杂度就越低。

如果我们希望从代码中来观察$\gamma$的作用，使用sklearn中传统的学习曲线等工具就比较困难了。来看下面这段代码，这是一段让参数$\gamma$在0~5之间均匀取值的学习曲线。其运行速度较缓慢并且曲线的效果匪夷所思，大家若感兴趣可以自己运行一下。

#======【TIME WARNING: 1 min】=======#
axisx = np.arange(0,5,0.05)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators=180,random_state=420,gamma=i)
    result = CVS(reg,Xtrain,Ytrain,cv=cv)
    rs.append(result.mean())
    var.append(result.var())
    ge.append((1 - result.mean())**2+result.var())
print(axisx[rs.index(max(rs))],max(rs),var[rs.index(max(rs))])
print(axisx[var.index(min(var))],rs[var.index(min(var))],min(var))
print(axisx[ge.index(min(ge))],rs[ge.index(min(ge))],var[ge.index(min(ge))],min(ge))
rs = np.array(rs)
var = np.array(var)*0.1
plt.figure(figsize=(20,5))
plt.plot(axisx,rs,c="black",label="XGB")
plt.plot(axisx,rs+var,c="red",linestyle='-.')
plt.plot(axisx,rs-var,c="red",linestyle='-.')
plt.legend()
plt.show()

以上这段代码的运行结果如下：

可以看到，我们完全无法从中看出什么趋势，偏差时高时低，方差时大时小，参数$\gamma$引起的波动远远超过其他参数（其他参数至少还有一个先上升再平稳的过程，而$\gamma$则是仿佛完全无规律）。在sklearn下XGBoost太不稳定，如果这样来调整参数的话，效果就很难保证$\gamma$。因此，为了调整 ，我们需要来引入新的工具，xgboost库中的类xgboost.cv。

xgboost.cv (params, dtrain, num_boost_round=10, nfold=3, stratified=False, folds=None, metrics=(), obj=None,feval=None, maximize=False, early_stopping_rounds=None, fpreproc=None, as_pandas=True, verbose_eval=None,show_stdv=True, seed=0, callbacks=None, shuffle=True)

import xgboost as xgb
#为了便捷，使用全数据
dfull = xgb.DMatrix(X,y)
#设定参数
param1 = {'silent':True,'obj':'reg:linear',"gamma":0}
num_round = 100
n_fold=5 #sklearn - KFold
#使用类xgb.cv
time0 = time()
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:610364
#看看类xgb.cv生成了什么结果？
cvresult1 #随着树不断增加，我们的模型的效果如何变化

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,101),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,101),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.legend()
plt.show()

#从这个图中，我们可以看出什么？
#怎样从图中观察模型的泛化能力？
#从这个图的角度来说，模型的调参目标是什么？

为了使用xgboost.cv，我们必须要熟悉xgboost自带的模型评估指标。xgboost在建库的时候本着大而全的目标，和sklearn类似，包括了大约20个模型评估指标，然而用于回归和分类的其实只有几个，大部分是用于一些更加高级的功能比如ranking。来看用于回归和分类的评估指标都有哪些：

param1 = {'silent':True,'obj':'reg:linear',"gamma":0,"eval_metric":"mae"}
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,101),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,101),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.legend()
plt.show()

来看看如果我们调整$\gamma$，会发生怎样的变化：

param1 = {'silent':True,'obj':'reg:linear',"gamma":0}
param2 = {'silent':True,'obj':'reg:linear',"gamma":20}
num_round = 180
n_fold=5

time0 = time()
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:01:083104
time0 = time()
cvresult2 = xgb.cv(param2, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:01:359378

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,181),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,181),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.plot(range(1,181),cvresult2.iloc[:,0],c="green",label="train,gamma=20")
plt.plot(range(1,181),cvresult2.iloc[:,2],c="blue",label="test,gamma=20")
plt.legend()
plt.show()

#从这里，你看出gamma是如何控制过拟合了吗？控制训练集上的训练 - 降低训练集上的表现

试一个分类的例子：

import xgboost as xgb
import matplotlib.pyplot as plt
from time import time
import datetime
from sklearn.datasets import load_breast_cancer
data2 = load_breast_cancer()

x2 = data2.data
y2 = data2.target

dfull2 = xgb.DMatrix(x2,y2)

param1 = {'silent':True,'obj':'binary:logistic',"gamma":0,"nfold":5
          ,"eval_metrics":"error"
         }
param2 = {'silent':True,'obj':'binary:logistic',"gamma":1,"nfold":5}
num_round = 100
time0 = time()
cvresult1 = xgb.cv(param1, dfull2, num_round,metrics=("error"))
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:271581
time0 = time()
cvresult2 = xgb.cv(param2, dfull2, num_round,metrics=("error")) 
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:443810
plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,101),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,101),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.plot(range(1,101),cvresult2.iloc[:,0],c="green",label="train,gamma=1")
plt.plot(range(1,101),cvresult2.iloc[:,2],c="blue",label="test,gamma=1")
plt.legend()
plt.show()

有了xgboost.cv这个工具，我们的参数调整就容易多了。这个工具可以让我们直接看到参数如何影响了模型的泛化能力。接下来，我们将重点讲解如何使用xgboost.cv这个类进行参数调整。到这里，所有关于XGBoost目标函数的原理就讲解完毕了，这个目标函数及这个目标函数所衍生出来的各种数学过程是XGB原理的重中之重，大部分XGB中基于原理的参数都集中在这个模块之中，到这里大家应该已经基本掌握。