從CART到經典GBDT到XGBoost——樹模型重溫

正如從AdaBoost到GBM到LambdaMART——Boosting方法重溫所提到的，GBM中應用最廣泛的模型是GBDT，即以CART回歸樹為基函數的梯度提升樹模型。XGBoost作為實現Gradient Boosting的主流框架之一，除了在計算性能上作優化之外，在演算法層面也作了一些改變。這裡挑若干值得一提的點稍作對比和總結。

減輕過擬合的手段

機器學習上的優化和一般優化問題最大的區別在於真實數據分佈未知，因此無法對真實的優化目標期望泛化誤差直接進行優化，只能從分佈中抽樣出訓練集，通過優化訓練集上的經驗風險函數來間接優化期望泛化誤差。這就意味著在經驗風險函數上的最優解不一定是期望泛化誤差的最優解，通常的情況是在經驗風險函數上表現良好，但基於測試集來做估計的泛化誤差高，也即所謂的過擬合。從優化的角度來看，可以通過一些約束將經驗風險函數上的最優解「往隔壁挪一挪」，讓它沒那麼完美，但在泛化誤差上表現更好，從而減輕過擬合。

決策樹本身是一個容易高方差的模型，因為它挑選切分特徵和切分點是完全基於訓練數據的分佈來的，也就是多次抽樣下訓練集之間的差異很容易導致得到的決策樹也顯著不同，因此方差理論上會比較大。

理論上GBDT作為多個決策樹的組合，會得到一個相對單棵決策樹更平滑的模型，因此能一定程度上降低方差，也就沒那麼容易過擬合。但這並非意味著GBDT不存在過擬合，畢竟無論如何它都只能基於訓練數據來學習。Friedman在提出GBM時給出了一些降低GBDT過擬合程度的方法，而XGBoost也相應地做了一些補充和改良。

首先將GBDT模型記為 $F(x)=sum_{m=1}^Malpha_mg_m(x)$ ，其中 $g_m(x)=sum_{t=1}^T c_{mt}I(x in R_{mt})$ . 那麼 $F(x)=sum_{m=1}^Malpha_mg_m(x)=sum_{m=1}^M sum_{t=1}^T alpha_mc_{mt}I(x in R_{mt})$ . 又記 $w_{mt}=alpha_mc_{mt}$ ，因此 $F(x)=sum_{m=1}^M sum_{t=1}^T w_{mt}I(x in R_{mt})$ .

GBDT模型化簡之後可以看作擬合負梯度方向的是決定 $I(x in R_{mt})$ 的一系列規則，而 $w_{mt}$ 是步長。由於 $R_{mt}$ 之間的互斥關係，一個對應於唯一一個前進方向和步長。

1）控制迭代次數和步長

由於Boosting的過程可以看作是在函數空間通過line search方法求最優值解，因此為了減輕過擬合，可以在Boosting層面上可以控制迭代次數和步長，分別對應於和 $w_{mt}$ . 的最優值可以由交叉檢驗實驗得出，而控制 $w_{mt}$ 的常見方法有兩種：

shrinkage

即指定收縮因子，使收縮後的模型變為 $F(x)=sum_{m=1}^M etasum_{t=1}^T w_{mt}I(x in R_{mt})$ 。

這種做法可以從多個角度進行解釋，比如從Boosting角度可以看成是通過引入學習率對步長進行收縮；從決策樹角度看是對葉節點輸出值作收縮；若是將 $sum_{t=1}^T w_{mt}I(x in R_{mt})$ 記為一棵新的決策樹， $F(x)=sum_{m=1}^M eta f_m(x)$ 又可以看作是降低每個基函數的影響力。

L2正則化

即在損失函數上加上 $w_{mt}$ 對應向量的L2模，使得第次boost的損失函數為 $sum_{i=1}^nL(y_i, F^{m-1}(x_i) +sum_{t=1}^T w_{mt}I(x in R_{mt})) + frac{1}{2}lambda||vec{w}||^2$ 。

