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