今回は、optimizerと言えば、ほぼadam一択だったけど、実は性質はいろいろというのを実感したので、まとめておく。
これをやろうと思ったのは、参考の記事を読んでなるほどと思うと同時にいくつかの疑問がおきたのでやってみました。
【参考】
①【前編】Pytorchの様々な最適化手法(torch.optim.Optimizer)の更新過程や性能を比較検証してみた！
###やったこと
・基本的な動作
・疑問
・詳細に見る
・自前optimizer
・optimizerの処理を比較し性質を考える
###・基本的な動作
以下の表は、上記の参考①と以下の参考②を参考に作成しています。
・【参考】
②確率的勾配降下法

名称	数式	備考	コード（リンク）
SGD (Stochastic Gradient Descent)	$W\leftarrow$ $W-\eta\frac{\partial L}{\partial W}$	関数$Ｌ$を最小化（極小化）するために、パラメータを勾配（坂道）を下る方向に動かす。物体運動なら$W$が位置、$L$がポテンシャルに相当し、ポテンシャルの等高線に垂直に運動する。（$W$：パラメータ、$η$：学習率、$L$：$W$の関数、$\frac{\partial L}{\partial W}$：勾配）Stochastic は学習データを毎回ランダムに変更する	pythorch解説 SOURCE CODE FOR TORCH.OPTIM.SGD
Momentum	$W\leftarrow W + v$ $v←αv−η\frac{∂L}{∂W}$	$v$は前回のパラメータの変化量；速度の$\alpha$倍をパラメータ変位に加味する	pythorch解説 SOURCE CODE FOR TORCH.OPTIM.Momentum
Adagrad	$W \leftarrow W - \frac{\eta}{\sqrt{h}} \frac{\partial L}{\partial W}$$h \leftarrow h + \left(\frac{\partial L}{\partial W}\right)^2$	学習率ηを$h$で規格化する	pythorch解説torch.optim.Adagrad(params, lr=0.01, lr_decay=0, weight_decay=0, initial_accumulator_value=0, eps=1e-10)SOURCE CODE FOR TORCH.OPTIM.ADAGRAD
RMSprop (Root Mean Square propagation)	$W \leftarrow W - \frac{\eta}{\sqrt{h + \epsilon}} \frac{\partial L}{\partial W}$ $h \leftarrow \alpha h + (1 - \alpha) \left( \frac{\partial L}{\partial W} \right)^2$	Adagradに$\alpha$で学習率の分母をhと勾配で寄与調整、$\epsilon$は0除算のため	pythorch解説torch.optim.RMSprop(params, lr=0.01, alpha=0.99, eps=1e-08, weight_decay=0, momentum=0, centered=False)SOURCE CODE FOR TORCH.OPTIM.RMSPROP
Adadelta	$W\leftarrow W -\frac{\sqrt{\Delta W + \epsilon}}{\sqrt{h + \epsilon}}\frac{\partial L}{\partial W}$ $h \leftarrow \alpha h + (1 - \alpha) \left( \frac{\partial L}{\partial W} \right)^2$ $\Delta W \leftarrow \beta \Delta W + (1 - \beta) (\Delta W)^2$	RMSpropの更新量$\Delta W = - \frac{\eta}{\sqrt{h + \epsilon}}\frac{\partial L}{\partial W}$、移動平均量$\Delta W$の一次二次の調整された$\Delta W$で次元調整;学習率の設定が不要	pythorch解説torch.optim.Adadelta(params, lr=1.0, rho=0.9, eps=1e-06, weight_decay=0)SOURCE CODE FOR TORCH.OPTIM.Adadelta
Adam	$W \leftarrow W - \frac{\eta \hat{m}}{\sqrt{\hat{v} + \epsilon}}$ $\hat{m} = \frac{m}{1 - \beta_1}$ $m \leftarrow \beta_1 m + (1 - \beta_1)\frac{\partial L}{\partial W}$ $\hat{v} = \frac{v}{1 - \beta_2}$ $v \leftarrow \beta_2 v + (1 - \beta_2)\left(\frac{\partial L}{\partial W}\right)^2$	$\hat{m};\textrm {momentum,}$ $\hat{v};\textrm{ RMSprop の採用}$ $初期値　m_0=0, v_0=0$ $\textrm{推奨params }\eta = 0.001 \beta_1 =0.9, \beta_2 = 0.999 \epsilon = 1e^{-8}$	pythorch解説;torch.optim.Adam(params, lr=0.001, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False)SOURCE CODE FOR TORCH.OPTIM.Adam
AdamW	$W \leftarrow W - \frac{\eta \hat{m}}{\sqrt{\hat{v} + \epsilon}}$ $\hat{m} = \frac{m}{1 - \beta_1}$ $m \leftarrow \beta_1 m + (1 - \beta_1)\frac{\partial L}{\partial W}$ $\hat{v} = \frac{v}{1 - \beta_2}$ $v \leftarrow \beta_2 v + (1 - \beta_2)\left(\frac{\partial L}{\partial W}\right)^2$	$\hat{m};\textrm {momentum,}$ $\hat{v};\textrm{ RMSprop の採用}$ $初期値　m_0=0, v_0=0$ $\textrm{推奨params }\eta = 0.001 \beta_1 =0.9, \beta_2 = 0.999 \epsilon = 1e^{-8}$ 上記にWeight decayを追加　weight_decay=0.01 AdamとAdamWのweight-decayは異なる	pythorch解説;torch.optim.AdamW(params, lr=0.001, betas=(0.9, 0.999), eps=1e-08, weight_decay=0.01, amsgrad=False) SOURCE CODE FOR TORCH.OPTIM.AdamW 論文Decoupled Weight Decay RegularizationリンクのkerasコードAdamW_Keras/AdamW.py

