最大熵模型学习的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;
}