relnewton.m

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

function [W,res]=relnewton(X,lam,niter,normg_tol,method)
% Relative Newton method for ICA 
%
% Use:
%
% W=relnewton(X)
% W=relnewton(X,lam)
% W=relnewton(X,lam,niter)
% W=relnewton(X,lam,niter,normg_tol)
% W=relnewton(X,lam,niter,normg_tol,method)
%[W,res]=relnewton(X,lam,niter,normg_tol,method)
%
%
%Input:
%
%  X         - matrix with the mixed signals (row-wise)
%
%  lam       - parameter of nonlinearity, default: lam=1e-2
%
%                 lam > 0 - use  smoothing of the absolute value function with parameter lambda,
%                                becomes more sharp and accurate  when lam -> 0;
%                                recommended for sparse/sepergaussian sources 
%                                (see the paper newt_ica_ica2003.pdf)
%
%                 lam < 0 - use  the power function |s|^|lam| (lam= -4 ... -20),  
%                                recommended for subgaussian sources
%                
%  niter     - Maximum number of Newton steps. Default: 200
%
%  normg_tol - norm of the gradient for the stopping criterion. Default: 1e-5
%
%  method    - 'fastnwt' - fast relative Newton (default)
%              'relnwt'  - standard relative Newton
%              'natgrad' - natural gradient (batch mode)
%              'diagprecond' -  diaginally preconditioned natural gradient (batch mode)
%
%Output:
%
%  W   - separating matrix:   recovered sources S = WX
%  res - report of the iterative process

%%%%%%%%%%%%%%%%%%%%%%%%  by Michael Zibulevsky, version: 30.10.2002 %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%  by Michael Zibulevsky, version: 03.12.2003 %%%%%%%%%%%%%%%%%%



pcaflag= 0;   % 1 - use PCA preprocessing (accelerates convergence, but you can loose 
              %                            superefficiency for sparse sources !!!)

armijoflag=1; % 1 - armijo step (recommended), 0 - cubic linesearch


fghpar.penpar=lam;  % parameter of smoothing of absolute value function
fghpar.pentype='smooth_abs'; % Type of non-linearity
fghpar.fghname='fghd_lagr'; % Don't ask about this parameter: just write

