集成学习—GBDT原理理解

GBDT，全称为梯度提升决策树，即Gradient Boosting Decision Tree。上一篇博客：集成学习—GBDT（论文研读）中已经对GBDT原始论文——Jerome H. Friedman的《Greedy function approximation: A gradient boosting machine》进行了翻译和精读，这篇博客主要是对GBDT的原理再啃啃，梳理一下，原理理解主要在1-5中体现。

总结

对于回归问题，给定NNN个带标签的训练集{(x1,y1),…,(xN,yN)}\{(x_1,y_1),\dots,(x_N,y_N)\}{(x1,y1),…,(xN,yN)}，这里先给出基于损失函数LLL的不带正则化的GBDT算法：

初始化F0(x)=arg⁡min⁡γ∑i=1NL(yi,γ)F_0(x)=\arg\min_\gamma\sum_{i=1}^NL(y_i,\gamma)F0(x)=argminγ∑i=1NL(yi,γ)
对m=1m=1m=1到MMM循环：
(1) 对i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N，计算：rmi=−[∂L(yi,F(xi))∂F(xi)]F=Fm−1r_{m}^i=-\bigg[{\partial L(y_i,F(x_i))\over \partial F(x_i)}\bigg]_{F=F_{m-1}}rmi=−[∂F(xi)∂L(yi,F(xi))]F=Fm−1(2) 对{xi,rmi}1N\{x_i,r_{m}^i\}_1^N{xi,rmi}1N拟合一个回归树，得到第mmm颗树的叶结点区域Rmj,j=1,2,…,JmR_m^{j},\ j=1,2,\dots,J_mRmj, j=1,2,…,Jm
(3) 对j=1,2,…,Jmj=1,2,\dots,J_mj=1,2,…,Jm，计算：γmj=arg⁡min⁡γ∑xi∈RmjL(yi,Fm−1(xi)+γ)\gamma_m^{j}=\arg\min_\gamma\sum_{x_i\in R_m^{j}}L(y_i,F_{m-1}(x_i)+\gamma)γmj=argγminxi∈Rmj∑L(yi,Fm−1(xi)+γ)(4)更新Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)F_m(x)=F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^{j}I(x\in R_m^{j})Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

这个算法也称为“多重加性回归树”（multiple additive regression trees），简称MART。之所以这样叫，我自己的理解是最后我们有这个算法得到的输出F^(x)\hat F(x)F^(x)，实际上写开了就是F^(x)=FM(x)=∑m=1M∑j=1JmγmjI(x∈Rmj)\hat F(x)=F_M(x)=\sum_{m=1}^M\sum_{j=1}^{J_m}\gamma_m^{j}I(x\in R_m^{j})F^(x)=FM(x)=∑m=1M∑j=1JmγmjI(x∈Rmj)，是一个多重的加性树的结果（一般说来，我们都把每棵树的大小设置为相同的即Jm=JJ_m=JJm=J）。

总的来说，GBDT就是把决策树作为基学习器，首先初始化使得最小化损失函数（损失函数可以由多种选择，例如最小二乘、最小绝对偏差、负二项对数似然、Huber损失等），一般说来这种初始化都很粗糙；然后利用前向分步算法，在每次迭代的时候用一个基学习器去拟合整个函数的最速下降方向，再求最小化当前损失的线性搜索步长，得到的基学习器和步长放入加性模型中就得到这次迭代的模型；这样经过MMM次迭代后最后求和出来的MMM个基学习器的组合就是最后的模型。此外，为了抵抗过拟合，对每次迭代的项还增加了一个收缩系数作为正则化参数，它和迭代次数一起调节最终模型的性能。

由于GBDT也是boosting家族中重要的一员，它是串行生成的，它的理解不同于Adaboost的M个基学习器加权投票，GBDT就像是一个优化算法，从初始化的函数开始出发，每步朝着损失函数下降最快的方向——也就是梯度方向优化，即添加一个函数并走相应的步长，最后走了M步才走到最优点，也就是离目标函数最近的地方，从而得到目标函数的逼近。

下面就按照这个思路，将GBDT一点点理解吃透。第1、2、3、4讲的是GBDT的前两个字母“GB”，后面第5、6讲的是后两个字母“DT”，第7讲的是过抵抗拟合的手段，第8讲利用GBDT可以得到特征的相对重要性，这个可以用来做特征选择或者解释整个模型，以及第9的GBDT的优缺点。

1. 最速下降

在理解GBDT之前，首先要理解下最速下降法（Steepest Descent），在欧式范数下，最速下降法和梯度下降法等价。最速下降法以负梯度的方向作为极小化算法的下降方向，是无约束最优化中最简单的方法。原理如下：

