原理进阶: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算法来说我们都有如下流程(对大部分算法来说是损失函数LL,对XGBoost来说是目标函数OO):


依据上一个弱评估器f(x)k1f(x)_{k-1}的结果,计算目标函数OO

并使用OO自适应地影响下一个弱评估器f(x)kf(x)_k的构建。集成模型输出的结果,受到整体所有弱评估器f(x)0f(x)_0 ~ f(x)Kf(x)_K的影响。

假设现有数据集NN,含有形如(xi,yi)(x_i,y_i)的样本MM个,ii为任意样本的编号,单一样本的损失函数为l(yi,H(xi))l(y_i,H(x_i)),其中H(xi)H(x_i)ii号样本在集成算法上的预测结果,整个算法的损失函数为L(y,H(x))L(y,H(x)),且总损失等于全部样本的损失之和:L(y,H(x))=il(yi,H(xi))L(y,H(x)) = \sum_i l(y_i,H(x_i))。目标函数中使用L2正则化(λ\lambda为0,α\alpha为0),并且γ\gamma不为0。

同时,弱评估器为回归树ff,总共学习KK轮(注意在GBDT当中我们使用的是大写字母T来表示迭代次数,由于在XGBoost当中字母T被用于表示目标函数中的叶子总量,因此我们在这里使用字母K表示迭代次数)。则XGBoost回归的基本流程如下所示:

  • 1) 初始化

初始化数据迭代的起点H0(x)H_0(x)。在应用xgboost时,我们可以指定任意数字来作为H0(x)H_0(x),但在xgboost原始论文当中,并未详细讨论如何计算迭代的初始值。考虑到XGBoost在许多方面继承了梯度提升树GBDT的思想,我们可以使用公式来计算XGBoost的H0H_0

H0(x)=argminCi=1Ml(yi,C)=argminCL(y,C)\begin{aligned} H_0(x) &= \mathop{argmin}_{C} \sum_{i=1}^M l(y_i,C)\\ \\ &= \mathop{argmin}_{C} L(y,C) \end{aligned}

其中yiy_i为真实标签,CC为任意常数。以上式子表示,找出令i=1Ml(yi,C)\sum_{i=1}^Ml(y_i,C)最小的常数CC值,并输出最小的i=1Ml(yi,C)\sum_{i=1}^Ml(y_i,C)作为H0(x)H_0(x)的值。需要注意的是,由于H0(x)H_0(x)是由全部样本的ll计算出来的,因此所有样本的初始值都是H0(x)H_0(x),不存在针对某一样本的单一初始值。

由于在初始的时候没有树结构,因此没有复杂度等信息,因此没有使用目标函数求初始值,而是使用了损失函数。在GBDT的数学过程当中,我们详细展示过如何求解令初始损失最小的CC(对损失求一阶导数并让一阶导数为0),并且我们详细证明过,当损失函数为MSE时,令整体初始损失最小的CC值就是yy的均值。对XGBoost来说这一切都成立,只不过在xgboost库中我们默认的初始值为0.5。

开始循环,for k in 1,2,3…K:

  • 2) 抽样


    在现有数据集NN中,抽样MM * subsample个样本,构成训练集NkN^k

  • 3) 求拟合项


    对任意一个样本ii,计算一阶导数gikg_{ik},二阶导数hikh_{ik},以及伪残差(pseudo-residuals)rikr_{ik},具体公式为:

    gik=l(yi,Hk1(xi))Hk1(xi)g_{ik} = \frac{\partial{l(y_i,H_{k-1}(x_i))}}{\partial{H_{k-1}(x_i)}}

    hik=2l(yi,Hk1(xi))Hk12(xi)h_{ik} = \frac{\partial^2{l(y_i,H_{k-1}(x_i))}}{\partial{H^2_{k-1}(x_i)}}

    rik=gikhikr_{ik} = -\frac{g_{ik}}{h_{ik}}

不难发现,伪残差是一个样本的一阶导数除以二阶导数并取负的结果,并且在进行第k次迭代、计算第k个导数时,我们使用的是前k-1次迭代后输出的集成算法结果。同时,我们是先对目标函数ll中的自变量H(x)H(x)求导,再令求导后的结果等于Ht1(xi)H_{t-1}(x_i)的值,并不是直接对Ht1(xi)H_{t-1}(x_i)这一常数求导。


