とてもシンプルなAdamを実装する（スカラー）

ディープラーニングでよく使用されるソルバー、**Adam (Kingma, Ba, 2015)**を理解するためにpythonでシンプルな実装を行った。
最適化の対象は関数$y=x^\theta$の$\theta$とし、$x, y, \theta$いずれもスカラーとした。
まだニューラルネットワークは登場しない。



## アルゴリズム
Adamの理論的な面について筆者は別に記事を書いている([俺はまだ本当のAdamを理解していない]
(https://qiita.com/exp/items/99145796a87cc6cd47e1))。
Adamは生の勾配ではなく、勾配の平均と(中心化前の)分散を使って、パラメータの更新幅が大きくなりすぎないように調節していることが特徴だった。また、この平均と分散を0にバイアスさせない目的で定数$\beta_1,\beta_2$による補正を行っている。

Adamのアルゴリズムを再掲する。

###アルゴリズム1
$g_t^2$は要素ごとの2乗$g_t \odot g_t $を表す。
ベクトル上のすべての操作は要素ごとに行う。
$\beta_1^t, \beta_2^t$は$\beta_1, \beta_2$の$t$乗のこと。

----
>
**Require:**$\quad \alpha \quad$ ステップサイズ 
**Require:**$\quad \beta_1, \beta_2 \in [0,1) \quad$モーメント推定に使う指数減衰率
**Require:**$\quad f(\theta) \quad$ パラメターを$\theta$とする確率的目的関数
**Require:**$\quad \theta_0 \quad$ 初期パラメータベクトル
$\quad m_0 \leftarrow 0 \quad$(最初の一次モーメントベクトル)
$\quad v_0 \leftarrow 0 \quad$(最初の二次モーメントベクトル)
$\quad t \leftarrow 0 \quad$(タイムステップ初期化)
$\quad$ **while** $\;\theta_t$が収束していない **do**
$\qquad t \leftarrow t+1$
$\qquad g_t \leftarrow \nabla_\theta f_t (\theta_{t-1}) \quad$ 目的関数の時刻tでの勾配 
$\qquad m_t \leftarrow \beta_1 \cdot m_{t-1} +(1-\beta_1) \cdot g_t \quad$ 一次のモーメントを更新（バイアスあり）
$\qquad v_t \leftarrow \beta_2 \cdot v_{t-1} + (1-\beta_2) \cdot g_t^2 \quad$　二次の生のモーメントを更新（バイアスあり）
$\qquad \hat{m_t} \leftarrow \frac{m_t}{(1-\beta_1^t)} \quad$ バイアス修正した一次のモーメント
$\qquad \hat{v_t} \leftarrow \frac{v_t}{(1-\beta_2^t)} \quad$　バイアス修正した二次の生のモーメント
$\qquad \theta_t \leftarrow  \theta_{t-1} - \alpha \cdot \frac{\hat{m_t}}{\sqrt{\hat{v_t}+\epsilon}}  \quad$　パラメータ更新
$\quad$**end while**
**return $\theta_t$**

----
アルゴリズムで使用される固定値の推奨設定：

```math
\alpha=0.001 \\
\beta_1=0.9  \\
\beta_2=0.999 \\
\epsilon=10^{-8}
```

##問題設定
冒頭にもある通り、最適化の対象は

```math
y=x^\theta
```
の$\theta$とする。最適な$\theta$は$\theta=2$になるようにする。
つまり、適当に$x$を選んで、$x^2$を正解としてデータセットを作り、ランダムな$\theta$から始めて、Adamで最適化する。

ちなみに、$y=x^2$のグラフはこんな感じ。
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/65afd6ff-55cb-7cff-fc0a-11bd16542cdf.png>

###y=x<sup>θ</sup>の問題点
####勾配
最適化のためには対象となる目的関数の勾配を計算する必要がある。勾配(勾配ベクトル)とは
>偏微分係数を並べたもの

であり、
>勾配ベクトルの向きは，今いる点からちょっと動いたときに関数の値が一番大きくなる向き

である。実用的には`大きくなる`か小さくなるかは符号をつけかえて調節すればよい。 **高校数学の美しい物語**様の[勾配ベクトルの意味と例題](https://mathtrain.jp/koubai)より引用した。

####偏微分
モデルは$y=x^\theta$であるが、説明のために$y$を$y'$と書き換える。つまり

```math
y'=x^\theta \tag{1}
```

今回の目的関数はこの出力$y'$に正解データ$y$との２乗誤差関数を適用した

```math
f=(y' - y)^2 \tag{2}
```

である。これを$\theta$に関して偏微分すると、

```math
\frac {\partial f}{\partial \theta}=2(y' - y)(\frac {\partial y'}{\partial \theta}) \tag{3}
```
となる。右辺の$\frac {\partial y'}{\partial \theta}$は$(1)$式の$\theta$に関する偏微分であり、

```math
\frac {\partial y'}{\partial \theta}=(log x)x^\theta \tag{4}
```
である。すぐわかる通り、$(4)$式が最適化処理内にあるとまずい。なぜなら$x \leq 0$の場合に$log x$が**計算できない。**
これを避けるため、本記事では $ x \geq 0 $ を満たすようにデータセットを作成した。$x \gt 0$でないのは、ごく小さい値($\epsilon=10^{-8}$)を$x$に加えることで$x=0$の場合も許すようにしたからである。

つまりデータセットは
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/14babebe-67d4-f299-7cf7-9391601cb3b3.png>
こんな感じである。
このグラフを見ていると、”これってx=0のときに最小になるんだよね！”と言いたくなるが、**このグラフはデータセットである。目的関数ではない。**

ちなみに、筆者は意識していなかったが、ニューラルネットワークを構成する関数でも、上記$log x$のようなものが含まれた場合は計算できない可能性が高くなるはず。（$x \leq 0$なんて入力はどこにでもある）
こういう関数を使わなくても、**+**と**\***だけで表現力十分ってことだろうか

## 実装

```python:adam.py
# -*- coding: utf-8 -*-
import math
import random

#定義
ALPHA=0.001 #ステップサイズ
BETA_1=0.9 #一次モーメントの指数減衰率
BETA_2=0.999 #二次モーメントの指数減衰率
THETA_0=random.random() #最適化するθの初期値
EPS=1e-08 #イプシロン
EPOCH=1000 #エポック数

#初期化------
m_0=0 #一次モーメント
m_t_old=m_0 #１時刻前の一次モーメント

v_0=0 #二次モーメント
v_t_old=v_0 #１時刻前の二次モーメント

t=0 #タイムステップ

#θ
theta_t=THETA_0
theta_t_old=THETA_0


#関数の定義------
#y=x^θ
def parabola(x, theta):
    result=pow(x,theta)
    return result
  
#上記の微分
def dy_dtheta(x,theta):
    result=math.log(x+EPS)*pow(x,theta)
    return result

#2乗誤差
def error_func(x, y):
    subs=x-y
    result=pow(subs, 2)
    return result

#２乗誤差の微分
def df_dy(f, y):
    result=2*(f-y)
    return result

#可視化用------
gts=[]
mts=[]
vts=[]
mthats=[]
vthats=[]
thetas=[]
losses=[]

#データセット------
DS_X=[]
DS_y=[]
NUM_DATA=1000
# yとして、x^2 を使う。
# 0<= x < 100
for d in range(NUM_DATA):
    x=d*0.1
    DS_X.append(x)
    DS_y.append(parabola(x, 2)) 

# 実行------
# 収束までではなくエポック数だけ回す
while  t<EPOCH :

    #データはランダムに送る
    idx=random.randint(0,NUM_DATA-1)
    x=DS_X[idx]
    y=DS_y[idx]

    #時刻増やす 
    t=t+1
    
    #目的関数の勾配g_t
    y_prime=parabola(x,theta_t_old)
    g_t=df_dy(y_prime,y)*dy_dtheta(x,theta_t_old)
    gts.append(g_t)

    #誤差
    losses.append(error_func(y_prime, y))

    # 一次モーメントの更新
    # beta_1が大きいと過去の値が支配的になる
    m_t=BETA_1*m_t_old + (1- BETA_1)* g_t 
    mts.append(m_t)
    # 二次モーメントの更新
    v_t=BETA_2*v_t_old + (1-BETA_2) * pow(g_t,2)
    vts.append(v_t)
    

    #一次モーメントのバイアス補正
    m_t_hat=m_t / (1-pow(BETA_1, t))
    mthats.append(m_t_hat)
    #二次モーメントのバイアス補正
    v_t_hat=v_t / (1-pow(BETA_2, t))
    vthats.append(v_t_hat)

    #パラメータθを更新
    theta_t=theta_t_old -ALPHA * m_t_hat / (math.sqrt(v_t_hat + EPS)) 
    #古いθとして保持
    theta_t_old=theta_t
    thetas.append(theta_t)


#結果
print("result:",theta_t)
```

##実行
**Google Colab**で実行した。
> result: 2.0003402224491365

こんな感じで、2.0に近い値が出てくれば成功。

##学習の可視化
次のようにグラフをプロットして確認した。横軸は時刻t。

```python
ts=range(NUM_DATA)
plt.plot(ts, losses,  lw=0.4, alpha=0.7) 
```
###誤差
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/3e7c1285-2139-ab24-cf28-57229152b07b.png>
一般的な誤差のグラフと比べると不安だが、ミニバッチ的なものを使っていないのでこのようになるのだと思う。およそ800エポック程度で収束している。

###パラメータθ
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/43a6207d-8d73-b81e-ac2e-860d10104cf5.png>
滑らかに2.0へ向かっている。

###勾配g<sub>t</sub>
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/8c12995d-8275-ea9c-c87d-73dc7261c138.png>
勾配がない状態から700エポック程度で最大になり、また減少する。
なぜこのようになるのだろうか？
論文中で**"真の勾配方向と一致しているか不確かなときはステップサイズが0に近づく"**といった説明があった。
逆に連続して真の勾配方向に近づいていれば勾配は増大する？計算式を見た感じだと特にそんな気はしないが・・・

###勾配の平均m<sub>t</sub>
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/3715fcb7-63eb-f558-7d64-64f084179142.png>
勾配$g_t$とほぼ同じ。$\beta_1=0.9$で、

```python
# beta_1が大きいと過去の値が支配的になる
m_t=BETA_1*m_t_old + (1- BETA_1)* g_t 
```
だったので、1時刻前の$g_t$とほぼ同じなのだろう。

###バイアス補正した勾配の平均m<sub>t</sub>^
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/a5247178-9ac6-b683-5cda-bf1f42f1f810.png>
丸を付けた個所以外はほぼ$m_t$と変わりない。
0で初期化された学習初期だけ補正が強く働いていることがわかる。

###勾配の分散v<sub>t</sub>
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/5a9fd1dd-346a-3476-7f61-e61964522777.png>
ピーク時点は$m_t$同様。$\beta_2=0.99$で、

```python
# 二次モーメントの更新
v_t=BETA_2*v_t_old + (1-BETA_2) * pow(g_t,2)
```
見た目にも$g_t$の2乗であることがわかりやすい($g_t$は下向きのグラフだった)

###バイアス補正した勾配の分散v<sub>t</sub>^
<img src=https://qiita-image-store.s3.amazonaws.com/0/296108/2f0ab9f2-9c87-88c3-5810-bddd8aacb85d.png>
$v_t$と変わりない形だが、丸を付けた個所がわかりやすい違い。
全体的に$v_t$より値が大きい。(分かりにくいが、このグラフの左上に$1e14$とあり、$v_t$は$1e13$。一桁程度の違いはある。)


##まとめ

* スカラーのパラメータ$\theta$に対して最適化を行うAdamをpythonで実装した
* 目的関数は$y=x^\theta$に2乗誤差を適用したもの
* パラメータ$\theta$の収束、Adamのバイアス補正の挙動の一部等を可視化で確認できた

##展望

* パラメータをベクトル化する
* $y=x^\theta$というモデルがあまり良くない(入力に対する制約が強い)
* Adam以外のソルバーとの比較がほしい
* ニューラルネットワークが出てこないぞ