###・疑問
上記参考は分かり易いです。
その結果、以下の疑問がわきました。
①parameterの妥当性；adamやadamWって、蛇行しているけど
②ポテンシャルが$f(x,y)=x^2 + y ^2$なのに、SGD以外はまっすぐ進んでいないけど、ある意味勾配小さいy軸方向に動いているように見える
③adamとadamWの原点近傍の振動？はなんだ
####①parameterの妥当性；adamやadamWって、蛇行しているけど
まず、同じような図を出力してみました。
確かにこの条件だと、同じ絵が出力されました。

しかし、普通に考えれば、AdamやAdamWの振動は学習率が大きすぎるために起こっている振動だと思われます。そこで、学習率Irを小さくして振動が止まる条件で出力させてみると、以下の通りとなりました。

####②SGD以外はまっすぐ進んでいない
ポテンシャルが$f(x,y)=x^2 + y ^2$なのに、SGD以外はまっすぐ進んでいないけど、ある意味勾配小さいy軸方向に動いているように見える。この現象はどこから来るのでしょう。上の表からSGD以外は、$\frac{\partial L}{\partial W}$の係数は、傾きの和で除しているので、等方的な係数にはなっていないことが分かります。つまり、傾きに素直に追従しない。むしろ、除しているので傾きの履歴が大きいものほど、大きな数字で除すことになるので、その方向のベクトルは小さくなりそうです。
ということで、ここではポテンシャルが異方性を持つ場合についてやってみました。
例えば、ポテンシャルの形状が$x^2 + 0.01y^2$の場合、以下の図が得られました。ここでもSGDは素直にまずx=0へ進んでから傾きが小さなy方向へ移動しています。一方、他の手法ではこのような場合でもy方向へ大きく移動してy=0を実現してからx=0へ進んでいきます。

###③adamとadamWの原点近傍の振動？はなんだ
これは、学習率が大きすぎて、振動して収束できないということです。ということで、小さな学習率でやってみた、上記のIr = 0.25の結果を見ると原点に収束しているのが分かります。
####変な挙動
実は、上記のポテンシャル形状を色々変化させてみてみると以下のような結果を得ました。以下は、ポテンシャル形状が$x^2+1e^{-8}y^2$の場合の各optimでの収束の様子です。SGDとAdadeltaは素直にx=0方向に動いており、収束していませんが、その他のoptimでは収束というかy=0を探り当てています。

ちなみに、$x^2+0y^2$の場合は、以下の通りAdamW以外は収束しませんでした。ある意味ほっとしましたが、AdamWは一体どうしたのでしょう。

さらに、$x^2-1e^{-8}y^2$の場合は、収束しませんでした。

###詳細に見る
$x^2 + 0.01y^2$の場合について、収束の様子を3dプロットしてみてみます。
なお、この部分のコードは以下のとおりです。
このコードから、10まで'blue'1000まで'red'そして、10000までは、'black'と'white'に色分けしています。これを見ると収束性のスピードが見えると思います。

