Mercurial > hg > smallbox

changeset 169:290cca7d3469 danieleb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added dictionary decorrelation functions and test script for ICASSP paper.
author Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
date Thu, 29 Sep 2011 09:46:52 +0100
parents ff866a412be5
children 68fb71aa5339
files DL/two-step DL/SMALL_two_step_DL.m DL/two-step DL/dico_decorr_symetric.m examples/SMALL_test_coherence.m examples/SMALL_test_coherence2.m util/classes/dictionaryMatrices/grassmannian.m util/classes/dictionaryMatrices/iterativeprojections.m util/classes/dictionaryMatrices/rotatematrix.m util/classes/util/regularisedictionary.m
diffstat 8 files changed, 411 insertions(+), 170 deletions(-) [+]
line wrap: on
line diff
--- a/DL/two-step DL/SMALL_two_step_DL.m	Tue Sep 20 15:52:33 2011 +0100
+++ b/DL/two-step DL/SMALL_two_step_DL.m	Thu Sep 29 09:46:52 2011 +0100
@@ -30,7 +30,7 @@
 
 % initialize the dictionary %
 
-if (isfield(DL.param,'initdict'))
+if (isfield(DL.param,'initdict')) && ~isempty(DL.param.initdict);
   if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
     dico = sig(:,DL.param.initdict(1:dictsize));
   else
@@ -108,7 +108,7 @@
 % main loop %
 
 for i = 1:iternum
-    disp([num2str(i) '/' num2str(iternum)]);
+    %disp([num2str(i) '/' num2str(iternum)]);
     Problem.A = dico;
     solver = SMALL_solve(Problem, solver);
     [dico, solver.solution] = dico_update(dico, sig, solver.solution, ...
@@ -119,7 +119,13 @@
                 dico = dico_decorr_symetric(dico, mu, solver.solution);
             case 'tropp'
                 [n m] = size(dico);
-                dico = grassmannian(n,m,[],[],[],dico,true);
+                dico = grassmannian(n,m,[],[],[],dico);
+			case 'barchiesi'
+				[n m] = size(dico);
+				params.nIter = 100;
+				dico = iterativeprojections(n,m,[],[],[],dico);
+				[~, ~, W] = rotatematrix(Problem.b,dico*solver.solution,'conjgradlie',params);
+				dico = W*dico;
             otherwise
         end
     
--- a/DL/two-step DL/dico_decorr_symetric.m	Tue Sep 20 15:52:33 2011 +0100
+++ b/DL/two-step DL/dico_decorr_symetric.m	Thu Sep 29 09:46:52 2011 +0100
@@ -26,8 +26,7 @@
     % several decorrelation iterations might be needed to reach global
     % coherence mu. niter can be adjusted to needs.
     niter = 1;
-    while niter < 5 && ...
-            max(max(abs(dico'*dico -eye(length(dico))))) > mu + eps
+    while max(max(abs(dico'*dico -eye(length(dico))))) > mu + 0.01
         % find pairs of high correlation atoms
         colors = dico_color(dico, mu);
         
--- a/examples/SMALL_test_coherence.m	Tue Sep 20 15:52:33 2011 +0100
+++ b/examples/SMALL_test_coherence.m	Thu Sep 29 09:46:52 2011 +0100
@@ -18,8 +18,8 @@
 % Dependent parameters
 nActiveAtoms = fix(blockSize/100*percActiveAtoms); %number of active atoms
 epsilon      = 1/sigma; %error constraint for sparse representation step (corresponds to noise applied to signals)
-%minCoherence = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1))); %target coherence (based on coherence lower bound)
-minCoherence = 0.4;		%target coherence
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1))); %target coherence (based on coherence lower bound)
+minCoherence = [0.1:0.1:1];
 
 % Initial dictionaries
 dctDict = dictionary('dct',blockSize,dictSize);
@@ -29,7 +29,7 @@
 
 %% Generate audio denoising problem with low noise (audio representation)
 SMALL.Problem = generateAudioDenoiseProblem(signal.s,[],blockSize,...
-    dictSize,overlap,sigma); % generate representation problem
+	dictSize,overlap,sigma); % generate representation problem
 SMALL.Problem.b1 = SMALL.Problem.b; % copy signals from training set b to test set b1 (needed for later functions)
 
 % omp2 sparse representation solver