fghname=fghpar.fghname;

  
  if nargin <2
    lam=1e-2;      % smoothing of the absolute value function with parametre lambda
  end


  if nargin <3
    niter=200;      % Maximum number of newton steps
  end
  
   if nargin <4
     normg_tol=1e-4; % norm of the gradient for the stopping criterion.
   end
   
   if nargin <5
     method='fastnwt';
   end
   
 
  [Nsn,TotT]=size(X);
  
  I=eye(Nsn);
  W=I;
 
  fold=1e10;
  normg=1e10;
  %res.df=[];
  res.normg=[];
  res.dw=[];
  res.f=[];
  res.time=[];
  time0=cputime;
  
  if pcaflag
    [X,B]=mypca(X);  % normalized principal components X=B*X
  else
    B=eye(Nsn);
  end
  

  for i=1:niter   %%%%%%%%%%%% Optimization loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     U=W*X;
     w=I(:);
     
     if strcmp(method,'relnwt')  %%%%%%%%%%%%%%%% STANDARD RELATIVE NEWTON %%%%%%%%%%
         [f,g,hess_diag,hess]=feval(fghname,w,U,fghpar); 
         %[f,g,hess]=fgh_stand(w,U,lam);   %hess
         %[f,g,hess_diag,hess]=fghdiag(w,U,lam);   %Approximate hess: diag + logdet''
         
         if 1   
             [CholFact,D]=mcholmz3(full(hess)); CholFact=CholFact'; 
             d= cholsub(CholFact,-g,D); % Cholessky substitution
         else
             d=zeros(Nsn^2,1);
             for j=1:Nsn      % Solve block-diagonal Newton system
                 k0=(j-1)*Nsn; kk= [(k0+1):(k0+Nsn)];
                 d(kk)= -hess(:,:,j)\g(kk);
             end
         end

     elseif strcmp(method,'fastnwt')  %%%%%%%%%%%%%%% FAST RELATIVE NEWTON %%%%%%%%%%%%%%%%%%
        
         [f,g,hess_diag]=feval(fghname,w,U,fghpar); 
         
         %%%%%%%%%% Fast solution of Newton matrix equation  D*X + X' = G %%%%%%%%%%
         
         D=reshape(hess_diag,Nsn,Nsn)';
         G=reshape(g,Nsn,Nsn)';
         dd=fast_newt_dir(D,-G); d=dd'; d=d(:);
         
         %dd= -(G.*D'-G') ./ (D.*D'-ones(Nsn)); d=dd'; d=d(:);
         % D.*dd+dd'+G; d1= -hess\g; dmd1=d-d1, keyboard
     
     elseif strcmp(method,'natgrad') 
        [f,g]=feval(fghname, w,U,fghpar);   %hess
        d=-g;
     elseif strcmp(method,'diagprecond') 
        [f,g,hess_diag]=feval(fghname, w,U,fghpar); 
        d= -g./(hess_diag+1); % Preconditioned gradient descent direction
     else   
        error('Unknown optimization method') 
     end

      
     tbest=1;
     
     if normg > 1e-5 | strcmp(method,'natgrad') %%%%%%%% Line search %%%%%%%%%%%%%%%%%

       if armijoflag, %% Armijo step (recommended)
	 
	      [wnew,fnew,gradnew,tbest]= armijostep(fghname,w,f,g,d,0.3,0.3, U,fghpar);
       
       else  %%%  cubic line search

	     b1=0;b2=1;EpsGold=0.4; dftol=1e-3; dxtol=1e-3; nlinsrchmax=10;
         [wnew,fnew,gradnew,tbest,resEpsGold,b1,b2] = cublinsrch1(fghname,...
             w,f,g,d,b1,b2,EpsGold,dftol,dxtol,nlinsrchmax,  U,fghpar);
       end

     else  % make full step when newton method is close to solution
       
       wnew=w+d;  
       [fnew,gradnew] = feval(fghname,wnew,U,fghpar);   % new func. and grad
     
     end  %%%%%%%%%%%%%%%%%%%%%%%% end of Line search %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     fprintf('tbest=%g\n',tbest);
      
     Wtmp=reshape(wnew,Nsn,Nsn)';

     GRAD=reshape(gradnew,Nsn,Nsn)';
     %NATGRAD=GRAD*Wtmp'*Wtmp;

     W=Wtmp*W;  % Update separation matrix
     
     %WW=W'; f0=feval(fghname,WW(:),X,fghpar);  df=fold-f0; fold=f0;   res.df=[res.df;df];

     
     dW=Wtmp-I;
     maxdW=max(abs(dW(:)));
     normg=norm(gradnew);
     res.normg=[res.normg;normg];
     res.dw=[res.dw;maxdW];
     res.f=[res.f;fnew];
     res.time=[res.time;cputime-time0];
     fprintf('iter %d: f=%g  normg=%.2e,  dW=%.2e    ',i, fnew,normg,maxdW);

     figure(200); semilogy([1:i],res.normg,'*'); %,[1:i],res.dw,'o'); 
     title('Relative Newton iterations: gradient norm'); xlabel('Iteration'); drawnow
      
     if normg < normg_tol, break; end
     
     if i==niter, fprintf('WARNING: MAXIMAL NUMBER OF ITERATIONX EXEEDED IN RELNEWTON.M !!!');end
     
  end    %%%%%%%%%%%% end of optimization loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  W=W*B; % Account for the PCA, which was done  at the beginning
    
return


%%%%%%%%%%%%%%%%%%%%%%% TEST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%% Generate the data %%%%%%%%
N=3;
T=1000;
A=rand(N);
S=sprandn(N,T,0.5);
X=A*S;

%%%%%% Setup the parameters

lam=1e-1; % parameter of smoothing of the absolute value function
niter=100
normg_tol=1e-12;
%method='relnwt';  
method= 'fastnwt'; 


%%%%%%% Separate with the relative NEWTON %%%%%%%%
W=relnewton(X,lam,niter,normg_tol,method)

%%%%%%%% COMPUTE QUALITY OF SEPARATION
WA=W*A
ISR=max(sep_quality(WA))