# HG changeset patch # User Daniele Barchiesi # Date 1317286012 -3600 # Node ID 290cca7d346912ea544acb1ed287e83b8ac3c07f # Parent ff866a412be51b694331e450a7d783c742923e20 Added dictionary decorrelation functions and test script for ICASSP paper. diff -r ff866a412be5 -r 290cca7d3469 DL/two-step DL/SMALL_two_step_DL.m --- 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 diff -r ff866a412be5 -r 290cca7d3469 DL/two-step DL/dico_decorr_symetric.m --- 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); diff -r ff866a412be5 -r 290cca7d3469 examples/SMALL_test_coherence.m --- 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); diff -r ff866a412be5 -r 290cca7d3469 examples/SMALL_test_coherence2.m --- /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 diff -r ff866a412be5 -r 290cca7d3469 util/classes/dictionaryMatrices/grassmannian.m --- 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) diff -r ff866a412be5 -r 290cca7d3469 util/classes/dictionaryMatrices/iterativeprojections.m --- /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,:)); diff -r ff866a412be5 -r 290cca7d3469 util/classes/dictionaryMatrices/rotatematrix.m --- /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 diff -r ff866a412be5 -r 290cca7d3469 util/classes/util/regularisedictionary.m --- 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