ともにゃん的データ分析ブログ

勉強したことの備忘録とかね

XGBoostの概要

XGBoostの凄さに最近気がついたので、もうちょっと詳しく知りたいと思って以下の論文を読みました。

XGBoost: A Scalable Tree Boosting System

せっかくなので、簡単にまとめてみたいと思います。。。と思っていたら結構な量になってしいました。
何か間違い等がありましたらコメントをしていただくか、@kefism へ連絡をしてくださると嬉しいです。


XGBoostとは

基本的に内部で行われていることは決定木を複数個作るということです。しかし、その作り方に特徴があります。

ここで記号を定義しておきましょう。
説明変数をm次元としてx_{i}=(x_{i1},\ldots, x_{im})^{T}、目的変数をy_{i}i=1,\ldots,nとします。ここでnはデータの数です。またデータx_{i}に対する予測値を\hat{y}_{i}とします。

まず決定木を1つ構築します。するとその決定木を使用して予測ができるようになります。1つ目の決定木から得られるデータx_{i}に対する予測値を\hat{y}_{i}^{(1)}としましょう。このときに実測値y_{i}と予測値\hat{y}_{i}^{(1)}の誤差は

 \epsilon_{i}^{(1)} = y_{i} - \hat{y}_{i}^{(1)}

となります。上記の誤差は何を表しているかというと、1つ目の決定木では学習することができなかった部分、つまり説明することができなかった部分になります。y_{i} - \hat{y}_{i}^{(1)}が0に近い値であれば実測と予測がほほ一致しているわけですから学習ができている、つまりi番目のデータについては1つの決定木で説明することができたということになりますし、0から離れた値になっていれば予測が外れているわけですから学習がうまくいっていないことになります。

さて、もしも誤差が予測できたらどうでしょうか。実測値 - 予測値 = 誤差 なわけですから、誤差がわかれば 実測値 = 予測値 + 誤差 ということで、予測値に誤差を足してあげれば実測値になります。そこでXGBoostでは、1つ目の決定木では予測できなかった部分、つまり誤差を目的変数として、2つ目の決定木を構築します。誤差を予測するモデルを作るわけです。

2つ目の決定木を構築したとします。するとその決定木を使用して誤差の予測ができるようになります。2つ目の決定木から得られるデータx_{i}に対する予測値を\hat{\epsilon}_{i}^{(2)}としましょう。つまり\hat{\epsilon}_{i}^{(2)}は、1つ目の決定木によって説明できなかった部分\epsilon_{i}^{(1)}の予測値です。このときに実測値y_{i}に対する予測値は、1つ目の決定木の予測値に2つ目の決定木の予測値、つまり予測された誤差を加えて次のようになります:

\hat{y}_{i}^{(2)} = \hat{y}_{i}^{(1)} + \hat{\epsilon}_{i}^{(2)}

したがって、2つ目の決定木を構築し終えた時点での実測値と予測値の誤差は


\begin{eqnarray}
\epsilon_{i}^{(2)} &=& y_{i} - \hat{y}_{i}^{(2)} \\
&=& y_{i} - (\hat{y}_{i}^{(1)} + \hat{\epsilon}_{i}^{(2)})
\end{eqnarray}

となります。ここで \epsilon_{i}^{(2)}は、2つ目の決定木を構築してもなお学習できなかった部分です。ということで、さらに \epsilon_{i}^{(2)}を目的変数として3つ目の決定木を構築します。繰り返しになりますが、3つ目の決定木を構築し終えた後のデータx_{i}に対する予測値は、1つ目の決定木で予測された値 + 1つ目の決定木では学習できなかった誤差の予測値(2つ目の決定木の結果) + 2つの決定木を構築してもなお学習できなかった誤差の予測値(3つめの決定木の結果)ということで、次のようになります:


\begin{eqnarray}
\hat{y}_{i}^{(3)} &=& \hat{y}_{i}^{(2)} + \hat{\epsilon}_{i}^{(3)} \\
&=& \hat{y}_{i}^{(1)} + \hat{\epsilon}_{i}^{(2)} + \hat{\epsilon}_{i}^{(3)}
\end{eqnarray}

