jdqz.m

数据降维工具箱

function varargout=jdqz(varargin)
%JDQZ computes a partial generalized Schur decomposition (or QZ
%  decomposition) of a pair of square matrices or operators.
%  
%  LAMBDA=JDQZ(A,B) and JDQZ(A,B) return K eigenvalues of the matrix pair
%  (A,B), where K=min(5,N) and N=size(A,1) if K has not been specified.
%  
%  [X,JORDAN]=JDQZ(A,B) returns the eigenvectors X and the Jordan
%  structure JORDAN:  A*X=B*X*JORDAN. The diagonal of JORDAN contains the
%  eigenvalues: LAMBDA=DIAG(JORDAN). JORDAN is an K by K matrix with the
%  eigenvalues on the diagonal and zero or one on the first upper diagonal
%  elements. The other entries are zero.
%  
%  [X,JORDAN,HISTORY]=JDQZ(A,B) returns also the convergence history.
%  
%  [X,JORDAN,Q,Z,S,T,HISTORY]=JDQZ(A,B) 
%  If between four and seven output arguments are required, then Q and Z
%  are N by K orthonormal, S and T are K by K upper triangular such that
%  they form a partial generalized Schur decomposition: A*Q=Z*S and
%  B*Q=Z*T. Then LAMBDA=DIAG(S)./DIAG(T) and X=Q*Y with Y the eigenvectors
%  of the pair (S,T): S*Y=T*Y*JORDAN (see also OPTIONS.Schur).
%  
%  JDQZ(A,B) 
%  JDQZ('Afun','Bfun')
%  The first input argument is either a square matrix (which can be full
%  or sparse, symmetric or nonsymmetric, real or complex), or a string
%  containing the name of an M-file which applies a linear operator to the
%  columns of a given matrix. In the latter case, the M-file, say Afun.m,
%  must return the dimension N of the problem with N = Afun([],'dimension').
%  For example, JDQZ('fft',...) is much faster than JDQZ(F,...), where F is
%  the explicit FFT matrix.
%  If another input argument is a square N by N matrix or the name of an
%  M-file, then B is this argument (regardless whether A is an M-file or a
%  matrix). If B has not been specified, then B is assumed to be the
%  identity unless A is an M-file with two output vectors of dimension N
%  with [AV,BV]=Afun(V), or with AV=Afun(V,'A') and BV=Afun(V,'B').
%  
%  The remaining input arguments are optional and can be given in
%  practically any order:
%  
%  [X,JORDAN,Q,Z,S,T,HISTORY] = JDQZ(A,B,K,SIGMA,OPTIONS)
%  [X,JORDAN,Q,Z,S,T,HISTORY] = JDQZ('Afun','Bfun',K,SIGMA,OPTIONS)
%  
%  where
%  
%      K         an integer, the number of desired eigenvalues.
%      SIGMA     a scalar shift or a two letter string.
%      OPTIONS   a structure containing additional parameters.
%  
%  If K is not specified, then K = MIN(N,5) eigenvalues are computed.
%  
%  If SIGMA is not specified, then the Kth eigenvalues largest in
%  magnitude are computed. If SIGMA is a real or complex scalar, then the
%  Kth eigenvalues nearest SIGMA are computed. If SIGMA is column vector
%  of size (L,1), then the Jth eigenvalue nearest to SIGMA(MIN(J,L))
%  is computed for J=1:K. SIGMA is the "target" for the desired eigenvalues.
%  If SIGMA is one of the following strings, then it specifies the desired 
%  eigenvalues.
%  
%    SIGMA            Specified eigenvalues
%  
%    'LM'             Largest Magnitude  
%    'SM'             Smallest Magnitude (same as SIGMA = 0)
%    'LR'             Largest Real part
%    'SR'             Smallest Real part
%    'BE'             Both Ends. Computes K/2 eigenvalues
%                     from each end of the spectrum (one more
%                     from the high end if K is odd.)
%  
%  If 'TestSpace' is 'Harmonic' (see OPTIONS), then SIGMA = 0 is the
%  default, otherwise SIGMA = 'LM' is the default.
%  
%  
%  The OPTIONS structure specifies certain parameters in the algorithm.
%  
%   Field name            Parameter                             Default
%  
%   OPTIONS.Tol           Convergence tolerance:                1e-8 
%                           norm(r) <= Tol/SQRT(K)   
%   OPTIONS.jmin          Minimum dimension search subspace V   K+5
%   OPTIONS.jmax          Maximum dimension search subspace V   jmin+5
%   OPTIONS.MaxIt         Maximum number of iterations.         100
%   OPTIONS.v0            Starting space                        ones+0.1*rand
%   OPTIONS.Schur         Gives schur decomposition             'no'
%                           If 'yes', then X and JORDAN are
%                           not computed and [Q,Z,S,T,HISTORY]
%                           is the list of output arguments.
%   OPTIONS.TestSpace     Defines the test subspace W           'Harmonic'
%                           'Standard':    W=sigma*A*V+B*V
%                           'Harmonic':    W=A*V-sigma*B*V
%                           'SearchSpace': W=V
%                            W=V is justified if B is positive
%                            definite.
%   OPTIONS.Disp          Shows size of intermediate residuals  'no'
%                           and the convergence history
%   OPTIONS.NSigma        Take as target for the second and     'no'
%                           following eigenvalues, the best  
%                           approximate eigenvalues from the 
%                           test subspace.  
%   OPTIONS.Pairs         Search for conjugated eigenpairs      'no'
%   OPTIONS.LSolver       Linear solver                         'GMRES'
%   OPTIONS.LS_Tol        Residual reduction linear solver      1,0.7,0.7^2,..
%   OPTIONS.LS_MaxIt      Maximum number it.  linear solver     5
%   OPTIONS.LS_ell        ell for BiCGstab(ell)                 4
%   OPTIONS.Precond       Preconditioner  (see below)           identity.
%   OPTIONS.Type_Precond  Way of using preconditioner           'left'
%  
%  For instance
%  
%    options=struct('Tol',1.0e-8,'LSolver','BiCGstab','LS_ell',4,'Precond',M);
%  
%  changes the convergence tolerance to 1.0e-8, takes BiCGstab as linear 
%  solver, and takes M as preconditioner (for ways of defining M, see below).
%
%
%  PRECONDITIONING. The action M-inverse of the preconditioner M (an 
%  approximation of A-lamda*B) on an N-vector V can be defined in the 
%  OPTIONS
%  
%     OPTIONS.Precond
%     OPTIONS.L_Precond     same as OPTIONS.Precond
%     OPTIONS.U_Precond
%     OPTIONS.P_Precond
%
%  If no preconditioner has been specified (or is []), then M\V=V (M is
%  the identity).
%  If Precond is an N by N matrix, say, K, then
%        M\V = K\V.
%  If Precond is an N by 2*N matrix, say, K, then
%        M\V = U\L\V, where K=[L,U], and L and U are N by N matrices.
%  If Precond is a string, say, 'Mi', then
%        if Mi(V,'L') and Mi(V,'U') return N-vectors 
%               M\V = Mi(Mi(V,'L'),'U')
%        otherwise 
%               M\V = Mi(V) or M\V=Mi(V,'preconditioner').
%  Note that Precond and A can be the same string.
%  If L_Precond and U_Precond are strings, say, 'Li' and 'Ui', 
%  respectively, then
%        M\V=Ui(Li(V)).
%  If (P_precond,) L_Precond, and U_precond are N by N matrices, say, 
%  (P,) L, and U, respectively, then
%        M\V=U\L\(P*V)      (P*M=L*U)
%
%     OPTIONS.Type_Precond
%  The preconditioner can be used as explicit left preconditioner
%  ('left', default), as explicit right preconditioner ('right') or 
%  implicitly ('impl').
%  
%
%  JDQZ without input arguments returns the options and its defaults.
%

