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