设函数f(x)f(x)f(x)在xkx_kxk附近连续可微，且梯度∇f(xk)≠0\nabla f(x_k)\neq0∇f(xk)=0，由泰勒展式把f(x)f(x)f(x)展开到一阶f(x)≈f(xk)+(x−xk)T∇f(xk)f(x)\approx f(x_k)+(x-x_k)^T\nabla f(x_k)f(x)≈f(xk)+(x−xk)T∇f(xk)可知，若记x−xk=αkdkx-x_k=\alpha_k d_kx−xk=αkdk，其中dkd_kdk为单位向量表示方向，αk>0\alpha_k>0αk>0为矢量表示线性搜索步长，我们希望f(x)f(x)f(x)的值变小，则满足dkT∇f(xk)<0d_k^T\nabla f(x_k)<0dkT∇f(xk)<0的那个方向dkd_kdk就是下降方向。当αk\alpha_kαk取定后，dkT∇f(xk)d_k^T\nabla f(x_k)dkT∇f(xk)的值越小，即−dkT∇f(xk)-d_k^T\nabla f(x_k)−dkT∇f(xk)的值越大，函数下降得越快。由Cauchy-Schwartz不等式∣dkT∇f(xk)∣⩽∥dk∥∥∇f(xk)∥|d_k^T\nabla f(x_k)|\leqslant \|d_k\| \|\nabla f(x_k)\|∣dkT∇f(xk)∣⩽∥dk∥∥∇f(xk)∥故当且仅当dk=−∇f(xk)d_k=-\nabla f(x_k)dk=−∇f(xk)时dkT∇f(xk)d_k^T\nabla f(x_k)dkT∇f(xk)的值越小，函数下降最快。因此，称负梯度方向为最速下降方向。
再求线性搜索步长αk\alpha_kαkαk=arg⁡min⁡αf(xk+αdk)\alpha_k=\arg\min_\alpha f(x_k+\alpha d_k)αk=argαminf(xk+αdk)这样就得到了一次最速下降步的更新：xk+1=xk+αkdk=xk−αk∇f(xk)x_{k+1}=x_k+\alpha_k d_k=x_k-\alpha_k \nabla f(x_k)xk+1=xk+αkdk=xk−αk∇f(xk)

下面就将最速下降的思想应用到boosting上。

2. 加性模型

设我们的损失函数为LLL，我们用一个加性模型FM(x)=∑m=1Mβmh(x;am)F_M(x)=\sum_{m=1}^M\beta_mh(x;a_m)FM(x)=∑m=1Mβmh(x;am)去逼近真实的目标函数，也就是求解min⁡FM(x)∑i=1NL(yi,FM(xi))\min_{F_M(x)}\sum_{i=1}^NL\big(y_i,F_M(x_i)\big)FM(x)mini=1∑NL(yi,FM(xi))这里这样假设加性模型的形式βmh(x;am)\beta_mh(x;a_m)βmh(x;am)是因为它适用于决策树、神经网络、径向基函数、支持向量机等一系列方法。对于求解上式，采取前向分步算法，在每次迭代贪婪地选择βm\beta_mβm和ama_mam来最小化当前损失，即对于m=1,2,…,Mm=1,2,\dots,Mm=1,2,…,M(βm,am)=arg⁡min⁡β,a∑i=1NL(yi,Fm−1(xi)+βh(xi;am))(\beta_m,a_m)=\arg\min_{\beta,a}\sum_{i=1}^NL\big(y_i,F_{m-1}(x_i)+\beta h(x_i;a_m)\big)(βm,am)=argβ,amini=1∑NL(yi,Fm−1(xi)+βh(xi;am))

将βmh(x;am)\beta_mh(x;a_m)βmh(x;am)记为fm(x)f_m(x)fm(x)，则FM(x)=∑m=1Mfm(x)F_M(x)=\sum_{m=1}^Mf_m(x)FM(x)=∑m=1Mfm(x)。利用最速下降法在第mmm次迭代时可以求得fm(x)=αmdmf_m(x)=\alpha_m d_mfm(x)=αmdm，其中N维数据空间中的最佳下降方向dmd_mdm为损失函数的负梯度方向：dm=−{∂∑i=1NL(yi,F(xi))∂F(xi)∣F=Fm−1(x)}d_m=-\bigg\{{\partial \sum_{i=1}^NL\big(y_i,F(x_i)\big)\over \partial F(x_i) }\bigg|_{F=F_{m-1}(x)}\bigg\}dm=−{∂F(xi)∂∑i=1NL(yi,F(xi))∣∣∣∣F=Fm−1(x)}线性搜索步长αm\alpha_mαm为：αm=arg⁡min⁡α∑i=1NL(yi,Fm−1(xi)+αdm)\alpha_m=\arg\min_\alpha \sum_{i=1}^NL\big(y_i,F_{m-1}(x_i)+\alpha d_m\big)αm=argαmini=1∑NL(yi,Fm−1(xi)+αdm)

3. 用基学习器拟合最速下降方向

最佳下降方向dmd_mdm是一个NNN维的向量(dm1,dm2,…,dmN)\big(d_m^1,d_m^2,\dots,d_m^N\big)(dm1,dm2,…,dmN)，其中dmi=−[∂L(yi,F(xi))∂F(xi)]F=Fm−1(x)d_m^{\ i}=-\bigg[{\partial L\big(y_i,F(x_i)\big)\over \partial F(x_i) }\bigg]_{F=F_{m-1}(x)}dm i=−[∂F(xi)∂L(yi,F(xi))]F=Fm−1(x)它不能表示在N维数据空间中的函数，因此GBDT的最奥妙的地方就在于每次迭代中，用基学习器βmh(x;am)\beta_m h(x;a_m)βmh(x;am)去拟合当前的最佳下降方向dmd_mdm，也就是把这个NNN维的向量当做新的“伪标签”，用最简单的最小二乘去拟合得到当前的一个基学习器，即最小化平方损失函数得到am=arg⁡min⁡β,a∑i=1N[dmi−βh(xi;a)]2a_m=\arg\min_{\beta,a} \sum_{i=1}^N\big[d_m^{\ i}-\beta h(x_i;a)\big]^2am=argβ,amini=1∑N[dm i−βh(xi;a)]2此时的线性搜索为βm=arg⁡min⁡β∑i=1NL(yi,Fm−1(xi)+βh(xi;am))\beta_m=\arg\min_\beta \sum_{i=1}^NL\big(y_i,F_{m-1}(x_i)+\beta h(x_i;a_m)\big)βm=argβmini=1∑NL(yi,Fm−1(xi)+βh(xi;am))这样就得到了第mmm次迭代时的基学习器βmh(x;am)\beta_m h(x;a_m)βmh(x;am)，以及当前对于目标函数的逼近Fm(x)=∑k=1mβkh(x;ak)F_m(x)=\sum_{k=1}^m\beta_k h(x;a_k)Fm(x)=∑k=1mβkh(x;ak)

