【最適化】準ニュートン法（BFGS公式とアルミホ条件による直線探索）C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法（反復法）である、準ニュートン法（BFGS公式とアルミホ条件による直線探索）のC++コードを公開します。例題として2変数関数

$\min f(x_1,x_2) = \frac{1}{2} (x_1-x_2^2)^2 + \frac{1}{2} (x_2-2)^2$

を考えます。最適解は $\boldsymbol{x}^*= (4, 2)$ です。反復法とは、適当な初期値 $\boldsymbol{x}^{0}$ を定め、

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

という漸化式によって値を更新していき、最終的に最適解 $\boldsymbol{x}^*$ へと収束させるような方法のことです。ここで、 $\boldsymbol{d}$ が探索方向、 $\alpha$ がステップのサイズです。つまり、反復法を構成しているのは探索方向とステップのサイズです。準ニュートン法の場合、ニュートン法とは異なりステップのサイズは直線探索を行います。

参考文献

参考文献は『最適化とその応用』pp.156-165、『最適化と変分法』pp.49-57、『最適化法』pp.119-122それと以下のPDFファイルです。今回は、特に『最適化とその応用』のp.158の例6が非常に助けになりました。準ニュートン法の1ステップを具体的な数値を出しながら説明してくれています。これを丁寧に追えば、自分のコードが正しいか判断できます。『最適化とその応用』は他の話題でも具体的な数値例が出てくるので大好きです。

http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

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

作者: 矢部博
出版社/メーカー: 数理工学社
発売日: 2006/04
メディア: 単行本
購入: 1人クリック: 10回
この商品を含むブログ (11件) を見る

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

作者: 寒野善博,土谷隆,東京大学工学教程編纂委員会
出版社/メーカー: 丸善出版
発売日: 2014/10/22
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

最適化法 (工系数学講座 17)

作者: 田村明久,村松正和
出版社/メーカー: 共立出版
発売日: 2002/04/01
メディア: 単行本
購入: 2人クリック: 27回
この商品を含むブログ (8件) を見る

準ニュートン法とは？

準ニュートン法とは、ニュートン法を改良した手法です。ニュートン法の説明は以下の記事を見て下さい。

ニュートン法はまとめると

$\nabla^2 f(\boldsymbol{x}) \boldsymbol{d} = -\nabla f(\boldsymbol{x})$

という連立一次方程式を解いて $\boldsymbol{d}$ を求め

$\boldsymbol{x}^{n+1} = \boldsymbol{x}^{n} + \boldsymbol{d}$

のように値を更新して収束させる手法です。ただし、ニュートン法はヘシアン $\nabla^2 f(\boldsymbol{x})$ が正定値であれば降下方向（目的関数の値が減少する方向）を与えるのですが、一般にはヘシアンは常には正定値になりません。そこで、常に正定値となるようなヘシアンを近似する行列 $B$ をヘシアンの代わりに使おう、というのが準ニュートン法です。つまり、方向は

$B \boldsymbol{d} = -\nabla f(\boldsymbol{x})$

を解いて決めます。今回は連立一次方程式の解法として直接法であるガウスの消去法を使います。以下にC++コードがあります。

方向の次はステップのサイズですが、これはアルミホの条件で決めます。以下の記事を見て下さい。

さて最後に残ったのがヘシアンの近似行列 $B$ です。これは、BFGS公式と呼ばれる式を用いて更新していきます。Broyden、Fletcher、Goldfarb、Shannoらによって独立に提案されました。以下のような式です。

$B^{n+1} = B^{n} -\frac{B^{n}\boldsymbol{s}^{n} (B^{n}\boldsymbol{s}^{n})^{T}}{(\boldsymbol{s}^{n})^{T} B^{n} \boldsymbol{s}^{n}} + \frac{\boldsymbol{y}^{n} (\boldsymbol{y}^{n})^{T}}{(\boldsymbol{s}^{n})^{T} \boldsymbol{y}^{n}}$

ここで、 $\boldsymbol{s}^{n} = \boldsymbol{x}^{n+1} - \boldsymbol{x}^{n}$ 、 $\boldsymbol{y}^{n} = \nabla f(\boldsymbol{x}^{n+1}) - \nabla f(\boldsymbol{x}^{n})$ です。大事なのは、この $B$ がセカント条件と正定値条件を満たすことです。これらは、ヘシアンが持っていて然るべき性質をあらわしています。詳細は『最適化法』pp.120-121を見て下さい。

準ニュートン法の詳細は上記の参考文献を見てください。

C++コード

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

using namespace std;

