optimize_without_constraint.c

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

		  double *pk,
		  double *para,
		  int (* get_gk)(double *para,int paranum,double *gk),
		  double (* func)(double *para,int paranum),
		  int paranum)
{
	double alpha=1,a=0,b=UnLimit;
	double funcval1,funcval2,*para2,*gk1,*gk2;
	int paracount;
	char *position="wolfe";

	Enter(position);
	//printf("sigma=%f,u=%f\n",sigma,u);
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;

	get_gk(para,paranum,gk1);//计算原先的梯度
	funcval1=func(para,paranum);//计算第一个函数值
	while(1)
	{
		for (paracount=0;paracount<paranum;paracount++)
			para2[paracount]=para[paracount]+alpha*pk[paracount];
		get_gk(para2,paranum,gk2);//计算新的梯度
		funcval2=func(para2,paranum);//计算第二个函数值
		if ((funcval1-funcval2)>=(-u*alpha*vector_inter_multiply(gk1,pk,paranum)))//满足条件1.6
		{
			if (vector_inter_multiply(gk2,pk,paranum)>=sigma*vector_inter_multiply(gk1,pk,paranum))//满足条件1.7
			{
				free_allocated;
				*ak=alpha;
				return SUCCESS;
			}
			//不满足1.7
			a=alpha;
			if (2*alpha>(alpha+b)/2)
				alpha=(alpha+b)/2;
			else
				alpha=2*alpha;
			muticount++;
			addcount++;
			continue;
		}
		//不满足1.6
		b=alpha;
		alpha=(a+alpha)/2;
		muticount++;
		addcount++;
	}
	return FAIL;
}


//功能说明:最速下降法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gx--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int fasterdown(double e,
			   int (* get_ak)(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),
			   int (* get_gk)(double *para,int paranum,double *gk),
			   double (* get_funcval)(double *para,int paranum),
			   double *para,int paranum)
{
	int stepcount=0,paracount;
	double ak=0,*gk,*pk,*para2;
	char *position="fasterdown";

	if (!(gk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;

	memcpy(para2,para,sizeof(double)*paranum);
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		stepcount++;
		get_gk(para2,paranum,gk);//取梯度值
		//取pk,即-gk
		for (paracount=0;paracount<paranum;paracount++)
			pk[paracount]=-gk[paracount];

		if (vector_val(gk,paranum)<e)//当前已到精度要求
		{
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount);			
			free(pk);
			free(para2);
			free(gk);
			return SUCCESS;
		}
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		get_ak(&ak,U,SIGMA,pk,para2,get_gk,get_funcval,paranum);//计算步长
		for (paracount=0;paracount<paranum;paracount++)
			para2[paracount]+=pk[paracount]*ak;
	}
	return FAIL;
}

//功能说明:阻尼Newton法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gk--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int newton(double e,
			   int (* get_ak)(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),
			   int (* get_gk)(double *para,int paranum,double *gk),
			   int (* get_G)(double *para,int paranum,double *G),
			   double (* get_funcval)(double *para,int paranum),
			   double *para,int paranum)
{
	int stepcount=0,paracount;
	double ak=0,*gk,*pk,*para2,*g,*reverseg;
	char *position="fasterdown";

	if (!(gk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(g=(double *)_malloc(position,sizeof(double)*paranum*paranum)))
		return FAIL;
	if (!(reverseg=(double *)_malloc(position,sizeof(double)*paranum*paranum)))
		return FAIL;

	memcpy(para2,para,sizeof(double)*paranum);
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		stepcount++;
		get_gk(para2,paranum,gk);//取梯度值

		if (vector_val(gk,paranum)<e)//当前已到精度要求
		{
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount);			
 			free(pk);
			free(para2);
			free(gk);
			free(reverseg);
			free(g);
			return SUCCESS;
		}
		//取G
		get_G(para2,paranum,g);
		if (!reverse_matrix(g,reverseg,paranum,paranum))//取逆阵不成功
			return FAIL;
		//取pk,即-G-1gk
		matrix_multiply(reverseg,paranum,paranum,gk,paranum,1,pk);
		//取pk,即-gk
		for (paracount=0;paracount<paranum;paracount++)
			pk[paracount]=-pk[paracount];
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		get_ak(&ak,U,SIGMA,pk,para2,get_gk,get_funcval,paranum);//计算步长
		for (paracount=0;paracount<paranum;paracount++)
			para2[paracount]+=pk[paracount]*ak;
	}
	return FAIL;
}

