読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

準ニュートン法による最適化Pythonサンプルコード

f:id:meison_amsl:20160503211305p:plain

目次

はじめに

以前、

ニュートン法による最適化サンプルコードの記事を書きました。

myenigma.hatenablog.com

 

上記のニュートン法は、

収束が非常に早いという特徴がありますが、

ヘッセ行列を計算しないといけないという問題があります。

 

複雑なモデルの場合、

各変数で二次微分した

ヘッセ行列を計算するのは

数式が複雑で難しいという問題があります。

 

そこで、

最適化手法の一つである、

ヘッセ行列を近似的に逐次計算する

準ニュートン法を使うことで、

ヘッセ行列が計算できないモデルにおいても、

高速な最適化が可能になります。

 

今回の記事では、

この準ニュートン法の概要と、

Pythonのサンプルコードを紹介したいと思います。

 

準ニュートン法

準ニュートン法は、

ニュートン法で使用するヘッセ行列を、

逐次的に推定するアルゴリズムです。

 

更新則としては、

ニュートン法と同じ下記の式を利用します。

ニュートン法との違いは、

上記のHkをヤコビ行列と最適化の履歴から、

逐次的に推定することです。  

Hkの推定方法はいくつか方式があります。

まず、下記のように勾配ベクトルの差をykとすると、

Hkの推定方法は下記のように複数ありますので、

いずれかを利用しましょう。

f:id:meison_amsl:20160503095647p:plain

 

いずれの方法も、

ヘッセ行列そのものは不要で、

Δxとyk、そして一つ前のHk-1のみで、

擬似ヘッセ行列を計算できます。

 

Pythonサンプルコード

下記が準ニュートン法による最適化Pythonサンプルコードです。

今回はHkの推定方法として、BSGS法を使いました。

(この方法を使った深い意味はありません)

#!/usr/bin/python
# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np
import random
import math

delta = 0.1
minXY=-5.0
maxXY=5.0
nContour=50
alpha=0.001

def Jacob(state):
    u"""
    jacobi matrix of Himmelblau's function
    """
    x=state[0,0]
    y=state[0,1]
    dx=4*x**3+4*x*y-44*x+2*x+2*y**2-14
    dy=2*x**2+4*x*y+4*y**3-26*y-22
    J=np.matrix([dx,dy]).T
    return J

def HimmelblauFunction(x,y):
    u"""
    Himmelblau's function
    see Himmelblau's function - Wikipedia, the free encyclopedia 
    http://en.wikipedia.org/wiki/Himmelblau%27s_function
    """
    return (x**2+y-11)**2+(x+y**2-7)**2

def CreateMeshData():
    x = np.arange(minXY, maxXY, delta)
    y = np.arange(minXY, maxXY, delta)
    X, Y = np.meshgrid(x, y)
    Z=[HimmelblauFunction(x,y) for (x,y) in zip(X,Y)]
    return(X,Y,Z)

def QuasiNewtonMethod(start,Jacob):
    u"""
    Quasi Newton Method Optimization
    """

    result=start
    x=start

    H= np.identity(2)
    preJ=None
    preG=None

    while 1:
        J=Jacob(x)

        sumJ=abs(np.sum(J))
        if sumJ<=0.01:
            print("OK")
            break

        grad=-np.linalg.inv(H)*J
        x+=alpha*grad.T
        
        result=np.vstack((result,np.array(x)))

        if preJ is not None:
            y=J-preJ
            H=H+(y*y.T)/(y.T*preG)-(H*preG*preG.T*H)/(preG.T*H*preG)

        preJ=J
        preG=(alpha*grad.T).T

    return result

# Main
start=np.matrix([random.uniform(minXY,maxXY),random.uniform(minXY,maxXY)])

result=QuasiNewtonMethod(start,Jacob)
(X,Y,Z)=CreateMeshData()
CS = plt.contour(X, Y, Z,nContour)

plt.plot(start[0,0],start[0,1],"xr");

optX=result[:,0]
optY=result[:,1]
plt.plot(optX,optY,"-r");

plt.show()

github.com

 

上記のシミュレーションを実行すると、

記事冒頭のような最適化結果のグラフが表示されます。

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

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

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

非線形最適制御入門 (システム制御工学シリーズ)

非線形最適制御入門 (システム制御工学シリーズ)