最急降下法(gradient descent)

勾配法の下るバージョンの最急降下法(gradient descent)を実装してみた. ちなみに上るバージョンは最急上昇法(gradient ascent)と呼ばれる. 最近書いている微分を近似するという内容の記事は実はこれの準備のためなのでした.笑
導関数や勾配を明示的に与えることの無い実装を,したくてしたくてたまらなかったのです.
このアルゴリズムにより極小値を求めることができる. 通常勾配を明示的に与える必要があるが,最近書いている勾配の近似で代用した. つまり明示的に勾配を与える必要がない.
実装は以下.

 1: # f(x1, ..., xn)のxiによる偏導関数を返す
 2: def partial(f, n, i, dx=0.001):
 3:     I = np.identity(n)
 4:     def part(x):
 5:         return (f(x+I[i]*dx) - f(x)) / dx
 6:         #return (f(x+I[i]*2*dx) - f(x)) / 2*dx
 7:     return part
 8: 
 9: # f(x1, ..., xn)の勾配を返す
10: def Gradient(f, n, dx=0.01):
11:     def g(x):
12:         return np.array([G[i](x) for i in np.arange(n)], dtype=float)
13:     G = [partial(f, n, i, dx=dx) for i in np.arange(n)]
14:     return g
15: 
16: # n変数関数fを最小にするxを勾配法により求める. 
17: def gradient_descent(f, n, h0=0.03, dx=1./10**2):
18:     fx = Gradient(f, n, dx=dx)
19:     pp = np.random.random_sample(n)
20:     g = fx(pp)
21:     p = pp - h0*g
22:     h = h0
23:     while (p-pp).dot(p-pp) > 1./10**12:
24:         ph = h; h = 2 * ph
25:         while f(p-h*g) < f(p-ph*g):
26:             ph = h; h *= 2
27:         h /= 2; pp = p
28:         g = fx(p); p = p - h*g; h = h0
29:     return p

実は,この本を読んでいます.

これなら分かる最適化数学―基礎原理から計算手法まで

この本の例題3.2に出てくる関数

\displaystyle  f(x,y) = x^{3} + y^{3} - 9xy + 27

に上記のアルゴリズムを適用します.

n = 2
f = lambda x: x[0]**3 + x[1]**3 - 9*x[0]*x[1] + 27
solution = gradient_descent(f, n)
print solution, f(solution)

実行結果は以下のようになりました.

[ 2.98999076  2.98999076] 0.000899659139463

解は(x,y) = (3, 3)なので,良い感じだと思われます.
この例題はニュートン法の問題ですが,実はニュートン法も実装してみたのでまた記事を書きたいと思います. いやあ,楽しいなぁ.

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中