annotate solvers/SMALL_cgp.m @ 5:f44689e95ea4

(none)
author idamnjanovic
date Mon, 22 Mar 2010 10:43:01 +0000
parents 01cad25206d6
children 524cc3fff5ac
rev   line source
idamnjanovic@3 1 function [A, resF]=SMALL_cgp(Dict,X, m, maxNumCoef, errorGoal, varargin)
idamnjanovic@1 2 %=============================================
idamnjanovic@1 3 % Sparse coding of a group of signals based on a given
idamnjanovic@1 4 % dictionary and specified number of atoms to use.
idamnjanovic@1 5 % input arguments: Dict - the dictionary
idamnjanovic@1 6 % X - the signals to represent
idamnjanovic@1 7 % m - number of atoms in Dictionary
idamnjanovic@1 8 % errorGoal - the maximal allowed representation error for
idamnjanovic@1 9 % each signal.
idamnjanovic@1 10 %
idamnjanovic@1 11 % optional: if Dict is function handle then Transpose Dictionary
idamnjanovic@1 12 % handle needs to be specified.
idamnjanovic@1 13 %
idamnjanovic@1 14 % output arguments: A - sparse coefficient matrix.
idamnjanovic@1 15 %
idamnjanovic@1 16 % based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
idamnjanovic@1 17 % with pinv() is changed with conjugate gradient method)
idamnjanovic@1 18 % Ivan Damnjanovic 2009
idamnjanovic@1 19 %=============================================
idamnjanovic@1 20 %%
idamnjanovic@1 21 %%
idamnjanovic@1 22
idamnjanovic@1 23 % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers
idamnjanovic@1 24 if isa(Dict,'float')
idamnjanovic@1 25 D =@(z) Dict*z;
idamnjanovic@1 26 Dt =@(z) Dict'*z;
idamnjanovic@1 27 elseif isobject(Dict)
idamnjanovic@1 28 D =@(z) Dict*z;
idamnjanovic@1 29 Dt =@(z) Dict'*z;
idamnjanovic@1 30 elseif isa(Dict,'function_handle')
idamnjanovic@1 31 try
idamnjanovic@1 32 DictT=varargin{1};
idamnjanovic@1 33 if isa(DictT,'function_handle');
idamnjanovic@1 34 D=Dict;
idamnjanovic@1 35 Dt=DictT;
idamnjanovic@1 36 else
idamnjanovic@1 37 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1 38 end
idamnjanovic@1 39 catch
idamnjanovic@1 40 error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.');
idamnjanovic@1 41 end
idamnjanovic@1 42 else
idamnjanovic@1 43 error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1 44 end
idamnjanovic@1 45 %%
idamnjanovic@1 46
idamnjanovic@1 47 [n,P]=size(X);
idamnjanovic@1 48
idamnjanovic@1 49
idamnjanovic@1 50
idamnjanovic@1 51 A = sparse(m,size(X,2));
idamnjanovic@1 52
idamnjanovic@1 53 for k=1:1:P,
idamnjanovic@1 54
idamnjanovic@1 55 x = X(:,k);
idamnjanovic@1 56 residual = x;
idamnjanovic@1 57 indx =[];
idamnjanovic@1 58 j=0;
idamnjanovic@1 59
idamnjanovic@1 60
idamnjanovic@1 61 currResNorm = residual'*residual/n;
idamnjanovic@3 62 errGoal=errorGoal*currResNorm;
idamnjanovic@1 63 a = zeros(m,1);
idamnjanovic@1 64 p = zeros(m,1);
idamnjanovic@1 65
idamnjanovic@1 66 while j < maxNumCoef,
idamnjanovic@1 67
idamnjanovic@1 68 j = j+1;
idamnjanovic@1 69 dir=Dt(residual);
idamnjanovic@1 70
idamnjanovic@1 71 [tmp__, pos]=max(abs(dir));
idamnjanovic@1 72 indx(j)=pos;
idamnjanovic@1 73
idamnjanovic@1 74 p(indx)=dir(indx);
idamnjanovic@1 75 Dp=D(p);
idamnjanovic@1 76 pDDp=Dp'*Dp;
idamnjanovic@1 77
idamnjanovic@1 78 if (j==1)
idamnjanovic@1 79 beta=0;
idamnjanovic@1 80 else
idamnjanovic@1 81 beta=-Dp'*residual/pDDp;
idamnjanovic@1 82 end
idamnjanovic@1 83
idamnjanovic@1 84 alfa=residual'*Dp/pDDp;
idamnjanovic@1 85 a=a+alfa*p;
idamnjanovic@1 86 p(indx)=dir(indx)+beta*p(indx);
idamnjanovic@1 87
idamnjanovic@1 88 residual=residual-alfa*Dp;
idamnjanovic@1 89
idamnjanovic@1 90 currResNorm = residual'*residual/n;
idamnjanovic@3 91 if currResNorm<errGoal
idamnjanovic@1 92 fprintf('\nFound exact representation! \n');
idamnjanovic@1 93 break;
idamnjanovic@1 94 end
idamnjanovic@1 95 end;
idamnjanovic@1 96 if (~isempty(indx))
idamnjanovic@3 97 resF(k)=currResNorm;
idamnjanovic@1 98 A(indx,k)=a(indx);
idamnjanovic@1 99 end
idamnjanovic@1 100 end;
idamnjanovic@1 101 return;
idamnjanovic@1 102
idamnjanovic@1 103
idamnjanovic@1 104