勾配ブースティング決定木を理解する

本記事では、機械学習コンペなどでよく見られる勾配ブースティング決定木(gradient boosting decision tree)を説明します。勾配ブースティング決定木は、MNISTデータに対して、ニューラルネットの最高精度と同等の精度を出したり、また高速な実装xgboostなどで有名な手法です。ライブラリを使用している方も多いと思いますが、意外とどのような構造になっているかを知らない人もいるかもしれません。

そこで、本記事では、決定木とは何か、というところから始めて、アンサンブル学習、勾配ブースティング決定木について見ていきます。

決定木

決定木(decision tree)は、データに対して一連の質問を与えることによって、目標に対する決定を下す予測モデルです。決定木の優れた点として、意味解釈性(interpretability)が挙げられます。予測モデルには、与えられたデータから得られる結果の過程がブラックボックス、すなわちどうしてその結果が導かれたのかがわかりづらいものが多く存在します。一方、決定木では、その木構造からどうやってその結果が得られたかを解釈しやすい、という人にとって非常に嬉しい性質を持ちます。

決定木の例

この図では、天気予報が晴れなら海に行く、曇りならランニング、雨なら映画を見るというように与えられたデータ(=天気情報)から次の行動(=その日の予定)を予測しています。これは、分類問題の例ですが、決定木は回帰問題、すなわち連続値を扱うような問題にも適用できます。例えば、身長であれば170cm以上、もしくはそれより下、というようにわけることで分類問題のように木を成長させることができます。分類問題と回帰問題、それぞれに対する決定木を分類木(classification tree)、回帰木(regression tree)と呼びます。

決定木は、その根から始めて、後述する情報利得(information gain) が最大となるようにデータを分割し、その葉が純粋(=葉の要素が全て同じクラスに属する)となるまで分割を行います。ただし、実際は葉が純粋になるまで決定木を成長させると、非常に構造が深くなり、過学習に陥りがちです。そこで、最大の深さを設定するなどにより、木が複雑になりすぎるのを防ぎます。

情報利得

決定木を分割する上で、はじめに考えなければいけないのはデータをどのように分割するか、ということです。この問いに一つの答えを与えてくれるのが情報利得であり、情報利得は分割された要素のばらつきの減少をあらわします。

例えば、クラスAに属するデータが40個、クラスBに属するデータが40個、計80個のデータを決定木を用いて分割することを考えます。

分割の例

さて、ここで決定木1と決定木2では、どちらの分割が、(語弊を恐れずに言えば)優れているのでしょうか。直感的には、決定木2であれば、右側の葉の要素はすべてクラスAに属しており、これ以上葉を増やさなくて済むので、良さそうに見えます。これを、式で裏付けていきます。

まず、決定木を二分決定木、すなわち親ノードは二つの子ノードまでしか持たない、を仮定します。このとき、I$I$を不純度、すなわち異なるクラスがどれだけ混じっているか、とすると、情報利得IG$IG$は以下の式であらわされます。

IG(Dp)=I(Dp)−NleftNpI(Dleft)−NrightNpI(Dright)$$IG(D_p)=I(D_p)-\frac{N_{left}}{N_p}I(D_{left})-\frac{N_{right}}{N_p}I(D_{right})$$

ここで、それぞれの変数は以下で定義されます。

• Dp$D_p$ : 親のデータセット
• Np$N_p$ : 親のデータセットの総数
• Dleft$D_{left}$ : 分割された左のデータセット
• Nleft$N_{left}$ : 分割された左のデータセットの総数
• Dright$D_{right}$ : 分割された右のデータセット
• Nright$N_{right}$ : 分割された右のデータセットの総数

すなわち、情報利得とは、「親の不純度」と「子の不純度の総和」の差であると言えます。

不純度として、代表的なものに以下が挙げられます。

• ジニ不純度(gini impurity) : IG=∑ci=1p(i|t)(1−p(i|t))$I_G=\sum _{i=1}^{c}p(i|t)(1-p(i|t))$
• エントロピー(entropy) : IH=−∑ci=1p(i|t)log2p(i|t)$I_H=-\sum _{i=1}^{c}p(i|t)lo{g}_{2}p(i|t)$
• 分類誤差(classification error) : IE=1−max(p(i|t))$I_E=1-max(p(i|t))$

