optimize_without_constraint.c

求解无约束最优化问题的各种算法集合!

//无约束最优化方法函数库
//Date: 	2004-11-24 
//Author:	lubangjian
//email:[e-mail]
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>

#define SUCCESS	1
#define FAIL	0
#define UnLimit	0xffffffff//用0xffffffff表示无穷大
#define SIGMA	0.5
#define U		0.1
#define false	0
#define true	1
float FirstStep=0.001;
//#define Enter(position)	printf("enter:%s\n",position)
long muticount=0,addcount=0;


#define Enter(position)

/*分配内存*/
void *_malloc(char *position,int memsize)
{
	void *mem;
	if ((mem=malloc(memsize))==NULL)
	{
		printf("%s:内存分配失败\n",position);
		return NULL;
	}
	return mem;
}

//输出最优结果
int output(double *para,
		   int paranum,
		   double (* get_funcval)(double *para,
							  int paranum))
{
	int paracount;

	printf("(");
	for (paracount=0;paracount<paranum;paracount++)
	{
		printf("%10.12f",para[paracount]);
		if (paracount<paranum-1)
			printf(",");
	}
	printf(")\t函数值:%10.12f\n",get_funcval(para,paranum));
	return SUCCESS;
}

int output_matrix(double *matrix,int row,int col)
{
	int rowcount,colcount;
	for (rowcount=0;rowcount<row;rowcount++)
	{
		for (colcount=0;colcount<col;colcount++)
			printf("%f ",matrix[rowcount*col+colcount]);
		printf("\n");
	}
	return 1;
}

//矩阵相乘
int matrix_multiply(double *matrix1,int row1,int col1,
					double *matrix2,int row2,int col2,
					double *res)
{
	int row1count,col1count,row2count,col2count;
	double tempval;

	if (col1!=row2)//检查是否可以相乘
		return FAIL;

	for (row1count=0;row1count<row1;row1count++)
	{
		for (col2count=0;col2count<col2;col2count++)
		{
			tempval=0.0;
			for (row2count=0;row2count<row2;row2count++)
			{
				col1count=row2count;
				tempval+=matrix1[row1count*col1+col1count]*matrix2[row2count*col2+col2count];
				muticount++;
				addcount++;
			}
			res[row1count*col2+col2count]=tempval;
		}
	}
	return SUCCESS;
}

#define ADD	0
#define SUB	1

//矩阵加减运算
int matrix_add_sub(double *matrix1,int row1,int col1,
				   double *matrix2,int row2,int col2,
				   double *matrix,
				   int flag)
{
	int row1count,col1count;

	if ((row1!=row2)||(col1!=col2))//检查是否可以进行运算
		return FAIL;

	for (row1count=0;row1count<row1;row1count++)
	{
		for (col1count=0;col1count<col1;col1count++)
		{
			if (flag==SUB)
				matrix[row1count*col1+col1count]=matrix1[row1count*col1+col1count]-matrix2[row1count*col1+col1count];
			else
				matrix[row1count*col1+col1count]=matrix1[row1count*col1+col1count]+matrix2[row1count*col1+col1count];
			addcount++;
		}
	}
	return SUCCESS;
}

#define MUL	0
#define DIV	1

//矩阵乘除一个数的运算
int matrix_mul_div(double *matrix1,int row1,int col1,
				   double *matrix,double num,
				   int flag)
{
	int row1count,col1count;

	if ((num==0.0)&&(flag==DIV))//检查是否可以相除
		return FAIL;

	for (row1count=0;row1count<row1;row1count++)
	{
		for (col1count=0;col1count<col1;col1count++)
		{
			if (flag==MUL)
				matrix[row1count*col1+col1count]=matrix1[row1count*col1+col1count]*num;
			else
				matrix[row1count*col1+col1count]=matrix1[row1count*col1+col1count]/num;
			muticount++;
		}
	}
	return SUCCESS;
}

