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

Gradient Boosting Machine#

複数の関数に重みをかけて足し合わせた関数

\[ \newcommand{\argmin}{\mathop{\rm arg~min}\limits} f(x) = f_0(x; \theta_0) + \beta_1 f_1(x; \theta_1) + \cdots + \beta_M f_M(x; \theta_M) \]

の形で予測モデルを構築することを考える。ここで\(\theta_0, \dots, \theta_M\)は関数を形づくるパラメータ(例えば線形回帰の重みや決定木の分岐の閾値)である。

このモデルは\(\beta_1, \beta_2, \dots, \beta_M\)\(\theta_0, \theta_1, \dots, \theta_M\)のパラメータを推定する必要がある。 今回はすべてのパラメータを一度に学習するのではなく、\(\beta_m f_m(x; \theta_m)\)を一つずつ学習していく方法を考える。具体的には次のように行う。

前向き段階的加法モデリング(forward stagewise additive modeling)

  1. \(f_0(x) = 0\)で初期化

  2. \(m = 1\)から\(M\)までについて、

    1. パラメータを推定する:\((\beta_m, \theta_m) = \argmin_{\beta, \theta} \sum^N_{i=1} L(y_i, f_{m-1}(x_i) + \beta f(x_i; \theta))\)

    2. 新たなモデルを足す:\(f_m(x) = f_{m-1}(x) + \beta_m f(x; \theta_m)\)

これを前向き段階的加法モデリング(forward stagewise additive modeling)という。ブースティングはこの方法でアンサンブル学習を行う。

誤差関数が二乗誤差\(L(y, f(x)) = (y - f(x))^2\)の場合、

\[\begin{split} \begin{align} L(y_i, f_{m-1}(x_i) + \beta f(x_i; \theta)) &= (y_i - f_{m-1}(x_i) - \beta f(x_i; \theta))^2 \\ &= (\text{residual}_{i,m-1} - \beta f(x_i; \theta))^2 \end{align} \end{split}\]

となり、\(m-1\)回目のモデルの残差\(\text{residual}_{i,m-1} = y_i - f_{m-1}(x_i)\)を近似するように\(m\)回目のモデル\(\beta f(x_i; \theta)\)を学習させていると捉えることができる。残差が大きければそれだけ訓練中に重視されるため「間違えた箇所を重点的に学習する手法」とも捉えることができる。

最適化の観点からの説明#

勾配降下法#

数理最適化において関数の最小化問題

