原理进阶:XGBoost的求解流程与数学原理
学习XGBoost的数学原理需要大量梯度提升树(GBDT)相关知识,本课程将假设你已经非常熟悉梯度提升树的原理与特点。如果你不熟悉梯度提升树,强烈建议你回顾之前的课程。
1 XGBoost的基本数学流程
作为Boosting算法的里程碑,XGBoost以它高度复杂的数学原理闻名。无论你是阅读XGBoost原始论文,还是寻找经典书籍作为资料,都会发现数学占了很大的篇幅。然而事实上,随着Boosting算法的数学过程越来越复杂,实际建模过程和实现代码却是越来越简单的。数学并没有让算法变得更复杂,而是让算法变得更简单。
XGBoost最初是由陈天奇及华盛顿大学团队在2014年提出,在提出时就已经具有了今天我们要证明的全套复杂数学理论。但在搭建xgboost算法库的过程中,研发团队不断吸收当代各种经典算法的元素,形成更加复杂的算法系统:例如,DART树是吸收了深度学习中的Dropout技术,我们没有详细讲解的基于直方图的估计贪婪算法是借鉴LightGBM的技术,参考以下论文:
《XGBoost:一种可拓展的提升树系统》
Chen,T.Q.; Geustrin, C. (2014). “XGBoost: A Scalable Tree Boosting System”
《DART树:当Dropout遇见自适应回归树》
Rashmi, K.V.; Ran Gilad-Bachrach.(2015) “DART: Dropouts meet Multiple Additive Regression Trees”
《Lightgbm:一种极度高效的梯度提升树》
Ke, G. et al.(2017).“Lightgbm: A highly efficient gradient boosting decision tree”
你可以在本周上传的课程中找到以上三篇论文。需要注意的是,在本节中我们将会专注于最原始的XGBoost论文,而不会涉及到其他两篇文献的内容。其中,XGBoost基本数学流程是总结自GBDT的流程与XGBoost论文中涉及到的伪代码,而所有证明过程都来自于XGBoost原始论文。同时,为了与之前的课程保持一致,我们将使用与原论文略有不同的符号,但数学过程将是完全相同的。在阅读原论文时,记得辨析不同的数学符号。
不要忘记,对任何Boosting算法来说我们都有如下流程(对大部分算法来说是损失函数L,对XGBoost来说是目标函数O):
依据上一个弱评估器f(x)k−1的结果,计算目标函数O,
并使用O自适应地影响下一个弱评估器f(x)k的构建。集成模型输出的结果,受到整体所有弱评估器f(x)0 ~ f(x)K的影响。
假设现有数据集N,含有形如(xi,yi)的样本M个,i为任意样本的编号,单一样本的损失函数为l(yi,H(xi)),其中H(xi)是i号样本在集成算法上的预测结果,整个算法的损失函数为L(y,H(x)),且总损失等于全部样本的损失之和:L(y,H(x))=∑il(yi,H(xi))。目标函数中使用L2正则化(λ为0,α为0),并且γ不为0。
同时,弱评估器为回归树f,总共学习K轮(注意在GBDT当中我们使用的是大写字母T来表示迭代次数,由于在XGBoost当中字母T被用于表示目标函数中的叶子总量,因此我们在这里使用字母K表示迭代次数)。则XGBoost回归的基本流程如下所示:
初始化数据迭代的起点H0(x)。在应用xgboost时,我们可以指定任意数字来作为H0(x),但在xgboost原始论文当中,并未详细讨论如何计算迭代的初始值。考虑到XGBoost在许多方面继承了梯度提升树GBDT的思想,我们可以使用公式来计算XGBoost的H0:
H0(x)=argminCi=1∑Ml(yi,C)=argminCL(y,C)
其中yi为真实标签,C为任意常数。以上式子表示,找出令∑i=1Ml(yi,C)最小的常数C值,并输出最小的∑i=1Ml(yi,C)作为H0(x)的值。需要注意的是,由于H0(x)是由全部样本的l计算出来的,因此所有样本的初始值都是H0(x),不存在针对某一样本的单一初始值。
由于在初始的时候没有树结构,因此没有复杂度等信息,因此没有使用目标函数求初始值,而是使用了损失函数。在GBDT的数学过程当中,我们详细展示过如何求解令初始损失最小的C(对损失求一阶导数并让一阶导数为0),并且我们详细证明过,当损失函数为MSE时,令整体初始损失最小的C值就是y的均值。对XGBoost来说这一切都成立,只不过在xgboost库中我们默认的初始值为0.5。
开始循环,for k in 1,2,3…K:
- 2) 抽样
在现有数据集N中,抽样M * subsample
个样本,构成训练集Nk
- 3) 求拟合项
对任意一个样本i,计算一阶导数gik,二阶导数hik,以及伪残差(pseudo-residuals)rik,具体公式为:gik=∂Hk−1(xi)∂l(yi,Hk−1(xi))
hik=∂Hk−12(xi)∂2l(yi,Hk−1(xi))
rik=−hikgik
不难发现,伪残差是一个样本的一阶导数除以二阶导数并取负的结果,并且在进行第k次迭代、计算第k个导数时,我们使用的是前k-1次迭代后输出的集成算法结果。同时,我们是先对目标函数l中的自变量H(x)求导,再令求导后的结果等于Ht−1(xi)的值,并不是直接对Ht−1(xi)这一常数求导。
对常数求导,以及对变量求导是两个概念,举例说明:
l=x2+x
对常数0求导:
∂0∂l=0∂(x2+x)=0
对变量x求导并让x=0:
∂x∂l=∂x∂(x2+x)=2x+1=2∗0+1=1
因此,gik标准的写法应该是:
gik=[∂H(xi)∂l(yi,H(xi))]H(xi)=Hk−1 (xi)
在实际推导过程中,为公式简洁,简写为上述流程中的写法。
在k=1时,所有求导计算过程中的Hk−1(xi)都等于初始H0(x),在k>1时,每个样本上的Hk−1(xi)都是不同的取值。
-
4) 建树
求解出伪残差后,在数据集(xi,rik)上按colsample_by*
规则进行抽样,再按照结构分数增益规则建立一棵回归树fk。注意在这个过程中,训练时拟合的标签为样本的伪残差rik,并且叶子节点j的结构分数和任意分枝时的结构分数增益的公式为:
Scorej=∑i∈jhi+λ(∑i∈jgi)2
Gain=21(∑i∈Lhi+λ(∑i∈Lgi)2+∑i∈Rhi+λ(∑i∈Rgi)2−∑i∈Phi+λ(∑i∈Pgi)2)−γ
建树过程不影响任何gik与hik的值。
-
5) 输出树上的结果
建树之后,依据回归树fk的结构输出叶子节点上的输出值(预测值)。对任意叶子节点j来说,输出值为:
假设样本i被分割到叶子j上,则有:
fk(xi)=wj
使用字母w表示叶子节点上的输出值是XGBoost论文所规定的,我们曾经见过一次w,你还记得在哪里吗?在我们介绍XGBoost的目标函数时,L2正则项的表达式为21λ∑j=1Twj2。我们曾说过w代表XGBoost中的叶子权重,实际上叶子权重就是叶子上的输出值。为不和其他权重混淆,之后我们统一称呼w为输出值。
不难发现,叶子节点上的输出值与结构分数很相似,只不过结构分数的分子上是平方,而输出值的分子上没有平方。在数学上我们可以证明,该输出值能让目标函数最快减小。
你可能注意到了,在迭代刚开始时我们已经知道了输出值式子中所需的所有g和h。为什么还要建树呢?只有当我们建立了决策树,我们才能够知道具体哪些样本i在叶子节点j上。因此树fk提供的是结构信息。
由于任意样本必然被分到任意叶子上,因此对整棵树fk来说,任意fk(xi)一定有对应的w。
- 6) 迭代
根据预测结果fk(xi)迭代模型,具体来说:
Hk(xi)=Hk−1(xi)+fk(xi)
假设输入的步长为η,则Hk(x)应该为:
Hk(xi)=Hk−1(xi)+ηfk(xi)
对整个算法则有:
Hk(x)=Hk−1(x)+ηfk(x)
- 7) 循环结束
输出HK(x)的值作为集成模型的输出值。
以上就是XGBoost的完整数学流程。不难发现,作为从GBDT改进而来的算法,XGBoost在基础数学流程上基本继承了GBDT的流程(7步走的流程与GBDT一模一样,同时也有继承伪残差等细节),但又在具体每个流程中都做出了改进,进一步简化了Boosting算法的运算流程——比如说,虽然整个算法持续再向降低目标函数的方向运行,但整个过程中不存在任何的求最优解的数学计算。除了建树流程以外,其他流程都是非常简单的按公式计算而已。
由于XGBoost原始论文中并不存在上述流程的完整说明,因此要梳理出该流程并不容易,但如果我们对GBDT的流程相对熟悉,那XGBoost的流程也并不难。对XGBoost来说,真正难度较大的部分并不是梳理以上算法流程,而是证明这一流程可以让模型向着目标函数最小化的方向运行。在这个流程中包括如下很明显的问题:
-
- 建树时拟合的rik=−hikgik究竟是什么?拟合它有什么意义?
-
- 结构分数和结构分数增益的公式是如何推导出来的?为什么这样建树可以提升模型的效果?
-
- 为什么叶子节点的输出值wj是−∑i∈jhik+λ(∑i∈jgik)?这样输出有什么意义?
-
- 课程的第一部分说XGBoost拟合的也是残差,残差在哪里?
接下来我们就将展示完整的证明来回答这些问题。
2 化简XGBoost的目标函数
首先,根据之前对目标函数的定义,XGBoost中目标函数是针对一棵树的目标函数,而不是针对一个样本或一整个算法的目标函数。并且,任意树的目标函数都包括三大部分:损失函数l、叶子数量T以及正则项。具体地来说:
假设单一树fk的目标函数为Ok,总共有T片叶子,该树上任意样本i的损失函数为l((yi,H(xi)),其中H(xi)是i号样本在集成算法上的预测结果。树上总共有M个样本,目标函数中使用L2正则化(λ不为0,α为0),并且γ不为0,则该树的目标函数为:
Ok=i=1∑Ml(yi,Hk(xi))+γT+21λj=1∑Twj2
我们的目标是令目标函数最小,并找出令目标函数最小的某个自变量。对使用普通损失函数的Boosting算法来说,算法的输出值H(x)是在迭代过程中不断变化的,损失函数l(y,H(x))也是在迭代中不断变小的:
Hk(xi)=Hk−1(xi)+fk(xi)
lk=l(yi,Hk−1(xi)+fk(xi))
当迭代到第k次时,损失函数中的yi与Hk−1(xi)都是常数,只有fk(xi)是变量,因此我们只需要在损失函数上对fk(xi)求导,并找到令整体损失函数最小的预测值fk(xi)即可。在GBDT当中,我们提到过,无论弱评估器fk是什么结构、什么规则、如何建立、如何拟合,只要其最终的输出值fk(xi)是令整体损失函数L最小化的fk(xi),那随着算法逐步迭代,损失函数必然会越来越小。因此,一个适合的fk(xi)不仅能保证损失持续减小,还可以指导单个评估器的建立。
在GBDT当中,我们证明了令GBDT整体损失函数最小化的fk(xi)就是损失函数的负梯度−gi(详见GBDT课程《四 4:证明:拟合伪残差可以令损失函数最快地减小》),也因此GBDT在建树时拟合负梯度。当损失函数为21MSE时,GBDT中的负梯度在数值上就等于残差,因此GBDT是拟合残差的算法。
在XGBoost当中,我们也可以对目标函数求导、并找出令目标函数最小的某个自变量,但问题在于,XGBoost的目标函数中存在多个自变量:
Ok=i=1∑Ml(yi,Hk(xi))+γT+21λj=1∑Twj2=i=1∑Ml(yi,Hk−1(xi)+fk(xi))+γT+21λj=1∑Twj2
其中,T是第k棵树上的叶子总量,fk(xi)与wj都是模型输出的预测值(叶子上的输出值),不过表现形式不同,对任意位于叶子j上的样本i而言,数值上fk(xi)=wj。对XGBoost来说,只能选择一个变量作为自变量,考虑到fk(xi)只与单个样本的精确程度有关,而T只与树结构有关,XGBoost论文最终选择了即与精确度有关、又与树结构有关的变量wj。同时,如果知道叶子的最佳输出值wj就可以引导树成长为合理的结构,但只知道叶子的总量T是无法指导建树的。
因此,求解XGBoost目标函数的第一步,就是将目标函数尽量整理成以wj表示的形式。
在GBDT当中,我们借助一阶泰勒展开化简了损失函数,借鉴于GBDT的思路,XGBoost也使用了泰勒展开。具体来看:
在数学中,泰勒级数使用无限个项的连加式来表示一个函数。实际应用当中,我们一般取有限项的连加式来逼近一个函数。当总共有N项时,连加式被叫做N阶泰勒展开(Nth-order Taylor approximation)。假设现在存在函数f(x),则有:
其中(x-a)是非常小的任意实数/复数,n!是n的阶乘,f(n)(a)是函数f(x)的n阶导数在a点的取值。当a为0时,泰勒级数也被叫做麦克劳思级数。
f(x)≈n=0∑2n!f(n)(a)(x−a)n≈f(a)+1!f′(a)(x−a)+2!f′′(a)(x−a)2
阶数越大,泰勒展开的值越接近$f(x)$的真实值。
在我们的目标函数$O_k$中,**可以被泰勒展开的是第一部分损失函数**$l$:
$$O_k = \sum_{i=1}^Ml \left( y_i,H_{k-1}(x_i) + f_k(x_i) \right) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^T w_j^2
由于损失函数l中只有唯一变量Hk−1(xi)+fk(xi),因此可以将函数简写为l(Hk−1(xi)+fk(xi))。
根据二阶泰勒展开,已知:
f(x)≈n=0∑2n!f(n)(a)(x−a)n≈f(a)+1!f′(a)(x−a)+2!f′′(a)(x−a)2
令泰勒展开中的x=Hk−1(xi)+fk(xi),令泰勒展开中的a=Hk−1(xi),则(x−a)=fk(xi)。据此,损失函数l(Hk−1(xi)+fk(xi))可以被表示为:
l(Hk−1(xi)+fk(xi))≈l(Hk−1(xi))+∂Hk−1(xi)∂l(Hk−1(xi))∗fk(xi)+2∂Hk−12(xi)∂2l(Hk−1(xi))∗fk2(xi)
在XGBoost中我们定义过损失函数的一阶导数与二阶导数:
gik=∂Ht−1(xi)∂l(yi,Hk−1(xi))
hik=∂Ht−12(xi)∂2l(yi,Hk−1(xi))
在XGBoost原论文中,为了公式简洁,gi和hi并没有呈现下标k,但我们已经很清楚:g与h是在每一轮迭代时需要被重新计算的。在这里我们也参照原论文中的做法去掉下标k。因此,经过泰勒展开后的式子可以化简为:
l(Hk−1(xi)+fk(xi))≈l(Hk−1(xi))+gifk(xi)+21hifk2(xi)≈常数+gifk(xi)+21hifk2(xi)
不难发现,该式子中Hk−1(xi)是常数,因此第一部分l(Ht−1(xi))也是一个常数,常数无法被最小化,因此我们可以将常数从该目标函数中剔除。经过泰勒展开,目标函数变为:
O~k=i=1∑M(gifk(xi)+21hifk2(xi))+γT+21λj=1∑Twj2=i=1∑Mgifk(xi)+21i=1∑Mhifk2(xi)+γT+21λj=1∑Twj2
现在目标函数的前两项分别代表所有样本的gifk(xi)之和,以及所有样本的hifk2(xi)之和乘1/2。别忘记,我们选择的唯一的自变量是wj,因此我们希望能够将fk以某种方式转化为wj。之前已经提到过多次,对任意位于叶子j上的样本i而言,数值上fk(xi)=wj,我们可以尝试着从一个样本开始进行转化:
对于单一样本i,假设这个样本位于叶子j上,应该有:
gifk(xi)=giwj
对于一片叶子j,我们可以计算这片叶子上所有样本的giwj之和:
i∈j∑giwj
而一片叶子上所有样本的wj都是一致的,因此一片叶子上的giwj之和可以转变为:
i∈j∑giwj=g1wj + g2wj + ... + gnwj,其中1,2...n是叶子j上的样本=wji∈j∑gi
假设现在一共有T片叶子,则整棵树上所有样本的giwj之和为:
j=1∑T(wji∈j∑gi)
所以:
i=1∑Mgifk(xi)=j=1∑T(wji∈j∑gi)
同理,单一样本i的hifk2(xi)也可以以相同方式转化。对单一样本:
hifk2(xi)=hiwj2
对一片叶子:
i∈j∑hiwj2=h1wj2 + h2wj2 + ... + hnwj2,其中1,2...n是叶子j上的样本=wj2i∈j∑hi
对整棵树:
i=1∑Mhifk2(xi)=j=1∑T(wj2i∈j∑hi)
因此对整个目标函数有:
O~k=i=1∑Mgifk(xi)+21i=1∑Mhifk2(xi)+γT+21λj=1∑Twj2=j=1∑T(wji∈j∑gi+21wj2i∈j∑hi)+γT+21λj=1∑Twj2
不难发现,现在正则项可以与原来损失函数的部分合并了:
=j=1∑T(wji∈j∑gi+21wj2i∈j∑hi+21λwj2)+γT=j=1∑T(wji∈j∑gi+21wj2(i∈j∑hi+λ))+γT
合并之后,整个目标函数变为两项,一项是所有叶子上的(损失+正则)之和,另一项是叶子总量。现在,我们可以开始求解最小目标函数以及对应的最优自变量wj了。
3 求解XGBoost的目标函数
首先,令目标函数中的叶子总量最小是不可能的,过度降低叶子总量会大幅度伤害模型的学习能力,因此我们只能考虑令所有叶子上的(损失+正则)之和最小。
其次,当树建好之后,叶子与叶子之间是相互独立的,因此每片叶子上的(损失+正则)也是相互独立的。我们只要令每片叶子的(损失+正则)都最小,就可以保证全部叶子的(损失+正则)之和最小。故而,我们要令式子中标注为红色的部分最小:
O~k=j=1∑Twji∈j∑gi+21wj2(i∈j∑hi+λ)+γT
将标注为红色的部分命名为μj,表示叶子j上的损失+正则。则有:
μj=wji∈j∑gi+21wj2(i∈j∑hi+λ)
现在,对叶子j而言,在μj上对唯一自变量wj求导,则有:
∂wj∂μj=∂wj∂wj∑i∈jgi+21wj2(∑i∈jhi+λ)=i∈j∑gi+wj(i∈j∑hi+λ)
令一阶导数为0,则有:
i∈j∑gi+wj(i∈j∑hi+λ)wj(i∈j∑hi+λ)wj=0=−i∈j∑gi=−∑i∈jhi+λ∑i∈jgi
你应该发现了,对一片叶子来说,令目标函数最小的wj就是我们之前提过的叶子权重,也就是XGBoost数学流程当中叶子上的输出值。如果要令叶子的输出非常接近叶子权重公式,那应该如何拟合每个样本呢?
对任意位于叶子j上的样本i来说:
μi=wjgi+21wj2hi
将一片叶子上的μj转变成μi时,原则上需要将μj中的每一项都转换为单个样本所对应的项,然而在转换正则项时则存在问题:与∑i∈jgi这样可以直接指向单个样本的项不同,λ是针对与一片叶子设置的值,如果要将λ转变为针对单一样本的正则项,则需要知道当前叶子上一共有多少样本。然而,拟合发生在建树之前,因此在这一时间点不可能知道一片叶子上的样本总量,因此在xgboost的实际实现过程当中,拟合每一片叶子时不涉及正则项,只有在计算结构分数与叶子输出值时才使用正则项。
对μi上唯一的自变量wj求导,则有:
∂wj∂μi=∂wj∂(wjgi+21wj2hi)=gi+wjhi
令一阶导数为0,则有:
gi+wjhiwjhiwj=0=−gi=−higi
对任意样本i而言,令目标函数最小的最优wj就是我们的伪残差ri,也就是XGBoost数学流程当中用于进行拟合的拟合值。
现在,我们把令目标函数最小的最优wj带回到μj中,查看μj如何变化:
μj=wji∈j∑gi+21wj2(i∈j∑hi+λ)=−∑i∈jhi+λ∑i∈jgi∗i∈j∑gi+21(−∑i∈jhi+λ∑i∈jgi)2∗i∈j∑hi+λ=−∑i∈jhi+λ(∑i∈jgi)2+21∑i∈jhi+λ(∑i∈jgi)2=−21∑i∈jhi+λ(∑i∈jgi)2
因此,目标函数(所有叶子上的损失)就可以变为:
O~k=j=1∑Twji∈j∑gi+21wj2(i∈j∑hi+λ)+γT=j=1∑T(−21∑i∈jhi+λ(∑i∈jgi)2)+γT
因此,一片叶子上的目标函数就是:
Oj=−21∑i∈jhi+λ(∑i∈jgi)2+γ
对任意一片叶子来说,目标函数可以衡量叶子的质量,其中γ是可以设定的超参数,21为常数,因此对任意叶子,我们希望标注为红色的部分越小越好:
Oj=21(−∑i∈jhi+λ(∑i∈jgi)2)+γ
故而,我们希望以下式子越大越好:
∑i∈jhi+λ(∑i∈jgi)2
这个式子,正是XGBoost用于分枝时的指标“结构分数”(Structure Score)。
当分枝的时候,我们希望目标函数越小越好,因此在分枝过程中,父节点的目标函数是大于子节点的目标函数的,因此我们可以使用(父节点目标函数 - 子节点目标函数之和)来衡量分枝的质量,则有:
Gain=Op−(Ol+Or)=−21∑i∈Phi+λ(∑i∈Pgi)2+γ−(−21∑i∈Lhi+λ(∑i∈Lgi)2+γ−21∑i∈Rhi+λ(∑i∈Rgi)2+γ)=−21∑i∈Phi+λ(∑i∈Pgi)2+γ+21∑i∈Lhi+λ(∑i∈Lgi)2−γ+21∑i∈Rhi+λ(∑i∈Rgi)2−γ=21(∑i∈Lhi+λ(∑i∈Lgi)2+∑i∈Rhi+λ(∑i∈Rgi)2−∑i∈Phi+λ(∑i∈Pgi)2)−γ=21(ScoreL+ScoreR−ScoreP)−γ
其中,γ是可以设定的超参数,21为常数,因此:
Gain=ScoreL+ScoreR−ScoreP
这就是我们在分枝时所使用的结构分数增益了。
现在你发现了,XGBoost流程中所使用的全部新公式(包括独特的拟合值、独特的分枝指标、独特的输出值)都是通过令目标函数最小而求解出来的。因此,XGBoost整个流程就保证了目标函数一定是向着最小化方向进行迭代的,新生成的每片叶子上的输出值wj都是会令目标函数最小化的输出值。现在,你可以回答以下的问题了:
-
- 建树时拟合的rik=−hikgik究竟是什么?拟合它有什么意义?
-
- 结构分数和结构分数增益的公式是如何推导出来的?为什么这样建树可以提升模型的效果?
-
- 为什么叶子节点的输出值wj 是−∑i∈jhik+λ(∑i∈jgik)?这样输出有什么意义?
-
- 课程的第一部分说XGBoost拟合的也是残差,残差在哪里?
唯一余留问题4其实我们早已解答过,回去查看第二部分中的《2.2 弱评估器的分枝》,你会发现其实我们早已经求出残差了。当目标函数为21MSE,负梯度−gi就等于残差,而hi=1,因此拟合项−higi自然也是残差本身了。因此,XGBoost也是拟合负梯度的算法,并且在特定损失函数下,XGBoost也拟合残差。在了解这个推导过程之后,再返回复习XGBoost的整体数学流程,你会发现数学真的让算法变得更简单,而不是更复杂。如果你已经学到了这里,你已经跨过了机器学习原理中几乎最难的一部分。在接下来的课程中,所有数学推导对你来说都不再是阻碍了。