4. Gradient Boosting框架

这样我们就得到了下面基本的Gradient Boosting框架：

初始化F0(x)=arg⁡min⁡c∑i=1NL(yi,c)F_0(x)=\arg\min_c\sum_{i=1}^NL(y_i,c)F0(x)=argminc∑i=1NL(yi,c)
对m=1m=1m=1到MMM：
(1) 对i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N，计算第mmm次的最速下降方向：dmi=−[∂L(yi,F(xi))∂F(xi)]F=Fm−1(x)d_{m}^{\ i}=-\bigg[{\partial L(y_i,F(x_i))\over \partial F(x_i)}\bigg]_{F=F_{m-1}(x)}dm i=−[∂F(xi)∂L(yi,F(xi))]F=Fm−1(x)(2) 用基学习器拟合下降方向，得到第mmm个基学习器的参数am=arg⁡min⁡β,a∑i=1N[dmi−βh(xi;a)]2a_m=\arg\min_{\beta,a} \sum_{i=1}^N\big[d_m^{\ i}-\beta h(x_i;a)\big]^2am=argβ,amini=1∑N[dm i−βh(xi;a)]2(3) 计算线性搜索：βm=arg⁡min⁡β∑i=1NL(yi,Fm−1(xi)+βh(xi;am))\beta_m=\arg\min_\beta \sum_{i=1}^NL\big(y_i,F_{m-1}(x_i)+\beta h(x_i;a_m)\big)βm=argβmini=1∑NL(yi,Fm−1(xi)+βh(xi;am))(4)更新Fm(x)=Fm−1(x)+βmh(x;am)F_m(x)=F_{m-1}(x)+\beta_mh(x;a_m)Fm(x)=Fm−1(x)+βmh(x;am)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

5. 决策树的形式

之前在集成学习—决策树（CART）中曾经详细讨论过CART树 (Classification And Regression Tree)，这里主要使用的是回归树，GBDT的第mmm个基学习器是一个JmJ_mJm个叶结点的回归树，于是相当于这棵树把特征空间划分为JmJ_mJm个互不相交的单元Rm1,Rm2⋯,RmJmR_m^1,R_m^2\cdots,R_m^{J_m}Rm1,Rm2⋯,RmJm，并且在每个单元RmjR_m^jRmj上有一个固定的输出值bmjb_m^jbmj，于是第mmm个基学习器可以表示为：
h(x;{bmj,Rmj}j=1Jm)=∑j=1JmbmjI(x∈Rmj)h(x;\{b_m^j,R_m^j\}_{j=1}^{J_m})=\sum_{j=1}^{J_m}b_m^j I(x\in R_m^j) h(x;{bmj,Rmj}j=1Jm)=j=1∑JmbmjI(x∈Rmj)这样，上述Gradient Boosting框架中用基学习器拟合下降方向就相当于是用回归树来拟合下降方向（伪标签）{dmi}i=1N\{d_{m}^{\ i}\}_{i=1}^N{dm i}i=1N，得到JmJ_mJm个单元Rm1,Rm2⋯,RmJmR_m^1,R_m^2\cdots,R_m^{J_m}Rm1,Rm2⋯,RmJm以及在每个单元上的系数bm1,bm2,…,bmJmb_m^1,b_m^2,\dots,b_m^{J_m}bm1,bm2,…,bmJm。

对于如何划分特征空间得到这样JmJ_mJm个单元，是根据MSE准则寻找最优切分变量和切分点，具体形式我们不关心，但我们关心每个单元上的系数的具体形式。CART回归树就是用MSE最小的准则来求解每个单元上的最优输出值，即在任意单元RmjR_m^jRmj中最小化∑xi∈Rmj(dmi−bmj)2\sum_{x_i\in R_m^j}(d_{m}^{\ i}-b_m^j)^2xi∈Rmj∑(dm i−bmj)2令它对bmjb_m^jbmj求导为零，可知最优的bmjb_m^jbmj是RmjR_m^jRmj上所有特征xix_ixi对应的伪标签dmid_{m}^{\ i}dm i的均值，即bmj=avexi∈Rmj(dmi)b_m^j=ave_{x_i\in R_m^j}(d_{m}^{\ i})bmj=avexi∈Rmj(dm i)

因此，我们可以把上述内容带入到Gradient Boosting框架，就得到了更新Fm(x)=Fm−1(x)+βm∑j=1JmbmjI(x∈Rmj)≜Fm−1(x)+∑j=1JmγmjI(x∈Rmj)\begin{aligned} F_m(x)&=F_{m-1}(x)+\beta_m\sum_{j=1}^{J_m}b_m^j I(x\in R_m^j)\\ &\triangleq F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^j I(x\in R_m^j) \end{aligned}Fm(x)=Fm−1(x)+βmj=1∑JmbmjI(x∈Rmj)≜Fm−1(x)+j=1∑JmγmjI(x∈Rmj)其中γmj=βmbmj\gamma_m^j =\beta_mb_m^jγmj=βmbmj。因此我们线性搜索的目标变为
{γmj}j=1J=arg⁡min⁡{γj}1J∑i=1NL(yi,Fm−1(xi)+∑j=1JmγjI(xi∈Rmj))\{\gamma_m^j \}_{j=1}^J=\arg\min_{\{\gamma^{j}\}_1^J}\sum_{i=1}^NL\bigg(y_i,F_{m-1}(x_i)+\sum_{j=1}^{J_m}\gamma^j I(x_i\in R_m^j)\bigg) {γmj}j=1J=arg{γj}1Jmini=1∑NL(yi,Fm−1(xi)+j=1∑JmγjI(xi∈Rmj))由于回归树产生的区域具有不相交的性质，这就简化为γmj=arg⁡min⁡γ∑xi∈RmjL(yi,Fm−1(xi)+γ)\gamma_{m}^j=\arg\min_{\gamma}\sum_{x_i\in R_{m}^j}L\bigg(y_i,F_{m-1}(x_i)+\gamma \bigg)γmj=argγminxi∈Rmj∑L(yi,Fm−1(xi)+γ)

