bfgs.m

这是在网上下的一个东东

%
% xstar=bfgs(func,delf,x0,tol,maxiter)
%
% Use steepest descent to minimize a function f(x).
%
%    func         name of the function f(x)
%    delf         name of the gradient del f(x)
%    x0           initial guess
%    tol          stopping tolerance
%    maxiter      maximum number of iterations allowed
%
function xstar=bfgs(func,grad,x0,tol,maxiter)
global FUN;
global X;
global P;
%
%  Figure out the size of the problem.  
%
n=size(x0,1);
%
%  Initialize x and oldx.  We use oldx=1000*x just to make sure that
%  it isn't too close to x.  
%
x=x0;
fx=feval(func,x);
oldfx=fx*10;
oldx=x*10;
%
% The main loop.  While the current solution isn't good enough, keep
% trying...  Stop after maxiter iterations in the worst case. 
%
%
%  Initialize the iteration count.
%
  iter=0;
%
% Initialize B and its Cholesky factorization.
%
  B=eye(n);
  [L,D]=ldlt(B);
%
% Get the gradient at the initial point.
%
  g=feval(grad,x);
%
%  The loop.
%
  while (iter <maxiter)
%
% Check the termination criteria.
%
    if ((norm(g,2)< sqrt(tol)*(1+abs(fx))) & ...
        (abs(oldfx-fx)<tol*(1+abs(fx))) & ...
        (norm(oldx-x,2)<sqrt(tol)*(1+norm(x,2))))
%      disp('Problem solved');
%      iter
      xstar=x;
      return;
    end
%
%  Compute the new search direction by solving:
%           Bp=-g
%
    u=L\(-g);
    v=D\u;
    p=L'\v;
%
%  Setup to minimize f(x+alpha*p)
%
    X=x;
    P=p;
    FUN=func;
%
%  Minimize along the line.
%
    [a,u,b]=bracket('subfun');
    [alphamin,falphamin,newa,newb]=gssearch('subfun',a,u,b,1.0e-4);
%
%  The new point is at X+alphamin*P.  Update oldx to.
%
    oldx=x;
    oldfx=fx;
    x=X+alphamin*P;
    fx=falphamin;
%
% Compute the gradient at the new point.
%
    oldg=g;
    g=feval(grad,x);
    y=g-oldg;
%
%  Update B.
%
    B=B+(oldg*oldg')/(oldg'*p)+y*y'/(alphamin*y'*p);
%
%  Update the factorization of B.
%
    v=oldg/sqrt(-oldg'*p);
    w=y/sqrt(alphamin*y'*p);
    [L,D]=ldltup(L,D,w);
    [L,D]=ldltdown(L,D,v);
%
%  If the matrix is close to singular, then restart with B=I.  Note that
%  If we have any NaN's in B, then we really need to restart too.
%
    if (max(diag(D))/min(diag(D)) > 1.0e15) | (sum(sum(isnan(B)))>0)
      B=eye(n);
      [L,D]=ldlt(B);
%      disp('restarting with B=I');
%      iter
    end
%
%  Update the iteration counter.
%
    iter=iter+1;
%
%  End of the loop.
%
  end
%
%  Save the results.
%
xstar=x;
%
%  Print out some information
%
disp('BFGS method failed');
disp('Iterations were');
iter