function [k,x,val]=bfgs(fun,gfun,x0,varargin)
%输入:fun,gfun 分别是目标函数与目标函数的梯度,x1
%是初始目标点,varargin是输入的可变参数变量
%输出:k迭代次数,x,val近似最优点和最优值
%N=500;
N=1000;
epsilon = 1.e-5;%给定容许误差
beta = 0.55;
sigma = 0.4;
n=length(x0);
Bk=eye(n);
k=0;
while(k<N)
gk = feval(gfun,x0);%计算梯度
if(norm(gk)<epsilon)
break;
end
dk=-Bk\gk;%解方程Bk*dk=-gk;,计算搜素方向
m=0;
mk = 0;
while(m<20)
newf = feval(fun,x0+beta^m*mk);
oldf = feval(fun,x0);
if(newf<=oldf+sigma*beta^m*gk'*dk)
mk = m;
break;
end
m=m+1;
end
%BFGS校正
size(dk)
x=x0+beta^mk*dk;
sk=x-x0;
yk = feval(gfun,x)-gk;
if(yk'*sk>0)
Bk = Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);
end
k=k+1;
x0 = x;
end
val=feval(fun,x0);
x0
val
k
end