annotate util/classes/dictionaryMatrices/grassmannian.m @ 175:9eb5f0d4c1a4 danieleb

added MOCOD dictionary update
author Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
date Thu, 17 Nov 2011 11:17:00 +0000
parents 68fb71aa5339
children
rev   line source
daniele@160 1 function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
daniele@160 2 % grassmanian attempts to create an n by m matrix with minimal mutual
daniele@160 3 % coherence using an iterative projection method.
daniele@160 4 %
daniele@160 5 % [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
daniele@160 6 %
daniele@166 7 % REFERENCE
daniele@166 8 % M. Elad, Sparse and Redundant Representations, Springer 2010.
daniele@166 9
daniele@160 10 %% Parameters and Defaults
daniele@160 11 error(nargchk(2,7,nargin));
daniele@160 12
daniele@160 13 if ~exist('verb','var') || isempty(verb), verb = false; end %verbose output
daniele@160 14 if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
daniele@170 15 if ~exist('dd2','var') || isempty(dd2), dd2 = 0.9; end %shrinking factor
daniele@160 16 if ~exist('dd1','var') || isempty(dd1), dd1 = 0.9; end %percentage of coherences to be shrinked
daniele@162 17 if ~exist('nIter','var') || isempty(nIter), nIter = 10; end %number of iterations
daniele@160 18
daniele@162 19 %% Main algo
daniele@160 20 A = normc(initA); %normalise columns
daniele@166 21 [Uinit Sigma] = svd(A);
daniele@162 22 G = A'*A; %gram matrix
daniele@160 23
daniele@162 24 muMin = sqrt((m-n)/(n*(m-1))); %Lower bound on mutual coherence (equiangular tight frame)
daniele@160 25 res = zeros(nIter,1);
daniele@162 26 if verb
daniele@162 27 fprintf(1,'Iter mu_min mu \n');
daniele@160 28 end
daniele@160 29
daniele@162 30 % optimise gram matrix
daniele@162 31 for iIter = 1:nIter
daniele@162 32 gg = sort(abs(G(:))); %sort inner products from less to most correlated
daniele@162 33 pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
daniele@162 34 G(pos) = G(pos)*dd2; %shrink large elements of gram matrix
daniele@162 35 [U S V] = svd(G); %compute new SVD of gram matrix
daniele@162 36 S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
daniele@162 37 G = U*S*V'; %update gram matrix
daniele@162 38 G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
daniele@162 39 if verb
daniele@162 40 Geye = G - eye(size(G));
daniele@162 41 fprintf(1,'%6i %12.8f %12.8f \n',iIter,muMin,max(abs(Geye(:))));
daniele@162 42 end
daniele@162 43 end
daniele@162 44
daniele@166 45 [V_gram Sigma_gram] = svd(G); %calculate svd decomposition of gramian
daniele@166 46 Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma); %calculate singular values of dictionary
daniele@169 47 A = Uinit*Sigma_new*V_gram'; %update dictionary