Xgboost光速推导

1.为什么要推导Xgboost?

  Xgboost是今年数据科学竞赛非常常用的boosting大杀器，尤其是与传统bagging型学习器（如RF）做融合时有往往会有奇效。虽然近两年新出了传说中比Xgboost“更强大”的lightBGM选手，但是Xgboost在比赛中的表现依然亮眼，在小体量的数据集上，可以说是和后者不相上下。

  毕竟，如霍元甲所言，天下算法本没有高低贵贱之分。

2.为什么要快速推导Xgboost?

  因为天下武功，唯快不破！

3.开始推导！

  首先介绍两种常见的损失函数，作为推导时理解记忆的辅助：
MSE（均方损失）:L(Θ)=∑in(yi−y^i)2MSE（均方损失）: L(\Theta )=\sum_{i}^{n}(y_{i}-\hat{y}_{i})^{2}MSE（均方损失）:L(Θ)=i∑n​(yi​−y^​i​)2
Logloss（对数损失）：L(Θ)=∑in[yiln(h(Θ))+(1−yi)ln(1−h(Θ))]Logloss（对数损失）：L(\Theta )=\sum_{i}^{n}[y_{i}ln(h_{( \Theta)})+(1-y_{i})ln(1-h_{(\Theta)})]Logloss（对数损失）：L(Θ)=i∑n​[yi​ln(h(Θ)​)+(1−yi​)ln(1−h(Θ)​)]
  由于Xgboost的计算思路是一棵树一棵树的叠加计算，所以在使用损失函数时不需要对其做平均处理（即除以n）。下面进入正题！

====================================================================================
开始计时（600s）

  设Xgboost由T棵树组成，用boosting的思想一颗一颗依次建立；
则有：y^=∑t=1tft(xi),其中ft∈F(F表示所有可能的CART树)\hat {y} = \sum_{t=1}^{t}f_{t}(x_{i}) , 其中f_t \in F (F表示所有可能的CART树)y^​=t=1∑t​ft​(xi​),其中ft​∈F(F表示所有可能的CART树)
  为找到算法模型的结构风险最小化形式，需建立一个目标函数obj(Θ)obj(\Theta)obj(Θ);
有:obj(Θ)=∑inl(yi,y^i)+∑inΩ(ft)(1)obj( \Theta)=\sum_{i}^{n}l(y_{i},\hat {y}_{i})+\sum_{i}^{n}\Omega(f_{t})(1)obj(Θ)=i∑n​l(yi​,y^​i​)+i∑n​Ω(ft​)(1)
  这时可以看出目标函数由损失函数及正则化项组成，由boosting的思想，我们可以通过前t-1棵训练好的树学习器的训练效果，来决定如何产生第t颗树的生成。因此为求得第t棵树的目标函数，可以用迭代的方法进行加法训练：
y^i(0)=0\hat y_{i}^{(0)}=0y^​i(0)​=0
y^i(1)=y^i(0)+f1(xi)\hat y_{i}^{(1)}=\hat y_{i}^{(0)}+f_{1}(x_{i})y^​i(1)​=y^​i(0)​+f1​(xi​)
……
y^i(t)=y^i(t−1)+ft(xi)\hat y_{i}^{(t)}=\hat y_{i}^{(t-1)} + f_{t}(x_{i})y^​i(t)​=y^​i(t−1)​+ft​(xi​)
  代入（1）式中可得：
obj(t)=∑inl(yi,y^i(t−1)+ft(xi))+Ω(ft)+constantobj^{(t)}= \sum_{i} ^ {n}l(y_{i},\hat {y}_{i}^{(t-1)}+f_{t}(x_{i}))+\Omega(f_{t})+constantobj(t)=i∑n​l(yi​,y^​i(t−1)​+ft​(xi​))+Ω(ft​)+constant
  其中前t-1棵树的正则项均为常数，故归入constant留在式后暂且不作考虑。将(2)式中的损失函数做二阶泰勒展开，可得：