這種方法相當於在原優化目標中加入 $||vec{w}||^2<C$ 的約束，其中常數與存在某種關聯。相比起shrinkage方法的人為指定收縮因子，這種方法對各個分量產生的影響是不均等的，比如說對絕對值大的分量 $w_{mt}$ 上的收縮力度會大一些。

2）控制決策樹的結構

以上這些方法能影響的是樹的數目、葉節點輸出值的平滑程度，若要通過調整樹的結構來降低過擬合程度，常用的方法是從隨機森林借鑒過來的row subsampling和column subsampling，分別控制每棵樹採用的樣本集和特徵集，從而使得樹的結構與使用全集時有所不同。

相對於經典CART樹上的後剪枝，以上這些策略更為高效，畢竟GBDT需要擬合多個CART樹，假如讓每棵CART樹都完全生長再剪枝，一是低效，二是讓單棵CART樹儘可能榨乾訓練數據的信息對於Boosting方法來說沒有必要。

雖然XGBoost也支持shrinkage和subsampling等方法，但使其與經典GBDT區別開來的還是它在經典GBDT的損失函數上加上L2正則項和決策樹大小，即損失函數變為 $sum_{i=1}^nL(y_i, F^{m-1}(x_i) +sum_{t=1}^T w_{mt}I(x in R_{mt})) + frac{1}{2}lambda||vec{w}||^2 + u T$ ，其中為當前這棵決策樹的葉節點數目。這也使得它的參數估計方法與經典GBDT有所不同。

正則化損失函數下的參數估計

傳統GBDT每一步的損失函數為 $sum_{i=1}^nL(y_i, F^{m-1}(x_i) +sum_{t=1}^T w_{mt}I(x in R_{mt}))$

在參數估計時只需要找到最優的 $R_{mt}$ 來擬合負梯度方向 $abla L |_{F^{m-1}(x)}$ ，再求出最優步長 $vec{w_m}^*=argmin_{vec{w_m}} sum_{i=1}^nL(y_i, F^{m-1}(x_i) +sum_{t=1}^T w_{mt}I(x in R_{mt}))$ 。由於 $R_{mt}$ 之間的互斥關係， $w_{mt}^*=argmin_{w_{mt}} sum_{x_i in R_{mt}}L(y_i, F^{m-1}(x_i) + w_{mt})$ ，即各分量的最優解僅與其對應的子區域上的樣本有關。

而正如上文所提到的，XGBoost使用正則化的目標函數 $sum_{i=1}^nL(y_i, F^{m-1}(x_i) +sum_{t=1}^T w_{mt}I(x in R_{mt})) + frac{1}{2}lambda||vec{w}||^2 + u T$

這個目標函數對並非連續可導的，因此求最優 $R_{mt}$ 無法直接通過擬合負梯度方向 $abla L |_{F^{m-1}(x)}$ 來完成。但實際上負梯度方向也是通過對損失函數作一階泰勒近似推導而來的，因此不妨直接通過對損失函數作二階泰勒近似推導出最優的 $w_{mt}$ 和 $R_{mt}$ 。

1）估計 $w_{mt}$

假設現在處於第次Boosting，當前的樹結構是，劃分出了個子區域。同樣由子區域之間的互斥關係得到在當前樹結構下 $w_{mt}$ 的最優解：

$w_{mt}^* = argmin_{w_{mt}}{sum_{x_i in R_{mt}} L(y_i, F^{m-1}(x_i) + w_{mt}) + frac{1}{2}lambda w_{mt}^2 + u}$

對 $L(y_i, F^{m-1}(x_i) + w_{mt})$ 作二階泰勒近似得到 $w_{mt}^* = argmin_{w_{mt}}{sum_{x_i in R_{mt}} L(y_i, F^{m-1}(x_i)) + w_{mt} * abla L + frac{1}{2}w_{mt}^2* abla^2L + frac{1}{2}lambda w_{mt}^2 + u}\=argmin_{w_{mt}}{sum_{x_i in R_{mt}} w_{mt} * abla L + frac{1}{2}w_{mt}^2* abla^2L + frac{1}{2}lambda w_{mt}^2}\=argmin_{w_{mt}} {sum_{x_i in R_{mt}} w_{mt} * abla L + frac{1}{2}( abla^2L+lambda)w_{mt}^2}$

