xgboost 实战以及源代码分析

1.序

　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

2.xgboost vs gbdt

　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

 

 

图1

 

　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

 

 

　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

3.原理

对于上面给出的目标函数，我们可以进一步化简

（1）定义树的复杂度

对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

 

 

定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

 

 

注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

 

 

这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

 

 

最终公式可以化简为

 

 

通过对求导等于0，可以得到

 

 

然后把最优解代入得到：

 

 

（2）打分函数计算示例

Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

（3）分裂节点

论文中给出了两种分裂节点的方法

（1）贪心法：

每一次尝试去对已有的叶子加入一个分割

 

 

对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

 

 

我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

下面是论文中的算法

 

 

（2）近似算法：

主要针对数据太大，不能直接进行计算

 

 

4.自定义损失函数（指定grad、hess）

（1）损失函数

（2）grad、hess推导

（3）官方代码

1. #!/usr/bin/python

2. import numpy as np

3. import xgboost as xgb

4. ###

5. # advanced: customized loss function

6. #

7. print ('start running example to used customized objective function')

8.  

9. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

10. dtest = xgb.DMatrix('../data/agaricus.txt.test')

11.  

12. # note: for customized objective function, we leave objective as default

13. # note: what we are getting is margin value in prediction

14. # you must know what you are doing

15. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

16. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

17. num_round = 2

18.  

19. # user define objective function, given prediction, return gradient and second order gradient

20. # this is log likelihood loss

21. def logregobj(preds, dtrain):

22. labels = dtrain.get_label()

23. preds = 1.0 / (1.0 + np.exp(-preds))

24. grad = preds - labels

25. hess = preds * (1.0-preds)

26. return grad, hess

27.  

28. # user defined evaluation function, return a pair metric_name, result

29. # NOTE: when you do customized loss function, the default prediction value is margin

30. # this may make builtin evaluation metric not function properly

31. # for example, we are doing logistic loss, the prediction is score before logistic transformation

32. # the builtin evaluation error assumes input is after logistic transformation

33. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

34. def evalerror(preds, dtrain):

35. labels = dtrain.get_label()

36. # return a pair metric_name, result

37. # since preds are margin(before logistic transformation, cutoff at 0)

38. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

39.  

40. # training with customized objective, we can also do step by step training

41. # simply look at xgboost.py's implementation of train

42. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42

5.Xgboost调参

由于xgboost的参数过多，这里介绍三种思路

（1）GridSearch

（2）Hyperopt

（3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

6.工程实现优化

（1）Column Blocks and Parallelization

 

 

（2）Cache Aware Access

• A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
• The main thread accumulates gradients statistics in the continuous buff er.

（3）System Tricks

• Block pre-fetching.
• Utilize multiple disks to parallelize disk operations.
• LZ4 compression(popular recent years for outstanding performance).
• Unrolling loops.
• OpenMP

7.代码走读

这块非常感谢杨军老师的无私奉献【4】

个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

• 1

• 1

实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

1. cli_main.cc:

2. main()

3. -> CLIRunTask()

4. -> CLITrain()

5. -> DMatrix::Load()

6. -> learner = Learner::Create()

7. -> learner->Configure()

8. -> learner->InitModel()

9. -> for (i = 0; i < param.num_round; ++i)

10. -> learner->UpdateOneIter()

11. -> learner->Save()

12. learner.cc:

13. Create()

14. -> new LearnerImpl()

15. Configure()

16. InitModel()

17. -> LazyInitModel()

18. -> obj_ = ObjFunction::Create()

19. -> objective.cc

20. Create()

21. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

22. LambdaRankObj(rank_obj.cc)/

23. RegLossObj(regression_obj.cc)/

24. PoissonRegression(regression_obj.cc)

25. -> gbm_ = GradientBooster::Create()

26. -> gbm.cc

27. Create()

28. -> GBTree(gbtree.cc)/

29. GBLinear(gblinear.cc)

30. -> obj_->Configure()

31. -> gbm_->Configure()

32. UpdateOneIter()

33. -> PredictRaw()

34. -> obj_->GetGradient()

35. -> gbm_->DoBoost()

36.  

37. gbtree.cc:

38. Configure()

39. -> for (up in updaters)

40. -> up->Init()

41. DoBoost()

42. -> BoostNewTrees()

43. -> new_tree = new RegTree()

44. -> for (up in updaters)

45. -> up->Update(new_tree)

46.  

47. tree_updater.cc:

48. Create()

49. -> ColMaker/DistColMaker(updater_colmaker.cc)/

50. SketchMaker(updater_skmaker.cc)/

51. TreeRefresher(updater_refresh.cc)/

52. TreePruner(updater_prune.cc)/

53. HistMaker/CQHistMaker/

54. GlobalProposalHistMaker/

55. QuantileHistMaker(updater_histmaker.cc)/

56. TreeSyncher(updater_sync.cc)

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56

从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
　　I). inter-feature exact parallelism （特征级精确并行） 
　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
　　III). intra-feature parallelism （特征内并行）

此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

1. updater_colmaker.cc:

2. ColMaker::Update()

3. -> Builder builder;

4. -> builder.Update()

5. -> InitData()

6. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

7. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

8. -> for (depth = 0; depth < 树的最大深度; ++depth)

9. -> FindSplit()

10. -> for (each feature) // 通过OpenMP获取

11. // inter-feature parallelism

12. -> UpdateSolution()

13. -> EnumerateSplit() // 每个执行线程处理一个特征，

14. // 选出每个特征的

15. // 最优split point

16. -> ParallelFindSplit()

17. // 多个执行线程同时处理一个特征，选出该特征

18. //的最优split point;

19. // 在每个线程里汇总各个线程内分配到的数据样

20. //本的统计量(grad/hess);

21. // aggregate所有线程的样本统计(grad/hess)，

22. //计算出每个线程分配到的样本集合的边界特征值作为

23. //split point的最优分割点;

24. // 在每个线程分配到的样本集合对应的特征值集合进

25. //行枚举作为split point，选出最优分割点

26. -> SyncBestSolution()

27. // 上面的UpdateSolution()/ParallelFindSplit()

28. //会为所有待扩展分割的叶结点找到特征维度的最优split

29. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

30. //的最优split point，OpenMP线程2会找到特征F2的最

31. //优split point，所以需要进行全局sync，找到叶结点A

32. //的最优split point。

33. -> 为需要进行分割的叶结点创建孩子结点

34. -> ResetPosition()

35. //根据上一步的分割动作，更新样本到树结点的映射关系

36. // Missing Value(i.e. default)和非Missing Value(i.e.

37. //non-default)分别处理

38. -> UpdateQueueExpand()

39. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

40. //起始基础

41. -> InitNewNode() // 为可用于split的树结点计算统计量

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41

8.python、R对于xgboost的简单使用

任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

【Python】

 

 

【R】

 

 

9.xgboost中比较重要的参数介绍

（1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

• “reg:linear” –线性回归。
• “reg:logistic” –逻辑回归。
• “binary:logistic” –二分类的逻辑回归问题，输出为概率。
• “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
• “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
• “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
• “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
• “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

（2）’eval_metric’ The choices are listed below，评估指标:

• “rmse”: root mean square error
• “logloss”: negative log-likelihood
• “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
• “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
• “mlogloss”: Multiclass logloss
• “auc”: Area under the curve for ranking evaluation.
• “ndcg”:Normalized Discounted Cumulative Gain
• “map”:Mean average precision
• “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
• “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

（3）lambda [default=0] L2 正则的惩罚系数

（4）alpha [default=0] L1 正则的惩罚系数

（5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

（6）eta [default=0.3] 
为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
取值范围为：[0,1]

（7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

（8）min_child_weight [default=1] 
孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
取值范围为: [0,∞]

10.DART

核心思想就是将dropout引入XGBoost

示例代码

1. import xgboost as xgb

2. # read in data

3. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

4. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

5. # specify parameters via map

6. param = {'booster': 'dart',

7. 'max_depth': 5, 'learning_rate': 0.1,

8. 'objective': 'binary:logistic', 'silent': True,

9. 'sample_type': 'uniform',

10. 'normalize_type': 'tree',

11. 'rate_drop': 0.1,

12. 'skip_drop': 0.5}

13. num_round = 50

14. bst = xgb.train(param, dtrain, num_round)

15. # make prediction

16. # ntree_limit must not be 0

17. preds = bst.predict(dtest, ntree_limit=num_round)

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17

更多细节可以阅读参考文献5

11.csr_matrix训练XGBoost

当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

1. import pandas as pd

2. import numpy as np

3. import seaborn as sns

4. import matplotlib.pyplot as plt

5. import os

6. from sklearn.preprocessing import LabelEncoder

7. from scipy.sparse import csr_matrix, hstack

8. import xgboost as xgb

9. from sklearn.cross_validation import StratifiedKFold

10. from sklearn.metrics import log_loss

11.  

12. datadir = '../input'

13. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

14. index_col='device_id')

15. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

16. index_col = 'device_id')

17. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

18. # Get rid of duplicate device ids in phone

19. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

20. events = pd.read_csv(os.path.join(datadir,'events.csv'),

21. parse_dates=['timestamp'], index_col='event_id')

22. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

23. usecols=['event_id','app_id','is_active'],

24. dtype={'is_active':bool})

25. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

26.  

27. gatrain['trainrow'] = np.arange(gatrain.shape[0])

28. gatest['testrow'] = np.arange(gatest.shape[0])

29.  

30. brandencoder = LabelEncoder().fit(phone.phone_brand)

31. phone['brand'] = brandencoder.transform(phone['phone_brand'])

32. gatrain['brand'] = phone['brand']

33. gatest['brand'] = phone['brand']

34. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

35. (gatrain.trainrow, gatrain.brand)))

36. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

37. (gatest.testrow, gatest.brand)))

38. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

39.  

40. m = phone.phone_brand.str.cat(phone.device_model)

41. modelencoder = LabelEncoder().fit(m)

42. phone['model'] = modelencoder.transform(m)

43. gatrain['model'] = phone['model']

44. gatest['model'] = phone['model']

45. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

46. (gatrain.trainrow, gatrain.model)))

47. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

48. (gatest.testrow, gatest.model)))

49. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

50.  

51. appencoder = LabelEncoder().fit(appevents.app_id)

52. appevents['app'] = appencoder.transform(appevents.app_id)

53. napps = len(appencoder.classes_)

54. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

55. .groupby(['device_id','app'])['app'].agg(['size'])

56. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

57. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

58. .reset_index())

59.  

60. d = deviceapps.dropna(subset=['trainrow'])

61. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

62. shape=(gatrain.shape[0],napps))

63. d = deviceapps.dropna(subset=['testrow'])

64. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

65. shape=(gatest.shape[0],napps))

66. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

67.  

68. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

69. applabels['app'] = appencoder.transform(applabels.app_id)

70. labelencoder = LabelEncoder().fit(applabels.label_id)

71. applabels['label'] = labelencoder.transform(applabels.label_id)

72. nlabels = len(labelencoder.classes_)

73.  

74. devicelabels = (deviceapps[['device_id','app']]

75. .merge(applabels[['app','label']])

76. .groupby(['device_id','label'])['app'].agg(['size'])

77. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

78. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

79. .reset_index())

80. devicelabels.head()

81.  

82. d = devicelabels.dropna(subset=['trainrow'])

83. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

84. shape=(gatrain.shape[0],nlabels))

85. d = devicelabels.dropna(subset=['testrow'])

86. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

87. shape=(gatest.shape[0],nlabels))

88. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

89.  

90. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

91. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

92. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

93.  

94. targetencoder = LabelEncoder().fit(gatrain.group)

95. y = targetencoder.transform(gatrain.group)

96.  

97. ########## XGBOOST ##########

98.  

99. params = {}

100. params['booster'] = 'gblinear'

101. params['objective'] = "multi:softprob"

102. params['eval_metric'] = 'mlogloss'

103. params['eta'] = 0.005

104. params['num_class'] = 12

105. params['lambda'] = 3

106. params['alpha'] = 2

107.  

108. # Random 10% for validation

109. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

110.  

111. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

112. ytr, yte = y[kf[0]], y[kf[1]]

113.  

114. print('Training set: ' + str(Xtr.shape))

115. print('Validation set: ' + str(Xte.shape))

116.  

117. d_train = xgb.DMatrix(Xtr, label=ytr)

118. d_valid = xgb.DMatrix(Xte, label=yte)

119.  

120. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

121.  

122. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

123.  

124. pred = clf.predict(xgb.DMatrix(Xtest))

125.  

126. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

127. pred.head()

128. pred.to_csv('sparse_xgb.csv', index=True)

129.  

130. #params['lambda'] = 1

131. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

132. # params['alpha'] = alpha

133. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

134. 1.序

     　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

     　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

     2.xgboost vs gbdt

     　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

      

      

     图1

      

     　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

      

     

      

     　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

     3.原理

     对于上面给出的目标函数，我们可以进一步化简

     （1）定义树的复杂度

     对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

      

     

      

     定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

      

     

      

     注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

     在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

      

     

      

     这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

      

     

      

     最终公式可以化简为

      

     

      

     通过对求导等于0，可以得到

      

     

      

     然后把最优解代入得到：

      

     

      

     （2）打分函数计算示例

     Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

     

     （3）分裂节点

     论文中给出了两种分裂节点的方法

     （1）贪心法：

     每一次尝试去对已有的叶子加入一个分割

      

     

      

     对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

      

     

      

     我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

     观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

     下面是论文中的算法

      

     

      

     （2）近似算法：

     主要针对数据太大，不能直接进行计算

      

     

      

     4.自定义损失函数（指定grad、hess）

     （1）损失函数

     

     （2）grad、hess推导

     

     （3）官方代码

135. #!/usr/bin/python

136. import numpy as np

137. import xgboost as xgb

138. ###

139. # advanced: customized loss function

140. #

141. print ('start running example to used customized objective function')

142.  

143. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

144. dtest = xgb.DMatrix('../data/agaricus.txt.test')

145.  

146. # note: for customized objective function, we leave objective as default

147. # note: what we are getting is margin value in prediction

148. # you must know what you are doing

149. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

150. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

151. num_round = 2

152.  

153. # user define objective function, given prediction, return gradient and second order gradient

154. # this is log likelihood loss

155. def logregobj(preds, dtrain):

