comparison solvers/SMALL_cgp.m @ 1:7750624e0c73 version0.5

(none)

+function [A]=SMALL_cgp(Dict,X, m,  maxNumCoef, errorGoal, varargin) 
+%=============================================
+% 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.
+%
+% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
+% with pinv() is changed with conjugate gradient method)
+% Ivan Damnjanovic 2009
+%=============================================
+%%
+%%
+
+% 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
+%%
+
+[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;
+    errorGoal=errorGoal*currResNorm;
+    a = zeros(m,1);
+    p = zeros(m,1);
+    
+    while  j < maxNumCoef,
+        
+        j = j+1;
+        dir=Dt(residual);
+        
+        [tmp__, pos]=max(abs(dir));
+        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<errorGoal
+            fprintf('\nFound exact representation! \n');
+            break;
+        end
+   end;
+   if (~isempty(indx))
+       A(indx,k)=a(indx);
+   end
+end;
+return;
+
+
+