很容易可以求出 $w_{mt}^*=-frac{sum_{x_i in R_{mt}} abla L}{sum_{x_i in R_{mt}} abla^2L+lambda}$

相比起無正則化約束時的經典GBDT的解 $-frac{1}{|R_{mt}|}sum_{x_i in R_{mt}} abla L$

$w_{mt}^*$ 可以直觀而不嚴謹地看作是用Hessian矩陣和正則化係數對負梯度方向作了修正。

2）估計 $R_{mt}$

由於目前是處於第次boost，因此 $L(y_i, F^{m-1}(x_i))$ 是常數，因此在考慮 $sum_{t=1}^T w_{mt}I(x in R_{mt})$ 的生成時可以忽略 $L(y_i, F^{m-1}(x_i))$ ，因此將 $w_{mt}^*$ 代入到 $w_{mt} * abla L + frac{1}{2}w_{mt}^2* abla^2L + frac{1}{2}lambda w_{mt}^2 + u$ 就得到結構下 $R_{mt}$ 這個子區域的損失函數極小值，個子區域的整體損失函數為 $sum_{t=1}^Tw_{mt}^* * abla L + frac{1}{2}w_{mt}^{*2}* abla^2L + frac{1}{2}lambda w_{mt}^{*2} + u T\=sum_{t=1}^T-frac{1}{2}frac{(sum_{x_i in R_{mt}} abla L)^2}{sum_{x_i in R_{mt}} abla^2L + lambda} + u T$

假設現在在結構的基礎上要對某一個葉節點 $R_{mt}$ 繼續進行切分成兩個新的葉節點 $R_{mtL}$ 和 $R_{mtR}$ ，那麼切分之後的整體損失函數的下降值為

$frac{1}{2}(frac{(sum_{x_i in R_{mt_L}} abla L)^2}{sum_{x_i in R_{mt_L}} abla^2L + lambda} +frac{(sum_{x_i in R_{mt_R}} abla L)^2}{sum_{x_i in R_{mt_R}} abla^2L + lambda} -frac{(sum_{x_i in R_{mt_R}} abla L)^2}{sum_{x_i in R_{mt_R}} abla^2L + lambda} )- u$

使此下降值最大的即是這一次生長的最優切分特徵和最優切分點。使得下降值最大顯然也等價於生長後的損失函數最小，這與CART的生成邏輯是一致的。

另一方面，XGBoost利用了損失函數的二階展開來迭代求解，理論上這會比傳統GBDT的收斂速度更快。

對缺失值和無序離散型特徵的處理方式

缺失值和離散特徵本身也是CART樹生成過程中的的重要話題。之所以要強調「無序」是因為有序的離散特徵在切分時的做法可以和連續特徵一樣，而無序離散特徵的難點在於它總共有這麼多種切分可能，其中為該特徵的取值個數。

1）缺失值

先以CART回歸樹為例，生成過程決定切分特徵和切分點的目標函數是 $min_{j,s}{min_{c_1}sum_{x_iin R_1(j,s)} (y_i-c_1)^2 + min_{c_2} sum_{x_i in R_2(j,s)} (y_i - c_2)^2 }$

其中和根據在第維特徵上的值與切分點的大小關係來定義。

在訓練時假如在第維上的特徵缺失，那麼目標函數將無法考慮到這個樣本的信息，相當於將它丟棄；而在預測時，第維上的特徵缺失的樣本也無法得到預測結果。

CART處理缺失值的常見方式有兩種：

將缺失值看作是這個特徵的新取值，填充為「missing」或「NaN」（這是lightGBM採用的方式）
採用代理切分特徵。仍然是以回歸樹為例，每次基於非缺失數據上的目標函數找到最優的之後，再從以外的特徵尋找和的「切分行為」很相似的若干個代理特徵。在針對維特徵缺失的樣本做預測時，在這些代理特徵裡面找一個上不缺失且與相近程度最高的代理特徵來完成切分