156. labels = dtrain.get_label()

157. preds = 1.0 / (1.0 + np.exp(-preds))

158. grad = preds - labels

159. hess = preds * (1.0-preds)

160. return grad, hess

161.  

162. # user defined evaluation function, return a pair metric_name, result

163. # NOTE: when you do customized loss function, the default prediction value is margin

164. # this may make builtin evaluation metric not function properly

165. # for example, we are doing logistic loss, the prediction is score before logistic transformation

166. # the builtin evaluation error assumes input is after logistic transformation

167. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

168. def evalerror(preds, dtrain):

169. labels = dtrain.get_label()

170. # return a pair metric_name, result

171. # since preds are margin(before logistic transformation, cutoff at 0)

172. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

173.  

174. # training with customized objective, we can also do step by step training

175. # simply look at xgboost.py's implementation of train

176. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

177. 1
178. 2
179. 3
180. 4
181. 5
182. 6
183. 7
184. 8
185. 9
186. 10
187. 11
188. 12
189. 13
190. 14
191. 15
192. 16
193. 17
194. 18
195. 19
196. 20
197. 21
198. 22
199. 23
200. 24
201. 25
202. 26
203. 27
204. 28
205. 29
206. 30
207. 31
208. 32
209. 33
210. 34
211. 35
212. 36
213. 37
214. 38
215. 39
216. 40
217. 41
218. 42
219. 1
220. 2
221. 3
222. 4
223. 5
224. 6
225. 7
226. 8
227. 9
228. 10
229. 11
230. 12
231. 13
232. 14
233. 15
234. 16
235. 17
236. 18
237. 19
238. 20
239. 21
240. 22
241. 23
242. 24
243. 25
244. 26
245. 27
246. 28
247. 29
248. 30
249. 31
250. 32
251. 33
252. 34
253. 35
254. 36
255. 37
256. 38
257. 39
258. 40
259. 41
260. 42

261. 5.Xgboost调参

     由于xgboost的参数过多，这里介绍三种思路

     （1）GridSearch

     （2）Hyperopt

     （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

     6.工程实现优化

     （1）Column Blocks and Parallelization

      

     

      

     （2）Cache Aware Access

262. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
263. The main thread accumulates gradients statistics in the continuous buff er.

264. （3）System Tricks

265. Block pre-fetching.
266. Utilize multiple disks to parallelize disk operations.
267. LZ4 compression(popular recent years for outstanding performance).
268. Unrolling loops.
269. OpenMP

270. 7.代码走读

     这块非常感谢杨军老师的无私奉献【4】

     个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

     find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

271. 1
272. 1

273. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

274. cli_main.cc:

275. main()

276. -> CLIRunTask()

277. -> CLITrain()

278. -> DMatrix::Load()

279. -> learner = Learner::Create()

280. -> learner->Configure()

281. -> learner->InitModel()

282. -> for (i = 0; i < param.num_round; ++i)

283. -> learner->UpdateOneIter()

284. -> learner->Save()

285. learner.cc:

286. Create()

287. -> new LearnerImpl()

288. Configure()

289. InitModel()

290. -> LazyInitModel()

291. -> obj_ = ObjFunction::Create()

292. -> objective.cc

293. Create()

294. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

295. LambdaRankObj(rank_obj.cc)/

296. RegLossObj(regression_obj.cc)/

297. PoissonRegression(regression_obj.cc)

298. -> gbm_ = GradientBooster::Create()

299. -> gbm.cc

300. Create()

301. -> GBTree(gbtree.cc)/

302. GBLinear(gblinear.cc)

303. -> obj_->Configure()

304. -> gbm_->Configure()

305. UpdateOneIter()

306. -> PredictRaw()

307. -> obj_->GetGradient()

308. -> gbm_->DoBoost()

309.  

310. gbtree.cc:

311. Configure()

312. -> for (up in updaters)

313. -> up->Init()

314. DoBoost()

315. -> BoostNewTrees()

316. -> new_tree = new RegTree()

317. -> for (up in updaters)

318. -> up->Update(new_tree)

319.  

320. tree_updater.cc:

321. Create()

322. -> ColMaker/DistColMaker(updater_colmaker.cc)/

323. SketchMaker(updater_skmaker.cc)/

324. TreeRefresher(updater_refresh.cc)/

325. TreePruner(updater_prune.cc)/

326. HistMaker/CQHistMaker/

327. GlobalProposalHistMaker/

328. QuantileHistMaker(updater_histmaker.cc)/

329. TreeSyncher(updater_sync.cc)

330. 1
331. 2
332. 3
333. 4
334. 5
335. 6
336. 7
337. 8
338. 9
339. 10
340. 11
341. 12
342. 13
343. 14
344. 15
345. 16
346. 17
347. 18
348. 19
349. 20
350. 21
351. 22
352. 23
353. 24
354. 25
355. 26
356. 27
357. 28
358. 29
359. 30
360. 31
361. 32
362. 33
363. 34
364. 35
365. 36
366. 37
367. 38
368. 39
369. 40
370. 41
371. 42
372. 43
373. 44
374. 45
375. 46
376. 47
377. 48
378. 49
379. 50
380. 51
381. 52
382. 53
383. 54
384. 55
385. 56
386. 1
387. 2
388. 3
389. 4
390. 5
391. 6
392. 7
393. 8
394. 9
395. 10
396. 11
397. 12
398. 13
399. 14
400. 15
401. 16
402. 17
403. 18
404. 19
405. 20
406. 21
407. 22
408. 23
409. 24
410. 25
411. 26
412. 27
413. 28
414. 29
415. 30
416. 31
417. 32
418. 33
419. 34
420. 35
421. 36
422. 37
423. 38
424. 39
425. 40
426. 41
427. 42
428. 43
429. 44
430. 45
431. 46
432. 47
433. 48
434. 49
435. 50
436. 51
437. 52
438. 53
439. 54
440. 55
441. 56

442. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

     I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
     II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
     III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
     　　I). inter-feature exact parallelism （特征级精确并行） 
     　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
     　　III). intra-feature parallelism （特征内并行）

     此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

     XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

     以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

443. updater_colmaker.cc:

444. ColMaker::Update()

445. -> Builder builder;

446. -> builder.Update()

447. -> InitData()

448. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

449. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

450. -> for (depth = 0; depth < 树的最大深度; ++depth)

451. -> FindSplit()

452. -> for (each feature) // 通过OpenMP获取

453. // inter-feature parallelism

454. -> UpdateSolution()

455. -> EnumerateSplit() // 每个执行线程处理一个特征，

456. // 选出每个特征的

457. // 最优split point

458. -> ParallelFindSplit()

459. // 多个执行线程同时处理一个特征，选出该特征

460. //的最优split point;

461. // 在每个线程里汇总各个线程内分配到的数据样

462. //本的统计量(grad/hess);

463. // aggregate所有线程的样本统计(grad/hess)，

464. //计算出每个线程分配到的样本集合的边界特征值作为

465. //split point的最优分割点;

466. // 在每个线程分配到的样本集合对应的特征值集合进

467. //行枚举作为split point，选出最优分割点

468. -> SyncBestSolution()

469. // 上面的UpdateSolution()/ParallelFindSplit()

470. //会为所有待扩展分割的叶结点找到特征维度的最优split

471. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

472. //的最优split point，OpenMP线程2会找到特征F2的最

473. //优split point，所以需要进行全局sync，找到叶结点A

474. //的最优split point。

475. -> 为需要进行分割的叶结点创建孩子结点

476. -> ResetPosition()

477. //根据上一步的分割动作，更新样本到树结点的映射关系

478. // Missing Value(i.e. default)和非Missing Value(i.e.

479. //non-default)分别处理

480. -> UpdateQueueExpand()

481. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

482. //起始基础

483. -> InitNewNode() // 为可用于split的树结点计算统计量

484. 1
485. 2
486. 3
487. 4
488. 5
489. 6
490. 7
491. 8
492. 9
493. 10
494. 11
495. 12
496. 13
497. 14
498. 15
499. 16
500. 17
501. 18
502. 19
503. 20
504. 21
505. 22
506. 23
507. 24
508. 25
509. 26
510. 27
511. 28
512. 29
513. 30
514. 31
515. 32
516. 33
517. 34
518. 35
519. 36
520. 37
521. 38
522. 39
523. 40
524. 41
525. 1
526. 2
527. 3
528. 4
529. 5
530. 6
531. 7
532. 8
533. 9
534. 10
535. 11
536. 12
537. 13
538. 14
539. 15
540. 16
541. 17
542. 18
543. 19
544. 20
545. 21
546. 22
547. 23
548. 24
549. 25
550. 26
551. 27
552. 28
553. 29
554. 30
555. 31
556. 32
557. 33
558. 34
559. 35
560. 36
561. 37
562. 38
563. 39
564. 40
565. 41

566. 8.python、R对于xgboost的简单使用

     任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

     【Python】

      

     

      

     【R】

      

     

      

     9.xgboost中比较重要的参数介绍

     （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

567. “reg:linear” –线性回归。
568. “reg:logistic” –逻辑回归。
569. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
570. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
571. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
572. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
573. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
574. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

575. （2）’eval_metric’ The choices are listed below，评估指标:

576. “rmse”: root mean square error
577. “logloss”: negative log-likelihood
578. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
579. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
580. “mlogloss”: Multiclass logloss
581. “auc”: Area under the curve for ranking evaluation.
582. “ndcg”:Normalized Discounted Cumulative Gain
583. “map”:Mean average precision
584. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
585. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

586. （3）lambda [default=0] L2 正则的惩罚系数

     （4）alpha [default=0] L1 正则的惩罚系数

     （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

     （6）eta [default=0.3] 
     为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
     取值范围为：[0,1]

     （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

     （8）min_child_weight [default=1] 
     孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
     取值范围为: [0,∞]

     10.DART

     核心思想就是将dropout引入XGBoost

     示例代码

587. import xgboost as xgb

588. # read in data

589. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

590. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

591. # specify parameters via map

592. param = {'booster': 'dart',

593. 'max_depth': 5, 'learning_rate': 0.1,

594. 'objective': 'binary:logistic', 'silent': True,

595. 'sample_type': 'uniform',

596. 'normalize_type': 'tree',

597. 'rate_drop': 0.1,

598. 'skip_drop': 0.5}

599. num_round = 50

600. bst = xgb.train(param, dtrain, num_round)

601. # make prediction

602. # ntree_limit must not be 0

603. preds = bst.predict(dtest, ntree_limit=num_round)

604. 1
605. 2
606. 3
607. 4
608. 5
609. 6
610. 7
611. 8
612. 9
613. 10
614. 11
615. 12
616. 13
617. 14
618. 15
619. 16
620. 17
621. 1
622. 2
623. 3
624. 4
625. 5
626. 6
627. 7
628. 8
629. 9
630. 10
631. 11
632. 12
633. 13
634. 14
635. 15
636. 16
637. 17

638. 更多细节可以阅读参考文献5

     11.csr_matrix训练XGBoost

     当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

     下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

639. import pandas as pd

640. import numpy as np

641. import seaborn as sns

642. import matplotlib.pyplot as plt

643. import os

644. from sklearn.preprocessing import LabelEncoder

645. from scipy.sparse import csr_matrix, hstack

646. import xgboost as xgb

647. from sklearn.cross_validation import StratifiedKFold

648. from sklearn.metrics import log_loss

649.  

650. datadir = '../input'

651. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

652. index_col='device_id')

653. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

654. index_col = 'device_id')

655. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

656. # Get rid of duplicate device ids in phone

657. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

658. events = pd.read_csv(os.path.join(datadir,'events.csv'),

659. parse_dates=['timestamp'], index_col='event_id')

660. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

661. usecols=['event_id','app_id','is_active'],

662. dtype={'is_active':bool})

663. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

664.  

665. gatrain['trainrow'] = np.arange(gatrain.shape[0])

666. gatest['testrow'] = np.arange(gatest.shape[0])

667.  

668. brandencoder = LabelEncoder().fit(phone.phone_brand)

669. phone['brand'] = brandencoder.transform(phone['phone_brand'])

670. gatrain['brand'] = phone['brand']

671. gatest['brand'] = phone['brand']

672. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

673. (gatrain.trainrow, gatrain.brand)))

674. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

675. (gatest.testrow, gatest.brand)))

676. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

677.  

678. m = phone.phone_brand.str.cat(phone.device_model)

679. modelencoder = LabelEncoder().fit(m)

680. phone['model'] = modelencoder.transform(m)

681. gatrain['model'] = phone['model']

682. gatest['model'] = phone['model']

683. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

684. (gatrain.trainrow, gatrain.model)))

685. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

686. (gatest.testrow, gatest.model)))

687. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

688.  

689. appencoder = LabelEncoder().fit(appevents.app_id)

690. appevents['app'] = appencoder.transform(appevents.app_id)

691. napps = len(appencoder.classes_)

692. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

693. .groupby(['device_id','app'])['app'].agg(['size'])

694. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

695. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

696. .reset_index())

697.  

698. d = deviceapps.dropna(subset=['trainrow'])

699. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

700. shape=(gatrain.shape[0],napps))

701. d = deviceapps.dropna(subset=['testrow'])

702. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

703. shape=(gatest.shape[0],napps))

704. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

705.  

706. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

707. applabels['app'] = appencoder.transform(applabels.app_id)

708. labelencoder = LabelEncoder().fit(applabels.label_id)

709. applabels['label'] = labelencoder.transform(applabels.label_id)

710. nlabels = len(labelencoder.classes_)

711.  

712. devicelabels = (deviceapps[['device_id','app']]

713. .merge(applabels[['app','label']])

714. .groupby(['device_id','label'])['app'].agg(['size'])

715. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

716. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

717. .reset_index())

718. devicelabels.head()

719.  

720. d = devicelabels.dropna(subset=['trainrow'])

721. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

722. shape=(gatrain.shape[0],nlabels))

723. d = devicelabels.dropna(subset=['testrow'])

724. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

725. shape=(gatest.shape[0],nlabels))

726. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

727.  

728. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

729. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

730. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

731.  

732. targetencoder = LabelEncoder().fit(gatrain.group)

733. y = targetencoder.transform(gatrain.group)

734.  

735. ########## XGBOOST ##########

736.  

737. params = {}

738. params['booster'] = 'gblinear'

