集成學習之Boosting —— XGBoost

集成學習之Boosting —— AdaBoost

集成學習之Boosting —— Gradient Boosting

集成學習之Boosting —— XGBoost

Gradient Boosting 可以看做是一個總體的算法框架，起始於Friedman 的論文 [Greedy Function Approximation: A Gradient Boosting Machine] 。XGBoost (eXtreme Gradient Boosting) 是於2015年提出的一個新的 Gradient Boosting 實現，由華盛頓大學的 陳天奇 等人開發，在速度和精度上都有顯著提升，因而近年來在 Kaggle 等各大數據科學比賽中都得到了廣泛應用。本文主要對其原理進行闡述，並將其與傳統的 GBDT 進行比較。

大體來看，XGBoost 在原理方面的改進主要就是在損失函數上作文章。一是在原損失函數的基礎上添加了正則化項產生了新的目標函數，這類似於對每棵樹進行了剪枝並限制了葉結點上的分數來防止過擬合。二是對目標函數進行二階泰勒展開，以類似牛頓法的方式來進行優化（事實上早在 [Friedman, J., Hastie, T. and Tibshirani, R., 1999] 中就已有類似方案，即利用二階導信息來最小化目標函數，陳天奇在論文中也提到了這一點）。

在上一篇文章中，了解到 Gradient Boosting 的思想可以理解為通過函數空間的梯度下降來最小化目標函數 \(L(f) = \sum\limits_{i=1}^NL(y_i,f_m(x_i))\)，其只使用了一階導信息，而 XGBoost 引入二階導的一大好處是可以推導出一種新的增益計算方法，事實證明采用新的增益計算方法在優化目標函數上更加有效，精確度上也勝過傳統的 GBDT。所以下面也從目標函數的優化入手進行推導。

XGBoost的目標函數

與 GBDT 一樣，XGBoost 同樣采用加法模型，設基學習器為\(f(x)\)，預測值為\(\hat{y}\)：

\[\hat{y} = \sum\limits_{m=1}^M f_m(x)\;, \quad f_m \in \mathcal{F} \tag{1.1} \]

在第m步，前m-1個基學習器是固定的，因而目標函數為

\[\mathcal{L}^{(m)} = \sum\limits_{i=1}^N L(y_i, \;\hat{y_i}^{(m-1)} + f_m(x_i)) + \Omega(f_m) \tag{1.2} \]

在傳統的 GBDT 中，是沒有正則化項 \(\Omega\) 的，而在XGBoost中采用的是

\[\Omega(f) = \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \]

其中 \(J\) 為葉結點數目，\(b_j\) 為各個葉結點上的值。該項限制了樹的復雜度，這樣相當於使葉結點的數目變小，同時限制葉結點上的分數，因為通常分數越大學得越快，就越容易過擬合。

接下來將 \(\Omega\) 代入\((1.2)\) 式並對目標函數在 \(\hat{y}_i^{(m-1)}\) 處進行二階泰勒展開：