对常数求导,以及对变量求导是两个概念,举例说明:

l=x2+xl = x^2+x

对常数0求导:

l0=(x2+x)0=0\frac{\partial{l}}{\partial{0}} = \frac{\partial{(x^2 + x)}}{0} = 0

对变量x求导并让x=0:

lx=(x2+x)x=2x+1=20+1=1\frac{\partial{l}}{\partial{x}} = \frac{\partial{(x^2 + x)}}{\partial x} = 2x + 1 = 2*0 + 1 = 1

因此,gikg_{ik}标准的写法应该是:

gik=[l(yi,H(xi))H(xi)]H(xi)=Hk1  (xi)g_{ik} = \big[\frac{\partial{l(y_i,H(x_i))}}{\partial{H(x_i)}}\big]_{H(x_i) = H_{k-1}\ \ (x_i)}

在实际推导过程中,为公式简洁,简写为上述流程中的写法。


在k=1时,所有求导计算过程中的Hk1(xi)H_{k-1}(x_i)都等于初始H0(x)H_0(x),在k>1时,每个样本上的Hk1(xi)H_{k-1}(x_i)都是不同的取值。

  • 4) 建树


    求解出伪残差后,在数据集(xi,rik)(x_i, r_{ik})上按colsample_by*规则进行抽样,再按照结构分数增益规则建立一棵回归树fkf_k。注意在这个过程中,训练时拟合的标签为样本的伪残差rikr_{ik},并且叶子节点jj的结构分数和任意分枝时的结构分数增益的公式为:

    Scorej=(ijgi)2ijhi+λScore_j = \frac{(\sum_{i \in j}g_i)^2}{\sum_{i \in j}h_i + \lambda}

    Gain=12((iLgi)2iLhi+λ+(iRgi)2iRhi+λ(iPgi)2iPhi+λ)γGain = \frac{1}{2} \left( \frac{(\sum_{i \in L}g_i)^2}{\sum_{i \in L}h_i + \lambda} + \frac{(\sum_{i \in R}g_i)^2}{\sum_{i \in R}h_i + \lambda} - \frac{(\sum_{i \in P}g_i)^2}{\sum_{i \in P}h_i + \lambda} \right) - \gamma

    建树过程不影响任何gikg_{ik}hikh_{ik}的值。

  • 5) 输出树上的结果


    建树之后,依据回归树fkf_k的结构输出叶子节点上的输出值(预测值)。对任意叶子节点jj来说,输出值为:

    假设样本ii被分割到叶子jj上,则有:

    fk(xi)=wjf_k(x_i) = w_j

使用字母ww表示叶子节点上的输出值是XGBoost论文所规定的,我们曾经见过一次ww,你还记得在哪里吗?在我们介绍XGBoost的目标函数时,L2正则项的表达式为12λj=1Twj2\frac{1}{2} \lambda \sum_{j=1}^T w_j^2。我们曾说过ww代表XGBoost中的叶子权重,实际上叶子权重就是叶子上的输出值。为不和其他权重混淆,之后我们统一称呼ww为输出值。


不难发现,叶子节点上的输出值与结构分数很相似,只不过结构分数的分子上是平方,而输出值的分子上没有平方。在数学上我们可以证明,该输出值能让目标函数最快减小


你可能注意到了,在迭代刚开始时我们已经知道了输出值式子中所需的所有gghh。为什么还要建树呢?只有当我们建立了决策树,我们才能够知道具体哪些样本ii在叶子节点jj上。因此树fkf_k提供的是结构信息。


由于任意样本必然被分到任意叶子上,因此对整棵树fkf_k来说,任意fk(xi)f_k(x_i)一定有对应的ww

  • 6) 迭代


    根据预测结果fk(xi)f_k(x_i)迭代模型,具体来说:

    Hk(xi)=Hk1(xi)+fk(xi)H_k(x_i) = H_{k-1}(x_i) + f_k(x_i)

    假设输入的步长为η\eta,则Hk(x)H_k(x)应该为:

    Hk(xi)=Hk1(xi)+ηfk(xi)H_k(x_i) = H_{k-1}(x_i) + \eta f_k(x_i)

    对整个算法则有:

    Hk(x)=Hk1(x)+ηfk(x)H_k(x) = H_{k-1}(x) + \eta f_k(x)

  • 7) 循环结束

    输出HK(x)H_K(x)的值作为集成模型的输出值。

