簡單易學的機器學習演算法——梯度提升決策樹GBDT
梯度提升決策樹（Gradient Boosting Decision Tree，GBDT）演算法是近年來被提及比較多的一個演算法，這主要得益於其演算法的效能，以及該演算法在各類資料探勘以及機器學習比賽中的卓越表現，有很多人對GBDT演算法進行了開原始碼的開發，比較火的是陳天奇的XGBoost和微軟的LightGBM。
一、監督學習
1、監督學習的主要任務
監督學習是機器學習演算法中重要的一種，對於監督學習，假設有mm個訓練樣本：
{(X(1),y(1)),(X(2),y(2)),⋯,(X(m),y(m))}{(X(1),y(1)),(X(2),y(2)),⋯,(X(m),y(m))}
其中，X(i)={x(i)1,x(i)2,⋯,x(i)m}X(i)={x1(i),x2(i),⋯,xm(i)}稱為第ii個樣本的特徵，y(1)y(1)稱為第ii個樣本的標籤，樣本標籤可以為離散值，如分類問題；也可以為連續值，如迴歸問題。在監督學習中，利用訓練樣本訓練出模型，該模型能夠實現從樣本特徵X(i)X(i)到樣本標籤y(1)y(1)的對映，即：
X(i)→Fy(i)X(i)→Fy(i)
為了能夠對對映F(X)F(X)進行求解，通常對模型設定損失函式L(y,F(X))L(y,F(X))，並求得在損失函式最小的情況下的對映為最好的對映：
F∗=argminF(X)L(y,F(X))F∗=argminF(X)L(y,F(X))
對於一個具體的問題，如線性迴歸問題，其對映函式的形式為：
F(X;W)=WX=w0+w1x1+w2x2+⋯+wnxnF(X;W)=WX=w0+w1x1+w2x2+⋯+wnxn
此時對於最優對映函式F(X;W)F(X;W)的求解，實質是對對映函式中的引數WW的求解。對於引數的求解方法有很多，如梯度下降法。
2、梯度下降法
梯度下降法（Gradient Descent，GD）演算法是求解最優化問題最簡單、最直接的方法。梯度下降法是一種迭代的優化演算法，對於優化問題：
minf(w)minf(w)
其基本步驟為：
- 隨機選擇一個初始點w0w0
- 重複以下過程：
- 決定下降的方向：di=−∂∂wf(w)∣widi=−∂∂wf(w)∣wi
- 選擇步長ρρ
- 更新：wi+1=wi+ρ⋅diwi+1=wi+ρ⋅di
- 直到滿足終止條件
梯度下降法的具體過程如下圖所示：
由以上的過程，我們可以看出，對於最終的最優解w∗w∗，是由初始值w0w0經過MM代的迭代之後得到的，在這裡，設w0=d0w0=d0，則w∗w∗為：
w∗=∑i=0Mρi⋅diw∗=∑i=0Mρi⋅di
3、在函式空間的優化
以上是在指定的函式空間中對最優函式進行搜尋，那麼，能否直接在函式空間（function space）中查詢到最優的函式呢？根據上述的梯度下降法的思路，對於模型的損失函式L(y,F(X))L(y,F(X))，為了能夠求解出最優的函式F∗(X)F∗(X)，首先，設定初始值為：
F0(X)=f0(X)F0(X)=f0(X)
以函式F(X)F(X)作為一個整體，對於每一個樣本X(i)X(i)，都存在對應的函式值F(X(i))F(X(i))。與梯度下降法的更新過程一致，假設經過MM代，得到最有的函式F∗(X)F∗(X)為：
F∗(X)=∑i=0Mfi(X)F∗(X)=∑i=0Mfi(X)
其中，fi(X)fi(X)為：
fi(X)=−ρigm(X)fi(X)=−ρigm(X)
其中，gm(X)=[∂L(y,F(X))∂F(X)]F(X)=Fm−1(X)gm(X)=[∂L(y,F(X))∂F(X)]F(X)=Fm−1(X)。
由上述的過程可以得到函式F(X)F(X)的更新過程：
Fm(X)=∑i=0mfi(X)Fm(X)=∑i=0mfi(X)
與上面類似，函式f(X)f(X)是由引數aa決定的，即：
f(X)=−ρ⋅h(X;a)f(X)=−ρ⋅h(X;a)
二、Boosting
1、整合方法之Boosting
Boosting方法是整合學習中重要的一種方法，在整合學習方法中最主要的兩種方法為Bagging和Boosting，在Bagging中，通過對訓練樣本重新取樣的方法得到不同的訓練樣本集，在這些新的訓練樣本集上分別訓練學習器，最終合併每一個學習器的結果，作為最終的學習結果，Bagging方法的具體過程如下圖所示：
在Bagging方法中，最重要的演算法為隨機森林Random Forest演算法。由以上的圖中可以看出，在Bagging方法中，bb個學習器之間彼此是相互獨立的，這樣的特點使得Bagging方法更容易並行。與Bagging方法不同，在Boosting演算法中，學習器之間是存在先後順序的，同時，每一個樣本是有權重的，初始時，每一個樣本的權重是相等的。首先，第11個學習器對訓練樣本進行學習，當學習完成後，增大錯誤樣本的權重，同時減小正確樣本的權重，再利用第22個學習器對其進行學習，依次進行下去，最終得到bb個學習器，最終，合併這bb個學習器的結果，同時，與Bagging中不同的是，每一個學習器的權重也是不一樣的。Boosting方法的具體過程如下圖所示：
在Boosting方法中，最重要的方法包括：AdaBoost和GBDT。
2、Gradient Boosting
由上圖所示的Boosting方法中，最終的預測結果為bb個學習器結果的合併：
f(X)=∑j=1bθjφj(X)f(X)=∑j=1bθjφj(X)
這與上述的在函式空間中的優化類似：
Fm(X)=∑i=0m−ρi⋅h(X;ai)Fm(X)=∑i=0m−ρi⋅h(X;ai)
根據如上的函式空間中的優化可知，每次對每一個樣本的訓練的值為：
y¯i=[∂L(yi,F(X(i)))∂F(X(i))]F(X)=Fm−1(X)y¯i=[∂L(yi,F(X(i)))∂F(X(i))]F(X)=Fm−1(X)
上建立模型，由於上述是一個求解梯度的過程，因此也稱為基於梯度的Boost方法，其具體過程如下所示：
三、Gradient Boosting Decision Tree
在上面簡單介紹了Gradient Boost框架，梯度提升決策樹Gradient Boosting Decision Tree是Gradient Boost框架下使用較多的一種模型，在梯度提升決策樹中，其基學習器是分類迴歸樹CART，使用的是CART樹中的迴歸樹。
1、分類迴歸樹CART
分類迴歸樹CART演算法是一種基於二叉樹的機器學習演算法，其既能處理迴歸問題，又能處理分類為題，在梯度提升決策樹GBDT演算法中，使用到的是CART迴歸樹演算法，對於CART樹演算法的更多資訊，可以參考簡單易學的機器學習演算法——分類迴歸樹CART。
對於一個包含了mm個訓練樣本的迴歸問題，其訓練樣本為：
{(X(1),y(1)),(X(2),y(2)),⋯,(X(m),y(m))}{(X(1),y(1)),(X(2),y(2)),⋯,(X(m),y(m))}
其中，X(i)X(i)為nn維向量，表示的是第ii個樣本的特徵，y(i)y(i)為樣本的標籤，在迴歸問題中，標籤y(i)y(i)為一系列連續的值。此時，利用訓練樣本訓練一棵CART迴歸樹：
- 開始時，CART樹中只包含了根結點，所有樣本都被劃分在根結點上：
此時，計算該節點上的樣本的方差（此處要乘以mm），方差表示的是資料的波動程度。那麼，根節點的方差的mm倍為：
s2⋅m=(y(1)−y¯)2+(y(2)−y¯)2+⋯+(y(m)−y¯)2s2⋅m=(y(1)−y¯)2+(y(2)−y¯)2+⋯+(y(m)−y¯)2
其中，y¯y¯為標籤的均值。此時，從nn維特徵中選擇第jj維特徵，從mm個樣本中選擇一個樣本的值：xjxj作為劃分的標準，當樣本ii的第jj維特徵小於等於xjxj時，將樣本劃分到左子樹中，否則，劃分到右子樹中，通過以上的操作，劃分到左子樹中的樣本個數為m1m1，劃分到右子樹的樣本的個數為m2=m−m1m2=m−m1，其劃分的結果如下圖所示：
那麼，什麼樣本的劃分才是當前的最好劃分呢？此時計算左右子樹的方差之和：s21⋅m1+s22⋅m2s12⋅m1+s22⋅m2：
s21⋅m1+s22⋅m2=∑X(i)∈left(y(i)−y¯1)2+∑X(j)∈right(y(j)−y¯2)2s12⋅m1+s22⋅m2=∑X(i)∈left(y(i)−y¯1)2+∑X(j)∈right(y(j)−y¯2)2
其中，y¯1y¯1為左子樹中節點標籤的均值，同理，y¯2y¯2為右子樹中節點標籤的均值。選擇其中s21⋅m1+s22⋅m2s12⋅m1+s22⋅m2最小的劃分作為最終的劃分，依次這樣劃分下去，直到得到最終的劃分，劃分的結果為：
注意：對於上述最優劃分標準的選擇，以上的計算過程可以進一步優化。
首先，對於s2⋅ms2⋅m：
s2⋅m=∑X(i)(y(i)−y¯)2=∑X(i)((y(i))2−2y(i)⋅y¯+(y¯)2)=∑X(i)(y(i))2−2m(∑X(i)y(i))2+1m(∑X(i)y(i))2=∑X(i)(y(i))2−1m(∑X(i)y(i))2s2⋅m=∑X(i)(y(i)−y¯)2=∑X(i)((y(i))2−2y(i)⋅y¯+(y¯)2)=∑X(i)(y(i))2−2m(∑X(i)y(i))2+1m(∑X(i)y(i))2=∑X(i)(y(i))2−1m(∑X(i)y(i))2
而對於s21⋅m1+s22⋅m2s12⋅m1+s22⋅m2：
s21⋅m1+s22⋅m2=∑X(i)∈left(y(i)−y¯1)2+∑X(j)∈right(y(j)−y¯2)2=∑X(i)∈left((y(i))2−2y(i)⋅y¯1+(y¯1)2)+∑X(j)∈right((y(j))2−2y(j)⋅y¯2+(y¯2)2)s12⋅m1+s22⋅m2=∑X(i)∈left(y(i)−y¯1)2+∑X(j)∈right(y(j)−y¯2)2=∑X(i)∈left((y(i))2−2y(i)⋅y¯1+(y¯1)2)+∑X(j)∈right((y(j))2−2y(j)⋅y¯2+(y¯2)2)
=∑X(i)(y(i))2−2m1⎛⎝∑X(i)∈lefty(i)⎞⎠2+1m1⎛⎝∑X(i)∈lefty(i)⎞⎠2−2m2⎛⎝∑X(j)∈righty(j)⎞⎠2+1m2⎛⎝∑X(j)∈righty(j)⎞⎠2=∑X(i)(y(i))2−1m1⎛⎝∑X(i)∈lefty(i)⎞⎠2−1m2⎛⎝∑X(j)∈righty(j)⎞⎠2=∑X(i)(y(i))2−2m1(∑X(i)∈lefty(i))2+1m1(∑X(i)∈lefty(i))2−2m2(∑X(j)∈righty(j))2+1m2(∑X(j)∈righty(j))2=∑X(i)(y(i))2−1m1(∑X(i)∈lefty(i))2−1m2(∑X(j)∈righty(j))2
通過以上的過程，我們發現，劃分前，記錄節點的值為：
1m(∑X(i)y(i))21m(∑X(i)y(i))2
當劃分後，兩個節點的值的和為：
1m1⎛⎝∑X(i)∈lefty(i)⎞⎠2+1m2⎛⎝∑X(j)∈righty(j)⎞⎠21m1(∑X(i)∈lefty(i))2+1m2(∑X(j)∈righty(j))2
最好的劃分，對應著兩個節點的值的和的最大值。
2、GBDT——二分類
在梯度提升決策樹GBDT中，通過定義不同的損失函式，可以完成不同的學習任務，二分類是機器學習中一類比較重要的分類演算法，在二分類中，其損失函式為：
L(y,F)=log(1+exp(−2yF)),y∈{−1,1}L(y,F)=log(1+exp(−2yF)),y∈{−1,1}
套用上面介紹的GB框架，得到下述的二分類GBDT的演算法：
在構建每一棵CART迴歸樹的過程中，對一個樣本的預測值應與y~y~儘可能一致，對於y~y~，其計算過程為：
y~(i)=−[∂L(y(i),F(X(i)))∂F(X(i))]F(X)=Fm−1(X)=−[∂log(1+exp(−2y(i)F(X(i))))∂F(X(i))]F(X)=Fm−1(X)=−[11+exp(−2y(i)F(X(i)))⋅exp(−2y(i)F(X(i)))⋅(−2y(i))]F(X)=Fm−1(X)y~(i)=−[∂L(y(i),F(X(i)))∂F(X(i))]F(X)=Fm−1(X)=−[∂log(1+exp(−2y(i)F(X(i))))∂F(X(i))]F(X)=Fm−1(X)=−[11+exp(−2y(i)F(X(i)))⋅exp(−2y(i)F(X(i)))⋅(−2y(i))]F(X)=Fm−1(X)
=2y(i)⋅exp(−2y(i)F(X(i)))1+exp(−2y(i)F(X(i)))F(X)=Fm−1(X)=2y(i)1+exp(2y(i)Fm−1(X(i)))=2y(i)⋅exp(−2y(i)F(X(i)))1+exp(−2y(i)F(X(i)))F(X)=Fm−1(X)=2y(i)1+exp(2y(i)Fm−1(X(i)))
在y~y~（通常有的地方稱為殘差，在這裡，更準確的講是梯度下降的方向）上構建CART迴歸樹。最終將每一個訓練樣本劃分到對應的葉子節點中，計算此時該葉子節點的預測值：
γjm=argminγ∑X(i)∈Rjmlog(1+exp(−2y(i)(Fm−1(X(i))+γ)))γjm=argminγ∑X(i)∈Rjmlog(1+exp(−2y(i)(Fm−1(X(i))+γ)))
由Newton-Raphson迭代公式可得：
γjm=∑X(i)∈Rjmy~(i)∑X(i)∈Rjm∣∣y~(i)∣∣(2−∣∣y~(i)∣∣)γjm=∑X(i)∈Rjmy~(i)∑X(i)∈Rjm|y~(i)|(2−|y~(i)|)
- GBDT訓練的主要程式碼為：
void GBDT::fit(Problem const &Tr, Problem const &Va)
{
bias = calc_bias(Tr.Y); //用於初始化的F
std::vector<float> F_Tr(Tr.nr_instance, bias), F_Va(Va.nr_instance, bias);
Timer timer;
printf("iter time tr_loss va_loss\n");
// 開始訓練每一棵CART樹
for(uint32_t t = 0; t < trees.size(); ++t)
{
timer.tic();
std::vector<float> const &Y = Tr.Y;
std::vector<float> R(Tr.nr_instance), F1(Tr.nr_instance); // 記錄殘差和F
#pragma omp parallel for schedule(static)
for(uint32_t i = 0; i < Tr.nr_instance; ++i)
R[i] = static_cast<float>(Y[i]/(1+exp(Y[i]*F_Tr[i]))); //計算殘差，或者稱為梯度下降的方向
// 利用上面的殘差值，在此函式中構造一棵樹
trees[t].fit(Tr, R, F1); // 分類樹的生成
double Tr_loss = 0;
// 用上面訓練的結果更新F_Tr，並計算log_loss
#pragma omp parallel for schedule(static) reduction(+: Tr_loss)
for(uint32_t i = 0; i < Tr.nr_instance; ++i)
{
F_Tr[i] += F1[i];
Tr_loss += log(1+exp(-Y[i]*F_Tr[i]));
}
Tr_loss /= static_cast<double>(Tr.nr_instance);
// 用上面訓練的結果預測測試集，列印log_loss
#pragma omp parallel for schedule(static)
for(uint32_t i = 0; i < Va.nr_instance; ++i)
{
std::vector<float> x = construct_instance(Va, i);
F_Va[i] += trees[t].predict(x.data()).second;
}
double Va_loss = 0;
#pragma omp parallel for schedule(static) reduction(+: Va_loss)
for(uint32_t i = 0; i < Va.nr_instance; ++i)
Va_loss += log(1+exp(-Va.Y[i]*F_Va[i]));
Va_loss /= static_cast<double>(Va.nr_instance);
printf("%4d %8.1f %10.5f %10.5f\n", t, timer.toc(), Tr_loss, Va_loss);
fflush(stdout);
}
}- CART迴歸樹的訓練程式碼為：
void CART::fit(Problem const &prob, std::vector<float> const &R, std::vector<float> &F1){
uint32_t const nr_field = prob.nr_field; // 特徵的個數
uint32_t const nr_sparse_field = prob.nr_sparse_field;
uint32_t const nr_instance = prob.nr_instance; // 樣本的個數
std::vector<Location> locations(nr_instance); // 樣本資訊
#pragma omp parallel for schedule(static)
for(uint32_t i = 0; i < nr_instance; ++i)
locations[i].r = R[i]; // 記錄每一個樣本的殘差
for(uint32_t d = 0, offset = 1; d < max_depth; ++d, offset *= 2){// d:深度
uint32_t const nr_leaf = static_cast<uint32_t>(pow(2, d)); // 葉子節點的個數
std::vector<Meta> metas0(nr_leaf); // 葉子節點的資訊
for(uint32_t i = 0; i < nr_instance; ++i){
Location &location = locations[i]; //第i個樣本的資訊
if(location.shrinked)
continue;
Meta &meta = metas0[location.tnode_idx-offset]; //找到對應的葉子節點
meta.s += location.r; //殘差之和
++meta.n;
}
std::vector<Defender> defenders(nr_leaf*nr_field); //記錄每一個葉節點的每一維特徵
std::vector<Defender> defenders_sparse(nr_leaf*nr_sparse_field);
// 針對每一個葉節點
for(uint32_t f = 0; f < nr_leaf; ++f){
Meta const &meta = metas0[f]; // 葉子節點
double const ese = meta.s*meta.s/static_cast<double>(meta.n); //該葉子節點的ese
for(uint32_t j = 0; j < nr_field; ++j)
defenders[f*nr_field+j].ese = ese;
for(uint32_t j = 0; j < nr_sparse_field; ++j)
defenders_sparse[f*nr_sparse_field+j].ese = ese;
}
std::vector<Defender> defenders_inv = defenders;
std::thread thread_f(scan, std::ref(prob), std::ref(locations),
std::ref(metas0), std::ref(defenders), offset, true);
std::thread thread_b(scan, std::ref(prob), std::ref(locations),
std::ref(metas0), std::ref(defenders_inv), offset, false);
scan_sparse(prob, locations, metas0, defenders_sparse, offset, true);
thread_f.join();
thread_b.join();
// 找出最佳的ese，scan裡是每個欄位的最佳ese，這裡是所有欄位的最佳ese，賦值給相應的tnode
for(uint32_t f = 0; f < nr_leaf; ++f){
// 對於每一個葉節點都找到最好的劃分
Meta const &meta = metas0[f];
double best_ese = meta.s*meta.s/static_cast<double>(meta.n);
TreeNode &tnode = tnodes[f+offset];
for(uint32_t j = 0; j < nr_field; ++j){
Defender defender = defenders[f*nr_field+j];//每一個葉節點都對應著所有的特徵
if(defender.ese > best_ese)
{
best_ese = defender.ese;
tnode.feature = j;
tnode.threshold = defender.threshold;
}
defender = defenders_inv[f*nr_field+j];
if(defender.ese > best_ese)
{
best_ese = defender.ese;
tnode.feature = j;
tnode.threshold = defender.threshold;
}
}
for(uint32_t j = 0; j < nr_sparse_field; ++j)
{
Defender defender = defenders_sparse[f*nr_sparse_field+j];
if(defender.ese > best_ese)
{
best_ese = defender.ese;
tnode.feature = nr_field + j;
tnode.threshold = defender.threshold;
}
}
}
// 把每個instance都分配給樹裡的一個葉節點下
#pragma omp parallel for schedule(static)
for(uint32_t i = 0; i < nr_instance; ++i){
Location &location = locations[i];
if(location.shrinked)
continue;
uint32_t &tnode_idx = location.tnode_idx;
TreeNode &tnode = tnodes[tnode_idx];
if(tnode.feature == -1){
location.shrinked = true;
}else if(static_cast<uint32_t>(tnode.feature) < nr_field){
if(prob.Z[tnode.feature][i].v < tnode.threshold)
tnode_idx = 2*tnode_idx;
else
tnode_idx = 2*tnode_idx+1;
}else{
uint32_t const target_feature = static_cast<uint32_t>(tnode.feature-nr_field);
bool is_one = false;
for(uint64_t p = prob.SJP[i]; p < prob.SJP[i+1]; ++p)
{
if(prob.SJ[p] == target_feature)
{
is_one = true;
break;
}
}
if(!is_one)
tnode_idx = 2*tnode_idx;
else
tnode_idx = 2*tnode_idx+1;
}
}
}
// 用於計算gamma
std::vector<std::pair<double, double>>
tmp(max_tnodes, std::make_pair(0, 0));
for(uint32_t i = 0; i < nr_instance; ++i)
{
float const r = locations[i].r;
uint32_t const tnode_idx = locations[i].tnode_idx;
tmp[tnode_idx].first += r;
tmp[tnode_idx].second += fabs(r)*(1-fabs(r));
}
for(uint32_t tnode_idx = 1; tnode_idx <= max_tnodes; ++tnode_idx)
{
double a, b;
std::tie(a, b) = tmp[tnode_idx];
tnodes[tnode_idx].gamma = (b <= 1e-12)? 0 : static_cast<float>(a/b);
}
#pragma omp parallel for schedule(static)
for(uint32_t i = 0; i < nr_instance; ++i)
F1[i] = tnodes[locations[i].tnode_idx].gamma;// 重新更新F1的值
}