%   Gerard Sleijpen.
%   Copyright (c) 2002
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.3b.
% The toolbox can be obtained from http://www.cs.unimaas.nl/l.vandermaaten
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you 
% maintain the name of the original author.
%
% (C) Laurens van der Maaten
% Maastricht University, 2007

global Qschur Zschur Sschur Tschur ...
       Operator_MVs Precond_Solves ...
       MinvZ QastMinvZ

if nargin==0
   possibilities, return,
end

%%% Read/set parameters
[n,nselect,Sigma,kappa,SCHUR,...
   jmin,jmax,tol0,maxit,V,AV,BV,TS,DISP,PAIRS,JDV0,FIX_tol,track,NSIGMA,...
   lsolver,LSpar] = ReadOptions(varargin{1:nargin});

Qschur = zeros(n,0);    Zschur=zeros(n,0);; 
MinvZ  = zeros(n,0);    QastMinvZ=zeros(0,0); 
Sschur = []; Tschur=[]; history = []; 

%%% Return if eigenvalueproblem is trivial
if n<2
  if n==1, Qschur=1; Zschur=1; [Sschur,Tschur]=MV(1); end
  if nargout == 0, Lambda=Sschur/Tschur, else
  [varargout{1:nargout}]=output(history,SCHUR,1,Sschur/Tschur); end, 
