# Оптимизация формы

## 4. Текст функции ConstrBFGS

    function [xmin,ymin,gmin,iter] = ...
      ConstrBFGS(fun,x0,A,b,tbeg,tol,dd,maxc,retone,varargin)
    % ConstrBFGS- minimization of function of several variables
    %             by Broyden-Fletcher-Goldfard-Shanno method
    %             with linear restrictions in kind Ax<=b.
    % fun       - name of function;
    %             the function must be in kind:
    %             [y,g] = myfunc(x,varargin);
    %             where:
    %             x  - the vector-column of arguments;
    %             y  - the value of function in point x;
    %             g  - the vector-column of gradient;
    % x0        - the initial arguments;
    % A, b      - the parameters of restrictions;
    % tbeg      - the initial length of the step;
    % tol       - the tolerance by arguments;
    % dd        - the diagonal element in Hessian;
    % maxc      - maximal number of cicles of n iterations;
    %             (the number of calculations of function =3*maxiter);
    % retone    - after retone iterations we return to diagonal Hessian;
    % Function return:
    % xmin - the arguments in point of minimum;
    % ymin - min value of function;
    % gmin - the gradient in point of minimum;
    % iter - real number of iterations.

    % Initial values;
    [m,n] = size(A); % number of variables and rextrictions;
    xmin=x0; % initial values of arguments;
    [ymin,gmin]=feval(fun,x0,varargin{:}); % function and gradient;
    H=eye(n)*dd; % the diagonal Hessian
    maxiter=maxc*n; % the maximal number of iterations;
    nupd=0; % if need update the Hessian;
    tcurr=tbeg; % current length of the step

    % We begin the iterations
    for iter=1:maxiter,
      fprintf('The iteration %d, objective function y=%16.8f\n',...
        iter,ymin);
      x=xmin % current values
      g=gmin;
      if norm(g)<tol, % tolerance by gradient
        disp('Gradient is low! Exit.');
        return;
      end;
      if nupd, % if need update the Hessian
        dx=x-xold; % the increment of the argument
        if norm(dx)<tol, % tolerance by argument;
          disp('The increment of the argument is low! We in min - exit.');
          return;
        end;
        dg=g-gold; % the increment of the gradient
        H=H+(dg*dg')/(dg'*dx)-(H*dx*dx'*H)/(dx'*H*dx); % BFGS-formula
      end; % if nupd
      s=-H\g; % the quazi-Newton direction
      s=s/norm(s); % normalize
      % we project vector of a direction on active restrictions
      nb=find((abs(b-A*x)<1e4*eps)&(A*s>0)); % numbers of active restrictions
      if ~isempty(nb), % active restrictions is present
        for k=1:length(nb), % for each active restriction
          na=A(nb(k),:); % row of matrix of rextrictions A
          s=s-((na*s)/(na*na'))*na'; % we is subtrahend a projection to a normal
        end
      end
      if norm(s)<1e4*eps, % we is in the angle point
        disp('Angle point is reached.');
        if ~nupd, % Hessian is diagonal
          disp('Hessian is diagonal. Exit.');
          return;
        end
        disp('We return to diagonal Hessian.');
        H=eye(n)*dd;
        nupd=0; % not need update the Hessian
      else
        s=s/norm(s); % normalize
        bi=(b-A*x); % right members
        aki=A*s; % factors
        nap=find(aki>1e4*eps); % numbers of positive coordinates
        if isempty(nap), % not the positive coordinates
          tmax=inf; % it is possible increase unrestrictebly
        else
          tmax=min(abs(bi(nap))./aki(nap)); % min relation
        end
        % 1D search
        tt=zeros(1,3); % 3 arguments
        ff=zeros(1,3); % 3 functions
        gg=zeros(n,3); % 3 gradients
        gs=zeros(1,3); % 3 projections of gradient
        ff(1)=ymin; % 1-st point - current (have calculated);
        gg(:,1)=g;
        gs(1)=gg(:,1)'*s;
        if gs(1)>0,
          disp('Bad direction.')
          if ~nupd, % Hessian is diagonal
            disp('Hessian is diagonal. Exit.');
            return;
          end
          disp('Return to diagonal Hessian.');
          H=eye(n)*dd;
          nupd=0; % the Hessian will not be necessary to update
        else
          tt(2)=min(tcurr,tmax); % 2-nd point
          xcurr=x+tt(2)*s;
          [ff(2),gg(:,2)]=feval(fun,xcurr,varargin{:}); % 2-nd point
          gs(2)=gg(:,2)'*s;
          while (gs(2)<0)&(ff(2)<ff(1)), % search along direction
            if tt(2)>tmax-1e4*eps, % only to nearest bound
              break;
            end
            tt(3)=tt(1);
            tt(1)=tt(2);
            ff(1)=ff(2);
            gg(:,1)=gg(:,2);
            gs(1)=gs(2);
            tt(2)=min(tmax,tt(2)+2*(tt(1)-tt(3)));
            tcurr=tt(2);
            xcurr=x+tt(2)*s;
            [ff(2),gg(:,2)]=feval(fun,xcurr,varargin{:}); % 2-nd point
            gs(2)=gg(:,2)'*s;
          end
          if tt(2)>tmax-1e4*eps, % in border
            tt(3)=tmax; % 3-rd point is on the border
          else % not border - cubic interpolation
            b1=gs(1)+gs(2)-3*(ff(1)-ff(2))/(tt(1)-tt(2));
            b2=(b1^2-gs(1)*gs(2))^0.5;
            tt(3)=tt(2)-(tt(2)-tt(1))*(gs(2)+b2-b1)/(gs(2)-gs(1)+2*b2);
          end
          [ff(3),gg(:,3)]=feval(fun,x+tt(3)*s,varargin{:}); % 3-rd point
          gs(3)=gg(:,3)'*s;
          [ymin,imin]=min(ff); % minimal value and it number
          if (imin==1)&(tt(1)==0), % minima is in initial point
            disp('The function not decrease.');
            if ~nupd, % Hessian is diagonal
              disp('Hessian is diagonal. Decrease step.');
              tcurr=tcurr/2;
              if tcurr<tol
                disp('Step is low. We in minima point - exit.');
                return;
              end;
            end
            disp('We return to diagonal Hessian.');
            H=eye(n)*dd;
            nupd=0; % not need update the Hessian
          else
            tcurr=tt(imin);
            xold=xmin; gold=gmin; yold=ymin; % we save old values
            xmin=x+tt(imin)*s; % argumets
            gmin=gg(:,imin); % gradients
            if round(rem(iter,retone))==0, % full cicle of retone iterations
              H=eye(n)*dd;
              disp('The full cicle is held - we return to a diadonal Hessian');
              nupd=0; % not need update the Hessian
            else
              nupd=1; % need update the Hessian
            end;
          end; % there is an advance to a minimum
        end;
      end;
    end % while
    return