Edited at

Gradient Boosted Tree (Xgboost) の取り扱い説明書

勾配ブースティング木は性能が高くてすごい、ということしかわかっていなかったので、ブースティングだったり木構造の学習アルゴリズムだったりの勉強を兼ねて内容をまとめてみました。


Introduction to Boosted Trees

http://xgboost.readthedocs.io/en/latest/model.html



Notationについて


  • n番目の入力データ $x_n$, その教師ラベル $y_n$

  • $T$ 番目の学習段階においての予測モデルの出力値 $\tilde{y}^{(T)}=\sum_{\tau=1}^T f^{(\tau)}$


考えていく問題


一般的な機械学習の問題設定

通常の機械学習のように入力 $x$ について、予測モデル $f(x)$ をいかにして設計するかということを考えます。


入力データ $x_n$ に対する正しい値を $y_n$ とおき、予測モデルの出力 $f(x_n)$ との差分を計算する関数をロス関数 $l(y,\tilde{y})$ として以下のように定義します。

$$

l\left( y_n,f(x_n) \right)

$$

また、この予測器は複数の関数の和として表されている(加法的モデル)であるとします。


加法的モデルを式にすると $f(x)=\sum_{i=1}^t f^{(i)}(x)$ となります。このそれぞれのモデル $f^{(i)}$ については反復を繰り返して求めていきます。

データ点に対する当てはまりだけを考えているとデータに過剰にフィットしてしまうので、これに正則化 $\Omega$ を加えた

$$

{\rm Obj}=\sum_{n=1}^N l(y_n,f(x_n))+ \Omega(f)

$$

を最小化する問題を考えていきます。ここで正則化の関数 $\Omega$ は引数に関数を取り、その複雑さを計算する関数です。一般線形回帰モデル $f(x)=g(w^{\mathrm T}x)$ に対して、予測モデルの内部で用いられている重みパラメータ $w$ の関数として表すものなどが有名です。

例えば、$\Omega(f)=\|w \|_2$ とすればRidge回帰に相当し、$\Omega(f)= \| w \|_1$ ととればLassoに相当する問題となります。


各反復での学習

各学習ステップ $t$ においては、それ以前に決定した関数をすべて固定して、 $f^{(t)}$ のみについて学習を行います。ですから、直前の反復までで求められた予測関数は固定です。なので、 $t-1$ 番目の反復での出力 $\tilde{y}^{(t-1)}$ は、 $t$ の反復の時変化しません。

この $\tilde{y}^{(t)}$ と各反復で求められて予測関数に加えられる関数 $f^{(t)}$ との関係を表したのが以下の等式です。

\begin{align}

\tilde{y}^{(0)}&= 0 \\
\tilde{y}^{(1)}&= \tilde{y}^{(0)} + f^{(1)} = f^{(0)}+ f^{(1)}\\
\tilde{y}^{(2)}&= \tilde{y}^{(1)} + f^{(2)} = f^{(0)}+ f^{(1)}+ f^{(2)}\\
\vdots \\
\tilde{y}^{(t)}&= \tilde{y}^{(t-1)} + f^{(t)} = \sum_{i=1}^{t-1}f^{(i)}+ f^{(t)}
\end{align}

これを目的関数に代入すると以下のようになります

\begin{align}

{\rm Obj}^{(t)}&=\sum_{n=1}^N l(y_n,\tilde{y}^{(t)}_n)+\sum_{i=1}^t\Omega(f^{(i)}) \\
&=\sum_{n=1}^N l(y_n,\tilde{y}^{(t-1)}_n+f^{(t)})+\Omega(f^{(t)}) + {\rm Const.}
\end{align}

正則化項の $i<t$ の部分に関しては $t$ ステップでは変更しませんから定数とみなせます(定数項 Const.に対応しています)

目的関数の中で、ロス関数の部分を $\tilde{y}$ に関してテイラー展開して二次の項までの近似を行うと、以下のようになります。

\begin{align}

&l(y_n,\tilde{y}^{(t-1)}_n+f^{(t)}) \\
=& l(y_n,\tilde{y}_n^{(t-1)}) + \frac{\partial}{\partial\tilde{y}}l(y_n,\tilde{y}_n^{(t-1)})f^{(t)} +
\frac{1}{2} \frac{\partial^2}{\partial\tilde{y}^2} l(y_n,\tilde{y}_n^{(t-1)}){f^{(t)}}^2 
\end{align}

これを先の目的関数に代入し、最適化の際に関係のない定数の部分を省くと、以下の関数が得られます。

{\rm Obj}^{(t)} = \sum_{n=1}^N \left(g_n f^{(t)}(x_n) + \frac{1}{2}h_n f^{(t)}(x_n)^2 \right) + \Omega(f^{(t)}) \\

