comparison util/classes/dictionaryMatrices/iterativeprojections.m @ 170:68fb71aa5339 danieleb

Added dictionary decorrelation functions and test script for Letters paper.

-function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
+function [dic mus srs] = iterativeprojections(dic,mu,Y,X,params)
 % grassmanian attempts to create an n by m matrix with minimal mutual
 % coherence using an iterative projection method.
 %
-% [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
+% REFERENCE
 %
-% REFERENCE
-% M. Elad, Sparse and Redundant Representations, Springer 2010.
 
 %% Parameters and Defaults
-error(nargchk(2,7,nargin));
+if ~nargin, testiterativeprojections; return; end
 
-if ~exist('verb','var')  || isempty(verb),  verb = false; end %verbose output
-if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
-if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.99; end %shrinking factor
-if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
-if ~exist('nIter','var') || isempty(nIter), nIter = 10; end %number of iterations
+if ~exist('params','var') || isempty(param), params = struct; end
+if ~isfield(params,'nIter'), params.nIter = 10; end %number of iterations
+if ~isfield(params,'eps'),	 params.eps = 1e-9;	 end %tolerance level
+[n m] = size(dic);
 
-%% Main algo
-A = normc(initA); %normalise columns
-[Uinit Sigma] = svd(A);
-G = A'*A; %gram matrix
+SNR = @(dic) snr(Y,dic*X);									 %SNR function
+MU  = @(dic) max(max(abs((dic'*dic)-diag(diag(dic'*dic))))); %coherence function
 
-muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence (equiangular tight frame)
-res = zeros(nIter,1);
-if verb
-	fprintf(1,'Iter		mu_min		mu \n');
+%% Main algorithm
+dic = normc(dic);	%normalise columns
+alpha = m/n;		%ratio between number of atoms and ambient dimension
+
+mus = zeros(params.nIter,1);		%coherence at each iteration
+srs = zeros(params.nIter,1);		%signal to noise ratio at each iteration
+iIter = 1;
+while iIter<=params.nIter && MU(dic)>mu
+	fprintf(1,'Iteration number %u\n', iIter);
+	% calculate snr and coherence
+	mus(iIter) = MU(dic);
+	srs(iIter) = SNR(dic);
+	
+	% calculate gram matrix
+	G = dic'*dic;
+	
+	% project into the structural constraint set
+	H = zeros(size(G));				%initialise matrix
+	ind1 = find(abs(G(:))<=mu);		%find elements smaller than mu
+	ind2 = find(abs(G(:))>mu);		%find elements bigger than mu
+	H(ind1) = G(ind1);				%copy elements belonging to ind1
+	H(ind2) = mu*sign(G(ind2));		%threshold elements belonging to ind2
+	H(1:m+1:end) = 1;				%set diagonal to one
+	
+	% project into spectral constraint set
+	[~ , S, V] = svd(H);
+	%G = alpha*(V(:,1:n)*V(:,1:n)');
+	G = V(:,1:n)*S(1:n,1:n)*V(:,1:n)';
+	
+	% calculate dictionary
+	[~, S V] = svd(G);
+	dic = sqrt(S(1:n,:))*V';
+	
+	% rotate dictionary
+	options = struct('nIter',100,'step',0.001);
+	[~, ~, W] = rotatematrix(Y,dic*X,'conjgradlie',options);
+	dic = W*dic;
+	
+	iIter = iIter+1;
 end
+if iIter<params.nIter
+	mus(iIter:end) = mus(iIter);
+	srs(iIter:end) = srs(iIter);
+end
+% Test function
+function testiterativeprojections
+clc
+%define parameters
+n = 256;							%ambient dimension
+m = 512;							%number of atoms
+N = 1024;							%number of signals
+mu_min = sqrt((m-n)/(n*(m-1)));		%minimum coherence
 
-% optimise gram matrix
-for iIter = 1:nIter
-	gg  = sort(abs(G(:))); %sort inner products from less to most correlated
-	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
-	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
-	[U S V] = svd(G);	%compute new SVD of gram matrix
-	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
-	G = U*S*V';			%update gram matrix
-	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
-	if verb
-		Geye = G - eye(size(G));
-		fprintf(1,'%6i    %12.8f  %12.8f  \n',iIter,muMin,max(abs(Geye(:))));
-	end
-end
+%initialise data
+X = sprandn(m,N,1);					%matrix of coefficients
+phi = normc(randn(n,m));			%dictionary
+temp = randn(n);
+W = expm(0.5*(temp-temp'));			%rotation matrix
+Y = W*phi*X;						%observed signals
 
-% [~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
+%optimise dictionary
+[~, mus srs] = iterativeprojections(phi,0.2,Y,X);
 
-% A = normc(A);					%normalise dictionary
+%plot results
+nIter = length(mus);
 
-[V_gram Sigma_gram] = svd(G);				%calculate svd decomposition of gramian
-Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma); %calculate singular values of dictionary
-A = Uinit*Sigma_new*V_gram';				%update dictionary
+figure, 
+subplot(2,1,1)
+plot(1:nIter,srs,'kd-');
+xlabel('nIter');
+ylabel('snr (dB)');
+grid on
 
-% param.step  = 0.01;
-% param.reg   = 0.01;
-% param.nIter = 20;
-% A = rotatematrix(initA,A,'linesearchlie',param);
+subplot(2,1,2), hold on
+plot(1:nIter,mus,'ko-');
+plot([1 nIter],[mu_min mu_min],'k')
+grid on
+legend('\mu','\mu_{min}');
 
-% %% Debug visualization function
-% function plotcart2d(A)
-% compass(A(1,:),A(2,:));