共役勾配法について
共役勾配法は、有名な最適化アルゴリズムの一つで、
最急降下法の勾配ベクトルを逐次補正することで、
最急降下法よりも少ない計算回数で、
非線形最適化を実現できるアルゴリズムです。
最急降下法による非線形最適化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
また、今回のシミュレーションでは、
学習率の決定に、黄金分割法による線形探索を用いました。
Pythonサンプルコード
#!/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
関連資料
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
myenigma.hatenablog.com
# MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。