greed_nomp.m

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Main algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
   display('Main iterations...') 
end
tic
t=0;
p=zeros(m,1);
DR=Pt(Residual);
[v I] = max(abs(DR));
if weakness ~= 1
    [vals inds] = sort(abs(DR),'descend');
    I=inds( find( vals >= alpha * v ) );
end
    
IN = union(IN,I);
if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
    IN=IN(1:STOPTOL);
end
MASK=zeros(size(DR));
pDDp=1;
done = 0;
iter=1;

while ~done

    
    % Select new element
    if isa(GradSteps,'char')
        if strcmp(GradSteps,'auto')
             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Iteration to automatic selection of the number of gradient steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%             finished=0;    
%             while ~finished
            % Update direction    
                 if iter==1
                     p(IN)=DR(IN);
                     Dp=P(p);
                 else
                     MASK(IN)=1;
                     PDR=P(DR.*MASK);
                     b=-Dp'*PDR/pDDp;
                     p(IN)=DR(IN) +b*p(IN);
                     Dp=PDR +b* Dp;
                 end
             % Step size
%                  Dp=P(p); % =P(DR(IN)) +b P(p(IN));
                 pDDp=Dp'*Dp;
                 a=Residual'*Dp/(pDDp);
             % Update coefficients   
                 s=s+a*p;
             % New Residual and inner products
                 Residual=Residual-a*Dp;
                 DR=Pt(Residual);
                 % select new element
                     [v I] = max(abs(DR));
                     if weakness ~= 1
                         [vals inds] = sort(abs(DR),'descend');
                         I=inds( find( vals >= alpha * v ) );
                     end
                     IN = union(IN,I);
                     if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
                        IN=IN(1:STOPTOL);
                     end

%                  % Only if we select new element do we leave the loop    
%                      if isempty(find (IN==I, 1))
%                         IN=[IN I];
%                         finished=1;
%                      end
%             end
        else
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Is option known?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            error('Undefined option for GradSteps, use ''auto'' or an integer.')
        end
    elseif isa(GradSteps,'numeric') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               Iteration for fixed number of gradient steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
                
        % Do GradSteps gradient steps
        count=1;
        while count<=GradSteps
            % Update direction    
                 if iter==1
                     p(IN)=DR(IN);
                     Dp=P(p);
                 else
                     MASK(IN)=1;
                     PDR=P(DR.*MASK);
                     b=-Dp'*PDR/pDDp;
                     p(IN)=DR(IN) +b*p(IN);
                     Dp=PDR +b* Dp;
                 end
             % Step size
%                  Dp=P(p);   
                 pDDp=Dp'*Dp;
                 a=Residual'*Dp/(pDDp);
             % Update coefficients   
                 s=s+a*p;
             % New Residual and inner products
                 Residual=Residual-a*Dp;
                 DR=Pt(Residual);
                  count=count+1;
        end
             % select new element
                 [v I] = max(abs(DR));
                 if weakness ~= 1
                     [vals inds] = sort(abs(DR),'descend');
                     I=inds( find( vals >= alpha * v ) );
                 end
                 IN = union(IN,I);
                 if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
                    IN=IN(1:STOPTOL);
                 end
                
     else
          error('Undefined option for GradSteps, use ''auto'' or an integer.')
     end


     ERR=Residual'*Residual/n;
     if comp_err
         err_mse(iter)=ERR;
     end
     
     if comp_time
         iter_time(iter)=toc;
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Are we done yet?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if strcmp(STOPCRIT,'M')
         if length(IN) >= STOPTOL
             done =1;
         elseif verbose && toc-t>10
            display(sprintf('Iteration %i. --- %i selected elements',iter ,length(IN))) 
            t=toc;
         end
    elseif strcmp(STOPCRIT,'mse')
         if comp_err
            if err_mse(iter)<STOPTOL;
                done = 1; 
            elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter))) 
                t=toc;
            end
         else
             if ERR<STOPTOL;
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,ERR)) 
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'mse_change') && iter >=2
         if comp_err && iter >=2
              if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) 
                t=toc;
             end
         else
             if ((oldERR - ERR)/sigsize < STOPTOL);
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize)) 
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'corr') 
          if max(abs(DR)) < STOPTOL;
             done = 1; 
          elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR)))) 
                t=toc;
          end
     end
     
    % Also stop if residual gets too small or maxIter reached
     if comp_err
         if err_mse(iter)<1e-16
             display('Stopping. Exact signal representation found!')
             done=1;
         end
     else


         if iter>1
             if ERR<1e-16
                 display('Stopping. Exact signal representation found!')
                 done=1;
             end
         end
     end

     if iter >= MAXITER
         display('Stopping. Maximum number of iterations reached!')
         done = 1; 
     end
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    If not done, take another round
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
     if ~done
        iter=iter+1;
        oldERR=ERR;
     end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  Only return as many elements as iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout >=2
    err_mse = err_mse(1:iter);
end
if nargout ==3
    iter_time = iter_time(1:iter);
end
if verbose
   display('Done') 
end


%
% Change history
% 9th of May: Set update direction to gradient for first iteration. 
%             Otherwise incorrect direction is used when initial vector is given. 
%
% 22 of November: Original code did not calculate b value correctly. This
%                 has now been fixed.
%
% 14 of December: Faster implementation using recursive calculation of Dp
%
%  4 of January: Implemented weak element selection. This allows the
%  algorithm to select more than a single element in each iteration.
%
% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.
%
% 29 of Februray: USed unon command to combine selected elements.