//function//
inline double f(double x, double y)
{
	return 0.5*(x-y*y)*(x-y*y)+0.5*(y-2.0)*(y-2.0);
}

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

//gradient y//
inline double fy(double x, double y)
{
	return -2.0*y*(x-y*y)+y-2.0;
}

//matix vector multiplication Ax=b
inline void MatVec(double A[],double x[],double b[], int size)
{
	for(int i=0;i<size;i++)
	{
		b[i]=0.0;
		for(int j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	double value=0.0;
	for(int i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//Newton method, Gaussian elimination//
inline double Newton(int n, double b[], double B[], double x[])
{	
	double A[n*n];
	
	//copy//
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			A[i*n+j]=B[i*n+j];
		}
	}
	
	//forward elimination//
	for(int i=0;i<n-1;i++)
	{
		for(int j=i+1;j<n;j++)
		{
			int m=A[j*n+i]/A[i*n+i];
			
			for(int k=i+1;k<n;k++)
			{
				A[j*n+k]=A[j*n+k]-m*A[i*n+k];
			}
			
			b[j]=b[j]-m*b[i];
			
			A[j*n+i]=0.;
		}
	}
	
	//backward substitution//
	x[n-1]=b[n-1]/A[(n-1)*n+n-1];
	
	for(int i=n-2;i>=0;i--)
	{
		for(int j=i+1;j<n;j++)
		{
			b[i]-=A[i*n+j]*x[j];
		}
		
		x[i]=b[i]/A[i*n+i];
	}
}

//dysdic, vector*vector=matrix, a*b=ab//
inline void dyadic(double a[], double b[], double ab[], int size)
{
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			ab[i*size+j]=a[i]*b[j];
		}
	}
}

inline void Bupdate(double Bs[], double dgrad[], double B[], double sBs, double sdgrad, int size)
{
	double BsBs[size*size];
	double dgraddgrad[size*size];
	dyadic(Bs,Bs,BsBs,size);
	dyadic(dgrad,dgrad,dgraddgrad,size);
	
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			B[i*size+j]=B[i*size+j]-BsBs[i*size+j]/sBs+dgraddgrad[i*size+j]/sdgrad;
		}
	}
}

//Armijo line search//
inline double Armijo(double x, double y, double d[])
{
	double alpha=1.0;//step size
	double rho=0.4;//contraction ratio
	double c1=0.001;//coefficient
	
	for(int i=0;i<10000;i++)
	{
		if(f(x+d[0]*alpha,y+d[1]*alpha)<=f(x,y)+c1*(fx(x,y)*d[0]+fy(x,y)*d[1])*alpha) break;
		else alpha=rho*alpha;
	}
	
	return alpha;
}

int main()
{
	int n=2;
	double err=pow(10.0,-10.0);//tolerance
	double x,y;//unknown values
	double xnew,ynew;//next time step values
	double d[n];//direction and step size
	double B[n*n];//approximate hessian
	double grad[n];//gradient
	double s[n];
	double dgrad[n];
	double alpha;//step size
	double Bs[n];
	double sdgrad;
	double sBs;
	
	//initial value//
	x=0.0;
	y=0.0;
	B[0*n+0]=1.0;
	B[0*n+1]=0.0;
	B[1*n+0]=0.0;
	B[1*n+1]=1.0;
	
	//iteration//
	for(int i=0;i<10000;i++)
	{
		cout<<i<<" "<<abs(fx(x,y))<<" "<<abs(fy(x,y))<<endl;
		
		//convergence check//
		if((abs(fx(x,y))<err)&&(abs(fy(x,y))<err)) break;
		
		grad[0]=-fx(x,y);
		grad[1]=-fy(x,y);
		
		//direction//
		Newton(n,grad,B,d);
		
		//step size//
		alpha=Armijo(x,y,d);
		
		//update for x//
		xnew=x+alpha*d[0];
		ynew=y+alpha*d[1];
		
		//update for B//
		s[0]=alpha*d[0];
		s[1]=alpha*d[1];
		dgrad[0]=fx(xnew,ynew)-fx(x,y);
		dgrad[1]=fy(xnew,ynew)-fy(x,y);
		MatVec(B,s,Bs,n);
		sdgrad=ip(s,dgrad,n);
		sBs=ip(s,Bs,n);
		
		Bupdate(Bs,dgrad,B,sBs,sdgrad,n);
		
		//system("pause");
		//update//
		x=xnew;
		y=ynew;
	}
	
	cout<<endl;
	cout<<x<<" "<<y<<endl;
	
    return 0;
}