\[ \newcommand{\b}[1]{ \boldsymbol{#1} } \min_\b{x} f(\b{x}) \]

を解く方法のひとつに勾配降下法(gradient descent method)あるいは最急降下法(steepest descent method)と呼ばれるものがある。これは目的関数の微分のベクトルである勾配

\[\begin{split} \nabla f(\b{x}) = \begin{bmatrix} \frac{ \partial f(x_1) }{ \partial x_1 }\\ \vdots \\ \frac{ \partial f(x_m) }{ \partial x_m } \end{bmatrix} \end{split}\]

を用いて

\[ \b{x}_{m} = \b{x}_{m-1} - \alpha_{m-1} \nabla f(\b{x}_{m-1}) \]

という値の更新を何度も繰り返して最適化を行っていく。ここで\(\alpha\)は学習率と呼ばれるパラメータで、値の更新量が多すぎると最適解を通り過ぎてしまうことがあるので小さめの値を乗じて更新幅を抑えるために用いられる。

最終的に\(M\)回反復して得た最適解\(x^*\)

\[ x^* = x_0 - \alpha_1 \frac{\partial f(x_{1})}{\partial x_{1}} - \alpha_2 \frac{\partial f(x_{2})}{\partial x_{2}} - \cdots - \alpha_M \frac{\partial f(x_{M})}{\partial x_{M}} \]

となり、ブースティングにより得られる予測モデル

\[ f(x) = f_0(x; \theta_0) + \beta_1 f_1(x; \theta_1) + \cdots + \beta_M f_M(x; \theta_M) \]

と同様に重み付き和の形になる。

ブースティング#

ブースティングは勾配降下法を機械学習で行っていると捉えることができる。

機械学習においては予測値\(f(x)\)と実測値\(y\)の誤差の最小化問題

\[ \min_{f(x)} L(y, f(x)) \]

を解きたいため、勾配は誤差関数の予測モデルによる微分\(\frac{ \partial L(y, f(x)) }{ \partial f(x) }\)によって得られる。

二乗誤差\(L(y, f(x)) = \frac{1}{2} (y - f(x))^2\)の場合、負の勾配は残差である

\[ - \frac{ \partial L(y, f(x)) }{ \partial f(x) } = y - f(x) = \text{residual} \]

前向き段階的加法モデルの節で「二乗誤差の場合は残差を近似するように学習している」と述べた。

\[ L(y_i, f_{m-1}(x_i) + \beta f(x_i; \theta)) = (\text{residual}_{i,m-1} - \beta f(x_i; \theta))^2 \]

これにより、学習されるモデル\(\beta f(x; \theta)\)は負の勾配を学習するようになり、最終的にそれらの和となるモデルは勾配降下法を解いた状態を近似することになる。

正則化つきGBDT#

\(n\)個の観測データがあり、\(m\)次元の特徴量があるとする。

\[ \mathcal{D} = \{ (x_i, y_i) \}, |D| = n, x_i \in \mathbb{R}^m, y \in \mathbb{R} \]

勾配ブースティング決定木のモデルは次のように表される

\[ \hat{y}_i = \phi(x_i) = \sum^K_{k=1} f_k(x_i) \]

ここで

  • \(f_k \in \mathcal{F}\)は予測器

  • \(\mathcal{F} = \{ f(x) = w_{q(x)} \}\)は回帰木の空間

    • \(q: \mathbb{R}^m \to T\)は入力データ\(x\)を木の各葉のインデックスに割り振る写像。\(T\)は各木の葉の数

    • \(w \in \mathbb{R}^T\)は葉の重み(weights)と呼ばれ、予測に使われた葉の出力値。予測値は\(w_{q(x)}\)の和となるので予測値のベースでもある

学習の際は正則化付き誤差関数

\[ \mathcal{L}(\phi) = \sum^n_{i=1} l(\hat{y}, y_i) + \sum_k \Omega(f_k) \]

を最小化する。ここで\(l\)は微分可能な凸関数である誤差関数で、\(\Omega\)は正則化項

\[ \Omega(f) = \gamma T + \frac{1}{2} \lambda ||w||^2 \]

学習は加法的に行うため\(t\)番目の誤差は次のようになる。

\[ \mathcal{L}^{(t)}(\phi) = \sum^n_{i=1} l(y_i, \hat{y}^{(t-1)} + f_t(x_i)) + \Omega(f_t) \]

テイラー展開による近似#

この誤差を二次近似したものを使うことで計算量を削減することもできることが知られている(Friedman et al., 2000

\[ \mathcal{L}^{(t)} \approx \sum^n_{i=1} [ l(y_i, \hat{y}^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) ] + \Omega(f_t) \]

ここで

\[\begin{split} g_i = \frac{ \partial l(y_i, \hat{y}^{(t-1)}) }{\partial \hat{y}^{(t-1)} }\\ h_i = \frac{ \partial^2 l(y_i, \hat{y}^{(t-1)}) }{\partial (\hat{y}^{(t-1)})^2 } \end{split}\]

定数項を省略すると

\[ \tilde{\mathcal{L}}^{(t)} = \sum^n_{i=1} [ g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) ] + \Omega(f_t) \]

\(j\)におけるインスタンス(サンプル)の集合を\(I_j = \{i|q(x_i) = j\}\)と表記すると、次のように書き換えることができる

\[\begin{split} \begin{align} \tilde{\mathcal{L}}^{(t)} &= \sum^n_{i=1} [ g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) ] + \gamma T + \frac{1}{2} \lambda \sum^T_{j=1} w_j^2 \\ &= \sum^T_{j=1} [ (\sum_{i\in I_j} g_i) w_j + \frac{1}{2} (\sum_{i\in I_j} h_i + \lambda ) w^2_j ] + \gamma T \end{align} \end{split}\]

固定した木の構造\(q(x)\)について、葉\(j\)の最適な重みは

\[ w^*_j = -\frac{ \sum_{i\in I_j} g_i }{ \sum_{i \in I_j} h_i + \lambda } \]

となる。

導出

\(j\)についての部分だけ取り出して導関数を0とおいて整理する

\[\begin{split} \frac{ \partial \tilde{\mathcal{L}}^{(t)}_j }{\partial w_j } = \sum_{i\in I_j} g_i + (\sum_{i\in I_j} h_i + \lambda ) w_j = 0 \\ \implies (\sum_{i\in I_j} h_i + \lambda ) w_j = -\sum_{i\in I_j} g_i \\ \implies w_j = -\frac{ \sum_{i\in I_j} g_i }{ \sum_{i\in I_j} h_i + \lambda } = w_j^* \end{split}\]

L1の場合

\[\begin{split} \begin{align} \tilde{\mathcal{L}}^{(t)} &= \sum^n_{i=1} [ g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) ] + \gamma T + \frac{1}{2} \lambda \sum^T_{j=1} w_j^2 + \alpha \sum^T_{j=1} |w_j| \\ &= \sum^T_{j=1} [ (\sum_{i\in I_j} g_i) w_j + \frac{1}{2} (\sum_{i\in I_j} h_i + \lambda ) w^2_j ] + \gamma T \end{align} \end{split}\]

L1の場合

\[\begin{split} \begin{align} \tilde{\mathcal{L}}^{(t)} &= \sum^n_{i=1} [ g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) ] + \gamma T + \frac{1}{2} \lambda \sum^T_{j=1} w_j^2 + \alpha \sum^T_{j=1} |w_j| \\ &= \sum^T_{j=1} [ (\sum_{i\in I_j} g_i) w_j + \alpha |w_j| + \frac{1}{2} (\sum_{i\in I_j} h_i + \lambda ) w^2_j ] + \gamma T \end{align} \end{split}\]

