bfgs 算法 拟牛顿

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