739. params['objective'] = "multi:softprob"

740. params['eval_metric'] = 'mlogloss'

741. params['eta'] = 0.005

742. params['num_class'] = 12

743. params['lambda'] = 3

744. params['alpha'] = 2

745.  

746. # Random 10% for validation

747. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

748.  

749. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

750. ytr, yte = y[kf[0]], y[kf[1]]

751.  

752. print('Training set: ' + str(Xtr.shape))

753. print('Validation set: ' + str(Xte.shape))

754.  

755. d_train = xgb.DMatrix(Xtr, label=ytr)

756. d_valid = xgb.DMatrix(Xte, label=yte)

757.  

758. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

759.  

760. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

761.  

762. pred = clf.predict(xgb.DMatrix(Xtest))

763.  

764. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

765. pred.head()

766. pred.to_csv('sparse_xgb.csv', index=True)

767.  

768. #params['lambda'] = 1

769. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

770. # params['alpha'] = alpha

771. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

772. 1.序

     　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

     　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

     2.xgboost vs gbdt

     　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

      

      

     图1

      

     　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

      

     

      

     　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

     3.原理

     对于上面给出的目标函数，我们可以进一步化简

     （1）定义树的复杂度

     对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

      

     

      

     定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

      

     

      

     注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

     在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

      

     

      

     这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

      

     

      

     最终公式可以化简为

      

     

      

     通过对求导等于0，可以得到

      

     

      

     然后把最优解代入得到：

      

     

      

     （2）打分函数计算示例

     Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

     

     （3）分裂节点

     论文中给出了两种分裂节点的方法

     （1）贪心法：

     每一次尝试去对已有的叶子加入一个分割

      

     

      

     对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

      

     

      

     我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

     观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

     下面是论文中的算法

      

     

      

     （2）近似算法：

     主要针对数据太大，不能直接进行计算

      

     

      

     4.自定义损失函数（指定grad、hess）

     （1）损失函数

     

     （2）grad、hess推导

     

     （3）官方代码

773. #!/usr/bin/python

774. import numpy as np

775. import xgboost as xgb

776. ###

777. # advanced: customized loss function

778. #

779. print ('start running example to used customized objective function')

780.  

781. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

782. dtest = xgb.DMatrix('../data/agaricus.txt.test')

783.  

784. # note: for customized objective function, we leave objective as default

785. # note: what we are getting is margin value in prediction

786. # you must know what you are doing

787. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

788. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

789. num_round = 2

790.  

791. # user define objective function, given prediction, return gradient and second order gradient

792. # this is log likelihood loss

793. def logregobj(preds, dtrain):

794. labels = dtrain.get_label()

795. preds = 1.0 / (1.0 + np.exp(-preds))

796. grad = preds - labels

797. hess = preds * (1.0-preds)

798. return grad, hess

799.  

800. # user defined evaluation function, return a pair metric_name, result

801. # NOTE: when you do customized loss function, the default prediction value is margin

802. # this may make builtin evaluation metric not function properly

803. # for example, we are doing logistic loss, the prediction is score before logistic transformation

804. # the builtin evaluation error assumes input is after logistic transformation

805. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

806. def evalerror(preds, dtrain):

807. labels = dtrain.get_label()

808. # return a pair metric_name, result

809. # since preds are margin(before logistic transformation, cutoff at 0)

810. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

811.  

812. # training with customized objective, we can also do step by step training

813. # simply look at xgboost.py's implementation of train

814. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

815. 1
816. 2
817. 3
818. 4
819. 5
820. 6
821. 7
822. 8
823. 9
824. 10
825. 11
826. 12
827. 13
828. 14
829. 15
830. 16
831. 17
832. 18
833. 19
834. 20
835. 21
836. 22
837. 23
838. 24
839. 25
840. 26
841. 27
842. 28
843. 29
844. 30
845. 31
846. 32
847. 33
848. 34
849. 35
850. 36
851. 37
852. 38
853. 39
854. 40
855. 41
856. 42
857. 1
858. 2
859. 3
860. 4
861. 5
862. 6
863. 7
864. 8
865. 9
866. 10
867. 11
868. 12
869. 13
870. 14
871. 15
872. 16
873. 17
874. 18
875. 19
876. 20
877. 21
878. 22
879. 23
880. 24
881. 25
882. 26
883. 27
884. 28
885. 29
886. 30
887. 31
888. 32
889. 33
890. 34
891. 35
892. 36
893. 37
894. 38
895. 39
896. 40
897. 41
898. 42

899. 5.Xgboost调参

     由于xgboost的参数过多，这里介绍三种思路

     （1）GridSearch

     （2）Hyperopt

     （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

     6.工程实现优化

     （1）Column Blocks and Parallelization

      

     

      

     （2）Cache Aware Access

900. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
901. The main thread accumulates gradients statistics in the continuous buff er.

902. （3）System Tricks

903. Block pre-fetching.
904. Utilize multiple disks to parallelize disk operations.
905. LZ4 compression(popular recent years for outstanding performance).
906. Unrolling loops.
907. OpenMP

908. 7.代码走读

     这块非常感谢杨军老师的无私奉献【4】

     个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

     find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

909. 1
910. 1

911. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

912. cli_main.cc:

913. main()

914. -> CLIRunTask()

915. -> CLITrain()

916. -> DMatrix::Load()

917. -> learner = Learner::Create()

918. -> learner->Configure()

919. -> learner->InitModel()

920. -> for (i = 0; i < param.num_round; ++i)

921. -> learner->UpdateOneIter()

922. -> learner->Save()

923. learner.cc:

924. Create()

925. -> new LearnerImpl()

926. Configure()

927. InitModel()

928. -> LazyInitModel()

929. -> obj_ = ObjFunction::Create()

930. -> objective.cc

931. Create()

932. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

933. LambdaRankObj(rank_obj.cc)/

934. RegLossObj(regression_obj.cc)/

935. PoissonRegression(regression_obj.cc)

936. -> gbm_ = GradientBooster::Create()

937. -> gbm.cc

938. Create()

939. -> GBTree(gbtree.cc)/

940. GBLinear(gblinear.cc)

941. -> obj_->Configure()

942. -> gbm_->Configure()

943. UpdateOneIter()

944. -> PredictRaw()

945. -> obj_->GetGradient()

946. -> gbm_->DoBoost()

947.  

948. gbtree.cc:

949. Configure()

950. -> for (up in updaters)

951. -> up->Init()

952. DoBoost()

953. -> BoostNewTrees()

954. -> new_tree = new RegTree()

955. -> for (up in updaters)

956. -> up->Update(new_tree)

957.  

958. tree_updater.cc:

959. Create()

960. -> ColMaker/DistColMaker(updater_colmaker.cc)/

961. SketchMaker(updater_skmaker.cc)/

962. TreeRefresher(updater_refresh.cc)/

963. TreePruner(updater_prune.cc)/

964. HistMaker/CQHistMaker/

965. GlobalProposalHistMaker/

966. QuantileHistMaker(updater_histmaker.cc)/

967. TreeSyncher(updater_sync.cc)

968. 1
969. 2
970. 3
971. 4
972. 5
973. 6
974. 7
975. 8
976. 9
977. 10
978. 11
979. 12
980. 13
981. 14
982. 15
983. 16
984. 17
985. 18
986. 19
987. 20
988. 21
989. 22
990. 23
991. 24
992. 25
993. 26
994. 27
995. 28
996. 29
997. 30
998. 31
999. 32
1000. 33
1001. 34
1002. 35
1003. 36
1004. 37
1005. 38
1006. 39
1007. 40
1008. 41
1009. 42
1010. 43
1011. 44
1012. 45
1013. 46
1014. 47
1015. 48
1016. 49
1017. 50
1018. 51
1019. 52
1020. 53
1021. 54
1022. 55
1023. 56
1024. 1
1025. 2
1026. 3
1027. 4
1028. 5
1029. 6
1030. 7
1031. 8
1032. 9
1033. 10
1034. 11
1035. 12
1036. 13
1037. 14
1038. 15
1039. 16
1040. 17
1041. 18
1042. 19
1043. 20
1044. 21
1045. 22
1046. 23
1047. 24
1048. 25
1049. 26
1050. 27
1051. 28
1052. 29
1053. 30
1054. 31
1055. 32
1056. 33
1057. 34
1058. 35
1059. 36
1060. 37
1061. 38
1062. 39
1063. 40
1064. 41
1065. 42
1066. 43
1067. 44
1068. 45
1069. 46
1070. 47
1071. 48
1072. 49
1073. 50
1074. 51
1075. 52
1076. 53
1077. 54
1078. 55
1079. 56

1080. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

      I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
      II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
      III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
      　　I). inter-feature exact parallelism （特征级精确并行） 
      　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
      　　III). intra-feature parallelism （特征内并行）

      此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

      XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

      以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

1081. updater_colmaker.cc:

1082. ColMaker::Update()

1083. -> Builder builder;

1084. -> builder.Update()

1085. -> InitData()

1086. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

1087. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

1088. -> for (depth = 0; depth < 树的最大深度; ++depth)

1089. -> FindSplit()

1090. -> for (each feature) // 通过OpenMP获取

1091. // inter-feature parallelism

1092. -> UpdateSolution()

1093. -> EnumerateSplit() // 每个执行线程处理一个特征，

1094. // 选出每个特征的

1095. // 最优split point

1096. -> ParallelFindSplit()

1097. // 多个执行线程同时处理一个特征，选出该特征

1098. //的最优split point;

1099. // 在每个线程里汇总各个线程内分配到的数据样

1100. //本的统计量(grad/hess);

1101. // aggregate所有线程的样本统计(grad/hess)，

1102. //计算出每个线程分配到的样本集合的边界特征值作为

1103. //split point的最优分割点;

1104. // 在每个线程分配到的样本集合对应的特征值集合进

1105. //行枚举作为split point，选出最优分割点

1106. -> SyncBestSolution()

1107. // 上面的UpdateSolution()/ParallelFindSplit()

1108. //会为所有待扩展分割的叶结点找到特征维度的最优split

1109. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

1110. //的最优split point，OpenMP线程2会找到特征F2的最

1111. //优split point，所以需要进行全局sync，找到叶结点A

1112. //的最优split point。

1113. -> 为需要进行分割的叶结点创建孩子结点

1114. -> ResetPosition()

1115. //根据上一步的分割动作，更新样本到树结点的映射关系

1116. // Missing Value(i.e. default)和非Missing Value(i.e.

1117. //non-default)分别处理

1118. -> UpdateQueueExpand()

1119. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

1120. //起始基础

1121. -> InitNewNode() // 为可用于split的树结点计算统计量

1122. 1
1123. 2
1124. 3
1125. 4
1126. 5
1127. 6
1128. 7
1129. 8
1130. 9
1131. 10
1132. 11
1133. 12
1134. 13
1135. 14
1136. 15
1137. 16
1138. 17
1139. 18
1140. 19
1141. 20
1142. 21
1143. 22
1144. 23
1145. 24
1146. 25
1147. 26
1148. 27
1149. 28
1150. 29
1151. 30
1152. 31
1153. 32
1154. 33
1155. 34
1156. 35
1157. 36
1158. 37
1159. 38
1160. 39
1161. 40
1162. 41
1163. 1
1164. 2
1165. 3
1166. 4
1167. 5
1168. 6
1169. 7
1170. 8
1171. 9
1172. 10
1173. 11
1174. 12
1175. 13
1176. 14
1177. 15
1178. 16
1179. 17
1180. 18
1181. 19
1182. 20
1183. 21
1184. 22
1185. 23
1186. 24
1187. 25
1188. 26
1189. 27
1190. 28
1191. 29
1192. 30
1193. 31
1194. 32
1195. 33
1196. 34
1197. 35
1198. 36
1199. 37
1200. 38
1201. 39
1202. 40
1203. 41

1204. 8.python、R对于xgboost的简单使用

      任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

      【Python】

       

      

       

      【R】

       

      

       

      9.xgboost中比较重要的参数介绍

      （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

1205. “reg:linear” –线性回归。
1206. “reg:logistic” –逻辑回归。
1207. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
1208. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
1209. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
1210. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
1211. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
1212. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

1213. （2）’eval_metric’ The choices are listed below，评估指标:

1214. “rmse”: root mean square error
1215. “logloss”: negative log-likelihood
1216. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
1217. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
1218. “mlogloss”: Multiclass logloss
1219. “auc”: Area under the curve for ranking evaluation.
1220. “ndcg”:Normalized Discounted Cumulative Gain
1221. “map”:Mean average precision
1222. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
1223. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

1224. （3）lambda [default=0] L2 正则的惩罚系数

      （4）alpha [default=0] L1 正则的惩罚系数

      （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

      （6）eta [default=0.3] 
      为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
      取值范围为：[0,1]

      （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

      （8）min_child_weight [default=1] 
      孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
      取值范围为: [0,∞]

      10.DART

      核心思想就是将dropout引入XGBoost

      示例代码

1225. import xgboost as xgb

1226. # read in data

1227. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

1228. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

1229. # specify parameters via map

1230. param = {'booster': 'dart',

1231. 'max_depth': 5, 'learning_rate': 0.1,

1232. 'objective': 'binary:logistic', 'silent': True,

1233. 'sample_type': 'uniform',

1234. 'normalize_type': 'tree',

1235. 'rate_drop': 0.1,

1236. 'skip_drop': 0.5}

1237. num_round = 50

1238. bst = xgb.train(param, dtrain, num_round)

1239. # make prediction

1240. # ntree_limit must not be 0

1241. preds = bst.predict(dtest, ntree_limit=num_round)

1242. 1
1243. 2
1244. 3
1245. 4
1246. 5
1247. 6
1248. 7
1249. 8
1250. 9
1251. 10
1252. 11
1253. 12
1254. 13
1255. 14
1256. 15
1257. 16
1258. 17
1259. 1
1260. 2
1261. 3
1262. 4
1263. 5
1264. 6
1265. 7
1266. 8
1267. 9
1268. 10
1269. 11
1270. 12
1271. 13
1272. 14
1273. 15
1274. 16
1275. 17

1276. 更多细节可以阅读参考文献5

      11.csr_matrix训练XGBoost

      当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

      下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

1277. import pandas as pd

1278. import numpy as np

1279. import seaborn as sns

