greed_omp.m

Sparse Estimation Algorithms by Blumensath and Davies min ||x||_0 subject to ||y - Ax||_2<e

            else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end 
        case {'stopTol'}
            if isa(Options{i+1},'numeric') ; STOPTOL     = Options{i+1};   
            else error('stopTol must be number. Exiting.'); end
        case {'P_trans'} 
            if isa(Options{i+1},'function_handle'); Pt = Options{i+1};   
            else error('P_trans must be function _handle. Exiting.'); end
        case {'maxIter'}
            if isa(Options{i+1},'numeric'); MAXITER     = Options{i+1};             
            else error('maxIter must be a number. Exiting.'); end
        case {'verbose'}
            if isa(Options{i+1},'logical'); verbose     = Options{i+1};   
            else error('verbose must be a logical. Exiting.'); end 
        case {'start_val'}
            if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ;
                s_initial     = Options{i+1};   
            else error('start_val must be a vector of length m. Exiting.'); end
         otherwise
            error('Unrecognised option. Exiting.') 
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Make P and Pt functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if          isa(A,'float')      P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isobject(A)         P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isa(A,'function_handle') 
    try
        if          isa(Pt,'function_handle'); P=A;
        else        error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end
    catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end
else        error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Random Check to see if dictionary is normalised 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        mask=zeros(m,1);
        mask(ceil(rand*m))=1;
        nP=norm(P(mask));
        if abs(1-nP)>1e-3;
            display('Dictionary appears not to have unit norm columns.')
        end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Select algorithm to use
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~any(strcmp(SOLVER,{'auto','qr','chol','cgp','cg','pinv','linsolve'})) 
    display('Unknown solver specified, using automatic selection')
    SOLVER = 'auto';
end
if strcmp(SOLVER,'auto')
    % determine problem size and fastest algo to use
    if strcmp(STOPCRIT,'M') 
        maxM=STOPTOL;
    else
        maxM=MAXITER;
    end
        
	if maxM*n < 1e7 && isa(A,'float')
        try TESTMAT1=zeros(n,maxM);TESTMAT2=zeros(maxM); 
        catch
            try TESTMAT2=zeros(maxM);
            catch
                SOLVER='cg';
                display('Memory requirements too large. Using full conjugate solver in each iteration.')
            end
            clear TESTMAT2
            SOLVER='chol';
            display('Memory requirements too large for QR. Using Cholesky solver.')
        end
            clear TESTMAT1 TESTMAT2
            SOLVER = 'qr';
            display ('Memory requirement acceptable. Using QR method')
	elseif maxM^2 < 1e7
        try TESTMAT2=zeros(maxM);
        catch
            SOLVER='cg';
            display('Memory requirements too large. Using full conjugate solver in each iteration. Reduce maxIter to use other solver.')
        end
        clear TESTMAT2
        SOLVER='chol';
        display('Memory requirements possibly too large for QR. Using Cholesky solver. Reduce maxIter to use QR.')
    else
        SOLVER='cg';
        display('Memory requirements too large. Using full conjugate solver in each iteration.')
    end

end
if strcmp(SOLVER,'qr')
    display('Trying to use QR OMP implementation.')
    try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_qr(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_qr(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_qr(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Error using QR algorithm. Trying Cholesky instead.')
         SOLVER = 'chol';
     end
end
if strcmp(SOLVER,'cgp') 
    display('Trying to use CGP implementation.')
     try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_cgp(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_cgp(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_cgp(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Error using CG algorithm. Trying Cholesky instead.')
         SOLVER = 'chol';
     end
end
if strcmp(SOLVER,'chol') 
    display('Trying to use Cholesky OMP implementation.')
     try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_chol(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_chol(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_chol(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Error using Cholesky algorithm. Problem instance probably too large. We recommend the use of greed_gp or greed_acgp.')
     end
end   
if strcmp(SOLVER,'cg') 
    display('Trying to use Full CG OMP implementation. This can take a while....')
     try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_cg(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_cg(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_cg(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Something wrong. Full CG did not work.')
     end
end 
if strcmp(SOLVER,'pinv') 
    display('Trying to use pinv OMP implementation. This can take a while....')
     try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_pinv(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_pinv(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_pinv(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Something wrong. pinv did not work.')
     end
end 
if strcmp(SOLVER,'linsolve') 
    display('Trying to use linsolve OMP implementation. This can take a while....')
     try 
        if comp_err
            if comp_time
                [s, err_norm, iter_time]    = greed_omp_linsolve(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});
            else 
                [s, err_norm]               = greed_omp_linsolve(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});  
            end
        else
                [s]                         = greed_omp_linsolve(x,P,m,{'start_val',s_initial,'stopCrit',STOPCRIT,'stopTol',STOPTOL,'P_trans',Pt,'maxIter',MAXITER,'verbose',verbose});           
        end
     catch 
         display('Something wrong. linsolve did not work.')
     end
end 

if ~exist('s','var')
    error('Something is wrong in the code. We should have an answer by now. Exiting.')
end
    
% Change history
%
% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.
%                Changed automatic selection of QR and Cholesky method.