\(\sum^T_{j=1} [\cdot]\)の内側を整理すると

\[ w^*_j = -\frac { \sum_{i\in I_j} g_i \pm \alpha } { \sum_{i \in I_j} h_i + \lambda } \]

最適な重み\(w_j^*\)を誤差関数に戻すと

\[ \tilde{\mathcal{L}}^{(t)}(q) = -\frac{1}{2} \sum^T_{j=1} \frac{ (\sum_{i\in I_j} g_i)^2 } { \sum_{i\in I_j} h_i + \lambda } + \gamma T \]

となり、これ木の構造\(q\)の品質をスコアリングする関数として使うことができる。

導出
\[\begin{split} \begin{align} \tilde{\mathcal{L}}^{(t)} &= \sum^T_{j=1} [ (\sum_{i\in I_j} g_i) (-\frac{ \sum_{i\in I_j} g_i }{ \sum_{i \in I_j} h_i + \lambda }) + \frac{1}{2} (\sum_{i\in I_j} h_i + \lambda ) (-\frac{ \sum_{i\in I_j} g_i }{ \sum_{i \in I_j} h_i + \lambda })^2 ] + \gamma T \\ &= \sum^T_{j=1} [ -\frac{ (\sum_{i\in I_j} g_i)^2 }{ \sum_{i \in I_j} h_i + \lambda } + \frac{1}{2} \frac{ (\sum_{i\in I_j} g_i)^2 }{ \sum_{i \in I_j} h_i + \lambda } ] + \gamma T \\ &= \sum^T_{j=1} [ - \frac{1}{2} \frac{ (\sum_{i\in I_j} g_i)^2 }{ \sum_{i \in I_j} h_i + \lambda } ] + \gamma T \end{align} \end{split}\]

参考#