obj(t)≈∑in[l(yi,y^i(t−1))+gift(xi)+12hift(xi)2]+Ω(ft)+constant(2)obj^{(t)}\approx\sum_{i} ^ {n}[l(y_{i},\hat {y}_{i}^{(t-1)})+g_{i}f_{t}(x_{i})+\frac{1}{2}h_{i}f_{t}(x_{i})^{2}]+\Omega(f_{t})+constant(2)obj(t)≈i∑n​[l(yi​,y^​i(t−1)​)+gi​ft​(xi​)+21​hi​ft​(xi​)2]+Ω(ft​)+constant(2)
  为简化公式，式中gi=∂y^i(t−1)l(yi,y^i(t−1))g_{i}=\partial _{\hat y_{i}^{(t-1)}}l(y_{i},\hat y_{i}^{(t-1)})gi​=∂y^​i(t−1)​​l(yi​,y^​i(t−1)​),hi=∂y^i(t−1)2l(yi,y^i(t−1))h_{i}=\partial^{2} _{\hat y_{i}^{(t-1)}}l(y_{i},\hat y_{i}^{(t-1)})hi​=∂y^​i(t−1)​2​l(yi​,y^​i(t−1)​)，此时分析可知在式中∑inl(yi,y^i(t−1)+ft(xi))\sum_{i} ^ {n}l(y_{i},\hat {y}_{i}^{(t-1)}+f_{t}(x_{i}))∑in​l(yi​,y^​i(t−1)​+ft​(xi​))为常数，可并入constant项（大于0）。为了使目标函数obj得到最小值，必将把constant项除去，由（2）可得：
objmin(t)=∑in[l(yi,y^i(t−1))+gift(xi)+12hift(xi)2]+Ω(ft)obj^{(t)}_{min}=\sum_{i} ^ {n}[l(y_{i},\hat {y}_{i}^{(t-1)})+g_{i}f_{t}(x_{i})+\frac{1}{2}h_{i}f_{t}(x_{i})^{2}]+\Omega(f_{t})objmin(t)​=i∑n​[l(yi​,y^​i(t−1)​)+gi​ft​(xi​)+21​hi​ft​(xi​)2]+Ω(ft​)
  此时设q(x)为一映射，表示x所分到的叶子节点，而w为某叶子节点上的分数。则有wq(xi)=ft(xi)w_{q(x_i)}=f_{t}(x_{i})wq(xi​)​=ft​(xi​)表示x所在叶子节点的分数，则可以定义正则化项为：
Ω(ft)=γT+12λ∑j=1Twj2(3)\Omega(f_{t})=\gamma T+ \frac {1}{2}\lambda \sum _{j=1}^{T}w_{j}^{2}(3)Ω(ft​)=γT+21​λj=1∑T​wj2​(3)
将（3）代入上式，即可得：
objmin(t)=∑in[giwq(xi)+12hiwq(xi)2]+12λ∑j=1Twj2+γTobj^{(t)}_{min}=\sum_{i} ^ {n}[g_{i}w_{q(x_{i})}+\frac{1}{2}h_{i}w_{q(x_{i})}^{2}]+\frac {1}{2}\lambda \sum _{j=1}^{T}w_{j}^{2}+\gamma Tobjmin(t)​=i∑n​[gi​wq(xi​)​+21​hi​wq(xi​)2​]+21​λj=1∑T​wj2​+γT
  令Gj=∑i∈IjgiG_{j}=\sum _{i\in{I_{j}}}g_{i}Gj​=∑i∈Ij​​gi​,Hj=∑i∈IjhiH_{j}=\sum _{i\in{I_{j}}}h_{i}Hj​=∑i∈Ij​​hi​,其中IjI_{j}Ij​是分到 j 叶子节点的所有样本的集合。
  易得简化后的目标函数：
objmin(t)=∑j=1T[Giwj+12(Hi+λ)wj2]+γTobj^{(t)}_{min}=\sum_{j=1} ^ {T}[G_{i}w_{j}+\frac{1}{2}(H_{i}+\lambda)w_{j}^{2}]+\gamma Tobjmin(t)​=j=1∑T​[Gi​wj​+21​(Hi​+λ)wj2​]+γT
进一步求最值，即对wjw_{j}wj​求最值，易得wj∗=−GjHj+λw_{j}^{*}=-\frac{G_{j}}{H_{j}+\lambda}wj∗​=−Hj​+λGj​​，obj∗=−12∑j=1TGj2Hj+λ+γTobj^{*}=-\frac {1}{2}\sum_{j=1}^{T}\frac{G_{j}^{2}}{H_{j}+\lambda}+ \gamma Tobj∗=−21​∑j=1T​Hj​+λGj2​​+γT

推导完成！

====================================================================================

4.一些细节的解释

（1）为什么要泰勒展开啊？为啥还是二阶？

  让我们回顾一下前文留下的伏笔！记不记得有一个MSE函数？这时候可以代进obj(t)obj^{(t)}obj(t)试一试：