ここでp(i|t)$p(i|t)$は、ノードt$t$においてクラスi$i$に属するサンプルの割合、c$c$はクラス数をあらわします。例えばクラス数が二つのときのあるノードtにおけるエントロピーは次のようになります。

IH=−∑i=1cp(i|t)log2p(i|t)=−(p(1|t)log2p(1|t)+p(0|t)log2p(0|t))$$I_H=-\sum _{i=1}^{c}p(i|t)lo{g}_{2}p(i|t)=-(p(1|t)lo{g}_{2}p(1|t)+p(0|t)lo{g}_{2}p(0|t))$$

ここで、t$t$におけるデータがすべてクラス1に属するとすると、エントロピーは次のようになります。

IH=−∑i=1cp(i|t)log2p(i|t)=−(1log21+0log20)=0$$I_H=-\sum _{i=1}^{c}p(i|t)lo{g}_{2}p(i|t)=-(1lo{g}_{2}1+0lo{g}_{2}0)=0$$

ただし、0log20=0$0lo{g}_{2}0=0$とします。

これより、あるノードにおいて、全ての要素が一つのクラスに属しているとき、エントロピーは最小となります。一方、各クラスが一様に分布している時エントロピーは最大となります。

これに基づいて先ほどの例の情報利得を求めてみます。

各決定木の情報利得

これより、決定木2の方が情報利得が大きい、すなわち"良い"分割であることがわかりました。

上の例では、不純度としてエントロピーを用いましたが、どれを用いるのが良いのでしょうか。それを確認するために、ニクラス分類における三つの不純度を比較してみます。図の横軸はあるノードにおける要素のクラス1に属する割合、縦軸は不純度をあらわします。

不純度の比較

上図からわかるように、ジニ不純度とエントロピーが近く、分類誤差が少し異なることがわかります。実際、先ほどの例で分類誤差を計算して見ると、面白いかもしれません。結論としては、不純度は特に悩まずジニ不純度かエントロピーを用いるのが無難だと思います。

以上より、"良い"分割方法がわかったので、情報利得が最大となるような質問を、葉が純粋になるまで繰り返すことで、決定木が得られます。

アンサンブル学習

アンサンブル学習とは、分類/回帰問題に対して、複数の予測器を組み合わせる手法を指します。アンサンブル学習の根幹にある考え方は、弱い学習器を組み合わせることにより、複雑な問題に対する強い学習器が得られる、というものです。アンサンブル学習には、代表的な手法としてバギングとブースティングが挙げられます。

バギング

訓練データからランダムに復元抽出をM回繰り返すことにより、M組の新たなデータ集合を得る手法をブートストラップ法(bootstrap)と呼びます。

復元抽出と非復元抽出

バギング(bagging)は、こうして得られたM組のデータに対して、それぞれ予測器を学習させて、予測器の重み付け和を用いて予測を行う手法です。 訓練データ集合D$D$からブートストラップ法で得られたM組のデータセットをD1$D_1$, ..., DM$D_M$とし、データ集合Di$D_i$で学習した予測器をhi$h_i$とします。 このとき、回帰問題では次のように予測を行います。

H(x)=1M∑i=1Mhi(x)$$H(x)=\frac{1}{M}\sum _{i=1}^M{h}_i(x)$$

分類問題では最多票を得たクラスを分類結果とします。

H(x)=argmax∑i=1MI(hi(x)=t)$$H(x)=argmax\sum _{i=1}^{M}I(h_i(x)=t)$$

ブースティング

バギングは各予測器を並列に学習できるという利点がありますが、独立に予測器を学習させるため多様性が生まれにくいという欠点があります。そこで、独立に学習するのではなく、逐次的に予測器を学習させる方法が提案されました。これをブースティング(boosting)と呼びます。すなわち、前段の分類器hi(x)$h_i(x)$で誤分類されたデータに大きな重みが与えられ、後段の分類器hi+1(x)$h_{i+1}(x)$の学習に利用されます。また、各分類器では、誤り率ϵi${ϵ}_i$が計算され、それに応じた信頼度αi${\alpha }_i$が計算されます。この信頼度により、重み付けされた各学習器の和を用いて分類を行います。

勾配ブースティング決定木

