XGBoostの論文を読んだのでGBDTについてまとめた

はじめに

今更ですが、XGboostの論文を読んだので、2章GBDT部分のまとめ記事を書こうと思います。*1
この記事を書くにあたって、できるだけ数式の解釈を書くように心がけました。数式の意味をひとつひとつ追っていくことは、実際にXGBoost（またはLightGBMやCatBoostなどのGBDT実装）を使う際にも役立つと考えています。たとえばハイパーパラメータがどこに効いているかを理解することでチューニングを効率化したり、モデルを理解することでよりモデルに合った特徴量のエンジニアリングができるのではないかと思います。

また、この記事に限りませんが、記述に間違いや不十分な点などあればご指摘頂ければ嬉しいです。

XGBoost論文

目的関数の設定

一般的な状況として、サンプルサイズがで特徴量の数がのデータに対する予測モデルを構築することを想定しましょう。*2
今回はツリーをアンサンブルした予測モデルを構築します。
のツリーを加法的に組み合わせた予測モデルは以下のように定式化できます。*3

\begin{align}
\hat{y}_{i} &= \phi\left(\mathbf{x}_{i}\right)=\sum_{k\in\mathcal{K}} f_{k}\left(\mathbf{x}_{i}\right),\\
\text{where}\quad f_{k} \in \mathcal{F} &= \left\{f(\mathbf{x})=w_{q(\mathbf{x})}\right\}\left(q : \mathbb{R}^{m} \rightarrow \mathcal{T},\; \mathcal{T} = \{1, \dots, T\}, \; w \in \mathbb{R}^{T}\right)
\end{align}

ここで、はひとつひとつのツリーを表しています。ツリーは特徴量 が与えられると、それをに従って各ノードに紐づけ、それぞれのノードに対応する予測値を返します。そして、ひとつひとつのツリーの予測値を足し合わせることで、最終的な予測結果とします。

では、具体的にツリーをどうやって作っていくかを決めるために、最小化したい目的関数を設定します。

\begin{align}
\mathcal{L}(\phi) &= \sum_{i\in \mathcal{I}} l\left(y_{i}, \hat{y}_{i}\right)+\sum_{k\in \mathcal{K}} \Omega\left(f_{k}\right), \\
\text{where}\quad\Omega(f) &= \gamma T+\frac{1}{2} \lambda\|w\|^{2} = \gamma T+\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2}
\end{align}

ここで、は損失関数で、たとえば二乗誤差になります。ただし、単に二乗誤差を最小化するのではなく、過適合を回避して汎化性能を上げるために正則化項が追加されています。なお、とはハイパーパラメータであり、交差検証などで最適な値を探索する必要があります。

• の第一項はツリーのノードの数に応じてペナルティが課されるようになっています。ハイパーパラメータを大きくするとよりノード数少ないツリーが好まれるようになります。ツリーの大きさに制限をかけることで過適合を回避することが目的です。
• の第二項は各ノードが返す値の大きさに対してペナルティがかかることを意味しています。ハイパーパラメータを大きくすると、（絶対値で見て）より小さいが好まれるようになります。が小さいということは最終的な出力を決める部分で足し合わされる値が小さくなるので、過適合を避けることに繋がります。

勾配ブースティング

さて、目的関数を最小化するような個のツリー構築したいわけですが、個ツリーを同時に構築して最適化するのではなく、個目のツリーを作る際には、個目までに構築したツリーを所与として、目的関数を最小化するようなツリーを作ることにしましょう。*4

\begin{align}
\min_{f_k}\;\mathcal{L}^{(k)}=\sum_{i\in\mathcal{I}} l\left(y_{i}, \hat{y}_i^{(k - 1)}+f_{k}\left(\mathbf{x}_{i}\right)\right)+\Omega\left(f_{k}\right)
\end{align}

このステップで作成する個目のツリーを合わせた予測値はであり、個目までのツリーではうまく予測できていない部分に対してフィットするように新しいツリーを構築すると解釈できます。このように残差に対してフィットするツリーを逐次的に作成していく手法をブースティングと呼びます。

さて、損失関数を直接最適化するのではなく、その2階近似を最適化することにしましょう。後にわかるように、2次近似によって解析的に解を求めることができます。の周りで2階のテイラー展開を行うと、目的関数は以下で近似できます。

\begin{align}
\mathcal{L}^{(k)} &\approx \sum_{i\in\mathcal{I}}\left[l\left(y_{i}, \hat{y}^{(k - 1)}\right)+g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right] + \Omega\left(f_{k}\right),\\
\text{where} \quad g_i &= \frac{\partial }{\partial \hat{y}^{(k - 1)}}l\left(y_{i}, \hat{y}^{(k - 1)}\right),\\
h_i &= \frac{\partial^2 }{\partial \left(\hat{y}^{(k - 1)}\right)^2}l\left(y_{i}, \hat{y}^{(k - 1)}\right)
\end{align}

