Python 上昇勾配法を使って極値を計算する

上昇勾配法とは

滑らかな曲線 $y=f(x)$ に対して、極値を求める方法。

ステップ1. ユーザーは初期値、曲線を与える
ステップ2. 初期値の点に対して曲線の傾きを計算し、初期値が曲線の上昇している方向（x軸方向）へ移動するような処理を行う。
ステップ3. ステップ2で初期値が移動した点に対して次のステップ4の評価を行い、再度ステップ2を行う。
ステップ4. 移動した距離が十分小さくなったとき、処理を終了する。

極値付近では曲線の傾きが非常に小さい。このときステップ2で移動する点の距離は非常に小さくなり、処理が終了する。

ステップ2については以下のような処理である。
$x_0$ を初期値とする。$x_1$ を次のように定義する。
$x_1 = x_0 + \theta \cdot f'(x_0)$
ここで $\theta$ は微小な量である。
これで $x_1$ は曲線の上昇しているx軸の方向へ移動する。
一般的に定義すれば次のようになる。
$x_{i + 1} = x_i + \theta \cdot f'(x_i) \ \ (i = 0, 1, \cdots)$

サンプルコード

epsilonが、ステップ4で移動した距離が十分小さいかどうかを判定する微小な量である。

step_sizeが、ステップ2で初期値を曲線の上昇している方向へ移動するために用いる微小な量である。

from sympy import *
import numpy as np
import matplotlib.pyplot as plt

def gradient(initial_value, f):
    # 変数を指定
    x = Symbol('x')
    # fを微分する
    f1 = diff(f, x)
    # 微小な量を定義する
    epsilon = 1e-6
    step_size = 1e-2
    x_old = initial_value
    x_new = x_old + step_size * f1.subs(x, x_old)
    cnt = 0
    # 変化を追うための各点を格納する配列
    x_array = np.empty(0)
    y_array = np.empty(0)
    # ステップ2とステップ4の処理
    while abs(x_old - x_new) > epsilon:
        x_array = np.append(x_array, x_old)
        y_array = np.append(y_array, f.subs(x, x_old))
        x_old = x_new
        x_new = x_old + step_size * f1.subs(x, x_old)
    # プロットする区間を定義する
    x_interval = np.linspace(-2, 2, 100)
    y_images = np.empty(0)
    # プロットする区間の各点xに対応するf(x)を計算する
    for x_val in x_interval:
        y_images = np.append(y_images, f.subs(x, x_val))
    # (x, f(x))をプロット
    plt.plot(x_interval, y_images, label='curve')
    # 上昇勾配法による変化をプロット
    plt.scatter(x_array, y_array, marker='o', s=50, color='lime')
    # 初期位置をプロット
    plt.scatter(
        [initial_value], [f.subs(x, initial_value)],
        label='initial position', marker='o', s=200, color='blue'
    )
    # 停止位置をプロット
    plt.scatter(
        [x_new], [f.subs(x, x_new)],
        label='stop position', marker='o', s=200, color='red'
    )
    print('勾配法で求めた極大値\nx = {} のとき極大値 {}\n'.format(
        x_new, f.subs(x, x_new)
    ))

    
def cal_extremes(curve, length):
    """
    微分を使って極値を計算するメソッド
    """
    x = Symbol('x')
    epsilon = 1e-10
    critical_points = solve(diff(curve, x))
    critical_points = [re(point) for point in critical_points
                       if abs(im(point)) < epsilon]
    for critical_point in critical_points:
        if im(critical_point) != 0\
           or length < critical_point or critical_point < -length:
            continue
        extreme_value = curve.subs(x, critical_point)
        diff1 = diff(curve, x)
        diff2 = diff(diff1, x)
        if diff2.subs(x, critical_point) < 0:
            print('微分で求めた極大値\nx = {} のとき極大値 {}\n'.format(
                critical_point.evalf(), extreme_value.evalf()
            ))


plt.figure(figsize=(14.0, 12.0))
initial_value = 0
x = Symbol('x')
f = x ** 3 - 2 * x
gradient(initial_value, f)
cal_extremes(f, 5)
plt.tick_params(labelsize=20)
plt.rc('legend', fontsize=15)
plt.legend()
plt.show()

実行結果

勾配法で求めた極大値
x = -0.816477683180498 のとき極大値 1.08866210702887

微分で求めた極大値
x = -0.816496580927726 のとき極大値 1.08866210790363

降下勾配法とは

上昇勾配法とは逆に、初期値が曲線の降下している方向へ移動するのが降下勾配法である。ステップ2については以下のように定義される。
$x_{i + 1} = x_i - \theta \cdot f'(x_i) \ \ (i = 0, 1, \cdots)$

次のように符号をマイナスにすればよい。

x_new = x_old - step_size * f1.subs(x, x_old)

失敗例

勾配法はステップ2で移動する距離が大きいと失敗する。移動する距離が大きく、目指している極値のx座標を超えてしまう。以下は失敗したときのグラフである。

解決策は
$x_{i + 1} = x_i + \theta \cdot f'(x_i) \ \ (i = 0, 1, \cdots)$
における $\theta$ を小さくすることである。

• 失敗時

step_size = 1e-2

この値をもっと小さくする。

step_size = 1e-4

参考文献

Amit Saha著 黒川 利明 訳 2016年05月 Pythonからはじめる数学入門