最大熵模型学习的BFGS算法c++实现

BFGS算法步骤：

注意：代码需要在VS中编译我用devc++编译可能不通过

#include <iostream>
//#include <iomanip.h>
#include <math.h>
using namespace std;
const double gama=2;
const double lamda=0.618033989;
const double h=0.001;
const double eps=10e-5;
double X0[4]={1,2,2,2};
double I[16]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};

class Matrix
{
public:
	double *elems;
	int row, col;
public:
	Matrix(int r, int c); 
	Matrix(double *m, int r, int c);  
	~Matrix(){}
	Matrix& operator = (Matrix &b); 
	Matrix trans();
	double norm(); 
	double MTod(); 
};
Matrix operator + (Matrix &a, Matrix &b); 
Matrix operator - (Matrix &a, Matrix &b); 
Matrix operator * (double a, Matrix &b); 
Matrix operator * (Matrix &a, Matrix &b); 


double f(double x1, double x2,double x3,double x4) 
{
	return (exp(x1)-x2)*(exp(x1)-x2)*(exp(x1)-x2)*(exp(x1)-x2)+100*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)+(tan(x3-x4))*(tan(x3-x4))*(tan(x3-x4))*(tan(x3-x4))+x1*x1*x1*x1*x1*x1*x1*x1;
}
double df1(double x1, double x2,double x3,double x4) 
{
	return 4*(exp(x1)-x2)*(exp(x1)-x2)*(exp(x1)-x2)*exp(x1)+8*x1*x1*x1*x1*x1*x1*x1;
}
double df2(double x1, double x2,double x3,double x4) 
{
	return -4*(exp(x1)-x2)*(exp(x1)-x2)*(exp(x1)-x2)+600*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3);
}
double df3(double x1, double x2,double x3,double x4)
{
	return -600*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)*(x2-x3)+4*tan(x3-x4)*tan(x3-x4)*tan(x3-x4)/(cos(x3-x4)*cos(x3-x4)*cos(x3-x4));
}
double df4(double x1, double x2,double x3,double x4)
{
	return -4*tan(x3-x4)*tan(x3-x4)*tan(x3-x4); 
}
Matrix::Matrix(int r, int c)
{
	if( r>0 && c>0 )
	{
		row=r; 
		col=c;
		elems=new double[r*c];
	}
	else
	{
		elems=NULL;
		row=col=0;
	}
}
Matrix::Matrix(double *m, int r, int c)
{
	if( r>0 && c>0 )
	{
		row=r; 
		col=c;
		elems=new double[r*c];
	}
	else
	{
		elems=NULL;
		row=col=0;
	}
	for (int i=0; i<r*c; i++)
		elems[i]=m[i];
}
Matrix operator + (Matrix &a, Matrix &b)
{
	Matrix *ans=new Matrix(a.row,a.col);
	for (int r=0; r<a.row; r++)
	{
		for (int c=0; c<a.col; c++)
			ans->elems[r*ans->col+c]=a.elems[r*a.col+c]+b.elems[r*b.col+c];
	}
	return *ans;
}
Matrix operator - (Matrix &a, Matrix &b)
{
	Matrix *ans=new Matrix(a.row,a.col);
	for (int r=0; r<a.row; r++)
	{
		for (int c=0; c<a.col; c++)
			ans->elems[r*ans->col+c]=a.elems[r*a.col+c]-b.elems[r*b.col+c];
	}
	return *ans;
}
Matrix operator * (Matrix &a, Matrix &b)
{
	Matrix *ans=new Matrix(a.row,b.col);
	double s;
	for (int r=0; r<a.row; r++)
	{
		for (int c=0; c<b.col; c++)
		{
			s=0;
			for (int k=0; k<a.col; k++)
			{
				s+=a.elems[r*a.col+k]*b.elems[k+c*b.row];
			}
			ans->elems[r*b.col+c]=s;
		}
	}
	return *ans;
}
Matrix operator * (double a, Matrix &b)
{
	Matrix *ans=new Matrix(b.row,b.col);
	for (int r=0; r<b.row; r++)
	{
		for (int c=0; c<b.col; c++)
			ans->elems[r*b.col+c]=a*b.elems[r*b.col+c];
	}
	return *ans;
}
Matrix & Matrix::operator = (Matrix &b)
{
	for (int r=0; r<row; r++)
	{
		for (int c=0; c<col; c++)
			elems[r*col+c]=b.elems[r*col+c];
	}
	return *this;
}
Matrix Matrix::trans()
{
	Matrix *ans=new Matrix(col,row);
	for (int r=0; r<row; r++)
	{
		for (int c=0; c<col; c++)
			ans->elems[c*row+r]=elems[r*col+c];
	}
	return *ans;
}
double Matrix::norm()
{
	int i;
	double sum=0;
	if (col==1)
	{
		for(i=0; i<row; i++)
			sum+=fabs(elems[i*col+0]);
	}
	return sum;
}
double Matrix::MTod()
{
	if(row==1 && col==1)
		return elems[0];
}