収束をプロットするコード

n=10000
for kin in range(6): #6
    fig = plt.figure(figsize=(12, 6))
    ax1 = fig.add_subplot(121 , projection='3d')
    ax2 = fig.add_subplot(1, 2, 2)
    hv = ax2.pcolormesh(X,Y,Z, cmap='jet')
    # 等高線図の生成
    cont=ax2.contour(X,Y,Z,  10, colors=['black']) #Vmax=1
    cont.clabel(fmt='%1.1f', fontsize=14)
    sc0_ = ax1.plot_wireframe(X, Y, Z, rstride=2, cstride=2)
    #sc1_ = ax1.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='jet') 
     
    key = classes[kin]
    j =0
    while j <n:
        if j<10:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c ='blue')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='blue')
            j +=1

        elif j < 1000:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c='red')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='red')
            j += 5
        elif j < 10000:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c='black')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='white')
            j += 200

    if key =='SGD':
        Ir = 0.1*Ir0
    else:
        Ir = Ir0            
    ax1.set_title('{}_Ir={}_{}x^2+{}y^2'.format(key,Ir,r0,r1), fontsize=16)
    ax1.set_xlabel('x', fontsize=16)
    ax1.set_ylabel('y', fontsize=16)    
    ax1.set_zlabel('z', fontsize=16)    
    fig.colorbar(hv, ax=ax2)

    plt.savefig('fig3d_optim_{}_Ir{}_r0{}_r1{}.png'.format(key,Ir, r0,r1))
    plt.pause(3)

###・自前optimizer
上記は、一応の傾向は見えるので、それなりに面白いですが、ある意味らちがあきません。
ということで、以下の参考と同じように自分で上記の表の関数を使って収束の様子を見たいと思います。
【参考】
③Gradient Descent Optimizations
####学習率による収束の様子；VGD； Vanila Gradient Descent

####長細いポテンシャル形状の場合のmomentumの効果
####ポテンシャル形状　
$L=x^2+100y^2$

#####Ir = 0.009 simple VGD

#####Ir = 0.009 with momentum

####長細いポテンシャル形状の場合の学習率による違い
$L=0.05x^2+y^2$
#####Ir = 0.9

#####Ir = 0.75

#####Ir = 0.3

###RMSprop, Adam, Adagrad, そしてAdadeltaの違い
これらを見ると、それぞれの収束性能が垣間見えます。
####L=x^2+100y^2_Ir0.1 with RMSprop

RMSpropのコード

def f2(x):
    return x[0]**2 + 100*x[1]**2

def grad2(x):
    return np.array([2*x[0], 200*x[1]])

def gd2_rmsprop(x, grad, alpha, beta=0.9, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    for i in range(max_iter):
        v = beta*v + (1-beta)*grad(x)**2
        x = x - alpha * grad(x) / (eps + np.sqrt(v))
        xs[i+1, :] = x
    return xs

alpha = 0.1
x0 = np.array([-1,-1])
xs = gd2_rmsprop(x0, grad2, alpha, beta=0.9, max_iter=10)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradient descent with RMSprop')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withRMSprop.png'.format(alpha))
plt.show()

![fig_y=x^2+100y^2_Ir0.1withRMSprop.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/233744/77088280-1688-43ee-8862-59a4275c3ea6.png) ####L=x^2+100y^2_Ir0.1 with Adam

Adamのコード