1280. import matplotlib.pyplot as plt

1281. import os

1282. from sklearn.preprocessing import LabelEncoder

1283. from scipy.sparse import csr_matrix, hstack

1284. import xgboost as xgb

1285. from sklearn.cross_validation import StratifiedKFold

1286. from sklearn.metrics import log_loss

1287.  

1288. datadir = '../input'

1289. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

1290. index_col='device_id')

1291. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

1292. index_col = 'device_id')

1293. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

1294. # Get rid of duplicate device ids in phone

1295. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

1296. events = pd.read_csv(os.path.join(datadir,'events.csv'),

1297. parse_dates=['timestamp'], index_col='event_id')

1298. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

1299. usecols=['event_id','app_id','is_active'],

1300. dtype={'is_active':bool})

1301. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

1302.  

1303. gatrain['trainrow'] = np.arange(gatrain.shape[0])

1304. gatest['testrow'] = np.arange(gatest.shape[0])

1305.  

1306. brandencoder = LabelEncoder().fit(phone.phone_brand)

1307. phone['brand'] = brandencoder.transform(phone['phone_brand'])

1308. gatrain['brand'] = phone['brand']

1309. gatest['brand'] = phone['brand']

1310. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

1311. (gatrain.trainrow, gatrain.brand)))

1312. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

1313. (gatest.testrow, gatest.brand)))

1314. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

1315.  

1316. m = phone.phone_brand.str.cat(phone.device_model)

1317. modelencoder = LabelEncoder().fit(m)

1318. phone['model'] = modelencoder.transform(m)

1319. gatrain['model'] = phone['model']

1320. gatest['model'] = phone['model']

1321. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

1322. (gatrain.trainrow, gatrain.model)))

1323. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

1324. (gatest.testrow, gatest.model)))

1325. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

1326.  

1327. appencoder = LabelEncoder().fit(appevents.app_id)

1328. appevents['app'] = appencoder.transform(appevents.app_id)

1329. napps = len(appencoder.classes_)

1330. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

1331. .groupby(['device_id','app'])['app'].agg(['size'])

1332. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

1333. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

1334. .reset_index())

1335.  

1336. d = deviceapps.dropna(subset=['trainrow'])

1337. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

1338. shape=(gatrain.shape[0],napps))

1339. d = deviceapps.dropna(subset=['testrow'])

1340. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

1341. shape=(gatest.shape[0],napps))

1342. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

1343.  

1344. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

1345. applabels['app'] = appencoder.transform(applabels.app_id)

1346. labelencoder = LabelEncoder().fit(applabels.label_id)

1347. applabels['label'] = labelencoder.transform(applabels.label_id)

1348. nlabels = len(labelencoder.classes_)

1349.  

1350. devicelabels = (deviceapps[['device_id','app']]

1351. .merge(applabels[['app','label']])

1352. .groupby(['device_id','label'])['app'].agg(['size'])

1353. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

1354. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

1355. .reset_index())

1356. devicelabels.head()

1357.  

1358. d = devicelabels.dropna(subset=['trainrow'])

1359. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

1360. shape=(gatrain.shape[0],nlabels))

1361. d = devicelabels.dropna(subset=['testrow'])

1362. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

1363. shape=(gatest.shape[0],nlabels))

1364. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

1365.  

1366. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

1367. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

1368. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

1369.  

1370. targetencoder = LabelEncoder().fit(gatrain.group)

1371. y = targetencoder.transform(gatrain.group)

1372.  

1373. ########## XGBOOST ##########

1374.  

1375. params = {}

1376. params['booster'] = 'gblinear'

1377. params['objective'] = "multi:softprob"

1378. params['eval_metric'] = 'mlogloss'

1379. params['eta'] = 0.005

1380. params['num_class'] = 12

1381. params['lambda'] = 3

1382. params['alpha'] = 2

1383.  

1384. # Random 10% for validation

1385. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

1386.  

1387. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

1388. ytr, yte = y[kf[0]], y[kf[1]]

1389.  

1390. print('Training set: ' + str(Xtr.shape))

1391. print('Validation set: ' + str(Xte.shape))

1392.  

1393. d_train = xgb.DMatrix(Xtr, label=ytr)

1394. d_valid = xgb.DMatrix(Xte, label=yte)

1395.  

1396. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

1397.  

1398. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

1399.  

1400. pred = clf.predict(xgb.DMatrix(Xtest))

1401.  

1402. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

1403. pred.head()

1404. pred.to_csv('sparse_xgb.csv', index=True)

1405.  

1406. #params['lambda'] = 1

1407. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

1408. # params['alpha'] = alpha

1409. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

1410. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

1411. #!/usr/bin/python

1412. import numpy as np

1413. import xgboost as xgb

1414. ###

1415. # advanced: customized loss function

1416. #

1417. print ('start running example to used customized objective function')

1418.  

1419. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

1420. dtest = xgb.DMatrix('../data/agaricus.txt.test')

1421.  

1422. # note: for customized objective function, we leave objective as default

1423. # note: what we are getting is margin value in prediction

1424. # you must know what you are doing

1425. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

1426. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

1427. num_round = 2

1428.  

1429. # user define objective function, given prediction, return gradient and second order gradient

1430. # this is log likelihood loss

1431. def logregobj(preds, dtrain):

1432. labels = dtrain.get_label()

1433. preds = 1.0 / (1.0 + np.exp(-preds))

1434. grad = preds - labels

1435. hess = preds * (1.0-preds)

1436. return grad, hess

1437.  

1438. # user defined evaluation function, return a pair metric_name, result

1439. # NOTE: when you do customized loss function, the default prediction value is margin

1440. # this may make builtin evaluation metric not function properly

1441. # for example, we are doing logistic loss, the prediction is score before logistic transformation

1442. # the builtin evaluation error assumes input is after logistic transformation

1443. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

1444. def evalerror(preds, dtrain):

1445. labels = dtrain.get_label()

1446. # return a pair metric_name, result

1447. # since preds are margin(before logistic transformation, cutoff at 0)

1448. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

1449.  

1450. # training with customized objective, we can also do step by step training

1451. # simply look at xgboost.py's implementation of train

1452. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

1453. 1
1454. 2
1455. 3
1456. 4
1457. 5
1458. 6
1459. 7
1460. 8
1461. 9
1462. 10
1463. 11
1464. 12
1465. 13
1466. 14
1467. 15
1468. 16
1469. 17
1470. 18
1471. 19
1472. 20
1473. 21
1474. 22
1475. 23
1476. 24
1477. 25
1478. 26
1479. 27
1480. 28
1481. 29
1482. 30
1483. 31
1484. 32
1485. 33
1486. 34
1487. 35
1488. 36
1489. 37
1490. 38
1491. 39
1492. 40
1493. 41
1494. 42
1495. 1
1496. 2
1497. 3
1498. 4
1499. 5
1500. 6
1501. 7
1502. 8
1503. 9
1504. 10
1505. 11
1506. 12
1507. 13
1508. 14
1509. 15
1510. 16
1511. 17
1512. 18
1513. 19
1514. 20
1515. 21
1516. 22
1517. 23
1518. 24
1519. 25
1520. 26
1521. 27
1522. 28
1523. 29
1524. 30
1525. 31
1526. 32
1527. 33
1528. 34
1529. 35
1530. 36
1531. 37
1532. 38
1533. 39
1534. 40
1535. 41
1536. 42

1537. 5.Xgboost调参

      由于xgboost的参数过多，这里介绍三种思路

      （1）GridSearch

      （2）Hyperopt

      （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

      6.工程实现优化

      （1）Column Blocks and Parallelization

       

      

       

      （2）Cache Aware Access

1538. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
1539. The main thread accumulates gradients statistics in the continuous buff er.

1540. （3）System Tricks

1541. Block pre-fetching.
1542. Utilize multiple disks to parallelize disk operations.
1543. LZ4 compression(popular recent years for outstanding performance).
1544. Unrolling loops.
1545. OpenMP

1546. 7.代码走读

      这块非常感谢杨军老师的无私奉献【4】

      个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

      find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

1547. 1
1548. 1

1549. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

1550. cli_main.cc:

1551. main()

1552. -> CLIRunTask()

1553. -> CLITrain()

1554. -> DMatrix::Load()

1555. -> learner = Learner::Create()

1556. -> learner->Configure()

1557. -> learner->InitModel()

1558. -> for (i = 0; i < param.num_round; ++i)

1559. -> learner->UpdateOneIter()

1560. -> learner->Save()

1561. learner.cc:

1562. Create()

1563. -> new LearnerImpl()

1564. Configure()

1565. InitModel()

1566. -> LazyInitModel()

1567. -> obj_ = ObjFunction::Create()

1568. -> objective.cc

1569. Create()

1570. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

1571. LambdaRankObj(rank_obj.cc)/

1572. RegLossObj(regression_obj.cc)/

1573. PoissonRegression(regression_obj.cc)

1574. -> gbm_ = GradientBooster::Create()

1575. -> gbm.cc

1576. Create()

1577. -> GBTree(gbtree.cc)/

1578. GBLinear(gblinear.cc)

1579. -> obj_->Configure()

1580. -> gbm_->Configure()

1581. UpdateOneIter()

1582. -> PredictRaw()

1583. -> obj_->GetGradient()

1584. -> gbm_->DoBoost()

1585.  

1586. gbtree.cc:

1587. Configure()

1588. -> for (up in updaters)

1589. -> up->Init()

1590. DoBoost()

1591. -> BoostNewTrees()

1592. -> new_tree = new RegTree()

1593. -> for (up in updaters)

1594. -> up->Update(new_tree)

1595.  

1596. tree_updater.cc:

1597. Create()

1598. -> ColMaker/DistColMaker(updater_colmaker.cc)/

1599. SketchMaker(updater_skmaker.cc)/

1600. TreeRefresher(updater_refresh.cc)/

1601. TreePruner(updater_prune.cc)/

1602. HistMaker/CQHistMaker/

1603. GlobalProposalHistMaker/

1604. QuantileHistMaker(updater_histmaker.cc)/

1605. TreeSyncher(updater_sync.cc)

1606. 1
1607. 2
1608. 3
1609. 4
1610. 5
1611. 6
1612. 7
1613. 8
1614. 9
1615. 10
1616. 11
1617. 12
1618. 13
1619. 14
1620. 15
1621. 16
1622. 17
1623. 18
1624. 19
1625. 20
1626. 21
1627. 22
1628. 23
1629. 24
1630. 25
1631. 26
1632. 27
1633. 28
1634. 29
1635. 30
1636. 31
1637. 32
1638. 33
1639. 34
1640. 35
1641. 36
1642. 37
1643. 38
1644. 39
1645. 40
1646. 41
1647. 42
1648. 43
1649. 44
1650. 45
1651. 46
1652. 47
1653. 48
1654. 49
1655. 50
1656. 51
1657. 52
1658. 53
1659. 54
1660. 55
1661. 56
1662. 1
1663. 2
1664. 3
1665. 4
1666. 5
1667. 6
1668. 7
1669. 8
1670. 9
1671. 10
1672. 11
1673. 12
1674. 13
1675. 14
1676. 15
1677. 16
1678. 17
1679. 18
1680. 19
1681. 20
1682. 21
1683. 22
1684. 23
1685. 24
1686. 25
1687. 26
1688. 27
1689. 28
1690. 29
1691. 30
1692. 31
1693. 32
1694. 33
1695. 34
1696. 35
1697. 36
1698. 37
1699. 38
1700. 39
1701. 40
1702. 41
1703. 42
1704. 43
1705. 44
1706. 45
1707. 46
1708. 47
1709. 48
1710. 49
1711. 50
1712. 51
1713. 52
1714. 53
1715. 54
1716. 55
1717. 56

1718. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

      I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
      II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
      III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
      　　I). inter-feature exact parallelism （特征级精确并行） 
      　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
      　　III). intra-feature parallelism （特征内并行）

      此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

      XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

      以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

1719. updater_colmaker.cc:

1720. ColMaker::Update()

1721. -> Builder builder;

1722. -> builder.Update()

1723. -> InitData()

1724. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

1725. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

1726. -> for (depth = 0; depth < 树的最大深度; ++depth)

1727. -> FindSplit()

1728. -> for (each feature) // 通过OpenMP获取

1729. // inter-feature parallelism

1730. -> UpdateSolution()

1731. -> EnumerateSplit() // 每个执行线程处理一个特征，

1732. // 选出每个特征的

1733. // 最优split point

1734. -> ParallelFindSplit()

1735. // 多个执行线程同时处理一个特征，选出该特征

1736. //的最优split point;

1737. // 在每个线程里汇总各个线程内分配到的数据样

1738. //本的统计量(grad/hess);

1739. // aggregate所有线程的样本统计(grad/hess)，

1740. //计算出每个线程分配到的样本集合的边界特征值作为

1741. //split point的最优分割点;

1742. // 在每个线程分配到的样本集合对应的特征值集合进

1743. //行枚举作为split point，选出最优分割点

1744. -> SyncBestSolution()

1745. // 上面的UpdateSolution()/ParallelFindSplit()

1746. //会为所有待扩展分割的叶结点找到特征维度的最优split

1747. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

1748. //的最优split point，OpenMP线程2会找到特征F2的最

1749. //优split point，所以需要进行全局sync，找到叶结点A

1750. //的最优split point。

1751. -> 为需要进行分割的叶结点创建孩子结点

1752. -> ResetPosition()

1753. //根据上一步的分割动作，更新样本到树结点的映射关系

1754. // Missing Value(i.e. default)和非Missing Value(i.e.

1755. //non-default)分别处理

1756. -> UpdateQueueExpand()

1757. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

1758. //起始基础

1759. -> InitNewNode() // 为可用于split的树结点计算统计量\

1760. 8.python、R对于xgboost的简单使用

      任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

      【Python】

       

      

       

      【R】

       

      

       

      9.xgboost中比较重要的参数介绍

      （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

1761. “reg:linear” –线性回归。
1762. “reg:logistic” –逻辑回归。
1763. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
1764. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
1765. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
1766. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
1767. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
1768. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

1769. （2）’eval_metric’ The choices are listed below，评估指标:

