daniele@171
|
1 function [dic mus] = shrinkgram(dic,mu,dd1,dd2,params)
|
daniele@171
|
2 % grassmanian attempts to create an n by m matrix with minimal mutual
|
daniele@171
|
3 % coherence using an iterative projection method.
|
daniele@171
|
4 %
|
daniele@171
|
5 % [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
|
daniele@171
|
6 %
|
daniele@171
|
7 % REFERENCE
|
daniele@171
|
8 % M. Elad, Sparse and Redundant Representations, Springer 2010.
|
daniele@171
|
9
|
daniele@171
|
10 %% Parameters and Defaults
|
daniele@171
|
11 if ~nargin, testshrinkgram; return; end
|
daniele@171
|
12
|
daniele@171
|
13 if ~exist('dd2','var') || isempty(dd2), dd2 = 0.9; end %shrinking factor
|
daniele@171
|
14 if ~exist('dd1','var') || isempty(dd1), dd1 = 0.9; end %percentage of coherences to be shrinked
|
daniele@171
|
15 if ~exist('params','var') || isempty(params), params = struct; end
|
daniele@171
|
16 if ~isfield(params,'nIter'), params.nIter = 100; end
|
daniele@171
|
17
|
daniele@171
|
18 %% Main algo
|
daniele@171
|
19 dic = normc(dic); %normalise columns
|
daniele@171
|
20 G = dic'*dic; %gram matrix
|
daniele@171
|
21 [n m] = size(dic);
|
daniele@171
|
22
|
daniele@171
|
23 MU = @(G) max(max(abs(G-diag(diag(G))))); %coherence function
|
daniele@171
|
24
|
daniele@171
|
25 mus = ones(params.nIter,1);
|
daniele@171
|
26 iIter = 1;
|
daniele@171
|
27 % optimise gram matrix
|
daniele@171
|
28 while iIter<=params.nIter && MU(G)>mu
|
daniele@171
|
29 mus(iIter) = MU(G); %calculate coherence
|
daniele@171
|
30 gg = sort(abs(G(:))); %sort inner products from less to most correlated
|
daniele@171
|
31 pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
|
daniele@171
|
32 G(pos) = G(pos)*dd2; %shrink large elements of gram matrix
|
daniele@171
|
33 [U S V] = svd(G); %compute new SVD of gram matrix
|
daniele@171
|
34 S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
|
daniele@171
|
35 G = U*S*V'; %update gram matrix
|
daniele@171
|
36 G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
|
daniele@171
|
37 iIter = iIter+1;
|
daniele@171
|
38 end
|
daniele@171
|
39 %if iIter<params.nIter
|
daniele@171
|
40 % mus(iIter:end) = mus(iIter-1);
|
daniele@171
|
41 %end
|
daniele@171
|
42
|
daniele@171
|
43 [V_gram Sigma_gram] = svd(G); %calculate svd decomposition of gramian
|
daniele@171
|
44 dic = sqrt(Sigma_gram(1:n,:))*V_gram'; %update dictionary
|
daniele@171
|
45
|
daniele@171
|
46 function testshrinkgram
|
daniele@171
|
47 clc
|
daniele@171
|
48 %define parameters
|
daniele@171
|
49 n = 256; %ambient dimension
|
daniele@171
|
50 m = 512; %number of atoms
|
daniele@171
|
51 N = 1024; %number of signals
|
daniele@171
|
52 mu_min = sqrt((m-n)/(n*(m-1))); %minimum coherence
|
daniele@171
|
53
|
daniele@171
|
54 %initialise data
|
daniele@171
|
55 phi = normc(randn(n,m)); %dictionary
|
daniele@171
|
56
|
daniele@171
|
57 %optimise dictionary
|
daniele@171
|
58 [~, mus] = shrinkgram(phi,0.2);
|
daniele@171
|
59
|
daniele@171
|
60 %plot results
|
daniele@171
|
61 nIter = length(mus);
|
daniele@171
|
62
|
daniele@171
|
63 figure, hold on
|
daniele@171
|
64 plot(1:nIter,mus,'ko-');
|
daniele@171
|
65 plot([1 nIter],[mu_min mu_min],'k')
|
daniele@171
|
66 grid on
|
daniele@171
|
67 legend('\mu','\mu_{min}');
|