以下、これと同じ操作をK回繰り返します。この繰り返す回数Kは分析者が決めるパラメータとなり、クロスバリデーション等でうまく決定してあげる必要があります。


このように、今までに学習したモデルの情報を使って、新たなモデルを構築することでデータの学習を進める方法をブースティング(Boosting)と呼びます。

この記事の冒頭で、XGBoostの内部で行われていることは決定木を複数個作るということと言いました。より正確にはブースティングする際に、つまり1つ前までの決定木の結果を利用して新たな決定木を作る際に、実測値と予測値との誤差がある意味で最小になるように決定木のアルゴリズムを構築したものがXGBoostです。ここで"ある意味で"といったのは、世の機械学習アルゴリズムは目的関数を最小(または最大)にするようにアルゴリズムを構築するわけですが、まだXGBoostにおける目的関数の説明をしていないからです。

ここで強調しなければならないことは、単純に誤差を目的変数にして新たな決定木を構築するのではなく、前回までの決定木の情報を利用して新たな決定木を構築するということです。というわけで、新たな決定木を構築する際に前回までの情報をどのように利用すればよいのかという話に繋がっていきます。

決定木は、説明変数を分割していって、わかりやすく言えば樹形図を構築するアルゴリズムです。このことから前回までの情報を決定木を構築する際に利用するということは、前回までの情報を利用して説明変数の分割を最適化するということなのです。そして、後で説明しますが、その前回までの情報というのは、目的関数の勾配(Gradient)として取り出すことができます。XGBoostのGはGradientのGです。


XGBoostの定式化

さて、これから説明しなければならないことは、「どのようにして前回までの情報を使った決定木を構築するか」です。そこで今はひとまずk個目の決定木そのものを抽象的に関数f_{k}とします。関数 f にデータ x_{i}を渡すと、そのデータに対する予測値を返してくれます。今までの例だと、f_{1}(x_{i})=\hat{y}_{i}^{(1)}, f_{2}(x_{i})=\hat{\epsilon}_{i}^{(2)}, f_{3}(x_{i})=\hat{\epsilon}_{i}^{(3)}です。そしてブースティングをK回行ったときの、つまりK個目の決定木を構築したときの予測値は次のようになります:

\hat{y}_{i}^{(K)} = \sum_{k=1}^{K}f_{k}(x_{i})

まだfがどのような形になっているか、つまり「前回までの情報を使った決定木はどういったアルゴリズムか」は分かりません。そこで、色々と計算をすることで、このfがどういった性質を満たしていれば良いのかというのを見ていくことにします。

今までの話を数式として表現しましょう。

損失関数 l(a,b) を導入します。損失関数として、例えば回帰の際は多くの場合で二乗誤差

二乗誤差 = (実測値 - 予測値)^{2}

が用いられますね。今回の説明では損失関数を二乗誤差と設定するのではなく、より一般的な(抽象的な)形で議論を進めていきます。損失関数は微分可能で、a(例えば、実測値)とb(例えば、予測値)が近ければ近いほど小さい値をとり、遠ければ遠いほど大きな値をとるように設定されます。したがって、今回の「前までのt-1個の決定木の情報を利用して、t個目の決定木をどのように構築すれば良いのか」という問題は次のように定式化されます:


\begin{eqnarray}
&& \min_{f_{t}}\ \sum_{i=1}^{n}l(y_{i}, \hat{y}_{i}^{(k)}) \\
&\iff& \\
&& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \sum_{k=1}^{t}f_{k}(x_{i})\right) \\
&\iff& \\
&& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \sum_{k=1}^{t-1}f_{k}(x_{i}) + f_{t}(x_{i})\right) \\
&\iff& \\
&& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \hat{y}_{i}^{(t-1)} + f_{t}(x_{i})\right) \\
\end{eqnarray}

この上記の式を最小にする決定木 f_{t} を構築すればいいということになります。

しかし、このままだと学習に用いたデータに過剰にフィッティングしすぎることによって、予測性能が下がる過学習(over fitting)という現象が起きるでしょう。とくに上記の定式化のままではこれが顕著に現れます。なぜなら、実際に手元に手に入っている y_{i},\ i=1,\ldots,n と予測の誤差が小さくなることだけを今は考えていて、未来のデータに対する予測のことを一切考慮していないからです。

