数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【最適化】最急降下法(アルミホ条件による直線探索)C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法(反復法)である、最急降下法(アルミホ条件による直線探索)のC++コードを公開します。例題として

 \min (x-2)^2

を考えます。もちろん、最適解は  x^*=2 です。反復法とは、適当な初期値  x^{0} を定め、

 x^{n+1} = x^{n} + \alpha d

という漸化式によって値を更新していき、最終的に最適解  x^* へと収束させるような方法のことです。ここで、 d が探索方向、 \alpha がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。今回は、探索方向としては最急降下法を、ステップのサイズとしてはアルミホ(Armijo)条件による直線探索を行います。

参考文献

参考文献は『最適化と変分法』pp.35-43、『最適化とその応用』pp.130-139、それと以下のPDFファイルです。どちらの本も非常に参考になります。いつも読んでいます。
http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

基礎系 数学 最適化と変分法 (東京大学工学教程)

基礎系 数学 最適化と変分法 (東京大学工学教程)

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

最急降下法とは?

最急降下法とは、勾配が最もきつい方向、すなわちgradのマイナス方向に進んでいく手法のことを言います。勾配のきつい方向に進んでいくと底の部分にたどり着けるというわけです。ただし、場合によってはこの方法は効率がよくありません。そのような場合にも適用できる手法として共役勾配法などがあります。詳細は上記の参考文献を見てください。

アルミホ条件による直線探索とは?

アルミホ(Armijo)条件による直線探索とは、スタート地点での傾きに1より小さい係数をかけた傾きを使って直線をひき、その直線よりも関数の値が下回るようなステップサイズを採用する、というものです。上記のPDFの図がわかりやすいです。詳細は上記の参考文献を見てください。

C++コード

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//function//
inline double f(double x)
{
	return (x-2.0)*(x-2.0);
}

//gradient//
inline double fx(double x)
{
	return 2.0*(x-2.0);
}

//steepest gradient//
inline double SteepGrad(double x)
{
	return -fx(x);
}

//Armijo line search//
inline double Armijo(double x, double d, double alpha, double c1, double rho)
{
	alpha=1.0;
	for(int i=0;i<10000;i++)
	{
		if(f(x+d*alpha)<=f(x)+c1*fx(x)*d*alpha) break;
		else alpha=rho*alpha;
	}
	
	return alpha;
}

int main()
{
	double err=pow(10.0,-10.0);//tolerance
	double x=10.0;//initial value
	double alpha;//step size
	double rho=0.4;//contraction ratio
	double c1=0.001;//coefficient
	double d;//direction
	
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x))<<endl;
		
		//convergence check//
		if(abs(fx(x))<err) break;
		
		//step size search//
		d=SteepGrad(x);//direction
		alpha=Armijo(x,d,alpha,c1,rho);//step size;
		
		//update//
		x=x+alpha*d;
	}
	
	cout<<endl;
	cout<<x<<endl;
	
    return 0;
}