1770. “rmse”: root mean square error
1771. “logloss”: negative log-likelihood
1772. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
1773. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
1774. “mlogloss”: Multiclass logloss
1775. “auc”: Area under the curve for ranking evaluation.
1776. “ndcg”:Normalized Discounted Cumulative Gain
1777. “map”:Mean average precision
1778. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
1779. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

1780. （3）lambda [default=0] L2 正则的惩罚系数

      （4）alpha [default=0] L1 正则的惩罚系数

      （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

      （6）eta [default=0.3] 
      为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
      取值范围为：[0,1]

      （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

      （8）min_child_weight [default=1] 
      孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
      取值范围为: [0,∞]

      10.DART

      核心思想就是将dropout引入XGBoost

      示例代码

1781. import xgboost as xgb

1782. # read in data

1783. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

1784. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

1785. # specify parameters via map

1786. param = {'booster': 'dart',

1787. 'max_depth': 5, 'learning_rate': 0.1,

1788. 'objective': 'binary:logistic', 'silent': True,

1789. 'sample_type': 'uniform',

1790. 'normalize_type': 'tree',

1791. 'rate_drop': 0.1,

1792. 'skip_drop': 0.5}

1793. num_round = 50

1794. bst = xgb.train(param, dtrain, num_round)

1795. # make prediction

1796. # ntree_limit must not be 0

1797. preds = bst.predict(dtest, ntree_limit=num_round

1798. 更多细节可以阅读参考文献5

      11.csr_matrix训练XGBoost

      当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

      下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

1799. import pandas as pd

1800. import numpy as np

1801. import seaborn as sns

1802. import matplotlib.pyplot as plt

1803. import os

1804. from sklearn.preprocessing import LabelEncoder

1805. from scipy.sparse import csr_matrix, hstack

1806. import xgboost as xgb

1807. from sklearn.cross_validation import StratifiedKFold

1808. from sklearn.metrics import log_loss

1809.  

1810. datadir = '../input'

1811. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

1812. index_col='device_id')

1813. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

1814. index_col = 'device_id')

1815. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

1816. # Get rid of duplicate device ids in phone

1817. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

1818. events = pd.read_csv(os.path.join(datadir,'events.csv'),

1819. parse_dates=['timestamp'], index_col='event_id')

1820. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

1821. usecols=['event_id','app_id','is_active'],

1822. dtype={'is_active':bool})

1823. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

1824.  

1825. gatrain['trainrow'] = np.arange(gatrain.shape[0])

1826. gatest['testrow'] = np.arange(gatest.shape[0])

1827.  

1828. brandencoder = LabelEncoder().fit(phone.phone_brand)

1829. phone['brand'] = brandencoder.transform(phone['phone_brand'])

1830. gatrain['brand'] = phone['brand']

1831. gatest['brand'] = phone['brand']

1832. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

1833. (gatrain.trainrow, gatrain.brand)))

1834. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

1835. (gatest.testrow, gatest.brand)))

1836. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

1837.  

1838. m = phone.phone_brand.str.cat(phone.device_model)

1839. modelencoder = LabelEncoder().fit(m)

1840. phone['model'] = modelencoder.transform(m)

1841. gatrain['model'] = phone['model']

1842. gatest['model'] = phone['model']

1843. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

1844. (gatrain.trainrow, gatrain.model)))

1845. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

1846. (gatest.testrow, gatest.model)))

1847. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

1848.  

1849. appencoder = LabelEncoder().fit(appevents.app_id)

1850. appevents['app'] = appencoder.transform(appevents.app_id)

1851. napps = len(appencoder.classes_)

1852. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

1853. .groupby(['device_id','app'])['app'].agg(['size'])

1854. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

1855. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

1856. .reset_index())

1857.  

1858. d = deviceapps.dropna(subset=['trainrow'])

1859. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

1860. shape=(gatrain.shape[0],napps))

1861. d = deviceapps.dropna(subset=['testrow'])

1862. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

1863. shape=(gatest.shape[0],napps))

1864. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

1865.  

1866. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

1867. applabels['app'] = appencoder.transform(applabels.app_id)

1868. labelencoder = LabelEncoder().fit(applabels.label_id)

1869. applabels['label'] = labelencoder.transform(applabels.label_id)

1870. nlabels = len(labelencoder.classes_)

1871.  

1872. devicelabels = (deviceapps[['device_id','app']]

1873. .merge(applabels[['app','label']])

1874. .groupby(['device_id','label'])['app'].agg(['size'])

1875. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

1876. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

1877. .reset_index())

1878. devicelabels.head()

1879.  

1880. d = devicelabels.dropna(subset=['trainrow'])

1881. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

1882. shape=(gatrain.shape[0],nlabels))

1883. d = devicelabels.dropna(subset=['testrow'])

1884. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

1885. shape=(gatest.shape[0],nlabels))

1886. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

1887.  

1888. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

1889. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

1890. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

1891.  

1892. targetencoder = LabelEncoder().fit(gatrain.group)

1893. y = targetencoder.transform(gatrain.group)

1894.  

1895. ########## XGBOOST ##########

1896.  

1897. params = {}

1898. params['booster'] = 'gblinear'

1899. params['objective'] = "multi:softprob"

1900. params['eval_metric'] = 'mlogloss'

1901. params['eta'] = 0.005

1902. params['num_class'] = 12

1903. params['lambda'] = 3

1904. params['alpha'] = 2

1905.  

1906. # Random 10% for validation

1907. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

1908.  

1909. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

1910. ytr, yte = y[kf[0]], y[kf[1]]

1911.  

1912. print('Training set: ' + str(Xtr.shape))

1913. print('Validation set: ' + str(Xte.shape))

1914.  

1915. d_train = xgb.DMatrix(Xtr, label=ytr)

1916. d_valid = xgb.DMatrix(Xte, label=yte)

1917.  

1918. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

1919.  

1920. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

1921.  

1922. pred = clf.predict(xgb.DMatrix(Xtest))

1923.  

1924. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

1925. pred.head()

1926. pred.to_csv('sparse_xgb.csv', index=True)

1927.  

1928. #params['lambda'] = 1

1929. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

1930. # params['alpha'] = alpha

1931. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

1932. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

1933. #!/usr/bin/python

1934. import numpy as np

1935. import xgboost as xgb

1936. ###

1937. # advanced: customized loss function

1938. #

1939. print ('start running example to used customized objective function')

1940.  

1941. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

1942. dtest = xgb.DMatrix('../data/agaricus.txt.test')

1943.  

1944. # note: for customized objective function, we leave objective as default

1945. # note: what we are getting is margin value in prediction

1946. # you must know what you are doing

1947. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

1948. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

1949. num_round = 2

1950.  

1951. # user define objective function, given prediction, return gradient and second order gradient

1952. # this is log likelihood loss

1953. def logregobj(preds, dtrain):

1954. labels = dtrain.get_label()

1955. preds = 1.0 / (1.0 + np.exp(-preds))

1956. grad = preds - labels

1957. hess = preds * (1.0-preds)

1958. return grad, hess

1959.  

1960. # user defined evaluation function, return a pair metric_name, result

1961. # NOTE: when you do customized loss function, the default prediction value is margin

1962. # this may make builtin evaluation metric not function properly

1963. # for example, we are doing logistic loss, the prediction is score before logistic transformation

1964. # the builtin evaluation error assumes input is after logistic transformation

1965. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

1966. def evalerror(preds, dtrain):

1967. labels = dtrain.get_label()

1968. # return a pair metric_name, result

1969. # since preds are margin(before logistic transformation, cutoff at 0)

1970. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

1971.  

1972. # training with customized objective, we can also do step by step training

1973. # simply look at xgboost.py's implementation of train

1974. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

1975. 1
1976. 2
1977. 3
1978. 4
1979. 5
1980. 6
1981. 7
1982. 8
1983. 9
1984. 10
1985. 11
1986. 12
1987. 13
1988. 14
1989. 15
1990. 16
1991. 17
1992. 18
1993. 19
1994. 20
1995. 21
1996. 22
1997. 23
1998. 24
1999. 25
2000. 26
2001. 27
2002. 28
2003. 29
2004. 30
2005. 31
2006. 32
2007. 33
2008. 34
2009. 35
2010. 36
2011. 37
2012. 38
2013. 39
2014. 40
2015. 41
2016. 42
2017. 1
2018. 2
2019. 3
2020. 4
2021. 5
2022. 6
2023. 7
2024. 8
2025. 9
2026. 10
2027. 11
2028. 12
2029. 13
2030. 14
2031. 15
2032. 16
2033. 17
2034. 18
2035. 19
2036. 20
2037. 21
2038. 22
2039. 23
2040. 24
2041. 25
2042. 26
2043. 27
2044. 28
2045. 29
2046. 30
2047. 31
2048. 32
2049. 33
2050. 34
2051. 35
2052. 36
2053. 37
2054. 38
2055. 39
2056. 40
2057. 41
2058. 42

2059. 5.Xgboost调参

      由于xgboost的参数过多，这里介绍三种思路

      （1）GridSearch

      （2）Hyperopt

      （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

      6.工程实现优化

      （1）Column Blocks and Parallelization

       

      

       

      （2）Cache Aware Access

2060. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
2061. The main thread accumulates gradients statistics in the continuous buff er.

2062. （3）System Tricks

2063. Block pre-fetching.
2064. Utilize multiple disks to parallelize disk operations.
2065. LZ4 compression(popular recent years for outstanding performance).
2066. Unrolling loops.
2067. OpenMP

2068. 7.代码走读

      这块非常感谢杨军老师的无私奉献【4】

      个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

      find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

2069. 1
2070. 1

2071. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

2072. cli_main.cc:

2073. main()

2074. -> CLIRunTask()

2075. -> CLITrain()

2076. -> DMatrix::Load()

2077. -> learner = Learner::Create()

2078. -> learner->Configure()

2079. -> learner->InitModel()

2080. -> for (i = 0; i < param.num_round; ++i)

2081. -> learner->UpdateOneIter()

2082. -> learner->Save()

2083. learner.cc:

2084. Create()

2085. -> new LearnerImpl()

2086. Configure()

2087. InitModel()

2088. -> LazyInitModel()

2089. -> obj_ = ObjFunction::Create()

2090. -> objective.cc

2091. Create()

2092. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

2093. LambdaRankObj(rank_obj.cc)/

2094. RegLossObj(regression_obj.cc)/

2095. PoissonRegression(regression_obj.cc)

2096. -> gbm_ = GradientBooster::Create()

2097. -> gbm.cc

2098. Create()

2099. -> GBTree(gbtree.cc)/

2100. GBLinear(gblinear.cc)

2101. -> obj_->Configure()

2102. -> gbm_->Configure()

2103. UpdateOneIter()

2104. -> PredictRaw()

2105. -> obj_->GetGradient()

2106. -> gbm_->DoBoost()

2107.  

2108. gbtree.cc:

2109. Configure()

2110. -> for (up in updaters)

2111. -> up->Init()

2112. DoBoost()

2113. -> BoostNewTrees()

2114. -> new_tree = new RegTree()

2115. -> for (up in updaters)

2116. -> up->Update(new_tree)

2117.  

2118. tree_updater.cc:

2119. Create()

2120. -> ColMaker/DistColMaker(updater_colmaker.cc)/

2121. SketchMaker(updater_skmaker.cc)/

2122. TreeRefresher(updater_refresh.cc)/

2123. TreePruner(updater_prune.cc)/

2124. HistMaker/CQHistMaker/

2125. GlobalProposalHistMaker/

2126. QuantileHistMaker(updater_histmaker.cc)/

2127. TreeSyncher(updater_sync.cc)

2128. 1
2129. 2
2130. 3
2131. 4
2132. 5
2133. 6
2134. 7
2135. 8
2136. 9
2137. 10
2138. 11
2139. 12
2140. 13
2141. 14
2142. 15
2143. 16
2144. 17
2145. 18
2146. 19
2147. 20
2148. 21
2149. 22
2150. 23
2151. 24
2152. 25
2153. 26
2154. 27
2155. 28
2156. 29
2157. 30
2158. 31
2159. 32
2160. 33
2161. 34
2162. 35
2163. 36
2164. 37
2165. 38
2166. 39
2167. 40
2168. 41
2169. 42
2170. 43
2171. 44
2172. 45
2173. 46
2174. 47
2175. 48
2176. 49
2177. 50
2178. 51
2179. 52
2180. 53
2181. 54
2182. 55
2183. 56
2184. 1
2185. 2
2186. 3
2187. 4
2188. 5
2189. 6
2190. 7
2191. 8
2192. 9
2193. 10
2194. 11
2195. 12
2196. 13
2197. 14
2198. 15
2199. 16
2200. 17
2201. 18
2202. 19
2203. 20
2204. 21
2205. 22
2206. 23
2207. 24
2208. 25
2209. 26
2210. 27
2211. 28
2212. 29
2213. 30
2214. 31
2215. 32
2216. 33
2217. 34
2218. 35
2219. 36
2220. 37
2221. 38
2222. 39
2223. 40
2224. 41
2225. 42
2226. 43
2227. 44
2228. 45
2229. 46
2230. 47
2231. 48
2232. 49
2233. 50
2234. 51
2235. 52
2236. 53
2237. 54
2238. 55
2239. 56

2240. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

      I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
      II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
      III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
      　　I). inter-feature exact parallelism （特征级精确并行） 
      　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
      　　III). intra-feature parallelism （特征内并行）

      此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

      XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

      以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

2241. updater_colmaker.cc:

2242. ColMaker::Update()

2243. -> Builder builder;

2244. -> builder.Update()

2245. -> InitData()

2246. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

2247. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

2248. -> for (depth = 0; depth < 树的最大深度; ++depth)

2249. -> FindSplit()

2250. -> for (each feature) // 通过OpenMP获取

2251. // inter-feature parallelism

2252. -> UpdateSolution()

2253. -> EnumerateSplit() // 每个执行线程处理一个特征，

2254. // 选出每个特征的

2255. // 最优split point

2256. -> ParallelFindSplit()

2257. // 多个执行线程同时处理一个特征，选出该特征

2258. //的最优split point;

2259. // 在每个线程里汇总各个线程内分配到的数据样

2260. //本的统计量(grad/hess);

2261. // aggregate所有线程的样本统计(grad/hess)，

2262. //计算出每个线程分配到的样本集合的边界特征值作为

2263. //split point的最优分割点;

2264. // 在每个线程分配到的样本集合对应的特征值集合进