ここまでの話から、勾配ブースティング決定木は、決定木を用いてアンサンブル学習、特にブースティングを用いる手法だということがわかると思います。では、どのように前段の決定木から後段の決定木を生成するのでしょうか。

そこで、まずは次の目的関数を考えます。

obj(θ)=L(θ)+Ω(θ)$$obj(\theta )=L(\theta )+\mathrm{\Omega }(\theta )$$

ここで、L$L$はロス関数、Ω$\mathrm{\Omega }$は正則化項をあらわします。ロス関数の例としては、次のようなものが挙げられます。

• 二乗誤差(squared error): L(θ)=∑i(yi−y^i)2$L(\theta )=\sum _i(y_i-{\hat{y}}_i{)}^2$
• ロジスティックロス(logistic loss): L(θ)=∑i[(yiln(1+e−y^i)+(1−yi)ln(1+ey^i)]$L(\theta )=\sum _i[(y_{i}ln(1+e^{-{\hat{y}}_i})+(1-y_i)ln(1+e^{{\hat{y}}_i})]$

ただし、yi$y_i$は真の値、y^i${\hat{y}}_i$は予測値をあらわします。

いま、K個の決定木fk$f_k$を用いて予測を行うとき、i$i$番目のデータxi$x_i$に対する予測値は次のようになります。

y^i=∑k=1Kfk(xi)$${\hat{y}}_i=\sum _{k=1}^K{f}_k(x_i)$$

このとき、目的関数はロス関数をl$l$とすると、次のようにあらわされます。

obj(θ)=∑i=1nl(yi,y^i)+∑k=1KΩ(fk)$$obj(\theta )=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i)+\sum _{k=1}^K\mathrm{\Omega }(f_k)$$

各反復での訓練

t$t$番目の決定木ft$f_t$を求めるために、t$t$番目の予測値y^(t)i${\hat{y}}_i^{(t)}$について考えます。

y^(0)i=0y^(1)i=f1(xi)=y^(0)i+f1(xi)... y^(t)i=∑k=1tfk(xi)=y^(t−1)i+ft(xi)$$\begin{gathered} {\hat{y}}_i^{(0)}=0 \\ {\hat{y}}_i^{(1)}=f_1(x_i)={\hat{y}}_i^{(0)}+f_1(x_i) \\ ... \\ {\hat{y}}_i^{(t)}=\sum _{k=1}^t{f}_k(x_i)={\hat{y}}_i^{(t-1)}+f_t(x_i) \end{gathered}$$

これより、t$t$番目の目的関数は次のようになります。

obj(t)=∑i=1nl(yi,y^(t−1)i+ft(xi))+Ω(ft)+Const$$ob{j}^{(t)}=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i^{(t-1)}+f_t(x_i))+\mathrm{\Omega }(f_t)+Const$$

ここで、t−1$t-1$番目までの決定木はfixなので、正則化項はt$t$番目の決定木のみ考えています。 ここで、ロス関数をテイラー展開し、二次の項までで打ち切った目的関数を考えます。

obj(t)=∑i=1n[l(yi,y^(t−1)i)+gift(xi)+12hif2t(xi)]+Ω(ft)+Const$$ob{j}^{(t)}=\sum _{i=1}^n\left[l(y_i,{\hat{y}}_i^{(t-1)})+g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}_t^2(x_i)\right]+\mathrm{\Omega }(f_t)+Const$$

ただし、gi$g_i$、hi$h_i$はそれぞれロス関数のy^(t−1)i${\hat{y}}_i^{(t-1)}$における一階微分と二階微分をあらわします。

以上より、t$t$番目の目的関数の定数部分を除いた式は、次のとおりです。

∑i=1n[gift(xi)+12hif2t(xi)]+Ω(ft)$$\sum _{i=1}^n\left[g_i{f}_t(x_i)+\frac{1}{2}{h}_i{f}_t^2(x_i)\right]+\mathrm{\Omega }(f_t)$$

これを最小化するようにt$t$番目の決定木f(t)$f^{(t)}$を構築します。 上式より、gi$g_i$、hi$h_i$が与えられれば、f(t)$f^{(t)}$を求めることができ、これがまさにXGBoostがサポートしているロス関数の条件となります。

モデルの複雑性