以上就是XGBoost的完整数学流程。不难发现,作为从GBDT改进而来的算法,XGBoost在基础数学流程上基本继承了GBDT的流程(7步走的流程与GBDT一模一样,同时也有继承伪残差等细节),但又在具体每个流程中都做出了改进,进一步简化了Boosting算法的运算流程——比如说,虽然整个算法持续再向降低目标函数的方向运行,但整个过程中不存在任何的求最优解的数学计算。除了建树流程以外,其他流程都是非常简单的按公式计算而已。

由于XGBoost原始论文中并不存在上述流程的完整说明,因此要梳理出该流程并不容易,但如果我们对GBDT的流程相对熟悉,那XGBoost的流程也并不难。对XGBoost来说,真正难度较大的部分并不是梳理以上算法流程,而是证明这一流程可以让模型向着目标函数最小化的方向运行。在这个流程中包括如下很明显的问题:

    1. 建树时拟合的rik=gikhikr_{ik} = -\frac{g_{ik}}{h_{ik}}究竟是什么?拟合它有什么意义?
    1. 结构分数和结构分数增益的公式是如何推导出来的?为什么这样建树可以提升模型的效果?
    1. 为什么叶子节点的输出值wjw_j(ijgik)ijhik+λ-\frac{(\sum_{i \in j} g_{ik})}{\sum_{i \in j} h_{ik} + \lambda}?这样输出有什么意义?
    1. 课程的第一部分说XGBoost拟合的也是残差,残差在哪里?

接下来我们就将展示完整的证明来回答这些问题。

2 化简XGBoost的目标函数

  • 定义目标函数与目标函数的自变量

首先,根据之前对目标函数的定义,XGBoost中目标函数是针对一棵树的目标函数,而不是针对一个样本或一整个算法的目标函数。并且,任意树的目标函数都包括三大部分:损失函数ll、叶子数量TT以及正则项。具体地来说:

