FriedmanのGradient Boosting

# はじめに
今ではお馴染みとなったGradient Boostingですが、きちんと説明を読んだことがなかったので調べてみました。

Boostingというと色々ありますが、今回まとめたのはFriedmanのGradient Boostingです。


# 参考
- https://en.wikipedia.org/wiki/Gradient_boosting  
- [Gradientboostingmachines,atutorial(PDF)](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=4&ved=0ahUKEwiAvL2V6avVAhXEgrwKHUKaBcYQFghBMAM&url=http%3A%2F%2Fjournal.frontiersin.org%2Farticle%2F10.3389%2Ffnbot.2013.00021%2Fpdf&usg=AFQjCNGVcZg5TVEnSDMAnoNYcSqb2I5QWQ)

# 手順
## 定義
説明変数を$\mathbf{x}_i(i=1,...,N)$とし、目的変数を$t_i (i=1,...,N)$とします。説明変数の次元は$d$とします。

予測器を$f$と書き、予測値を$y$で表します。我々が求めたいものは、損失関数$\Psi$を最小にする$\hat f$です。

```math
\begin{align}
y_i &= \hat f(\mathbf{x_i}) \\
\hat f &= \arg \min_{f \in \mathcal{F}} \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}))
\end{align}
```

ここで$\mathcal{F}$は関数$f$の空間です。

関数推定の問題を扱いやすくするために、関数$f$を関数空間がパラメータ$\theta$で張られる族$f(\mathbf{x}, \theta)$に制限します。すると上の式は次のように書き直せます。

```math
\begin{align}
y_i &= \hat f(\mathbf{x_i}) \\
\hat f(\mathbf{x}) &= f(\mathbf{x}, \hat\theta) \\
\hat \theta &= \arg \min_{\theta \in D} \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}, \theta))
\end{align}
```

## 勾配降下法
$\theta$を導入したことで全データの損失関数を次のように表せます。

```math
J(\theta) = \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}, \theta))
```

$J$が最小になる$\theta$が$\hat \theta$です。解析的には$J$のgradを計算し、大域的に極小となる点を求めればよいです。しかし、損失関数の形によっては解析的に解を求めるのが難しいこともあります。そこで、$J$の微分を計算し、$J$が小さくなる方向に$\theta$を変化させる、ということを逐次的に行います。

```math
\theta_{\mathrm{new}} =\theta - \eta \nabla_\theta J(\theta)
```

上記の式では更新ごとにすべての訓練用データが含まれていますが、データ一つごとに$\theta$を更新すればオンライン学習に、少量ごとに更新すればミニバッチ法になります。ちなみに、訓練データをすべて使い切るのに必要な回数をイテレーション数と呼びます。

```math
\mathrm{イテレーション数} = \mathrm{データ数} / \mathrm{バッチサイズ}
```

## Gradient Boosting
前章での$\theta$の更新ですが、勾配法との比較のために次のように書き直してみましょう。

イテレーションを$M$回行うとし、各イテレーションを$m=1,...,M$でラベルします。更新の式を次のように書き換えます。

```math
\begin{align}
\theta^{(m)} &= \theta^{(m-1)} + \delta\theta^{(m)}\\
\delta \theta^{(m)} &= - \eta \nabla_\theta J(\theta^{(m-1)})
\end{align}
```

ここで$\theta^{(0)}$は初期値で、適当に与えます。すると$M$回のイテレーションで$\theta$が次のように更新されることが分かります。

```math
\theta^{(M)} = \theta^{(0)} + \sum_{m=1}^{M}\delta\theta^{(m)}
```

式からわかるように、初期パラメータ$\theta^{(0)}$に修正$\delta\theta^{(m)}$を加えていくことで最適な$\hat\theta$を求める仕組みになっています。

そこで、この発想を関数空間に適用することを考えます。すなわちパラメータを修正するのではなく、初期予測器に対して修正予測器を加えていく仕組みを考えます。

```math
f^{(M)} = f^{(0)} + \sum_{m=1}^{M}\delta f^{(m)}
```

勾配法とのアナロジーで更新は次のようにすれば良いとわかります。

```math
\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) + \delta f^{(m)}(\mathbf{x})\\
\delta f^{(m)} &= -\rho \nabla_f J(f^{(m-1)})
\end{align}
```

上の式は素朴なアナロジーでしたが、関数空間の場合はこの更新をブラッシュアップできます。

どういうことかというと、$\theta$の空間の場合、どれくらい$\theta$をシフトさせれば良いかわかりませんでしたので、学習率$\eta$を手動でチューニングする必要がありました。ところが、関数空間の場合は手元に予測器があるので、予測がより良くなるように、すなわち損失関数が小さくなるように関数のシフトを決めることができます。

