XGboost: A Scalable Tree Boosting System论文及源码导读
这篇论文一作为陈天齐,XGBoost是从竞赛pk中脱颖而出的算法,目前开源在github,和传统gbdt方式不同,XGBoost对loss function进行了二阶的泰勒展开,并增加了正则项,用于权衡目标函数的下降和模型的复杂度[12]。罗列下优势:
- 可扩展性强
- 为稀疏数据设计的决策树训练方法
- 理论上得到验证的加权分位数略图法
- 并行和分布式计算
- 设计高效核外计算,进行cache-aware数据块处理
分布式训练树模型boosting方法已有[1,2,3]。
整体目标
L(ϕ)=∑il(yi,ˆyi)+∑kΩ(fk)L(ϕ)=∑il(yi,^yi)+∑kΩ(fk)其中L(⋅)L(⋅)为目标函数,l(⋅)l(⋅)是损失函数,通常是凸函数,用于 刻画预测值ˆyi^yi和真实值yiyi的差异,第二项Ω(⋅)Ω(⋅)为模型的正则化项, 用于降低模型的复杂度,减轻过拟合问题,类似的正则化方法可以在引文[4]看到。模型目标是最小化目标函数。
L(⋅)L(⋅)为函数空间上的表达,我们可以将其转换为下面这张gradient boosting的方式,记ˆy(t)i^y(t)i为第ii个样本第tt轮迭代:
L(t)=n∑i=1l(yi,ˆyi(t−1)+ft(xi))+Ω(ft)L(t)=n∑i=1l(yi,^yi(t−1)+ft(xi))+Ω(ft)
对该函数在ˆy(t)i^y(t)i位置进行二阶泰勒展开,可以加速优化过程,我们得到目标函数的近似:
L(t)≃n∑i=1[l(yi,ˆy(t−1))+gift(xi)+12hif2t(xi)]+Ω(ft)L(t)≃n∑i=1[l(yi,^y(t−1))+gift(xi)+12hif2t(xi)]+Ω(ft)
泰勒展开的推导部分,可以参考思考1,其中第一项是常数项,删除可得:
˜L(t)=n∑i=1[gift(xi)+12hif2t(xi)]+Ω(ft)−−(1)~L(t)=n∑i=1[gift(xi)+12hif2t(xi)]+Ω(ft)−−(1)
下面对正则化项进行参数化定义。延续前文GBDT的概念,引入分裂节点jj定义的区域记作Ij={i|q(xi)=j}Ij={i|q(xi)=j}
那么ΩΩ展开原目标函数改写为:
˜L(t)=n∑i=1[gift(xi)+12hif2t(xi)]+γT+12λT∑j=1w2j=T∑j=1[(∑i∈Ijgi)wj+12(∑i∈Ijhi+λ)w2j]+γT−−(2)~L(t)=n∑i=1[gift(xi)+12hif2t(xi)]+γT+12λT∑j=1w2j=T∑j=1⎡⎣⎛⎝∑i∈Ijgi⎞⎠wj+12⎛⎝∑i∈Ijhi+λ⎞⎠w2j⎤⎦+γT−−(2)
对于固定的树结构q(x)q(x),对wjwj求导得解析解w∗jw∗j:
w∗j=∑i∈Ijgi∑i∈Ijhi+λw∗j=∑i∈Ijgi∑i∈Ijhi+λ
代入到(1)(1)式中,可得
˜L(t)(q)=−12(∑i∈Ijgi)2∑i∈Ijhi+λ+γT−−(3)~L(t)(q)=−12(∑i∈Ijgi)2∑i∈Ijhi+λ+γT−−(3)
公式(3)(3)可以作为分裂节点的打分,形式上很像CART树纯度打分的计算,区别在于它是从目标函数中推导而得。图中显示目标函数值的计算:
实践中,很难去穷举每一颗树进行打分,再选出最好的。通常采用贪心的方式,逐层选择最佳的分裂节点。假设ILIL和IRIR为分裂节点的左右节点,记I=IL∪IRI=IL∪IR。
则选择此节点分裂的收益为:
Lsplit=12[(∑i∈ILgi)2∑i∈ILhi+λ+(∑i∈IJgi)2∑i∈IJhi+λ−(∑i∈Igi)2∑i∈Ihi+λ]−γ−−(4)Lsplit=12⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣(∑i∈ILgi)2∑i∈ILhi+λ+(∑i∈IJgi)2∑i∈IJhi+λ−(∑i∈Igi)2∑i∈Ihi+λ⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦−γ−−(4)
[补充]
1.作为GB类方法,也可以采用shrinkage策略
2.随机森林[5]的特征降采样(subsampling)克服过拟合代码,xgboost用了类似的技术。和传统的降采样不同,xgboost按列进行降采样,在并行化有加速作用。(待研究)
查找分裂节点(split finding)
贪心算法
贪心算法是最基本的方法,前面介绍的时候有提过。具体做法:遍历所有特征中可能的分裂点位置,根据公式(4)(4)找到最合适的位置。Split_Finding()Split_Finding()算法流程如下
Split_Finding():Input:I,instancesetofcurrentnodeInput:d,featuredimensiongain←0G←∑i∈Igi,H←∑i∈Ihifork=1tomdoGL←0,HL←0forjinsorted(I,byxjk)doGL←GL+gj,HL←HL+hjGR←G−GL,HR←H−HLscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endendOutput:split_valuewithmaxscoreSplit_Finding():Input:I,instancesetofcurrentnodeInput:d,featuredimensiongain←0G←∑i∈Igi,H←∑i∈Ihifork=1tomdoGL←0,HL←0forjinsorted(I,byxjk)doGL←GL+gj,HL←HL+hjGR←G−GL,HR←H−HLscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endendOutput:split_valuewithmaxscore
近似算法
如果数据不能一次读入内存,使用贪心算法效率较低。近似算法在过去也有应用[6, 7, 8]。具体描述为,对于某个特征xkxk,找到该特征若干值域分界点{sk1,sk2,...,skl}{sk1,sk2,...,skl}。根据特征的值对样本进行分桶,对于每个桶内的样本统计值GG、HH进行累加(两个统计量含义同贪心算法),记为分界点vv的统计量,vv满足{xkj=skv}{xkj=skv}。最后在分界点集合上调用Split_Finding()Split_Finding()进行贪心查找,得到的结果为最佳分裂点的近似。具体如下:
fork=1tomdoProposeSk={sk1,sk2,...,skl}bypercentileonfeaturekProposecanbedonepertree(global),orpersplit(local)endfork=1tomdoGkv←=∑j∈{j|sk,v≥xjk>sk,v−1}gjHkv←=∑j∈{j|sk,v≥xjk>sk,v−1}hjendcallSplit_Finding()fork=1tomdoProposeSk={sk1,sk2,...,skl}bypercentileonfeaturekProposecanbedonepertree(global),orpersplit(local)endfork=1tomdoGkv←=∑j∈{j|sk,v≥xjk>sk,v−1}gjHkv←=∑j∈{j|sk,v≥xjk>sk,v−1}hjendcallSplit_Finding()
下面介绍找特征值域分界点{sk1,sk2,...,skl}{sk1,sk2,...,skl}的方法,加权分位数略图。为了尽可能地逼近最佳分裂点,我们需要保证采样后数据分布同原始数据尽可能一致。记Dk={(x1k,h1),(x2k,h2)⋯(xnk,hn)}Dk={(x1k,h1),(x2k,h2)⋯(xnk,hn)}表示 每个训练样本的第kk维特征值和对应二阶导数。接下来定义排序函数为rk(⋅):R→[0,+∞)rk(⋅):R→[0,+∞)
rk(z)=1∑(x,h)∈Dkh∑(x,h)∈Dk,x<zhrk(z)=1∑(x,h)∈Dkh∑(x,h)∈Dk,x<zh
函数表示特征的值小于zz的样本分布占比,其中二阶导数hh可以视为权重,后面论述。在这个排序函数下,我们找到一组点{sk1,sk2,...,skl}{sk1,sk2,...,skl} ,满足:
|rk(sk,j)−rk(sk,j+1)|<ε∣∣rk(sk,j)−rk(sk,j+1)∣∣<ε
其中,sk1=minixik,skl=maxixiksk1=minixik,skl=maxixik。εε为采样率,直观上理解,我们最后会得到1/ε1/ε个分界点。
讨论下为何hihi表示权重。从目标函数(1)(1)出发,配方得到n∑i=1[12(ft(xi)−gi/hi)2]+Ω(ft)+constantn∑i=1[12(ft(xi)−gi/hi)2]+Ω(ft)+constant
具体证明过程见原作附录。这是一个关于标签为gi/higi/hi和权重为hihi的平方误差形式。当每个权重相同的时候,退化为普通的分位数略图[9, 10]
接下来介绍另一个论文的亮点,寻找分裂点的过程中,如何克服数据稀疏。稀疏数据可能来自于missing value、大量的0值、或者特征工程例如采用one-hot表示带来的。为了解决这个问题,设定一个默认指向,当发生特征缺失的时候,将样本分类到默认分支,如下图:
默认方向由训练集中non-missing value学习而得,把不存在的值也当成missing value进行学习和处理,如下:
Sparsity_Split_Finding():Input:I,instancesetofcurrentnodeInput:d,featuredimensiongain←0G←∑i∈Igi,H←∑i∈Ihifork=1tomdoGL←0,HL←0forjinsorted(Ik,ascentorderbyxjk)doGL←GL+gj,HL←HL+hjGR←G−GL,HR←H−HLscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endGR←0,HR←0forjinsorted(Ik,ascentorderbyxjk)doGR←GR+gj,HR←HR+hjGL←G−GR,HL←H−HRscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endendOutput:split_valueanddefaultdirectionswithmaxscoreSparsity_Split_Finding():Input:I,instancesetofcurrentnodeInput:d,featuredimensiongain←0G←∑i∈Igi,H←∑i∈Ihifork=1tomdoGL←0,HL←0forjinsorted(Ik,ascentorderbyxjk)doGL←GL+gj,HL←HL+hjGR←G−GL,HR←H−HLscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endGR←0,HR←0forjinsorted(Ik,ascentorderbyxjk)doGR←GR+gj,HR←HR+hjGL←G−GR,HL←H−HRscore←max(score,G2LHL+λ+G2RHR+λ−G2H+λ)endendOutput:split_valueanddefaultdirectionswithmaxscore
论文后续内容为系统部分(包括并行化、cache-aware加速、数据块预处理,需要结合代码研读)、实验(在不同的数据集上进行了实验,包括分类、LTR、CTR预估)详见论文[11]
代码阅读:
代码易用,注释较完备,支持多种语言,目前仍在持续集成和更新,其中稀疏矩阵存储为CSR格式(见附录思考[2])。
很多核心的机器学习函数复用了作者的另一个机器学习库DMLC(Distributed Machine Learning Common Codebase),和参考[12]中介绍代码版本对比,运用较多通用的工厂类方法(见附录思考[3]),草草介绍下主体代码组成,下一篇会进行详细拆解:
|
|
附录
引文
[1] Panda, Biswanath, et al. “PLANET: massively parallel learning of tree ensembles with MapReduce.” Proceedings of the Vldb Endowment 2.2(2009):1426-1437.
[2] Tyree, Stephen, et al. “Parallel boosted regression trees for web search ranking.” International Conference on World Wide Web ACM, 2011:387-396.
[3] Ye, Jerry, et al. “Stochastic gradient boosted distributed decision trees.” ACM Conference on Information and Knowledge Management, CIKM 2009, Hong Kong, China, November 2009:2061-2064.
[4] Johnson, Rie, and Z. Tong. “Learning Nonlinear Functions Using Regularized Greedy Forest.” IEEE Transactions on Pattern Analysis & Machine Intelligence 36.5(2013):942-54.
[5] Friedman, Jerome H., and B. E. Popescu. “Importance Sampled Learning Ensembles.” (2003).
[6] Li, Ping, C. J. C. Burges, and Q. Wu. “McRank: Learning to Rank Using Multiple Classification and Gradient Boosting.” Advances in Neural Information Processing Systems (2007):897-904.
[7] Bekkerman, Ron, M. Bilenko, and J. Langford. “Scaling up machine learning: parallel and distributed approaches.” ACM SIGKDD International Conference Tutorials ACM, 2011:1-1.
[8] Tyree, Stephen, et al. “Parallel boosted regression trees for web search ranking.” International Conference on World Wide Web ACM, 2011:387-396.
[9] Greenwald, Michael. “Space-efficient online computation of quantile summaries.” Acm Sigmod Record 30.2(2001):58–66.
[10] Zhang, Qi, and W. Wang. “A Fast Algorithm for Approximate Quantiles in High Speed Data Streams.” (2007):29-29.
[11]Chen, Tianqi, and C. Guestrin. “XGBoost: A Scalable Tree Boosting System.” (2016).
[12] xgboost导读和实战(百度云地址),王超,陈帅华
[13] xgboost代码参数说明, @zc02051126译
[14]xgboost: 速度快效果好的boosting模型,何通
[15]Introduction to Boosted Trees, 官网
思考
[1].不能直接看明白,手痒推导一把就通顺多了。函数f(x)f(x)在x0x0处的二阶展开式为:f(x)=f(x0)+f′(x0)(x−x0)+f″(x0)(x−x0)2
我们对l(yi,x)在ˆy(t−1)i处进行二阶展开得到
l(yi,x)=l(yi,ˆy(t−1))+∂l(yi,ˆy(t−1))∂ˆy(t−1)(x−ˆy(t−1))+∂2l(yi,ˆy(t−1))∂ˆy(t−1)(x−ˆy(t−1))2
令x=ˆy(t−1)+ft(xi),且记一阶导为gi=∂l(yi,ˆy(t−1))∂ˆy(t−1),二阶导为hi=∂2l(yi,ˆy(t−1))∂ˆy(t−1)。
我们得到l(yi,ˆy(t−1)+ft(xi))的二阶泰勒展开:
l(yi,ˆy(t−1))+gift(xi)+12hif2t(xi)
带入目标函数可得:
L(t)≃n∑i=1[l(yi,ˆy(t−1))+gift(xi)+12hif2t(xi)]+Ω(ft)
[2]:csr格式是什么?类似的格式有什么,它们之间的区别是什么?
参见:稀疏矩阵存储格式总结+存储效率对比:COO,CSR,DIA,ELL,HYB
CSR有行偏移,列下标,值三种元素。行偏移数组大小为(总行数目+1),行偏移表示某一行的第一个元素在values里面的起始偏移位置。数值和列号与COO一致,表示一个元素以及其列号。如上图中,第一行元素1是0偏移,第二行元素2是2偏移,第三行元素5是4偏移,第4行元素6是7偏移。在行偏移的最后补上矩阵总的元素个数,本例中是9。
[3]思考:通用的工厂类,@tenfyzhong
ShawnXiao@baidu