Mercurial > hg > smallbox
view solvers/SMALL_pcgp.m @ 217:8b3c71bb44eb luisf_dev
Removed "clear all" from example scripts (subs by "clear" instead)
author | luisf <luis.figueira@eecs.qmul.ac.uk> |
---|---|
date | Thu, 22 Mar 2012 14:41:04 +0000 |
parents | 855025f4c779 |
children |
line wrap: on
line source
function [A, resF]=SMALL_pcgp(Dict,X, m, maxNumCoef, errorGoal, varargin) %% Partial Conjugate Gradient Pursuit Multiple vectors sparse representation % % Sparse coding of a group of signals based on a given % dictionary and specified number of atoms to use. % input arguments: Dict - the dictionary % X - the signals to represent % m - number of atoms in Dictionary % errorGoal - the maximal allowed representation error for % each signal. % % optional: if Dict is function handle then Transpose Dictionary % handle needs to be specified. % % output arguments: A - sparse coefficient matrix. % % Centre for Digital Music, Queen Mary, University of London. % This file copyright 2009 Ivan Damnjanovic. % % This program is free software; you can redistribute it and/or % modify it under the terms of the GNU General Public License as % published by the Free Software Foundation; either version 2 of the % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. % %% % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers if isa(Dict,'float') D =@(z) Dict*z; Dt =@(z) Dict'*z; elseif isobject(Dict) D =@(z) Dict*z; Dt =@(z) Dict'*z; elseif isa(Dict,'function_handle') try DictT=varargin{1}; if isa(DictT,'function_handle'); D=Dict; Dt=DictT; else error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. '); end catch error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.'); end else error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.'); end %% %positivity=1; [n,P]=size(X); A = sparse(m,size(X,2)); for k=1:1:P, x = X(:,k); residual = x; indx =[]; j=0; currResNorm = residual'*residual/n; errGoal=errorGoal*currResNorm; a = zeros(m,1); p = zeros(m,1); while j < maxNumCoef, j = j+1; dir=Dt(residual); if exist('positivity','var')&&(positivity==1) [tmp__, pos]=max(dir); else [tmp__, pos]=max(abs(dir)); end indx(j)=pos; p(indx)=dir(indx); Dp=D(p); pDDp=Dp'*Dp; if (j==1) beta=0; else beta=-Dp'*residual/pDDp; end alfa=residual'*Dp/pDDp; a=a+alfa*p; p(indx)=dir(indx)+beta*p(indx); residual=residual-alfa*Dp; currResNorm = residual'*residual/n; if currResNorm<errGoal fprintf('\nFound exact representation! \n'); break; end end; if (~isempty(indx)) resF(k)=currResNorm; A(indx,k)=a(indx); end end; return;