そこで、次のような罰則項を定義します:

\Omega(f_{t})=\gamma T + \frac{1}{2} \lambda ||w||^{2}

そして、上記の罰則項を損失関数に加えたものを \mathcal{L}^{(t)}(f_{t}) とし、最終的に解きたい問題を以下のように設定します:


\begin{eqnarray}
\min_{f_{t}} \mathcal{L}^{(t)}(f_{t}) &=& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \hat{y}_{i}^{(t-1)} + f_{t}(x_{i})\right) + \Omega(f_{t})\\
&=& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \hat{y}_{i}^{(t-1)} + f_{t}(x_{i})\right) + \gamma T + \frac{1}{2} \lambda ||w||^{2}
\end{eqnarray}

ここでTは決定木を構築したときの最終ノードの数(つまり、木の大きさ)です。木の大きさTを大きく設定すればするほど、過学習を起こしやすくなるのでlの部分は小さくなるので目的関数 \mathcal{L}(f_{t}) の値を小さくするように働きますが、\gamma Tの部分は \mathcal{L}(f_{t})が大きくなるように働きます。結果として、あまりTの大きな木は好まれなくなります。\gammaTの大きさに対するペナルティーで、\gamma=0のときは過学習なんてどうでもいいから、とにかくlの部分が小さくできればいい!!ということを表していて、\gamma > 0を設定してやることで、\gammaの値が大きければ大きいほどTの値が小さい、つまり小さな木が好まれるようになります。したがって \gamma は分析者が事前に決めてやる必要のあるパラメータです。

また、w は決定木 f_{t} が返すことのできる値のベクトルです。例えば下図の木の場合はT=4で、返すことのできる値はw=(w_{1},w_{2},w_{3},w_{4})の4つということになります。つまり、w\hat{y}_{i}^{(t)} がとることのできる値の集合ということになります。

f:id:kefits:20170610234136p:plain:w300

\lambda はこの決定木f_{t} が返すことのできる値wの大きさ調整するパラメータです。大きな\lambdaを設定すると、小さなwが好まれるようになります。t個目の決定木を構築し終えた後の予測値は、t-1個目までの結果の和にt個目の結果を加えた\sum_{k=1}^{t-1}f_{k}(x_{i}) + f_{t}(x_{i})ですから、 f_{t}(x_{i}) の値が小さいということは、実測値と予測値の差をちびちびと更新していくような感じになります。一方で小さな\lambdaを設定すると、実測値と予測値の誤差が早く小さくなるように大胆に更新をおこなっていきます。つまり、\lambda は新たな決定木を構築したときの更新の幅を調整するパラメータになります。ちびちびと更新すればKの数が大きくなる、つまりたくさんの木を構築することになるので学習に時間がかかりますが、その分様々なパターンを表現できる一方、過学習の要因ともなります。大胆に更新をすれば早く実測値と予測値の差が小さくなっていきますが、早く差が小さくなるということは、最終的に構築される決定木の数Kが少なくなるということなので入力データの色々なパターンをとらえられません。バランスが大切ということで、\lambdaもクロスバリデーション等で前もって分析者が決めてやる必要があります。

以上でXGBoostの定式化が完了しました。次からは、具体的に

「前までのt-1個の結果を使って、\mathcal{L}^{(t)}(f_{t})を最小にするようにt個目の決定木f_{t}をどのように構築すればいいのか?」

というところを見ていきます。


XGBoostのアルゴリズムの導出

t 個目の決定木が返すべき値

目的関数\mathcal{L}^{(t)}(f_{t})の第1項、lf_{t}に関して0の周りで2次のテーラー展開を行います。


\begin{eqnarray}
\min_{f_{t}} \mathcal{L}^{(t)}(f_{t}) &=& \min_{f_{t}}\ \sum_{i=1}^{n}l\left(y_{i}, \hat{y}_{i}^{(t-1)} + f_{t}(x_{i})\right) + \Omega(f_{t})\\
&\approx& \min_{f_{t}}\ \sum_{i=1}^{n}\left( l\left(y_{i}, \hat{y}_{i}^{(t-1)}\right) +g_{i}f_{t}(x_{i}) + \frac{1}{2}h_{i}f_{t}^{2}(x_{i}) \right)+ \Omega(f_{t})
\end{eqnarray}