至此，就得到了最开始总结里的GBDT的框架。

6. 损失函数（LS/LAD/M/logistic）

下面具体将不同的损失函数带入GBDT的框架：

LS_TreeBoost

LS即最小二乘（Least-squares），这里的损失函数为平方损失L(y,F)=(y−F)2/2L(y,F)=(y-F)^2/2L(y,F)=(y−F)2/2，因此初始化为数据标签的均值，即F0(x)=yˉF_0(x)=\bar yF0(x)=yˉ，最速下降方向的第iii个分量为dmi=yi−Fm−1(xi)d_{m}^{\ i}=y_i-F_{m-1}(x_i)dm i=yi−Fm−1(xi)，从形式上看，这个量是数据标签和当前拟合值的残差（residuals），因此也有很多人说GBDT是用基学习器去拟合当前的残差，所以最开始给出的最速下降方向是以rmir_{m}^irmi而不是dmid_{m}^{\ i}dm i的形式（对此，我个人觉得这个是比较狭义的看法，如果换一个损失函数，不一定就是这么标准的残差形式了，所以宁愿称之为最速下降方向）。由于CART回归树也是使用最小化平方损失，因此在用残差拟合回归树时，除了可以得到划分单元RmjR_m^{j}Rmj外，还能得到与每个单元上的系数bmjb_m^{j}bmj相等的线性搜索。

最后得到LS_TreeBoost算法为：

初始化F0(x)=yˉF_0(x)=\bar yF0(x)=yˉ
对m=1m=1m=1到MMM循环：
(1) rmi=yi−Fm−1(xi)r_{m}^i=y_i-F_{m-1}(x_i)rmi=yi−Fm−1(xi)，i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N
(2) 对{xi,rmi}1N\{x_i,r_{m}^i\}_1^N{xi,rmi}1N拟合一个回归树，得到第mmm颗树的参数(γmj,Rmj)=arg⁡min⁡Rj,γj∑xi∈Rj(rmi−γj)2,j=1,2,…,Jm(\gamma_m^{j},R_m^j)=\arg\min_{R^j,\gamma^j}\sum_{x_i\in R^j}(r_{m}^{\ i}-\gamma^j)^2,\ \ j=1,2,\dots,J_m(γmj,Rmj)=argRj,γjminxi∈Rj∑(rm i−γj)2, j=1,2,…,Jm (3) 更新Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)F_m(x)=F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^j I(x\in R_m^j)Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

LAD_TreeBoost

LAD即最小绝对偏差（Least absolute deviation），这里的损失函数为绝对偏差损失L(y,F)=∣y−F∣L(y,F)=|y-F|L(y,F)=∣y−F∣，因此初始化为数据标签的中位数，即F0(x)=arg⁡min⁡γ∑i=1N∣yi−γ∣=median{yi}1NF_0(x)=\arg\min_\gamma\sum_{i=1}^N|y_i-\gamma|=median\{y_i\}_1^NF0(x)=argγmini=1∑N∣yi−γ∣=median{yi}1N

最速下降方向的第iii个分量为dmi=−[∂L(yi,F(xi))∂F(xi)]F=Fm−1=sign(yi−Fm−1(xi))d_{m}^{\ i}=-\bigg[{\partial L(y_i,F(x_i))\over \partial F(x_i)}\bigg]_{F=F_{m-1}}=sign\big(y_i-F_{m-1}(x_i)\big)dm i=−[∂F(xi)∂L(yi,F(xi))]F=Fm−1=sign(yi−Fm−1(xi))也就是说基学习器需要去拟合当前残差的符号函数。在得到划分单元后，线性搜索为
γmj=arg⁡min⁡γ∑xi∈RmjL(yi,Fm−1(xi)+γ)=arg⁡min⁡γ∑xi∈Rmj∣yi−Fm−1(xi)−γ∣=medianxi∈Rmj{yi−Fm−1(xi)}1N\begin{aligned} \gamma_m^{j}&=\arg\min_\gamma\sum_{x_i\in R_m^{j}}L(y_i,F_{m-1}(x_i)+\gamma)\\ &=\arg\min_\gamma\sum_{x_i\in R_m^{j}}|y_i-F_{m-1}(x_i)-\gamma|\\ &=median_{x_i\in R_{m}^j}\{y_i-F_{m-1}(x_i)\}_1^N \end{aligned}γmj=argγminxi∈Rmj∑L(yi,Fm−1(xi)+γ)=argγminxi∈Rmj∑∣yi−Fm−1(xi)−γ∣=medianxi∈Rmj{yi−Fm−1(xi)}1N最后得到LAD_TreeBoost算法为：

