| [s, err_mse, iter_time]=block_gp(x,A,group,varargin) |
function [s, err_mse, iter_time]=block_gp(x,A,group,varargin)
% block_gp: Block Gradient Pursuit algorithm (modification from [1])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage
% [s, err_mse, iter_time]=block_gp(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)
% Also requires specification of P_trans option.
% 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
%--------------------------------------------------------------------------
% 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
% GradSteps | 'auto' or integer | 'auto'
%
% Available stopping criteria :
% 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 of allowed iterations.
%
% verbose: Logical value to allow algorithm progress to be displayed.
%
% start_val: Allows algorithms to start from partial solution.
%
% GradSteps: Number of gradient optimisation steps per iteration.
% 'auto' uses inner products to decide if more gradient steps
% are required.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs
% s Solution vector
% err_mse Vector containing mse of approximation error for each
% iteration
% iter_time Vector containing times for each iteration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% References
% [1] T. Blumensath and M.E. Davies, "Gradient Pursuits", submitted, 2007
%
%
% Original code by Thomas Blumensath 2007
% Modified for Block Sparsity by Angshul Majumdar 2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values and initialisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = size(A,2);
[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 = [];
STOPCRIT = 'M';
STOPTOL = ceil(n/4);
MAXITER = n;
verbose = false;
s_initial = zeros(m,1);
GradSteps = 'auto';
if verbose
display('Initialising...')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 {'stopCrit'}
if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact'));
STOPCRIT = Options{i+1};
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};
initial_given=1;
else error('start_val must be a vector of length m. Exiting.'); end
case {'GradSteps'}
if isa(Options{i+1},'numeric') || strcmp(Options{i+1},'auto') ;
GradSteps = Options{i+1};
else error('start_val must be a vector of length m. Exiting.'); end
otherwise
error('Unrecognised option. Exiting.')
end
end
if strcmp(STOPCRIT,'M')
maxM=STOPTOL;
else
maxM=MAXITER;
end
if nargout >=2
err_mse = zeros(maxM,1);
end
if nargout ==3
iter_time = zeros(maxM,1);
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do we start from zero or not?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if initial_given ==1;
IN = find(s_initial);
Residual = x-P(s_initial);
s = s_initial;
oldERR = Residual'*Residual/n;
else
IN = [];
Residual = x;
s = s_initial;
sigsize = x'*x/n;
oldERR = sigsize;
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
c = max(group);
for j = 1:c
g{j} = find(group == j);
end
tic
t=0;
p=zeros(m,1);
DR=Pt(Residual);
[v L]=max(abs(DR));
I = g{group(L)};
IN=[IN I'];
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
p(IN)=DR(IN);
% Step size
Dp=P(p);
a=Residual'*Dp/(Dp'*Dp);
% Update coefficients
s=s+a*p;
% New Residual and inner products
Residual=Residual-a*Dp;
DR=Pt(Residual);
% select new element
[v L]=max(abs(DR));
I = g{group(L)};
if isempty(intersect(IN,I))
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
p(IN)=DR(IN);
% Step size
Dp=P(p);
a=Residual'*Dp/(Dp'*Dp);
% 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 L]=max(abs(DR));
I = g{group(L)};
IN=[IN I'];
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 iter >= STOPTOL
done =1;
elseif verbose && toc-t>10
display(sprintf('Iteration %i. --- %i iterations to go',iter ,STOPTOL-iter))
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
%
% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.
|
|
