annotate solvers/SMALL_cgp.m @ 62:92a9b18c3bd3

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