```math
\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) - \rho^{(m)}\nabla_f J\left(f^{(m-1)}(\mathbf{x})\right) \\
\rho^{(m)} &=  \arg \min_{\rho} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) - \rho \nabla_f J\left(f^{(m-1)}(\mathbf{x_i})\right)\right)
\end{align}
```

この発想で考えれば、もはやシフトに使用する関数も自由に選べます。したがってより一般には次のように書けます。

```math
\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) +g^{(m)}(\mathbf{x}) \\
g^{(m)} &=  \arg \min_{g\in \mathcal{G}} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + g(\mathbf{x_i})\right)
\end{align}
```

このようにしてGradient Boostingでは、ひとつ前のイテレーションの予測器に新たな予測器を加えていくことで、経験損失を小さくしていきます。

## Friedmanの論文の説明について補足
[Friedmanの論文](https://statweb.stanford.edu/~jhf/ftp/stobst.pdf)は実装のことも考えて書かれており、説明が少し異なるので、これについて補足を残します。

前章の最後で、関数のシフトに任意の$g$を使用しました。Friedmanの論文では、この$g$（と$f^{(0)}$）にパラメータ$\theta$を持つ弱学習器$h(\mathbf{x}, \theta)$を考えます。すなわち、前章の最後の式は次のようになります。

```math
\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) + \beta^{(m)} h(\mathbf{x}, \theta^{(m)}) \\
\left(\beta^{(m)}, \theta^{(m)}\right) &=  \arg \min_{\beta, \theta} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \beta h(\mathbf{x}_i, \theta)\right)
\end{align}
```

$g$をこのように置くのは、任意の損失関数に対応するためです。損失関数が複雑な場合、$\beta^{(m)}, \theta^{(m)}$を求めるのは大変です。そこで、gradient boosting近似を行います。

まず、損失関数の勾配に予測器をフィットさせます。

```math
\begin{align}
\theta^{(m)} &= \arg \min_{\rho, \theta}\left[ \sum_{i=1}^N\left(\tilde y_{i}^{(m)} - \rho h(\mathbf{x}_i, \theta) \right)^2\right] \\
\tilde y_{i}^{(m)} &= -\left[\frac{\partial\Psi(t_i, f)}{\partial f}\right]_{f=f^{(m-1)}, \mathbf{x}=\mathbf{x}_i}
\end{align}
```

次に、この弱学習機の係数を最適化します。

```math
\beta^{(m)} =  \arg \min_{\beta} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \beta h(\mathbf{x}_i, \theta^{(m)})\right)
```

この手順によって、最適化問題を最二乗問題と損失関数$\Psi$の1変数最適化問題に分解しました。

予測器の具体例として、弱学習器である決定木を考えます。この決定木は葉を$L$個持つとします。また、この決定木は回帰木です。イテレーション$m$の決定木により、説明変数$\mathbf{x}$の空間が$L$個に分割されます。この分割された領域を$R_m^l (m=1,...,M;l=1,...,L)$でラベルします。gradient boosting近似によって、決定木は次の予測値を出力するように学習します。

```math
\begin{align}
h(\mathbf{x}, \{R_m^l\}_{l=1}^L) &= \sum_{l=1}^L \bar y_l^{(m)} 1(\mathbf{x} \in R_m^l)\\
\bar y_l^{(m)} &= \frac{1}{\#(\mathbf{x}_i \in R_m^l)}\sum_{\mathbf{x}_i \in R_m^l}\tilde y_{i}^{(m)}
\end{align}
```

あとは$\beta$を決めるのですが、領域ごとに$\beta$を最適化できます。予測器が出力する値は$R_m^l$の各領域内で定数なので、$\beta$との積にしてしまって、これを$\gamma$とすると次のように書けます。

```math
\gamma_{l}^{(m)} = \arg \min_{\gamma} \sum_{\mathbf{x}_i\in R_m^l} \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \gamma\right)
```

これによって学習器が更新されます。

```math
f^{(m)}(\mathbf{x}) = f^{(m-1)}(\mathbf{x}) +  \gamma_{l}^{(m)}1(\mathbf{x} \in R_m^l)
```

オーバーフィッティングを避けるために、学習率$\nu(0<\nu<1)$を導入して次のようにする場合もあります。

```math
f^{(m)}(\mathbf{x}) = f^{(m-1)}(\mathbf{x}) +  \nu \gamma_{l}^{(m)}1(\mathbf{x} \in R_m^l)
```

これをShrinkageと呼びます。Friedmanによると$\nu<0.1$とすると汎化性能が上がるそうです。

## おわりに
Boostingアルゴリズムは残差を予測する予測器を加えていくというアルゴリズムですが、他にもいくつか手法があり、それぞれアルゴリズムが結構違うので調べてみると面白いです。

- XGBoost -> [XGBoostの概要](http://kefism.hatenablog.com/entry/2017/06/11/182959)
- AdaBoost -> [AdaBoostについて調べたのでまとめる](https://www.st-hakky-blog.com/entry/2017/08/08/025846)