obj(t)=∑in[yi−(y^(t−1)+ft(xi))]2+Ω(ft)+constantobj^{(t)}=\sum_{i}^{n}[y_{i}-(\hat{y}^{(t-1)}+f_{t}(x_{i}))]^{2}+\Omega(f_{t})+constantobj(t)=i∑n​[yi​−(y^​(t−1)+ft​(xi​))]2+Ω(ft​)+constant
=∑in[(yi−y^(t−1))−ft(xi)]2+Ω(ft)+constant=\sum_{i}^{n}[(y_{i}-\hat{y}^{(t-1)})-f_{t}(x_{i})]^{2}+\Omega(f_{t})+constant=i∑n​[(yi​−y^​(t−1))−ft​(xi​)]2+Ω(ft​)+constant
=∑in[yi2−2yiy^(t−1)+(y^(t−1))2−2(yi−y^(t−1))ft(xi)+ft2(xi)]+Ω(ft)+constant=\sum_{i}^{n}[y_{i}^{2}-2y_{i}\hat{y}^{(t-1)}+(\hat{y}^{(t-1)})^{2}-2(y_{i}-\hat{y}^{(t-1)})f_{t}(x_{i})+f_{t}^{2}(x_{i})]+\Omega(f_{t})+constant=i∑n​[yi2​−2yi​y^​(t−1)+(y^​(t−1))2−2(yi​−y^​(t−1))ft​(xi​)+ft2​(xi​)]+Ω(ft​)+constant
将常数项合并：
=∑in[2(y^(t−1)−yi)ft(xi)+ft2(xi)]+Ω(ft)+constant=\sum_{i}^{n}[2(\hat{y}^{(t-1)}-y_{i})f_{t}(x_{i})+f_{t}^{2}(x_{i})]+\Omega(f_{t})+constant=i∑n​[2(y^​(t−1)−yi​)ft​(xi​)+ft2​(xi​)]+Ω(ft​)+constant
  上式中可以看出有非常整齐的一次项及二次项，这也是Xgboost对于传统boosting方式（如GBDT）的优点，考虑了二阶的损失函数，相当于在梯度下降的过程中不仅考虑了下降的速度，还考虑了下降的加速度，这就使得学习器在训练过程中的视野更远，大局观更强，相对不那么容易陷入到局部最优的陷阱当中。
  但这种干净形式不是所有损失函数都可以直接得出，比如伏笔中的logloss（由于计算太麻烦以至于笔者也不想算= =）。所以就需要创在一种通用的损失函数构造方法，使目标函数obj可以得到上述的一阶、二阶形式。
  这个时候泰勒展开就站了出来！
f(x+Δx)=f(x)+f(x)′Δx+12f(x)′′Δx+O(x3)f(x+\Delta x )=f(x)+{f(x)}' \Delta x+\frac{1}{2}{f(x)}''\Delta x+O(x^{3})f(x+Δx)=f(x)+f(x)′Δx+21​f(x)′′Δx+O(x3)
而原式中的l(yi,y^i(t−1)+ft(xi))l(y_{i},\hat {y}_{i}^{(t-1)}+f_{t}(x_{i}))l(yi​,y^​i(t−1)​+ft​(xi​))就等价于f(x+Δx)f(x+\Delta x )f(x+Δx)，舍去高阶小项代入即可得到推导中的（2）式。
obj(t)≈∑in[l(yi,y^i(t−1))+gift(xi)+12hift(xi)2]+Ω(ft)+constantobj^{(t)}\approx\sum_{i} ^ {n}[l(y_{i},\hat {y}_{i}^{(t-1)})+g_{i}f_{t}(x_{i})+\frac{1}{2}h_{i}f_{t}(x_{i})^{2}]+\Omega(f_{t})+constantobj(t)≈i∑n​[l(yi​,y^​i(t−1)​)+gi​ft​(xi​)+21​hi​ft​(xi​)2]+Ω(ft​)+constant

（2）Xgboost的并行体现在哪里？

  使用过Xgboost的同学肯定会记得，在训练XGB的时候电脑CPU占用率会高达百分之百！（不能一边GridSearch一边看直播了 😦 ）。这部分体现了Xgboost速度来源——多线程并行训练，那么它的并行具体体现在哪里呢？

1. 可以发现在计算一阶导和二阶导的时候，采用了并行处理。在GiG_{i}Gi​、HiH_{i}Hi​的计算中，样本与样本之间的运算是相互独立的，因而可以多线程地批量计算G、H的值，最终一同代入目标函数进行求解。（体现在同一棵树上的并行）

2. 在特征粒度上的也存在并行。决策树的学习最耗时的一个步骤就是对特征的值进行排序（因为要确定最佳分割点），xgboost在训练之前，预先对数据进行了排序，然后保存为block结构，后面的迭代中重复地使用这个结构，大大减小计算量（当然这个步骤也使得先期会消耗一段时间，在lightGBM中对这一点有所改善）。这个block结构也使得并行成为了可能，在进行节点的分裂时，需要计算每个特征的增益，最终选增益最大的那个特征去做分裂，那么各个特征的增益计算就可以开多线程进行。