ここで、とはそれぞれ損失関数の1階と2階の勾配情報になります。勾配情報を使ったブースティングなので勾配ブースティングと呼ばれています。*5
今回を動かすことで目的関数を最小化するので、と関係ない第一項は無視できます。

\begin{align}
\tilde{\mathcal{L}}^{(k)} &=\sum_{i\in\mathcal{I}}\left[g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right]+\gamma T+\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2} \\
&= \sum_{t \in \mathcal{T}}\left[\sum_{i \in \mathcal{I}_{t}} g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} \sum_{i \in \mathcal{I}_{t}} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right]+\gamma T +\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2}\\
&=\sum_{t \in \mathcal{T}}\left[\left(\sum_{i \in \mathcal{I}_{t}} g_{i}\right) w_{t}+\frac{1}{2}\left(\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}\right) w_{t}^{2}\right]+\gamma T
\end{align}

• 1行目の式では、と関係ない第一項を取り除き、の中身を書き下しました。
• 2行目への変換ですが、全てのについて足し合わせている部分を、まずノードに所属する部分を足し合わせてから、全てのノードについて足し合わせるように分解しています。
• 3行目への変換では、同じノードに所属するは全てを返すというツリーの性質を利用しています。また、の共通部分をくくっています。

さて、はに関しての2次式なので、解析的に解くことができます。

\begin{align}
w_{t}^{*}=-\frac{\sum_{i \in \mathcal{I}_{t}} g_{i}}{\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}}
\end{align}

以上で、個目のツリーに関して、各ノードが返すべき値が解析的に求まりました。*6
この式からもハイパーパラメータを大きくするとが（絶対値で見て）小さくなることが見て取れます。このを元の目的関数に代入してあげることで

\begin{align}
\tilde{\mathcal{L}}^{(k)}(q)=-\frac{1}{2} \sum_{t\in\mathcal{T}} \frac{\left(\sum_{i \in \mathcal{I}_{t}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}}+\gamma T
\end{align}

を得ます。あとはツリーの構造、言い換えれば特徴量の分割ルールを決める必要があります。たとえば、一番シンプルなケースとして、全く分割を行わない場合（）と一度だけ分割を行う場合（に分割）を比較しましょう。分割による目的関数の値の減少分は

\begin{align}
\mathcal{L}_{\text{split}}&= -\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}} h_{i}}+\gamma - \left(-\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}_L} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_L} h_{i}}-\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}_R} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_R} h_{i}}+2\gamma \right)\\
&= \frac{1}{2}\left(\frac{\left(\sum_{i \in \mathcal{I}_L} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_L} h_{i}}+\frac{\left(\sum_{i \in \mathcal{I}_R} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_R} h_{i}} - \frac{\left(\sum_{i \in \mathcal{I}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}} h_{i}}\right) - \gamma
\end{align}

であり、これがプラスなら分割を行い、マイナスなら分割を行わないということになります。上式からも、ハイパーパラメータ を大きくするとより分割が行われなくなることが見て取れます。

ところで、そもそもどの特徴量のどの値で分割するべきかなのでしょうか？一番ナイーブな考え方は、全ての変数に対して全ての分割点を考慮して、一番目的関数の値を減少させるような分割を選ぶというものがあります。ただし、この方法は膨大な計算量が必要になるため、XGBoostでは近似手法が提供されており、3章に記述されています。さらに、4章移行では並列計算や比較実験などが記されています。

まとめ

XGboostの論文を読んだので、自身の理解を深めるために2章GBDT部分のまとめ記事を書きました。
今までなんとなく使っていたハイパーパラメータが具体的にどの部分に効いているのか学ぶことができて、とても有意義だったと感じています。

参考文献

[1603.02754] XGBoost: A Scalable Tree Boosting System
Overview of tree algorithms from decision tree to xgboost
NIPS2017読み会 LightGBM: A Highly Efficient Gradient Boosting Decision T…
XGBoostの概要 - ともにゃん的データ分析ブログ

*1:なお、元論文のノーテーションがおかしい/統一的でないように感じたので、一部表記を変更しています。

*2:元論文ではですが、 とが混じるとややこしい気もするのでにしました。

*3:元論文のとなっているのですが、はインデックスを返す関数なので、タイポかと思われます。

*4:元論文の添字がツリーの数とややこしいのでのまま進めることにしました。

*5:たぶん

*6:論文にはこれがimpuityみたいなものと書かれているのですが、不勉強で理解できませんでした。