//采用部分主元法的高斯消去法求方阵A的逆矩阵B
//A 方阵  (IN)
//n 方阵的阶数 (IN)
//B 方阵  (OUT)
//返回true 说明正确计算出逆矩阵
//返回false说明矩阵A不存在逆矩阵
int reverse_matrix(double *A,double *B,int row,int col)
{
	int i,j,k;
	double lMax,temp;
	//临时矩阵存放A
	double *tempm;

	if (!(tempm=_malloc("gaussj",sizeof(double)*row*col)))
		return false;
	memcpy(tempm,A,sizeof(double)*row*col);

	//初始化B为单位阵
	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
		{
			if(i!=j)
				B[i*col+j] = 0;
			else
				B[i*col+j] = 1;
		}
	}
	for(i=0;i<row;i++)
	{
		//寻找主元
		lMax = tempm[i*col+i];
		k = i;
		for(j=i+1;j<row;j++) //扫描从i+1到n的各行
		{
			if( fabs(tempm[j*col+i]) > fabs(lMax))
			{
				lMax = tempm[j*col+i];
				k = j;
			}
		}
		//如果主元所在行不是第i行，进行行交换
		if(k!=i)
		{
			for(j=0;j<col;j++)
			{
				temp = tempm[i*col+j] ;
				tempm[i*col+j] = tempm[k*col+j];
				tempm[k*col+j] = temp;
				//B伴随计算
				temp = B[i*col+j];
				B[i*col+j] = B[k*col+j];
				B[k*col+j] = temp;
			}
		}
		//判断主元是否是0，如果是，则矩阵A不是满秩矩阵，不存在逆矩阵
		if(tempm[i*col+i] == 0)
			return false;
		//消去A的第i列除去i行以外的各行元素
		temp = tempm[i*col+i];
		for(j=0;j<col;j++)
		{
			tempm[i*col+j] = tempm[i*col+j] / temp; //主对角线上元素变成1
			B[i*col+j] = B[i*col+j] / temp; //伴随计算
			muticount+=2;
		}

		for(j=0;j<row;j++) //行 0 -> n
		{
			if(j!=i)  //不是第i行
			{
				temp = tempm[j*col+i];
				for(k=0;k<col;k++) // j行元素 - i行元素* j行i列元素
				{
					tempm[j*col+k] = tempm[j*col+k] - tempm[i*col+k] * temp;
					B[j*col+k] = B[j*col+k] - B[i*col+k] * temp;
					muticount+=2;
					addcount+=2;
				}
			}
		}
	}
	return true;
}

//功能说明:计算两个向量的内积
//入参说明:vector1--向量1
//		   vector2--向量2
//		   elemnum--向量元素个数
//返 回 值:向量的模
double vector_inter_multiply(double *vector1,double *vector2,int elemnum)
{
	double res=0;
	int elemcount;

	Enter("vector_inter_multiply");
	for (elemcount=0;elemcount<elemnum;elemcount++)
	{
		res+=vector1[elemcount]*vector2[elemcount];
		muticount++;
		addcount++;
	}

	return res;
}

//功能说明:计算向量的模
//入参说明:vector--向量
//		   elemnum--向量元素个数
//返 回 值:向量的模
#define vector_val(vector,elemnum) sqrt(vector_inter_multiply(vector,vector,elemnum))


#define free_a_allocated(mem)	if (mem) free(mem)
#define free_allocated	free_a_allocated(para2);\
						free_a_allocated(gk1);\
						free_a_allocated(gk2)


//采用进退法进行精确一维搜索
#define ReturnRes		free(para1);\
		free(para2);\
		free(para0);\
		return SUCCESS;
//和wolfe统一,其中只有几个参数可用para,paranum,func
int forwardandback(double *ak,
		  double u,
		  double sigma,
		  double *pk,
		  double *para,
		  int (* get_gk)(double *para,int paranum,double *gk),
		  double (* func)(double *para,int paranum),
		  int paranum)
{
	double alpha=FirstStep,*para0,*para1,*para2,res0,res1,res2;//初始步长
	int paracount;

	if (!(para2=_malloc("forwardandback",sizeof(double)*paranum)))
		return FAIL;
	if (!(para1=_malloc("forwardandback",sizeof(double)*paranum)))
		return FAIL;
	if (!(para0=_malloc("forwardandback",sizeof(double)*paranum)))
		return FAIL;
	memcpy(para0,para,sizeof(double)*paranum);
	res0=func(para0,paranum);
	for (paracount=0;paracount<paranum;paracount++)
		para1[paracount]=para0[paracount]+alpha*pk[paracount];
	res1=func(para1,paranum);
	if (res0>=res1)
	{
		while(1)
		{
			alpha*=2;
			for (paracount=0;paracount<paranum;paracount++)
				para2[paracount]=para1[paracount]+alpha*pk[paracount];
			res2=func(para2,paranum);
			if (res1<=res2)
			{
				*ak=alpha/2;
				ReturnRes;
			}
			else
			{
/*
				alpha/=2;
				memcpy(para1,para2,sizeof(double)*paranum);
				for (paracount=0;paracount<paranum;paracount++)
					para1[paracount]=para0[paracount]+alpha*pk[paracount];
				res1=func(para1,paranum);
				if (res0>res1)
				{
					*ak=alpha*2;
					ReturnRes;
				}
				*/
				memcpy(para0,para1,sizeof(double)*paranum);
				memcpy(para1,para2,sizeof(double)*paranum);
				res0=res1;
				res1=res2;
			}
		}
	}
	else
	{
		while(1)
		{
			alpha*=-2;
			for (paracount=0;paracount<paranum;paracount++)
				para2[paracount]=para0[paracount]+alpha*pk[paracount];
			res2=func(para2,paranum);
			if (res0<=res2)
			{
				*ak=-alpha/2;
				ReturnRes;
			}
			else
			{
				memcpy(para1,para0,sizeof(double)*paranum);
				memcpy(para0,para2,sizeof(double)*paranum);
				res1=res0;
				res0=res2;
			}
		}
	}
}
		
//用wolfe求迭代步长
int wolfe(double *ak,
		  double u,
		  double sigma,