greed_omp.m

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

function [s, err_norm, iter_time]=greed_omp(x,A,m,varargin)
% greed_omp: Orthogonal Matching Pursuit using a range of implementations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage
%	[s, cost, iter_time ]=greed_omp(x,P,m,'option_name','option_value')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input
%   Mandatory:
%               x   Observation vector to be decomposed
%               P   Either:
%                       1) An nxm matrix (n must be dimension of x)
%                       2) A function handle (type "help function_format" 
%                          for more information)
%                       3) An object handle (type "help object_format" for 
%                          more information)
%               m   length of s 
%
%   Possible additional options:
%   (specify as many as you want using 'option_name','option_value' pairs)
%   See below for explanation of options:
%__________________________________________________________________________
%   option_name    |     available option_values                | default
%--------------------------------------------------------------------------
%   solver         | auto, qr, chol, cgp, cg, pinv, linsolve    | auto
%   stopCrit       | M, corr, mse, mse_change                   | M
%   stopTol        | number (see below)                         | n/4
%   P_trans        | function_handle (see below)                | 
%   maxIter        | positive integer (see below)               | n
%   verbose        | true, false                                | false
%   start_val      | vector of length m                         | zeros
%
%
% Explanation of possible option
%
%   solver: different implementations of OMP are available. The fastest
%           method is the QR algorithm, followed by the Cholesky based 
%           approach. The QR based algorithm requires storage of an nxM
%           matrix and an MxM triangular matrix, for long signals or 
%           large M, this can be too much. The Cholesky based approach
%           requires storage of an upper triangular matrix of size MxM,
%           however, it also requires repeated solutions of inverse systems
%           involving this matrix, which for large M can become slow.
%           Available options are:
%               auto        -   selects from qr, chol or cg depending on
%                               problem size
%               qr          -   uses QR based method
%               chol        -   Uses Cholesky based method
%               cgp         -   Uses Conjugate Gradient Pursuit 
%                               implementation (See [1] for details)
%               cg          -   Solves the inverse problem in each
%                               iteration using Conjugate Gradient
%                               algorithm. This is the only viable option
%                               to solve OMP for very large problems. 
%               pinv        -   Uses matlab PINV command to solve inverse
%                               problem in OMP iteration. (For
%                               reference only.)
%               linsolve    -   Uses matlab LINSOLVE command to solve 
%                               inverse problem in OMP iteration. (For
%                               reference only.)
%
%   stopCrit: Stopping criterion for the algorithm.
%               M           -   Extracts exactly M = stopTol elements.
%               corr        -   Stops when maximum correlation between
%                               residual and atoms is below stopTol value.
%               mse         -   Stops when mean squared error of residual 
%                               is below stopTol value.
%               mse_change  -   Stops when the change in the mean squared 
%                               error falls below stopTol value.
%
%   stopTol: Value for stopping criterion.
%
%   P_trans: If P is a function handle, then P_trans has to be specified and 
%            must be a function handle. 
%
%   maxIter: Maximum number of allowed iterations.
%
%   verbose: Logical value to allow algorithm progress to be displayed.
%
%   start_val: Allows algorithms to start from partial solution.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs
%    s              Solution vector 
%    err_norm       Vector containing norm of approximation error for each 
%                   iteration
%    iter_time      Vector containing computation times for each iteration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description
%   greed_omp performs a greedy signal decomposition. 
%   In each iteration a new element is selected depending on the inner
%   product between the current residual and columns in P.
%   The non-zero elements of s are approximated by orthogonally projecting 
%   x onto the selected elements in each iteration.
%   Different algorithms are possible to solve this projection. 
%   greed_omp has access to the following implementations (see options):
%   1) QR decomposition based algorithm.
%   2) Cholesky factorisation based algorithm.
%   3) Single Step Conjugate Gradient implementation as described in [1].
%   4) Full Conjugate Gradient solver in each iteration.
%   4) Pseudo Inverse solution in each iteration. (Not recommended!)
%   5) Use of matlab linsolve in each iteration. (Not recommended!)
%
% References
%   [1] T. Blumensath and M.E. Davies, "Gradient Pursuits", submitted, 2007
%
% See Also
%   greed_mp, greed_gp, greed_nomp, greed_omp, greed_pcgp
%
% Copyright (c) 2007 Thomas Blumensath
%
% The University of Edinburgh
% Email: [e-mail]
% Comments and bug reports welcome
%
% This file is part of sparsity Version 0.1
% Created: April 2007
%
% Part of this toolbox was developed with the support of EPSRC Grant
% D000246/1
%
% Please read COPYRIGHT.m for terms and conditions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Default values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n1 n2]=size(x);
if n2 == 1
    n=n1;
elseif n1 == 1
    x=x';
    n=n2;
else
   error('x must be a vector.');
end
    
sigsize     = x'*x/n;
initial_given=0;
err_mse     = [];
iter_time   = [];
SOLVER      = 'auto';
STOPCRIT    = 'M';
STOPTOL     = ceil(n/4);
MAXITER     = n;
verbose     = false;
s_initial   = zeros(m,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Output variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch nargout 
    case 3
        comp_err=true;
        comp_time=true;
    case 2 
        comp_err=true;
        comp_time=false;
    case 1
        comp_err=false;
        comp_time=false;
    case 0
        error('Please assign output variable.')
    otherwise
        error('Too many output arguments specified')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Look through options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OS
    if isa(varargin{i},'cell')
        CellSize=length(varargin{i});
        ThisCell=varargin{i};
        for j=1:CellSize
            Options{c}=ThisCell{j};
            c=c+1;
        end
    else
        Options{c}=varargin{i};
        c=c+1;
    end
end
OS=length(Options);
if rem(OS,2)
   error('Something is wrong with argument name and argument value pairs.') 
end
for i=1:2:OS
   switch Options{i}
        case {'solver'}  
            if isa(Options{i+1},'char'); SOLVER      = Options{i+1};
            else error('solver must be char string [auto, qr, chol, cgp, cg, pinv, linsolve]. Exiting.'); end
        case {'stopCrit'}
            if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact'));
                STOPCRIT    = Options{i+1};