eigifp.m

Frequently used algorithms for numerical analysis and signal processing

function [EVals, EVecs] = eigifp(varargin) 
%
%
%   EIGIFP (version 2.1.1):  computes a few (algebraically) smallest or largest 
%       eigenvalues/eigenvectors of the symmetric matrix eigenvalue problem:
%
%             A x = lambda x   or    A x = lambda B x
%
%       where A and B are symmetric and B is positive definite.
%
%   [d, x] = eigifp(A)      => smallest eigenvalue/eigenvector of A;
%   [d, X] = eigifp(A,k)    => k smallest eigenvalues/eigenvectors of A;
%
%   [d, x] = eigifp(A,B)    => smallest eigenvalue/eigenvector of (A,B);
%   [d, X] = eigifp(A,B,k)  => k smallest eigenvalue/eigenvector of (A,B);
%
%   A (and/or B) can be either a sparse matrix or a function that 
%   returns A*x for x. In the latter case, the dimension of A has 
%   to be input through options.size (see below). It is also desirable  
%   to have an estimate of a norm of A (or B) inputed through 
%   options.NORMA (or options.NORMB). 
%
%   To compute the largest eigenvalues, one can call as above by adding 
%   the option options.MAXEIG = 1 (see below). Alternatively, if A is 
%   numeric, one can do the following.
%
%   [d, x] = eigifp(-A); d=-d;     => largest eigenvalue/eigenvector of A;
%   [d, X] = eigifp(-A,k); d=-d;   => k largest eigenvalues/eigenvectors of A;
%
%   [d, x] = eigifp(-A,B); d=-d;   => largest eigenvalue/eigenvector of (A,B);
%   [d, X] = eigifp(-A,B,k); d=-d; => k largest eigenvalues/eigenvectors of (A,B);
%
%   Other optional inputs can be provided to use any a priori information,
%   to control computational cost, and to optimize performance.
%  
%   eigifp(A,opt), eigifp(A,k opt), eigifp(A,B,opt), eigifp(A,B,k,opt)
%   take the following optional inputs defined by opt:  
%
%     opt.SIZE:     
%       the size of A. 
%     opt.NORMA:    
%       an estimate of the 1-norm of A 
%     opt.NORMB:    
%       an estimate of the 1-norm of B
%     opt.INITIALVEC or opt.V0: 
%       a matrix whose i-th column is the i-th initial eigenvector; 
%     opt.TOLERANCE or opt.TOL: 
%       termination tolerance for the 1-norm of residual; 
%       default: 10*eps*sqrt(n)*(||A||+|lambda|*||B||)
%     opt.MAXITERATION or opt.MAXIT: 
%       set the maximum number of (outer) iterations; default: 500;
%     opt.INNERITERATION or opt.INNERIT: 
%       set a fixed inner iteration to control the memory requirement; 
%       Default: between 1 and 128 as adaptively determined.   
%     opt.USEPRECON:    
%       set preconditioning off or completely on:  
%           1.  set opt.USEPRECON='NO' to switch off preconditioning; 
%           2.  set opt.USEPRECON = an approximate eigenvalue 
%               to start preconditioned iterations with a preconditioner 
%           computed using the eigenvalue as the shift;              
%     opt.ILUTHRESH or opt.ILU:     
%       threshold between 0 and 1 used in incomplete LU for computing  
%           preconditioner (i.e. used in luinc.m); default: 1e-4; 
%       Special cases: 0 = exact LU; 1 = 0-level ILU;  
%     opt.PRECONDITIONER or opt.PRECON:  
%       user supplied preconditioner. It can be a matrix or a function.
%     opt.DISP:
%       Set to 0 to disable on-screen display of output, and to any other 
%       numerical value to enable display. Default: 1
%     opt.MAXEIG:
%       Set to 1 to compute the largest eigenvalues instead of the 
%       smallest eigenvalue of (A,B). Any other value causes the code to 
%       compute the smallest eigenvalue. Default: 0
%
%   
%   EIGIFP is based on the algorithm in: G. Golub and Q. Ye, An Inverse 
%   Free Preconditioned Krylov Subspace Method for Symmetric Generalized
%       Eigenvalue Problems,  SIAM J. Sci. Comput., 24(2002):312-334. It also
%       incorporates several enhancements. 
%   
%       It is particularly suitable for problems where 
%       a)  factorization of B is difficult; or
%       b)  factorization of A-lambda_k B is difficult; or 
%       c)  a shift near eigenvalues sought is not known. 
%
%
%   EIGIFP is developed by Qiang Ye ([e-mail]) with the programming 
%   assistance provided by James Money. This work was supported by NSF  
%   under Grant CCR-0098133. 
%   
%   This program is provided for research and educational use and  
%   is distributed through http://www.ms.uky.edu/~qye/software. 
%   Neither redistribution nor commercial use is permitted without 
%   consent of the author. 
%   
%   Copyright (c): Qiang Ye
%   Version 2.1.1 is dated Oct. 2004
%   Version 2.1 is dated July 2004
%   Version 2.0 is dated Aug. 2003
%   Version 1 is dated Sept. 2002
%

                    % Process inputs 