ここで g_{i}=\partial_{\hat{y}_{i}^{(t-1)}}\ l\left(y_{i}, \hat{y}_{i}^{(t-1)}\right), h_{i}=\partial_{\hat{y}_{i}^{(t-1)}}^{2}\ l\left(y_{i}, \hat{y}_{i}^{(t-1)}\right) です。

最適化に関係のない項、つまりf_{t}が関わっていない項を取り除いたものを \tilde{\mathcal{L}}^{(t)} とすると


\begin{eqnarray}
\min_{f_{t}} \tilde{\mathcal{L}}^{(t)}(f_{t}) = \min_{f_{t}}\ \sum_{i=1}^{n}\left(g_{i}f_{t}(x_{i}) + \frac{1}{2}h_{i}f_{t}^{2}(x_{i}) \right)+ \Omega(f_{t})
\end{eqnarray}

となります。
ここで、最終ノードj にあるデータの集合を I_{j} とします。どういう事かというと、例えば I_{4}

f:id:kefits:20170611132136p:plain:w300

上記の図の「ノード4」と書かれたノードの中にあるデータ x_{i} たちのことです。図では添字の i を省略しました。

このときに、さらに \tilde{\mathcal{L}}^{(t)} を以下のように書き直します:


\begin{eqnarray}
\min_{f_{t}} \tilde{\mathcal{L}}^{(t)}(f_{t}) &=& \min_{f_{t}}\ \sum_{i=1}^{n}\left(g_{i}f_{t}(x_{i}) + \frac{1}{2}h_{i}f_{t}^{2}(x_{i}) \right)+ \Omega(f_{t})\\
&=& \min_{f_{t}}\ \sum_{i=1}^{n}\left(g_{i}f_{t}(x_{i}) + \frac{1}{2}h_{i}f_{t}^{2}(x_{i}) \right)+ \gamma T + \frac{1}{2} \lambda ||w||^{2}\\
&=& \min_{f_{t}}\ \sum_{j=1}^{T} \left( \sum_{i \in I_{j}} g_{i}f_{t}(x_{i}) + \frac{1}{2}\sum_{i \in I_{j}}h_{i}f_{t}^{2}(x_{i}) \right)+ \gamma T + \frac{1}{2} \lambda \sum_{k=1}^{T}w_{j}^{2}\\
&=& \min_{f_{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}w_{j}^{2} \right)+ \gamma T + \frac{1}{2} \lambda \sum_{k=1}^{T}w_{j}^{2}\\
&=& \min_{f_{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}w_{j}^{2} + \frac{1}{2}\lambda w_{j}^{2} \right)+ \gamma T\\
&=& \min_{f_{t}}\ \sum_{j=1}^{T} \left( \sum_{i \in I_{j}} g_{i}w_{j} + \frac{1}{2}\left(\sum_{i \in I_{j}}h_{i} + \lambda \right) w_{j}^{2} \right)+ \gamma T
\end{eqnarray}

ここで、2行目から3行目にかけて\sum_{i=1}^{n}\sum_{j=1}^{T}\sum_{i \in I_{j}} となったのは、最終ノードごとについて足し上げてた後に、それらをさらに足し上げた結果と、最初から全データ(がごちゃまぜになっているときに)ついて一度に足し上げたときの結果は同じであることからわかると思います。また、同じく2行目から3行目にかけて ||w||^2 をきちんと書き下しています。

さらに、3行目から4行目にかけてf_{t}(x_{i})w_{j} となっているのは、データをノード j にあるデータのみに制限したとき、つまり、例えばデータx_{i}がノード j に来たときに決定木 f_{t} が返す結果は f_{t}(x_{i}) = w_{j} であることを思い出すと理解できると思います。

あとは共通項でくくれるものはくくって、目的関数の式変形はおしまいです。