def gd2_adam(x, grad, alpha, beta1=0.9, beta2=0.999, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    m = 0
    v = 0
    for i in range(max_iter):
        m = beta1*m + (1-beta1)*grad(x)
        v = beta2*v + (1-beta2)*grad(x)**2
        mc = m/(1+beta1**(i+1))
        vc = v/(1+beta2**(i+1))
        x = x - alpha * m / (eps + np.sqrt(vc))
        xs[i+1, :] = x
    return xs

####L=x^2+100y^2_Ir0.1 with Adagrad

Adagradのコード

def gd2_adagrad(x, grad, alpha, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    for i in range(max_iter):
        v = v + grad(x)**2
        x = x - alpha * grad(x) / (eps + np.sqrt(v))
        xs[i+1, :] = x
    return xs

####L=x^2+100y^2_Ir0.01 with Adadelta
このコードはうまく収束していますが、今一つ理論が不明です。

Adadeltaのコード

def gd2_adadelta(x, grad, alpha, beta1=0.9, beta2 = 0.999, beta3 = 0.999, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    h = 0
    for i in range(max_iter):
        h = beta1*h + (1-beta1)*grad(x)*grad(x)
        dw = -alpha*grad(x) / (eps + np.sqrt(h))
        dw = beta2*dw + (1-beta2)*dw*dw
        v = beta3*v + (1-beta3)*grad(x)**2
        x = x - np.sqrt(eps + dw) * grad(x) / np.sqrt(eps + v)
        xs[i+1, :] = x
    return xs

###・optimizerの処理を比較し性質を考える
最後にoptimizerの間の関係を見るために、処理部分を比較したいと思います。
まず、基本のVGD

VGD.py

x = x - alpha * grad(x)

VGDにｖを導入して、$\sqrt{v+grad(x)^2}$で除して規格化して等方性を導入したのがAdagradです。すなわち、もしgrad(x)**2が大きなところに来ると、ほぼx = x -alphaで更新することになり、傾きが急激なのに反映されません。また、vが大きくなると学習率が実質的に小さくなることとなります。

Adagrad.py

v = v + grad(x)**2
x = x - alpha * grad(x) / (eps + np.sqrt(v))

RPSpropは、Adagradのv履歴の寄与と現在の傾きの寄与の割合を調整できるようにしたものである。betaにより、学習率の調整が可能となっている。

RMSprop.py

v = beta*v + (1-beta)*grad(x)**2
x = x - alpha * grad(x) / (eps + np.sqrt(v))

そして、Adamは、Adaptive Moment Estimationの略であり、moment項mと上記のRMSpropのいいとこ取りをしているモデルです。moment項はそれまでの傾きの履歴を保持するような項になっています。

Adam.py

m = beta1*m + (1-beta1)*grad(x)
v = beta2*v + (1-beta2)*grad(x)**2
mc = m/(1+beta1**(i+1))
vc = v/(1+beta2**(i+1))
x = x - alpha * m / (eps + np.sqrt(vc))

###まとめ
・optimizerを変更して遊んでみた
・図に表すことにより収束性能の一端が垣間見えた
・学習率の影響は大きいので、今後は慎重に決めたいと思う
・今回の結果ではAdamとAdamWの収束性能がよさそうである

・途中でoptimizerを変更するのもポテンシャル形状によっては有効な場合がありそうである
・図にするとイメージはしやすいがそれぞれの特徴を見出すまでは至っていない
###おまけ

torch.optimizerのコード

import torch
from torch import optim
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D 
    
def func(x, y, r0,r1):
    return r0*x*x + r1*y*y #+x*y

#各データの格納リスト
classes = ['SGD', 'Adagrad', 'RMSprop', 'Adadelta', 'Adam', 'AdamW']
x_list = {}
y_list = {}
f_list = {}
optimizers = {}

#学習率（SGDのみ0.1としています。）
#lr = 1
Ir0= 0.01 #1 0.1 0.25 0.5 
r0 = 0.0005  #1 10 0.1
r1 = 1  #2 1.0.1 0.01 0.0 0.00000001 -0.00000001 0.0001
for key in classes:
    x_list[key] = []
    y_list[key] = []
    f_list[key] = []

#各種optimizer毎に関数の最適化を行う
for key in classes:
    print(key)
    #各optimizerの定義

    #更新パラメータ（x, y）の初期値
    x = torch.tensor(-7.0, requires_grad=True) #-75
    y = torch.tensor(2.0, requires_grad=True) #-10
    params = [x, y]

    optimizers['SGD'] = optim.SGD(params, lr=Ir0*0.1, momentum=0.0)
    optimizers['Adagrad'] = optim.Adagrad(params, lr=Ir0)
    optimizers['RMSprop'] = optim.RMSprop(params, lr=Ir0)
    optimizers['Adadelta'] = optim.Adadelta(params, lr=Ir0)
    optimizers['Adam'] = optim.Adam(params, lr=Ir0)
    optimizers['AdamW'] = optim.AdamW(params, lr=Ir0)

    #関数の出力に対して勾配を求め、optimizerで更新を繰り返す
    for i in range(10000):
        #print('Step : ', i+1)
        optimizers[key].zero_grad()
        outputs = func(x, y, r0, r1)
        outputs.backward()
        optimizers[key].step()
        x_list[key].append(x.item())
        y_list[key].append(y.item())
        f_list[key].append(outputs.item())

#一覧グラフ作成
fig = plt.figure(figsize=(15, 10))
i = 1
for key in classes:
    ax = fig.add_subplot(2, 3, i)
    ax.scatter(x_list[key], y_list[key])
    ax.plot(-7, 2, '+') #-75,-10
    if key =='SGD':
        Ir = 0.1*Ir0
    else:
        Ir = Ir0
    ax.set_title('{}_Ir={}_{}x^2+{}y^2'.format(key, Ir, r0,r1), fontsize=16)
    ax.set_xlabel('x', fontsize=16)
    ax.set_ylabel('y', fontsize=16)
    ax.set_xlim(-10, 10) #-80,80
    ax.set_ylim(-10, 10) #-15,15
    i += 1
 
plt.tight_layout()
plt.savefig('fig_optim_Ir{}_r0{}_r1{}.png'.format(Ir0,r0,r1))
plt.pause(1)
plt.close()

X = np.arange(-10, 10, 1) #-80, 15, 5  
Y = np.arange(-10, 10, 1) #-15, 15, 1
X, Y = np.meshgrid(X, Y)
Z = func(X,Y,r0,r1) #r0*X*X + r1*Y*Y
#print(Z)


fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
hs = ax1.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='jet') # cm.viridis
fig.colorbar(hs, ax=ax1)
ax1.set_title('potential_{}x^2+{}y^2'.format(r0,r1), fontsize=16)

ax2 = fig.add_subplot(1, 2, 2)
hv = ax2.pcolormesh(X,Y,Z, cmap='jet')
fig.colorbar(hv, ax=ax2)
plt.pause(1)
plt.savefig('orijinal_potential2_r0={}_r1={}.png'.format(r0,r1))
plt.close()

n=10000
for kin in range(6): #6
    fig = plt.figure(figsize=(12, 6))
    ax1 = fig.add_subplot(121 , projection='3d')
    ax2 = fig.add_subplot(1, 2, 2)
    hv = ax2.pcolormesh(X,Y,Z, cmap='jet')
    # 等高線図の生成
    cont=ax2.contour(X,Y,Z,  10, colors=['black']) #Vmax=1
    cont.clabel(fmt='%1.1f', fontsize=14)
    sc0_ = ax1.plot_wireframe(X, Y, Z, rstride=2, cstride=2)
    #sc1_ = ax1.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='jet') 
     
    key = classes[kin]
    j =0
    while j <n:
        if j<10:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c ='blue')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='blue')
            j +=1

        elif j < 1000:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c='red')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='red')
            j += 5
        elif j < 10000:
            zk = func(x_list[key][j], y_list[key][j], r0, r1)
            sc = ax1.scatter(x_list[key][j], y_list[key][j],zs=zk,zdir='z', s=20, c='black')
            sc1 = ax2.scatter(x_list[key][j], y_list[key][j], s=10, c ='white')
            j += 200

    if key =='SGD':
        Ir = 0.1*Ir0
    else:
        Ir = Ir0            
    ax1.set_title('{}_Ir={}_{}x^2+{}y^2'.format(key,Ir,r0,r1), fontsize=16)
    ax1.set_xlabel('x', fontsize=16)
    ax1.set_ylabel('y', fontsize=16)    
    ax1.set_zlabel('z', fontsize=16)    
    fig.colorbar(hv, ax=ax2)

    plt.savefig('fig3d_optim_{}_Ir{}_r0{}_r1{}.png'.format(key,Ir, r0,r1))
    plt.pause(3)