\begin{eqnarray}
\left\{
\begin{array}{l}
g_n = \frac{\partial}{\partial y}l(y_n,\tilde{y}_n^{(t-1)}) \\
h_n = \frac{\partial^2}{\partial y^2}l(y_n,\tilde{y}_n^{(t-1)})
\end{array}
\right.
\end{eqnarray}


  • 毎回学習するべき目的関数は、ロス関数の二次近似


    • 誤差関数が二乗和ロス以外の場合は近似誤差が生じる



  • 勾配 $(g_n)$ およびヘシアン $(h_n)$ を直前の予測値 $\tilde{y}^{(t-1)}$ で評価したものからなる


    • 最適化の際には出力自体は必要ではないが代わりに勾配情報とヘシアン情報が必要



  • $g_n,h_n$ 以外の部分は、毎ステップでとかなければならない関数は全く同じ

以下では簡略化のために何番目のステップかという表示を消して考えます。なので考える関数は以下のようになります。

$$

\min_f{\rm Obj}^{(t)} = \sum_{n=1}^N\left(g_n f(x_n) + \frac{1}{2}h_n f(x_n)^2\right) + \Omega(f)

$$


どのようにしてこれを解くのか

目的関数が決定したのであとはこれを解いていく作業に入ります。

問題を解くためにはまず予測器のクラスを制限しないといけません。今の定式化では予測器がどういう形をしているか、ということに対する制限がありません。すなわちどのような形の関数でも大丈夫な形式になっているため、調べるべき関数の形が膨大になり、数学的に関数を決定するのが難しくなります。

そこで、様々な関数のクラスがある中でも今回は木構造を考えていきます。


木構造はデータ点を根から枝に向かう途中で条件によって分類して、末端のノード(葉)にたどり着いたらそのノードに与えられた値を予測として返すという予測器です。これを定式化していきます。


木構造の予測器

末端のノード(葉)の数を $K$ として、ノード $k\in\{1,2,...,T\}$ の時の予測値をベクトル $\omega$ の第 $k$ 成分として $\omega_k$ とします。また入力に対してどの末端ノードに割り振るかという関数を $q(x)\in\{1,2,...T\}$ と定義すると、木による予測 $f(x)$ は

$$

f(x)=\omega_{q(x)}

$$

と書くことができます。

カジュアルに言えば木構造は


  1. 適当な分類規則にしたがってデータ点を $T$ 個のクラスタに分類して

  2. それぞれのクラスタに分類されたデータ点のすべてに対して $w_j$ の値を与えるもの

と言えます。


目的関数の定式化

これを用いて先の関数を変形すると

\begin{align}

{\rm Obj}
&= \sum_{n=1}^N \left(g_n \omega_{q(x_n)}+\frac{1}{2}h_n\omega_{q(x_n)}^2 \right) + \Omega(f) \\
&= \sum_{j=1}^T \left\{ \left( \sum_{n\in I_j}g_n\right)\omega_j + \frac{1}{2}\left( \sum_{n\in I_j}h_n\right)\omega_j^2 \right\} +\Omega(f)
\end{align}

ここで $j$ に分類されるデータ番号の集合を $I_j =\left\{ n|q(x_n)=j \right\}$ と定義しました。

この変形はデータnの和で表されていたものをクラスごとの和に書きなおしたものです。ですから勾配の部分とヘシアンの部分がそれぞれそのクラスに所属するデータの和に変化しています。

このように書いた時、ロス関数部分(目的関数のデータへの当てはまり部分)のみに注目すると、 $\omega_j(j=1,2,...,T)$ に対してのただの二次関数になっています。なので単なる二次関数の最小化問題です。


正則化項の定式化

目的関数部分が求めるべきパラメータ $w$ の二次関数になっていますが、この値 $w$ が大きな値を取れれば、一つの関数 $f^{(t)}$ だけで、より大胆に目的の値 $y$ へと近づけることが可能になります。

今回考えているのは加法モデルなので、あまり大きな $w$ を取らなくともそれが有用なパラメータであればステップが進むに連れて自然と値は目的の値 $y$ に近づいていくことが期待されますから、あまり大きな $w$ を取ってしまうとデータへと過剰にフィットすることに寄与してしまいます。なのであまり大きな値を取らないように $\frac{1}{2}\lambda \sum_{j=1}^T w_j^2$ を正則化項に加えます。

また、加法モデルに加える関数があまりに複雑な関数だと、これもまた過剰なデータへの当てはまりに寄与します。よって複雑度に対するペナルティとして、葉の数 $T$ に比例する $\gamma T$ を正則化として加えます。

以上をまとめると正則化項は

$$

\Omega(f) = \frac{1}{2} \lambda \sum_{j=1}^T w_j^2 + \gamma T

$$

となります。このうち目的変数は各ノードの値 $w$ と葉の数 $T$ です。


 目的関数の最適化

以上をまとめると以下のように変形できます。

\begin{align}

{\rm Obj}
&= \sum_{n=1}^N \left(g_n \omega_{q(x_n)}+\frac{1}{2}h_n\omega_{q(x_n)}^2 \right) + \Omega(f) \\
&= \sum_{j=1}^T \left\{ \left( \sum_{n\in I_j}g_n\right)\omega_j + \frac{1}{2}\left( \sum_{n\in I_j}h_n\right)\omega_j^2 \right\} + \frac{1}{2} \lambda \sum_{j=1}^T w_j^2 + \gamma T \\
&=\sum_{j=1}^T \left\{ \left(\sum_{n\in I_j}g_n\right)\omega_j+ \frac{1}{2}\left( \sum_{n\in I_j}h_n +\lambda \right)\omega_j^2 \right\} +\gamma T
\end{align}

これは目的変数 $w_j$ に関する二次関数なので高校数学を思い出してやれば最適解 $w^{\ast}$ 及びその時の目的関数の値 ${ \rm Obj}^{\ast}$ は $\sum_{n\in I_j}g_n=G_j,\sum_{n\in I_j}h_n=H_j$ とおくと

\begin{align}

w_j^* &= - \frac{G_j}{H_j + \lambda} \\
{\rm Obj}^*&= -\frac{1}{2}\sum_{j=1}^T \frac{G_j^2}{H_j + \lambda} +\gamma T
\end{align}

となります。これはまだ葉の数 $T$ に対して最適化が行われていないことに注意してください。言い方を変えると、木の形については何も言っていないけれど、形の決まった木に対してもっとも良い葉ノードの重みはすぐに計算できるよ、ということです。


したがって、残りの仕事は木の構造を決定する、言い方を変えれば分割する規則をどう決定するのか、という問題になります。


木の構造を決定するには


簡単に思いつく方法

「木の構造が決まれば、その最適な葉ノードの重みは容易に計算できる」ことはわかりました。なので簡単に思いつく方法として


  1. 適当に可能性のある木の構造を列挙

  2. それ全てに対して最適な目的関数値 ${\rm Obj}^*$ を計算

  3. その中で最も小さい値を取る構造を選ぶ

とういう方法が考えられます。しかしこれは現実問題としてうまく機能しません。

例えば、入力 $x$ に連続変数が存在している時を考えてみることにします。この時ノードの分割の規則を決める際、いくらでも可能性のある木を列挙できてしまいます。なのでこの方法では必ずすべての可能性を列挙できません。他の方法を考えましょう。


成長戦略

なので、一気に構造を作るのは諦めて、今得られている木からより良い木を作っていく、すなわち木をどんどんと成長させていく方法を取ります。

木が成長すると何が起こるでしょうか。

木が成長するということは、どこかの末端のノードが2つに分割されるということを意味しています。仮にそのノードが $T$ 番目のノードだったとしましょう。


そのノードには $I_T$ のデータ点が分類されています。これが $R,L$ の2つの集合へと分割されるというケースについて考えてみます。すると

\begin{align}

G_T &= \sum_{i\in I_T}g_i \\
&= \sum_{i\in I_R}g_i + \sum_{i\in I_L}g_i \\
&= G_R + G_L
\end{align}

が成り立ちます。ここで新しく $G_R,G_L$ を定義しました。これと同様にして $H_T = H_R + H_L$ も成立します。

ここで、先ほどの最適解の時の目的関数の値、最適値を振り返ってみます。木を成長させる前の最適値 ${\rm Obj^{\ast}}$ は

{\rm Obj^{\ast}}=-\frac{1}{2}\left(\sum_{j=1}^{T-1}\frac{G_j^2}{H_j+\lambda} + \frac{G_T^2}{H_T+\lambda} \right) + \gamma T

一方で $T$ 番目のノードを2つに分割し、一段階成長させた木に対する最適値 ${\rm Obj'}$ は

{\rm Obj'}=-\frac{1}{2}\left(\sum_{j=1}^{T-1}\frac{G_j^2}{H_j+\lambda} + \frac{G_R^2}{H_R+\lambda} + \frac{G_L^2}{H_L+\lambda} \right) + \gamma \left(T+1 \right)

となります。目的関数が大きく減少するように関数を選びたい、ようするに木を成長させた時に、目的関数が大きく減れば良いのです。すなわち目的関数の値の差 (Gain) を見れば良いことになります。よって

\begin{align}

{\rm Gain} &= {\rm Obj^*} - {\rm Obj'} \\
&= \frac{1}{2}\left( \frac{G_R^2}{H_R+\lambda} + \frac{G_L^2}{H_L+\lambda} - \frac{G_T^2}{H_T+\lambda} \right) - \gamma \\
&= \underbrace{ \frac{1}{2}\left( \frac{G_R^2}{H_R+\lambda} + \frac{G_L^2}{H_L+\lambda} - \frac{(G_R + G_L)^2}{H_R + H_L+\lambda} \right) }_{トレーニングデータへの当てはまりの改善度}- \underbrace{\gamma}_{木を成長させることに対するペナルティ} \\
\end{align}

のようになり


  • データへの当てはまりの部分

  • 木を成長させることに対するペナルティの部分

の2つで構成されていることがわかります。

これにはもはや分割するノード以外の情報は含まれていません( $T$ に含まれるデータ以外の情報が入っていません). したがって分割した際の GAIN の計算は並列に実行することが可能です。

また勾配、ヘシアンの計算も直前の予測値 $\tilde{y}^{(t-1)}$ での評価があればよく、それ以前 $f^{(t-2)},f^{(t-3)},...$ などのパラメータなどは必要ありませんから、計算に必要なメモリも少なくて済むことがわかります。


最終的なアルゴリズム

ざっくりと書くと以下のような感じになります。


  • すべての末端ノードに対して、それに含まれるデータ点で以下の反復を繰り返す


    • ある特徴量でソートをする

    • linear scanをかけて最善の分割点を見つける

    • Gainの値を記録し、これをすべての特徴量に対して行う

    • すべての特徴の中で、最もGainが大きくなった特徴を選択する



  • 木が成長しきったらそれを加法モデルに加える

  • 以上を繰り返す

またこれに加えて


  • 各反復で得られる木をそのまま加法モデル加えず、0から1の間を取る値をかけてからモデルに加える

  • 各木の深さ(場合分けの数)が一定の数 max_depth に達したらそこで成長をやめる

  • 一定の確率で木を縮小させる(枝刈)

などのテクニックも使うことでよりオーバフィットしにくい予測モデルを生成しています。


まとめ


  • 解いている問題は通常の機械学習の定式化と同様「ロス関数+正則化項」で表される

  • 予測関数は「加法的モデル」


    • ひとつひとつのモデルが「木」であわらされる → 勾配ブースティング”木”

    • 毎回の反復では「それ以前の出力」+「今のモデル」ができるだけ正解ラベルに近づくように学習

    • 「それ以前の出力」は固定



  • 実際に解く際にはロス関数を二次近似


    • 二乗ロス以外のロス関数では近似誤差が生じる

    • 最適化の文脈で言えば逐次二次計画に近い戦略



  • 正則化項は重みの二乗和と木の深さ(場合分けの数)

  • 木を成長させた時にもっとも目的関数が小さくなるように分割の条件を探索

  • 木構造および boosting の性質から計算量的に有利


思ったこと


  • アルゴリズムにおいて入力するのは正解の値と $G_i,H_i$ の2つ


    • ロス関数の勾配とヘシアンを計算さえすればどんなモデルでも適応可能

    • xgboostではユーザーが勾配とヘシアンを定義して任意のロス関数に対して学習を行うことが可能



  • 勾配ブースティングにおける正則化は目的関数自体を変更させている


    • 早期打ち切りとか枝刈りなどの、木構造での分類における過学習を抑える手法とは根本的に発想が違う

    • 普通の機械学習と同じ発想の正則化



  • 全体がひとつの木で構成されない


    • 予測器がなりうる可能性が単一の木に比べ広がっている

    • 単一の木の平均を使う手法と比べると、探索できる空間が異なっている

    • 予測性能が高いのもそのためかもしれない



  • 正則化の与え方が自然


    • それに加えて枝狩り等の早期打ち切り策による過学習抑制も使えるために、汎化性も良い?



  • 深さに対する正則化 $\gamma T$ はパラメータ数に対するペナルティに対応する?


    • 木が進行するに連れて、分類Tの数は増えていくことが予想される

    • Tが増える=使う変数の数が増えるということ


    • max_depth=9に設定すればそれは最大で9変数に関する分割をするということ(同じ変数で分割が入るかもしれないのであくまで最大で)

    • ステップが進んでいけば、段々とすべての変数がにたようなばらつきになっていくはずで、そうなるとmat_depthに設定した変数の数だけの特徴量を使って分類をするようになっていくはず