前言
xgboost是有监督学习的利器，相较于与决策树、随机森林、gbdt等都有着明显优势，也许是各个算法有自己适用的领域吧。xgboost厉害之处在于同时控制了残差和复杂度，在迭代优化过程中，会同时考量这两个因素，本文会按照我自己的学习时的思路来介绍该算法，包括陈天奇(天才少年)的C++源码解读。
带着问题走
由于之前也是使用了挺多的算法，包括dtree, random forest, 还有各种pair wise、list wise算法，但是听到同事的介绍还是产生了兴趣，我本人介入机器学习领域时间较短，开始在听完xgboost的介绍时，就迫切的想了解两个东西:
xgboost如何并行化？
xgboost是如何控制复杂度的?
xgboost如何选择最优分裂?
这三个问题是我继续了解的动力，另外后面的阅读和学习中，又知道了关于xgboost的最优化split、最优化Gain等详细知识，另外还有各种objective 灵活设置，非常方便。有些材料仅供参考:
原理介绍?
基本的概念如tree, mart, residual等就不介绍了，直接进入主题，xgboost的目标函数:
$Obj(\theta) = L(\theta) + \Omega(\theta)$
$\theta$是我们要求的最优解，可以对应xgboost上各个tree上最有split和最有leaf value，$L$是残差，即模型预测值和实际训练样本的差值，当然是越小越好，$\Omega$是模型复杂度，这个也很重要，防止过拟合以及严重的抖动，一般模型需要做bias和varience的校验，一个对应$L$，一个对应$\Omega$。xgboost是mart，预测值是所有的tree的加权相加:
$\hat{y_i}=\sum_{i=1}^{K}f_k(x_i),f_k \in{\mathcal{F}}$
以上是mart建立之后如何预测，tree建立在分裂时，会有目标，Objective细化为
$Obj=\sum_{i=1}^{n} l(y_i, \hat{y_i})+\sum_{k=1}^{K} \Omega(f_i)
=\sum_{i=1}^{n}l(y_i, \hat y_i^{t-1} + f_t(x_i)) + \Omega(f_t) + constant $
$f_t$是第t轮迭代，即第t颗树，由于t轮之前的参数已经确定，所以最后归结为constant，目标也是寻找最优的$f_t$，让Obj最小。
使用二阶泰勒展开的方式
$f(x+\Delta{x}) \approx f(x) + f’(x)\Delta{x} + \frac{1}{2}f’’(x)\Delta{x^2}$
做一下转化:
$l(y_i, \hat y_i^{t-1} + f_t(x_i)) => f(x+\Delta{x}) $
$g_i=\delta_{\hat y_i^{t-1}}{l(y_i, \hat{y}^{(t-1)})}$
$h_i=\delta_{\hat y_i^{t-1}}^2 {l(y_i, \hat{y}^{(t-1)})}$
则Obj在t轮迭代的目标为
$Obj=\sum_{i=1}^{n}l(y_i, \hat y_i^{t-1} + f_t(x_i)) + \Omega(f_t) + constant
=\sum_{i=1}^{n}[l(y_i, \hat y_i^{t-1}) + g_if_t(x_i) + \frac{1}{2}h_if_t^2(x_i)] + \Omega(f_t) + constant$
模型复杂度$\Omega(f_t)$的定义：
$\Omega(f_t)=\gamma{T} + \frac{1}{2}\lambda\sum_{j=1}^{T}w_j{^2}$
$T$是叶子节点的个数，$w_j$是叶子节点对叶的leaf_score,进一步化简$Obj(f_t)$:
$Obj=\sum_{i=1}^{n}[g_if_t(x_i) + \frac{1}{2}h_if_t^2(x_i)] + \Omega(f_t)
=\sum_{i=1}^{n}[g_iw_{q(x_i)} + \frac{1}{2}h_iw^2_{q(x_i)}] + \gamma{T} + \frac{1}{2}\lambda\sum_{j=1}^{T}w_j{^2}$
$=\sum_{j=1}^{T}[(\sum_{i \in I_j}g_i)w_j + \frac{1}{2}(\sum_{i \in I_j}h_i + \lambda)w^2_j] + \gamma T$
这样终于看出GBDT计算收益的最终公式了：
$(1) w_j = - \frac{G_j}{H_j+\lambda}$
$(2) Obj =- \frac{1}{2}\sum_{j=1}^{T}\frac{G^2_j}{H_j + \lambda} + \gamma T$
$(3) Gain = \frac{1}{2}[\frac{G^2_L}{H_L+\lambda} + \frac{G^2_R}{H_R+\lambda} - \frac{(G_L + G_R)^2}{H_L + H_R + \lambda}] - \gamma$
$\gamma$就是传说中的一阶正则项，$\gamma$是后面提到的参数min_split_loss，剪枝使用，太小的收益不必分裂。树节点分裂的收益达到某个值后再分裂，而二阶正则项$\lambda$主要是为了防止过拟合, 增强模型的稳定性。
如何并行化？
xgboost相对于gbdt是有一个并行的优势，怎么并行的呢，下面会读数据并行化和树分裂时的并行化分两部分介绍
读数据并行
开始的输入数据是libsvm格式，形如”1 1:1.2 2:2.1 3:3.5 …”一行一行的数据，训练集合可能会非常大，大家都知道xgboost会枚举所有的特征值，选用最佳分割，所以开始是需要按列把各个feature排序好，这一步在xgboost里是并行做的。。。xgboost使用DMatrix存储训练数据，每次训练前会调用InitColAccess()将行训练数据支持按列取的数据，这里是通过MakeOneBatch多线程来实现的，这么多列属性的排序被n个线程平分，最后将数据转化成 key,value的方式，key是列id，value是rowid, val的pair，rowid是该列所对应的行编号，用于拿出整行训练数据的信息，val是该行该列下的具体数字，这个信息最后会存入DMatrix的SparsePage中，可知道最后是存了feature_num个实际train_set，用于训练阶段方便的取数据，不用做重复的排序工作了。
分裂选择并行化
tree进行分裂的时候会考量每个feature的每个value，如果串行计算，肯定会比较慢，xgboost使用了并行化，多个线程综合考虑每个列的情况，主要在FindSplit()里做的
● this->UpdateSolution(iter->Value(), gpair, *p_fmat); 主函数，枚举每一个feature，的每一个value，尝试收益
● this->SyncBestSolution(qexpand); 多线程merger，当前node的最好的split
● sync solution; 同步成果，给tree建node，左右儿子等
每个线程会枚举自己的feature的各个split value，将最佳分裂保存到当前线程的结构中，最后所有的线程会根据各自的best split做merge，从而找到当前tree的最佳分裂，这里也用到了并行化，共享了读数据阶段产生的SparsePage
这些并行化和共享让xgboost的训练速度非常快。以上基本上解决了我开始的三个疑问，下面说一些xgboost的细节。
细节1: tree如何维护当前的数据集？
这个是通过postion这个vector来实现的，vector是训练机的大小，每个值代表当前训练样本属于的树节点编号，开始都是0，即开始所有的样本都是在根节点上。expand这个vector保存当前叶子节点，即即将开始下一轮分裂的叶子节点编号，当然开始也是初始化为0，postion里的所有值都是取自于expand
细节2: 模型复杂度
主要说下$\gamma,\alpha, \lambda$这三个，这三个控制这模型复杂度。
- $\gamma$是
min_split_loss，控制树的分裂深度，收益必须大于该值才分裂。 $\alpha$是一阶正则项，也是控制复杂度的，不过比较弱，一般训练设置为0，控制这Grad(一阶梯度)不至于过大或者过小
123456if (p.reg_alpha == 0.0f) {return Sqr(sum_grad) / (sum_hess + p.reg_lambda);} else {return Sqr(ThresholdL1(sum_grad, p.reg_alpha)) /(sum_hess + p.reg_lambda);}$\lambda$是二阶正则项，具体使用上面的推到也已经明了了
Lambda Object
非常重要的应用，相对于gbdt，或许这个才算是最重要的区别，即可以自定义目标函数(loss)，传统的Objective loss即预测值和标签值的差异，在实际工程中的应用显然是不够的，很多指标有实际的公式，如果表示成残差，梯度，这个活就交给lambda来做了，如下要素:
● 组的概念
● 选择两个组内两个标签不同的一对
● $P_{ij} = \frac{1}{1+e^{-\sigma(s_i-s_j)}}$
● $C = -\bar P_{ij} log P_{ij} - (1 - \bar P_{ij})log(1 - P_{ij})$
● $C=log(1+e^{-\sigma(s_i-s_j)})$
● $\delta{w_k} = -\eta\sum_i\lambda_{i}\frac{\partial s_i}{\partial w_k}$
● $\lambda_i = \sum_{j:{i, j}\in I} \lambda_{ij} - \sum_{j:{j, i}\in I} \lambda_{ij}$
● $\lambda_{ij} = \frac{\partial C(s_i - s_j)}{\partial{s_i}} = \frac{-\sigma}{1+e^{\sigma(s_i-s_j)}} * \Delta$
这样，成功的讲工程上的知道公式产生的残差计算成功的引入到了梯度中去,$\Delta$可以是NDCG以及MAP等等。如何从pointwise的思路转到pairwise listwise，$\lambda$的应用是关键，其中的转化思路如上面的推理所示，具体代码在src/objective/rank_obj.cc中，定义了基本的类LambdaRankObj，在此基础上实现了PairwiseRankObj和LambdaRankObjNDCG以及LambdaRankObjMAP等。LambdaRankObj主要的函数就是
该函数的流程如下
- 遍历group
- 组内排序，按预测值放入lst, 按照标签值放入rec
- rec内按照不同的标签组成pairs, 采用uniform_int_distribution的方式组
- 调用GetLambdaWeight获取$\Delta$, pairwise的什么也没干，默认是1
- 根据预测值和标签，算Sigmod值，然后乘$\Delta$, 算出的g,h分别放入 out_gpair 的 i, j中，i > j, 一个是加，一个是减1234gpair[pos.rindex].grad += g * w;gpair[pos.rindex].hess += 2.0f * w * h;gpair[neg.rindex].grad -= g * w;gpair[neg.rindex].hess += 2.0f * w * h;
NDCG的GetLambdaWeight怎么算？这里就引入了计算公式。
计算流程：
- 按照labels，算出最大的maxNDCG作为分母
- 算original ndcg以及交换之后的ndcg，做差值，处maxNDCG
问题1:如何做到好的结果按照NDCG更快排上去的?
好的结果产生的$\Delta$在NDCG计算时会更大，weight会更高，所以获得加权会更高
问题2:好的结果每次都会加权，最后预测值会很大？
不会的, 如果好的记过得到的score很高了，即$s_i$很高，会拉开与其他结果($s_j$)的差距，所以系数$\frac{-\sigma}{1+e^{\sigma(s_i-s_j)}}$会很小，这个系数表示目前模型的好坏，如果很好了即使$\Delta$很大也是没有多少梯度的
应用参数
首先推荐官方的API，这个没有更权威的了，如果使用，一定要仔细阅读，我是在python中使用，所以会看python的API)。首先说一下命令行版本也就是github拉下来之后直接编译出的二进制，直接使用”./xgboost config.ini”的方式即可，注意配置的格式，其实参考”src/cli_main.cc”也可以大概看出需要配置那些参数，我贴一下自己的训练配置和预测配置:
源码中是使用dmlc-core来做的参数处理，十分方便，我是做搜索排序相关，所以目标使用ndcg, reg_lambda对应上面所说的$\lambda$，min_split_loss是树分裂所达到的最小收益值，即上面公式(3)中的Gain，一阶正则项$\alpha＝0$。
不过还是推荐使用python版本，python版本在训练时，可以看到每轮之后test集的准确度(ndcg值)，可以直观的看到效果，主要是param的设置，如下:
最后贴一张训练截图