A=varargin{1};
n=[]; 
if ( isnumeric(A) )
    n=size(A,1);
    if (any(size(A) ~= n))
        error('Input Error: Matrix A is not square');
    end
elseif ( ~isstr(A) ) 
    error('Input Error: A must be either a matrix or a function');
end 

B=[];
if ( nargin>1 )
    if ( (~ isstruct(varargin{2})) ) 
        if( isstr(varargin{2}) | (any(size(varargin{2}) ~= 1)))
        B=varargin{2};
        if ( isnumeric(B) )
            if ( isempty(n) )
                n=size(B,1);
            end 
            if (any(size(B) ~= n))
                    error('Input Error: Matrix B must be square and of the size of A');
            end
        elseif ( ~isstr(B) ) 
            error('Input Error: B must be either a matrix or a function');
        end         
        end
    end
end

                    % set initial or default values 
m=0;
k=1;
X=[];
L=[];
eta=1e-4;
usePrecon=1;
iterMax=500;
ksetbyuser=0;
tolerance=[];
normA=[]; 
normB=[]; 
outputYes=1;
EVals=[];
EVecs=[];
findMax=0;

if ( nargin > 1 + ~isempty(B) )
    if ( (~ isstruct(varargin{2+~isempty(B)})) & all(size(varargin{2+~isempty(B)}) == 1))
        k=varargin{2+~isempty(B)};
        ksetbyuser=1;
    end