自前optimizerのコード

import matplotlib.pyplot as plt
import numpy as np

"""
n = 50
x = np.arange(n) * np.pi
y = np.cos(x) * np.exp(x/100) - 10*np.exp(-0.01*x)

def ewa(y, beta):
    #Exponentially weighted average.

    zs = np.zeros(len(y))
    z = 0
    for i in range(n):
        z = beta*z + (1 - beta)*y[i]
        zs[i] = z
    return zs
    
def ewabc(y, beta):
    #Exponentially weighted average with hias correction.

    zs = np.zeros(len(y))
    z = 0
    for i in range(n):
        z = beta*z + (1 - beta)*y[i]
        zc = z/(1 - beta**(i+1))
        zs[i] = zc
    return zs
    
beta = 0.9

plt.plot(x, y, 'o-')
plt.plot(x, ewa(y, beta), c='red', label='EWA')
plt.plot(x, ewabc(y, beta), c='orange', label='EWA with bias correction')
plt.legend()
pass

plt.show()
"""
def f(x):
    return x**2

def grad(x):
    return 2*x

def gd(x, grad, alpha, max_iter=10):
    xs = np.zeros(1 + max_iter)
    xs[0] = x
    for i in range(max_iter):
        x = x - alpha * grad(x)
        xs[i+1] = x
    return xs

