ニュートン法による非線形最適化MATLAB＆Pythonサンプルプログラム

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

ニュートン法について

ニュートン法は、最適化アルゴリズムの代表的な方法の一つです。

以前紹介した、最急降下法や共役勾配法よりも収束性が高く、

最適化を試みる際には、一番最初に候補に上がる方法です。

共役勾配法による非線形最適化MATLABサンプルプログラム - MY ENIGMA

最急降下法による非線形最適化MATLABサンプルプログラム - MY ENIGMA

ニュートン法は、

最急降下法のようにモデルの一次微分情報(ヤコビ行列)だけでなく、

二次微分情報(ヘッセ行列)も利用して最適化を行うため、

高速な最適化が実現できます。

Newton法

更新式

しかし、モデルや評価関数が複雑な時には、

二次微分を求めることができないことも多いため、

その際には最急降下法や共役勾配法、

または、二次微分情報であるヘッセ行列を逐次推定する

準ニュートン法などが利用されるようです。

準ニュートン法 - Wikipedia

MATLABサンプルプログラム

% -------------------------------------------------------------------------
%
% File : NewtonMethod.m
%
% Discription : Newton Method Method
% 
% Environment : Matlab
%
% Author : Atsushi Sakai
%
% Copyright (c): 2014 Atsushi Sakai
%
% -------------------------------------------------------------------------

function [] = NewtonMethod()
close all;
clear all;

disp('Optimization simulation with Newton method starts!!')
%シミュレーション空間の最大最小範囲
maxxy=5;
minxy=-5;

%メッシュデータの生成
[x1range,x2range,yrange]=CreateMeshData(minxy,maxxy);

%初期パラメータの生成[x1 x2 y]
param=InitSearchParameter(minxy,maxxy);

finalParam=NewtonMethodOptimization(param,x1range,x2range,yrange);

function finalParam=NewtonMethodOptimization(param,x1range,x2range,yrange)
%最急降下法による最適化

iterationMax=50;%イタレーションの最大回数
alpha=1;%学習率
result=[];%パラメータの履歴
for i=1:iterationMax
    
    %ヤコビ行列の取得
    J=CalcJacobi(param(1),param(2));
    %ヘッセ行列の取得
    H=CalcHessian(param(1),param(2));
    
    %ニュートン法によるパラメータの更新
    param(1:2)=param(1:2)-(alpha*inv(H)*J')';
    param(3)=f(param(1),param(2));
    
    if sum(abs(alpha*J))<=0.01
        disp('Converge');
        i
        break;
    end
    
    %シミュレーション結果の描画
    figure(1)
    hold off;
    contour(x1range,x2range,yrange,100); hold on;
    result=[result;param];
    plot3(result(:,1),result(:,2),result(:,3),'.-k','markersize',10);
    xlabel('x1');
    ylabel('x2');
    view(2)
    drawnow;
    pause(0.5);
     
end

finalParam=param(end,:);

function param=InitSearchParameter(minxy,maxxy)
%初期パラメータを作成する関数

x1=maxxy - (maxxy-minxy).*rand(1,1);
x2=maxxy - (maxxy-minxy).*rand(1,1);
y=f(x1,x2);
param=[x1' x2' y'];

function [x1range,x2range,yrange]=CreateMeshData(minxy,maxxy)
%シミュレーション空間のメッシュデータを作成する関数

%メッシュデータの作成
[x1range,x2range]=meshgrid(minxy:0.3:maxxy);
yrange=f(x1range,x2range);

function y = f(x1,x2)
% Himmelblau's function
% see Himmelblau's function - Wikipedia, the free encyclopedia 
% http://en.wikipedia.org/wiki/Himmelblau%27s_function
y=(x1.^2+x2-11).^2+(x1+x2.^2-7).^2;

function J= CalcJacobi(x,y)
% jacobi matrix of Himmelblau's function
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=[dx dy];

function H= CalcHessian(x,y)
% Hessian matrix of Himmelblau's function
dxx=12*x^2+4*y-42;
dxy=4*x+4*y;
dyy=4*x+12*y^2-26;
H=[dxx dxy;
   dxy dyy];

MATLABのサンプルプログラムは下記のGitHubリポジトリで公開しています。

MATLABRobotics/NewtonMethod.m at master · AtsushiSakai/MATLABRobotics

Pythonサンプルコード

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

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

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

def Hessian(state):
    u"""
    Hessian matrix of Himmelblau's function
    """
    x=state[0]
    y=state[1]
    dxx=12*x**2+4*y-42;
    dxy=4*x+4*y
    dyy=4*x+12*y**2-26
    H=np.array([[dxx,dxy],[dxy,dyy]])
    return H
    

def Jacob(state):
    u"""
    jacobi matrix of Himmelblau's function
    """
    x=state[0]
    y=state[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=[dx,dy]
    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 NewtonMethod(start,Jacob):
    u"""
    Newton Method Optimization
    """

    result=start
    x=start

    while 1:
        J=Jacob(x)
        H=Hessian(x)
        sumJ=sum([abs(alpha*j) for j in J])
        if sumJ<=0.01:
            print("OK")
            break

        grad=-np.linalg.inv(H).dot(J) 
        print(grad)

        x=x+[alpha*j for j in grad]
        
        result=np.vstack((result,x))

    return result

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

result=NewtonMethod(start,Jacob)
(X,Y,Z)=CreateMeshData()
CS = plt.contour(X, Y, Z,nContour)
#  plt.clabel(CS, inline=1, fontsize=10)
#  plt.title('Simplest default with labels')

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

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

plt.show()