最急降下法とニュートン法による最適化

# はじめに
最急降下法は、名前のとおり曲線の極小値の探索に用いられる。具体的には、曲線の微分の傾きが小さくなる方に値を移動させるものである。一方で、ニュートンラフソン法（以下ニュートン法）では、極小値では微分係数が0となることを利用する。具体的には、そのような方程式をまずは解くことで極小値を求める。したがって、最急降下法では1階微分を必要として、ニュートン法では2階微分を必要とする。そこで、今回は2つのアルゴリズムについて紹介してその収束度について調査することを目的とする。

以下、最急降下法での結果である。

![二次関数の最小値探索.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/a85e52bd-e3b3-d631-3292-4724b430970c.png)

![誤差関数と積算回数の関係.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/6d2ebf7c-95a1-6aeb-b3b6-3bc8eec4c71f.png)


# アルゴリズム
## 最急降下法
$z=f(x,y)$の最小値を求めることを考える。
解空間の最小値を取る極値が一つ解空間の谷が一つの場合に有効な最小値の推定手法である。$(a,b)$に初期値を与え、以下のような式を用いて$(a,b)$を更新していく。つまり、$(a_k,b_k)$を算出する。


``` math
a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a_{i}}f(a_{k},b_{k}))
```

``` math
b_{k+1}=b_{k}-\beta(\frac{\delta }{\delta b_{k}}f(a_k,b_k))
```
ただし$\alpha,\beta$は$(a,b)$の学習率である。

これは、以下のグラフのようにボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。

## ニュートン法

ニュートン法は、最急降下法を発展させたもので、一次元での場合は以下のようなアルゴリズムで値を更新する。

例.$y=f(x)$の最小値を求めたい場合

```math
x_{n+1}=x_n-\frac{f''(x)}{f'(x)}
```

これを3次元に拡張する。
$z=f(x,y)$とするとその関数の最小値を求める。
```math
\textbf{x}_{n+1}=\textbf{x}_n-\textbf{H}_f^{-1}\nabla f
```
ただし、

```math
\textbf{H}_f=

\begin{pmatrix}
   \frac{\delta^2 f}{\delta x^2} & \frac{\delta^2 f}{\delta xy}  \\
   \frac{\delta^2 f}{\delta yx} & \frac{\delta^2 f}{\delta y^2} 
\end{pmatrix}

```

# プログラム
## 前提
ベンチマークの関数として以下の2次関数を用いた。

``` math
f(x,y)=(x-2)^2+(y-2)^2+1
```

これの最小値は$(x=2,2)$のとき$1$となる。

この関数が最小になる条件を探索する。

## 最急降下法
``` python gradient_descent.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

n=100
x=np.linspace(-2,2,n)
y=np.linspace(-2,2,n)

n=100
L=10
x=np.linspace(-L,L,n)
y=np.linspace(-L,L,n)
X,Y=np.meshgrid(x,y)
Z=np.zeros((n,n))
h=1.0*10**-5

def fun(a,b):
  function_minimize=(a-2)**2+(b-2)**2+1
  return function_minimize

for i in range(n):
  for j in range(n):
    Z[i][j]=fun(X[i][j],Y[i][j])
#等高線の表示(二次関数の描写)
# plt.contourf(X,Y,Z,cmap="jet",alpha=0.8)
# plt.colorbar(label="z")
# plt.xlabel('x')
# plt.ylabel('y')
# plt.savefig("gauss1.png")
# plt.show()

iterations=100
h=1.0*10**-5
#初期値
a=10
b=10

#最急降下法
#学習率
alpha=0.2
beta=0.2

a_pre=[]
b_pre=[]
fun_pre=[]
i_ary=[]

for i in range(iterations):
  i_ary.append(i+1)
  d_err_a=(fun(a+h,b)-fun(a,b))/h
  d_err_b=(fun(a,b+h)-fun(a,b))/h
  a=a-alpha*d_err_a
  b=b-beta*d_err_b
  a_pre.append(a)
  b_pre.append(b)
  fun_pre.append(fun(a,b))
#誤差関数と積算回数の関係
# plt.plot(i_ary,fun_pre)
# plt.xlabel('積算回数')
# plt.ylabel('二次関数の値')
# plt.savefig("誤差関数と積算回数の関係.png")
# plt.show()

#等高線の表示
plt.contourf(X,Y,Z,cmap="jet",alpha=0.5)
plt.scatter(a_pre,b_pre,c=fun_pre,cmap="jet")
plt.colorbar(label="z")
plt.xlabel('x')
plt.ylabel('y')
plt.savefig("二次関数の最小値探索.png")
plt.show()
```



