cublinsrch1.m

该代码为盲源分离中的相对牛顿法

function [x,f,grad,t,resEpsGold,b1,b2] = cublinsrch1(fg,...
        x0,f0,grad0,dx,b1,b2,EpsGold,gtol,xtol,nitermax,varargin)

%cublinsrch1 - cubic linesearch with bisection safeguard and early
%               stopping by Goldstain criterion
%
%               We  minimize a function f(x) in direction dx:
%
%                     min_t f(x0 + t*dx)
%CALLED AS:
%==========
%
%        [x,f,grad,t,resEpsGold,b1,b2] = cublinsrch1(fg,...
%        x0,f0,grad0,dx,b1,b2,EpsGold,dir_drv_tol,delta_t_tol,...
%        nitermax,varargin)
%
%INPUT:
%======
%
% fg    - name of matlab function, which computes function f and 
%         gradient grad at the point x:
%
%                   [f,grad]=fg(x,varargin)
%
%                   varargin - an arbitrary number of additional arguments   
%          
%
% x0    - initial value of the optimized variable x
% f0    - value  of the optimized function at the point x0
% grad0 - value  of the gradient at the point x0
% dx    - direction of search
% b1    - initial lower bound on t (usually b1=0)
% b2    - initial upper bound on t (usually b2=1)
%
% Stopping occures if one of 4 criteria is achieved:
%
% EpsGold     - Goldstain stopping criterion. Stop if
%                     abs((f-f0)/(t*f'(t))-0.5) < EpsGold
% dir_drv_tol - stopping ctiterion in the value of directional derivative f'(t)
% delta_t_tol - stopping ctiterion in the  size of unsertain interval in t
% nitermax    - max number of iterations
%
%
% varargin - an arbitrary number of additional arguments,which are 
%            passed to fg(...)
%
%OUTPUT:
%=======
% x          - resulting x
% f          - resulting f(x)
% grad       - resulting gradient of f(x)
% t          - resulting search parameter: x = x0 + t*dx
% resEpsGold - resulting value of Goldstain criterion
% b1         - resulting lower bound on t
% b2         - resulting upper bound on t

  %EpsGold=0; gtol=0; xtol=0 %for debugging

  %Some parameters:
  Delta=1e-8;
  alpha= 0.7; % required minimal progress of cubic step;
  mu=0.2;    % artificial increase in cubic step to provide V-shape

  resEpsGold=0;

  b1=0;

  f1=f0; grad1=grad0;
  normdx= norm(dx);
  df0=grad0'*dx; 
  if df0>0, 
    disp('cublinsrch1: Direction dx is of increase, changing it to the opposite !!!');
    dx=-dx; df0=-df0; 
  end
  df1=df0;

  t=b2;
  x=x0+t*dx;
  [f,grad] = feval(fg,x,varargin{:});
  f2=f;df2=grad'*dx; df=df2;
  if df2<=0, return; end
    
  
  for i=1:nitermax,
    %fprintf('cublinsrch1 ');keyboard
    fprintf('-');
    %fprintf('1: b1=%g,b2=%g,t=%g,f=%g,df=%g\n', b1,b2,t,f,df); 

    resEpsGold=abs((f-f0)/(t*df0-1.765786e-80)-0.5);
    if(resEpsGold<EpsGold) break; end %Goldstain criterion
    if (abs(df/df0)<gtol | b2-b1<xtol) & f<f0,break;end

    tcub=b1+mycubici1(f2,f1,df2,df1,b2-b1);            %cubic interpolation 

    db=b2-b1; db_old=db; 
     
    
    if tcub>b1 & b2>tcub, %cubic step

      bisect_flag=0;
      
      if (tcub-b1 < b2-tcub), t =tcub+mu*(tcub-b1);
      else                    t =tcub-mu*(b2-tcub);
      end
      
      x=x0+t*dx;
      [f,grad] = feval(fg,x,varargin{:});
      df=grad'*dx;
      
      %fprintf('linesearch ratio %g\n',(t-b1)/db);
      
      if df<0, b1=t; f1=f; df1=df;
      else,    b2=t; f2=f; df2=df;
      end


      %fprintf('2: b1=%g,b2=%g,t=%g,f=%g,df=%g\n', b1,b2,t,f,df); 
      
      resEpsGold=abs((f-f0)/(t*df0-1.765786e-80)-0.5);
      if(resEpsGold<EpsGold) break; end %Goldstain criterion
      if (abs(df/df0)<gtol | b2-b1<xtol) & f<f0,break;end

      db=b2-b1;
      if(db>alpha*db_old), bisect_flag=1; end
    else
      bisect_flag=1;
    end

    if bisect_flag,

      fprintf('b');
      %disp('bisection step');

      t=(b1+b2)/2;

      x=x0+t*dx;
      [f,grad] = feval(fg,x,varargin{:});
      df=grad'*dx;


      if df<0, b1=t; f1=f; df1=df;
      else,    b2=t; f2=f; df2=df;
      end

    end

  end % end_for
 

function [x,f,grad,tmid] = bis(fg,x0,dx,t1,t2,tacc,nitermax,varargin),

 x=x0+t1*dx;
 [f,grad] = feval(fg,x,varargin{:});
 f0=f; f1=f; df1=grad'*dx;

 if (df1<=0)
   tmid=t1;
   disp('bis: bisection, direction of increase, zero step is taken!');
   return
 end

 x=x0+t2*dx;
 [f,grad] = feval(fg,x,varargin{:});
 f2=f;df2=grad'*dx;
 
 if (df2<=0)
   tmid=t2;
   return
 end
 
 for j=1:150,
   tmid=(t1+t2)/2.;

   x=x0+tmid*dx;
   [f,grad] = feval(fg,x,varargin{:});
   dfmid=grad'*dx;

   if(dfmid*df1>0),  t1=tmid; df1=fmid;
   else              t2=tmid; df2=fmid;
   end

   %fprintf('j=%d,dfmid=%g',j,dfmid);
   
   if(j>=jmax & f<f0), break; end

   %ngradbisglob = ngradbisglob+1;
 end

function [xbest,fbest,gradbest,tbest,resEpsGold] = cubbislinsrch(fg,...
        x0,f0,grad0,dx,b1,b2,EpsGold,gtol,xtol,nitermax,varargin),
  
  
  [xbest,fbest,gradbest,tbest,resEpsGold,b1,b2] = cublinsrch(fg,...
        x0,f0,grad0,dx,b1,b2,EpsGold,gtol,xtol,nitermax,varargin);

  if resEpsGold > 0.495,  %if_ Goldstain criterion not satis. then bisect
    [x,f,grad,tmid] = bis(fg,x0,dx,b1,b2,tacc,nitermax,varargin);
  end


function r=mycubici1(fnew,fold,graddnew,graddold,stepsize)
%CUBICI1 Cubicly interpolates 2 points and gradients to estimate minimum.
%
%   This function uses cubic interpolation and the values of two 
%   points and their gradients in order to estimate the minimum of a 
%   a function along a line.

if(abs(stepsize)<eps) stepsize= sign(stepsize)*eps;end
if isinf(fnew), fnew=1/eps; end
z=3*(fold-fnew)/stepsize+graddold+graddnew;
w=real(sqrt(z*z-graddold*graddnew));
r=stepsize*((z+w-graddold)/(graddnew-graddold+2*w));
return


%%%%%%%%%%%%%%% test mycubici1

y=inline( '5*x + 3*x.^2+ 0.3*x.^3');
dy=inline( '5 + 6*x+ 0.9*x.^2');
z=[-10:0.1:10];figure;plot(z,y(z));grid 

t1=-3; t2=5; t3= t1+mycubici1(y(t2),y(t1),dy(t2),dy(t1),t2-t1)
         

ii=0; 
ss=b1:0.001:1;
%ss=b1:0.02*db:b2;
ff=ss;dff=ss;
for s=ss,
  ii=ii+1;
  x=x0+s*dx;
  %[ff(ii),g]=feval(fg,x,varargin{:});
  dff(ii)=g'*dx;
end
figure;plot(ss,ff);grid
figure;plot(ss,dff);grid