return, end

%---------- SET PARAMETERS & STRINGS FOR OUTPUT -------------------------

if     TS==0, testspace='sigma(1)''*Av+sigma(2)''*Bv';
elseif TS==1, testspace='sigma(2)*Av-sigma(1)*Bv';
elseif TS==2, testspace='v'; 
elseif TS==3, testspace='Bv';
elseif TS==4, testspace='Av';
end

String=['\r#it=%i #MV=%3i, dim(V)=%2i, |r_%2i|=%6.1e  '];

%------------------- JDQZ -----------------------------------------------

% fprintf('Scaling with kappa=%6.4g.',kappa)

k=0; nt=0; j=size(V,2); nSigma=size(Sigma,1);
it=0; extra=0; Zero=[]; target=[]; tol=tol0/sqrt(nselect);

INITIATE=1;  JDV=0; 
rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0;

time=clock;
if TS ~=2
while (k<nselect & it<maxit)

   %%% Initialize target, test space and interaction matrices
   if INITIATE, % set new target
      nt=min(nt+1,nSigma); sigma = Sigma(nt,:); nlit=0; lit=0;  
      if j<2
        [V,AV,BV]=Arnoldi(V,AV,BV,sigma,jmin,nselect,tol);
        rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0; target=[];
        j=min(jmin,n-k);
      end
      if DETECTED & NSIGMA
         [Ur,Ul,St,Tt] = SortQZ(WAV,WBV,sigma,kappa);
         y=Ur(:,1); q=V*y; Av=AV*y; Bv=BV*y; 
         [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa);
         sigma=ScaleEig(theta);
         USE_OLD=NSIGMA; rKNOWN=1; lit=10;
      end   
      NEWSHIFT= 1; 
      if DETECTED & TS<2, NEWSHIFT= ~min(target==sigma); end
      target=sigma; ttarget=sigma;
      if ischar(ttarget), ttrack=0; else, ttrack=track; end
      if NEWSHIFT 
         v=V; Av=AV; Bv=BV; W=eval(testspace);
         %%% V=RepGS(Qschur,V); [AV,BV]=MV(V); %%% more stability??
         %%% W=RepGS(Zschur,eval(testspace));  %%% dangerous if sigma~lambda
         if USE_OLD, W(:,1)=V(:,1); end, 
         W=RepGS(Zschur,W); WAV=W'*AV;  WBV=W'*BV;
      end
      INITIATE=0; DETECTED=0; JDV=0;

   end % if INITIATE

   %%% Solve the preconditioned correction equation
   if rKNOWN,
      if JDV, z=W; q=V; extra=extra+1; 
         if DISP,  fprintf('  %2i-d proj.\n',k+j-1), end 
      end
      if FIX_tol*nr>1 & ~ischar(target), theta=target; else, FIX_tol=0; end
      t=SolvePCE(theta,q,z,r,lsolver,LSpar,lit); 
      nlit=nlit+1; lit=lit+1; it=it+1;
      EXPAND=1; rKNOWN=0; JDV=0;
   end % if rKNOWN
    
   %%% Expand the subspaces and the interaction matrices
   if EXPAND
      [v,zeta]=RepGS([Qschur,V],t);
      V=[V,v]; 
      [Av,Bv]=MV(v); AV=[AV,Av]; BV=[BV,Bv]; 
      w=eval(testspace); w=RepGS([Zschur,W],w);
      WAV=[WAV,W'*Av;w'*AV]; WBV=[WBV,W'*Bv;w'*BV]; W=[W,w];
      j=j+1; EXPAND=0;

      %%% Check for stagnation
      if abs(zeta(size(zeta,1),1))/norm(zeta)<0.06, JDV=JDV0; end

   end % if EXPAND
 
   %%% Solve projected eigenproblem
   if USE_OLD
      [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin,y); 
   else
      [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin); 
   end

   %%% Compute approximate eigenpair and residual
   y=Ur(:,1); q=V*y; Av=AV*y; Bv=BV*y; 
   [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa); 
   %%%=== an alternative, but less stable way of computing z =====
   % beta=Tt(1,1); alpha=St(1,1); theta=[alpha,beta];
   % r=RepGS(Zschur,beta*Av-alpha*Bv,0); nr=norm(r); z=W*Ul(:,1);
   rKNOWN=1; if nr<ttrack, ttarget=ScaleEig(theta); end

         if DISP,                                  %%% display history
            fprintf(String,it,Operator_MVs,j,nlit,nr), 
         end 
         history=[history;nr,it,Operator_MVs];    %%% save history

   %%% check convergence 
   if (nr<tol)
      %%% EXPAND Schur form
      Qschur=[Qschur,q]; Zschur=[Zschur,z];
      Sschur=[[Sschur;zeros(1,k)],Zschur'*Av]; 
      Tschur=[[Tschur;zeros(1,k)],Zschur'*Bv];  Zero=[Zero,0];
      k=k+1; 
      if ischar(target), Target(k,:)=[nt,0,0];
      else, Target(k,:)=[0,target]; end
      if DISP, ShowEig(theta,target,k); end
      if (k>=nselect), break; end;
      %%% Expand preconditioned Schur matrix MinvZ=M\Zschur
      UpdateMinvZ;
      J=[2:j]; j=j-1; Ur=Ur(:,J); Ul=Ul(:,J); 
      V=V*Ur; AV=AV*Ur; BV=BV*Ur; W=W*Ul; 
      WAV=St(J,J); WBV=Tt(J,J);
  
      rKNOWN=0; DETECTED=1;  USE_OLD=0;

      %%% check for conjugate pair
      if PAIRS & (abs(imag(theta(1)/theta(2)))>tol) 
         t=ImagVector(q); % t=conj(q); t=t-q*(q'*t);
         if norm(t)>tol, t=RepGS([Qschur,V],t,0); 
            if norm(t)>200*tol
               target=ScaleEig(conj(theta));
               EXPAND=1; DETECTED=0; 
               if DISP, fprintf('--- Checking for conjugate pair ---\n'), end
            end
         end
      end
    
      INITIATE = ( j==0 & DETECTED);

   elseif DETECTED %%% To detect whether another eigenpair is accurate enough 
      INITIATE=1; 
   end % if (nr<tol)
   
   %%% restart if dim(V)> jmax
   if j==jmax
      j=jmin; J=[1:j]; 
      Ur=Ur(:,J); Ul=Ul(:,J); 
      V=V*Ur; AV=AV*Ur; BV=BV*Ur; W=W*Ul; 
      WAV=St(J,J); WBV=Tt(J,J); 
   end % if j==jmax

end % while k
end % if TS~=2


if TS==2
Q0=Qschur; ZastQ=[];
% WAV=V'*AV; WBV=V'*BV;
while (k<nselect & it<maxit)

   %%% Initialize target, test space and interaction matrices
   if INITIATE & ( nSigma>k | NSIGMA), % set new target
      nt=min(nt+1,nSigma); sigma = Sigma(nt,:); nlit=0; lit=0;                
      if j<2
        [V,AV,BV]=Arnoldi(V,AV,BV,sigma,jmin,nselect,tol);
        rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0; target=[];
        j=min(jmin,n-k);; 
      end
      if DETECTED & NSIGMA
         [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,sigma,kappa,1); 
         q=RepGS(Zschur,V*Ur(:,1)); [Av,Bv]=MV(q); 
         [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa); 
         sigma=ScaleEig(theta); 
         USE_OLD=NSIGMA; rKNOWN=1; lit=10;
      end   
      target=sigma; ttarget=sigma;
      if ischar(ttarget), ttrack=0; else, ttrack=track; end
      if ~DETECTED
         %%% additional stabilisation. May not be needed
         %%% V=RepGS(Zschur,V); [AV,BV]=MV(V); 
         %%% end add. stab.
         WAV=V'*AV; WBV=V'*BV;
      end
      DETECTED=0; INITIATE=0; JDV=0; 
   end % if INITIATE

   %%% Solve the preconditioned correction equation
   if rKNOWN,
      if JDV, z=V; q=V; extra=extra+1; 
         if DISP,  fprintf('  %2i-d proj.\n',k+j-1), end 
      end
      if FIX_tol*nr>1 & ~ischar(target), theta=target; else, FIX_tol=0; end
      t=SolvePCE(theta,q,z,r,lsolver,LSpar,lit); 
      nlit=nlit+1; lit=lit+1; it=it+1;
      EXPAND=1; rKNOWN=0; JDV=0;
   end % if rKNOWN

   %%% expand the subspaces and the interaction matrices
   if EXPAND
      [v,zeta]=RepGS([Zschur,V],t); [Av,Bv]=MV(v); 
      WAV=[WAV,V'*Av;v'*AV,v'*Av]; WBV=[WBV,V'*Bv;v'*BV,v'*Bv];
      V=[V,v]; AV=[AV,Av]; BV=[BV,Bv]; 
      j=j+1; EXPAND=0;

      %%% Check for stagnation
      if abs(zeta(size(zeta,1),1))/norm(zeta)<0.06, JDV=JDV0; end 

   end % if EXPAND
 
   %%% compute approximate eigenpair
   if USE_OLD
      [Ur,Ul]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin,Ur(:,1)); 
   else
      [Ur,Ul]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin);
   end
     
   %%% Compute approximate eigenpair and residual
   q=V*Ur(:,1); Av=AV*Ur(:,1); Bv=BV*Ur(:,1);  
   [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa); 
   rKNOWN=1; if nr<ttrack, ttarget=ScaleEig(theta); end

          if DISP,                                 %%% display history
             fprintf(String,it,Operator_MVs, j,nlit,nr), 
          end  
          history=[history;nr,it,Operator_MVs];   %%% save history
   
   %%% check convergence 
   if (nr<tol)
      %%% expand Schur form
      [q,a]=RepGS(Q0,q); a1=a(k+1,1); a=a(1:k,1);
      %%% ZastQ=Z'*Q0  
      Q0=[Q0,q]; %%% the final Qschur
      ZastQ=[ZastQ,Zschur'*q;z'*Q0]; Zschur=[Zschur,z]; Qschur=[Qschur,z]; 
      Sschur=[[Sschur;Zero],a1\(Zschur'*Av-[Sschur*a;0])];