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)k−1f(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)=argminC∑i=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,Hk−1(xi))∂Hk−1(xi)g_{ik} = \frac{\partial{l(y_i,H_{k-1}(x_i))}}{\partial{H_{k-1}(x_i)}}

    hik=∂2l(yi,Hk−1(xi))∂Hk−12(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)求导，再令求导后的结果等于Ht−1(xi)H_{t-1}(x_i)的值，并不是直接对Ht−1(xi)H_{t-1}(x_i)这一常数求导。


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

l=x2+xl = x^2+x

对常数0求导：

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

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

∂l∂x=∂(x2+x)∂x=2x+1=2∗0+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)=Hk−1  (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时，所有求导计算过程中的Hk−1(xi)H_{k-1}(x_i)都等于初始H0(x)H_0(x)，在k>1时，每个样本上的Hk−1(xi)H_{k-1}(x_i)都是不同的取值。

  • 4) 建树


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

    Scorej=(∑i∈jgi)2∑i∈jhi+λScore_j = \frac{(\sum_{i \in j}g_i)^2}{\sum_{i \in j}h_i + \lambda}

    Gain=12((∑i∈Lgi)2∑i∈Lhi+λ+(∑i∈Rgi)2∑i∈Rhi+λ−(∑i∈Pgi)2∑i∈Phi+λ)−γ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为输出值。


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


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


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

  • 6) 迭代


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

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

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

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

    对整个算法则有：

    Hk(x)=Hk−1(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是−(∑i∈jgik)∑i∈jhik+λ-\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)=Hk−1(xi)+fk(xi)H_k(x_i) = H_{k-1}(x_i) + f_k(x_i)

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

