MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

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


共役勾配法について

共役勾配法は、有名な最適化アルゴリズムの一つで、

最急降下法の勾配ベクトルを逐次補正することで、

最急降下法よりも少ない計算回数で、

非線形最適化を実現できるアルゴリズムです。

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


冒頭の最適化のシミュレーション結果を見ると分かる通り、

少ないイタレーション数で、

最適値を求めることができます。


目的関数がある条件を満たす場合(線形関数など)、

ある一定回数の最適化で最適化可能であることが

数学的に証明されています。


アルゴリズムの導出に関しては、

末尾の関連資料を参考にしていただきたいですが、

イメージ的には、

過去の勾配ベクトルの値を考慮した上で、

現在の勾配ベクトルを計算することで、

最急降下法の問題である、

ジグザグに最適化を実施してしまうという問題を

解決できることが多くなります。


ロボティクスの分野では、

ニューラルネットやSLAMなどの最適化にしばしば利用されているようです。

Subgraph-preconditioned Conjugate Gradients for Large Scale SLAM

モデルが複雑で、ヘッセ行列までは計算できないが、

最急降下法よりも収束性を上げたい時に、

この共役勾配法が使用されています。


また、下記のwikiのページにある通り、

共役勾配法の勾配補正係数であるβの計算方法は複数あります。

今回のMATLABシミュレーションでは、

なんとなくPolak–Ribièreを使いました。

(特別この方法がよかったわけではありません)

 

共役勾配法の更新則

共役勾配法の更新則は下記の通りです。

MATLABサンプルプログラム

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

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

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

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

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

finalParam=ConjugateGradientMethodOptimization(param,x1range,x2range,yrange);
disp('Final Param is ');
finalParam

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

iterationMax=50;%イタレーションの最大回数
result=param;%パラメータの履歴

%共役勾配法用パラメータ初期化
d=-CalcJacobi(param(1),param(2));%勾配ベクトルの初期化
preJ=CalcJacobi(param(1),param(2));

for i=1:iterationMax
    
    %ヤコビ行列の取得
    J=CalcJacobi(param(1),param(2));
    
    %勾配のベクトルの更新
    beta=max(0,((J-preJ)*J')/(preJ*preJ'));
    d=-J+beta*d;
    
    %学習率の線形探索
    alpha=GoldenSection(0,0.5,param,d);
    
    % 共役勾配の計算
    param(1:2)=param(1:2)+alpha*d;
    param(3)=f(param(1),param(2));%評価関数の計算
    preJ=J;%一つ前のヤコビ行列をストア
    
    %シミュレーション結果の描画
    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);
    
    %収束判定
    if sum(abs(alpha*J))<=0.01
        disp(['Converge:',i]);
        break;
    end
end

finalParam=param(end,:);

function [minX]=GoldenSection(a,b,param,d)
%黄金分割法による線形探索関数

%黄金率
GOLDEN_RATIO = 1.6180339887498948482045868343656;
 
%内分点の計算
x1 = (a-b)/(GOLDEN_RATIO + 1.0) + b;
x2 = (a-b)/GOLDEN_RATIO + b;
%評価関数を両方計算するのは最初だけ
paramTmp=param(1:2)+x1*d;
f1=f(paramTmp(1),paramTmp(2));
paramTmp=param(1:2)+x2*d;
f2=f(paramTmp(1),paramTmp(2));
while 1
    %ループを回して両点を更新
    if f1 < f2
        a = x2;
        x2 = x1;
        f2 = f1;
        x1 = (a - b)/(GOLDEN_RATIO + 1.0) + b;
        paramTmp=param(1:2)+x1*d;
        f1=f(paramTmp(1),paramTmp(2));
    else
        b = x1;
        x1 = x2;
        f1 = f2;
        x2 = (a - b)/GOLDEN_RATIO + b;
        paramTmp=param(1:2)+x2*d;
        f2=f(paramTmp(1),paramTmp(2));
    end
    %収束判定
    if abs(a-b)<=10^-3
        minX=(a+b)/2;
        break
    end
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/ConjugateGradientMethod.m at master · AtsushiSakai/MATLABRobotics


また、今回のシミュレーションでは、

学習率の決定に、黄金分割法による線形探索を用いました。

myenigma.hatenablog.com

 

Pythonサンプルコード

f:id:meison_amsl:20160502210029p:plain

#!/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]
    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=np.array([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 ConjugateGradientMethod(start,Jacob):
    u"""
    Conjugate Gradient Method Optimization
    """

    result=start
    x=start
    preJ=None

    while 1:
        J=Jacob(x)

        #convergence check
        sumJ=sum([abs(alpha*j) for j in J])
        if sumJ<=0.01:
            print("OK")
            break

        if preJ is not None:
            beta=np.linalg.norm(J)**2/np.linalg.norm(preJ)**2
            grad=-1.0*J+beta*grad

        else:
            grad=-1.0*J
            
        x=x+[alpha*g for g in grad]
        result=np.vstack((result,x))
        #  print(x)

        if math.isnan(x[0]):
            print("nan")
            break
        

        preJ=-1.0*J


    return result

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

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

こちらのサンプルコードは、

黄金分割法を利用していません。

https://github.com/AtsushiSakai/PythonRobotics/blob/master/scripts/optimization/ConjugateGradientMethod/ConjugateGradientMethod.pygithub.com