@@ -37,147 +37,144 @@
 solver = SMALL_init_solver('ompbox','omp2',ompParam,false); %solver structure
 
 
-%% Test ksvd dictionary update
+%% Test
 name = dicUpdate{1}; %use ksvd update
-SMALL.DL(1:9) = SMALL_init_DL(toolbox,name); %create dictionary learning structures
+SMALL.DL(1:36) = SMALL_init_DL(toolbox,name); %create dictionary learning structures
 
 % learn with random initialisation and no decorrelation
 SMALL.DL(1).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none'); %parameters for the dictionary learning
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','none'); %parameters for the dictionary learning
 SMALL.DL(1) = SMALL_learn(SMALL.Problem,SMALL.DL(1)); %learn dictionary
+%save('SMALL','SMALL');
 
 % learn with random initialisation and mailhe decorrelation
-SMALL.DL(2).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence); %parameters for the dictionary learning
-SMALL.DL(2) = SMALL_learn(SMALL.Problem,SMALL.DL(2)); %learn dictionary
-
-% learn with random initialisation and tropp decorrelation
-SMALL.DL(3).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence); %parameters for the dictionary learning
-SMALL.DL(3) = SMALL_learn(SMALL.Problem,SMALL.DL(3)); %learn dictionary
-
-% Learn with dct initialisation and no decorrelation
-SMALL.DL(4).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none','initdict',dctDict); %parameters for the dictionary learning
-SMALL.DL(4) = SMALL_learn(SMALL.Problem,SMALL.DL(4)); %learn dictionary
-
-% learn with dct initialisation and mailhe decorrelation
-SMALL.DL(5).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
-SMALL.DL(5) = SMALL_learn(SMALL.Problem,SMALL.DL(5)); %learn dictionary
-
-% learn with dct initialisation and tropp decorrelation
-SMALL.DL(6).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
-SMALL.DL(6) = SMALL_learn(SMALL.Problem,SMALL.DL(6)); %learn dictionary
-
-% Learn with gabor initialisation and no decorrelation
-SMALL.DL(7).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none','initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(7) = SMALL_learn(SMALL.Problem,SMALL.DL(7)); %learn dictionary
-
-% learn with gabor initialisation and mailhe decorrelation
-SMALL.DL(8).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(8) = SMALL_learn(SMALL.Problem,SMALL.DL(8)); %learn dictionary
-
-% learn with gabor initialisation and tropp decorrelation
-SMALL.DL(9).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(9) = SMALL_learn(SMALL.Problem,SMALL.DL(9)); %learn dictionary
-
-%% Test mailhe dictionary update
-name = dicUpdate{2}; %use mailhe update
-SMALL.DL(10:18) = SMALL_init_DL(toolbox,name); %create dictionary learning structure
-
-% learn with random initialisation and no decorrelation
-SMALL.DL(10).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none'); %parameters for the dictionary learning
-SMALL.DL(10) = SMALL_learn(SMALL.Problem,SMALL.DL(10)); %learn dictionary
-
-% learn with random initialisation and mailhe decorrelation
-SMALL.DL(11).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence); %parameters for the dictionary learning
-SMALL.DL(11) = SMALL_learn(SMALL.Problem,SMALL.DL(11)); %learn dictionary
+for iMu=1:10
+	SMALL.DL(1+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu)); %parameters for the dictionary learning
+	SMALL.DL(1+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(1+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
 
 % learn with random initialisation and tropp decorrelation
 SMALL.DL(12).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence); %parameters for the dictionary learning
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','tropp','coherence',minCoherence); %parameters for the dictionary learning
 SMALL.DL(12) = SMALL_learn(SMALL.Problem,SMALL.DL(12)); %learn dictionary
 
 % Learn with dct initialisation and no decorrelation
 SMALL.DL(13).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none','initdict',dctDict); %parameters for the dictionary learning
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','none','initdict',dctDict); %parameters for the dictionary learning
 SMALL.DL(13) = SMALL_learn(SMALL.Problem,SMALL.DL(13)); %learn dictionary
 
 % learn with dct initialisation and mailhe decorrelation
-SMALL.DL(14).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
-SMALL.DL(14) = SMALL_learn(SMALL.Problem,SMALL.DL(14)); %learn dictionary
+for iMu=1:10
+	SMALL.DL(13+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu),'initdict',dctDict); %parameters for the dictionary learning
+	SMALL.DL(13+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(13+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
 
 % learn with dct initialisation and tropp decorrelation
-SMALL.DL(15).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
-SMALL.DL(15) = SMALL_learn(SMALL.Problem,SMALL.DL(15)); %learn dictionary
+SMALL.DL(24).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','tropp','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(24) = SMALL_learn(SMALL.Problem,SMALL.DL(24)); %learn dictionary
+
 
 % Learn with gabor initialisation and no decorrelation
-SMALL.DL(16).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','none','initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(16) = SMALL_learn(SMALL.Problem,SMALL.DL(16)); %learn dictionary
+SMALL.DL(25).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','none','initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(25) = SMALL_learn(SMALL.Problem,SMALL.DL(25)); %learn dictionary
 
 % learn with gabor initialisation and mailhe decorrelation
-SMALL.DL(17).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','mailhe','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(17) = SMALL_learn(SMALL.Problem,SMALL.DL(17)); %learn dictionary
+for iMu=1:10
+	SMALL.DL(25+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu),'initdict',gaborDict); %parameters for the dictionary learning
+	SMALL.DL(25+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(25+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
 
 % learn with gabor initialisation and tropp decorrelation
-SMALL.DL(18).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
-    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
-    'decFcn','tropp','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
-SMALL.DL(18) = SMALL_learn(SMALL.Problem,SMALL.DL(18)); %learn dictionary
+SMALL.DL(36).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+	'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+	'decFcn','tropp','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(36) = SMALL_learn(SMALL.Problem,SMALL.DL(36)); %learn dictionary
 
 %% Evaluate coherence and snr of representation for the various methods
-sigNoiseRatio = zeros(18,1);
-mu = zeros(18,1);
-for i=1:18
-    SMALL.Problem.A = SMALL.DL(i).D;
-    tempSolver = SMALL_solve(SMALL.Problem,solver);
-    sigNoiseRatio(i) = snr(SMALL.Problem.b,SMALL.DL(i).D*tempSolver.solution);
-    dic(i) = dictionary(SMALL.DL(i).D);
-    mu(i) = dic(i).coherence;
+sigNoiseRatio = zeros(36,1);
+mu = zeros(36,1);
+for i=1:36
+	SMALL.Problem.A = SMALL.DL(i).D;
+	tempSolver = SMALL_solve(SMALL.Problem,solver);
+	sigNoiseRatio(i) = snr(SMALL.Problem.b,SMALL.DL(i).D*tempSolver.solution);
+	dic(i) = dictionary(SMALL.DL(i).D);
+	mu(i) = dic(i).coherence;
 end
 
 %% Plot results
 minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));
-maxSNR = max(sigNoiseRatio);
 
-figure, subplot(1,2,1)
-snrMat = buffer(sigNoiseRatio(1:9),3);
-bar(snrMat');
-title('Signal to noise ratio')
-xlabel('Initial dictionary')
-ylabel('SNR (dB)')
-set(gca,'XTickLabel',{'data','dct','gabor'});
-legend('none','Mailhe','Tropp')
-grid on
+figure,
+%subplot(3,1,1)
+hold on, grid on
+title('Data Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(1),sigNoiseRatio(1),'ks');
+plot(mu(12),sigNoiseRatio(12),'kd');
+plot(mu(2:11),sigNoiseRatio(2:11),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
 
+figure
+%subplot(3,1,2)
+hold on, grid on
+title('DCT Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(13),sigNoiseRatio(13),'ks');
+plot(mu(24),sigNoiseRatio(24),'kd');
+plot(mu(14:23),sigNoiseRatio(14:23),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
+
+figure
+%subplot(3,1,3)
+hold on, grid on
+title('Gabor Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(25),sigNoiseRatio(25),'ks');
+plot(mu(36),sigNoiseRatio(36),'kd');
+plot(mu(26:35),sigNoiseRatio(26:35),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
+
+% minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));
+% maxSNR = max(sigNoiseRatio);
+%
+% figure, subplot(2,2,1)
+% snrMat = buffer(sigNoiseRatio(1:9),3);
+% bar(snrMat');
+% title('Signal to noise ratio')
+% xlabel('Initial dictionary')
+% ylabel('SNR (dB)')
+% set(gca,'XTickLabel',{'data','dct','gabor'});
+% legend('none','Mailhe','Tropp')
+% grid on
+%
 % subplot(2,2,2), grid on
 % snrMat = buffer(sigNoiseRatio(10:18),3);
 % bar(snrMat');
@@ -187,18 +184,18 @@
 % set(gca,'XTickLabel',{'data','dct','gabor'},'YLim',[0 maxSNR+1]);
 % legend('none','Mailhe','Tropp')
 % grid on
-
-subplot(1,2,2), hold on, grid on
-title('Coherence')
-muMat = buffer(mu(1:9),3);
-line([0.5 3.5],[1 1],'Color','r');
-bar(muMat');
-line([0.5 3.5],[minMu minMu],'Color','k');
-set(gca,'XTick',1:3,'XTickLabel',{'data','dct','gabor'},'YLim',[0 1.05])
-legend('\mu_{max}','none','Mailhe','Tropp','\mu_{min}')
-ylabel('\mu')
-xlabel('Initial dictionary')
-
+%
+% subplot(2,2,3), hold on, grid on
+% title('Coherence')
+% muMat = buffer(mu(1:9),3);
+% line([0.5 3.5],[1 1],'Color','r');
+% bar(muMat');
+% line([0.5 3.5],[minMu minMu],'Color','k');
+% set(gca,'XTick',1:3,'XTickLabel',{'data','dct','gabor'},'YLim',[0 1.05])
+% legend('\mu_{max}','none','Mailhe','Tropp','\mu_{min}')
+% ylabel('\mu')
+% xlabel('Initial dictionary')
+%
 % subplot(2,2,4), hold on, grid on
 % title('Coherence - Mailhe Update')
 % muMat = buffer(mu(10:18),3);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/SMALL_test_coherence2.m	Thu Sep 29 09:46:52 2011 +0100
@@ -0,0 +1,119 @@
+clc, clear, close all
+
+%% Parameteres
+% Dictionary learning parameters
+toolbox   = 'TwoStepDL';							%dictionary learning toolbox
+dicUpdate = {'ksvd'};								%dictionary learning updates
+dicDecorr = {'none','mailhe','tropp','barchiesi'};	%dictionary decorrelation methods
+minCoherence = linspace(0.1,1,1);					%coherence levels
+minCoherence = 0.4;
+%dicDecorr = {'barchiesi'};
+
+iterNum   = 20;				%number of iterations
+epsilon   = 1e-6;			%tolerance level
+dictSize  = 512;			%number of atoms in the dictionary
+percActiveAtoms = 5;		%percentage of active atoms
+
+% Test signal parameters
+signal    = audio('music03_16kHz.wav'); %audio signal
+blockSize = 256;						%size of audio frames
+overlap   = 0.5;						%overlap between consecutive frames
+
+
+% Dependent parameters
+nActiveAtoms = fix(blockSize/100*percActiveAtoms); %number of active atoms
+
+% Initial dictionaries
+dctDict = dictionary('dct',blockSize,dictSize);
+dctDict = dctDict.phi;
+gaborParam = struct('N',blockSize,'redundancyFactor',2,'wd',@rectwin);
+gaborDict = Gabor_Dictionary(gaborParam);
+initDicts = {[],dctDict,gaborDict};
+initDicts = {[]};
+
+%% Generate audio approximation problem
+signal			 = buffer(signal,blockSize,blockSize*overlap,@rectwin);
+SMALL.Problem.b  = signal.S;
+SMALL.Problem.b1 = SMALL.Problem.b; % copy signals from training set b to test set b1 (needed for later functions)
+
+% omp2 sparse representation solver
+ompParam = struct('X',SMALL.Problem.b,'epsilon',epsilon,'maxatoms',nActiveAtoms); %parameters
+solver	 = SMALL_init_solver('ompbox','omp2',ompParam,false); %solver structure
+
+
+%% Test
+nDicUpdates = length(dicUpdate);		%number of dictionary updates
+nDecorrAlgs = length(dicDecorr);		%number of decorrelation algorithms
+nCorrLevels = length(minCoherence);		%number of coherence levels
+nInitDicts  = length(initDicts);		%number of initial dictionaries
+
+SMALL.DL(nInitDicts,nCorrLevels,nDecorrAlgs,nDicUpdates) = SMALL_init_DL(toolbox); %create dictionary learning structures
+for iInitDicts=1:nInitDicts
+	for iCorrLevels=1:nCorrLevels
+		for iDecorrAlgs=1:nDecorrAlgs
+			for iDicUpdates=1:nDicUpdates
+				SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).toolbox = toolbox;
+				SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).name = dicUpdate{iDicUpdates};
+				SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).profile = true;
+				SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).param = ...
+					struct( 'data',SMALL.Problem.b,...
+					'Tdata',nActiveAtoms,...
+					'dictsize',dictSize,...
+					'iternum',iterNum,...
+					'memusage','high',...
+					'solver',solver,...
+					'decFcn',dicDecorr{iDecorrAlgs},...
+					'coherence',minCoherence(iCorrLevels),...
+					'initdict',initDicts(iInitDicts));
+				
+				SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates) = ...
+					SMALL_learn(SMALL.Problem,SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates));
+			end
+		end
+	end
+end
+
+%% Evaluate coherence and snr of representation for the various methods
+sr = zeros(size(SMALL.DL));				%signal to noise ratio
+mu = zeros(size(SMALL.DL));				%coherence
+dic(size(SMALL.DL)) = dictionary;		%initialise dictionary objects
+for iInitDict=1:nInitDicts
+	for iCorrLevels=1:nCorrLevels
+		for iDecorrAlgs=1:nDecorrAlgs
+			for iDicUpdates=1:nDicUpdates
+				%Sparse representation
+				SMALL.Problem.A = SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).D;
+				tempSolver = SMALL_solve(SMALL.Problem,solver);
+				%calculate snr
+				sr(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates) = ...
+					snr(SMALL.Problem.b,SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).D*tempSolver.solution);
+				%calculate mu
+				dic(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates) = ...
+					dictionary(SMALL.DL(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).D);
+				mu(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates) = ...
+					dic(iInitDicts,iCorrLevels,iDecorrAlgs,iDicUpdates).coherence;
+			end
+		end
+	end
+end
+
+%% Plot results
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));	%lowe bound on coherence
+initDictsNames = {'Data','DCT','Gabor'};						
+dicDecorrNames = {'K-SVD','INK-SVD','Grassmannian','New'};
+lineStyles     = {'ks-','kd-','ko-','k*-'};
+for iInitDict=1:nInitDicts
+	figure, hold on, grid on
+	title([initDictsNames{iInitDict} ' Initialisation']);
+	plot([1 1],[0 25],'k-');
+	for iDecorrAlgs=1:nDecorrAlgs
+		plot(mu(iInitDicts,:,iDecorrAlgs,1),sr(iInitDicts,:,iDecorrAlgs,1),...
+			lineStyles{iDecorrAlgs});
+	end
+	plot([minMu minMu],[0 25],'k--')
+	
+	set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+	legend([{'\mu_{max}'},dicDecorrNames,{'\mu_{min}'}]);
+	xlabel('\mu');
+	ylabel('SNR (dB)');
+end
--- a/util/classes/dictionaryMatrices/grassmannian.m	Tue Sep 20 15:52:33 2011 +0100
+++ b/util/classes/dictionaryMatrices/grassmannian.m	Thu Sep 29 09:46:52 2011 +0100
@@ -48,7 +48,12 @@
 
 [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
+A = Uinit*Sigma_new*V_gram';				%update dictionary
+
+% param.step  = 0.01;
+% param.reg   = 0.01;
+% param.nIter = 20;
+% A = rotatematrix(initA,A,'linesearchlie',param);
 
 % %% Debug visualization function
 % function plotcart2d(A)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/iterativeprojections.m	Thu Sep 29 09:46:52 2011 +0100
@@ -0,0 +1,60 @@
+function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
+% 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
+% M. Elad, Sparse and Redundant Representations, Springer 2010.
+
+%% Parameters and Defaults
+error(nargchk(2,7,nargin));
+
+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
+
+%% Main algo
+A = normc(initA); %normalise columns
+[Uinit Sigma] = svd(A);
+G = A'*A; %gram matrix
+
+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');
+end
+
+% 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
+
+% [~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
+
+% A = normc(A);					%normalise dictionary
+
+[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
+
+% param.step  = 0.01;
+% param.reg   = 0.01;
+% param.nIter = 20;
+% A = rotatematrix(initA,A,'linesearchlie',param);
+
+% %% Debug visualization function
+% function plotcart2d(A)
+% compass(A(1,:),A(2,:));
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/rotatematrix.m	Thu Sep 29 09:46:52 2011 +0100
@@ -0,0 +1,96 @@
+function [Dhat cost W] = rotatematrix(D,Phi,method,param)
+%
+%
+%
+% REFERENCE
+% M.D. Plumbley, Geometrical Methods for Non-Negative ICA: Manifolds, Lie
+% Groups and Toral Subalgebra, Neurocomputing
+if ~nargin, testrotatematrix; return, end
+
+
+if ~exist('method','var') || isempty(method), method = 'unconstrained'; end
+
+J = @(W) 0.5*norm(D-W*Phi,'fro');
+cost = zeros(param.nIter,1);
+
+W   = eye(size(Phi,1));
+t = 0;
+Gprev = 0;
+Hprev = 0;
+for i=1:param.nIter
+	cost(i) = J(W);
+	grad = (W*Phi-D)*Phi';
+	switch method
+		case 'unconstrained'		% gradient descent
+			eta = param.step;
+			W = W - eta*grad;		% update W by steepest descent
+		case 'tangent'				% self correcting tangent
+			eta = param.step;
+			mu  = param.reg;
+			W = W - 0.5*eta*(grad - W*grad'*W + mu*W*(W'*W-eye(size(W))));
+		case 'steepestlie'
+			eta = param.step;
+			B = 2*skew(grad*W');	% calculate gradient in lie algebra
+			W = expm(-eta*B)*W;		% update W by steepest descent
+		case 'linesearchlie'
+			B = 2*skew(grad*W');	% calculate gradient in lie algebra
+			H = -B;					% calculate direction as negative gradient
+			t = searchline(J,H,W,t);% line search in one-parameter lie subalgebra
+			W = expm(t*H)*W;		% update W by line search
+		case 'conjgradlie'
+			G = 2*skew(grad*W');	% calculate gradient in lie algebra
+			H = -G + polakRibiere(G,Gprev)*Hprev; %calculate conjugate gradient direction
+			t = searchline(J,H,W,t);% line search in one-parameter lie subalgebra
+			W = expm(t*H)*W;		% update W by line search
+			Hprev = H;				% % save search direction
+			Gprev = G;				% % save gradient
+	end
+end
+Dhat = W*Phi;
+end
+% function C = matcomm(A,B)
+% %Matrix commutator
+% C = A*B-B*A;
+
+function gamma = polakRibiere(G,Gprev)
+gamma = G(:)'*(G(:)-Gprev(:))/(norm(Gprev(:))^2);
+if isnan(gamma) || isinf(gamma)
+	gamma = 0;
+end
+end
+
+function t = searchline(J,H,W,t)
+t = fminsearch(@(x) J(expm(x*H)*W),t);
+end
+
+function B = skew(A)
+B = 0.5*(A - A');
+end
+
+
+function testrotatematrix
+clear, clc, close all
+n = 256;
+m = 512;
+disp('A random matrix...');
+Phi = randn(n,m);
+disp('And its rotated mate...');
+Qtrue = expm(skew(randn(n)));
+D = Qtrue*Phi; 
+disp('Now, lets try to find the right rotation...');
+param.nIter = 1000;
+param.step  = 0.001;
+
+cost = zeros(param.nIter,4);
+[~, cost(:,1)] = rotatematrix(D,Phi,'unconstrained',param);
+[~, cost(:,2)] = rotatematrix(D,Phi,'steepestlie',param);
+[~, cost(:,3)] = rotatematrix(D,Phi,'linesearchlie',param);
+[~, cost(:,4)] = rotatematrix(D,Phi,'conjgradlie',param);
+
+figure, plot(cost)
+set(gca,'XScale','log','Yscale','log')
+legend({'uncons','settpestlie','linesearchlie','conjgradlie'})
+grid on
+xlabel('number of iterations')
+ylabel('J(W)')
+end
--- a/util/classes/util/regularisedictionary.m	Tue Sep 20 15:52:33 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-function Dreg = regularisedictionary(D)
-% REGULARISEDICTIONARY regularise dictionary for sparse representation
-%
-% Example
-% Dreg = regularisedictionary(D)
-%
-% Input
-% - D: initial dictionary matrix.
-%
-% Output
-% - Dreg: regularised dictionary
-%
-% References:
-%
-% See also
-
-% Author(s): Daniele Barchiesi
-% Copyright 2011-2011
-
-[n m] = size(D);	%dimension and number of atoms
-[U S V] = svd(D);	%SVD decomposition of the dictionary
-G = V*(S'*S)*V';	%gram matrix (equivalent to D'*D)
-I = eye(m);			%identity matrix
-lambda = 1;		    %regularisation coefficient
-
-% Optimise the gram matrix
-cvx_begin							%start cvx optimisation problem
-	variable Greg(m,m) symmetric;	%declare optimisation variable to be an mxm symmetric matrix
-	expression residual(m,m)		%declare residual as an intermediate calculation step
-	residual = G-Greg;				%define residual
-	minimise (max(max(abs(Greg-I))));	%define objective function
-	subject to						%declare constraints
-		Greg == semidefinite(m);	%positive semidefinite cone membership
-		diag(Greg) == ones(m,1);	%unit diagonal
-		norm(residual,'fro') <= lambda;	%
-cvx_end
-
-[~, Sgram Vgram] = svd(Greg);		%SVD decomposition of gramian
-Snew = sqrt(Sgram(1:n,:)).*sign(S); %calculate singular values of the regularised dictionary
-Dreg = U*Snew*Vgram';				%calculate regularised dictionary
-Dreg = normc(Dreg);					%normalise columns of regularised dictionary