end
if ( nargin > 1 + ~isempty(B) + ksetbyuser )
    options=varargin{2 + ~isempty(B) + ksetbyuser:nargin};
    if ( ~isstruct(options) ) 
        error('Input Error: too many input parameters.'); 
    end
    names=fieldnames(options);
    
    I=strmatch('DISP',upper(names),'exact');
    if (~isempty(I))
        outputYesStr=getfield(options,names{I});
        if (isnumeric(outputYesStr))
            if (outputYesStr==0)
                outputYes=0;
            else
                outputYes=1;
            end
        end
    end
    I=strmatch('SIZE',upper(names),'exact');
    if (~isempty(I))
        n=getfield(options,names{I});
        if( isnumeric(A) & n ~= size(A, 1) )
            error('Input Error: options.size not equal to the size of A');
        end 
        if( ~isempty(B) & isnumeric(B) & n ~= size(B, 1) )
            error('Input Error: options.size not equal to the size of B');
        end 
    end
    
    I=strmatch('NORMA',upper(names),'exact');
    if (~isempty(I))
        normA=getfield(options,names{I});
    end

    I=strmatch('NORMB',upper(names),'exact');
    if (~isempty(I))
        normB=getfield(options,names{I});
    end

    I=strmatch('INITIALVEC',upper(names),'exact');
    if (isempty(I))
        I=strmatch('V0',upper(names),'exact');
    end
    if (~isempty(I))
        X=getfield(options,names{I});
        if( isempty(n) )
            n=size(X, 1); 
        elseif ( size(X, 1) ~=n )
            error('Input Error: incorrect size of initial vectors');
        end 
        if( size(X, 2) ~=k )
            fprintf(1,'Warning: Number of initial vectors not equal to k\n');
        end 
    end
    
    I=strmatch('TOLERANCE',upper(names),'exact');
    if (isempty(I))
        I=strmatch('TOL',upper(names),'exact');
    end
    if (~isempty(I))
        tolerance=getfield(options,names{I});
        if ( tolerance <= 0 )
                error('Input Error: Invalid tolerance input.');
            end 
    end
    
    if ( isempty(n) )
        error('Input Error: need to input options.size = (dimension of A)');
    end
     
    I=strmatch('INNERITERATION',upper(names),'exact');
    if (isempty(I))
        I=strmatch('INNERIT',upper(names),'exact');
    end
    if (~isempty(I))
        m=getfield(options,names{I});
        if ( m < 0 | m >= n )
                error('Input Error: Invalid inner iteration number');
            end     
    end
    
    I=strmatch('ILUTHRESH',upper(names),'exact');
    if (isempty(I))
        I=strmatch('ILU',upper(names),'exact');
    end 
    if (~isempty(I))
        eta=getfield(options,names{I});
        if ( eta < 0 )
                error('Input Error: Invalid incomplete factorization threshold.');
            end             
    end
    
    I=strmatch('PRECONDITIONER',upper(names),'exact');
    if (isempty(I))
        I=strmatch('PRECON',upper(names),'exact');
    end
    if (~isempty(I))
        L=getfield(options,names{I});
        if ( isnumeric(L) & any(size(L) ~= n)) 
                error('Input Error: Invalid preconditioner; must be square.');
            end             
    end
    
    I=strmatch('MAXEIG',upper(names),'exact');
    if (~isempty(I))
        findMax=getfield(options,names{I});
        if ( ~isnumeric(findMax) | ~((findMax==0) | (findMax==1)))
            fprintf(1,'Warning: opt.MAXEIG should be either 0 or 1. Reset to 0.\n');
        end             
    end
    
    I=strmatch('USEPRECON',upper(names),'exact');
    if (~isempty(I))
        usePreconIn=getfield(options,names{I});
        if ( ~isempty(usePreconIn) &  ~isempty(L) )
            error('Input Error: input opt.useprecon while a preconditioner has been provided.'); 
        end
        if (isstr(usePreconIn))
            if (strcmp(upper(usePreconIn),'NO'))
                usePrecon=0;                    
            end
        elseif (~isempty(usePreconIn))      
            if ( any(size(usePreconIn) > 1) ) 
                    error('Input Error: Invalid shift; must be a real number.');
                end         
                if ( isnumeric(A) & isnumeric(B) )
                if(~issparse(A)) 
                    A=sparse(A);
                    fprintf(1,'Warning: A is not in the sparse format.\n');             
                end  
                    if (outputYes ~=0) 
                    fprintf(1,'Computing threshold ILU factorization using the shift provided.\n');
                end
                if ( ~isempty(B) )
                    if(~issparse(B)) 
                        B=sparse(B);
                        fprintf(1,'Warning: B is not in the sparse format.\n');             
                    end  
                    L=ildlte(A-usePreconIn*B,eta);
                else
                    n=size(A,1); 
                    L=ildlte(A-usePreconIn*speye(n),eta);
                end
            else
                fprintf(1,'Warning: A (or B) is a function; can not compute preconditioner using opt.useprecon.\n');            
            end
        end
    end
    I=strmatch('MAXITERATION',upper(names),'exact');
    if (isempty(I))
        I=strmatch('MAXIT',upper(names),'exact');
    end
    if (~isempty(I))
        iterMax=getfield(options,names{I});     
        if (  iterMax <= 0 )
                error('Input Error: Invalid maximum iteration.');
            end             

    end
end

if ( isempty(normA) )
    if ( isnumeric(A) )
        normA=norm(A,1); 
    else
        normA=norm(feval(A,ones(n,1)), 1)/n; 
        fprintf(1,'Warning: options.Anorm (estimate of ||A||_1) is not provided.\n');   
    end
end 
if ( isempty(normB) )
    if ( isempty(B) )
        normB=1;
    elseif ( isstr(B) )
        normB=norm(feval(B,ones(n,1)), 1)/n; 
        fprintf(1,'Warning: no options.Bnorm (estimate of ||B||_1) is not provided.\n');    
    elseif ( isnumeric(B) )
        normB=norm(B,1);
    end
end 

                        % set default tolerance 
toleranceA=10*eps*normA*sqrt(n);
if(~isempty(B)) 
   toleranceB=10*eps*normB*sqrt(n);
else 
   toleranceB=0; 
end 
if ( ~isempty(tolerance) )
    toleranceVec=[tolerance, 0];
    if ( tolerance < toleranceA )  
        fprintf(1,'Warning: Tolerance may be set too low. Suggested value: %E.\n',toleranceA);
    end 