2265. //行枚举作为split point，选出最优分割点

2266. -> SyncBestSolution()

2267. // 上面的UpdateSolution()/ParallelFindSplit()

2268. //会为所有待扩展分割的叶结点找到特征维度的最优split

2269. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

2270. //的最优split point，OpenMP线程2会找到特征F2的最

2271. //优split point，所以需要进行全局sync，找到叶结点A

2272. //的最优split point。

2273. -> 为需要进行分割的叶结点创建孩子结点

2274. -> ResetPosition()

2275. //根据上一步的分割动作，更新样本到树结点的映射关系

2276. // Missing Value(i.e. default)和非Missing Value(i.e.

2277. //non-default)分别处理

2278. -> UpdateQueueExpand()

2279. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

2280. //起始基础

2281. -> InitNewNode() // 为可用于split的树结点计算统计量

2282. 1
2283. 2
2284. 3
2285. 4
2286. 5
2287. 6
2288. 7
2289. 8
2290. 9
2291. 10
2292. 11
2293. 12
2294. 13
2295. 14
2296. 15
2297. 16
2298. 17
2299. 18
2300. 19
2301. 20
2302. 21
2303. 22
2304. 23
2305. 24
2306. 25
2307. 26
2308. 27
2309. 28
2310. 29
2311. 30
2312. 31
2313. 32
2314. 33
2315. 34
2316. 35
2317. 36
2318. 37
2319. 38
2320. 39
2321. 40
2322. 41
2323. 1
2324. 2
2325. 3
2326. 4
2327. 5
2328. 6
2329. 7
2330. 8
2331. 9
2332. 10
2333. 11
2334. 12
2335. 13
2336. 14
2337. 15
2338. 16
2339. 17
2340. 18
2341. 19
2342. 20
2343. 21
2344. 22
2345. 23
2346. 24
2347. 25
2348. 26
2349. 27
2350. 28
2351. 29
2352. 30
2353. 31
2354. 32
2355. 33
2356. 34
2357. 35
2358. 36
2359. 37
2360. 38
2361. 39
2362. 40
2363. 41

2364. 8.python、R对于xgboost的简单使用

      任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

      【Python】

       

      

       

      【R】

       

      

       

      9.xgboost中比较重要的参数介绍

      （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

2365. “reg:linear” –线性回归。
2366. “reg:logistic” –逻辑回归。
2367. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
2368. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
2369. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
2370. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
2371. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
2372. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

2373. （2）’eval_metric’ The choices are listed below，评估指标:

2374. “rmse”: root mean square error
2375. “logloss”: negative log-likelihood
2376. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
2377. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
2378. “mlogloss”: Multiclass logloss
2379. “auc”: Area under the curve for ranking evaluation.
2380. “ndcg”:Normalized Discounted Cumulative Gain
2381. “map”:Mean average precision
2382. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
2383. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

2384. （3）lambda [default=0] L2 正则的惩罚系数

      （4）alpha [default=0] L1 正则的惩罚系数

      （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

      （6）eta [default=0.3] 
      为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
      取值范围为：[0,1]

      （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

      （8）min_child_weight [default=1] 
      孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
      取值范围为: [0,∞]

      10.DART

      核心思想就是将dropout引入XGBoost

      示例代码

2385. import xgboost as xgb

2386. # read in data

2387. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

2388. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

2389. # specify parameters via map

2390. param = {'booster': 'dart',

2391. 'max_depth': 5, 'learning_rate': 0.1,

2392. 'objective': 'binary:logistic', 'silent': True,

2393. 'sample_type': 'uniform',

2394. 'normalize_type': 'tree',

2395. 'rate_drop': 0.1,

2396. 'skip_drop': 0.5}

2397. num_round = 50

2398. bst = xgb.train(param, dtrain, num_round)

2399. # make prediction

2400. # ntree_limit must not be 0

2401. preds = bst.predict(dtest, ntree_limit=num_round)

2402.  

2403. 更多细节可以阅读参考文献5

      11.csr_matrix训练XGBoost

      当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

      下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

2404. import pandas as pd

2405. import numpy as np

2406. import seaborn as sns

2407. import matplotlib.pyplot as plt

2408. import os

2409. from sklearn.preprocessing import LabelEncoder

2410. from scipy.sparse import csr_matrix, hstack

2411. import xgboost as xgb

2412. from sklearn.cross_validation import StratifiedKFold

2413. from sklearn.metrics import log_loss

2414.  

2415. datadir = '../input'

2416. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

2417. index_col='device_id')

2418. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

2419. index_col = 'device_id')

2420. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

2421. # Get rid of duplicate device ids in phone

2422. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

2423. events = pd.read_csv(os.path.join(datadir,'events.csv'),

2424. parse_dates=['timestamp'], index_col='event_id')

2425. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

2426. usecols=['event_id','app_id','is_active'],

2427. dtype={'is_active':bool})

2428. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

2429.  

2430. gatrain['trainrow'] = np.arange(gatrain.shape[0])

2431. gatest['testrow'] = np.arange(gatest.shape[0])

2432.  

2433. brandencoder = LabelEncoder().fit(phone.phone_brand)

2434. phone['brand'] = brandencoder.transform(phone['phone_brand'])

2435. gatrain['brand'] = phone['brand']

2436. gatest['brand'] = phone['brand']

2437. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

2438. (gatrain.trainrow, gatrain.brand)))

2439. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

2440. (gatest.testrow, gatest.brand)))

2441. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

2442.  

2443. m = phone.phone_brand.str.cat(phone.device_model)

2444. modelencoder = LabelEncoder().fit(m)

2445. phone['model'] = modelencoder.transform(m)

2446. gatrain['model'] = phone['model']

2447. gatest['model'] = phone['model']

2448. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

2449. (gatrain.trainrow, gatrain.model)))

2450. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

2451. (gatest.testrow, gatest.model)))

2452. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

2453.  

2454. appencoder = LabelEncoder().fit(appevents.app_id)

2455. appevents['app'] = appencoder.transform(appevents.app_id)

2456. napps = len(appencoder.classes_)

2457. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

2458. .groupby(['device_id','app'])['app'].agg(['size'])

2459. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

2460. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

2461. .reset_index())

2462.  

2463. d = deviceapps.dropna(subset=['trainrow'])

2464. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

2465. shape=(gatrain.shape[0],napps))

2466. d = deviceapps.dropna(subset=['testrow'])

2467. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

2468. shape=(gatest.shape[0],napps))

2469. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

2470.  

2471. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

2472. applabels['app'] = appencoder.transform(applabels.app_id)

2473. labelencoder = LabelEncoder().fit(applabels.label_id)

2474. applabels['label'] = labelencoder.transform(applabels.label_id)

2475. nlabels = len(labelencoder.classes_)

2476.  

2477. devicelabels = (deviceapps[['device_id','app']]

2478. .merge(applabels[['app','label']])

2479. .groupby(['device_id','label'])['app'].agg(['size'])

2480. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

2481. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

2482. .reset_index())

2483. devicelabels.head()

2484.  

2485. d = devicelabels.dropna(subset=['trainrow'])

2486. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

2487. shape=(gatrain.shape[0],nlabels))

2488. d = devicelabels.dropna(subset=['testrow'])

2489. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

2490. shape=(gatest.shape[0],nlabels))

2491. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

2492.  

2493. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

2494. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

2495. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

2496.  

2497. targetencoder = LabelEncoder().fit(gatrain.group)

2498. y = targetencoder.transform(gatrain.group)

2499.  

2500. ########## XGBOOST ##########

2501.  

2502. params = {}

2503. params['booster'] = 'gblinear'

2504. params['objective'] = "multi:softprob"

2505. params['eval_metric'] = 'mlogloss'

2506. params['eta'] = 0.005

2507. params['num_class'] = 12

2508. params['lambda'] = 3

2509. params['alpha'] = 2

2510.  

2511. # Random 10% for validation

2512. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

2513.  

2514. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

2515. ytr, yte = y[kf[0]], y[kf[1]]

2516.  

2517. print('Training set: ' + str(Xtr.shape))

2518. print('Validation set: ' + str(Xte.shape))

2519.  

2520. d_train = xgb.DMatrix(Xtr, label=ytr)

2521. d_valid = xgb.DMatrix(Xte, label=yte)

2522.  

2523. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

2524.  

2525. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

2526.  

2527. pred = clf.predict(xgb.DMatrix(Xtest))

2528.  

2529. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

2530. pred.head()

2531. pred.to_csv('sparse_xgb.csv', index=True)

2532.  

2533. #params['lambda'] = 1

2534. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

2535. # params['alpha'] = alpha

2536. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

2537. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

2538. #!/usr/bin/python

2539. import numpy as np

2540. import xgboost as xgb

2541. ###

2542. # advanced: customized loss function

2543. #

2544. print ('start running example to used customized objective function')

2545.  

2546. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

2547. dtest = xgb.DMatrix('../data/agaricus.txt.test')

2548.  

2549. # note: for customized objective function, we leave objective as default

2550. # note: what we are getting is margin value in prediction

2551. # you must know what you are doing

2552. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

2553. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

2554. num_round = 2

2555.  

2556. # user define objective function, given prediction, return gradient and second order gradient

2557. # this is log likelihood loss

2558. def logregobj(preds, dtrain):

2559. labels = dtrain.get_label()

2560. preds = 1.0 / (1.0 + np.exp(-preds))

2561. grad = preds - labels

2562. hess = preds * (1.0-preds)

2563. return grad, hess

2564.  

2565. # user defined evaluation function, return a pair metric_name, result

2566. # NOTE: when you do customized loss function, the default prediction value is margin

2567. # this may make builtin evaluation metric not function properly

2568. # for example, we are doing logistic loss, the prediction is score before logistic transformation

2569. # the builtin evaluation error assumes input is after logistic transformation

2570. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

2571. def evalerror(preds, dtrain):

2572. labels = dtrain.get_label()

2573. # return a pair metric_name, result

2574. # since preds are margin(before logistic transformation, cutoff at 0)

2575. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

2576.  

2577. # training with customized objective, we can also do step by step training

2578. # simply look at xgboost.py's implementation of train

2579. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

2580.  

2581. 5.Xgboost调参

      由于xgboost的参数过多，这里介绍三种思路

      （1）GridSearch

      （2）Hyperopt

      （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

      6.工程实现优化

      （1）Column Blocks and Parallelization

       

      

       

      （2）Cache Aware Access

2582. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
2583. The main thread accumulates gradients statistics in the continuous buff er.

2584. （3）System Tricks

2585. Block pre-fetching.
2586. Utilize multiple disks to parallelize disk operations.
2587. LZ4 compression(popular recent years for outstanding performance).
2588. Unrolling loops.
2589. OpenMP

2590. 7.代码走读

      这块非常感谢杨军老师的无私奉献【4】

      个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

      find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

2591. 1
2592. 1

2593. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

2594. cli_main.cc:

2595. main()

2596. -> CLIRunTask()

2597. -> CLITrain()

2598. -> DMatrix::Load()

2599. -> learner = Learner::Create()

2600. -> learner->Configure()

2601. -> learner->InitModel()

2602. -> for (i = 0; i < param.num_round; ++i)

2603. -> learner->UpdateOneIter()

2604. -> learner->Save()

2605. learner.cc:

2606. Create()

2607. -> new LearnerImpl()

2608. Configure()

2609. InitModel()

2610. -> LazyInitModel()

2611. -> obj_ = ObjFunction::Create()

2612. -> objective.cc

2613. Create()

2614. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

2615. LambdaRankObj(rank_obj.cc)/

2616. RegLossObj(regression_obj.cc)/

2617. PoissonRegression(regression_obj.cc)

2618. -> gbm_ = GradientBooster::Create()

2619. -> gbm.cc

2620. Create()

2621. -> GBTree(gbtree.cc)/

2622. GBLinear(gblinear.cc)

2623. -> obj_->Configure()

2624. -> gbm_->Configure()

2625. UpdateOneIter()

2626. -> PredictRaw()

2627. -> obj_->GetGradient()

2628. -> gbm_->DoBoost()

2629.  

2630. gbtree.cc:

2631. Configure()

2632. -> for (up in updaters)

2633. -> up->Init()

2634. DoBoost()

2635. -> BoostNewTrees()

2636. -> new_tree = new RegTree()

2637. -> for (up in updaters)

2638. -> up->Update(new_tree)

2639.  

2640. tree_updater.cc:

2641. Create()

2642. -> ColMaker/DistColMaker(updater_colmaker.cc)/

2643. SketchMaker(updater_skmaker.cc)/

2644. TreeRefresher(updater_refresh.cc)/

2645. TreePruner(updater_prune.cc)/

2646. HistMaker/CQHistMaker/

2647. GlobalProposalHistMaker/

2648. QuantileHistMaker(updater_histmaker.cc)/

2649. TreeSyncher(updater_sync.cc)

2650.  

2651. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

      I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
      II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
      III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
      　　I). inter-feature exact parallelism （特征级精确并行） 
      　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
      　　III). intra-feature parallelism （特征内并行）

      此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

      XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

      以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

2652. updater_colmaker.cc:

2653. ColMaker::Update()

2654. -> Builder builder;

2655. -> builder.Update()

2656. -> InitData()

2657. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

2658. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

2659. -> for (depth = 0; depth < 树的最大深度; ++depth)

2660. -> FindSplit()

2661. -> for (each feature) // 通过OpenMP获取

2662. // inter-feature parallelism

2663. -> UpdateSolution()

2664. -> EnumerateSplit() // 每个执行线程处理一个特征，

2665. // 选出每个特征的

2666. // 最优split point

2667. -> ParallelFindSplit()

2668. // 多个执行线程同时处理一个特征，选出该特征

2669. //的最优split point;

2670. // 在每个线程里汇总各个线程内分配到的数据样

2671. //本的统计量(grad/hess);

2672. // aggregate所有线程的样本统计(grad/hess)，

2673. //计算出每个线程分配到的样本集合的边界特征值作为

2674. //split point的最优分割点;

2675. // 在每个线程分配到的样本集合对应的特征值集合进

2676. //行枚举作为split point，选出最优分割点

2677. -> SyncBestSolution()

2678. // 上面的UpdateSolution()/ParallelFindSplit()

2679. //会为所有待扩展分割的叶结点找到特征维度的最优split