## ニュートン法

``` python newton.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

n=100
L=10
x=np.linspace(-L,L,n)
y=np.linspace(-L,L,n)
X,Y=np.meshgrid(x,y)
Z=np.zeros((n,n))
h=1.0*10**-5

def fun(a,b):

  function_minimize=(a-2)**2+(b-2)**2+1
  return function_minimize

for i in range(n):
  for j in range(n):
    Z[i][j]=fun(X[i][j],Y[i][j])

def fun_dif(a,b):
  d_err_a=(fun(a+h,b)-fun(a,b))/h
  d_err_b=(fun(a,b+h)-fun(a,b))/h
  dif_f_vec=np.array([d_err_a,d_err_b])
  return dif_f_vec
 
def fun_ddif(a,b):
  dd_err_a=(fun_dif(a+h,b)-fun_dif(a,b))/h
  dd_err_b=(fun_dif(a,b+h)-fun_dif(a,b))/h
  J=np.array([dd_err_a,dd_err_b])
  return J

iterations=10
#初期値
a=10
b=10

#ニュートンラフソン法
a_pre=[]
b_pre=[]
fun_pre=[]
i_ary=[]

XX=np.array([a,b])

for i in range(iterations):
  i_ary.append(i+1)
  J=fun_ddif(a,b)
  J_inv=np.linalg.inv(J)
  d_err_a,d_err_b=fun_dif(a,b)

  a=XX[0]
  b=XX[1]
  a_pre.append(a)
  b_pre.append(b)
  fun_pre.append(fun(a,b))
  XX=XX-J_inv@fun_dif(a,b)
  
plt.plot(i_ary,fun_pre)
plt.xlabel('積算回数')
plt.ylabel('二次関数の値')
plt.savefig("誤差関数と積算回数の関係_ニュートンラフソン法.png")
plt.show()

#等高線の表示
# plt.contourf(X,Y,Z,cmap="jet",alpha=0.5)
# plt.scatter(a_pre,b_pre,c=fun_pre,cmap="jet")
# plt.colorbar(label="z")
# plt.xlabel('x')
# plt.ylabel('y')
# plt.savefig("二次関数の最小値探索.png")
# plt.show()

print(a,b)
```
# 結果
上記のプログラムを実行すると以下のようなグラフが出力される。（コメント文を消したりするとできる）

## 最急降下法
最小値の探索状態について以下のようになった。
![二次関数の最小値探索.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/a85e52bd-e3b3-d631-3292-4724b430970c.png)

このように、最小値に向かってプロット（値）が移動しているが、収束に時間がかかっているようである。

![誤差関数と積算回数の関係.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/6d2ebf7c-95a1-6aeb-b3b6-3bc8eec4c71f.png)


## ニュートン法
最小値の探索状態について以下のようになった。

![二次関数の最小値探索.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/8c3943f1-c1cc-7aeb-617b-c18b47bc77ce.png)


![誤差関数と積算回数の関係_ニュートンラフソン法.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/2e618a34-cb2c-596d-06ec-6a5f77fc6a66.png)

最急降下法と比較して、10倍程度早く収束するということが分かった。


# まとめ

今回は、単純な二次関数の最適化を最急降下法とニュートン法を用いることで行い、収束速度の比較を行った。結果、ニュートン法の方が、計算の収束が速くなること可能性があること分かった。しかし、ニュートン法では、二階微分を計算する場合などにおいて値が発散してしまい上手く計算することができなくなることもあることから、適切な初期値を用いることが大切である。