データ開発本部 データサイエンティストの佐藤です。
こちらの記事では、XGBoostやLightGBMなどの勾配ブースティング木の分割についてご紹介いたします。
はじめに
データサイエンティストとして仕事を始めてから、「勾配ブースティング木（Gradient Boosting Decision Tree, GBDT）」に触れる機会が増えました。
GBDTは機械学習の世界で非常に広く使われているアルゴリズムの一つで、分類問題でも回帰問題でも高い精度を発揮することで知られています。[１]
実際、XGBoost や LightGBM といったライブラリを使えば、GBDTの実装は驚くほど簡単です。学習を行うコードはなんと数行で完了することができます。
model = xgb.train(params, dtrain)
実際にライブラリを使って動かしているだけだと、アルゴリズム内部で何が行われているかを深く理解することは少し難しい部分もあります。特に「GBDTの木がどのような基準で分割されるのか」という仕組みや数式は、実装の便利さに隠れてあまり見えないポイントではないでしょうか？
そこで今回は少し立ち止まり、XGBoost: A scalable tree boosting system[２]という元論文を参考に、GBDTの分割基準がどのように導き出されるのかを数式にフォーカスして整理することにしました。図やグラフを用いた説明は他の方がわかりやすくまとめているので[３]、この記事では数式に絞っていきます。
勾配ブースティング木の学習イメージ。[３]より引用。
木の分割の公式の導出
$n$個のサンプルと$m$個の特徴量からなるデータセットは
$$D=\{(x_i,y_i )\},\: ( |D|=n,\: x_i∈\mathbb{R}^m, \: y_i∈\mathbb{R} )$$
と書けます。$\mathbb{R}$は実数の集合です。この時、GBDTを使った予測モデルは
$$y ̂_i=ϕ(x_i )=\sum_{k=1}^K f_k (x_i), f_k∈\mathscr{F}$$
です。
ここでは、
$T$：葉の枚数、$q$：木の構造、$ω$：葉ごとのスコア（重み）
$f_k$：木の構造とその葉の重みをもつ関数（弱学習器）
とします。
$y ̂_i$は$y_i$の予測値です。予測においては、なるべく簡単なモデルで真の値と予測値の差をなくすことが望ましいです。なので以下$L$を最小化することを考えます。
$$L≔\sum_i l(y_i,y ̂_i ) +\sum_k Ω(f_k)$$
ここで$l(y_i,y ̂_i )$は損失関数、$Ω(f_k)$は木の複雑さに対する正則化項（ペナルティ）です。
$Ω(f_k)$は
$$Ω(f_k )=γT+\frac{1}{2}λ\sum_{j=1}^T ω_j^2$$
であり、$γT$は葉の数に対するペナルティ、$\frac{1}{2}λ\sum_{j=1}^T ω_j^2$は葉の重みに対するペナルティ（$L_2$正則化）です。
ここで、葉$t-1$枚での予測値から新しい木$f_t$を追加することで予測値が改善するとき、
$$y ̂_i^{(t)}=y ̂_i^{(t-1)}+f_t (x_i)$$
と書くことができます。これを$L$に代入すると
$$L=\sum_i l(y_i,y ̂_i^{(t-1)}+f_t (x_i) )+Ω(f_T )$$
です。ちなみに$l(y_i,y ̂_i^{(t-1)}+f_t (x_i) )$をテイラー展開すると
$$l(y_i,y ̂_i^{(t-1)}+f_t (x_i) )\fallingdotseq l(y_i,y ̂_i^{(t-1)})+g_i f_t (x_i )+\frac{1}{2}h_i f_t (x_i )^2$$
$$g_i=\frac{∂}{∂y ̂_i^{(t-1)}} l(y_i,y ̂_i^{(t-1)}),h_i=\frac{∂^2}{∂y ̂_i^{(t-1)^2}} l(y_i,y ̂_i^{(t-1)})$$
であるので
$$L\fallingdotseq\sum_i [l(y_i,y ̂_i^{(t-1)})+g_i f_t (x_i )+\frac{1}{2}h_i f_t (x_i )^2] +Ω(f_T )$$
の最小化を考えます。また、$Ω(f_k )=γT+\frac{1}{2}λ\sum_{j=1}^T ω_j^2$なので
$$L \fallingdotseq\sum_i [l(y_i,y ̂_i^{(t-1)})+g_i f_t (x_i )+\frac{1}{2}h_i f_t (x_i )^2] +γT+\frac{1}{2}λ\sum_{j=1}^T ω_j^2$$
$I_j={i|q(x_i )=j}$(どの葉に属するかを表す)、$f_t (x_i )=ω_i$なので、
$$\sum_{j=1}^T[(\sum_{i∈I_j} g_i ) ω_j+\frac{1}{2} (\sum_{i∈I_j} h_i +λ) ω_j^2)] +γT$$
これを最小化する$ω_i^*$は、
$$ω_i^*=-\frac{\sum_{i∈I_j} g_i}{\sum_{i∈I_j} h_i +λ}$$
となり、つまり最適値は
$$L=-\frac{1}{2} \sum_{j=1}^T \frac{(\sum_{i∈I_j} g_i)^2}{\sum_{i∈I_j} h_i +λ} +γT$$
分割後の集合$I_R$,$I_L$を考えると損失の減少は
$$L^{(t)}-L ^{(t-1)}=ΔL=\frac{1}{2}[\frac{(\sum_{i∈I_L} g_i)^2}{\sum_{i∈I_L} h_i +λ}+\frac{(\sum_{i∈I_R} g_i)^2}{\sum_{i∈I_R} h_i +λ}-\frac{(\sum_{i∈I_L∪I_R} g_i)^2}{\sum_{i∈I_L∪I_R} h_i +λ}]-γ$$
であり、これを最大化するような分割点を見つければそれが最良の木の分割となります。
分割ロジックからわかるより深い使い方
この導出の注目するべきポイントは、最終式の右辺に含まれる$g_i$も$h_i$も損失関数$l$から得られることです。つまり、条件さえ満たせば損失関数$l$を自由に定義してよいということです。
例えば「工場の温度管理で、一定以上の温度を超えると全て台無しになる」ような場面で、予測が上振れした場合コストを指数関数的に増大させたい場合は以下のような損失関数を考えます。例えば
$$l(y,y ̂ )=\exp(c \cdot (y-y ̂))-c \cdot (y-y ̂)-1$$
を考えれば、そのような挙動が再現できます。
def my_loss_func(y_pred, dtrain) -> [np.array, np.array]:
"""
XGBoostに渡すカスタム目的関数。
予測値と実際の値から、勾配(1階微分)とヘシアン(2階微分)を計算して返します。
"""
y_true = dtrain.get_label()
c = #{cの値}
x = y_pred - y_true
x_clipped = np.clip(x, a_min=-100, a_max=100)
#1階微分: g_i
grad = c * np.exp(c * x_clipped) - c
#2階微分: h_i
hess = (c ** 2) * np.exp(c * x_clipped)
return grad, hess
####上記で定義したのち以下のように用いる####
model = xgb.train(
params,
dtrain,
obj = my_loss_func #←ここを追記
)
以上、勾配ブースティング木（GBDT）の分割基準がどのように導出されるのかを、論文を基に数式ベースで整理しました。
このことにより損失関数を自由にカスタムしビジネスルールに合わせることも可能である事がわかります。
まとめ
GBDTは簡単に使えますが、その中身をより深く理解すればビジネス要件に沿ったより効果的なモデル構築や問題解決も可能になります。今後も本記事のように基礎理論を実務につなげられるような発信していけたらと思います。
参考文献
[1] A Comprehensive Survey on Prediction Models and the Impact of XGBoost
