Python NumPy サンプルコード: Newton 法の黄金探索法によるステップ幅の最適化

昨日 に続き、制約なし最適化問題の数値計算法である Newton 法についての投稿です。
今回は、ステップ幅を 黄金探索法 で更新する方法を紹介します。

黄金探索法によるステップ幅の最適化

直線探索においては、解の探索能力、収束速度はステップ幅 \(\alpha\) に強く依存します。
\(\alpha\) は小さな正の定数( \(0 \lt \alpha \leq {1}\) )ですが、\(\alpha\) が大きすぎると発散することがあり、逆に小さすぎると収束が遅くなってしまいます。

通常の Newton 法では \(\alpha=1\) なので収束が速いですが、場合によっては発散することがあるので、 黄金探索法で \(\alpha\) を最適化しながら収束解を探索したりします。
計算時間が増えてしまうので \(\alpha=1\) で十分収束するなら使う必要はありません。
黄金探索法の実装としては以下のようになります。

# coding: utf-8
import numpy as np


def objectiveFunction(x):
    # 目的関数
    # f(x,y) = 3*x1^2 + 2*x1*x2 + x2^2
    f = 3*x[0]**2 + 2*x[0]*x[1] + x[1]**2
    return f


def gradientHesssian(x):
    # 勾配ベクトルとヘッセ行列
    g = [6*x[0] + 2*x[1], 2*x[0] + 2*x[1]]
    h = [[6, 2],
	 [2, 2]]
    return g, h


def goldenRatioSearch(function, range, tolerance):
    # 黄金探索法によるステップ幅の最適化
    gamma = (-1+np.sqrt(5))/2
    a = range[0]
    b = range[1]
    p = b-gamma*(b-a)
    q = a+gamma*(b-a)
    Fp = function(p)
    Fq = function(q)
    width = 1e8
    while width > tolerance:
	if Fp <= Fq:
	    b = q
	    q = p
	    Fq = Fp
	    p = b-gamma*(b-a)
	    Fp = function(p)
	else:
	    a = p
	    p = q
	    Fp = Fq
	    q = a+gamma*(b-a)
	    Fq = function(q)
	 width = abs(b-a)/2
    alpha = (a+b)/2
    return alpha


def newtonLineSearch(function, ghFunction, x, tolerance, flag):
    # Newton 法の直線探索
    residual = 1e12
    fOld = 1e12
    while residual > tolerance:
	f = function(x)
	g, h = ghFunction(x)
	delta = -np.linalg.solve(h, g)
	if flag == 1:
	    alpha = goldenRatioSearch(lambda alpha: function(x+alpha*delta), [1e-2, 1.0], 1e-2)
	else:
	    alpha = 1
	x = x + alpha*delta
	residual = abs(f-fOld)
	fOld = f
    return x, f


x0 = [10.0, 10.0]
x, f = newtonLineSearch(objectiveFunction, gradientHesssian, x0, 1e-12, 1)
print(x, f)

コメント

Comments powered by Disqus
書籍更新情報
2017-01-16
Pythonによる科学技術計算 基礎編
1.3版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストアで販売中です
Pythonによる科学技術計算 基礎編
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share