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

MyEnigma

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

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

Robot MATLAB

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

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

最急降下法について

最急降下法は、関数の微分情報を使った非線形最適化の手法の中で

最もシンプルなものの一つです。


非線形モデルの偏微分行列であるヤコビ行列を使って、

勾配方向にパラメータを更新することで、

最適化を実施します。


(ヤコビ行列に関しては、下記の記事を参考してもらえると良いと思います。

ロボティクスにおける線形代数 - MY ENIGMA)


パラメータの更新式は下記の通りです。

 

ちなみにαは学習パラメータと呼ばれ、

良く使用される方法としては、

事前に決めた一意のパラメータを利用するか、

黄金分割法のような直線探索法を使って最適なαを逐次決定します。

myenigma.hatenablog.com

最急降下法の注意点

通常の場合、

最急降下法で最適化解を得ることができますが、

下記のグラフのように、

ある方向では極小値だが、

ある方向では極大値である

いわゆる鞍点(あんてん)の場合、

この点に収束してしまう可能性があります。

f:id:meison_amsl:20160502124048p:plain

 

この鞍点に収束しているかどうかを

確認する場合は、

この点におけるヘッセ行列が正定か準正定かどうかを、

確認すればOKです。

 

ちなみに、正定かどうかの確認は、

ヘッセ行列Hの固有値がすべて正か、

または、下記の二次形式の式がすべてのxにおいて、

正であることが条件になります。

 

MATLABサンプルプログラム

% -------------------------------------------------------------------------
%
% File : SteepestDescentMethod.m
%
% Discription : Non-Linear optimizion with Steepest Descent Method
% 
% Environment : Matlab
%
% Author : Atsushi Sakai
%
% Copyright (c): 2014 Atsushi Sakai
%
% License : Modified BSD Software License Agreement
% -------------------------------------------------------------------------

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

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

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

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

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

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

iterationMax=50;%イタレーションの最大回数
alpha=0.01;%学習率
result=[];%パラメータの履歴
for i=1:iterationMax
    
    %ヤコビ行列の取得
    J=CalcJacobi(param(1),param(2));
    
    %最急降下法によるパラメータの更新
    param(1:2)=param(1:2)-alpha*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];


MATLABのサンプルプログラムは

下記のGitHubリポジトリで公開されています。

MATLABRobotics/SteepestDescentMethod.m at master · AtsushiSakai/MATLABRobotics

 

Pythonサンプルプログラム

f:id:meison_amsl:20160502160838p:plain

下記が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 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 SteepestDescentMethod(start,Jacob):
    u"""
    Steepest Descent Method Optimization
    """

    result=start
    x=start

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

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

    return result

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

result=SteepestDescentMethod(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()

下記のリポジトリでも公開されています。

github.com


シミュレーションで使用したコスト関数

今回のシミュレーションでは、Himmelblau関数を使いました。

Himmelblau's function - Wikipedia, the free encyclopedia


Himmelblau関数のヤコビ行列は下記の通りになります。


参考資料

myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com

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

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

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

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