mathyw.cpp

C++最优化方法源程序

//mathyw.cpp    ：最优化计算模块的实现文件，应用时需加到工程中去.
#include "iostream.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"

//////////////////////////////////////////////////////////////////////////////////////
//compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
//    计算多元函数在某一点的梯度值.
//参数值:
//    pf: 所求函数的指针.
//    x: 函数的自变量(用数组指针来表示).
//    n: 自变量的个数.
//    point: 所求的点(用数组指针来表示).
//    grad: 存放结果的数组指针.
//////////////////////////////////////////////////////////////////////////////////////
void compute_grad( double ( * pf)(double * x ), int n , double * point, double * grad )
{
	double h = 1E-3;
	int i;
	double * temp = new double[n];
	for( i=0; i<n; i++ )
		temp[i] = point[i];
	for( i=0; i<n; i++ )
	{
		temp[i] += 0.5*h;
		grad[i] = 4* pf( temp )/ ( 3* h );
		temp[i] -= h;
		grad[i] -= 4* pf( temp )/ ( 3* h );
		temp[i] += 3* h/ 2;
		grad[i] -= pf( temp )/ ( 6* h );
		temp[i] -= 2* h;
		grad[i] += pf( temp )/ ( 6* h );
		temp[i] = point[i];
	}
	delete []temp;
}
//////////////////////////////////////////////////////////////////////////////////////
//hessian( double ( * pf)( double * x ), int n, double * point, double * H )
//     计算函数的二阶导数距阵
//参数值:
//     pf: 所求函数的指针.
//     x: 自变量数组指针.
//     n: 自变量个数.
//     point: 所求的点(用数组指针来表示).
//     H: 存放结果的数组指针.
//////////////////////////////////////////////////////////////////////////////////////
void hessian( double ( * pf)( double * x ), int n, double * point, double * H )
{
	int i, j;
	double h = 1E-3;
	double * temp = new double[n];
	for( i=0; i<n; i++ )
		temp[i] = point[i];
	for( i=0; i<n; i++ ) 
	{
		temp[i] += h/2;
		H[i* n+ i] = 16* pf( temp );
		temp[i] -= h;
		H[i* n+ i] += 16* pf( temp );
		temp[i] += 3* h/ 2;
		H[i* n+ i] -= pf( temp );
		temp[i] -= 2* h;
		H[i* n+ i] -= pf( temp );
		H[i* n+ i] -= 30* pf( point );
		H[i* n+ i] /= 3* h* h;
		temp[i] = point[i];
	}
	for( i=0; i<n; i++ ) 
	    for( j = i+1; j<n; j++ ) 
		{
			temp[i] += h/ 2;
			temp[j] += h/ 2;
			H[i* n+ j] = 16* pf( temp );
			temp[i] -= h;
			temp[j] -= h;
			H[i* n+ j] += 16* pf( temp );
			temp[i] += 3* h/ 2;
			temp[j] += 3* h/ 2;
			H[i* n+ j] -= pf( temp );
			temp[i] -= 2* h;
			temp[j] -= 2* h;
			H[i* n+ j] -= pf( temp );
			H[i* n+ j] -= 30* pf( point );
			H[i* n+ j] /= ( 3* h* h );
			H[i* n+ j] -= H[i* n+ i];
			H[i* n+ j] -= H[j* n+ j];
			H[i* n+ j] /= 2;
			temp[i] = point[i];
			temp[j] = point[j];
		}
	for( i=0; i<n; i++ ) 
	    for( j = i+1; j<n; j++ ) 
			H[j* n+ i] = H[i* n+ j];
	delete []temp;
}
///////////////////////////////////////////////////////////////////////////////////////
//bool cholesky( double * old_matrix, int n, double * triangle_low, double * diagonal )
//      进行Cholesky分解的函数
//参数值:
//      old_matrix: 输入的对称矩阵，所得矩阵与原矩阵相差的对角真
//      n: 矩阵的维数，
//      triangle_low: 记录下三角矩阵结果，
//      diagonal: 记录对角矩阵结果。
//返回值：
//      能正定分解则返回TRUE，否则返回FALSE。
///////////////////////////////////////////////////////////////////////////////////////
bool cholesky( double * matrix, int n, double * triangle_low, double * diagonal )
{
	memset( triangle_low, 0, n* n* sizeof( double ) );
	for( int i=0; i<n; i++ )
		triangle_low[i* n+ i] = 1;
	double * c = new double[n* n];

	for( int j=0; j<n; j++ )
	{
		diagonal[j] = matrix[j* n+ j];
		for( int r=0; r<j; r++ )
		{
			int t = j* n+ r;
			diagonal[j] -= c[t]* c[t]/ diagonal[r];
		}
		if( diagonal[j] <= 0 )
		{
			cout<< " 这不是正定矩阵!"<< endl;
			return false;
		}
		for( int i=j+1 ; i<n; i++ )
		{
			int t = i* n+ j;
			c[t] = matrix[t];
			for( int r=0; r<j; r++ )
			{
			    c[t] -= c[j* n+ r]* c[i* n+ r]/ diagonal[r];
			}
			triangle_low[t] = c[t]/ diagonal[j];
		}
	} 
	delete[] c;
	return true;
}
///////////////////////////////////////////////////////////////////////////////////////
//force_cholesky( double * matrix, int n, double * triangle_low, double * diagnoal )
//      进行强迫正定Cholesky分解的函数
//参数值:
//      matrix: 输入的对称矩阵，
//      n: 矩阵的维数，
//      triangle_low: 记录下三角矩阵结果，
//      diagnoal: 记录对角矩阵结果，
///////////////////////////////////////////////////////////////////////////////////////
void force_cholesky( double * matrix, int n, double * triangle_low, double * diagnoal )
{
	double precision = 1e-15;
	memset( triangle_low, 0, n* n* sizeof( double ) );
	for( int i=0; i<n; i++ )
		triangle_low[i* n+ i] = 1;

	double * c = new double[n* n];
	double * diag = new double[n];
	double x = fabs( matrix[0] ), y = fabs( matrix[n] ), z;
	double alternate, k;

	for( i=0; i<n; i++ )
	{
		alternate = matrix[i* n+ i];
		if( alternate > x)
			x = alternate; 
		for( int j= i+1; j<n; j++ )
		{
			alternate = matrix[j* n+ i];
			if( alternate > y )
				y = alternate;
		}
	}
	z = ( x > y/n ) ? x : y/n;

	for( int j=0; j<n; j++ )
	{
		diag[j] = matrix[j* n+ j];
		for( int r=0; r<j; r++ )
		{
			int t = j* n+ r;
			diag[j] -= c[t]* c[t]/ diagnoal[r];
		}
		diag[j] = fabs( diag[j] );
		if( diag[j] < precision )
			diag[j] = precision;
		for( int i=j+1 ; i<n; i++ )
		{
			int t = i* n+ j;
			c[t] = matrix[t];
			for( int r=0; r<j; r++ )
			{
			    c[t] -= c[j* n+ r]* c[i* n+ r]/ diagnoal[r];
			}
		}
		k = fabs( c[(j+1)* n+ j] );
		for( i=j+1; i<n; i++ )
		{
			alternate = fabs( c[i* n+ j] );
			if( k < alternate )
				k = alternate;
		}
		k *= k/ z;
		diagnoal[j] = ( diag[j] > k ) ? diag[j] : k;
		for( i=j+1; i<n; i++ )
			triangle_low[i* n+ j] = c[i* n+ j] / diagnoal[j];
	} 

	delete[] c;
	delete[] diag;
}
//////////////////////////////////////////////////////////////////////////////////
//double line_search1( double (*pf)( double * x ), int n, double * start, double * direction )
//       直接法的一维收索模块。
//参数值：
//     pf: 所求函数的指针.
//     x: 自变量数组指针.
//     n: 自变量个数.
//     start: 收索的起始点(用数组指针来表示).
//     direction: 收索的方向.
//返回值：最优步长。
//////////////////////////////////////////////////////////////////////////////////
double line_search1( double (*pf)( double * x ), int n, double * start, double * direction )
{
	double * point_l = new double[n];
	double * point_r = new double[n];
	double * point3  = new double[n];
	double left, right, left_in, right_in, value_l, value_r;
	double preportion = ( sqrt(5)- 1 )/ 2;
	double value, length;
	double precision = 1E-14;

	//找到包含极小点的区间
	left = 0, length = 0.1;
	value_l = pf( start );
	bool flag;
	do
	{
		flag = false;
		right_in = left+ length;
		for( int i=0; i<n; i++ )
		{
			point_r[i] = start[i]+ right_in* direction[i];
		}
		value_r = pf( point_r );
		if( value_r- value_l >= precision )
		{
			length = -length;
			right_in = left +length;
			for( int i=0; i<n; i++ )
			{
				point_r[i] = start[i]+ right_in* direction[i];
			}
			value_r = pf( point_r );
			if( value_r- value_l >= precision )
			{
				length /= 2;
				if( fabs( length ) -precision < precision )  return 0;
				flag = true;
			}
		}
	}while( flag );
	do
	{
		flag = false;
		length *= 3;
		right = right_in +length;
		for( int i=0; i<n; i++ )
		{
			point3[i] = start[i]+ right* direction[i];
		}
		value = pf( point3 );
		if( value -value_r <= precision )
		{
			flag = true;
			left = right_in;
			right_in = right;
			value_l = value_r;
			value_r = value;
		}
	}while( flag );
//*
	if( length < 0 )
	{
		value = right;
		right = left;
		left  = value;
	}
//*/
	value = preportion* ( right- left );
	left_in = right- value;
	right_in = left+ value;
	
	//根据步长求出各对应点
	for( int i=0; i<n; i++ )
	{
		point_l[i] = start[i]+ left_in* direction[i];
		point_r[i] = start[i]+ right_in* direction[i];
	}
	value_l = pf( point_l );
	value_r = pf( point_r );

	//进行收索
	do
	{
		if( value_l - value_r < 1E-15 )
		{
			right = right_in;
			value = right- left;
			if( fabs( value -precision ) > precision )
			{
				right_in = left_in;
				value_r = value_l;
				left_in = right- preportion* value;
				for( i=0; i<n; i++ )
					point_l[i] = start[i]+ left_in* direction[i];
				value_l = pf( point_l );
			}
		}
		else 
		{
			left = left_in;
			value = right- left;
			if( fabs( value -precision ) > precision )
			{
				left_in = right_in;
				value_l = value_r;
				right_in = left+ preportion* value;
				for( i=0; i<n; i++ )
					point_r[i] = start[i]+ right_in* direction[i];
				value_r = pf( point_r );
			}
		}
	}while( fabs( value -precision ) > precision );
	return ( right+ left )/ 2;

	delete[] point3;
	delete[] point_r;
	delete[] point_l;
}

