【最適化】多変数ニュートン法　C++コード付き

はじめに

今回は無制約最小化問題に対する数値解法（反復法）である、多変数ニュートン法のC++コードを公開します。例題として2変数関数

$\min f(x,y) = x^3 + y^3 -9xy +27$

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

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

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

参考文献

参考文献は『これならわかる最適化数学』pp.88-93、『最適化と変分法』pp.44-49、『最適化とその応用』pp.151-155、それと以下のPDFファイルです。『これならわかる最適化数学』は特にビギナー向けで読みやすいです。
http://www-optima.amp.i.kyoto-u.ac.jp/~nobuo/Ryukoku/2002/course7.pdf

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

作者: 金谷健一
出版社/メーカー: 共立出版
発売日: 2005/09/01
メディア: 単行本
購入: 29人クリック: 424回
この商品を含むブログ (42件) を見る

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

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

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

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

ニュートン法とは？

ニュートン法とは、目的関数（最小化しようとする対象）に対する1次の最適性条件

$\nabla f(\boldsymbol{x}) = \boldsymbol{0}$

に対して、非線型方程式に対するニュートン法を適用したもの、と言えます。詳しく説明しましょう。現在 $\boldsymbol{x}^{n}$ にいて $\boldsymbol{d}$ だけ進んだ時に1次の最適性条件を満たすとしましょう。つまり、

$\nabla f(\boldsymbol{x}+\boldsymbol{d}) = \boldsymbol{0}$

です。我々が知りたいのは $\boldsymbol{d}$ です。上の式を1次までテイラー展開すると

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

となります。 $\nabla^2 f(\boldsymbol{x})$ はヘシアンと呼ばれています。ここから $\boldsymbol{d}$ について解くと

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

となります。これが多変数ニュートン法における探索方向とステップサイズです。一度に求まります。多変数の場合は行列になるので逆行列をかける必要があります。これが1変数ニュートン法との違いです。この $\boldsymbol{d}$ を用いて

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

のように値を更新して収束させればよいのです。ただし、上記のように逆行列を求めるのは計算コストが高いので、実際には

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

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

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

C++コード

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

using namespace std;

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

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

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

//hessian xy//
inline double fxy(double x, double y)
{
	return -9.0;
}

//hessian fxx//
inline double fxx(double x, double y)
{
	return 6.0*x;
}

//hessian fyy//
inline double fyy(double x, double y)
{
	return 6.0*y;
}

//Newton method, Gaussian elimination//
inline double Newton(int n, double b[], double A[], double x[])
{	
	//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];
	}
}

int main()
{
	int n=2;
	double err=pow(10.0,-10.0);//tolerance
	double x,y;//unknown
	double d[n];//direction and step size
	double Hess[n*n];//hessian
	double grad[n];//gradient;
	
	//initial value//
	x=5.0;
	y=7.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);
		
		Hess[0*n+0]=fxx(x,y);
		Hess[0*n+1]=fxy(x,y);
		Hess[1*n+0]=fxy(x,y);
		Hess[1*n+1]=fyy(x,y);
		
		Newton(n,grad,Hess,d);
		
		//update//
		x=x+d[0];
		y=y+d[1];
	}
	
	cout<<endl;
	cout<<x<<" "<<y<<endl;
	
    return 0;
}