else
    toleranceVec=[toleranceA,toleranceB];
end

if ( isstr(A) | isstr(B) ) 
    usePrecon=0;
end

if ( isempty(L) & (usePrecon ~= 0) ) 
    if(~issparse(A)) 
        A=sparse(A);
        fprintf(1,'Warning: A is not in the sparse format.\n');             
    end  
    if ( ~isempty(B) & ~issparse(B)) 
        B=sparse(B);
        fprintf(1,'Warning: B is not in the sparse format.\n');             
    end
end  
                    
if (k>n)
    error('Error: # of eigenvalues sought is greater than the matrix size');
end

                        % call the main function ifree.m

if (isnumeric(A))       %if given a matrix, use (A,B) or (-A,B) if maxeig =1
    if (findMax==1)
        if (outputYes)
            fprintf(1,'Computing the smallest eigenvalues of (-A,B) first.\n');
        end
        [EVals, EVecs] = ifree(-A,B,n,m,toleranceVec,iterMax,k,X,eta,L,usePrecon,normA,normB,outputYes);
        if (outputYes)
            fprintf(1,'  Negating the eigenvalues of (-A, B) => the largest eigenvalues of (A,B):\n');
            EVals=-EVals;
            fprintf(1,'  %e\n',EVals);
        end        
    else
        [EVals, EVecs] = ifree(A,B,n,m,toleranceVec,iterMax,k,X,eta,L,usePrecon,normA,normB,outputYes);
    end
else                    % if A is a function 
    if (findMax==1)     % here we use minusA and negate the results when done
        if (outputYes)
            fprintf(1,'Computing the smallest eigenvalues of (-A,B) first.\n');
        end
        saveAfunc=0;
        if (exist('Afunc123')==1)   % we don't want to overwrite the user's global, so we save it
            saveAfunc=1;
            tempFunc=Afunc;
            clear Afunc123;
        end
        global Afunc123;
        Afunc123=A;     % this is the function it calls in minusA()
        [EVals, EVecs] = ifree('minusA',B,n,m,toleranceVec,iterMax,k,X,eta,L,usePrecon,normA,normB,outputYes);
        if (saveAfunc==1)
            Afunc123=saveAfunc; %restore the value if the user already defined it
        end
        if (outputYes)
            fprintf(1,'  Negating the eigenvalues of (-A, B) => the largest eigenvalues of (A,B):\n');
            EVals=-EVals;
            fprintf(1,'  %e\n',EVals);
        end      
    else
        [EVals, EVecs] = ifree(A,B,n,m,toleranceVec,iterMax,k,X,eta,L,usePrecon,normA,normB,outputYes);
    end
end





function [lambda, x, r, bx, diff, rDiff, bDiff, lam_max, res] = ...
     parnob(A, B, L, lambda, x, r, bx, diff, rDiff, bDiff, m, tol, EVecs, BEVecs, lambda_1, lambda_max)
%  
%   generate B-othonormal basis V of preconditioned (by L) Krylov subspace
%   by m steps of Arnoldi algorithm and then compute the Ritz value  
%   and Ritz vectors 
%   
%   BEVecs = B*EVecs  with EVecs =  converged eigenvectors
%

                    % initialization and normalization 
n=size(x,1);                     
V = zeros(n,m);
Wr=V; 
Am = zeros(m,m);
temp=bx'*x;
if (temp <= 0)
    if ( isnumeric(B) )     % double check on x'Bx
        bx=B*x;
    else
        bx=feval(B, x);
    end
    temp=bx'*x;
    if (temp <= 0)
        error('Error: B is not positive definite.');
    end
end
temp=sqrt(temp);
V(:,1) = x / temp;
r=r/temp; 
Wr(:,1) = r;
bx=bx/temp;
if ( ~ isempty(B) )
        Wb=V; 
    Wb(:,1) = bx;
end
Am(1,1)=0;

                    % Loop for Arnoldi iteration
for i = 2:m-1,                    
                                    % Apply preconditioner if given 
    if ( isnumeric(L) )         
        r=L\r;
        r=(L')\r;    
    else 
        r=feval(L,r); 
    end
                    % generate new basis vector 
    for k = 1:(i-1)
        if ( ~ isempty(B) )