さて、目的関数を2次形式で書くことができました。したがって、 \tilde{\mathcal{L}}^{(t)}w_{j}微分したものを0とおいて解くと最適解 w_{j}^{\ast} は、つまり目的関数を最小にする t 個目の決定木の最終ノード j が返すべき結果


\begin{eqnarray}
w_{j}^{\ast} = -\frac{\sum_{i \in I_{j}} g_{i}}{\sum_{i \in I_{j}} h_{i} + \lambda}
\end{eqnarray}

となります。

ここで g_{i}=\partial_{\hat{y}_{i}^{(t-1)}}\ l\left(y_{i}, \hat{y}_{i}^{(t-1)}\right), h_{i}=\partial_{\hat{y}_{i}^{(t-1)}}^{2}\ l\left(y_{i}, \hat{y}_{i}^{(t-1)}\right) であったことを思い出します。これらはそれぞれ損失関数の1次、2次の勾配であり、また前までの t-1 個の決定木の結果と手元にある学習データの実測値 y_{i} を使って計算できる値です。t-1 個の決定木の情報が勾配という形で t 個目の決定木を構築する際に利用されています。これがXGBoostです。


勾配情報を利用した木の分岐方法

今までの話は、データが最終ノードに辿り着いた後の話でした。どのように説明変数を分岐させていけばよいのかはまだ説明していません。ですのでこの節では、勾配情報をどのように利用して説明変数を分岐していけば良いのかという話をします。

さて、前節で求めた w_{j}^{\ast}\tilde{\mathcal{L}}^{(t)} に代入すると \tilde{\mathcal{L}}^{(t)} は次のようになります:


\begin{eqnarray}
\tilde{\mathcal{L}}^{(t)} = -\frac{1}{2}\sum_{j=1}^{T} \frac{\left( \sum_{i \in I_{j}} g_{i} \right)^{2}}{\sum_{i \in I_{j}} h_{i} + \lambda} + \gamma T
\end{eqnarray}

この式を眺めてみると、各最終ノードにおける


\frac{\left( \sum_{i \in I_{j}} g_{i} \right)^{2}}{\sum_{i \in I_{j}} h_{i} + \lambda}

の値がが大きければ大きいほど、目的関数 \tilde{\mathcal{L}}^{(t)} の値はより小さくなります。つまり、あるノードを左(L)と右(R)に分岐するときに、最も目的関数 \tilde{\mathcal{L}}^{(t)} の値が小さくなるように分岐をするためには、分岐前と分岐後の \frac{\left( \sum_{i \in I_{j}} g_{i} \right)^{2}}{\sum_{i \in I_{j}} h_{i} + \lambda}の値を比べてあげればよいということになります。

どういった比較をすれば良いのか具体的にみていきましょう。下図のように現時点で T 個のノードができていて、あるノード s を更に分割することを考えます。

f:id:kefits:20170611173519p:plain:w300

ノード s を分岐する前の目的関数を


\tilde{\mathcal{L}}_{\mbox{before}}^{(t)} = -\frac{1}{2}\sum_{j \neq s}^{T} \frac{\left( \sum_{i \in I_{j}} g_{i} \right)^{2}}{\sum_{i \in I_{j}} h_{i} + \lambda} - \frac{1}{2}\frac{\left( \sum_{i \in I_{s}} g_{i} \right)^{2}}{\sum_{i \in I_{s}} h_{i} + \lambda} + \gamma T

ノード r を分岐した後に左側のノードに入るデータの集合を I_{L}, 右側のノードに入るデータの集合を I_{R}としたときの (I_{s} = I_{L} \cup I_{R}) 分岐後の目的関数を


\tilde{\mathcal{L}}_{\mbox{after}}^{(t)} = -\frac{1}{2}\sum_{j \neq s}^{T} \frac{\left( \sum_{i \in I_{j}} g_{i} \right)^{2}}{\sum_{i \in I_{j}} h_{i} + \lambda} - \frac{1}{2}\frac{\left( \sum_{i \in I_{L}} g_{i} \right)^{2}}{\sum_{i \in I_{L}} h_{i} + \lambda} - \frac{1}{2}\frac{\left( \sum_{i \in I_{R}} g_{i} \right)^{2}}{\sum_{i \in I_{R}} h_{i} + \lambda} + \gamma T