\[\mathcal{L}^{(m)} \simeq \sum\limits_{i=1}^N \left[ L(y_i, \hat{y}_i^{(m-1)}) + g_if_m(x_i) + \frac12 h_i f_m^2(x_i)\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \tag{1.3} \]

其中 \(g_i = \frac{\partial L(y_i, \,\hat{y}_i^{m-1})}{\partial \, \hat{y}^{(m-1)}}\;\;,\;\; h_i = \frac{\partial^2 L(y_i, \,\hat{y}_i^{(m-1)})}{(\partial\, \hat{y}^{(m-1)})^2}\)

在上一篇 Gradient Boosting 中提到決策樹將特征空間划分為各個獨立區域，每個樣本只屬於其中一個區域，因而單顆決策樹可表示為 \(f(x) = \sum\limits_{j=1}^J b_j I(x \in R_j)\)，如果將所有樣本加起來，則可表示為 \(\sum\limits_{i=1}^N f(x_i) = \sum\limits_{j=1}^J \sum\limits _{x \in R_j} b_j\) ，代入 \((1.3)\) 式並將常數項 \(L(y_i, \hat{y}_i^{(m-1)})\) 移去：

\[\begin{align*} \mathcal{L}^{(m)} &= \sum\limits_{i=1}^N \left[ g_if_m(x_i) + \frac12 h_i f_m^2(x_i)\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \\ &= \sum\limits_{j=1}^J \left[ \sum\limits_{x_i \in R_j}g_ib_j + \frac12 \sum\limits_{x_i \in R_j} h_i b_j^2\right] + \gamma \, J + \frac12 \lambda \sum\limits_{j=1}^J b_j^2 \\ &= \sum\limits_{j=1}^J \left[ \sum\limits_{x_i \in R_j}g_ib_j + \frac12 (\sum\limits_{x_i \in R_j} h_i + \lambda) \,b_j^2 \right] + \gamma \, J \\ & = \sum\limits_{j=1}^J(G_j b_j + \frac12 (H_j + \lambda) \,b_j^2) + \gamma\, J \tag{1.4} \end{align*} \]

其中 \(G_j = \sum\limits_{x_i \in R_j} g_i \;\; , \;\; H_j = \sum\limits_{x_i \in R_j} h_i\)

XGBoost的增益計算和樹分裂方法

經過一系列推導后，現在我們的目標變為最小化 \((1.4)\) 式，並求得相應的樹結構 \(\{R_j\}^J_1\) 和葉結點上的值 \(\{b_j\}^J_1\) 。

精確地划分 \(\{R_j\}^J_1\) 是一個 NP hard 問題，現實中常使用貪心法，遍歷每個特征的每個取值，計算分裂前后的增益，並選擇增益最大的進行分裂。


對於具體增益的衡量標准，在幾種決策樹算法中，ID3 采用了信息增益：

\[I(D,A) = H(D) - H(D|A) = H(D) - \sum\limits_{v=1}^V\frac{|D^v|}{|D|} H(D^v) \]

其中 V 表示特征 A 有 V 個可能的取值，\(D^v\) 表示第 v 個取值上的樣本數量。

C4.5 采用了信息增益比：

\[I_R(D,A) = \frac{I(D, A)}{H_A(D)} \]

其中 \(H_A(D) = -\sum\limits_{v=1}^V \frac{|D^v|}{|D|} log_2 \frac{|D^v|}{|D|}\) 。

CART 分類樹采用了基尼系數：

\[Gini(D,A) = \sum\limits_{v=1}^V \frac{|D^v|}{|D|} Gini(D^v) = \sum\limits_{v=1}^V \frac{|D^v|}{|D|} \left(1 - \sum\limits_{k=1}^K \left\lbrace\frac{|C_k|}{|D|}\right\rbrace^2\right) \]

其中 K 為類別個數，\(|C_k|\) 為 \(D\) 中屬於第k類的樣本數量。

CART 回歸樹采用了均方誤差：

\[\mathop{min}_{A,s}\Bigg[\mathop{min}_{c_1}\sum\limits_{x_i \in R_1(A,s)}(y_i - c_1)^2 + \mathop{min}_{c_2}\sum\limits_{x_i \in R_2(A,s)}(y_i - c_2)^2\Bigg] \]

其中 \(c_1\)為特征 A 的切分點 s 划分出來的 \(R_1\) 區域的樣本輸出均值，\(c_1 = mean(y_i | x_i \in R_1(A,s))\)，\(c_2\)為 \(R_2\) 區域的樣本輸出均值，\(c_2 = mean(y_i | x_i \in R_2(A,s))\) 。

而XGBoost則提出了一種新的增益計算方法。

如果已經確定了樹的結構 \(\{R_j\}^J_1\) ，則直接對 \((1.4)\) 式求導，最優值為：

\[b^*_j = - \frac{G_j}{H_j + \lambda} \]

再代入\((1.4)\)式得到此時最小的為目標函數為：

\[\mathcal{L^{(m)}} = -\frac12 \sum\limits_{j=1}^J \frac{G_j^2}{H_j+ \lambda} + \gamma\, J \tag{1.5} \]

式 \((1.5)\) 可以認為是 XGBoost 的打分函數，該值越小，說明樹的結構越好，下圖示例了該式的計算方法：

該式的優點是除了能作為樹分裂的衡量標准外，還能使 XGboost 適用於各種不同的損失函數，所以 XGBoost 包中支持自定義損失函數，但前提是一階和二階可導。

從另一角度看， \((1.5)\) 式就類似於傳統決策樹中的不純度指標，在決策樹中我們希望一次分裂后兩邊子樹的不純度越低越好，對應到XGBoost中則是希望 \((1.5)\)式 越小越好，即 \(\frac{G_j^2}{H_j+ \lambda}\) 越大越好，這樣分裂前后的增益定義為：

\[Gain = \frac{G_L^2}{H_L+ \lambda} + \frac{G_R^2}{H_R+ \lambda} - \frac{(G_L + G_R)^2}{H_L+ H_R + \lambda} - \gamma \tag{1.6} \]

\(\frac{G_L^2}{H_L+ \lambda}\) 為分裂后左子樹的分數，\(\frac{G_R^2}{H_R+ \lambda}\) 為分裂后右子樹的分數，$\frac{(G_L + G_R)^2}{H_L+ H_R + \lambda} $ 為分裂前的分數，\(Gain\) 越大說明越值得分裂。當然 \((1.6)\) 式中 \(Gain\) 可能會變成負的，這個時候分裂后的目標函數不會減少，但這並不意味着不會分裂 。事實上 XGBoost 采用的是后剪枝的策略，建每棵樹的時候會一直分裂到指定的最大深度(max_depth)，然后遞歸地從葉結點向上進行剪枝，對之前每個分裂進行考察，如果該分裂之后的 \(Gain \leqslant 0\)，則咔咔掉。 \(\gamma\) 是一個超參數，具有雙重含義，一個是在 \((1.3)\) 式中對葉結點數目進行控制的參數；另一個是 \((1.6)\) 式中分裂前后 \(Gain\) 增大的閾值，當然二者的目的是一樣的，即限制樹的規模防止過擬合。

接下來考察決策樹建立的過程。如果是使用貪心法，就是遍歷一個葉結點上的所有候選特征和取值，分別計算 \(Gain\) ，選擇 \(Gain\) 最大的候選特征和取值進行分裂，如下樹分裂算法流程 (注意這是單個葉結點的分裂流程)：

有了上述單個葉結點上的分裂流程后，我們可以總結下整個 XGBoost 的學習算法：

1. 初始化 \(f_0(x)\)
2. for m=1 to M :
   (a) 計算損失函數在每個訓練樣本點的一階導數 \(g_i = \frac{\partial L(y_i, \,\hat{y}_i^{m-1})}{\partial \, \hat{y}^{(m-1)}}\) 和二階導數 \(h_i = \frac{\partial^2 L(y_i, \,\hat{y}_i^{(m-1)})}{(\partial\, \hat{y}^{(m-1)})^2}\) , $ i = 1,2 \cdots N$
   (b) 遞歸地運用樹分裂算法生成一顆決策樹 \(f_m(x)\)
   (c) 把新生成的決策樹添加到模型中， $\hat{y_i}^{(m)} = \hat{y_i}^{(m-1)} + f_m(x) $

如果把上述 XGBoost 的學習算法和上一篇中傳統 GBDT 的學習算法作比較，XGBoost 的主要優勢是在損失函數中加入正則化項后使得學習出來的樹更加簡單，防止過擬合，但除此以外並不能體現出 XGBoost 的速度優勢。XGBoost 之所以快的一大原因是在工程上實現了 Column Block 方法，使得並行訓練成為了可能。

Column Block for Parallel Learning

對於決策樹來說，對連續值特征進行划分通常比較困難，因為連續值特征往往取值眾多。通常的做法是先按特征取值對樣本排序，再按順序計算增益選擇划分點。若每次分裂都要排一次序，那是非常耗時的，所以 XGBoost 的做法是在訓練之前，預先按特征取值對樣本進行了排序，然后保存為block結構，采用CSC格式存儲，每一列(一個特征列)均升序存放，這樣，一次讀入數據並排好序后，以后均可重復使用，大大減小計算量。

由於已經預先排過序，所以在樹分裂算法中，通過一次線性掃描 (linear scan) 就能找出一個特征的最優分裂點，如下圖所示，同時可以看到缺失值並沒有保存：

陳天奇論文里有關鍵的一句話：Collecting statistics for each columns can be parallelized 。由於已經預先保存為block 結構，所以在對葉結點進行分裂時，每個特征的增益計算就可以開多線程進行，訓練速度也由此提升了很多。而且這種 block 結構也支持列抽樣，只要每次從所有 block 特征中選擇一個子集作為候選分裂特征就可以了，據我的使用經驗，列抽樣大部分時候都比行抽樣的效果好。

近似分裂算法

當數據量非常大難以被全部加載進內存時或者在分布式環境下時，上文的樹分裂算法將不再適用，因而作者提出了近似分裂算法，不僅解決了這個問題，也同時能提升訓練速度，具體流程見下圖：

近似分裂算法其實就是對連續性特征進行離散化，對於某個特征 \(k\)，算法首先根據特征分布找到 \(l\) 個分位點的候選集合 \(S_k = \{s_{k1}, s_{k2}, ... ,s_{kl} \}\) ，然后根據集合 \(S_k\) 中的分位點將相應樣本划分到桶(bucket)中。在遍歷該特征時，只需遍歷各個分位點，對每個桶內的樣本統計值 \(g\)、\(h\) 進行累加統計，尋找最佳分裂點進行分裂。該算法可分為 global 近似和 local 近似，global 就是在新生成一棵樹之前就對各個特征計算分位點並划分樣本，之后在每次分裂過程中都重復利用這些分位點進行划分，而 local 就是每一次結點分裂時都要重新計算分位點。

總結

最后總結一下 XGBoost 與傳統 GBDT 的不同之處：

• 傳統 GBDT 在優化時只用到一階導數信息，XGBoost 則對目標函數進行了二階泰勒展開，同時用到了一階和二階導數。另外 XGBoost 工具支持自定義損失函數，只要函數可一階和二階求導。

• XGBoost 在損失函數中加入了正則化項，用於控制模型的復雜度，防止過擬合，從而提高模型的泛化能力。

• 傳統 GBDT 采用的是均方誤差作為內部分裂的增益計算指標（因為用的都是回歸樹），而 XGBoost 使用的是經過優化推導后的式子，即式 \((1.6)\) 。

• XGBoost 借鑒了隨機森林的做法，支持列抽樣，不僅能降低過擬合，還能減少計算量，這也是 XGBoost 異於傳統 GBDT 的一個特性。

• XGBoost 添加了對稀疏數據的支持，在計算分裂增益時不會考慮帶有缺失值的樣本，這樣就減少了時間開銷。在分裂點確定了之后，將帶有缺失值的樣本分別放在左子樹和右子樹，比較兩者分裂增益，選擇增益較大的那一邊作為默認分裂方向。

• 並行化處理：由於 Boosting 本身的特性，無法像隨機森林那樣樹與樹之間的並行化。XGBoost 的並行主要體現在特征粒度上，在對結點進行分裂時，由於已預先對特征排序並保存為block 結構，每個特征的增益計算就可以開多線程進行，極大提升了訓練速度。

• 傳統 GBDT 在損失不再減少時會停止分裂，這是一種預剪枝的貪心策略，容易欠擬合。XGBoost采用的是后剪枝的策略，先分裂到指定的最大深度 (max_depth) 再進行剪枝。而且和一般的后剪枝不同， XGBoost 的后剪枝是不需要驗證集的。 不過我並不覺得這是“純粹”的后剪枝，因為一般還是要預先限制最大深度的呵呵。

說了這么多 XGBoost 的優點，其當然也有不完美之處，因為要在訓練之前先對每個特征進行預排序並將結果存儲起來，對於空間消耗較大。另外雖然相比傳統的 GBDT 速度是快了很多，但和后來的 LightGBM 比起來還是慢了不少，不知以后還會不會出現更加快的 Boosting 實現。

---

Reference

1. Tianqi Chen and Carlos Guestrin. XGBoost: A Scalable Tree Boosting System
2. Tianqi Chen. Introduction to Boosted Trees
3. https://zhuanlan.zhihu.com/p/34534004
4. https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/
5. https://stackoverflow.com/questions/52672116/what-is-xgboost-pruning-step-doing

附錄： 泰勒公式

XGBoost 推導的關鍵一步是二階泰勒展開，這里作一下拓展。泰勒公式的簡單解釋就是用多項式函數去逼近原函數，因為用多項式函數往往求解更加容易，而多項式中各項的系數則為原函數在某一點的n階導數值除以n階乘。這里力薦 3Blue1Brown 關於泰勒公式的視頻 微積分的本質 - 泰勒級數 ，講得非常形象 。

已知函數 \(f(x)\) 在 \(x=x_0\) 處n階可導，那么：

\[\begin{aligned} f(x) &= f(x_0) + f'(x_0)(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 + \cdots + \frac{f^{(n)}(x_0)}{n!}(x-x_0)^n \\ & = \sum\limits_{n=0}^{\infty}\frac{f^{(n)}x_0}{n!}(x - x_0)^n \end{aligned} \]

例如，\(x_0 = 0, \;\; f(x) = e^x\)時，\(e^x = \sum\limits_{n=0}^{\infty}\frac{x^n}{n!} = 1+x+\frac{x^2}{2!} + \frac{x^3}{3!} + \cdots\)

在機器學習中泰勒公式的一個典型應用就是牛頓法。在牛頓法中，將損失函數 \(L(\theta)\) 在 \(\theta^{k}\) 處進行二階泰勒展開(考慮 \(\theta\) 為標量)：

\[L(\theta) =L(\theta^k) + L'(\theta^k)(\theta - \theta^k) + \frac{L''(\theta^k)}{2}(\theta - \theta^k)^2 \]

將一階導和二階導分別記為 \(g_k\) 和 \(h_k\)，為使上式最小化，令

\[\frac{\partial\left[g(\theta - \theta^k) + \frac{h}{2}(\theta - \theta^k)^2\right]}{\partial \,\theta} = 0, \qquad 求得\; \theta-\theta^k = -\frac{g_k}{h_k}\;, \]

則參數更新公式為 \(\theta^{k+1} = \theta^k - \frac{g_k}{h_k} \;\)，推廣到向量形式則為 \(\boldsymbol{\theta}^{k+1} = \boldsymbol{\theta}^k - \mathbf{H}^{-1}_k\mathbf{g}_k\;\)，\(\mathbf{H}_k\) 為海森矩陣(Hessian matrix)，\(\mathbf{g}_k\) 為梯度向量：

\[\nabla_\theta L = \mathbf{g}_k = \begin {bmatrix} \frac{\partial L}{\partial \theta_1}\\ \frac{\partial L}{\partial\theta_2}\\ \vdots \\ \frac{\partial L}{\partial\theta_n} \end {bmatrix} \;\;,\;\; \nabla_{\theta}^2 L = \mathbf{H}_k = \begin {bmatrix} \frac{\partial^2 L}{\partial\, \theta_1^2} & \frac{\partial^2 L}{\partial\,\theta_1 \partial\,\theta_2} & \cdots & \frac{\partial^2 L}{\partial\,\theta_1 \partial\,\theta_n} \\ \frac{\partial^2 L}{\partial\,\theta_2 \partial\,\theta_1} & \frac{\partial^2 L}{\partial\, \theta_2^2} & \cdots & \frac{\partial^2 L}{\partial\,\theta_2 \partial\,\theta_n} \\ \vdots & \vdots & \ddots \\ \frac{\partial^2 L}{\partial\,\theta_n \partial\,\theta_1} & \frac{\partial^2 L}{\partial\,\theta_n \partial\,\theta_2} & \cdots & \frac{\partial^2 L}{\partial\, \theta_n^2} \end {bmatrix} \]

/