armijostep.m

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

function [x,f,grad,t] = armijostep(fg,x0,f0,grad0,dx,alpha,beta,varargin)

%armijostep - Descent of the function f(x) in direction dx:
%
%                     f(x0 + t*dx) <= f(x0) - alpha*|f'_dx (x0)|

%CALLED AS:
%==========
%
%        [x,f,grad,t,res_alpha] = armijostep(fg,x0,f0,grad0,dx,alpha,beta,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
% alpha - required descent rate:  f(x0 + t*dx) <= f(x0) - alpha*|f'_dx (x0)|
% beta  - factor of step decreasing: t_new = beta*t_old
%
% 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


  %Some parameters:
 
  nitermax=30;   % - max number of iterations
  t=1;           %  - initial stepsize


  df0=grad0'*dx; 
  if df0>0, 
    disp('Armijo search: Direction dx is of increase, changing it to the opposite !!!');
    dx=-dx; df0=-df0; 
  end

  x=x0+t*dx;
  [f,grad] = feval(fg,x,varargin{:});
  df=grad'*dx; 
  %if df<=0, return; end
    
  for i=1:nitermax,
    %fprintf('armijostep ');keyboard
    
    if (f <= f0 + alpha*df0*t) | ((df < 0) & (f < f0 + 1e-9 * max(1,abs(f0)))), return; end;
    %    if (f <= f0 + alpha*df0*t) | (df<0 & f<f0), return; end;
  
    fprintf('-');
    t= beta*t;
    x=x0+t*dx;
    [f,grad]=feval(fg,x,varargin{:}); 
    df=grad'*dx; 
    %if t<1e-6, fprintf('\n Armijostep.m); keyboard,end
  end    

  if f>f0, disp('Armijostep.m : f>f0 !!!');keyboard;end

return

%%%%%%%%%%%%%%% Garbage
ff=[];
ttt=[-1:0.01:1];
for tt= ttt
    x=x0+t*dx;
    ff=[ff  feval(fg,x,varargin{:})];
end
figure;plot(ttt,ff);