整合學習之Boosting —— XGBoost
ofollow,noindex" target="_blank"> 整合學習之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.readthedocs.io/en/latest/index.html" rel="nofollow,noindex" target="_blank"> 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 的學習演算法：
- 初始化 \(f_0(x)\)
- 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
- Tianqi Chen and Carlos Guestrin. XGBoost: A Scalable Tree Boosting System
- Tianqi Chen. Introduction to Boosted Trees
- https://zhuanlan.zhihu.com/p/34534004
- https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/
- 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} \]
/