初始化F0(x)=median{yi}1NF_0(x)=median\{y_i\}_1^NF0(x)=median{yi}1N
对m=1m=1m=1到MMM循环：
(1) rmi=sign(yi−Fm−1(xi))r_{m}^i=sign\big(y_i-F_{m-1}(x_i)\big)rmi=sign(yi−Fm−1(xi))，i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N
(2) 对{xi,rmi}1N\{x_i,r_{m}^i\}_1^N{xi,rmi}1N拟合回归树，得到第mmm颗树的划分单元{Rmj}j=1Jm\{R_m^j\}_{j=1}^{J_m}{Rmj}j=1Jm
(3) 线性搜索为γmj=medianxi∈Rmj{yi−Fm−1(xi)}1N\gamma_m^{j}=median_{x_i\in R_{m}^j}\{y_i-F_{m-1}(x_i)\}_1^Nγmj=medianxi∈Rmj{yi−Fm−1(xi)}1N
(4) 更新Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)F_m(x)=F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^j I(x\in R_m^j)Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

M_TreeBoost

这里的损失函数为Huber损失函数L(y,F)={12(y−F)2,∣y−F∣⩽δδ(∣y−F∣−δ/2),∣y−F∣>δL(y,F)=\begin{cases}{1\over2}(y-F)^2, &|y-F|\leqslant \delta\\\delta(|y-F|-\delta/2), &|y-F|>\delta \end{cases}L(y,F)={21(y−F)2,δ(∣y−F∣−δ/2),∣y−F∣⩽δ∣y−F∣>δ由于这里是和绝对值项∣y−F∣|y-F|∣y−F∣有关，因此初始化可以设为数据标签的中位数，即F0(x)=median{yi}1NF_0(x)=median\{y_i\}_1^NF0(x)=median{yi}1N

最速下降方向的第iii个分量（也就是第iii个伪标签的值）为dmi=−[∂L(yi,F(xi))∂F(xi)]F(x)=Fm−1(x)={yi−Fm−1(xi),∣yi−Fm−1(xi)∣⩽δδ⋅sign(yi−Fm−1(xi)),∣yi−Fm−1(xi)∣>δ\begin{aligned} d_{m}^{\ i}&=-\bigg[{\partial L\big(y_i,F(x_i)\big)\over\partial F(x_i)}\bigg]_{F(x)=F_{m-1}(x)}\\ &=\begin{cases} y_i-F_{m-1}(x_i), & |y_i-F_{m-1}(x_i)|\leqslant\delta\\ \delta\cdot sign(y_i-F_{m-1}(x_i)), & |y_i-F_{m-1}(x_i)|>\delta \end{cases} \end{aligned}dm i=−[∂F(xi)∂L(yi,F(xi))]F(x)=Fm−1(x)={yi−Fm−1(xi),δ⋅sign(yi−Fm−1(xi)),∣yi−Fm−1(xi)∣⩽δ∣yi−Fm−1(xi)∣>δ可以看出，δ\deltaδ的值受绝对损失控制，它的最优值取决于y−F∗(x)y-F^*(x)y−F∗(x)的分布，其中F∗F^*F∗是真正的目标函数。选择δ\deltaδ值的常见的做法是选为分布∣y−F∗(x)∣|y-F^*(x)|∣y−F∗(x)∣的α\alphaα分位数，这里(1−α)(1-\alpha)(1−α)把样本分为一部分，这部分可以任意修改而不会降低结果的质量。但由于F∗(x)F^*(x)F∗(x)是未知的，我们在第mmm次迭代中使用当前估计Fm−1(x)F_{m-1}(x)Fm−1(x)作为近似。∣y−Fm−1(x)∣|y-F_{m-1}(x)|∣y−Fm−1(x)∣的分布由当前的残差来估计，这样δm=quantileα{∣yi−Fm−1(xi)∣}1N\delta_m=quantile_\alpha\{|y_i-F_{m-1}(x_i)|\}_1^Nδm=quantileα{∣yi−Fm−1(xi)∣}1N

因此，基学习器需要去拟合的伪标签也是和当前残差y−Fm−1(x)y-F_{m-1}(x)y−Fm−1(x)有关。使用回归树作为基学习器，在得到划分单元RmjR_{m}^jRmj后，将Huber损失函数带入线性搜索得
γmj=arg⁡min⁡γ∑xi∈RmjL(yi,Fm−1(xi)+γ)=arg⁡min⁡γ∑xi∈Rmj{12(yi−Fm−1(xi)−γ)2,∣yi−Fm−1(xi)−γ∣⩽δmδm(∣yi−Fm−1(xi)−γ∣−δm/2),∣yi−Fm−1(xi)−γ∣>δm\begin{aligned} \gamma_m^{j}&=\arg\min_\gamma\sum_{x_i\in R_m^{j}}L(y_i,F_{m-1}(x_i)+\gamma)\\ &=\arg\min_\gamma\sum_{x_i\in R_m^{j}} \begin{cases} {1\over2}(y_i-F_{m-1}(x_i)-\gamma)^2,& |y_i-F_{m-1}(x_i)-\gamma|\leqslant\delta_m\\ \delta_m(|y_i-F_{m-1}(x_i)-\gamma|-\delta_m/2),& |y_i-F_{m-1}(x_i)-\gamma|>\delta_m\\ \end{cases}\\ \end{aligned}γmj=argγminxi∈Rmj∑L(yi,Fm−1(xi)+γ)=argγminxi∈Rmj∑{21(yi−Fm−1(xi)−γ)2,δm(∣yi−Fm−1(xi)−γ∣−δm/2),∣yi−Fm−1(xi)−γ∣⩽δm∣yi−Fm−1(xi)−γ∣>δm

因为平方最小求出来大概形如avexi∈Rmj{yi−Fm−1(xi)}1Nave_{x_i\in R_{m}^j}\{y_i-F_{m-1}(x_i)\}_1^Navexi∈Rmj{yi−Fm−1(xi)}1N，绝对值最小求出来大概形如medianxi∈Rmj{yi−Fm−1(xi)}1Nmedian_{x_i\in R_{m}^j}\{y_i-F_{m-1}(x_i)\}_1^Nmedianxi∈Rmj{yi−Fm−1(xi)}1N，但由于后面带有一个约束条件，所以不这么简单。因此，令r~mj=medianxi∈Rmj{rm−1(xi)}1N\tilde{r}_m^j=median_{x_i\in R_{m}^j}\{r_{m-1}(x_i)\}_1^Nr~mj=medianxi∈Rmj{rm−1(xi)}1N，其中rm−1(xi)=yi−Fm−1(xi)r_{m-1}(x_i)=y_i-F_{m-1}(x_i)rm−1(xi)=yi−Fm−1(xi)是当前残差，则上面复杂的线性搜索γmj\gamma_m^{j}γmj可以从中位数r~mj\tilde{r}_m^jr~mj开始优化近似得到γmj=r~mj+1Nmj∑xi∈Rmjsign(rm−1(xi)−r~mj)⋅min⁡(δm,abs(rm−1(xi)−r~mj))\gamma_m^{j}=\tilde{r}_m^{j}+{1\over N_m^{j}}\sum_{x_i\in R_m^{j}}sign\big(r_{m-1}(x_i)-\tilde{r}_m^{j}\big)\cdot \min\big(\delta_m,abs(r_{m-1}(x_i)-\tilde{r}_m^{j})\big)γmj=r~mj+Nmj1xi∈Rmj∑sign(rm−1(xi)−r~mj)⋅min(δm,abs(rm−1(xi)−r~mj))

最后得到M_TreeBoost算法为：

初始化F0(x)=median{yi}1NF_0(x)=median\{y_i\}_1^NF0(x)=median{yi}1N
对m=1m=1m=1到MMM循环：
(1) rm(xi)=yi−Fm−1(xi)r_{m}(x_i)=y_i-F_{m-1}(x_i)rm(xi)=yi−Fm−1(xi)，i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N
(2) δm=quantileα{∣rm(xi)∣}1N\delta_m=quantile_\alpha\{|r_{m}(x_i)|\}_1^Nδm=quantileα{∣rm(xi)∣}1N
(3) i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,Nrmi={rm(xi),∣rm(xi)∣⩽δmδm⋅sign(rm(xi)),∣rm(xi)∣>δm\begin{aligned} r_{m}^{\ i} &=\begin{cases} r_{m}(x_i), & |r_{m}(x_i)|\leqslant\delta_m\\ \delta_m\cdot sign(r_{m}(x_i)), & |r_{m}(x_i)|>\delta_m \end{cases} \end{aligned}rm i={rm(xi),δm⋅sign(rm(xi)),∣rm(xi)∣⩽δm∣rm(xi)∣>δm (4) 对{xi,rmi}1N\{x_i,r_{m}^i\}_1^N{xi,rmi}1N拟合回归树，得到第mmm颗树的划分单元{Rmj}j=1Jm\{R_m^j\}_{j=1}^{J_m}{Rmj}j=1Jm
(5) r~mj=medianxi∈Rmj{rm−1(xi)}1N,j=1,2,…,Jm\tilde{r}_m^j=median_{x_i\in R_{m}^j}\{r_{m-1}(x_i)\}_1^N,\ \ j=1,2,\dots,J_mr~mj=medianxi∈Rmj{rm−1(xi)}1N, j=1,2,…,Jm
(6) 线性搜索为γmj=r~mj+1Nmj∑xi∈Rmjsign(rm−1(xi)−r~mj)⋅min⁡(δm,abs(rm−1(xi)−r~mj)),j=1,2,…,Jm\gamma_m^{j}=\tilde{r}_m^{j}+{1\over N_m^{j}}\sum_{x_i\in R_m^{j}}sign\big(r_{m-1}(x_i)-\tilde{r}_m^{j}\big)\cdot \min\big(\delta_m,abs(r_{m-1}(x_i)-\tilde{r}_m^{j})\big),\ \ j=1,2,\dots,J_mγmj=r~mj+Nmj1xi∈Rmj∑sign(rm−1(xi)−r~mj)⋅min(δm,abs(rm−1(xi)−r~mj)), j=1,2,…,Jm(7) 更新Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)F_m(x)=F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^j I(x\in R_m^j)Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

上述三种不同损失的比较

在正态误差下LS_TreeBoost（平方损失）表现较好；
在异常值较多的情况下LAD_TreeBoost（绝对损失）表现较好；
M_TreeBoost（Huber损失）在这两种情况下都接近最好表现。

因此，在实际应用中，如果数据的清洁度不能保证或数据存在较多离群值，那么M_TreeBoost的相对较高的精度、一致的性能和鲁棒性可能是一个很大的优势。

LkL_kLk_TreeBoost

这里的损失函数为负的二项对数似然（FHT00）L(y,F)=log(1+exp(−2yF)),y∈{−1,1}L(y,F)=log(1+exp(-2yF)),\ \ y\in\{-1,1\}L(y,F)=log(1+exp(−2yF)), y∈{−1,1}其中F(x)=12log[Pr(y=1∣x)Pr(y=−1∣x)]F(x)={1\over2}log\bigg[{Pr(y=1|x)\over Pr(y=-1|x)}\bigg]F(x)=21log[Pr(y=−1∣x)Pr(y=1∣x)]因此初始化为F0(x)=12log1+yˉ1−yˉF_0(x)={1\over2}log{1+\bar y\over 1-\bar y}F0(x)=21log1−yˉ1+yˉ

第mmm次迭代的最速下降方向的第iii个分量（也就是第mmm次迭代的第iii个伪标签的值）为dmi=−[∂L(yi,F(xi))∂F(xi)]F(x)=Fm−1(x)=2yi/(1+exp(2yiFm−1(x)))d_{m}^{\ i}=-\bigg[{\partial L(y_i,F(x_i))\over \partial F(x_i)}\bigg]_{F(x)=F_{m-1}(x)}=2y_i\big/\bigg(1+exp\big(2y_iF_{m-1}(x)\big)\bigg)dm i=−[∂F(xi)∂L(yi,F(xi))]F(x)=Fm−1(x)=2yi/(1+exp(2yiFm−1(x)))用回归树作为基学习器，在得到划分单元RmjR_{m}^jRmj后，将负二项对数似然损失函数带入线性搜索得：γmj=arg⁡min⁡γ∑xi∈Rmjlog(1+exp(−2yi(Fm−1(xi)+γ)))\gamma_{m}^j=\arg\min_\gamma\sum_{x_i\in R_{m}^j}log\bigg(1+exp\big(-2y_i(F_{m-1}(x_i)+\gamma)\big)\bigg)γmj=argγminxi∈Rmj∑log(1+exp(−2yi(Fm−1(xi)+γ)))上式没有闭式解，因此用牛顿迭代法来求解得到γmj=∑xi∈Rmjdmi/∑xi∈Rmj∣dmi∣(2−∣dmi∣)\gamma_{m}^j=\sum_{x_i\in R_{m}^j}d_{m}^{\ i}\bigg/\sum_{x_i\in R_{m}^j}|d_{m}^{\ i}|(2-|d_{m}^{\ i}|)γmj=xi∈Rmj∑dm i/xi∈Rmj∑∣dm i∣(2−∣dm i∣)

最后得到LkL_kLk_TreeBoost算法为：

初始化F0(x)=12log1+yˉ1−yˉF_0(x)={1\over2}log{1+\bar y\over 1-\bar y}F0(x)=21log1−yˉ1+yˉ
对m=1m=1m=1到MMM循环：
(1) rmi=2yi/(1+exp(2yiFm−1(x)))r_{m}^i=2y_i\big/\bigg(1+exp\big(2y_iF_{m-1}(x)\big)\bigg)rmi=2yi/(1+exp(2yiFm−1(x)))，i=1,2,…,Ni=1,2,\dots,Ni=1,2,…,N
(2) 对{xi,rmi}1N\{x_i,r_{m}^i\}_1^N{xi,rmi}1N拟合回归树，得到第mmm颗树的划分单元{Rmj}j=1Jm\{R_m^j\}_{j=1}^{J_m}{Rmj}j=1Jm
(3) 线性搜索为γmj=∑xi∈Rmjrmi/∑xi∈Rmj∣rmi∣(2−∣rmi∣)\gamma_m^{j}=\sum_{x_i\in R_{m}^j}r_{m}^{\ i}\big/\sum_{x_i\in R_{m}^j}|r_{m}^{\ i}|(2-|r_{m}^{\ i}|)γmj=∑xi∈Rmjrm i/∑xi∈Rmj∣rm i∣(2−∣rm i∣)
(4) 更新Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)F_m(x)=F_{m-1}(x)+\sum_{j=1}^{J_m}\gamma_m^j I(x\in R_m^j)Fm(x)=Fm−1(x)+∑j=1JmγmjI(x∈Rmj)
输出F^(x)=FM(x)\hat F(x)=F_M(x)F^(x)=FM(x)

最后得到的FM(x)F_M(x)FM(x)通过F(x)=12log[Pr(y=1∣x)Pr(y=−1∣x)]F(x)={1\over2}log\bigg[{Pr(y=1|x)\over Pr(y=-1|x)}\bigg]F(x)=21log[Pr(y=−1∣x)Pr(y=1∣x)]和对数几率联系起来，这样可以转化得到概率估计p+(x)=Pr^(y=1∣x)=1/(1+e−2FM(x))p−(x)=Pr^(y=−1∣x)=1/(1+e2FM(x))p_+(x)=\hat{Pr}(y=1|x)=1/(1+e^{-2F_M(x)})\\ p_-(x)=\hat{Pr}(y=-1|x)=1/(1+e^{2F_M(x)})p+(x)=Pr^(y=1∣x)=1/(1+e−2FM(x))p−(x)=Pr^(y=−1∣x)=1/(1+e2FM(x))因此GBDT又可以用于分类，预测的标签为y^(x)=2⋅I[c(−1,1)p+(x)>c(1,−1)p−(x)]−1\hat{y}(x)=2\cdot I[c(-1,1)p_+(x)>c(1,-1)p_-(x)]-1y^(x)=2⋅I[c(−1,1)p+(x)>c(1,−1)p−(x)]−1其中c(y^,y)c(\hat y,y)c(y^,y)是当真实值为yyy，预测值为y^\hat yy^时有关的损失。

类似地，将损失函数改为多项分布的对数似然函数，并利用对称的多重逻辑变换就可以得到多分类逻辑回归和分类。这里就不讨论了。

7. 正则化

Shrinkage（收缩）的思想认为，每次走一小步逐渐逼近结果的效果，要比每次迈一大步很快逼近结果的方式更容易避免过拟合。即它不完全信任每一个棵残差树，它认为每棵树只学到了真理的一小部分，累加的时候只累加一小部分，希望通过多学几棵树来弥补不足。

在实际应用中，我们可以通过控制迭代次数来控制过拟合，但这隐含了一种先验信念，即“稀疏”近似的项数越少，预测效果越好。然而，人们经常发现，通过收缩实现的正则化比通过限制项的数量实现的正则化效果更好[Copas(1983)]。因此，Friedman提出在模型更新时添加一个“学习率”参数ν\nuν，GBDT的更新就变为Fm(x)=Fm−1(x)+ν⋅∑j=1JmγmjI(x∈Rmj),0<ν⩽1F_m(x)=F_{m-1}(x)+\nu\cdot\sum_{j=1}^{J_m}\gamma_m^{j}I(x\in R_m^{j}),\ \ \ 0<\nu\leqslant 1Fm(x)=Fm−1(x)+ν⋅j=1∑JmγmjI(x∈Rmj), 0<ν⩽1每次更新都是简单地根据“学习率”参数ν\nuν的值进行缩放。这样就达到了收缩的目的，减少了对单独的每棵树的信任，但一定程度上需要更多的树来达到一个很好的拟合效果。

因此，这里就有两种正则化参数，学习速率ν\nuν和迭代次数MMM。设置一个较小的ν\nuν值会增加MMM的最佳值，也就是增加迭代次数。经验表明，ν⩽0.1\nu\leqslant0.1ν⩽0.1会得到一个更好的泛化性能[Friedman(2002)]，但还应该考虑到计算的可行性，迭代次数MMM应该在计算方便或者可行的时候尽可能大。

8. 特征重要性

特征重要性就是看单个特征xjx_jxj对最后估计出的模型F(x)F(x)F(x)的相对影响IjI_jIj。Breiman、Friedman、Olshen和Stone(1983)提出用下面这个式子来度量特征xjx_jxj在单颗决策树TTT上对模型F(x)F(x)F(x)的相对影响Ij(T)I_j(T)Ij(T)Ij2(T)=∑t=1J−1lt2I(vt=j)I_j^2(T)=\sum_{t=1}^{J-1} l_t^2I(v_t=j)Ij2(T)=t=1∑J−1lt2I(vt=j)由于CART树特殊的生成规则，即只有一个孩子能产生孩子的二叉树结构，因此一个拥有JJJ个叶子结点的决策树TTT有且仅有J−1J-1J−1个内部结点，而分裂只在内部结点上进行，因此这里对t=1t=1t=1到J−1J-1J−1求和；vtv_tvt是和结点ttt有关的分裂变量，I(vt=j)I(v_t=j)I(vt=j)表示分裂变量是xjx_jxj时的示性函数；lt2l_t^2lt2是类似于Gini指数或信息熵的一个度量分裂准则，值越大表示左孩子和右孩子差别越大，因此用当前特征来做分割的影响更大。评估一个当前区域RRR可能分裂为两个子区域(Rl,Rr)(R_l,R_r)(Rl,Rr)的具体形式如下l2(Rl,Rr)=wlwrwl+wr(yˉl−yˉr)2l^2(R_l,R_r)={w_lw_r\over w_l+w_r}(\bar y_l-\bar y_r)^2l2(Rl,Rr)=wl+wrwlwr(yˉl−yˉr)2其中yˉl,yˉr\bar y_l,\bar y_ryˉl,yˉr分别为左、右两边的子标签的均值，wl,wrw_l,w_rwl,wr为对应的权重之和。所以这样的度量Ij2(T)I_j^2(T)Ij2(T)就相当于是把这棵树所有分裂结点取特征xjx_jxj时的影响全都加和起来。

因此，对于一个由Boosting得到的决策树集{Tm}1M\{T_m\}_1^M{Tm}1M，我们可以这个集合中所有树的平均值来定义这个模型中特征xjx_jxj在决策树集{Tm}1M\{T_m\}_1^M{Tm}1M上对模型F(x)F(x)F(x)的相对影响IjI_jIj，即Ij2=1M∑m=1MIj2(Tm)I_j^2={1\over M}\sum_{m=1}^{M} I_j^2(T_m)Ij2=M1m=1∑MIj2(Tm)这样就得到了每个特征对最后估计出的模型F(x)F(x)F(x)的相对影响了，这个值越大，说明对模型影响越大，这个特征也就越重要，因此是特征之间的一个相对重要性度量。我们可以用它来解释模型或者做变量选择等。

9. 优缺点

优点：

预测精度高；
在损失函数的选择上有很大的灵活性；
对特征的单调变换具有不变性，例如，用xj,log⁡xj,exjx_j,\log x_j,e^{x_j}xj,logxj,exj或xjax^a_jxja作为第jjj个输入变量都会得到相同的结果；
可以灵活处理各种类型的数据，包括连续值和离散值；
可以做分类也可以做回归；
继承了决策树可以很好处理缺失值的特性；
对长尾分布和异常值不敏感，使用一些稳健的损失函数，比如Huber损失函数等，对异常值的鲁棒性非常强；
可以作为模型预处理工具，先用较少的计算量快速建模，可以在先前停止的地方重新启动过程，不会造成计算损失。

缺点：

由于弱学习器之间存在依赖关系，难以并行训练数据。不过可以通过自采样的SGBT来达到部分并行；
如果数据维度较高时会加大算法的计算复杂度。

参考资料

Friedman J H . Greedy Function Approximation: A Gradient Boosting Machine[J]. The Annals of Statistics, 2001, 29(5):1189-1232.
Minasny B . Trevor Hastie, Robert Tishirani, Jerome Friedman, The Elements of Statistical Learning Second Edition (2009) Springer Series in Statistics, ISBN 0172-7397, 745 pp.[J]. Geoderma, 2009, 153(1-2):1.