daniele@171: function [dic mus] = shrinkgram(dic,mu,dd1,dd2,params)
daniele@171: % grassmanian attempts to create an n by m matrix with minimal mutual
daniele@171: % coherence using an iterative projection method.
daniele@171: %
daniele@171: % [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
daniele@171: %
daniele@171: % REFERENCE
daniele@171: % M. Elad, Sparse and Redundant Representations, Springer 2010.
daniele@171: 
daniele@171: %% Parameters and Defaults
daniele@171: if ~nargin, testshrinkgram; return; end
daniele@171: 
daniele@171: if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.9; end %shrinking factor
daniele@171: if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
daniele@171: if ~exist('params','var') || isempty(params), params = struct; end
daniele@171: if ~isfield(params,'nIter'), params.nIter = 100; end
daniele@171: 
daniele@171: %% Main algo
daniele@171: dic = normc(dic);	%normalise columns
daniele@171: G = dic'*dic;		%gram matrix
daniele@171: [n m] = size(dic);
daniele@171: 
daniele@171: MU  = @(G) max(max(abs(G-diag(diag(G))))); %coherence function
daniele@171: 
daniele@171: mus   = ones(params.nIter,1);
daniele@171: iIter = 1;
daniele@171: % optimise gram matrix
daniele@171: while iIter<=params.nIter && MU(G)>mu
daniele@171: 	mus(iIter) = MU(G);		%calculate coherence
daniele@171: 	gg  = sort(abs(G(:)));	%sort inner products from less to most correlated
daniele@171: 	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
daniele@171: 	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
daniele@171: 	[U S V] = svd(G);		%compute new SVD of gram matrix
daniele@171: 	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
daniele@171: 	G = U*S*V';				%update gram matrix
daniele@171: 	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
daniele@171: 	iIter = iIter+1;
daniele@171: end
daniele@171: %if iIter<params.nIter
daniele@171: %	mus(iIter:end) = mus(iIter-1);
daniele@171: %end
daniele@171: 
daniele@171: [V_gram Sigma_gram] = svd(G);				%calculate svd decomposition of gramian
daniele@171: dic = sqrt(Sigma_gram(1:n,:))*V_gram';		%update dictionary
daniele@171: 
daniele@171: function testshrinkgram
daniele@171: clc
daniele@171: %define parameters
daniele@171: n = 256;							%ambient dimension
daniele@171: m = 512;							%number of atoms
daniele@171: N = 1024;							%number of signals
daniele@171: mu_min = sqrt((m-n)/(n*(m-1)));		%minimum coherence
daniele@171: 
daniele@171: %initialise data
daniele@171: phi = normc(randn(n,m));			%dictionary
daniele@171: 
daniele@171: %optimise dictionary
daniele@171: [~, mus] = shrinkgram(phi,0.2);
daniele@171: 
daniele@171: %plot results
daniele@171: nIter = length(mus);
daniele@171: 
daniele@171: figure, hold on
daniele@171: plot(1:nIter,mus,'ko-');
daniele@171: plot([1 nIter],[mu_min mu_min],'k')
daniele@171: grid on
daniele@171: legend('\mu','\mu_{min}');