https://xgboost.readthedocs.io/en/latest//tutorials/index.html
https://arxiv.org/abs/1603.02754
主要参考了上面两个文献,第一个是XGBoost的官方文档,第二个是XGBoost的论文原文。
树提升算法回顾
- 学习目标的正则化
通常来讲,树提升算法是串行地搭建若干的树,每个树的训练会考虑到前边树的训练结果,最终把所有树的计算结果加在一起得到预测值。
y ^ i = ϕ ( x i ) = ∑ k = 1 K f k ( x i ) \hat y_i=\phi(x_i)=\sum _{k=1} ^K f_k(x_i) y^i=ϕ(xi)=k=1∑Kfk(xi)
其中, f k f_k fk表示第 k k k个树
同时,为了缓解过拟合问题,我们在定义损失函数的同时也会添加一定的正则项
L = ∑ i l ( y ^ i , y i ) + ∑ k Ω ( f k ) w h e r e Ω ( f ) = γ T + 1 2 λ ∣ ∣ w ∣ ∣ 2 L=\sum_i l(\hat y_i,y_i)+\sum_k\Omega(f_k) \\ where \ \Omega(f) =\gamma T + \frac 1 2 \lambda ||w||^2 L=i∑l(y^i,yi)+k∑Ω(fk)where Ω(f)=γT+21λ∣∣w∣∣2
T表示树的node个数, w w w表示叶子结点的值
- 梯度提升树
梯度提升树中,每个树的训练目标是前边的树留下的残差。我们的损失函数/优化目标可以写为
L ( t ) = ∑ i = 1 n l ( y i , y ^ i ( t − 1 ) + f t ( x i ) ) + Ω ( f t ) L^{(t)}=\sum_{i=1}^nl(y_i,\hat y_i ^{(t-1)}+f_t(x_i))+\Omega(f_t) L(t)=i=1∑nl(yi,y^i(t−1)+ft(xi))+Ω(ft)
我们对这个式子做二阶的泰勒展开,得到
L ( t ) ≃ ∑ i = 1 n l ( y i , y ^ i ( t − 1 ) ) + g i f t ( x i ) + 1 2 h i f t ( x i ) 2 + Ω ( f t ) L^{(t)} \simeq\sum_{i=1}^nl(y_i,\hat y_i^{(t-1)})+g_if_t(x_i)+ \frac 1 2 h_if_t(x_i)^2+\Omega(f_t) L(t)≃i=1∑nl(yi,y^i(t−1))+gift(xi)+21hift(xi)2+Ω(ft)
剔除优化目标中的常数项,
L ( t ) ≃ ∑ i = 1 n g i f t ( x i ) + 1 2 h i f t ( x i ) 2 + Ω ( f t ) L^{(t)} \simeq\sum_{i=1}^ng_if_t(x_i)+ \frac 1 2 h_if_t(x_i)^2+\Omega(f_t) L(t)≃i=1∑ngift(xi)+21hift(xi)2+Ω(ft)
我们定义 I j = { i ∣ q ( x i ) = j } I_j=\{i|q(\mathbf x_i)=j\} Ij={i∣q(xi)=j}作为被划分到叶子结点 j j j的样本index集合,那么我们可以把我们的优化目标重新写为
L ( t ) = ∑ j [ ( ∑ i ∈ I j g i ) w j + 1 2 ( ∑ i ∈ I j h i ) w j 2 ] + γ T + 1 2 λ w j 2 = ∑ j [ ( ∑ i ∈ I j g i ) w j + 1 2 ( ∑ i ∈ I j h i + λ ) w j 2 ] + γ T \begin{aligned} L^{(t)} &=\sum_{j}[ ( \sum_{i\in I_j}g_i )w_j + \frac 1 2 ( \sum_{i\in I_j}h_i )w_j^2]+ \gamma T + \frac 1 2 \lambda w_j^2 \\ &=\sum_{j}[ ( \sum_{i\in I_j}g_i )w_j + \frac 1 2 ( \sum_{i\in I_j}h_i +\lambda )w_j^2]+ \gamma T \end{aligned} L(t)=j∑[(i∈Ij∑gi)wj+21(i∈Ij∑hi)wj2]+γT+21λwj2=j∑[(i∈Ij∑gi)wj+21(i∈Ij∑hi+λ)wj2]+γT
对于具体的一个 j j j与 w j w_j wj来讲,上边这个优化目标是个二次函数
L j ( t ) = ( ∑ i ∈ I j g i ) w j + 1 2 ( ∑ i ∈ I j h i + λ ) w j 2 \begin{aligned} L^{(t)}_j = ( \sum_{i\in I_j}g_i )w_j + \frac 1 2 ( \sum_{i\in I_j}h_i +\lambda )w_j^2 \end{aligned} Lj(t)=(i∈Ij∑gi)wj+21(i∈Ij∑hi+λ)wj2
求导置零即可得到极值点,#损失函数一般是凸函数,所以二阶导大于等于0,因此 h i ≥ 0 h_i\geq 0 hi≥0,而正则项通常大于零,所以这是个开口朝上的二次函数,极值点即为最小值
w j ∗ = − ∑ i ∈ I j g i ∑ i ∈ I j h i + λ w_j^* =- \frac {\sum_{i\in I_j}g_i } {\sum_{i\in I_j}h_i +\lambda } wj∗=−∑i∈Ijhi+λ∑i∈Ijgi
我们把极值点的计算结果返回到优化目标中,便可以用 g i , h i g_i,h_i gi,hi表示出一棵树最小值。
L ( t ) ( q ) = − 1 2 ∑ j ( ∑ i ∈ I j g i ) 2 ∑ i ∈ I j h i + λ + γ T L^{(t)}(q)= -\frac 1 2 \sum_{j} \frac { {(\sum_{i\in I_j}g_i })^2} {\sum_{i\in I_j}h_i +\lambda } +\gamma T L(t)(q)=−21j∑∑i∈Ijhi+λ(∑i∈Ijgi)2+γT
我们为啥要这么费劲儿地用 g i , h i g_i,h_i gi,hi表达最小值呢?考虑到损失函数在一定程度上反映着叶子结点的优劣,loss越低越好,于是我们可以把损失函数的值当作结点分裂的判定依据。
L s p l i t = 1 2 [ ( ∑ i ∈ I L g i ) 2 ∑ i ∈ I L h i + λ + ( ∑ i ∈ I R g i ) 2 ∑ i ∈ I R h i + λ − ( ∑ i ∈ I g i ) 2 ∑ i ∈ I h i + λ ] − γ L_{split}=\frac 1 2 [ \frac { {(\sum_{i\in I_L}g_i })^2} {\sum_{i\in I_L}h_i +\lambda } +\frac { {(\sum_{i\in I_R}g_i })^2} {\sum_{i\in I_R}h_i +\lambda } -\frac { {(\sum_{i\in I}g_i })^2} {\sum_{i\in I}h_i +\lambda } ] -\gamma Lsplit=21[∑i∈ILhi+λ(∑i∈ILgi)2+∑i∈IRhi+λ(∑i∈IRgi)2−∑i∈Ihi+λ(∑i∈Igi)2]−γ
这里,我们用父结点的loss减去分裂后的两个叶子结点的loss得到了划分判定。#注意其中的 γ \gamma γ正则项 - 缩减与列降采样
决策树太深的时候很容易过拟合。前人们总结了很多的防止过拟合的措施。
一种是像梯度下降一样设置一个学习率 η \eta η
另一种是对属性进行降采样,像随机森林
划分查找算法
- Basic Exact Greedy Algorithm
在比较小的数据集中,对于某个具体的特征,我们可以把样本进行排序,然后逐个样本地去搜索最优的划分点。 - Approximate Algorithm
由于上一个算法要逐点确认,在比较大的数据集中会造成很大的时间开销,因此我们可以根据分位数来预先确定几个候选划分点,然后逐个检验它们的 L s p l i t L_{split} Lsplit
- Weighted Quantile Sketch
这里有了一个新的问题,就是如何去产生这些划分点,我们重新考察损失函数
L ( t ) ≃ ∑ i = 1 n g i f t ( x i ) + 1 2 h i f t ( x i ) 2 + Ω ( f t ) L^{(t)} \simeq\sum_{i=1}^ng_if_t(x_i)+ \frac 1 2 h_if_t(x_i)^2+\Omega(f_t) L(t)≃i=1∑ngift(xi)+21hift(xi)2+Ω(ft)
将其改写为
∑ i = 1 n 1 2 h i ( f ( x i ) + g i h i ) 2 + Ω ( f t ) + c o n s t \sum_{i=1}^n \frac 1 2 h_i(f(x_i) + \frac {g_i} {h_i})^2+\Omega(f_t)+const i=1∑n21hi(f(xi)+higi)2+Ω(ft)+const
可以将其看作用 h i h_i hi加权的平方和,因此 h i h_i hi也代表了每个样本的权重值
定义
r k ( z ) = 1 ∑ ( x , h ) h ∑ ( x , h ) ∈ D k , x < z h r_k(z)=\frac 1 {\sum_{(x,h)}h} \sum_{(x,h)\in D_k, x\lt z} h rk(z)=∑(x,h)h1(x,h)∈Dk,x<z∑h
表示在特征 k k k上的值小于 z z z的权重占比。
我们的目标是在找到满足如下关系的划分点集合 { s k 1 , s k 2 , . . . , s k l } \{s_{k1},s_{k2},...,s_{kl}\} {sk1,sk2,...,skl},
∣ r k ( s k , j ) − r k ( s k , j + 1 ) ∣ < ϵ , s k 1 = min i x i k , s k l = max i x i k |r_k(s_{k,j})-r_k(s_{k,j+1})|<\epsilon,s_{k1}=\min_i \mathbf x_{ik},s_{kl}=\max_i \mathbf x_{ik} ∣rk(sk,j)−rk(sk,j+1)∣<ϵ,sk1=iminxik,skl=imaxxik
#没搞清楚这里为啥是小于号 - Sparsity-aware Split Finding
在实际任务中,经常会出现一些缺省的情况。对于缺省的属性值,XGBoost会把选择一个分支作为默认分支,将缺省该属性值的样本点划分过去。
系统设计
- Column Block for Parallel Learning
对于一棵树的生成过程,最耗时的在于划分结点时对样本进行按照属性值大小的排序。
为了缓解这个问题,我们提前对数据集进行排序处理,存放在内存中,称为block
然后我们对每个属性值顺序读取计算,即得到划分的计算结果。
对于近似算法,block也有很大的作用。我们用多个block来存放数据,每个block存放数据集的若干行。这样分位数的查找就编程了一个线性的扫描过程。对每一列(每个特征)的统计也可以并行地进行,从而提高效率。
时间复杂度分析: 假设树的最大深度是 d d d,树的个数是 K K K,样本个数是 n n n- Exact Greedy Algorithm中,对于第m层,我们对其中的一个结点进行排序的复杂度是
O
(
n
/
2
m
−
1
log
(
n
/
2
m
−
1
)
)
O(n/2^{m-1}\log(n/2^{m-1}))
O(n/2m−1log(n/2m−1)),所有结点加在一起,变为
O
(
n
log
(
n
/
2
m
−
1
)
)
O(n\log(n/2^{m-1}))
O(nlog(n/2m−1)),通常来讲
2
m
−
1
2^{m-1}
2m−1与
n
n
n相比较小,所以我们可以近似认为复杂度为
O
(
n
log
(
n
)
)
O(n\log(n))
O(nlog(n))。考虑到缺省值,排序的复杂度改写为
O
(
x
log
(
n
)
)
O(x\log(n))
O(xlog(n)),其中
x
x
x表示未缺省的样本数量。再考虑到树的深度与个数,复杂度就变为
O
(
K
d
x
log
(
n
)
)
O(Kdx\log(n))
O(Kdxlog(n))。
如果我们使用块结构,我们首先需要对数据进行排序,复杂度为 O ( x log ( n ) ) O(x\log(n)) O(xlog(n)),后续的划分操作计算则是线性复杂度 O ( x ) O(x) O(x),那么总的复杂度为 O ( K d x + x log ( n ) ) O(Kdx+x\log(n)) O(Kdx+xlog(n)) - Approximate Algorithm中,对于第
m
m
m层中的一个结点,给出的一组划分点,我们需要找出该结点中的每个样本对应的属性值在哪个区间中,这可以用二分查找完成,复杂度为
O
(
x
/
(
2
m
−
1
)
log
q
)
O(x/(2^{m-1})\log q)
O(x/(2m−1)logq),其中
q
q
q是候选的划分点的个数,对于整个层,则为
O
(
x
log
q
)
O(x\log q)
O(xlogq)。具体实现中我们可以采取桶的思想,
q
q
q个分裂点对应着
q
+
1
q+1
q+1个区间,我们每次做二分查找,找出样本应该放在哪个区间中,最终再按照Algorithm 2进行统计。这样,对于整个树,我们的算法复杂度为
O
(
K
d
x
log
q
)
O(Kdx\log q)
O(Kdxlogq)
如果我们使用了block结构,那么每个样本不需要进行二分的查找,而是只需要线性扫描即可。初始化时,我们只需要在块内进行初始化,那么我们的复杂度可以写为 O ( K d x + x log B ) O(Kdx+x\log B) O(Kdx+xlogB), B B B为最大的块的行数,也即我们在块内进行了排序。
- Exact Greedy Algorithm中,对于第m层,我们对其中的一个结点进行排序的复杂度是
O
(
n
/
2
m
−
1
log
(
n
/
2
m
−
1
)
)
O(n/2^{m-1}\log(n/2^{m-1}))
O(n/2m−1log(n/2m−1)),所有结点加在一起,变为
O
(
n
log
(
n
/
2
m
−
1
)
)
O(n\log(n/2^{m-1}))
O(nlog(n/2m−1)),通常来讲
2
m
−
1
2^{m-1}
2m−1与
n
n
n相比较小,所以我们可以近似认为复杂度为
O
(
n
log
(
n
)
)
O(n\log(n))
O(nlog(n))。考虑到缺省值,排序的复杂度改写为
O
(
x
log
(
n
)
)
O(x\log(n))
O(xlog(n)),其中
x
x
x表示未缺省的样本数量。再考虑到树的深度与个数,复杂度就变为
O
(
K
d
x
log
(
n
)
)
O(Kdx\log(n))
O(Kdxlog(n))。
- Cache-aware Access
从上边的分析中,我们可以看出,在block中对属性值进行排序能够大大地降低复杂度。但是,如果我们对每一个属性列进行了排序,那么 g i , h i g_i,h_i gi,hi就会和属性值不对应,参考下图。这种排序导致我们对 g i , h i g_i,h_i gi,hi的访问变得很麻烦。
假设我们机器的内存比较有限,每次只能读取若干行数据集,那么就可能导致上图中的 g [ p t r [ i ] ] g[ptr[i]] g[ptr[i]]不在内存中,出现cache miss,进而使得算法变慢。为了解决这个问题,我们先把 g i , h i g_i,h_i gi,hi预读取到内存当中(如果装不下的话则适当多地读取),而对于属性值,我们采取mini-batch的方式逐次读取,这样就能有效地提升命中率。
基于这种方法,原文进行了进一步的分析。当block size比较大时,上述的方法仍会有比较多的cache miss(因为 g g g和 h h h不能全读进去);而当block size比较小时,则会导致每个线程的负载太小,从而并行化的效率不高。 - Blocks for Out-of-core Computation
在计算过程中,单独地开辟一个线程用于将磁盘上的block读取到内存中。原文用两种方法改进系统I/O- Block Compression. 对数据块进行压缩,读取后进行解压缩
- Block Sharding. 把block分到多个磁盘上,每个磁盘会分配一个线程进行读取,当手里的磁盘比较多的时候可以用这个方法来提升吞吐率。