假设单一树fkf_k的目标函数为OkO_k,总共有TT片叶子,该树上任意样本ii的损失函数为l((yi,H(xi))l((y_i,H(x_i)),其中H(xi)H(x_i)ii号样本在集成算法上的预测结果。树上总共有M个样本,目标函数中使用L2正则化(λ\lambda不为0,α\alpha为0),并且γ\gamma不为0,则该树的目标函数为:

Ok=i=1Ml(yi,Hk(xi))+γT+12λj=1Twj2O_k = \sum_{i=1}^Ml(y_i,H_k(x_i)) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^Tw_j^2

我们的目标是令目标函数最小,并找出令目标函数最小的某个自变量。对使用普通损失函数的Boosting算法来说,算法的输出值H(x)H(x)是在迭代过程中不断变化的,损失函数l(y,H(x))l(y,H(x))也是在迭代中不断变小的:

Hk(xi)=Hk1(xi)+fk(xi)H_k(x_i) = H_{k-1}(x_i) + f_k(x_i)

lk=l(yi,Hk1(xi)+fk(xi))l_k = l(y_i,H_{k-1}(x_i) + f_k(x_i))

当迭代到第kk次时,损失函数中的yiy_iHk1(xi)H_{k-1}(x_i)都是常数,只有fk(xi)f_k(x_i)是变量,因此我们只需要在损失函数上对fk(xi)f_k(x_i)求导,并找到令整体损失函数最小的预测值fk(xi)f_k(x_i)即可。在GBDT当中,我们提到过,无论弱评估器fkf_k是什么结构、什么规则、如何建立、如何拟合,只要其最终的输出值fk(xi)f_k(x_i)是令整体损失函数LL最小化的fk(xi)f_k(x_i),那随着算法逐步迭代,损失函数必然会越来越小。因此,一个适合的fk(xi)f_k(x_i)不仅能保证损失持续减小,还可以指导单个评估器的建立。

在GBDT当中,我们证明了令GBDT整体损失函数最小化的fk(xi)f_k(x_i)就是损失函数的负梯度gi-g_i(详见GBDT课程《四 4:证明:拟合伪残差可以令损失函数最快地减小》),也因此GBDT在建树时拟合负梯度。当损失函数为12MSE\frac{1}{2}MSE时,GBDT中的负梯度在数值上就等于残差,因此GBDT是拟合残差的算法。

在XGBoost当中,我们也可以对目标函数求导、并找出令目标函数最小的某个自变量,但问题在于,XGBoost的目标函数中存在多个自变量:

Ok=i=1Ml(yi,Hk(xi))+γT+12λj=1Twj2=i=1Ml(yi,Hk1(xi)+fk(xi))+γT+12λj=1Twj2\begin{aligned} O_k &= \sum_{i=1}^Ml(y_i,H_k(x_i)) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^Tw_j^2 \\ &= \sum_{i=1}^M l \left( y_i,H_{k-1}(x_i) + \boldsymbol{\color{red}{f_k(x_i)}} \right) + \gamma \boldsymbol{\color{red}T} + \frac{1}{2}\lambda\sum_{j=1}^T\boldsymbol{\color{red}{w_j}}^2 \end{aligned}

其中,TT是第kk棵树上的叶子总量,fk(xi)f_k(x_i)wjw_j都是模型输出的预测值(叶子上的输出值),不过表现形式不同,对任意位于叶子jj上的样本ii而言,数值上fk(xi)=wjf_k(x_i) = w_j。对XGBoost来说,只能选择一个变量作为自变量,考虑到fk(xi)f_k(x_i)只与单个样本的精确程度有关,而TT只与树结构有关,XGBoost论文最终选择了即与精确度有关、又与树结构有关的变量wjw_j。同时,如果知道叶子的最佳输出值wjw_j就可以引导树成长为合理的结构,但只知道叶子的总量TT是无法指导建树的。

因此,求解XGBoost目标函数的第一步,就是将目标函数尽量整理成以wjw_j表示的形式。

  • 对目标函数进行泰勒展开

在GBDT当中,我们借助一阶泰勒展开化简了损失函数,借鉴于GBDT的思路,XGBoost也使用了泰勒展开。具体来看:

在数学中,泰勒级数使用无限个项的连加式来表示一个函数。实际应用当中,我们一般取有限项的连加式来逼近一个函数。当总共有N项时,连加式被叫做N阶泰勒展开(Nth-order Taylor approximation)。假设现在存在函数f(x)f(x),则有:

  • 泰勒级数(无限项)

    f(x)=n=0f(n)(a)n!(xa)nf(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n

其中(x-a)是非常小的任意实数/复数,n!n!是n的阶乘,f(n)(a)f^{(n)}(a)是函数f(x)f(x)的n阶导数在a点的取值。当a为0时,泰勒级数也被叫做麦克劳思级数。

  • 一阶泰勒展开

    f(x)n=01f(n)(a)n!(xa)nf(a)+f(a)1!(xa)\begin{aligned} f(x) &\approx \sum_{n=0}^{1}\frac{f^{(n)}(a)}{n!}(x-a)^n \\ &\approx f(a) + \frac{f'(a)}{1!}(x-a) \end{aligned}

  • 二阶泰勒展开

f(x)n=02f(n)(a)n!(xa)nf(a)+f(a)1!(xa)+f(a)2!(xa)2\begin{aligned} f(x) &\approx \sum_{n=0}^{2}\frac{f^{(n)}(a)}{n!}(x-a)^n \\ &\approx f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 \end{aligned}

  • N阶泰勒展开

    f(x)n=0Nf(n)(a)n!(xa)nf(a)+f(a)1!(xa)+f(a)2!(xa)2+f(a)3!(xa)3+...\begin{aligned} f(x) &\approx \sum_{n=0}^{N}\frac{f^{(n)}(a)}{n!}(x-a)^n \\ &\approx f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \frac{f'''(a)}{3!}(x-a)^3 + ... \end{aligned}

阶数越大,泰勒展开的值越接近$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

由于损失函数ll中只有唯一变量Hk1(xi)+fk(xi)H_{k-1}(x_i) + f_k(x_i),因此可以将函数简写为l(Hk1(xi)+fk(xi))l(H_{k-1}(x_i) + f_k(x_i))

根据二阶泰勒展开,已知:

f(x)n=02f(n)(a)n!(xa)nf(a)+f(a)1!(xa)+f(a)2!(xa)2\begin{aligned} f(x) &\approx \sum_{n=0}^{2}\frac{f^{(n)}(a)}{n!}(x-a)^n \\ &\approx f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 \end{aligned}

令泰勒展开中的x=Hk1(xi)+fk(xi)x = H_{k-1}(x_i) + f_k(x_i),令泰勒展开中的a=Hk1(xi)a = H_{k-1}(x_i),则(xa)=fk(xi)(x-a) = f_k(x_i)。据此,损失函数l(Hk1(xi)+fk(xi))l(H_{k-1}(x_i) + f_k(x_i))可以被表示为:

l(Hk1(xi)+fk(xi))l(Hk1(xi))+l(Hk1(xi))Hk1(xi)fk(xi)+2l(Hk1(xi))2Hk12(xi)fk2(xi)\begin{aligned} l(H_{k-1}(x_i) + f_k(x_i)) &\approx l(H_{k-1}(x_i)) + \frac{\partial{l(H_{k-1}(x_i))}}{\partial{H_{k-1}(x_i)}} * f_k(x_i) + \frac{\partial^2{l(H_{k-1}(x_i))}}{2\partial{H^2_{k-1}(x_i)}} * f^2_k(x_i)\\ \end{aligned}

在XGBoost中我们定义过损失函数的一阶导数与二阶导数:

gik=l(yi,Hk1(xi))Ht1(xi)g_{ik} = \frac{\partial{l(y_i,H_{k-1}(x_i))}}{\partial{H_{t-1}(x_i)}}

hik=2l(yi,Hk1(xi))Ht12(xi)h_{ik} = \frac{\partial^2{l(y_i,H_{k-1}(x_i))}}{\partial{H^2_{t-1}(x_i)}}

在XGBoost原论文中,为了公式简洁,gig_ihih_i并没有呈现下标kk,但我们已经很清楚:gghh是在每一轮迭代时需要被重新计算的。在这里我们也参照原论文中的做法去掉下标kk。因此,经过泰勒展开后的式子可以化简为:

l(Hk1(xi)+fk(xi))l(Hk1(xi))+gifk(xi)+12hifk2(xi)常数+gifk(xi)+12hifk2(xi)\begin{aligned}l(H_{k-1}(x_i) + f_k(x_i)) &\approx l(H_{k-1}(x_i)) + g_if_k(x_i) + \frac{1}{2}h_if^2_k(x_i) \\ &\approx 常数 + g_if_k(x_i) + \frac{1}{2}h_if^2_k(x_i) \end{aligned}

不难发现,该式子中Hk1(xi)H_{k-1}(x_i)是常数,因此第一部分l(Ht1(xi))l(H_{t-1}(x_i))也是一个常数,常数无法被最小化,因此我们可以将常数从该目标函数中剔除。经过泰勒展开,目标函数变为:

O~k=i=1M(gifk(xi)+12hifk2(xi))+γT+12λj=1Twj2=i=1Mgifk(xi)+12i=1Mhifk2(xi)+γT+12λj=1Twj2\begin{aligned} \tilde{O}_k &= \sum_{i=1}^M\left(g_if_k(x_i) + \frac{1}{2}h_if^2_k(x_i)\right) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^T w_j^2 \\ &= \sum_{i=1}^Mg_if_k(x_i) + \frac{1}{2}\sum_{i=1}^Mh_if^2_k(x_i) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^T w_j^2\end{aligned}

  • 统一自变量

现在目标函数的前两项分别代表所有样本的gifk(xi)g_if_k(x_i)之和,以及所有样本的hifk2(xi)h_if^2_k(x_i)之和乘1/2。别忘记,我们选择的唯一的自变量是wjw_j,因此我们希望能够将fkf_k以某种方式转化为wjw_j。之前已经提到过多次,对任意位于叶子jj上的样本ii而言,数值上fk(xi)=wjf_k(x_i) = w_j,我们可以尝试着从一个样本开始进行转化:

对于单一样本ii,假设这个样本位于叶子jj上,应该有:

gifk(xi)=giwjg_if_k(x_i) = g_iw_j

对于一片叶子jj,我们可以计算这片叶子上所有样本的giwjg_iw_j之和:

ijgiwj\sum_{i \in j} g_iw_j

而一片叶子上所有样本的wjw_j都是一致的,因此一片叶子上的giwjg_iw_j之和可以转变为:

ijgiwj=g1wj + g2wj + ... + gnwj,其中1,2...n是叶子j上的样本=wjijgi\begin{aligned}\sum_{i \in j} g_iw_j &= g_1w_j \ + \ g_2w_j \ + \ ... \ + \ g_nw_j,其中1,2...n是叶子j上的样本 \\ &= w_j\sum_{i \in j} g_i\end{aligned}

假设现在一共有TT片叶子,则整棵树上所有样本的giwjg_iw_j之和为:

j=1T(wjijgi)\sum_{j=1}^T \left( w_j\sum_{i \in j} g_i \right)

所以:

i=1Mgifk(xi)=j=1T(wjijgi)\sum_{i=1}^Mg_if_k(x_i) = \sum_{j=1}^T \left( w_j\sum_{i \in j} g_i \right)

同理,单一样本iihifk2(xi)h_if^2_k(x_i)也可以以相同方式转化。对单一样本:

hifk2(xi)=hiwj2h_if^2_k(x_i) = h_iw^2_j

对一片叶子:

ijhiwj2=h1wj2 + h2wj2 + ... + hnwj2,其中1,2...n是叶子j上的样本=wj2ijhi\begin{aligned}\sum_{i \in j}h_iw^2_j &= h_1w^2_j \ + \ h_2w^2_j \ + \ ... \ + \ h_nw^2_j,其中1,2...n是叶子j上的样本 \\ &= w^2_j\sum_{i \in j} h_i \end{aligned}

对整棵树:

i=1Mhifk2(xi)=j=1T(wj2ijhi)\sum_{i=1}^Mh_if^2_k(x_i) = \sum_{j=1}^T \left( w^2_j\sum_{i \in j} h_i \right)

因此对整个目标函数有:

O~k=i=1Mgifk(xi)+12i=1Mhifk2(xi)+γT+12λj=1Twj2=j=1T(wjijgi+12wj2ijhi)+γT+12λj=1Twj2\begin{aligned} \tilde{O}_k &= \sum_{i=1}^Mg_if_k(x_i) + \frac{1}{2}\sum_{i=1}^Mh_if^2_k(x_i) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^T w_j^2 \\ &=\sum_{j=1}^T \left( w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j\sum_{i \in j} h_i \right) + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^T w_j^2\end{aligned}

不难发现,现在正则项可以与原来损失函数的部分合并了:

=j=1T(wjijgi+12wj2ijhi+12λwj2)+γT=j=1T(wjijgi+12wj2(ijhi+λ))+γT\begin{aligned} &= \sum_{j=1}^T \left( w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j\sum_{i \in j} h_i + \frac{1}{2}\lambda w_j^2 \right) + \gamma T \\ &= \sum_{j=1}^T \left( w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda) \right) + \gamma T\end{aligned}

合并之后,整个目标函数变为两项,一项是所有叶子上的(损失+正则)之和,另一项是叶子总量。现在,我们可以开始求解最小目标函数以及对应的最优自变量wjw_j了。

3 求解XGBoost的目标函数

首先,令目标函数中的叶子总量最小是不可能的,过度降低叶子总量会大幅度伤害模型的学习能力,因此我们只能考虑令所有叶子上的(损失+正则)之和最小。

其次,当树建好之后,叶子与叶子之间是相互独立的,因此每片叶子上的(损失+正则)也是相互独立的。我们只要令每片叶子的(损失+正则)都最小,就可以保证全部叶子的(损失+正则)之和最小。故而,我们要令式子中标注为红色的部分最小:

O~k=j=1T(wjijgi+12wj2(ijhi+λ))+γT\begin{aligned} \tilde{O}_k &= \sum_{j=1}^T \left( \boldsymbol{\color{red}{w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda)}} \right) + \gamma T\end{aligned}

  • 叶子权重wjw_j

将标注为红色的部分命名为μj\mu_j,表示叶子jj上的损失+正则。则有:

μj=wjijgi+12wj2(ijhi+λ)\mu_j = w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda)

现在,对叶子jj而言,在μj\mu_j上对唯一自变量wjw_j求导,则有:

μjwj=wjijgi+12wj2(ijhi+λ)wj=ijgi+wj(ijhi+λ)\begin{aligned}\frac{\partial{\mu_j}}{\partial w_j} &= \frac{\partial{w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda)}}{\partial w_j} \\ \\ &= \sum_{i \in j} g_i + w_j(\sum_{i \in j} h_i + \lambda)\end{aligned}

令一阶导数为0,则有:

ijgi+wj(ijhi+λ)=0wj(ijhi+λ)=ijgiwj=ijgiijhi+λ\begin{aligned} \sum_{i \in j} g_i + w_j(\sum_{i \in j} h_i + \lambda) &= 0 \\ \\ w_j(\sum_{i \in j} h_i + \lambda) &= -\sum_{i \in j} g_i \\ \\ w_j &= -\frac{\sum_{i \in j} g_i}{\sum_{i \in j} h_i + \lambda}\end{aligned}

你应该发现了,对一片叶子来说,令目标函数最小的wjw_j就是我们之前提过的叶子权重,也就是XGBoost数学流程当中叶子上的输出值。如果要令叶子的输出非常接近叶子权重公式,那应该如何拟合每个样本呢?

  • 拟合值

对任意位于叶子jj上的样本ii来说

μi=wjgi+12wj2hi\mu_i = w_jg_i + \frac{1}{2}w^2_jh_i

将一片叶子上的μj\mu_j转变成μi\mu_i时,原则上需要将μj\mu_j中的每一项都转换为单个样本所对应的项,然而在转换正则项时则存在问题:与ijgi\sum_{i \in j} g_i这样可以直接指向单个样本的项不同,λ\lambda是针对与一片叶子设置的值,如果要将λ\lambda转变为针对单一样本的正则项,则需要知道当前叶子上一共有多少样本。然而,拟合发生在建树之前,因此在这一时间点不可能知道一片叶子上的样本总量,因此在xgboost的实际实现过程当中,拟合每一片叶子时不涉及正则项,只有在计算结构分数与叶子输出值时才使用正则项。

μi\mu_i上唯一的自变量wjw_j求导,则有:

μiwj=(wjgi+12wj2hi)wj=gi+wjhi\begin{aligned}\frac{\partial{\mu_i}}{\partial w_j} &= \frac{\partial{\left( w_jg_i + \frac{1}{2}w^2_jh_i \right)}}{\partial w_j} \\ \\ &= g_i + w_jh_i\end{aligned}

令一阶导数为0,则有:

gi+wjhi=0wjhi=giwj=gihi\begin{aligned} g_i + w_jh_i &= 0 \\ \\ w_jh_i &= - g_i \\ \\ w_j &= -\frac{g_i}{h_i} \end{aligned}

对任意样本ii而言,令目标函数最小的最优wjw_j就是我们的伪残差rir_i,也就是XGBoost数学流程当中用于进行拟合的拟合值

  • 结构分数

现在,我们把令目标函数最小的最优wjw_j带回到μj\mu_j中,查看μj\mu_j如何变化:

μj=wjijgi+12wj2(ijhi+λ)=ijgiijhi+λijgi+12(ijgiijhi+λ)2ijhi+λ=(ijgi)2ijhi+λ+12(ijgi)2ijhi+λ=12(ijgi)2ijhi+λ\begin{aligned} \mu_j &= w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda) \\ &= -\frac{\sum_{i \in j} g_i}{\sum_{i \in j} h_i + \lambda} * \sum_{i \in j} g_i + \frac{1}{2}(-\frac{\sum_{i \in j} g_i}{\sum_{i \in j} h_i + \lambda})^2 * {\sum_{i \in j} h_i + \lambda}\\ &= -\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda} + \frac{1}{2}\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda} \\ &= - \frac{1}{2}\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda} \end{aligned}

因此,目标函数(所有叶子上的损失)就可以变为:

O~k=j=1T(wjijgi+12wj2(ijhi+λ))+γT=j=1T(12(ijgi)2ijhi+λ)+γT\begin{aligned} \tilde{O}_k &= \sum_{j=1}^T \left( \boldsymbol{\color{red}{w_j\sum_{i \in j} g_i + \frac{1}{2}w^2_j(\sum_{i \in j} h_i + \lambda)}} \right) + \gamma T \\ \\ &= \sum_{j=1}^T \left( -\frac{1}{2}\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda} \right) + \gamma T \end{aligned}

因此,一片叶子上的目标函数就是:

Oj=12(ijgi)2ijhi+λ+γO_j = -\frac{1}{2}\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda} + \gamma

对任意一片叶子来说,目标函数可以衡量叶子的质量,其中γ\gamma是可以设定的超参数,12\frac{1}{2}为常数,因此对任意叶子,我们希望标注为红色的部分越小越好:

Oj=12((ijgi)2ijhi+λ)+γO_j = \frac{1}{2}\left( \boldsymbol{\color{red}{-\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda}}} \right)+ \gamma

故而,我们希望以下式子越大越好:

(ijgi)2ijhi+λ\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda}

这个式子,正是XGBoost用于分枝时的指标“结构分数”(Structure Score)

  • 结构分数的增益

当分枝的时候,我们希望目标函数越小越好,因此在分枝过程中,父节点的目标函数是大于子节点的目标函数的,因此我们可以使用(父节点目标函数 - 子节点目标函数之和)来衡量分枝的质量,则有:

Gain=Op(Ol+Or)=12(iPgi)2iPhi+λ+γ(12(iLgi)2iLhi+λ+γ12(iRgi)2iRhi+λ+γ)=12(iPgi)2iPhi+λ+γ+12(iLgi)2iLhi+λγ+12(iRgi)2iRhi+λγ=12((iLgi)2iLhi+λ+(iRgi)2iRhi+λ(iPgi)2iPhi+λ)γ=12(ScoreL+ScoreRScoreP)γ\begin{aligned} Gain &= O_p - (O_l + O_r) \\ \\ &= -\frac{1}{2}\frac{(\sum_{i \in P} g_i)^2}{\sum_{i \in P} h_i + \lambda} + \gamma - (-\frac{1}{2}\frac{(\sum_{i \in L} g_i)^2}{\sum_{i \in L} h_i + \lambda} + \gamma -\frac{1}{2}\frac{(\sum_{i \in R} g_i)^2}{\sum_{i \in R} h_i + \lambda} + \gamma) \\ \\ &= -\frac{1}{2}\frac{(\sum_{i \in P} g_i)^2}{\sum_{i \in P} h_i + \lambda} + \gamma + \frac{1}{2}\frac{(\sum_{i \in L} g_i)^2}{\sum_{i \in L} h_i + \lambda} - \gamma + \frac{1}{2}\frac{(\sum_{i \in R} g_i)^2}{\sum_{i \in R} h_i + \lambda} - \gamma \\ \\ &= \frac{1}{2}\left( \frac{(\sum_{i \in L} g_i)^2}{\sum_{i \in L} h_i + \lambda} + \frac{(\sum_{i \in R} g_i)^2}{\sum_{i \in R} h_i + \lambda} - \frac{(\sum_{i \in P} g_i)^2}{\sum_{i \in P} h_i + \lambda} \right) - \gamma \\ \\ &= \frac{1}{2} (Score_L + Score_R - Score_P) - \gamma \end{aligned}

其中,γ\gamma是可以设定的超参数,12\frac{1}{2}为常数,因此:

Gain=ScoreL+ScoreRScorePGain = Score_L + Score_R - Score_P

这就是我们在分枝时所使用的结构分数增益了

现在你发现了,XGBoost流程中所使用的全部新公式(包括独特的拟合值、独特的分枝指标、独特的输出值)都是通过令目标函数最小而求解出来的。因此,XGBoost整个流程就保证了目标函数一定是向着最小化方向进行迭代的,新生成的每片叶子上的输出值wjw_j都是会令目标函数最小化的输出值。现在,你可以回答以下的问题了:

    1. 建树时拟合的rik=gikhikr_{ik} = -\frac{g_{ik}}{h_{ik}}究竟是什么?拟合它有什么意义?
    1. 结构分数和结构分数增益的公式是如何推导出来的?为什么这样建树可以提升模型的效果?
    1. 为什么叶子节点的输出值wjw_j(ijgik)ijhik+λ-\frac{(\sum_{i \in j} g_{ik})}{\sum_{i \in j} h_{ik} + \lambda}?这样输出有什么意义?
    1. 课程的第一部分说XGBoost拟合的也是残差,残差在哪里?

唯一余留问题4其实我们早已解答过,回去查看第二部分中的《2.2 弱评估器的分枝》,你会发现其实我们早已经求出残差了。当目标函数为12MSE\frac{1}{2}MSE,负梯度gi-g_i就等于残差,而hi=1h_i = 1,因此拟合项gihi-\frac{g_i}{h_i}自然也是残差本身了。因此,XGBoost也是拟合负梯度的算法,并且在特定损失函数下,XGBoost也拟合残差。在了解这个推导过程之后,再返回复习XGBoost的整体数学流程,你会发现数学真的让算法变得更简单,而不是更复杂。如果你已经学到了这里,你已经跨过了机器学习原理中几乎最难的一部分。在接下来的课程中,所有数学推导对你来说都不再是阻碍了。