Matrix tiduf(double x1,double x2,double x3,double x4) 
{
	Matrix *ans=new Matrix(4,1);
	double df[4]={df1(x1,x2,x3,x4),df2(x1,x2,x3,x4),df3(x1,x2,x3,x4),df4(x1,x2,x3,x4)};
	*ans=Matrix(df,4,1);
	return *ans;
}
void ssqj(double &a, double &b, Matrix x, Matrix d) 
{
	double a0=0.1, g1=gama, g2=gama, a1=a0+h, t, a2;
	if(f((x+a0*d).elems[0],(x+a0*d).elems[1],(x+a0*d).elems[2],(x+a0*d).elems[3])>f((x+a1*d).elems[0], (x+a1*d).elems[1],(x+a1*d).elems[2],(x+a1*d).elems[3]) )
	{
		a2=a1+g1*h;
		while(f((x+a2*d).elems[0],(x+a2*d).elems[1],(x+a2*d).elems[2],(x+a2*d).elems[3])<=f((x+a1*d).elems[0], (x+a1*d).elems[1],(x+a1*d).elems[2], (x+a1*d).elems[3]) )
		{
			g1=g1*gama;
			a0=a1;
			t=a2;
			a2=a1+g1*h;
			a1=t;
		}
		a=a0;
		b=a2;
	}
	else
	{
		a2=a0-g2*h;
		while(f((x+a2*d).elems[0],(x+a2*d).elems[1],(x+a2*d).elems[2],(x+a2*d).elems[3])<=f((x+a0*d).elems[0], (x+a0*d).elems[1], (x+a0*d).elems[2], (x+a0*d).elems[3]) )
		{
			g2=g2*gama;
			a1=a0;
			t=a2;
			a2=a0-g2*h;
			a0=t;
		}
		a=a2;
		b=a1;
	}
}
double hjfg(double a, double b, Matrix x, Matrix d) 
{
	double a1, a2;
	a1=b-lamda*(b-a);
	a2=a+lamda*(b-a);
	while ( (b-a)>=eps )
	{
		if(f((x+a1*d).elems[0],(x+a1*d).elems[1],(x+a1*d).elems[2],(x+a1*d).elems[3])<=f((x+a2*d).elems[0], (x+a2*d).elems[1],(x+a2*d).elems[2], (x+a2*d).elems[3]) )
		{
			b=a2;
			a2=a1;
			a1=b-lamda*(b-a);
		}
		else
		{
			a=a1;
			a1=a2;
			a2=a+lamda*(b-a);
		}
	}
	if(f((x+a1*d).elems[0],(x+a1*d).elems[1],(x+a1*d).elems[2],(x+a1*d).elems[3])<f((x+a2*d).elems[0], (x+a2*d).elems[1],(x+a2*d).elems[2], (x+a2*d).elems[3]) ) 
		return a1;
	else 
		return a2;
}
double alphak(Matrix x, Matrix d)
{
	double a, b;
	ssqj(a,b,x,d); 
	return hjfg(a,b,x,d);
}
int main()
{
	Matrix x0(X0,4,1), H0(I,4,4), x1(4,1), H1(4,4), d(4,1), delta(4,1), y(4,1);
	double alpha;
	while (1)
	{
		d=(-1)*H0*tiduf(x0.elems[0],x0.elems[1],x0.elems[2],x0.elems[3]);
		alpha=alphak(x0, d);
		x1=x0+alpha*d;
		if(tiduf(x0.elems[0],x0.elems[1],x0.elems[2],x0.elems[3]).norm()<eps||(x1-x0).norm()<eps )
		{
			cout << "该问题的最优解为:"<<endl;
			cout<<"x=("<<x1.elems[0] <<' '<< x1.elems[1] <<' '<<x1.elems[2]<<' '<< x1.elems[3]<<")\nf(x)=" << f(x1.elems[0],x1.elems[1],x1.elems[2],x1.elems[3]) << endl;
			return 0;
		}
		else
		{
			delta=x1-x0;
			y=tiduf(x1.elems[0],x1.elems[1],x1.elems[2],x1.elems[3])-tiduf(x0.elems[0],x0.elems[1],x0.elems[2],x0.elems[3]);
			H1=H0+(1+(y.trans()*H0*y).MTod()/(y.trans()*delta).MTod())/(y.trans()*delta).MTod()*delta*delta.trans()-1/(y.trans()*delta).MTod()*(H0*y*delta.trans()+(H0*y*delta.trans()).trans());
			H0=H1;
			x0=x1;
		}
	}
	return 0;
}