daniele@160: function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
daniele@160: % grassmanian attempts to create an n by m matrix with minimal mutual
daniele@160: % coherence using an iterative projection method.
daniele@160: %
daniele@160: % [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
daniele@160: %
daniele@166: % REFERENCE
daniele@166: % M. Elad, Sparse and Redundant Representations, Springer 2010.
daniele@166: 
daniele@160: %% Parameters and Defaults
daniele@160: error(nargchk(2,7,nargin));
daniele@160: 
daniele@160: if ~exist('verb','var')  || isempty(verb),  verb = false; end %verbose output
daniele@160: if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
daniele@166: if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.99; end %shrinking factor
daniele@160: if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
daniele@162: if ~exist('nIter','var') || isempty(nIter), nIter = 10; end %number of iterations
daniele@160: 
daniele@162: %% Main algo
daniele@160: A = normc(initA); %normalise columns
daniele@166: [Uinit Sigma] = svd(A);
daniele@162: G = A'*A; %gram matrix
daniele@160: 
daniele@162: muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence (equiangular tight frame)
daniele@160: res = zeros(nIter,1);
daniele@162: if verb
daniele@162: 	fprintf(1,'Iter		mu_min		mu \n');
daniele@160: end
daniele@160: 
daniele@162: % optimise gram matrix
daniele@162: for iIter = 1:nIter
daniele@162: 	gg  = sort(abs(G(:))); %sort inner products from less to most correlated
daniele@162: 	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
daniele@162: 	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
daniele@162: 	[U S V] = svd(G);	%compute new SVD of gram matrix
daniele@162: 	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
daniele@162: 	G = U*S*V';			%update gram matrix
daniele@162: 	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
daniele@162: 	if verb
daniele@162: 		Geye = G - eye(size(G));
daniele@162: 		fprintf(1,'%6i    %12.8f  %12.8f  \n',iIter,muMin,max(abs(Geye(:))));
daniele@162: 	end
daniele@162: end
daniele@162: 
daniele@162: % [~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
daniele@166: 
daniele@162: % A = normc(A);					%normalise dictionary
daniele@162: 
daniele@166: [V_gram Sigma_gram] = svd(G);				%calculate svd decomposition of gramian
daniele@166: Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma); %calculate singular values of dictionary
daniele@169: A = Uinit*Sigma_new*V_gram';				%update dictionary
daniele@169: 
daniele@169: % param.step  = 0.01;
daniele@169: % param.reg   = 0.01;
daniele@169: % param.nIter = 20;
daniele@169: % A = rotatematrix(initA,A,'linesearchlie',param);
daniele@162: 
daniele@162: % %% Debug visualization function
daniele@162: % function plotcart2d(A)
daniele@162: % compass(A(1,:),A(2,:));