def gd_momentum(x, grad, alpha, beta=0.9, max_iter=10):
    xs = np.zeros(1 + max_iter)
    xs[0] = x
    v = 0
    for i in range(max_iter):
        v = beta*v + (1-beta)*grad(x)
        vc = v/(1+beta**(i+1))
        x = x - alpha * vc
        xs[i+1] = x
    return xs

alpha = 0.1
x0 = 1
xs = gd(x0, grad, alpha)
xp = np.linspace(-1.2, 1.2, 100)
plt.title('fig_y=x^2_Ir{}'.format(alpha))
plt.plot(xp, f(xp))
plt.plot(xs, f(xs), 'o-', c='red')
for i, (x, y) in enumerate(zip(xs, f(xs)), 1):
    plt.text(x, y+0.2, i,
             bbox=dict(facecolor='yellow', alpha=0.5), fontsize=14)
pass
plt.savefig('fig_y=x^2_Ir{}.png'.format(alpha))
plt.show()

alpha = 0.95
xs = gd(1, grad, alpha)
xp = np.linspace(-1.2, 1.2, 100)
plt.title('fig_y=x^2_Ir{}'.format(alpha))
plt.plot(xp, f(xp))
plt.plot(xs, f(xs), 'o-', c='red')
for i, (x, y) in enumerate(zip(xs, f(xs)), 1):
    plt.text(x*1.2, y, i,
             bbox=dict(facecolor='yellow', alpha=0.5), fontsize=14)
pass
plt.savefig('fig_y=x^2_Ir{}.png'.format(alpha))
plt.show()

alpha = 0.95
xs = gd_momentum(1, grad, alpha, beta=0.9)
xp = np.linspace(-1.2, 1.2, 100)
plt.title('fig_y=x^2_Ir{} with momentum'.format(alpha))
plt.plot(xp, f(xp))
plt.plot(xs, f(xs), 'o-', c='red')
for i, (x, y) in enumerate(zip(xs, f(xs)), 1):
    plt.text(x, y+0.2, i,
             bbox=dict(facecolor='yellow', alpha=0.5), fontsize=14)
pass
plt.savefig('fig_y=x^2_Ir{}withmomentum.png'.format(alpha))
plt.show()

def f2(x):
    return x[0]**2 + 100*x[1]**2

def grad2(x):
    return np.array([2*x[0], 200*x[1]])

def grad3(x):
    return np.array([2*0.05*x[0], 2*x[1]])


x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
pass

plt.show()

def gd2(x, grad, alpha, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0,:] = x
    for i in range(max_iter):
        x = x - alpha * grad(x)
        xs[i+1,:] = x
    return xs

def gd2_momentum(x, grad, alpha, beta=0.9, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    for i in range(max_iter):
        v = beta*v + (1-beta)*grad(x)
        vc = v/(1+beta**(i+1))
        x = x - alpha * vc
        xs[i+1, :] = x
    return xs

alpha = 0.009 #0.01
x0 = np.array([-1,-1])
xs = gd2(x0, grad2, alpha, max_iter=75)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Vanilla gradient descent')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withVGD.png'.format(alpha))
plt.show()

def f3(x):
    return 0.05*x[0]**2 + x[1]**2

alpha = 0.9 #0.01
x0 = np.array([-1,-1])
xs = gd2(x0, grad3, alpha, max_iter=75)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.01,0.1,0.2,0.4,0.9, 1.6, 2.5, 3.6, 4.9, 6.4, 8.1, 10.0]
Z = 0.05*x**2 + Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Vanilla gradient descent_2')
pass
plt.savefig('fig_y=0.05x^2+y^2_Ir{}withVGD.png'.format(alpha))

plt.show()

