[Machine Learning] XGBoost

1. XGBoost介绍

XGBoost模型即是一些“串联”树结构的组合，最终预测结果由多棵树共同决定。

模型公式：$y=f(z)=\sum_{k=1}^K f_k(z)$
模型预测：$\hat y_i = \sum_{k=1}^K \hat f_k(x_i)$

类似模型预测，第$t$步的估计可以表示为 $\hat y_i = \sum_{k=1}^t \hat f_k(x_i)$。如何有效的估计每棵树，且如何重组多棵树使得XGBoost成为如今竞赛的常胜将军？我们需要解决以下三个问题：

问题0：树的组合方式
问题1：每棵树的结构-分枝
问题2：每个枝节点的预测值

下面分别按步骤依次解决上面三个问题：

1. 解决方案0：树的组合方式是“串联”，即每一个树是建立在前面所有树预测残差的基础上优化得到。所以后续树的结构和枝节点的预测值也将按照“并联”方式逐个估计。

2. 解决方案2：常规的解决方案是先解决树的结构，再确定每个节点的估计值，但XGBoost的推导过程是先假设树的结构的基础上先估计每个节点预测值，进而完善树的结构。任何优化问题之前，损失函数的确定是重中之重。XGBoost采用的是惩罚式损失$Obj = L+\Omega$，其中$L$为损失函数，$\Omega$为惩罚项。且惩罚函数为：$\Omega=\gamma T + \lambda\sum_{j=1}^T \omega_j^2$，其中$T$是枝节点的个数，$\omega_j$是$j$枝上的得分。假设已得到$t-1$棵树，我们来估计第$t$棵树的结构和节点估计。首先将数据和估计带入损失函数得：

   $$Obj^{(t)}= \sum_{i=1}^n \ell(y_i, \hat y_i^{(t-1)}+f_i) + \gamma T +\lambda \sum_{j=1}^T \omega_j^2$$ 将以上损失函数部分进行二次Taylor展开得：
   $$\begin{align}Obj^{(t)}&\sim\sum_{i=1}^n \Big(\ell(y_i, \hat y_i^{(t-1)})+g_i f_i + \frac{1}{2} h_i f_i^2\Big) + \gamma T +\lambda \sum_{j=1}^T \omega_j^2\nonumber\\ &\sim \sum_{i=1}^n \Big(g_i f_i + \frac{1}{2} h_i f_i^2\Big) + \gamma T +\lambda \sum_{j=1}^T \omega_j^2\nonumber\\ &\sim\sum_{j=1}^T \sum_{i \in I_j} \Big( g_i f_i + \frac{1}{2} h_i f_i^2 \Big)+ \gamma T + \lambda\sum_{j=1}^T \omega_j^2\nonumber\nonumber\\ &=\sum_{j=1}^T \Big( \big(\sum_{i \in I_j} g_i\big) \omega_j + \frac{1}{2} \big(\sum_{i \in I_j}h_i\big) \omega_j^2 \Big)+ \gamma T + \lambda\sum_{j=1}^T \omega_j^2\nonumber\\ &= \sum_{j=1}^T \Big( \big(\sum_{i \in I_j} g_i\big) \omega_j + \frac{1}{2} \big(\sum_{i \in I_j}h_i +\lambda\big) \omega_j^2 \Big)+ \gamma T\\ &= \sum_{j=1}^T \frac{1}{2} \Big(\sum_{i \in I_j}h_i +\lambda\Big) \Bigg(\omega_j^2 + 2\frac{\sum_{i \in I_j} g_i}{\sum_{i \in I_j}h_i +\lambda}\omega_j\Bigg) + \gamma T\nonumber\\ &= \sum_{j=1}^T \frac{1}{2} \Big(\sum_{i \in I_j}h_i +\lambda\Big)\Bigg(\omega_j+\frac{\sum_{i \in I_j} g_i}{\sum_{i \in I_j}h_i +\lambda}\Bigg)^2-\frac{1}{2} \sum_{j=1}^T\frac{\big(\sum_{i \in I_j} g_i\big)^2}{\sum_{i \in I_j}h_i +\lambda}+\gamma T=0\end{align}$$
   其中$g_i = \partial_{y^{(t-1)}}\ell(y_i, \hat y_i^{(t-1)})$和$h_i = \partial^2_{y^{(t-1)}}\ell(y_i, \hat y_i^{(t-1)})$，且$I_j$为$j$个节点包含的数据指标集。将上面(84)式对$\omega_j$求导得 $$\frac{d Obj^{(t)}}{d \omega_j} =\sum_{i \in I_j} g_i + \big(\sum_{i \in I_j}h_i +\lambda\big) \omega_j=0$$ 或从上式(85)的二次函数重组可以得到
   $$\omega_j^* = - \frac{\sum_{i \in I_j} g_i}{\sum_{i \in I_j}h_i +\lambda}$$ 且此时的目标函数可表达为： $$\begin{aligned}Obj^{(t)} &= -\frac{1}{2} \sum_{j=1}^T\frac{\big(\sum_{i \in I_j} g_i\big)^2}{\sum_{i \in I_j}h_i +\lambda}+ \gamma T\\ &= -\frac{1}{2}\sum_{j=1}^T \Bigg(\frac{\big(\sum_{i \in I_j} g_i\big)^2}{\sum_{i \in I_j}h_i +\lambda}+\gamma\Bigg)\\ &= -\frac{1}{2} \sum_{j=1}^T \Bigg(\frac{G}{H+\lambda}-2\gamma\Bigg)\end{aligned}$$

3. 解决方案1：如何得到树的结构呢？关键是确定每棵树的节点变量和如何划分节点。首先方案是：a. 对于每个节点，遍历全部变量 b. 对于每个变量，将所有的观测进行排序并进行扫描式，得到最好的切分点。如何判定切分点的好坏呢？根据以上增加该切分点，目标函数得变化量： $$\begin{aligned}\Delta &= Obj^{(t)} - \Bigg(Obg^{(t)}_L + Obj^{(t)}_R\Bigg)\\ &= \frac{1}{2} \Bigg(\frac{G_L}{H_L+\lambda}+ \frac{G_R}{H_R+\lambda}-\frac{G}{H+\lambda}\Bigg)-\gamma\end{aligned}$$ 达到最大。

2. R简单实现

# Load packages
library(xgboost)

# Read data from xgboost package
data('agaricus.train', package = 'xgboost')
data('agaricus.test', package = 'xgboost')

tr <- agaricus.train
te <- agaricus.test

# Investigate data set
str(tr$data)   ## dgCMatrix: Sparse matrix class from package Matrix

# Train xgboost model
## Input features: dense or sparse matrix both are ok
## Target variable: numeric vector(0-n for classification and real values for regression)
## Objective: 'reg:linear' for regression or 'binary:logistic'
## Number of iteration: number of trees added to the model


cv.res <- xgb.cv(data = tr$data, nfold = 5, label = tr$label,
               nround = 10, objective = 'binary:logistic')     # Cross Validation 
cv.res <- xgb.cv(data = tr$data, nfold = 5, label = tr$label,
              nround = 10, objective = 'binary:logistic',
              eval_metric = 'auc')                             # Cross Validation

bst <- xgboost(data = tr$data, label = tr$label,
               nround = 2, objective = 'binary:logistic',
               eval_metric = 'auc')
pred <- predict(bst, te$data)
table(pred, te$label)