確率的勾配降下法とは何か、をPythonで動かして解説する

#勾配降下法は何に使う？#
勾配降下法は統計学や機械学習で多く使われています。特に機械学習というのは基本的に何かしらの関数を最小化（最大化）する問題を数値解析的に解くことに帰結する場合が多いです。（e.g. 最小二乗法 → 誤差の二乗和を最小化（[参考](http://qiita.com/kenmatsu4/items/8b4e908d7c93d046110d)）、ニューラルネットワークのパラメータ決定 etc...）

なので、基本的にはひたすら微分して0となるところを探す問題ですね、微分して0。で、その微分して0となる値は何か、をプログラムで解く場合に重要になるのがこの勾配降下法です。幾つか勾配法にも種類がありますがここでは最急降下法、確率的勾配降下法の２つを扱います。まずはイメージをつかむために１次元のグラフで確認していきたいと思います。

##1次元の場合##
1次元の場合は、確率的という概念はなく、ただの勾配降下法になります。
（どういうことか、はのちほど）
1次元の例は、正規分布をマイナスにしたものを使ってみます。$x_{init} = 0$と$x_{init} = 6$のグラフを例として作っています。対象となるグラフは正規分布の密度関数をマイナスして下向きにしたものを使いました。


**初期値 = 0 のグラフ**
グラフの端から始まって、ちゃんと一番小さいところ、最小値で収束していることがわかりますね。青い線のグラフが正規分布のマイナス版で、緑の線のグラフがその微分した値のグラフです。
![normdist_decent_0.0_anim-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/1c51dc9e-c469-e97e-851e-a9a3ddf451fc.gif)

**初期値 = 6 のグラフ**
![normdist_decent_6.0_anim-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/0c80d559-a353-012c-4dd0-27ed5b763775.gif)

各点での傾きをベースに次の$x$の値を決め、またその位置での傾きを調べ次に降下して行く、という事を繰り返して最小値（本当は極小値）にたどり着きます。上の例だと、最小値に近ずくにつれ、傾きもなだらかになっていくので、1回のステップで進む量がだんだん小さくなり、それがある基準より小さくなるとストップする、というロジックを組んでいます。

これら２つのグラフを書いたコードは[こちら](https://gist.github.com/matsuken92/f92f8e8e4a1f924e6d99#file-1-dim-steepest-descent-py)です。

###2次元の場合(2次関数)###
多次元の場合は次のステップに進む向きも考慮して降下していきます。ターゲットとなる関数により動きの特徴が変わるのですが、最初に2次関数 $Z=-3x^2-5y^2+6xy$ という関数を例にグラフを書いてみました。なんだかジグザグ動きながら頂点に向かって行く姿をみることができます。
![quadratic_decent_anim-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/dd0f96bf-a434-a495-3fe7-59c900a1e90f.gif)


コードは[こちら](https://gist.github.com/matsuken92/f92f8e8e4a1f924e6d99#file-2-dim-steepest-descent_quadrant-py)

##2次元の場合##

###2次元の場合(2次元正規分布)###
2次元正規分布（にマイナスをかけたもの）に対して最急降下法を適用するとジグザグ動かず、綺麗に最小値に降下していくことがわかりました。このあたりターゲット関数による収束の様子の違いについては今度調べてみたいと思います。
![bi-normdist_decent_anim-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/b9b22d2b-edca-bbd3-6659-091e559ce2db.gif)


コードは[こちら](https://gist.github.com/matsuken92/f92f8e8e4a1f924e6d99#file-2-dim-steepest-descent_normdist-py)

###最急降下法の仕組み###
最急降下法の仕組みは、まず最初にスタート地点を選び、その点での一番降下する傾きが大きいベクトル、つまり

```math
\nabla f = \frac{d f({\bf x})}{d {\bf x}} = \left[ \begin{array}{r} \frac{\partial f}{\partial x_1} \\ ... \\ \frac{\partial f}{\partial x_2} \end{array} \right]
```
を利用して次のステップに行く方法です。上記のグラフでは、この勾配ベクトルに学習率$\eta$をかけたものを直線で表しています。
なので、

```math
x_{i+1} = x_i - \eta \nabla f
```
というステップを収束するまで繰り返す、という処理を行っています。

##確率的勾配降下法とは##
さて、次に確率的勾配降下法に移りたいと思います。
まず、対象とする関数ですが、以前[回帰分析の説明](http://qiita.com/kenmatsu4/items/1e37da1d55292035d985)で使った回帰直線の誤差関数を使って説明したいと思います。ちょっとその解釈を変えている部分(最尤法を使用)もあるので[文末](http://qiita.com/kenmatsu4/items/d282054ddedbd68fecb0#http://qiita.com/kenmatsu4/items/d282054ddedbd68fecb0#補足回帰直線のパラメータを最尤法で求める)に補足として記載しておきます。

データは、

```math
y_n = \alpha x_n + \beta + \epsilon_n　　　(\alpha=1, \beta=0)\\
\epsilon_n ∼ N(0, 2)　　　　　　　　　　　　　　　　
```

というN=1000個のデータセットを準備して、その回帰直線を求める例を使用します。散布図は以下の通りです。

![scatter-compressor.png](https://qiita-image-store.s3.amazonaws.com/0/50670/adf9be8d-1732-91fb-e1df-76c400bfef61.png)

誤差関数、

```math

E({\bf w})=\sum_{n=1}^{N} E_n({\bf w}) = \sum_{n=1}^{N} (y_n -\alpha x_n - \beta)^2 \\
{\bf w} = (\alpha, \beta)^{\rm T}　　　　　　　　　　　　　　　　　　　　　　　
```
の最小値を求めていきますが、まず勾配を計算するため、それぞれ$\alpha, \beta$で偏微分した勾配ベクトルを導出します。

```math

\frac{\partial E({\bf w})}{\partial \alpha} = \sum_{n=1}^{N} \frac{\partial E_n({\bf w})}{\partial \alpha} = 
\sum_{n=1}^N (2 x_n^2  \alpha +  2 x_n \beta - 2 x_n y_n )\\

\frac{\partial E({\bf w})}{\partial \beta} = \sum_{n=1}^{N}\frac{\partial E_n({\bf w})}{\partial \alpha} =
\sum_{n=1}^N (2\beta + 2  x_n \alpha - 2y_n)　　　
```
なので、勾配ベクトルは

```math
\nabla E({\bf w}) = \frac{d E({\bf w}) }{d {\bf w}} = \left[ \begin{array}{r} \frac{\partial E({\bf w})}{\partial \alpha} \\ \frac{\partial E({\bf w})}{\partial \beta} \end{array} \right] = \sum_{n=1}^{N}
\nabla E_n({\bf w}) = \sum_{n=1}^{N} \left[ \begin{array}{r} \frac{\partial E_n({\bf w})}{\partial \alpha} \\ \frac{\partial E_n({\bf w})}{\partial \beta} \end{array} \right]=
\left[ \begin{array}{r} 
2 x_n^2  \alpha +  2 x_n \beta - 2 x_n y_n \\ 
2\beta + 2  x_n \alpha - 2y_n 　
\end{array} \right]
```
となります。

単回帰直線を求める問題なので、方程式として解いても良いのですが、これを確率的勾配降下法を用いて解くことを行っていきます。ここでN=1000個の全てのデータをすべて使うのではなく、ここからランダムサンプリングして取り出したデータに対して勾配を計算し、次のステップを決めるのがこの手法のキモです。


まずこのデータセットから20個ランダムに取り出してグラフにしてみます。横軸が$\alpha$、縦軸が$\beta$、Z軸が誤差$E(\alpha,\beta)$です。

![SGD_4-compressor.png](https://qiita-image-store.s3.amazonaws.com/0/50670/50b5884b-23a0-8d1e-9597-3b3f3b83fac4.png)

なんとなくグラフの最小値が$\alpha=1, \beta=0$を中心に散らばっているかと思います。今回は1000個のデータからランダムに3つのデータをピックアップしてそれを元にこのグラフを書いています。
初期値を$\alpha=-1.5, \beta=-1.5$、学習率をイテレーション回数$t$の逆数に比例させ、さらにその地点での$E(\alpha, \beta)$の逆数もかけたものを使ってみました。この学習率と初期値の決め方について試行錯誤するしかないようなのですが、何か良い探し方をご存知の方がいれば教えてもらえると嬉しいです。ちょっと間違えるとあっという間に点が枠外に飛んで行って戻ってこなくなります（笑）

勾配を決める誤差関数が乱数に依存しているので毎回変化していることが見て取れます。回帰直線も最初は相当暴れていますが、だんだん大人しくなって収束していく様がわかると思います。

<img src="https://qiita-image-store.s3.amazonaws.com/0/50670/3080df2b-872d-90cf-ae12-61fda5feaaf1.gif" width="512" height="512" ></img>
コードは[こちら](https://github.com/matsuken92/Qiita_Contents/blob/master/General/SGD-lm.ipynb)。

正直、上記のアニメーションの例は収束が良い方のものでして、下記に10000回繰り返した際の$\alpha$と$\beta$の収束具合をグラフにしたものを載せていますが、$\alpha$は真の値1に近づいているのですが、$\beta$は0.5のあたりで収束してしまい、真の値の0にこれ以上近づきそうにありません。

![result_10000r-compressor.png](https://qiita-image-store.s3.amazonaws.com/0/50670/8bbcccd6-21d5-90e1-de4f-941fe68495cb.png)

これはそもそもの今回の回帰直線の誤差関数の特性にあるようですが、下記のグラフを見てもらうとわかるように、$\alpha$については$\alpha=1$のあたりで大体曲線の最小値あたりに落ちています。しかし、縦のラインについてはとても緩やかな変化しかないため、なかなか縦の動きが出る勾配が出てこないのです。これが通常の最急降下法等であればこの曲面自体は変化しないので、谷に沿って緩やかに最小値を目指すのですが、そもそも毎回この曲面が変わってしまうことが原因で$\beta$方向の動きが生まれないようです。

![graph2-compressor.png](https://qiita-image-store.s3.amazonaws.com/0/50670/649f7efe-b7f2-90ef-d86e-e43c1bc02e2d.png)

さて、確率的勾配降下法についてイメージが湧きましたでしょうか？私はこのアニメーショングラフを描くことでイメージがつかめました。
この確率的勾配降下法のアニメーションを描くpythonコードは[こちら](https://gist.github.com/matsuken92/f92f8e8e4a1f924e6d99#file-2-dim-stocastic-descent_quadrant-py)にありますので、よければ試してみてください。

###(補足)回帰直線のパラメータを最尤法で求める###

```math
y_n = f(x_n, \alpha, \beta) + \epsilon_n　　　　(n =1,2,...,N)
```
つまり、誤差は

```math
\epsilon_n = y_n - f(x_n, \alpha, \beta) 　　　　(n =1,2,...,N)  \\
\epsilon_n = y_n -\alpha x_n - \beta　　　　　　　　　　　　　　　　　　　
```
で表すことができます。この誤差$\epsilon_i$を平均=0, 分散=$\sigma^2$の正規分布に従うとすると、誤差関数の尤度は

```math

L = \prod_{n=1}^{N} \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{\epsilon_n^2}{2\sigma^2} \right)
````

であり、対数をとると

```math

\log L = -\frac{N}{2} \log (2\pi\sigma^2) -\frac{1}{2\sigma^2} \sum_{n=1}^{N} \epsilon_n^2
```

となり、パラメーター$\alpha, beta$に依存する項のみを取り出すと

```math
l(\alpha, \beta) = \sum_{n=1}^{N} \epsilon_n^2 = \sum_{n=1}^{N} (y_n -\alpha x_n - \beta)^2
```
を最小化すれば良いことがわかります。
これを誤差関数として

```math

E({\bf w})=\sum_{n=1}^{N} E_n({\bf w}) = \sum_{n=1}^{N} (y_n -\alpha x_n - \beta)^2 \\
{\bf w} = (\alpha, \beta)^{\rm T}　　　　　　　　　　　　　　　　　　　　　　　
```
と表現します。

<!--
**以下、当初確率的勾配降下法の私の解釈に誤りがあることがいただいたコメントで分かりましたので訂正をしております。同じ勘違

て、書き直しをしています。本来ランダムピックアップするのはデータ点なのですが、下記の説明ではm次元の勾配から１つ選ぶ、という誤解がありました。すみません・・・**


<S>
最急降下法では各点で一番降下する傾きが大きいベクトルを偏微分することで求めて次のステップに行く方法ですが、ここで扱っていた2次元程度の関数であればさほど大きな計算量にならないのですが、数百、数千、もっと多い次元のインプットがあった場合、相当の計算量になってしまい収束しなくなってしまうケースが出てきてしまいます。（私がよく扱っている手書き数字データの場合でも、784次元のデータがインプットになるので、1ステップの計算だけで2次元の392倍！）</S>

<S>なので、</S>

```math
\nabla f = \frac{d f({\bf x})}{d {\bf x}} = \left[ \begin{array}{r} \frac{\partial f}{\partial x_1} \\ ... \\ \frac{\partial f}{\partial x_2} \end{array} \right]
```

<S>のように全方向の勾配ベクトルを使うのではなく、この中からランダムに１つ

```math
\nabla f = \frac{\partial f}{\partial x_i}
```

<S>を選んで（複数選ぶケースもある）その方向のみに降下する事を行います。こうすると、全方向の偏微分を求めなくて済むので計算量がかなり減ります。ランダムに選ぶ、というのが”確率的”と名が付いている由来ですね。</S>

###グラフに書いてみる###
<S>"確率的"というだけあって、毎回、最小値に至るパスが違います。下にグラフを２つ書いてみたので見てみましょう。</S>

![quadratic_decent_stchastic_anim_1-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/0e5ccd94-930c-c1b1-34dc-871518c86f72.gif)
![quadratic_decent_stchastic_anim_2-compressor-2.gif](https://qiita-image-store.s3.amazonaws.com/0/50670/f3f16c2c-fc0e-4ec2-37d7-738425747349.gif)


<S>コードは[こちら](https://gist.github.com/matsuken92/f92f8e8e4a1f924e6d99#file-2-dim-stocastic-descent_quadrant-py)</S>

<S>どちらも正直精度が良い感じがしないですね（笑）ただ、今回は図示するため2次元関数を扱いましたので、この計算量だったら素直に最急降下法を使ったほうが良いということですね。上にある2次関数の最急降下法の場合、n=35回n=で収束基準を満たし終わっていますが、確率的勾配降下法をの場合はn=40回で0.3以上の誤差を残し打ち止めとなっています。（Qiitaのgifアニメの容量制限1MBにより、これ以上の回数をアニメーションで表示できませんでした・・・）</S>
<S>ただし、計算回数で考えると、</S>

<S>・最急降下法</S>
<S>　⇒ 2(次元) × 35(回で収束)  = 70回の計算</S>

<S>・確率的勾配降下法（ランダムに1次元選択する）</S>
<S>　⇒ 1(次元) × 40(回で収束)  = 40回の計算</S>

<S>なので、計算回数が少ないのです。特に次元数が多い場合は劇的に回数が変わってきます。次元数が多い場合は確率的勾配降下法を使ったほうが良いケースがありそうですね。</S>
-->
<!--
```math

f({\bf x}; {\bf \mu}, \Sigma) = \frac{1}{(2\pi)^{D/2}} \frac{1}{|\Sigma|^{1/2} }  \exp{ \left( ({\bf x - \mu})^{{\rm T}}\Sigma^{-1}(x)({\bf x - \mu}) \right)}

```

```math

\frac{d f({\bf x}; {\bf \mu}, \Sigma)}{d{\bf x}} = 2\Sigma^{-1}{\bf x}
\frac{1}{(2\pi)^{D/2}} \frac{1}{|\Sigma|^{1/2} }  \exp{ \left( ({\bf x - \mu})^{{\rm T}}\Sigma^{-1}(x)({\bf x - \mu}) \right)}

```
-->