mathyw.cpp

C++最优化方法源程序

		{
			a = b;
			diver_a = diver_b;
			step = 2* step;
		}
	}while( diver_b <= 1E-15 );
	for( i=1; i<=n; i++ )
		temp_point[i-1] = start[i-1]+ a* direction[i-1];
	value_a = (* pf)( temp_point );
	for( i=1; i<=n; i++ )
		temp_point[i-1] = start[i-1]+ b* direction[i-1];
	value_b = (* pf)( temp_point );
	do
	{
		s = 3* ( value_b - value_a )/ ( b - a );
		z = s - diver_a - diver_b;
		w = sqrt( z* z- diver_a* diver_b );
		t = a+ ( w- z- diver_a )* ( b- a )/ (diver_b- diver_a+ 2* w );
		value_b = (* pf)( temp_point );
		for( i=1; i<=n; i++ )
			temp_point[i-1] = start[i-1]+ t* direction[i-1];
		value_t = (* pf)( temp_point );
		compute_grad( pf, n, temp_point, grad );
		diver_t = 0;
		for( i=1; i<=n; i++ )
			diver_t = diver_t+ grad[i-1]* direction[i-1];
		if( diver_t > 1E-6 )
		{
			b = t;
			value_b = value_t;
			diver_b = diver_t;
		}
		else if( diver_t < -1E-6 )
		{
			a = t;
			value_a = value_t;
			diver_a = diver_t;
		}
		else break;
	}while( fabs( diver_t ) >= 1E-6 && ( fabs(b- a ) > 1E-6 ) );
	delete[] grad;
	delete[] temp_point;
	return(t);
}
/////////////////////////////////////////////////////////////////////////////////////
//void conjugate_search( double (* pf)( double * x), int n, double * start, double * point )
//         共轭梯度法求极小点的函数。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果，
//返回值：满足精度的极小点的一个近似值。
//函数调用：
//    梯度函数void compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
//    一维收索函数double line_search2( double (*pf)( double * x ), int n, double * start, double * direction )
/////////////////////////////////////////////////////////////////////////////////////
void conjugate_search( double (* pf)( double * x), int n, double * start, double * point )
{
	double * direction = new double[n];
	double * grad1 = new double[n];
	double * grad2 = new double[n];
	double length, value, data;
	double precision = 5E-3;
	
	compute_grad( pf, n, start, grad1 );
	do
	{
		for( int i=0; i<n; i++ )
			direction[i] = - grad1[i];
		length = line_search2( pf, n, start, direction );
		for( i=0; i<n; i++ )
			point[i] = start[i]+ length* direction[i];
	    for( int k=1; k<n; k++ )
		{
		    compute_grad( pf, n, point, grad2 );
			data = 0;
			value = 0;
			for( i=0; i<n; i++ )
			{
				value += grad2[i]* ( grad2[i]- grad1[i] );
				data += grad1[i]* grad1[i];
				data = sqrt( data );
			}
			value /= data;
			for( i=0; i<n; i++ )
			 	direction[i] = value* direction[i]- grad2[i];
			length = line_search2( pf, n, point, direction );
			for( i=0; i<n; i++ )
				point[i] += length* direction[i];
		}
		compute_grad( pf, n, point, grad1 );
		value = 0;
		for( i=0; i<n; i++ )
		{
			value += grad1[i]* grad1[i];
			value = sqrt( value );
		}
	}while( value > precision );
}
/////////////////////////////////////////////////////////////////////////////////////
//void newton_gill_search( double (* pf)( double * x), int n, double * start, double * point )
//         Gill-Murray改进牛顿法求极小点的函数。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果，
//函数调用：
//   梯度函数void compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
//   一维收索函数double line_search2( double (*pf)( double * x ), int n, double * start, double * direction )
//   强迫正定Cholesky分解函数void force_cholesky( double * matrix, int n, double * triangle_low, double * diagnoal )
//   计算函数的二阶导数距阵void hessian( double ( * pf)( double * x ), int n, double * point, double * H )
/////////////////////////////////////////////////////////////////////////////////////
void newton_gill_search( double (* pf)( double * x), int n, double * start, double * point )
{
	double precision = 1E-5;
	double * hesan = new double[n* n];
	double * grad = new double[n];
	double * diag = new double[n];
	double * low = new double[n* n];
	double * direction = new double[n];
	double * alternate = new double[n];
	double value, length;

	for( int i=0; i<n; i++ )
		point[i] = start[i];  //得到初始点
	do
	{
		compute_grad( pf, n, point, grad );
		hessian( pf, n, point, hesan );
		force_cholesky( hesan, n, low, diag );
		value = 0;
		for( i=0; i<n; i++ )
			value += grad[i]* grad[i];

		int j;
		if( sqrt( value ) >= precision )
		{
			//解线性方程组
			for( i=0; i<n; i++ )
			{
				alternate[i] = -grad[i];
				for( j=0; j<i; j++ )
				    alternate[i] -= low[i* n+ j]* alternate[j];
			}
			for( i=n-1; i>=0; i-- )
			{
				direction[i] = alternate[i]/ diag[i];
				for( int j=i+1; j<n; j++ )
				    direction[i] -= low[j* n+ i]* direction[j];
			}
			alternate[0] = -1;   //让程序继续收索下一点
		}
		else
		{
			for( i=0; i<n; i++ )
			{
				alternate[i] = diag[i];
				for( j=0; j<i; j++ )
				    alternate[i] += low[i* n+ j]* low[i* n+ j]* diag[j];
				alternate[i] -= hesan[i* n+ i];
				alternate[i] = diag[i]- alternate[i];
			}
			j=0;
			for( i=0; i<n; i++ )
				if( alternate[i] < alternate[0] )
				{
					alternate[0] = alternate[i];
					j = i;
				}
			if( alternate[0] < 0 )
			{
				memset( direction, 0, n* sizeof( double ) );
				direction[j] = 1;
				for( i=j-1; i>=0; i-- )
				{
					for( int k=j; k>i; k-- )
						direction[i] -= low[i* n+ k]* direction[k];
				}
				value = 0;
				for( i=0; i<=j; i++ )
					value += grad[i]* direction[i];
				if( value >= 0 )
				{
					for( i=0; i<=j; i++ ) 
						direction[i] = -direction[i];
				}	
			}
		}
		if( alternate[0] < 0 )
		{
			length = line_search2( pf, n, point, direction );
			for( i=0; i<n; i++ )
				point[i] += length* direction[i];
		}
	}while( alternate[0] < 0 );

	delete[] hesan;
	delete[] grad;
	delete[] diag;
	delete[] low;
	delete[] direction;
	delete[] alternate;
}
/////////////////////////////////////////////////////////////////////////////////////
//void newton_region_search( double (* pf)( double * x), int n, double * start, double * point )
//         信赖域法求极小点的函数。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果，
//返回值：满足精度的极小点的一个近似值。
//函数调用：
//   梯度函数void compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
//   Cholesky分解函数void cholesky( double * matrix, int n, double * triangle_low, double * diagnoal )
//   计算函数的二阶导数距阵void hessian( double ( * pf)( double * x ), int n, double * point, double * H )
/////////////////////////////////////////////////////////////////////////////////////
void newton_region_search( double (* pf)( double * x), int n, double * start, double * point )
{
	double precision = 1E-3;
	double * hesan = new double[n* n];
	double * grad = new double[n];
	double * diag = new double[n];
	double * low = new double[n* n];
	double * alternate = new double[n];
	double * shift = new double[n];

	double value, solution, length;

	for( int i=0; i<n; i++ )
		point[i] = start[i];  //得到初始点
	do
	{
		compute_grad( pf, n, point, grad );
		hessian( pf, n, point, hesan );
		value = 0;
		for( i=0; i<n; i++ )
			value += grad[i]* grad[i];		
		if( sqrt( value )- precision >= 1E-15 )
		{
			length = 0.1;
			do
			{
				for( i=0; i<n; i++ )
					hesan[i* n+ i] += length;
				length *= 4;
			}while( !cholesky( hesan, n, low, diag ) );
			do
			{		
				//解线性方程组
				for( i=0; i<n; i++ )
				{
					alternate[i] = -grad[i];
					for( int j=0; j<i; j++ )
					    alternate[i] -= low[i* n+ j]* alternate[j];
				}
				for( i=n-1; i>=0; i-- )
				{
					shift[i] = alternate[i]/ diag[i];
					for( int j=i+1; j<n; j++ )
					    shift[i] -= low[j* n+ i]* shift[j];
				}

				for( i=0; i<n; i++ )
					alternate[i] = point[i]+ shift[i];
				solution = 0;
				value = 0;
				for( i=0; i<n; i++ )
				{
					value = grad[i]* shift[i];
					for( int j=0; j<n; j++ )
						solution += shift[i]* hesan[i* n+ j]* shift[j]/ 2;
				}
				solution += value;
				value = pf( point )- pf( alternate );
				if( value< 0 )
					cholesky( hesan, n, low, diag );
			}while( value< 0 );

			for( i=0; i<n; i++ )
				point[i] += shift[i];
			if( value> 0.25* solution && value<= 0.75* solution )
				length /= 4;
			if( value> 0.75* solution )
				length /= 2;
			alternate[0] = -1;    //让程序继续收索下一点
		}
		else
		{
			for( i=0; i<n; i++ )
			{
				alternate[i] = diag[i];
				for( int j=0; j<i; j++ )
				    alternate[i] += low[i* n+ j]* low[i* n+ j]* diag[j];
				alternate[i] -= hesan[i* n+ i];
				alternate[i] = diag[i]- alternate[i];
			}
			int j=0;
			for( i=0; i<n; i++ )
				if( alternate[i] < alternate[0] )
				{
					alternate[0] = alternate[i];
					j = i;
				}
			if( alternate[0]< 0 )
			{
				memset( shift, 0, n* sizeof( double ) );
				shift[j] = 1;
				for( i=j-1; i>=0; i-- )
				{
					for( int k=j; k>i; k-- )
						shift[i] -= low[i* n+ k]* shift[k];
				}
				value = 0;
				for( i=0; i<=j; i++ )
					value += grad[i]* shift[i];
				if( value >= 0 )
				{
					for( i=0; i<=j; i++ ) 
						shift[i] = -shift[i];
				}
				for( i=0; i<n; i++ )
					alternate[i] = point[i]+ shift[i];
				while( pf( point )- pf( alternate ) <= 0 )
				{
					for( i=0; i<n; i++ )
						shift[i] /= 2;
					alternate[i] = point[i]+ shift[i];
				}
				for( i=0; i<n; i++ )
					point[i] = alternate[i];
			}
		}
	}while( alternate[0] < 0 );

	delete[] hesan;
	delete[] grad;
	delete[] diag;
	delete[] low;
	delete[] shift;
	delete[] alternate;
}
/////////////////////////////////////////////////////////////////////////////////////
//void DFP_search( double (* pf)( double * x), int n, double * start, double * point )
//         变尺度DFP法求极小点的函数。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果，
//函数调用：
//   梯度函数void compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
//   一维收索函数double line_search2( double (*pf)( double * x ), int n, double * start, double * direction )
/////////////////////////////////////////////////////////////////////////////////////
void DFP_search( double (* pf)( double * x), int n, double * start, double * point )
{
	double precision = 1E-5;
	double * f_hesan = new double[n* n];
	double * grad = new double[n];
	double * matrix = new double[n* n];
	double * direction = new double[n];
	double * dg = new double[n];
	double * dp = new double[n];
	double value, length, xandg, g_h_g;

	memset( f_hesan, 0, n* n* sizeof( double ) );
	for( int i=0; i<n; i++ )
	{
		point[i] = start[i];  //得到初始点
		f_hesan[i* n+ i] = 1;    //初始 H 为 I
	}
	value = 0;
	compute_grad( pf, n, point, grad );
	for( i=0; i<n; i++ )
	{
		value += grad[i]* grad[i];     //计算||grad||
	}
	if( sqrt( value ) >= precision )
	{
		for( i=0; i<n; i++ )
		{