2680. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

2681. //的最优split point，OpenMP线程2会找到特征F2的最

2682. //优split point，所以需要进行全局sync，找到叶结点A

2683. //的最优split point。

2684. -> 为需要进行分割的叶结点创建孩子结点

2685. -> ResetPosition()

2686. //根据上一步的分割动作，更新样本到树结点的映射关系

2687. // Missing Value(i.e. default)和非Missing Value(i.e.

2688. //non-default)分别处理

2689. -> UpdateQueueExpand()

2690. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

2691. //起始基础

2692. -> InitNewNode() // 为可用于split的树结点计算统计量

2693.  

2694. 8.python、R对于xgboost的简单使用

      任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

      【Python】

       

      

       

      【R】

       

      

       

      9.xgboost中比较重要的参数介绍

      （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

2695. “reg:linear” –线性回归。
2696. “reg:logistic” –逻辑回归。
2697. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
2698. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
2699. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
2700. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
2701. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
2702. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

2703. （2）’eval_metric’ The choices are listed below，评估指标:

2704. “rmse”: root mean square error
2705. “logloss”: negative log-likelihood
2706. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
2707. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
2708. “mlogloss”: Multiclass logloss
2709. “auc”: Area under the curve for ranking evaluation.
2710. “ndcg”:Normalized Discounted Cumulative Gain
2711. “map”:Mean average precision
2712. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
2713. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

2714. （3）lambda [default=0] L2 正则的惩罚系数

      （4）alpha [default=0] L1 正则的惩罚系数

      （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

      （6）eta [default=0.3] 
      为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
      取值范围为：[0,1]

      （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

      （8）min_child_weight [default=1] 
      孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
      取值范围为: [0,∞]

      10.DART

      核心思想就是将dropout引入XGBoost

      示例代码

2715. import xgboost as xgb

2716. # read in data

2717. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

2718. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

2719. # specify parameters via map

2720. param = {'booster': 'dart',

2721. 'max_depth': 5, 'learning_rate': 0.1,

2722. 'objective': 'binary:logistic', 'silent': True,

2723. 'sample_type': 'uniform',

2724. 'normalize_type': 'tree',

2725. 'rate_drop': 0.1,

2726. 'skip_drop': 0.5}

2727. num_round = 50

2728. bst = xgb.train(param, dtrain, num_round)

2729. # make prediction

2730. # ntree_limit must not be 0

2731. preds = bst.predict(dtest, ntree_limit=num_round)

2732.  

2733. 更多细节可以阅读参考文献5

      11.csr_matrix训练XGBoost

      当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

      下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

2734. import pandas as pd

2735. import numpy as np

2736. import seaborn as sns

2737. import matplotlib.pyplot as plt

2738. import os

2739. from sklearn.preprocessing import LabelEncoder

2740. from scipy.sparse import csr_matrix, hstack

2741. import xgboost as xgb

2742. from sklearn.cross_validation import StratifiedKFold

2743. from sklearn.metrics import log_loss

2744.  

2745. datadir = '../input'

2746. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

2747. index_col='device_id')

2748. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

2749. index_col = 'device_id')

2750. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

2751. # Get rid of duplicate device ids in phone

2752. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

2753. events = pd.read_csv(os.path.join(datadir,'events.csv'),

2754. parse_dates=['timestamp'], index_col='event_id')

2755. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

2756. usecols=['event_id','app_id','is_active'],

2757. dtype={'is_active':bool})

2758. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

2759.  

2760. gatrain['trainrow'] = np.arange(gatrain.shape[0])

2761. gatest['testrow'] = np.arange(gatest.shape[0])

2762.  

2763. brandencoder = LabelEncoder().fit(phone.phone_brand)

2764. phone['brand'] = brandencoder.transform(phone['phone_brand'])

2765. gatrain['brand'] = phone['brand']

2766. gatest['brand'] = phone['brand']

2767. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

2768. (gatrain.trainrow, gatrain.brand)))

2769. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

2770. (gatest.testrow, gatest.brand)))

2771. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

2772.  

2773. m = phone.phone_brand.str.cat(phone.device_model)

2774. modelencoder = LabelEncoder().fit(m)

2775. phone['model'] = modelencoder.transform(m)

2776. gatrain['model'] = phone['model']

2777. gatest['model'] = phone['model']

2778. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

2779. (gatrain.trainrow, gatrain.model)))

2780. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

2781. (gatest.testrow, gatest.model)))

2782. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

2783.  

2784. appencoder = LabelEncoder().fit(appevents.app_id)

2785. appevents['app'] = appencoder.transform(appevents.app_id)

2786. napps = len(appencoder.classes_)

2787. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

2788. .groupby(['device_id','app'])['app'].agg(['size'])

2789. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

2790. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

2791. .reset_index())

2792.  

2793. d = deviceapps.dropna(subset=['trainrow'])

2794. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

2795. shape=(gatrain.shape[0],napps))

2796. d = deviceapps.dropna(subset=['testrow'])

2797. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

2798. shape=(gatest.shape[0],napps))

2799. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

2800.  

2801. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

2802. applabels['app'] = appencoder.transform(applabels.app_id)

2803. labelencoder = LabelEncoder().fit(applabels.label_id)

2804. applabels['label'] = labelencoder.transform(applabels.label_id)

2805. nlabels = len(labelencoder.classes_)

2806.  

2807. devicelabels = (deviceapps[['device_id','app']]

2808. .merge(applabels[['app','label']])

2809. .groupby(['device_id','label'])['app'].agg(['size'])

2810. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

2811. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

2812. .reset_index())

2813. devicelabels.head()

2814.  

2815. d = devicelabels.dropna(subset=['trainrow'])

2816. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

2817. shape=(gatrain.shape[0],nlabels))

2818. d = devicelabels.dropna(subset=['testrow'])

2819. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

2820. shape=(gatest.shape[0],nlabels))

2821. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

2822.  

2823. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

2824. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

2825. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

2826.  

2827. targetencoder = LabelEncoder().fit(gatrain.group)

2828. y = targetencoder.transform(gatrain.group)

2829.  

2830. ########## XGBOOST ##########

2831.  

2832. params = {}

2833. params['booster'] = 'gblinear'

2834. params['objective'] = "multi:softprob"

2835. params['eval_metric'] = 'mlogloss'

2836. params['eta'] = 0.005

2837. params['num_class'] = 12

2838. params['lambda'] = 3

2839. params['alpha'] = 2

2840.  

2841. # Random 10% for validation

2842. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

2843.  

2844. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

2845. ytr, yte = y[kf[0]], y[kf[1]]

2846.  

2847. print('Training set: ' + str(Xtr.shape))

2848. print('Validation set: ' + str(Xte.shape))

2849.  

2850. d_train = xgb.DMatrix(Xtr, label=ytr)

2851. d_valid = xgb.DMatrix(Xte, label=yte)

2852.  

2853. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

2854.  

2855. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

2856.  

2857. pred = clf.predict(xgb.DMatrix(Xtest))

2858.  

2859. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

2860. pred.head()

2861. pred.to_csv('sparse_xgb.csv', index=True)

2862.  

2863. #params['lambda'] = 1

2864. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

2865. # params['alpha'] = alpha

2866. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

2867. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

2868. #!/usr/bin/python

2869. import numpy as np

2870. import xgboost as xgb

2871. ###

2872. # advanced: customized loss function

2873. #

2874. print ('start running example to used customized objective function')

2875.  

2876. dtrain = xgb.DMatrix('../data/agaricus.txt.train')

2877. dtest = xgb.DMatrix('../data/agaricus.txt.test')

2878.  

2879. # note: for customized objective function, we leave objective as default

2880. # note: what we are getting is margin value in prediction

2881. # you must know what you are doing

2882. param = {'max_depth': 2, 'eta': 1, 'silent': 1}

2883. watchlist = [(dtest, 'eval'), (dtrain, 'train')]

2884. num_round = 2

2885.  

2886. # user define objective function, given prediction, return gradient and second order gradient

2887. # this is log likelihood loss

2888. def logregobj(preds, dtrain):

2889. labels = dtrain.get_label()

2890. preds = 1.0 / (1.0 + np.exp(-preds))

2891. grad = preds - labels

2892. hess = preds * (1.0-preds)

2893. return grad, hess

2894.  

2895. # user defined evaluation function, return a pair metric_name, result

2896. # NOTE: when you do customized loss function, the default prediction value is margin

2897. # this may make builtin evaluation metric not function properly

2898. # for example, we are doing logistic loss, the prediction is score before logistic transformation

2899. # the builtin evaluation error assumes input is after logistic transformation

2900. # Take this in mind when you use the customization, and maybe you need write customized evaluation function

2901. def evalerror(preds, dtrain):

2902. labels = dtrain.get_label()

2903. # return a pair metric_name, result

2904. # since preds are margin(before logistic transformation, cutoff at 0)

2905. return 'error', float(sum(labels != (preds > 0.0))) / len(labels)

2906.  

2907. # training with customized objective, we can also do step by step training

2908. # simply look at xgboost.py's implementation of train

2909. bst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)

2910.  

2911. 5.Xgboost调参

      由于xgboost的参数过多，这里介绍三种思路

      （1）GridSearch

      （2）Hyperopt

      （3）老外写的一篇文章，操作性比较强，推荐学习一下。地址

      6.工程实现优化

      （1）Column Blocks and Parallelization

       

      

       

      （2）Cache Aware Access

2912. A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
2913. The main thread accumulates gradients statistics in the continuous buff er.

2914. （3）System Tricks

2915. Block pre-fetching.
2916. Utilize multiple disks to parallelize disk operations.
2917. LZ4 compression(popular recent years for outstanding performance).
2918. Unrolling loops.
2919. OpenMP

2920. 7.代码走读

      这块非常感谢杨军老师的无私奉献【4】

      个人看代码用的是SourceInsight，由于xgboost有些文件是cc后缀名，可以通过以下命令修改下（默认的识别不了）

      find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp

2921. 1
2922. 1

2923. 实际上，对XGBoost的源码进行走读分析之后，能够看到下面的主流程：

2924. cli_main.cc:

2925. main()

2926. -> CLIRunTask()

2927. -> CLITrain()

2928. -> DMatrix::Load()

2929. -> learner = Learner::Create()

2930. -> learner->Configure()

2931. -> learner->InitModel()

2932. -> for (i = 0; i < param.num_round; ++i)

2933. -> learner->UpdateOneIter()

2934. -> learner->Save()

2935. learner.cc:

2936. Create()

2937. -> new LearnerImpl()

2938. Configure()

2939. InitModel()

2940. -> LazyInitModel()

2941. -> obj_ = ObjFunction::Create()

2942. -> objective.cc

2943. Create()

2944. -> SoftmaxMultiClassObj(multiclass_obj.cc)/

2945. LambdaRankObj(rank_obj.cc)/

2946. RegLossObj(regression_obj.cc)/

2947. PoissonRegression(regression_obj.cc)

2948. -> gbm_ = GradientBooster::Create()

2949. -> gbm.cc

2950. Create()

2951. -> GBTree(gbtree.cc)/

2952. GBLinear(gblinear.cc)

2953. -> obj_->Configure()

2954. -> gbm_->Configure()

2955. UpdateOneIter()

2956. -> PredictRaw()

2957. -> obj_->GetGradient()

2958. -> gbm_->DoBoost()

2959.  

2960. gbtree.cc:

2961. Configure()

2962. -> for (up in updaters)

2963. -> up->Init()

2964. DoBoost()

2965. -> BoostNewTrees()

2966. -> new_tree = new RegTree()

2967. -> for (up in updaters)

2968. -> up->Update(new_tree)

2969.  

2970. tree_updater.cc:

2971. Create()

2972. -> ColMaker/DistColMaker(updater_colmaker.cc)/

2973. SketchMaker(updater_skmaker.cc)/

2974. TreeRefresher(updater_refresh.cc)/

2975. TreePruner(updater_prune.cc)/

2976. HistMaker/CQHistMaker/

2977. GlobalProposalHistMaker/

2978. QuantileHistMaker(updater_histmaker.cc)/

2979. TreeSyncher(updater_sync.cc)

2980.  

2981. 从上面的代码主流程可以看到，在XGBoost的实现中，对算法进行了模块化的拆解，几个重要的部分分别是：

      I. ObjFunction：对应于不同的Loss Function，可以完成一阶和二阶导数的计算。 
      II. GradientBooster：用于管理Boost方法生成的Model，注意，这里的Booster Model既可以对应于线性Booster Model，也可以对应于Tree Booster Model。 
      III. Updater：用于建树，根据具体的建树策略不同，也会有多种Updater。比如，在XGBoost里为了性能优化，既提供了单机多线程并行加速，也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现，按并行策略的不同，包括： 
      　　I). inter-feature exact parallelism （特征级精确并行） 
      　　II). inter-feature approximate parallelism（特征级近似并行，基于特征分bin计算，减少了枚举所有特征分裂点的开销） 
      　　III). intra-feature parallelism （特征内并行）

      此外，为了避免overfit，还提供了一个用于对树进行剪枝的updater(TreePruner)，以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher)，这样设计，关于建树的主要操作都可以通过Updater链的方式串接起来，比较一致干净，算是Decorator设计模式[4]的一种应用。

      XGBoost的实现中，最重要的就是建树环节，而建树对应的代码中，最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

      以ColMaker（单机版的inter-feature parallelism，实现了精确建树的策略）为例，其建树操作大致如下：

2982. updater_colmaker.cc:

2983. ColMaker::Update()

2984. -> Builder builder;

2985. -> builder.Update()

2986. -> InitData()

2987. -> InitNewNode() // 为可用于split的树结点（即叶子结点，初始情况下只有一个

2988. // 叶结点，也就是根结点) 计算统计量，包括gain/weight等

2989. -> for (depth = 0; depth < 树的最大深度; ++depth)

2990. -> FindSplit()

2991. -> for (each feature) // 通过OpenMP获取

2992. // inter-feature parallelism

2993. -> UpdateSolution()

2994. -> EnumerateSplit() // 每个执行线程处理一个特征，

2995. // 选出每个特征的

2996. // 最优split point

2997. -> ParallelFindSplit()

2998. // 多个执行线程同时处理一个特征，选出该特征

2999. //的最优split point;

3000. // 在每个线程里汇总各个线程内分配到的数据样