//////////////////////////////////////////////////////////////////////////////////
//double line_search2( double (*pf)( double * x ), int n, double * start, double * direction )
//       解析法的一维收索模块。
//参数值：
//     pf: 所求函数的指针.
//     x: 自变量数组指针.
//     n: 自变量个数.
//     start: 收索的起始点(用数组指针来表示).
//     direction: 收索的方向.
//返回值：最优步长。
//调用函数：
//     梯度计算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 )
{
	int i;
	double step = 0.001;
	double a = 0, value_a, diver_a;
	double b, value_b, diver_b;
	double t, value_t, diver_t;
	double s, z, w;
	double * grad, * temp_point;

	grad = new double[n];
	temp_point = new double[n];
	compute_grad( pf, n, start, grad );
	diver_a = 0;
	for( i=1; i<=n; i++ )
		diver_a = diver_a+ grad[i-1]* direction[i-1];
	do
	{
		b = a+ step;
		for( i=1; i<=n; i++ )
			temp_point[i-1] = start[i-1]+ b* direction[i-1];
		compute_grad( pf, n, temp_point, grad );
		diver_b = 0;
		for( i=1; i<=n; i++ )
			diver_b = diver_b+ grad[i-1]* direction[i-1];
		if( fabs(diver_b) < 1E-10 )
		{
			delete[] grad;
			delete[] temp_point;
			return(b);
		}
		if( diver_b < -1E-15 )