alpha = 0.75 #0.01
x0 = np.array([-1,-1])
xs = gd2(x0, grad3, alpha, max_iter=75)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.01,0.1,0.2,0.4,0.9, 1.6, 2.5, 3.6, 4.9, 6.4, 8.1, 10.0]
Z = 0.05*x**2 + Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Vanilla gradient descent_2')
pass
plt.savefig('fig_y=0.05x^2+y^2_Ir{}withVGD.png'.format(alpha))
plt.show()

alpha = 0.3 #0.01
x0 = np.array([-1,-1])
xs = gd2(x0, grad3, alpha, max_iter=75)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.01,0.1,0.2,0.4,0.9, 1.6, 2.5, 3.6, 4.9, 6.4, 8.1, 10.0]
Z = 0.05*x**2 + Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Vanilla gradient descent_2')
pass
plt.savefig('fig_y=0.05x^2+y^2_Ir{}withVGD.png'.format(alpha))
plt.show()

alpha = 0.009 #0.01
x0 = np.array([-1,-1])
xs = gd2_momentum(x0, grad2, alpha, beta=0.9, max_iter=75)


x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradieent descent with momentum')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withmomentum.png'.format(alpha))
plt.show()

def gd2_rmsprop(x, grad, alpha, beta=0.9, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    for i in range(max_iter):
        v = beta*v + (1-beta)*grad(x)**2
        x = x - alpha * grad(x) / (eps + np.sqrt(v))
        xs[i+1, :] = x
    return xs

alpha = 0.1
x0 = np.array([-1,-1])
xs = gd2_rmsprop(x0, grad2, alpha, beta=0.9, max_iter=10)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradient descent with RMSprop')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withRMSprop.png'.format(alpha))
plt.show()

def gd2_adam(x, grad, alpha, beta1=0.9, beta2=0.999, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    m = 0
    v = 0
    for i in range(max_iter):
        m = beta1*m + (1-beta1)*grad(x)
        v = beta2*v + (1-beta2)*grad(x)**2
        mc = m/(1+beta1**(i+1))
        vc = v/(1+beta2**(i+1))
        x = x - alpha * m / (eps + np.sqrt(vc))
        xs[i+1, :] = x
    return xs

alpha = 0.1
x0 = np.array([-1,-1])
xs = gd2_adam(x0, grad2, alpha, beta1=0.9, beta2=0.9, max_iter=10)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradient descent with Adam')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withAdam.png'.format(alpha))
plt.show()

def gd2_adagrad(x, grad, alpha, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    for i in range(max_iter):
        v = v + grad(x)**2
        x = x - alpha * grad(x) / (eps + np.sqrt(v))
        xs[i+1, :] = x
    return xs

alpha = 0.1
x0 = np.array([-1,-1])
xs = gd2_adagrad(x0, grad2, alpha, max_iter=10)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradient descent with Adagrad')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withAdagrad.png'.format(alpha))

plt.show()

def gd2_adadelta(x, grad, alpha, beta1=0.9, beta2 = 0.999, beta3 = 0.999, eps=1e-8, max_iter=10):
    xs = np.zeros((1 + max_iter, x.shape[0]))
    xs[0, :] = x
    v = 0
    h = 0
    #dw = -alpha*grad(x) / (eps + np.sqrt(h))
    for i in range(max_iter):
        h = beta1*h + (1-beta1)*grad(x)*grad(x)
        dw = -alpha*grad(x) / (eps + np.sqrt(h))
        dw = beta2*dw + (1-beta2)*dw*dw
        v = beta3*v + (1-beta3)*grad(x)**2
        x = x - np.sqrt(eps + dw) * grad(x) / np.sqrt(eps + v)
        xs[i+1, :] = x
    return xs

alpha = 0.01 #0.001 #0.01 #0.1発散
x0 = np.array([-1,-1])
xs = gd2_adadelta(x0, grad2, alpha, beta1=0.9, beta2=0.9, beta3 =0.9, max_iter=10)

x = np.linspace(-1.2, 1.2, 100)
y = np.linspace(-1.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
levels = [0.1,1,2,4,9, 16, 25, 36, 49, 64, 81, 100]
Z = x**2 + 100*Y**2
c = plt.contour(X, Y, Z, levels)
plt.plot(xs[:, 0], xs[:, 1], 'o-', c='red')
plt.title('Gradient descent with Adadelta')
pass
plt.savefig('fig_y=x^2+100y^2_Ir{}withAdadelta.png'.format(alpha))

plt.show()