とします。最も目的関数 \tilde{\mathcal{L}}^{(t)} の値が小さくなるように分岐をするためには、分岐前と分岐後の目的関数の差をとって、もっとも差が大きくなるところで分岐を行えば良いということになります。つまり、


\begin{eqnarray}
\tilde{\mathcal{L}}_{\mbox{before}}^{(t)} - \tilde{\mathcal{L}}_{\mbox{after}}^{(t)} &=& \frac{1}{2}\frac{\left( \sum_{i \in I_{L}} g_{i} \right)^{2}}{\sum_{i \in I_{L}} h_{i} + \lambda} + \frac{1}{2}\frac{\left( \sum_{i \in I_{R}} g_{i} \right)^{2}}{\sum_{i \in I_{R}} h_{i} + \lambda} - \frac{1}{2}\frac{\left( \sum_{i \in I_{s}} g_{i} \right)^{2}}{\sum_{i \in I_{s}} h_{i} + \lambda}\\
&=& \frac{1}{2}\left(\frac{\left( \sum_{i \in I_{L}} g_{i} \right)^{2}}{\sum_{i \in I_{L}} h_{i} + \lambda} + \frac{\left( \sum_{i \in I_{R}} g_{i} \right)^{2}}{\sum_{i \in I_{R}} h_{i} + \lambda} - \frac{\left( \sum_{i \in I_{s}} g_{i} \right)^{2}}{\sum_{i \in I_{s}} h_{i} + \lambda} \right)
\end{eqnarray}

が最大になるような分岐をおこなえば良いです。より詳しく言うと、すべての説明変数 x_{i}=(x_{i1},\ldots, x_{im})^{T} のあらゆる分割の候補のうち、最も\tilde{\mathcal{L}}_{\mbox{before}}^{(t)} - \tilde{\mathcal{L}}_{\mbox{after}}^{(t)} を大きくする候補を1つ選んで、そこで分岐を行うということです。

このことをアルゴリズムとしてまとめると、次の図のようになります(論文より引用)。

f:id:kefits:20170611173443p:plain:w400

しかし、すべての説明変数のあらゆる分割の候補を考えることは時間の制約もあり非常に大変なことです。そこで上図のアルゴリズムを近似したものが提案されており、精度としても問題なく使えるアルゴリズムとなっているようです。その近似アルゴリズムは次の図のようになります(論文より引用)。

f:id:kefits:20170611174245p:plain:w400

今回はこのアルゴリズムを詳しく説明することはありません。是非論文を読んでみてください。


以上でXGBoostのアルゴリズムの説明は終了になります。お疲れ様でした。


XGBoostで学習するときのテクニック

XGBoostでは、過学習を防ぐために正則化以外の様々なテクニックが利用されています。

Shrinkage

XGBoostにおいて、ある決定木から返されるノード j 結果は w_{j} でした。この w_{j}0 < \eta \leq 1 の値をかけた \eta w_{j} を返すようにすることを返り値のShrinkageと呼びます。パラメータ \lambda と似たような働きがありますが、\lambda は学習時に決定木が返す値の大きさを調整するパラメータでしたが、 \eta は学習後に直接返り値 w_{j} を修正するパラメータになります。

特徴量のsubsampling

ランダムフォレストでは、複数の決定木を作る際にすべての説明変数を使用するのではなく、ランダムに決められた割合で説明変数の数を選んで決定木を構築します。XGBoostでも同様に、ブースティングをする際の決定木を構築するときに、すべての説明変数を使用するのではなく、ランダムに決められた割合で説明変数の数を選んで決定木を構築します。特徴量のsubsamplingと呼ばれます。1つの決定木を構築する際に使用する特徴量の数が減るので、アルゴリズムが回るスピードも早くなります。

データのsubsampling

ブースティングをする際に n 個すべてのデータを使用して決定木を構築するのではなく、ランダムに決められた割合でデータの数を選んで決定木を構築します。これをデータのsubsamplingと呼びます。1つの決定木を構築する際に使用するデータの数が減るので、アルゴリズムが回るスピードも早くなります。


最後に

簡単にまとめるつもりだったのに大作になってしまった。

疲れた