changeset 163:855025f4c779 ivand_dev

renaiming small_cgp to small_pcgp
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Wed, 07 Sep 2011 14:16:50 +0100
parents f42aa8bcb82f
children 4205744092e6
files examples/ALPS solvers tests/SMALL_solver_test_ALPS.m examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m examples/SMALL_solver_test.m examples/SMALL_solver_test_Audio.m solvers/SMALL_pcgp.m
diffstat 6 files changed, 118 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/examples/ALPS solvers tests/SMALL_solver_test_ALPS.m	Wed Aug 31 12:02:19 2011 +0100
+++ b/examples/ALPS solvers tests/SMALL_solver_test_ALPS.m	Wed Sep 07 14:16:50 2011 +0100
@@ -88,7 +88,7 @@
 % SMALL Conjugate Gradient test 
 SMALL.solver(i)=SMALL_init_solver;
 SMALL.solver(i).toolbox='SMALL';    
-SMALL.solver(i).name='SMALL_cgp';
+SMALL.solver(i).name='SMALL_pcgp';
 
 % In the following string all parameters except matrix, measurement vector
 % and size of solution need to be specified. If you are not sure which
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m	Wed Aug 31 12:02:19 2011 +0100
+++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m	Wed Sep 07 14:16:50 2011 +0100
@@ -9,10 +9,7 @@
 %   -   KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient
 %              Implementation of the K-SVD Algorithm using Batch Orthogonal
 %              Matching Pursuit", Technical Report - CS, Technion, April 2008.
-%   -   RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares
-%       Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on,
-%       vol.58, no.4, pp.2121-2130, April 2010
-%
+
 
 
 %   Centre for Digital Music, Queen Mary, University of London.
--- a/examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m	Wed Aug 31 12:02:19 2011 +0100
+++ b/examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m	Wed Sep 07 14:16:50 2011 +0100
@@ -89,7 +89,7 @@
 %   Defining the parameters needed for sparse representation
 
 SMALL.solver(1).toolbox='SMALL';
-SMALL.solver(1).name='SMALL_cgp';
+SMALL.solver(1).name='SMALL_pcgp';
 
 %   Here we use mexLasso mode=2, with lambda=2, lambda2=0 and positivity
 %   constrain (type 'help mexLasso' for more information about modes):
--- a/examples/SMALL_solver_test.m	Wed Aug 31 12:02:19 2011 +0100
+++ b/examples/SMALL_solver_test.m	Wed Sep 07 14:16:50 2011 +0100
@@ -81,7 +81,7 @@
 % SMALL Conjugate Gradient test 
 SMALL.solver(i)=SMALL_init_solver;
 SMALL.solver(i).toolbox='SMALL';    
-SMALL.solver(i).name='SMALL_cgp';
+SMALL.solver(i).name='SMALL_pcgp';
 
 % In the following string all parameters except matrix, measurement vector
 % and size of solution need to be specified. If you are not sure which
--- a/examples/SMALL_solver_test_Audio.m	Wed Aug 31 12:02:19 2011 +0100
+++ b/examples/SMALL_solver_test_Audio.m	Wed Sep 07 14:16:50 2011 +0100
@@ -61,7 +61,7 @@
 
 SMALL.solver(i)=SMALL_init_solver;
 SMALL.solver(i).toolbox='SMALL';    
-SMALL.solver(i).name='SMALL_cgp';
+SMALL.solver(i).name='SMALL_pcgp';
 
 % In the following string all parameters except matrix, measurement vector
 % and size of solution need to be specified. If you are not sure which
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/solvers/SMALL_pcgp.m	Wed Sep 07 14:16:50 2011 +0100
@@ -0,0 +1,113 @@
+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;
+
+
+