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@160: %
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@160: if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.95; end %shrinking factor
daniele@160: if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
daniele@160: if ~exist('nIter','var') || isempty(nIter), nIter = 5; end %number of iterations
daniele@160: 
daniele@160: %% Compute svd and gramian
daniele@160: A = normc(initA); %normalise columns
daniele@160: [Uinit Sigma] = svd(A); %calculate svd of the matrix
daniele@160: G = A'*A; %gramian matrix
daniele@160: 
daniele@160: muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence
daniele@160: res = zeros(nIter,1);
daniele@160: for iIter = 1:nIter
daniele@160:    gg  = sort(abs(G(:))); %sort inner products from less to ost correlated
daniele@160:    pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); 
daniele@160:    G(pos) = G(pos)*dd2;
daniele@160:    [U S V] = svd(G);
daniele@160:    S(n+1:end,1+n:end) = 0;
daniele@160:    G = U*S*V';
daniele@160:    G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G))));
daniele@160:    gg = sort(abs(G(:)));
daniele@160:    pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6);
daniele@160:    res(iIter) = max(abs(G(pos)));
daniele@160:    if verb
daniele@160:     fprintf(1,'%6i  %12.8f   %12.8f  %12.8f \n',...
daniele@160:             [iIter,muMin,mean(abs(G(pos))),max(abs(G(pos)))]);
daniele@160:    end
daniele@160: end
daniele@160: 
daniele@160: [~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
daniele@160: Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma);
daniele@160: A = Uinit*Sigma_new*V_gram';
daniele@160: A = normc(A);