勾配法を用いた最適化計算

テキスト「これなら分かる最適化数学」の81ページにある勾配法(gradient method)アルゴリズムをpythonで実装した。勾配法とは、関数f(x)を最大にするxの値に近いと思われる値を初期値とし、必ず関数値f(x)が増加するようにし、かつステップ幅をなるべく大きくなるようその幅を調整することにより、できるだけ高速に関数を最大にするxを求めるというものである。

なお、今回扱う勾配法においては、変数の値がとりうる範囲について、関数の微分係数が0になるのは1箇所であると仮定する(複数極値があると、どこに向かうか分からないため)。また、微分できない箇所がある関数では不都合が生じるので、これも考えないことにする(例えばf(x)=|x|のx=0の点のように、尖っている状況下などがある)。

関数の最大値を求めるにあたっては、微分係数が0になるという性質を利用している。このため、微分という操作を行う必要であるが、この問題はsympyモジュールを利用することによって解決できる。sympyモジュールはこちらからダウンロードできる。なお、記事執筆時点のsympyのバージョンは0.7.1である。

以下、「2次関数f(x)」と「正規分布の確率変数と関係がある部分のみ取り出した関数g(x)」の2つの例を用いて、勾配法で正しく関数を最大化するようなxを推定できるか試してみた。結果としては、(制約が強いこともあるが)ほぼ正確に推定できているということになった。

# -*- encoding: utf-8 -*-

from sympy import *

# 符号関数
def sign(value):
    if value > 0: return 1
    if value < 0: return -1
    return 0

# 勾配法を計算するメインルーチン
# initは初期推定値、funcは最大値を求めたい関数、stepは初期移動量、
# epsは収束した(ほぼf'(x)が0となった)と判定する閾値である。
# 返り値は最大値となるxの値に相当する。
def gradient_method(init, func, step = 10**-7, eps = 10**-6):
    # diffは微分の関数、これで第一引数の関数をxで微分したことになる
    # 事前にxをSymbol関数でシンボル化しておく必要がある
    df = diff(func, x)
    xx, h = init, step
    while True:
        h = sign(limit(df, x, xx)) * abs(h)
        X, Xdash = xx, xx + h
        # 関数に値を入れるという方法が見つからないので
        # 極限の値を求めている
        if limit(func, x, X) < limit(func, x, Xdash):
            while limit(func, x, X) < limit(func, x, Xdash):
                h *= 2.0
                X, Xdash = Xdash, X + h
            xx, h = X, h / 2.0
        else:
            while limit(func, x, X) > limit(func, x, Xdash):
                h /= 2.0
                Xdash -= h
            xx, h = Xdash, 2.0 * h
        if abs(limit(df, x, xx)) <= eps:
            break
    return xx

x = Symbol('x') # xを記号として扱う
f = -x**2 + 2*x + 3
g = exp(-(x-2)**2)
print gradient_method(0, f) # 1.0000002、なお正解は1である
print gradient_method(0.5, g) # 2.0000002、なお正解は2である