当迭代到第kk次时，损失函数中的yiy_i与Hk−1(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,Hk−1(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=0∞f(n)(a)n!(x−a)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!(x−a)n≈f(a)+f′(a)1!(x−a)\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!(x−a)n≈f(a)+f′(a)1!(x−a)+f′′(a)2!(x−a)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!(x−a)n≈f(a)+f′(a)1!(x−a)+f′′(a)2!(x−a)2+f′′′(a)3!(x−a)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中只有唯一变量Hk−1(xi)+fk(xi)H_{k-1}(x_i) + f_k(x_i)，因此可以将函数简写为l(Hk−1(xi)+fk(xi))l(H_{k-1}(x_i) + f_k(x_i))。

根据二阶泰勒展开，已知：

f(x)≈∑n=02f(n)(a)n!(x−a)n≈f(a)+f′(a)1!(x−a)+f′′(a)2!(x−a)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=Hk−1(xi)+fk(xi)x = H_{k-1}(x_i) + f_k(x_i)，令泰勒展开中的a=Hk−1(xi)a = H_{k-1}(x_i)，则(x−a)=fk(xi)(x-a) = f_k(x_i)。据此，损失函数l(Hk−1(xi)+fk(xi))l(H_{k-1}(x_i) + f_k(x_i))可以被表示为：

l(Hk−1(xi)+fk(xi))≈l(Hk−1(xi))+∂l(Hk−1(xi))∂Hk−1(xi)∗fk(xi)+∂2l(Hk−1(xi))2∂Hk−12(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,Hk−1(xi))∂Ht−1(xi)g_{ik} = \frac{\partial{l(y_i,H_{k-1}(x_i))}}{\partial{H_{t-1}(x_i)}}

hik=∂2l(yi,Hk−1(xi))∂Ht−12(xi)h_{ik} = \frac{\partial^2{l(y_i,H_{k-1}(x_i))}}{\partial{H^2_{t-1}(x_i)}}

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

l(Hk−1(xi)+fk(xi))≈l(Hk−1(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}

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

O~k=∑i=1M(gifk(xi)+12hifk2(xi))+γT+12λ∑j=1Twj2=∑i=1Mgifk(xi)+12∑i=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之和：

∑i∈jgiwj\sum_{i \in j} g_iw_j

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

∑i∈jgiwj=g1wj + g2wj + ... + gnwj，其中1,2...n是叶子j上的样本=wj∑i∈jgi\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(wj∑i∈jgi)\sum_{j=1}^T \left( w_j\sum_{i \in j} g_i \right)

所以：

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

同理，单一样本ii的hifk2(xi)h_if^2_k(x_i)也可以以相同方式转化。对单一样本：

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

对一片叶子：

∑i∈jhiwj2=h1wj2 + h2wj2 + ... + hnwj2，其中1,2...n是叶子j上的样本=wj2∑i∈jhi\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(wj2∑i∈jhi)\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)+12∑i=1Mhifk2(xi)+γT+12λ∑j=1Twj2=∑j=1T(wj∑i∈jgi+12wj2∑i∈jhi)+γ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(wj∑i∈jgi+12wj2∑i∈jhi+12λwj2)+γT=∑j=1T(wj∑i∈jgi+12wj2(∑i∈jhi+λ))+γ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(wj∑i∈jgi+12wj2(∑i∈jhi+λ))+γ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=wj∑i∈jgi+12wj2(∑i∈jhi+λ)\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求导，则有：

∂μj∂wj=∂wj∑i∈jgi+12wj2(∑i∈jhi+λ)∂wj=∑i∈jgi+wj(∑i∈jhi+λ)\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，则有：

∑i∈jgi+wj(∑i∈jhi+λ)=0wj(∑i∈jhi+λ)=−∑i∈jgiwj=−∑i∈jgi∑i∈jhi+λ\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中的每一项都转换为单个样本所对应的项，然而在转换正则项时则存在问题：与∑i∈jgi\sum_{i \in j} g_i这样可以直接指向单个样本的项不同，λ\lambda是针对与一片叶子设置的值，如果要将λ\lambda转变为针对单一样本的正则项，则需要知道当前叶子上一共有多少样本。然而，拟合发生在建树之前，因此在这一时间点不可能知道一片叶子上的样本总量，因此在xgboost的实际实现过程当中，拟合每一片叶子时不涉及正则项，只有在计算结构分数与叶子输出值时才使用正则项。

对μi\mu_i上唯一的自变量wjw_j求导，则有：

∂μi∂wj=∂(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=wj∑i∈jgi+12wj2(∑i∈jhi+λ)=−∑i∈jgi∑i∈jhi+λ∗∑i∈jgi+12(−∑i∈jgi∑i∈jhi+λ)2∗∑i∈jhi+λ=−(∑i∈jgi)2∑i∈jhi+λ+12(∑i∈jgi)2∑i∈jhi+λ=−12(∑i∈jgi)2∑i∈jhi+λ\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(wj∑i∈jgi+12wj2(∑i∈jhi+λ))+γT=∑j=1T(−12(∑i∈jgi)2∑i∈jhi+λ)+γ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(∑i∈jgi)2∑i∈jhi+λ+γ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(−(∑i∈jgi)2∑i∈jhi+λ)+γ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

故而，我们希望以下式子越大越好：

(∑i∈jgi)2∑i∈jhi+λ\frac{(\sum_{i \in j} g_i)^2}{\sum_{i \in j} h_i + \lambda}

这个式子，正是XGBoost用于分枝时的指标“结构分数”（Structure Score）。

  • 结构分数的增益

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

Gain=Op−(Ol+Or)=−12(∑i∈Pgi)2∑i∈Phi+λ+γ−(−12(∑i∈Lgi)2∑i∈Lhi+λ+γ−12(∑i∈Rgi)2∑i∈Rhi+λ+γ)=−12(∑i∈Pgi)2∑i∈Phi+λ+γ+12(∑i∈Lgi)2∑i∈Lhi+λ−γ+12(∑i∈Rgi)2∑i∈Rhi+λ−γ=12((∑i∈Lgi)2∑i∈Lhi+λ+(∑i∈Rgi)2∑i∈Rhi+λ−(∑i∈Pgi)2∑i∈Phi+λ)−γ=12(ScoreL+ScoreR−ScoreP)−γ\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+ScoreR−ScorePGain = Score_L + Score_R - Score_P

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

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

    1. 建树时拟合的rik=−gikhikr_{ik} = -\frac{g_{ik}}{h_{ik}}究竟是什么？拟合它有什么意义？
    1. 结构分数和结构分数增益的公式是如何推导出来的？为什么这样建树可以提升模型的效果？
    1. 为什么叶子节点的输出值wjw_j 是−(∑i∈jgik)∑i∈jhik+λ-\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的整体数学流程，你会发现数学真的让算法变得更简单，而不是更复杂。

这里需要提醒：本科期间在CS专业学到的数学已经够用了，需要什么，学什么即可。如果你对其它的数学感兴趣，可以自己下去花时间钻研。但是若因为钻研纯数而放弃了CS的学习，到头来竹篮打水，这是笔者不建议的行为。