3001. //本的统计量(grad/hess);

3002. // aggregate所有线程的样本统计(grad/hess)，

3003. //计算出每个线程分配到的样本集合的边界特征值作为

3004. //split point的最优分割点;

3005. // 在每个线程分配到的样本集合对应的特征值集合进

3006. //行枚举作为split point，选出最优分割点

3007. -> SyncBestSolution()

3008. // 上面的UpdateSolution()/ParallelFindSplit()

3009. //会为所有待扩展分割的叶结点找到特征维度的最优split

3010. //point，比如对于叶结点A，OpenMP线程1会找到特征F1

3011. //的最优split point，OpenMP线程2会找到特征F2的最

3012. //优split point，所以需要进行全局sync，找到叶结点A

3013. //的最优split point。

3014. -> 为需要进行分割的叶结点创建孩子结点

3015. -> ResetPosition()

3016. //根据上一步的分割动作，更新样本到树结点的映射关系

3017. // Missing Value(i.e. default)和非Missing Value(i.e.

3018. //non-default)分别处理

3019. -> UpdateQueueExpand()

3020. // 将待扩展分割的叶子结点用于替换qexpand_，作为下一轮split的

3021. //起始基础

3022. -> InitNewNode() // 为可用于split的树结点计算统计量

3023.  

3024. 8.python、R对于xgboost的简单使用

      任务：二分类，存在样本不均衡问题（scale_pos_weight可以一定程度上解读此问题）

      【Python】

       

      

       

      【R】

       

      

       

      9.xgboost中比较重要的参数介绍

      （1）objective [ default=reg:linear ] 定义学习任务及相应的学习目标，可选的目标函数如下：

3025. “reg:linear” –线性回归。
3026. “reg:logistic” –逻辑回归。
3027. “binary:logistic” –二分类的逻辑回归问题，输出为概率。
3028. “binary:logitraw” –二分类的逻辑回归问题，输出的结果为wTx。
3029. “count:poisson” –计数问题的poisson回归，输出结果为poisson分布。 在poisson回归中，max_delta_step的缺省值为0.7。(used to safeguard optimization)
3030. “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题，同时需要设置参数num_class（类别个数）
3031. “multi:softprob” –和softmax一样，但是输出的是ndata * nclass的向量，可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
3032. “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

3033. （2）’eval_metric’ The choices are listed below，评估指标:

3034. “rmse”: root mean square error
3035. “logloss”: negative log-likelihood
3036. “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
3037. “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
3038. “mlogloss”: Multiclass logloss
3039. “auc”: Area under the curve for ranking evaluation.
3040. “ndcg”:Normalized Discounted Cumulative Gain
3041. “map”:Mean average precision
3042. “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
3043. “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

3044. （3）lambda [default=0] L2 正则的惩罚系数

      （4）alpha [default=0] L1 正则的惩罚系数

      （5）lambda_bias 在偏置上的L2正则。缺省值为0（在L1上没有偏置项的正则，因为L1时偏置不重要）

      （6）eta [default=0.3] 
      为了防止过拟合，更新过程中用到的收缩步长。在每次提升计算之后，算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
      取值范围为：[0,1]

      （7）max_depth [default=6] 数的最大深度。缺省值为6 ，取值范围为：[1,∞]

      （8）min_child_weight [default=1] 
      孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中，这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
      取值范围为: [0,∞]

      10.DART

      核心思想就是将dropout引入XGBoost

      示例代码

3045. import xgboost as xgb

3046. # read in data

3047. dtrain = xgb.DMatrix('demo/data/agaricus.txt.train')

3048. dtest = xgb.DMatrix('demo/data/agaricus.txt.test')

3049. # specify parameters via map

3050. param = {'booster': 'dart',

3051. 'max_depth': 5, 'learning_rate': 0.1,

3052. 'objective': 'binary:logistic', 'silent': True,

3053. 'sample_type': 'uniform',

3054. 'normalize_type': 'tree',

3055. 'rate_drop': 0.1,

3056. 'skip_drop': 0.5}

3057. num_round = 50

3058. bst = xgb.train(param, dtrain, num_round)

3059. # make prediction

3060. # ntree_limit must not be 0

3061. preds = bst.predict(dtest, ntree_limit=num_round)

3062.  

3063. 更多细节可以阅读参考文献5

      11.csr_matrix训练XGBoost

      当数据规模比较大、较多列比较稀疏时，可以使用csr_matrix训练XGBoost模型，从而节约内存。

      下面是Kaggle比赛中TalkingData开源的代码，可以学习一下，详见参考文献6。

3064. import pandas as pd

3065. import numpy as np

3066. import seaborn as sns

3067. import matplotlib.pyplot as plt

3068. import os

3069. from sklearn.preprocessing import LabelEncoder

3070. from scipy.sparse import csr_matrix, hstack

3071. import xgboost as xgb

3072. from sklearn.cross_validation import StratifiedKFold

3073. from sklearn.metrics import log_loss

3074.  

3075. datadir = '../input'

3076. gatrain = pd.read_csv(os.path.join(datadir,'gender_age_train.csv'),

3077. index_col='device_id')

3078. gatest = pd.read_csv(os.path.join(datadir,'gender_age_test.csv'),

3079. index_col = 'device_id')

3080. phone = pd.read_csv(os.path.join(datadir,'phone_brand_device_model.csv'))

3081. # Get rid of duplicate device ids in phone

3082. phone = phone.drop_duplicates('device_id',keep='first').set_index('device_id')

3083. events = pd.read_csv(os.path.join(datadir,'events.csv'),

3084. parse_dates=['timestamp'], index_col='event_id')

3085. appevents = pd.read_csv(os.path.join(datadir,'app_events.csv'),

3086. usecols=['event_id','app_id','is_active'],

3087. dtype={'is_active':bool})

3088. applabels = pd.read_csv(os.path.join(datadir,'app_labels.csv'))

3089.  

3090. gatrain['trainrow'] = np.arange(gatrain.shape[0])

3091. gatest['testrow'] = np.arange(gatest.shape[0])

3092.  

3093. brandencoder = LabelEncoder().fit(phone.phone_brand)

3094. phone['brand'] = brandencoder.transform(phone['phone_brand'])

3095. gatrain['brand'] = phone['brand']

3096. gatest['brand'] = phone['brand']

3097. Xtr_brand = csr_matrix((np.ones(gatrain.shape[0]),

3098. (gatrain.trainrow, gatrain.brand)))

3099. Xte_brand = csr_matrix((np.ones(gatest.shape[0]),

3100. (gatest.testrow, gatest.brand)))

3101. print('Brand features: train shape {}, test shape {}'.format(Xtr_brand.shape, Xte_brand.shape))

3102.  

3103. m = phone.phone_brand.str.cat(phone.device_model)

3104. modelencoder = LabelEncoder().fit(m)

3105. phone['model'] = modelencoder.transform(m)

3106. gatrain['model'] = phone['model']

3107. gatest['model'] = phone['model']

3108. Xtr_model = csr_matrix((np.ones(gatrain.shape[0]),

3109. (gatrain.trainrow, gatrain.model)))

3110. Xte_model = csr_matrix((np.ones(gatest.shape[0]),

3111. (gatest.testrow, gatest.model)))

3112. print('Model features: train shape {}, test shape {}'.format(Xtr_model.shape, Xte_model.shape))

3113.  

3114. appencoder = LabelEncoder().fit(appevents.app_id)

3115. appevents['app'] = appencoder.transform(appevents.app_id)

3116. napps = len(appencoder.classes_)

3117. deviceapps = (appevents.merge(events[['device_id']], how='left',left_on='event_id',right_index=True)

3118. .groupby(['device_id','app'])['app'].agg(['size'])

3119. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

3120. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

3121. .reset_index())

3122.  

3123. d = deviceapps.dropna(subset=['trainrow'])

3124. Xtr_app = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.app)),

3125. shape=(gatrain.shape[0],napps))

3126. d = deviceapps.dropna(subset=['testrow'])

3127. Xte_app = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.app)),

3128. shape=(gatest.shape[0],napps))

3129. print('Apps data: train shape {}, test shape {}'.format(Xtr_app.shape, Xte_app.shape))

3130.  

3131. applabels = applabels.loc[applabels.app_id.isin(appevents.app_id.unique())]

3132. applabels['app'] = appencoder.transform(applabels.app_id)

3133. labelencoder = LabelEncoder().fit(applabels.label_id)

3134. applabels['label'] = labelencoder.transform(applabels.label_id)

3135. nlabels = len(labelencoder.classes_)

3136.  

3137. devicelabels = (deviceapps[['device_id','app']]

3138. .merge(applabels[['app','label']])

3139. .groupby(['device_id','label'])['app'].agg(['size'])

3140. .merge(gatrain[['trainrow']], how='left', left_index=True, right_index=True)

3141. .merge(gatest[['testrow']], how='left', left_index=True, right_index=True)

3142. .reset_index())

3143. devicelabels.head()

3144.  

3145. d = devicelabels.dropna(subset=['trainrow'])

3146. Xtr_label = csr_matrix((np.ones(d.shape[0]), (d.trainrow, d.label)),

3147. shape=(gatrain.shape[0],nlabels))

3148. d = devicelabels.dropna(subset=['testrow'])

3149. Xte_label = csr_matrix((np.ones(d.shape[0]), (d.testrow, d.label)),

3150. shape=(gatest.shape[0],nlabels))

3151. print('Labels data: train shape {}, test shape {}'.format(Xtr_label.shape, Xte_label.shape))

3152.  

3153. Xtrain = hstack((Xtr_brand, Xtr_model, Xtr_app, Xtr_label), format='csr')

3154. Xtest = hstack((Xte_brand, Xte_model, Xte_app, Xte_label), format='csr')

3155. print('All features: train shape {}, test shape {}'.format(Xtrain.shape, Xtest.shape))

3156.  

3157. targetencoder = LabelEncoder().fit(gatrain.group)

3158. y = targetencoder.transform(gatrain.group)

3159.  

3160. ########## XGBOOST ##########

3161.  

3162. params = {}

3163. params['booster'] = 'gblinear'

3164. params['objective'] = "multi:softprob"

3165. params['eval_metric'] = 'mlogloss'

3166. params['eta'] = 0.005

3167. params['num_class'] = 12

3168. params['lambda'] = 3

3169. params['alpha'] = 2

3170.  

3171. # Random 10% for validation

3172. kf = list(StratifiedKFold(y, n_folds=10, shuffle=True, random_state=4242))[0]

3173.  

3174. Xtr, Xte = Xtrain[kf[0], :], Xtrain[kf[1], :]

3175. ytr, yte = y[kf[0]], y[kf[1]]

3176.  

3177. print('Training set: ' + str(Xtr.shape))

3178. print('Validation set: ' + str(Xte.shape))

3179.  

3180. d_train = xgb.DMatrix(Xtr, label=ytr)

3181. d_valid = xgb.DMatrix(Xte, label=yte)

3182.  

3183. watchlist = [(d_train, 'train'), (d_valid, 'eval')]

3184.  

3185. clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

3186.  

3187. pred = clf.predict(xgb.DMatrix(Xtest))

3188.  

3189. pred = pd.DataFrame(pred, index = gatest.index, columns=targetencoder.classes_)

3190. pred.head()

3191. pred.to_csv('sparse_xgb.csv', index=True)

3192.  

3193. #params['lambda'] = 1

3194. #for alpha in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:

3195. # params['alpha'] = alpha

3196. # clf = xgb.train(params, d_train, 1000, watchlist, early_stopping_rounds=25)

3197. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3198. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3199. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3200. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3201. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3202. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码

3203. 1.序

      　　距离上一次编辑将近10个月，幸得爱可可老师（微博）推荐，访问量陡增。最近毕业论文与xgboost相关，于是重新写一下这篇文章。

      　　关于xgboost的原理网络上的资源很少，大多数还停留在应用层面，本文通过学习陈天奇博士的PPT、论文、一些网络资源，希望对xgboost原理进行深入理解。（笔者在最后的参考文献中会给出地址）

      2.xgboost vs gbdt

      　　说到xgboost，不得不说gbdt，两者都是boosting方法（如图1所示），了解gbdt可以看我这篇文章 地址。

       

       

      图1

       

      　　如果不考虑工程实现、解决问题上的一些差异，xgboost与gbdt比较大的不同就是目标函数的定义。

       

      

       

      　　注：红色箭头指向的l即为损失函数；红色方框为正则项，包括L1、L2；红色圆圈为常数项。xgboost利用泰勒展开三项，做一个近似，我们可以很清晰地看到，最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

      3.原理

      对于上面给出的目标函数，我们可以进一步化简

      （1）定义树的复杂度

      对于f的定义做一下细化，把树拆分成结构部分q和叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去，而w给定了每个索引号对应的叶子分数是什么。

       

      

       

      定义这个复杂度包含了一棵树里面节点的个数，以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式，不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

       

      

       

      注：方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ，gamma

      在这种新的定义下，我们可以把目标函数进行如下改写，其中I被定义为每个叶子上面样本集合 ,g是一阶导数，h是二阶导数

       

      

       

      这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

       

      

       

      最终公式可以化简为

       

      

       

      通过对求导等于0，可以得到

       

      

       

      然后把最优解代入得到：

       

      

       

      （2）打分函数计算示例

      Obj代表了当我们指定一个树的结构的时候，我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

      

      （3）分裂节点

      论文中给出了两种分裂节点的方法

      （1）贪心法：

      每一次尝试去对已有的叶子加入一个分割

       

      

       

      对于每次扩展，我们还是要枚举所有可能的分割方案，如何高效地枚举所有的分割呢？我假设我们要枚举所有x < a 这样的条件，对于某个特定的分割a我们要计算a左边和右边的导数和。

       

      

       

      我们可以发现对于所有的a，我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

      观察这个目标函数，大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好，因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝， 当引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。大家可以发现，当我们正式地推导目标的时候，像计算分数和剪枝这样的策略都会自然地出现，而不再是一种因为heuristic（启发式）而进行的操作了。

      下面是论文中的算法

       

      

       

      （2）近似算法：

      主要针对数据太大，不能直接进行计算

       

      

       

      4.自定义损失函数（指定grad、hess）

      （1）损失函数

      

      （2）grad、hess推导

      

      （3）官方代码