//功能说明:fr法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gx--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//		   needreturn--是否n步重新开始
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int fr(double e,
	   int (* get_ak)(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),
	   int (* get_gk)(double *para,int paranum,double *gk),
	   double (* get_funcval)(double *para,int paranum),
	   double *para,int paranum,
	   int needreturn)
{
	int steepcount=0,paracount,stepcount=1;
	double ak=0,*gk1,*gk2,*pk,*para2,b=0;
	char *position="fr";

	if (!(gk1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;

	memcpy(para2,para,sizeof(double)*paranum);

	for (paracount=0;paracount<paranum;gk1[paracount]=0.0,paracount++);
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		get_gk(para2,paranum,gk2);//取梯度值
		
		if (vector_val(gk2,paranum)<e)//当前已到精度要求
		{
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount+1);			
 			free(gk1);
			free(gk2);
			free(pk);
			free(para2);
			return SUCCESS;
		}
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		//取b的值
		if (stepcount++<2)
			b=0;
		else
			b=vector_inter_multiply(gk2,gk2,paranum)/vector_inter_multiply(gk1,gk1,paranum);
		//取pk,即-gk
		if ((needreturn)&&(((stepcount+1)%(paranum+1))==0))
		{
			for (paracount=0;paracount<paranum;paracount++)
				pk[paracount]=-gk2[paracount];
		}
		else
		{
			for (paracount=0;paracount<paranum;paracount++)
				pk[paracount]=-gk2[paracount]+b*pk[paracount];
		}
		memcpy(gk1,gk2,sizeof(double)*paracount);
		memset(gk2,0,sizeof(double)*paracount);
		get_ak(&ak,U,SIGMA,pk,para2,get_gk,get_funcval,paranum);//计算步长
		for (paracount=0;paracount<paranum;paracount++)
			para2[paracount]+=pk[paracount]*ak;
	}
	return FAIL;
}

//功能说明:prp法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gx--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//		   needreturn--是否n步重新开始
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int prp(double e,
	   int (* get_ak)(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),
	   int (* get_gk)(double *para,int paranum,double *gk),
	   double (* get_funcval)(double *para,int paranum),
	   double *para,int paranum,
	   int needreturn)
{
	int steepcount=0,paracount,stepcount=1;
	double ak=0,*gk1,*gk2,*pk,*tempgk,*para2,b=0;
	char *position="fr";

	if (!(gk1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(tempgk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;

	memcpy(para2,para,sizeof(double)*paranum);

	for (paracount=0;paracount<paranum;gk1[paracount]=0.0,paracount++);
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		get_gk(para2,paranum,gk2);//取梯度值
		
		if (vector_val(gk2,paranum)<e)//当前已到精度要求
		{
			free(gk1);
			free(gk2);
			free(pk);
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount+1);			
 			return SUCCESS;
		}
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		//取b的值
		if (stepcount++<2)
			b=0;
		else
		{
			for (paracount=0;paracount<paracount;paracount++)
				tempgk[paracount]=gk2[paracount]-gk1[paracount];
			b=vector_inter_multiply(gk2,tempgk,paranum)/vector_inter_multiply(gk1,gk1,paranum);
		}
		//取pk,即-gk
		if ((needreturn)&&(((stepcount+1)%(paranum+1))==0))
		{
			for (paracount=0;paracount<paranum;paracount++)
				pk[paracount]=-gk2[paracount];
		}
		else
		{
			for (paracount=0;paracount<paranum;paracount++)
				pk[paracount]=-gk2[paracount]+b*pk[paracount];
		}
		memcpy(gk1,gk2,sizeof(double)*paracount);
		memset(gk2,0,sizeof(double)*paracount);
		get_ak(&ak,U,SIGMA,pk,para2,get_gk,get_funcval,paranum);//计算步长
		for (paracount=0;paracount<paranum;paracount++)
			para2[paracount]+=pk[paracount]*ak;