annotate solvers/SMALL_cgp.m @ 247:ecce33192fcc tip

Added tag ver_2.1 for changeset cef4500b936f
author luisf <luis.figueira@eecs.qmul.ac.uk>
date Wed, 31 Oct 2012 12:24:44 +0000
parents 8e660fd14774
children
rev   line source
idamnjanovic@3 1 function [A, resF]=SMALL_cgp(Dict,X, m, maxNumCoef, errorGoal, varargin)
ivan@128 2 %% Partial Conjugate Gradient Pursuit Multiple vectors sparse representation
ivan@128 3 %
ivan@128 4 % Sparse coding of a group of signals based on a given
ivan@128 5 % dictionary and specified number of atoms to use.
ivan@128 6 % input arguments: Dict - the dictionary
ivan@128 7 % X - the signals to represent
ivan@128 8 % m - number of atoms in Dictionary
ivan@128 9 % errorGoal - the maximal allowed representation error for
ivan@128 10 % each signal.
ivan@128 11 %
ivan@128 12 % optional: if Dict is function handle then Transpose Dictionary
ivan@128 13 % handle needs to be specified.
ivan@128 14 %
ivan@128 15 % output arguments: A - sparse coefficient matrix.
ivan@128 16
idamnjanovic@22 17 %
idamnjanovic@22 18 % Centre for Digital Music, Queen Mary, University of London.
idamnjanovic@22 19 % This file copyright 2009 Ivan Damnjanovic.
idamnjanovic@22 20 %
idamnjanovic@22 21 % This program is free software; you can redistribute it and/or
idamnjanovic@22 22 % modify it under the terms of the GNU General Public License as
idamnjanovic@22 23 % published by the Free Software Foundation; either version 2 of the
idamnjanovic@22 24 % License, or (at your option) any later version. See the file
idamnjanovic@22 25 % COPYING included with this distribution for more information.
idamnjanovic@1 26 %
idamnjanovic@1 27 %%
idamnjanovic@1 28
idamnjanovic@1 29 % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers
idamnjanovic@1 30 if isa(Dict,'float')
idamnjanovic@1 31 D =@(z) Dict*z;
idamnjanovic@1 32 Dt =@(z) Dict'*z;
idamnjanovic@1 33 elseif isobject(Dict)
idamnjanovic@1 34 D =@(z) Dict*z;
idamnjanovic@1 35 Dt =@(z) Dict'*z;
idamnjanovic@1 36 elseif isa(Dict,'function_handle')
idamnjanovic@1 37 try
idamnjanovic@1 38 DictT=varargin{1};
idamnjanovic@1 39 if isa(DictT,'function_handle');
idamnjanovic@73 40 D=Dict;
idamnjanovic@73 41 Dt=DictT;
idamnjanovic@1 42 else
idamnjanovic@1 43 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1 44 end
idamnjanovic@1 45 catch
idamnjanovic@1 46 error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.');
idamnjanovic@1 47 end
idamnjanovic@1 48 else
idamnjanovic@1 49 error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1 50 end
idamnjanovic@1 51 %%
idamnjanovic@62 52 %positivity=1;
idamnjanovic@1 53 [n,P]=size(X);
idamnjanovic@1 54
idamnjanovic@1 55
idamnjanovic@1 56
idamnjanovic@1 57 A = sparse(m,size(X,2));
idamnjanovic@1 58
idamnjanovic@1 59 for k=1:1:P,
idamnjanovic@1 60
idamnjanovic@1 61 x = X(:,k);
idamnjanovic@1 62 residual = x;
idamnjanovic@1 63 indx =[];
idamnjanovic@1 64 j=0;
idamnjanovic@1 65
idamnjanovic@1 66
idamnjanovic@1 67 currResNorm = residual'*residual/n;
idamnjanovic@3 68 errGoal=errorGoal*currResNorm;
idamnjanovic@1 69 a = zeros(m,1);
idamnjanovic@1 70 p = zeros(m,1);
idamnjanovic@1 71
idamnjanovic@1 72 while j < maxNumCoef,
idamnjanovic@1 73
idamnjanovic@1 74 j = j+1;
idamnjanovic@1 75 dir=Dt(residual);
idamnjanovic@22 76 if exist('positivity','var')&&(positivity==1)
idamnjanovic@22 77 [tmp__, pos]=max(dir);
idamnjanovic@22 78 else
idamnjanovic@22 79 [tmp__, pos]=max(abs(dir));
idamnjanovic@22 80 end
idamnjanovic@1 81 indx(j)=pos;
idamnjanovic@1 82
idamnjanovic@1 83 p(indx)=dir(indx);
idamnjanovic@1 84 Dp=D(p);
idamnjanovic@1 85 pDDp=Dp'*Dp;
idamnjanovic@1 86
idamnjanovic@1 87 if (j==1)
idamnjanovic@1 88 beta=0;
idamnjanovic@1 89 else
idamnjanovic@1 90 beta=-Dp'*residual/pDDp;
idamnjanovic@1 91 end
idamnjanovic@1 92
idamnjanovic@1 93 alfa=residual'*Dp/pDDp;
idamnjanovic@1 94 a=a+alfa*p;
idamnjanovic@1 95 p(indx)=dir(indx)+beta*p(indx);
idamnjanovic@1 96
idamnjanovic@1 97 residual=residual-alfa*Dp;
idamnjanovic@1 98
idamnjanovic@1 99 currResNorm = residual'*residual/n;
idamnjanovic@3 100 if currResNorm<errGoal
idamnjanovic@1 101 fprintf('\nFound exact representation! \n');
idamnjanovic@1 102 break;
idamnjanovic@1 103 end
idamnjanovic@1 104 end;
idamnjanovic@1 105 if (~isempty(indx))
idamnjanovic@3 106 resF(k)=currResNorm;
idamnjanovic@1 107 A(indx,k)=a(indx);
idamnjanovic@1 108 end
idamnjanovic@1 109 end;
idamnjanovic@1 110 return;
idamnjanovic@1 111
idamnjanovic@1 112
idamnjanovic@1 113