Python NumPy サンプルコード: Newton 法の直線探索

今回は制約なし最適化問題の数値計算法でメジャーな Newton 法 の直線探索を紹介します。

Newton 法の直線探索

制約なし最適化問題の数値計算法はいくつかありますが、数式の情報も使って素早く最適解を得るのが Newton 法です。 機械屋さんはエネルギーの最小化計算とかでよく使っています。

最適化の関数は SciPy にもたくさん用意されていますが、簡単なので初心者は勉強がてら自前で定義して使ってみては如何がでしょうか。 以下のプログラムでは Newton 法では勾配ベクトルとヘッセ行列から次の探索方向を計算しているのが分かるかと思います。 詳しい理屈の説明は Web にもいっぱい情報があるのでそちらに譲りますが、使い方は目的関数と、その勾配とヘッセ行列を定義してやるだけなので難しくありません。

では Matlab のドキュメント の例にあわせて、関数 \(f(\boldsymbol{x})=3x_1^2+2x_1x_2+x_2^2\) を最小化する例を示します。 収束解は単純に \(\boldsymbol{x} = (0, 0)\) で \(0\) です。

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 gradientHessian(x):
    # 勾配ベクトルとヘッセ行列
    g = [6*x[0] + 2*x[1], 2*x[0] + 2*x[1]]
    h = [[6, 2],
         [2, 2]]
    return g, h


def newtonLineSearch(function, ghFunction, x, tolerance):
    # Newton 法の直線探索
    residual = 1e12
    fOld = 1e12
    alpha = 1
    while residual > tolerance:
        f = function(x)
        g, h = ghFunction(x)
        delta = -np.linalg.solve(h, g)
        x = x + alpha*delta
        residual = abs(f-fOld)
        fOld = f
    return x, f


x0 = [10.0, 10.0]
x, f = newtonLineSearch(objectiveFunction, gradientHessian, x0, 1e-12)
print(x, f)
[ 0.  0.] 0.0

こんな風に簡単に最適解が求められるのは本当に便利です。

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編