最急降下法(Gradient Descent)とステップ幅の直線探索

勾配ベクトルを用いた漸化式的に表される反復法を用いて最適解を計算する手法を最急降下法(Gradient Descent)といいます。当記事では最急降下法の数式と、ステップ幅の適応的計算にあたって用いられる直線探索について取りまとめを行いました。
「数理計画入門(朝倉書店)」の$4.4$節の「最急降下法」やWikipediaの内容を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

新版　数理計画入門

福島雅夫

3,520円(10/04 09:41時点)

Amazon

数式の確認

最急降下法

最急降下法の数式

最急降下法の漸化式は下記のように表される。
$$
\large
\begin{align}
\mathbf{x}_{k+1} = \mathbf{x}_{k} – \alpha_{k} \nabla f(\mathbf{x}_{k}) \quad (1)
\end{align}
$$

最急降下法の解釈

$(1)$式における$\nabla f(\mathbf{x}_{k})$は$\mathbf{x}=\mathbf{x}_{k}$における勾配ベクトルであり、最急降下法は勾配の計算結果に基づいてベクトル$\mathbf{x}$の位置を更新すると解釈できる。

ステップ幅の決定

前節$(1)$式の$\alpha_{k}$の決定にあたっては主に「①定数を用いる方法」と、「②適応的計算を用いる方法」がある。以下、「定数を用いる方法」と都度計算を行う「適応的計算」の$1$例である「直線探索」について取りまとめを行なった。

定数を用いる

$(1)$式の$\alpha_{k}$には定数$\alpha$を用いることが多い。定数$\alpha$の値が大きい場合は$\mathbf{x}$が発散し、小さい場合は収束までに時間がかかることは抑えておくと良い。

直線探索

直線探索は「勾配ベクトルの方向に沿って目的関数が最小となるようなステップ幅$\alpha_{k}$を用いる手法」である。

Pythonを用いた最急降下法の実行

問題設定

「数理計画入門(朝倉書店)」の$4.4$節の「最急降下法」と同様に下記の関数の最適化にあたって、最急降下法のプログラムの作成を行う。
$$
\large
\begin{align}
f(\mathbf{x}) = f(x_1,x_2) = (x_1-1)^2 + 10(x_1^2-x_2)^2
\end{align}
$$

勾配ベクトルの計算

勾配ベクトル$\nabla f(\mathbf{x})$は下記のように計算できる。
$$
\large
\begin{align}
\nabla f(\mathbf{x}) &= \left(\begin{array}{c} \displaystyle \frac{\partial f}{\partial x_1} \\ \displaystyle \frac{\partial f}{\partial x_2} \end{array} \right) \\
&= \left(\begin{array}{c} 2(x_1-1) + 20(x_1^2-x_2) \cdot 2x_1 \\ 20(x_1^2-x_2) \cdot (-1) \end{array} \right) \\
&= \left(\begin{array}{c} 40x_1^3-40x_1x_2+2x_1-2 \\ -20(x_1^2-x_2) \end{array} \right)
\end{align}
$$

最急降下法の実行

下記を実行することで最急降下法による最適化を行うことができる。

import numpy as np

def calc_func_value(x1,x2):
    return (x1-1)**2 + 10*(x1**2-x2)**2

def calc_gradient(x1,x2):
    return np.array([40*x1**3 - 40*x1*x2 + 2*x1 - 2, -20*(x1**2-x2)])

x = np.array([0., 1.])
gamma = np.arange(-0.5, 0.5, 0.0000001)

display_flag = [1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500]

for k in range(500):
    grad = calc_gradient(x[0],x[1])
    x1 = x[0]+gamma*grad[0]
    x2 = x[1]+gamma*grad[1]
    f = calc_func_value(x1,x2)
    x[0] = x1[np.argmin(f)]
    x[1] = x2[np.argmin(f)]
    if k+1 in display_flag:
        print("Iter: {:.0f}, x: {}".format(k+1,x))

・実行結果

Iter: 1, x: [ 0.0998848  0.001152 ]
Iter: 2, x: [ 0.36070906  0.02723477]
Iter: 3, x: [ 0.35167097  0.11761506]
Iter: 4, x: [ 0.44424751  0.12687296]
Iter: 5, x: [ 0.43824528  0.18689394]
Iter: 10, x: [ 0.57217842  0.28642287]
Iter: 20, x: [ 0.6772309   0.43291189]
Iter: 30, x: [ 0.7394933   0.52795907]
Iter: 40, x: [ 0.78289396  0.59811139]
Iter: 50, x: [ 0.81554999  0.65307488]
Iter: 100, x: [ 0.90619007  0.81570209]
Iter: 200, x: [ 0.96896765  0.93720982]
Iter: 300, x: [ 0.98869242  0.97691039]
Iter: 400, x: [ 0.99575609  0.99130582]
Iter: 500, x: [ 0.9983905   0.99669874]

上記の結果は「数理計画入門(朝倉書店)」の表$4.1$の結果に基本的には一致する。「数理計画入門(朝倉書店)」に直線探索の設定がないので、厳密な再現ではないことに注意が必要である。

参考

新版　数理計画入門

福島雅夫

3,520円(10/04 09:41時点)

Amazon