はじめに

以前、

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

myenigma.hatenablog.com

上記のニュートン法は、

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

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

複雑なモデルの場合、

各変数で二次微分した

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

数式が複雑で難しかったり、正則でない場合は、

逆行列を計算することができず、

そもそもヘッセ行列も計算できないという問題があります。

また、計算したヘッセ行列が正定値でない場合は、

そもそも勾配方向が正確に計算できない場合もあります。

そこで、

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

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

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

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

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

今回の記事では、

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

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

準ニュートン法

準ニュートン法は、

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

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

更新則としては、

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

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

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

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

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

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

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

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

(引用: 準ニュートン法 - Wikipedia)

いずれの方法も、

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

Δ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()