次に正則化項の定義を与えます。正則化項の役割としては、モデルが複雑になりすぎることを防ぐ、すなわち汎化性能を高めることにあります。

まず、葉の数がT$T$個であると仮定します。ここでx$x$が与えられたとき、それをどの葉に割り当てるかを決定する関数をq:Rd→$q:R^d\to${1,2,...,T$1,2,...,T$}とします。さらに、葉のスコアのベクトルをw$w$とすると、t$t$番目の決定木のスコアは次のようにあらわされます。

ft(x)=wq(x)$$f_t(x)=w_{q(x)}$$

XGBoostでは、これを用いて正則化項を次のように定義します。

Ω(f)=γT+12λ∑j=1Twj$$\mathrm{\Omega }(f)=\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j$$

これにより、葉の数(木構造の深さ)T$T$と各ノードの値w$w$が過剰に大きくなりすぎることにペナルティを与えることができます。

最適解の導出

ロス関数と正則化項より、目的関数は次のようにあらわされます。(Const$Const$は無視)

Obj(t)=∑i=1n[giwq(xi)+12hiw2q(xi)]+γT+12λ∑j=1Twj$$Ob{j}^{(t)}=\sum _{i=1}^n\left[g_i{w}_{q(x_i)}+\frac{1}{2}{h}_i{w}_{q(x_i)}^2\right]+\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j$$

これを整理すると、

∑j=1T⎡⎣(∑i∈Ijgi)wj+12(∑i∈Ijhi+λ)w2j⎤⎦+γT$$\sum _{j=1}^T\left[(\sum _{i\in I_j}{g}_i)w_j+\frac{1}{2}(\sum _{i\in I_j}{h}_i+\lambda )w_j^2\right]+\gamma T$$

ただし、Ij=$I_j=${i|q(xi)=j$i|q(x_i)=j$}となるインデックス集合です。 さらに、

Gj=∑i∈Ijgi,Hj=∑i∈Ijhi$$G_j=\sum _{i\in I_j}{g}_i,H_j=\sum _{i\in I_j}{h}_i$$

とおくと、以下のようになります。

Obj(t)=∑j=1T[Gjwj+12(Hj+λ)w2j]+γT$$Ob{j}^{(t)}=\sum _{j=1}^T\left[G_j{w}_j+\frac{1}{2}(H_j+\lambda )w_j^2\right]+\gamma T$$

ここで、wj$w_j$は独立であり、二次関数の形から最適解とそのときの値は次のように与えられます。

w∗j=−GjHj+λ$$w_j^{\ast }=-\frac{G_j}{H_j+\lambda }$$

Obj∗=−12∑j=1TG2jHj+λ+γT$$Ob{j}^{\ast }=-\frac{1}{2}\sum _{j=1}^T\frac{G_j^2}{H_j+\lambda }+\gamma T$$

それぞれの値のイメージは以下となります。

決定木の構造(引用[3])

実装例

以上のアルゴリズムを実装した例がこちら( GitHub - nyk510/gradient-boosted-decision-tree: 勾配ブースティングのpythonによる実装 )にあるので、実行してみます。実装もわかりやすいので、アルゴリズムを細かく追いたい人は、オススメです。

分類問題の訓練データは、[1,1]$[1,1]$、[−1,−1]$[-1,-1]$を中心としたガウス分布からのサンプリングとなっています。

分類問題に対する予測

回帰問題の訓練データはsin関数にガウス分布からのノイズを付与したものとなります。

回帰問題に対する予測

分類、回帰ともに訓練データに対する予測ができていることがわかります。

まとめ

少し飛ばし気味のところもありますが、勾配ブースティング決定木について見てきました。より詳しく知りたい方は、参考の方をどうぞ。

この記事を書くきっかけとなったのが、現在開催中のkaggleの某コンペなので、そちらも参加される人は頑張りましょう。

不定期更新ですがよろしくお願いします。 何かあればコメントやツイッター(@hiyoko9t)でお気軽にお声がけください。

参考

[1] Sebastian Raschka, "Python 機械学習プログラミング", インプレス, 2016.
[2] nyk510, GitHub - nyk510/gradient-boosted-decision-tree: 勾配ブースティングのpythonによる実装
[3] Introduction to Boosted Trees — xgboost 0.6 documentation, 2015.