Mercurial > hg > smallbox
changeset 155:b14209313ba4 ivand_dev
Integration of Majorization Minimisation Dictionary Learning
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/AudioDicoLearning/AudioDicoLearning_Demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,9 @@ +M = '256'; +N = '512'; +L = '8192'; % N*16; +IT = '100'; % For a full convergence at least '10000' iterations; +Lambda = '.01'; % Lagrangian multiplier +admiss = 'un'; % admissible set for dictionaries 'bn' = bounded norm, 'un' = unit norm +tic +DLMM_Audio(M,N,L,IT,Lambda,admiss) +toc \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/AudioDicoLearning/DLMM_Audio.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,39 @@ +function [] = DLMM_Audio(m,n,l,it,Lambda,admiss) +M = str2num(m); +N = str2num(n); +L = str2num(l); +IT = str2num(it); +lambda = str2num(Lambda); +load(['Param',num2str(M),'X',num2str(N),'X',num2str(L),'kADCT2X01.mat']) +method = ['bn';'un']; +res = 1; +while method(res)~=admiss + res = res+1; +end +[PhiN,PhiM] = size(Phio); +Phi = randn(M,N); +Phi = Phi*diag(1./sqrt(sum(Phi.^2))); +[PhiN,L] = size(x); +phim = PhiM; +unhat = zeros(PhiM,L); +maxIST = 10; +maxIDU = 10; +it = 1; +ert(1) = norm(Phi*unhat-x,'fro')^2 + lambda*(sum(sum(abs(unhat)))); +while it <= IT + it +% %%%%%% Iterative Soft-Thresholding for Sparse Approximation %%%%% + tau2 = .1+(svds(Phi,1))^2; + eps = 3*10^-4; + map = 0; % Projecting on the selected space (0=no,1=yes) + [unhat,er] = mm1(Phi,x,unhat,tau2,lambda,maxIST,eps,map); + er + %%%%%% Majorisation Minimisation for Dictionary Update %%%%%% + c = .1 + svds(unhat,1)^2; + eps = 10^-7; + [Phi,unhat] = dict_update_REG_cn(Phi,x,unhat,maxIDU,eps,'un'); + %%%%% + ert(it+1) = norm(Phi*unhat-x,'fro')^2 + lambda*(sum(sum(abs(unhat)))); + it = it+1; +end +save(['DLMMl1',num2str(PhiN),'X',num2str(PhiM),'X',num2str(L),'iK',method(res,:),'Lambda',num2str(lambda),'Audio.mat'],'Phi','x','unhat','ert','lambda') \ No newline at end of file
Binary file DL/Majorization Minimization DL/AudioDicoLearning/Param256X512X8192kADCT2X01.mat has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/Demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,45 @@ +clear +M = 20; % Signal length +N = 40; % Coefficient Space Dimension +L = 32*N; % Number of Training Samples +R = 3; % Sparsity +IT = 1000; % Number of alternating sparse approximation and dictionary update +map = 1; % Debiasing. 0 = No, 1 = Yes +maxIT = 1000; % Inner-loop maximum iteration number. +lambda = 2*.2; % Lagrangian multiplier. +epsx = 10^-7; % Stopping criterion for iterative softthresholding +epsd = 10^-7; % Stopping criterion for MM dictionary update +cvset = 0; % Dictionary constraint. 0 = Non convex ||d|| = 1, 1 = Convex ||d||<=1 +Tre = .99; % Threshold for accepting too atoms identical +%%%% Generative Dictionaries +Do = randn(M,N); % Generative Dictionary +Do = Do*(diag(sum((Do'*Do).*eye(length(Do))).^-.5)); % Normalization +%%%% Sparse signal generation %%%%% +Xo = zeros(N,L); % Original Sparse Coefficients +for l = 1:L + r = 1; + while r<=R + ind = fix(rand(1)*N)+ones(1); + a = rand(1); + if Xo(ind)==0 + Xo(ind,l) = (.8*rand(1)+.2)*((a>=.5)-(a<.5)); + r = r+1; + end + end +end +Y = Do*Xo; % Sparse Signals +%%%% Algorithm initialization +D = randn(M,N); % Initial Dictionary +D = D*(diag(sum((D'*D).*eye(length(D))).^-.5)); % Normalization +X = ones(size(Xo)); % Initial coefficients +for it = 1:IT, + it + to = .1+svds(D,1); + [X,cost(it)] = mm1(D,Y,X,to,lambda,maxIT,epsx,map); + plot(cost); + [D,X] = dict_update_REG_cn(D,Y,X,maxIT,epsd,cvset); +end +%%% +success = sum(max(abs((Do'*D)))>=Tre); +display([' ------------------']) +display([' ',num2str(success),'% of the atoms successfully recovered after ',num2str(IT),' iterations.']); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/ExactRec_Demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,49 @@ +clear +IT = 1000; % Number of iterations +R = 3:7; % Sparsity of the signals +SN = 5; % Number of trials +M = 20; % Signal space dimension +N = 40; % Ambient space dimension +L = 32*N; % Number of training signals + +ParamGen(M,N,L,R,SN) % Generating generic problems + +%%%% Dictionary recovery algorithms +for sn = 1:SN + for r = R + mmdlcn_exactRec_demo(IT,r,sn,'un') + mapcn_exactRec_demo(IT,r,sn,'bn') + ksvd_exactRec_demo(IT,r,sn) + modcn_exactRec_demo(IT,r,sn) + end +end +%%%%%%% +Tre = .99; +for k = R + for i=1:SN + load(['RDLl1',num2str(M),'t',num2str(IT),'ikiun',num2str(k),'v1d',num2str(i),'.mat']) + nrRDLu(i,k-R(1)+1) = sum(max(abs((Phid'*Phi)))>=Tre); + %%%%%% + load(['MAPl1',num2str(M),'t',num2str(IT),'ikiun',num2str(k),'v1d',num2str(i),'.mat']) + nrMAP(i,k-R(1)+1) = sum(max(abs((Phid'*Phi)))>=Tre); + %%%%%% + load(['KSVDl1',num2str(M),'t',num2str(IT),'ikiun',num2str(k),'v1d',num2str(i),'.mat']) + nrKSVD(i,k-R(1)+1) = sum(max(abs((Phid'*Phi)))>=Tre); + %%%%%% + load(['MODl1',num2str(M),'t',num2str(IT),'ikiun',num2str(k),'v1d',num2str(i),'.mat']) + nrMOD(i,k-R(1)+1) = sum(max(abs((Phid'*Phi)))>=Tre); + end +end +clf +errorbar(R+.01,10*mean(nrRDLu,1)/4,std(nrRDLu,0,1),'k-.') +hold on +errorbar(R-.01,10*mean(nrKSVD,1)/4,std(nrKSVD,0,1),'r*-') +errorbar(R+.01,10*mean(nrMOD,1)/4,std(nrMOD,0,1),'b-v') +errorbar(R-.01,10*mean(nrMAP,1)/4,std(nrMAP,0,1),'b-^') +axis([2.5 7.5 15 105]); +title('Constrained column-norm dictionaries') +xlabel('Sparsity (# of non-zero elements in each coefficient vector)') +ylabel(['Average percents of exact recovery',sprintf('\n'),'after ',num2str(IT),' iterations']) +grid on +legend('MM Unit norm','K-SVD','MOD','MAPbased-DL'); +axis([2.8 7.2 -5 105]) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/dict_update_KSVD_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,22 @@ +function[Phiout,unhatnz] = dict_update_KSVD_cn(Phi,x,unhat) +%% K-SVD Dictionary Learning with the constraint on the column norms %%%%% + +%% +unhatn = unhat; +rPerm = randperm(size(Phi,2)); +%% for l = 1:size(Phi,2), +for l = rPerm + unhatns = unhat; + unhatns(l,:) = zeros(1,size(unhat,2)); + E = x-Phi*unhatns; + ER = E(:,(abs(unhat(l,:))>0)); + [U,S,V] = svd(ER,0); + Phi(:,l) = U(:,1) + unhatn(l,:) = zeros(1,size(unhat,2)); + unhatn(l,(abs(unhat(l,:))>0)) = S(1,1)*V(:,1); + unhat = unhatn; +end +%% +unhatnz = unhat; +Phiout = Phi; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/dict_update_MAP_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,27 @@ +function [Phiout,unhatnz] = dict_update_MAP_cn(Phi,x,unhat,mu,maxIT,eps,cvset) +%% Maximum A Posteriori Dictionary Update with the constraint on the column norms %%%%% + + + +%% + +K = Phi; +B = zeros(size(Phi,1),size(Phi,2)); +i = 1; +%% +while (sum(sum((B-K).^2))>eps)&&(i<=maxIT) + B = K; + for j = 1:size(K,2), + K(:,j) = K(:,j) - mu*(eye(size(K,1))-K(:,j)*K(:,j)')*(K*unhat-x)*unhat(j,:)'; + end + i = i+1; +end + +%% depleted atoms cancellation %%% + +[Y,I] = sort(sum(K.^2),'descend'); +RR = sum(Y>=.01); +Phiout = K(:,I(1:RR)); +unhatnz = unhat(I(1:RR),:); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/dict_update_MOD_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,20 @@ +function [Phiout,unhatnz] = dict_update_MOD_cn(x,unhat,cvset) +%% Method of Optimized Direction (MOD) Dictionary Update followed by the +%% constraint on the column norms %% + +%% +eps = .01; +K = x*unhat'*inv(unhat*unhat'+eps*eye(size(unhat,1))); +if cvset == 1, + K = K*diag(min(sum(K.^2).^(-.5),1)); % Projecting on the convex constraint set +else + K = K*diag(sum(K.^2).^(-.5)); % Projecting on the Unit-Norm constraint set +end + +%% depleted atoms cancellation %%% +[Y,I] = sort(sum(K.^2),'descend'); +RR = sum(Y>=.01); +Phiout = K(:,I(1:RR)); +unhatnz = unhat(I(1:RR),:); + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/ksvd_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,25 @@ +% K-SVD algorithm for Dictionary Learning +% Y = input data (M X L matrix) +% Phi = initial dictionary (M X N), e.g. random dictionary or first N data samples +% lambda = regularization coefficient (||Phi*X-Y||_F)^2 + lambda*||X||_1 +% IT = number of iterations +function [Phiout,X,ert] = ksvd_cn(Y,Phi,lambda,IT) +maxIT = 1000; +[PhiN,PhiM] = size(Phi); +RR1 = PhiM; +%%%%%%%%%%%%%% +% [PhiM,L] = size(ud); +[PhiN,L] = size(Y); +X = ones(PhiM,L); +for it = 1:IT + to = .1+svds(Phi,1); + [PhiN,PhiM] = size(Phi); + %%%% + eps = 3*10^-4; + map = 1; % Projecting on the selected space (0=no,1=yes) + [X,l1err] = mm1(Phi,Y,X,to,lambda,maxIT,eps,map); %% Sparse approximation with Iterative Soft-thresholding + ert(it) = l1err; + %%% + [Phi,X] = dict_update_KSVD_cn(Phi,Y,X); +end +Phiout = Phi; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/ksvd_exactRec_demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,20 @@ +function ksvd_exactRec_demo(it,k,sn) +tic +IT = it; +K = k; +SN = sn; +if SN<10, + samnum = ['0',num2str(SN)]; +else + samnum = num2str(SN); +end +load(['Param',num2str(K),'kS',samnum,'.mat']) +lambda = 2*.2; % 2 * Smallest coefficients (Soft Thresholding) +%%%%%%%%% +Phi = Phio; +M = size(Phi,1); +[Phi,unhat,ert] = ksvd_cn(x,Phi,lambda,IT); + +save(['KSVDl1',num2str(M),'t',num2str(IT),'ikiun',num2str(K),'v1d',num2str(SN),'.mat'],'Phi','Phid','x','ud','unhat','ert') + +toc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mapcn_exactRec_demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,29 @@ +function mapcn_exactRec_demo(it,k,sn,cs) +tic +IT = it; +K = k; +SN = sn; +if SN<10, + samnum = ['0',num2str(SN)]; +else + samnum = num2str(SN); +end +load(['Param',num2str(K),'kS',samnum,'.mat']) +lambda = 2*.2; % 2 * Smallest coefficients (Soft Thresholding) +if cs == 'bn' + res = 1; +elseif cs == 'un' + res = 2; +else + disp('Undefined dictionary admissible set.') +end +method = ['bn';'un']; +res = 2; % 1 for bounded-norm, 2 for unit-norm +%%%%%%%%%%%%%% +Phi = Phio; +M = size(Phi,1); +[Phi,unhat,ert] = mapdl_cn(x,Phi,lambda,IT,res); + +save(['MAPl1',num2str(M),'t',num2str(IT),'iki',method(res,:),num2str(K),'v1d',num2str(SN),'.mat'],'Phi','Phid','x','ud','unhat','ert') + +toc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mapdl_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,28 @@ +% Maximum A Posteriori Estimation for Dictionary Learning +% Y = input data (M X L matrix) +% Phi = initial dictionary (M X N), e.g. random dictionary or first N data samples +% lambda = regularization coefficient (||Phi*X-Y||_F)^2 + lambda*||X||_1 +% IT + number of iterations +% res = dictionary constraint. 'un' = unit colomn norm, 'bn' = bounded colomn norm +function [Phiout,X,ert] = mapdl_cn(Y,Phi,lambda,IT,res) +maxIT = 1000; +[PhiN,PhiM] = size(Phi); +RR1 = PhiM; +%%%%%%%%%%%%%% +% [PhiM,L] = size(ud); +[PhiN,L] = size(Y); +X = ones(PhiM,L); +for it = 1:IT + to = .1+svds(Phi,1); + [PhiN,PhiM] = size(Phi); + %%%% + eps = 3*10^-4; + map = 1; % Projecting on the selected space (0=no,1=yes) + [X,l1err] = mm1(Phi,Y,X,to,lambda,maxIT,eps,map); %% Sparse approximation with Iterative Soft-thresholding + ert(it) = l1err; + %%% + eps = 10^-7; + mu = 10^-4; + [Phi,X] = dict_update_MAP_cn(Phi,Y,X,mu,maxIT,eps,res); +end +Phiout = Phi; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mapfn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,16 @@ +%% Maximum A Posteriori Dictionary Learning with the constraint on the column norms %%%%% +function [Phiout,unhatnz] = mapfn(Phi,x,unhat,mu,maxIT,eps,phim,res) +K = Phi; +B = zeros(size(Phi,1),size(Phi,2)); +i = 1; +while (sum(sum((B-K).^2))>eps)&&(i<=maxIT) + B = K; + E = x-K*unhat; + K = K+mu*(phim*E*unhat'-trace(unhat*E'*K)*K); + i = i+1; +end +%%% depleted atoms cancellation %%% +[Y,I] = sort(sum(K.^2),'descend'); +RR = phim; +Phiout = K(:,I(1:RR)); +unhatnz = unhat(I(1:RR),:); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mmdl_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,26 @@ +% Majorization Minimization for Dictionary Learning +% Y = input data (M X L matrix) +% Phi = initial dictionary (M X N), e.g. random dictionary or first N data samples +% lambda = regularization coefficient (||Phi*X-Y||_F)^2 + lambda*||X||_1 +% IT = number of iterations +% res = dictionary constraint. 'un' = unit colomn norm, 'bn' = bounded colomn norm +function [Phiout,X,ert] = mmdl_cn(Y,Phi,lambda,IT,res) +maxIT = 1000; +[PhiN,PhiM] = size(Phi); +RR1 = PhiM; +%%%%%%%%%%%%%% +[PhiN,L] = size(Y); +X = ones(PhiM,L); +for it = 1:IT + to = .1+svds(Phi,1); + [PhiN,PhiM] = size(Phi); + %%%% + eps = 3*10^-4; + map = 1; % Projecting on the selected space (0=no,1=yes) + [X,l1err] = mm1(Phi,Y,X,to,lambda,maxIT,eps,map); %% Sparse approximation with Iterative Soft-thresholding + ert(it) = l1err; + %%% + eps = 10^-7; + [Phi,X] = dict_update_REG_cn(Phi,Y,X,maxIT,eps,res); +end +Phiout = Phi; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mmdlcn_exactRec_demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,29 @@ +function mmdlcn_exactRec_demo(it,k,sn,cs) +tic +IT = it; +K = k; +SN = sn; +if SN<10, + samnum = ['0',num2str(SN)]; +else + samnum = num2str(SN); +end +load(['Param',num2str(K),'kS',samnum,'.mat']) +lambda = 2*.2; % 2 * Smallest coefficients (Soft Thresholding) +if cs == 'bn' + res = 1; +elseif cs == 'un' + res = 2; +else + disp('Undefined dictionary admissible set.') +end +method = ['bn';'un']; +res = 2; % 1 for bounded-norm, 2 for unit-norm +%%%%%%%%%%%%%% +Phi = Phio; +M = size(Phi,1); +[Phi,unhat,ert] = mmdl_cn(x,Phi,lambda,IT,res); + +save(['RDLl1',num2str(M),'t',num2str(IT),'iki',method(res,:),num2str(K),'v1d',num2str(SN),'.mat'],'Phi','Phid','x','ud','unhat','ert') + +toc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mod_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,26 @@ +% MOD method for Dictionary Learning +% Y = input data (M X L matrix) +% Phi = initial dictionary (M X N), e.g. random dictionary or first N data samples +% lambda = regularization coefficient (||Phi*X-Y||_F)^2 + lambda*||X||_1 +% IT = number of iterations +% res = dictionary constraint. 'un' = unit colomn norm, 'bn' = bounded colomn norm +function [Phiout,X,ert] = mod_cn(Y,Phi,lambda,IT) +maxIT = 1000; +[PhiN,PhiM] = size(Phi); +RR1 = PhiM; +%%%%%%%%%%%%%% +[PhiN,L] = size(Y); +X = ones(PhiM,L); +for it = 1:IT + to = .1+svds(Phi,1); + [PhiN,PhiM] = size(Phi); + %%%% + eps = 3*10^-4; + map = 1; % Projecting on the selected space (0=no,1=yes) + [X,l1err] = mm1(Phi,Y,X,to,lambda,maxIT,eps,map); %% Sparse approximation with Iterative Soft-thresholding + ert(it) = l1err; + %%% + res = 2; % Fixed column norm dictionary constraint. + [Phi,X] = dict_update_MOD_cn(Y,X,res); +end +Phiout = Phi; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/mod_exactRec.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,39 @@ +%%% MOD (||ki||<=1) %%%% +function mod_exactRec(it,k,sn) +tic +IT = str2num(it); +K = str2num(k); +SN = str2num(sn); +if SN<10, + samnum = ['0',num2str(SN)]; +else + samnum = num2str(SN); +end +load(['Param',num2str(K),'kS',samnum,'.mat']) +method = ['bn';'un']; +res = 2; % 1 for bounded-norm, 2 for unit-norm +lambda = 2*.2; % 2 * Smallest coefficients (Soft Thresholding) +% lambda = 2*.2^2; % 2 * Smallest coefficients (Hard Thresholding) +%%%%%%%%%%%%%% +Phi = Phio; +[PhiN,PhiM] = size(Phi); +RR1 = PhiM; +%%%%%%%%%%%%%% +[PhiM,L] = size(ud); +unhat = ones(PhiM,L); +for it = 1:IT + it + to = .1+svds(Phi,1); + [PhiN,PhiM] = size(Phi); + %%%% +% eps = 10^-7; + eps = 3*10^-4; + maxIT = 1000; + map = 0; + [unhat,l1err] = mm1(Phi,x,unhat,to,lambda,maxIT,eps,map); %% Sparse approximation with Iterative Soft-thresholding + ert(it) = l1err; + %%% + [Phi,unhat] = modcn(x,unhat,res); +end +save(['MODl120t',num2str(IT),'iki',method(res,:),num2str(K),'v2d',num2str(SN),'.mat'],'Phi','Phid','x','ud','unhat','ert') +toc \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/ExactDicoRecovery/modcn_exactRec_demo.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,20 @@ +function modcn_exactRec_demo(it,k,sn) +tic +IT = it; +K = k; +SN = sn; +if SN<10, + samnum = ['0',num2str(SN)]; +else + samnum = num2str(SN); +end +load(['Param',num2str(K),'kS',samnum,'.mat']) +lambda = 2*.2; % 2 * Smallest coefficients (Soft Thresholding) +%%%%%%%%%%%%%% +Phi = Phio; +M = size(Phi,1); +[Phi,unhat,ert] = mod_cn(x,Phi,lambda,IT); + +save(['MODl1',num2str(M),'t',num2str(IT),'ikiun',num2str(K),'v1d',num2str(SN),'.mat'],'Phi','Phid','x','ud','unhat','ert') + +toc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/README Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,31 @@ +Majorization Minimization Dictionary Learning (MMDL) Package (Beta version) + +This package includes the codes needed to generate some of the plots in the paper, "Dictionary Learning for Sparse Approximations with the Majorization Method" by M. Yaghoobi, T. Blumensath and M. Davies, IEEE Transaction on Signal Processing, Vol. 57, No. 6, pp 2178-2191, 2009. + +The package also includes the authors' implementations of some other well-known dictionary learning methods. These have been used for comparison reason. + +The package includes two parts: + +1- Building blocks of the MMDL: + + a) mm1 : Iterative Softthresholding for sparse approximation + b) rdlcn : MMDL with column norm constraint + c) rdlfn : MMDL with Frobenius norm constraint + d) softthreshold : Softthresholding operator + e) Demo : A demo to demonstrate how to use MMDL and recover a generic dictionary + +2- MMDL in action: + a) ExactDicoRecovery : Exact generic atom recovery experiment. It includes a comparison with other methods + b) Dictionary Learning for Audio : Dictionary learning with MMDL method applied to some classic music recorded from BBC radio 3 + +Installation: + +You just need to decompress the package into your local disk and add the path of main directory to the Matlab path. + +Citation: + +This package is free to use. However, any publication based on using this package, in any kind, should cite the mentioned paper. + +Copyright: + +The use of this package does not need any farther permission, if it is restricted to the academic and industrial researches. Any military affiliated use of this package is strictly prohibited.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/dict_update_REG_cn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,48 @@ +function [Phiout,unhatnz] = dict_update_REG_cn(Phi,x,unhat,maxIT,eps,cvset) +%% Regularized Dictionary Learning with the constraint on the column norms %%%%% +% Phi = Normalized Initial Dictionary +% x = Signal(x). This can be a vector or a matrix +% unhat = Initial guess for the coefficients +% to = 1/(step size) . It is larger than spectral norm of coefficient matrix x +% eps = Stopping criterion for iterative softthresholding and MM dictionary update +% cvset = Dictionary constraint. 0 = Non convex ||d|| = 1, 1 = Convex ||d||<=1 +% Phiout = Updated Dictionary +% unhatnz Updated Coefficients (the same as input in this version) + +%% +B = Phi; +K = zeros(size(Phi,1),size(Phi,2)); +c = .1 + svds(unhat,1)^2; %.1 +c3 = (1/c)*eye(size(B,2)); +c1 = x*unhat'*c3; +c2 = (c*eye(size(B,2))-unhat*unhat')*c3; + +%% + +for i=1:maxIT +% if i>1 +% B = K; +% end + + K = c1 + B*c2; + + if cvset == 1, + K = K*diag(min(sum(K.^2).^(-.5),1)); % with convex constraint set + else + % Mehrdad original - + % K = K*diag(sum(K.^2).^(-.5)); % with fixed-norm constraint set + K = normc(K); + end + + if (sum(sum((B-K).^2)) < eps) + break; + end + + B = K; +end +%% depleted atoms cancellation %%% +[Y,I] = sort(sum(K.^2),'descend'); +RR = sum(Y>=.01); +Phiout = K(:,I(1:RR)); +unhatnz = unhat(I(1:RR),:); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/dict_update_REG_fn.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,39 @@ +function [Phiout,unhatnz] = dict_update_REG_fn(Phi,x,unhat,maxIT,eps,cvset) +%% Regularized Dictionary Learning with the constraint on the matrix frobenius norms %%%%% +% Phi = Normalized Initial Dictionary +% x = Signal(x). This can be a vector or a matrix +% unhat = Initial guess for the coefficients +% to = 1/(step size) . It is larger than spectral norm of coefficient matrix x +% eps = Stopping criterion for iterative softthresholding and MM dictionary update +% cvset = Dictionary constraint. 0 = Non convex ||D|| = N, 1 = Convex ||D||<=N +% Phiout = Updated Dictionary +% unhatnz Updated Coefficients (the same as input in this version) + +%% +B = Phi; +phim = norm(Phi, 'fro'); +K = zeros(size(Phi,1),size(Phi,2)); +c = .1 + svds(unhat,1)^2; + +%% +i = 1; +while (sum(sum((B-K).^2)) > eps)&&(i<=maxIT) + if i>1 + B = K; + end + K = 1/c *(x*unhat' + B*(c*eye(size(B,2))-unhat*unhat')); + Kfn = sum(sum(K.^2)); + if cvset == 1, + K = min(1,phim/Kfn)*K; % with convex constraint set + else + K = (phim/Kfn)*K; % with fixed-norm constraint set + end + i = i+1; +end + +%% depleted atoms cancellation %%% +[Y,I] = sort(sum(K.^2),'descend'); +RR = sum(Y>=0.0001); +Phiout = K(:,I(1:RR)); +unhatnz = unhat(I(1:RR),:); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/mm1.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,50 @@ +function [unhat,er] = mm1(Phi,x,u0,to,lambda,maxIT,eps,map) +%% Iterative Soft Thresholding (with optional debiasing) +% +% Phi = Normalized Dictionary +% x = Signal(x). This can be a vector or a matrix +% u0 = Initial guess for the coefficients +% to = 1/(step size) . It is larger than spectral norm of dictionary Phi +% lambda = Lagrangian multiplier. (regulates shrinkage) +% eps = Stopping criterion for iterative softthresholding and MM dictionary update +% map = Debiasing. 0 = No, 1 = Yes +% unhat = Updated coefficients +% er = Objective cost +%% +cont = 1; +in = 1; +% un = zeros(size(u0,1),size(u0,2)); +un = u0; +c1 = (1/to^2)*Phi'*x; +c2 = (1/to^2)*(Phi'*Phi); +%%%% +while (cont && (in<=maxIT)) + unold = un; + %%%%%% Soft Thresholding %%%%%%% + alphap = (un + c1 - c2*un); + un = (alphap-(lambda/(2*to^2))*sign(alphap)).*(abs(alphap)>=(lambda/(2*to^2))); + in = in+1; + cont = sum(sum((unold-un).^2))>eps; +end +%%%%%%%%%% +if map == 1, + %% Mapping on the selected space %%%% + [uN,uM] = size(un); + unhat = zeros(uN,uM); + for l = 1:uM, + unz = (abs(un(:,l))>0); + M = diag(unz); + PhiNew = Phi*M; + PhiS = PhiNew(:,unz); + unt = inv(PhiS'*PhiS+.0001*eye(sum(unz)))*PhiS'*x(:,l); + unhat(unz,l) = unt; + end +else + unhat = un; +end +%%% Cost function calculation +if map == 1, + er = sum(sum((Phi*unhat-x).^2))+lambda*(sum(sum(abs(unhat)>0))); %% l_0 Cost function +else + er = sum(sum((Phi*unhat-x).^2))+lambda*(sum(sum(abs(unhat)))); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/softthreshold.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,6 @@ +% Soft threshold operator for matrices +% [Y] = softthreshold(X,theta) +% X : input matrix +% theta : thresholding parameter +function [Y] = softthreshold(X,theta) +Y = X - min(theta,abs(X)).*sign(X); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/wrapper_mm_DL.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,185 @@ +function DL = wrapper_mm_DL(Problem, DL) + +% determine which solver is used for sparse representation % + +solver = DL.param.solver; + +% determine which type of udate to use +% (Mehrdad Yaghoobi implementations: 'MM_cn', MM_fn', 'MOD_cn', +% 'MAP_cn', 'KSVD_cn') + +typeUpdate = DL.name; + +sig = Problem.b; + +% determine dictionary size % + +if (isfield(DL.param,'initdict')) + if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:)))) + dictsize = length(DL.param.initdict); + else + dictsize = size(DL.param.initdict,2); + end +end + +if (isfield(DL.param,'dictsize')) % this superceedes the size determined by initdict + dictsize = DL.param.dictsize; +end + +if (size(sig,2) < dictsize) + error('Number of training signals is smaller than number of atoms to train'); +end + + +% initialize the dictionary % + +if (isfield(DL.param,'initdict')) + if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:)))) + dico = sig(:,DL.param.initdict(1:dictsize)); + else + if (size(DL.param.initdict,1)~=size(sig,1) || size(DL.param.initdict,2)<dictsize) + error('Invalid initial dictionary'); + end + dico = DL.param.initdict(:,1:dictsize); + end +else + data_ids = find(colnorms_squared(sig) > 1e-6); % ensure no zero data elements are chosen + perm = randperm(length(data_ids)); + dico = sig(:,data_ids(perm(1:dictsize))); +end + + +% number of iterations (default is 40) % + +if isfield(DL.param,'iternum') + iternum = DL.param.iternum; +else + iternum = 40; +end + +% number of iterations (default is 40) % + +if isfield(DL.param,'iterDictUpdate') + maxIT = DL.param.iterDictUpdate; +else + maxIT = 1000; +end + +% Stopping criterion for MM dictionary update (default = 1e-7) + +if isfield(DL.param,'epsDictUpdate') + epsD = DL.param.epsDictUpdate; +else + epsD = 1e-7; +end + +% Dictionary constraint - 0 = Non convex ||d|| = 1, 1 = Convex ||d||<=1 +% (default cvset is o) % + +if isfield(DL.param,'cvset') + cvset = DL.param.cvset; +else + cvset = 0; +end + +% determine if we should do decorrelation in every iteration % + +if isfield(DL.param,'coherence') + decorrelate = 1; + mu = DL.param.coherence; +else + decorrelate = 0; +end + +% show dictonary every specified number of iterations + +if isfield(DL.param,'show_dict') + show_dictionary = 1; + show_iter = DL.param.show_dict; +else + show_dictionary = 0; + show_iter = 0; +end + +% This is a small patch that needs to be resolved in dictionary learning we +% want sparse representation of training set, and in Problem.b1 in this +% version of software we store the signal that needs to be represented +% (for example the whole image) +if isfield(Problem,'b1') + tmpTraining = Problem.b1; + Problem.b1 = sig; +end +if isfield(Problem,'reconstruct') + Problem = rmfield(Problem, 'reconstruct'); +end +solver.profile = 0; + +% main loop % + +for i = 1:iternum + Problem.A = dico; + + solver = SMALL_solve(Problem, solver); + + switch lower(typeUpdate) + case 'mm_cn' + [dico, solver.solution] = ... + dict_update_REG_cn(dico, sig, solver.solution, maxIT, epsD, cvset); + case 'mm_fn' + [dico, solver.solution] = ... + dict_update_REG_fn(dico, sig, solver.solution, maxIT, epsD, cvset); + case 'mod_cn' + [dico, solver.solution] = dict_update_MOD_cn(sig, solver.solution, cvset); + case 'map_cn' + if isfield(DL.param,'muMAP') + muMAP = DL.param.muMAP; + else + muMAP = 1e-4; + end + [dico, solver.solution] = ... + dict_update_MAP_cn(dico, sig, solver.solution, muMAP, maxIT, epsD, cvset); + case 'ksvd_cn' + [dico, solver.solution] = dict_update_KSVD_cn(dico, sig, solver.solution); + otherwise + error('Dictionary update is not defined'); + end + + % Set previous solution as the best initial guess + % for the next iteration of iterative soft tresholding + + if (strcmpi(solver.toolbox, 'MMbox')) + solver.param.initcoeff = solver.solution; + end + + % Optional decorrelation of athoms - this is from Boris Mailhe and + % we need to test how it preforms with Mehrdad's updates + + if (decorrelate) + dico = dico_decorr(dico, mu, solver.solution); + end + + if ((show_dictionary)&&(mod(i,show_iter)==0)) + dictimg = SMALL_showdict(dico,[8 8],... + round(sqrt(size(dico,2))),round(sqrt(size(dico,2))),'lines','highcontrast'); + figure(2); imagesc(dictimg);colormap(gray);axis off; axis image; + pause(0.02); + end +end +if isfield(Problem,'b1') + Problem.b1 = tmpTraining; +end +DL.D = dico; + +end + +function Y = colnorms_squared(X) + +% compute in blocks to conserve memory +Y = zeros(1,size(X,2)); +blocksize = 2000; +for i = 1:blocksize:size(X,2) + blockids = i : min(i+blocksize-1,size(X,2)); + Y(blockids) = sum(X(:,blockids).^2); +end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/wrapper_mm_solver.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,74 @@ +function [X , cost] = wrapper_mm_solver(b, A, param) +%% SMALL wrapper for Majorization Maximization toolbos solver +% +% Function gets as input +% b - measurement vector +% A - dictionary +% param - structure containing additional parameters +% Output: +% x - sparse solution +% cost - Objective cost + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + +% Initial guess for the coefficients + +if (isfield(param, 'initcoeff')) + initX = param.initcoeff; +else + initX = zeros(size(A,2),size(b,2)); +end + +% to - 1/(step size) . It is larger than spectral norm of dictionary A + +if isfield(param, 'to') + to = param.to; +else + to = .1+svds(A,1); +end + +% lambda - Lagrangian multiplier. (regulates shrinkage) + +if isfield(param, 'lambda') + lambda = param.lambda; +else + lambda = 2*.2; +end + +% Inner-loop maximum iteration number. + +if isfield(param, 'iternum') + maxIT = param.iternum; +else + maxIT = 1000; +end + +% Stopping criterion for iterative softthresholding + +if isfield(param, 'epsilon') + epsilon = param.epsilon; +else + epsilon = 1e-7; +end + +% Debiasing. 0 = No, 1 = Yes + +if isfield(param, 'map') + map = param.map; +else + map = 1; +end + + +[X, cost] = mm1(A,b,initX,to,lambda,maxIT,epsilon,map); +cost +end \ No newline at end of file
--- a/Problems/generateAudioDeclippingProblem.m Fri Aug 12 11:17:47 2011 +0100 +++ b/Problems/generateAudioDeclippingProblem.m Mon Aug 22 11:46:35 2011 +0100 @@ -103,7 +103,7 @@ data.fs = x.fs; data.nbits = x.nbits; - +data.Upper_Limit = max(solutiondata.XClean); [data.m, data.n] = size(x_clip); data.p = windowSize*redundancyFactor; %number of dictionary elements
--- a/examples/AudioInpainting/Audio_Declipping_Example.m Fri Aug 12 11:17:47 2011 +0100 +++ b/examples/AudioInpainting/Audio_Declipping_Example.m Mon Aug 22 11:46:35 2011 +0100 @@ -49,7 +49,7 @@ % Defining the Problem structure -SMALL.Problem = generateAudioDeclippingProblem('', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2); +SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2); for idxSolver = 1:4
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,270 @@ +%% Dictionary Learning for Automatic Music Transcription - KSVD vs SPAMS +% +% +% This file contains an example of how SMALLbox can be used to test diferent +% dictionary learning techniques in Automatic Music Transcription problem. +% It calls generateAMT_Learning_Problem that will let you to choose midi, +% wave or mat file to be transcribe. If file is midi it will be first +% converted to wave and original midi file will be used for comparison with +% results of dictionary learning and reconstruction. +% The function will generarte the Problem structure that is used to learn +% Problem.p notes spectrograms from training set Problem.b using +% dictionary learning technique defined in DL structure. +% Two dictionary learning techniques were compared: +% +% - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient +% Implementation of the K-SVD Algorithm using Batch Orthogonal +% Matching Pursuit", Technical Report - CS, Technion, April 2008. +% +% - MMDL - M. Yaghoobi, T. Blumensath and M. Davies, "Dictionary Learning +% for Sparse Approximations with the Majorization Method", IEEE +% Trans. on Signal Processing, Vol. 57, No. 6, pp 2178-2191, +% 2009. + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +%% + +clear; + + +% Defining Automatic Transcription of Piano tune as Dictionary Learning +% Problem + +SMALL.Problem = generateAMT_Learning_Problem('',2048,0.75); + +%% +% Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in +% SMALL.Problem.p) using sparsity constrain only + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL(); + +% Defining fields needed for dictionary learning + +SMALL.DL(1).toolbox = 'KSVD'; +SMALL.DL(1).name = 'ksvd'; +% Defining the parameters for KSVD +% In this example we are learning 88 atoms in 100 iterations, so that +% every frame in the training set can be represented with maximum Tdata +% dictionary elements. Type help ksvd in MATLAB prompt for more options. + +SMALL.DL(1).param=struct(... + 'Tdata', 5,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 50); + +% Learn the dictionary + +SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); + +% Set SMALL.Problem.A dictionary and reconstruction function +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(1).D; +SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values +% As an example, SPAMS (Julien Mairal 2009) implementation of LARS +% algorithm is used for representation of training set in the learned +% dictionary. + +SMALL.solver(1)=SMALL_init_solver; + +% Defining the parameters needed for sparse representation + +SMALL.solver(1).toolbox='SMALL'; +SMALL.solver(1).name='SMALL_cgp'; + +% Here we use mexLasso mode=2, with lambda=2, lambda2=0 and positivity +% constrain (type 'help mexLasso' for more information about modes): +% +% min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2 + +SMALL.solver(1).param='20, 1e-2'; +% struct(... +% 'lambda', 2,... +% 'pos', 1,... +% 'mode', 2); + +% Call SMALL_soolve to represent the signal in the given dictionary. +% As a final command SMALL_solve will call above defined reconstruction +% function to reconstruct the training set (Problem.b) in the learned +% dictionary (Problem.A) + +SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +%% +% Analysis of the result of automatic music transcription. If groundtruth +% exists, we can compare transcribed notes and original and get usual +% True Positives, False Positives and False Negatives measures. + +if ~isempty(SMALL.Problem.notesOriginal) + AMT_res(1) = AMT_analysis(SMALL.Problem, SMALL.solver(1)); +end + + + +%% +% % Here we solve the same problem using non-negative sparse coding with +% % SPAMS online dictionary learning (Julien Mairal 2009) +% % +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values +% As an example, SPAMS (Julien Mairal 2009) implementation of LARS +% algorithm is used for representation of training set in the learned +% dictionary. + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for sparse representation + +SMALL.solver(2).toolbox='SPAMS'; +SMALL.solver(2).name='mexLasso'; + +% Here we use mexLasso mode=2, with lambda=3, lambda2=0 and positivity +% constrain (type 'help mexLasso' for more information about modes): +% +% min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2 + +SMALL.solver(2).param=struct('lambda', 3, 'pos', 1, 'mode', 2); + + +% You can also test ALPS, IST from MMbox or any other solver, but results +% are not as good as SPAMS +% +% % Initialising solver structure +% % Setting solver structure fields (toolbox, name, param, solution, +% % reconstructed and time) to zero values +% +% SMALL.solver(2)=SMALL_init_solver; +% +% % Defining the parameters needed for image denoising +% +% SMALL.solver(2).toolbox='ALPS'; +% SMALL.solver(2).name='AlebraicPursuit'; +% +% SMALL.solver(2).param=struct(... +% 'sparsity', 10,... +% 'memory', 1,... +% 'mode', 6,... +% 'iternum', 100,... +% 'tau',-1,... +% 'tolerance', 1e-14',... +% 'verbose',1); + +% % Initialising Dictionary structure +% % Setting Dictionary structure fields (toolbox, name, param, D and time) +% % to zero values +% % Initialising solver structure +% % Setting solver structure fields (toolbox, name, param, solution, +% % reconstructed and time) to zero values +% +% SMALL.solver(2)=SMALL_init_solver; +% +% % Defining the parameters needed for image denoising +% +% SMALL.solver(2).toolbox='MMbox'; +% SMALL.solver(2).name='mm1'; +% SMALL.solver(2).param=struct(... +% 'lambda',50,... +% 'iternum',1000,... +% 'map',0); + +SMALL.DL(2)=SMALL_init_DL('MMbox', 'MM_cn', '', 1); + + +% Defining the parameters for Majorization Minimization dictionary update +% +% In this example we are learning 88 atoms in 200 iterations, so that + + +SMALL.DL(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.A,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 200,... + 'iterDictUpdate', 1000,... + 'epsDictUpdate', 1e-7,... + 'cvset',0,... + 'show_dict', 0); + + +% Learn the dictionary + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Set SMALL.Problem.A dictionary and reconstruction function +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(2).D; +SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem); + + +% Call SMALL_soolve to represent the signal in the given dictionary. +% As a final command SMALL_solve will call above defined reconstruction +% function to reconstruct the training set (Problem.b) in the learned +% dictionary (Problem.A) + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + + +% Analysis of the result of automatic music transcription. If groundtruth +% exists, we can compare transcribed notes and original and get usual +% True Positives, False Positives and False Negatives measures. + +if ~isempty(SMALL.Problem.notesOriginal) + AMT_res(2) = AMT_analysis(SMALL.Problem, SMALL.solver(2)); +end + + +% Plot results and save midi files + +if ~isempty(SMALL.Problem.notesOriginal) + figAMT = SMALL_AMT_plot(SMALL, AMT_res); +else + figAMT = figure('Name', 'Automatic Music Transcription KSVD vs SPAMS'); + subplot(2,1,1); plot(SMALL.solver(1).reconstructed.notes(:,5), SMALL.solver(1).reconstructed.notes(:,3), 'kd '); + title (sprintf('%s dictionary in %.2f s', SMALL.DL(1).name, SMALL.DL(1).time)); + xlabel('Time'); + ylabel('Note Number'); + subplot(2,1,2); plot(SMALL.solver(2).reconstructed.notes(:,5), SMALL.solver(2).reconstructed.notes(:,3), 'b* '); + title (sprintf('%s dictionary in %.2f s', SMALL.DL(2).name, SMALL.DL(2).time)); + xlabel('Time'); + ylabel('Note Number'); +end + +FS=filesep; + +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'results']); + +[filename,pathname] = uiputfile({' *.mid;' },'Save KSVD result midi'); +if filename~=0 writemidi(SMALL.solver(1).reconstructed.midi, [pathname,FS,filename]);end + +[filename,pathname] = uiputfile({' *.mid;' },'Save SPAMS result midi'); +if filename~=0 writemidi(SMALL.solver(2).reconstructed.midi, [pathname,FS,filename]);end + +[filename,pathname] = uiputfile({' *.fig;' },'Save KSVD vs SPAMS AMT figure'); +if filename~=0 saveas(figAMT, [pathname,FS,filename]);end + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_ImgDenoise_DL_test_KSVDvsMajorizationMinimization.m Mon Aug 22 11:46:35 2011 +0100 @@ -0,0 +1,215 @@ +%% Dictionary Learning for Image Denoising - KSVD vs Recursive Least Squares +% +% This file contains an example of how SMALLbox can be used to test different +% dictionary learning techniques in Image Denoising problem. +% It calls generateImageDenoiseProblem that will let you to choose image, +% add noise and use noisy image to generate training set for dictionary +% learning. +% Two dictionary learning techniques were compared: +% +% - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient +% Implementation of the K-SVD Algorithm using Batch Orthogonal +% Matching Pursuit", Technical Report - CS, Technion, April 2008. +% +% - MMDL - M. Yaghoobi, T. Blumensath and M. Davies, "Dictionary Learning +% for Sparse Approximations with the Majorization Method", IEEE +% Trans. on Signal Processing, Vol. 57, No. 6, pp 2178-2191, 2009. + + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + + + +% If you want to load the image outside of generateImageDenoiseProblem +% function uncomment following lines. This can be useful if you want to +% denoise more then one image for example. +% Here we are loading test_image.mat that contains structure with 5 images : lena, +% barbara,boat, house and peppers. +clear; +TMPpath=pwd; +FS=filesep; +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'images']); +load('test_image.mat'); +cd(TMPpath); + +% Deffining the noise levels that we want to test + +noise_level=[10 20 25 50 100]; + +% Here we loop through different noise levels and images + +for noise_ind=2:2 +for im_num=1:1 + +% Defining Image Denoising Problem as Dictionary Learning +% Problem. As an input we set the number of training patches. + +SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); +SMALL.Problem.name=int2str(im_num); + +Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain; +maxatoms = floor(prod(SMALL.Problem.blocksize)/2); + +% results structure is to store all results + +results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr; + +%% +% Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL(); + +% Defining the parameters needed for dictionary learning + +SMALL.DL(1).toolbox = 'KSVD'; +SMALL.DL(1).name = 'ksvd'; + +% Defining the parameters for KSVD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (Edata) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(1).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'exact', 1, ... + 'iternum', 20,... + 'memusage', 'high'); + +% Learn the dictionary + +SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(1).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(1)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(1).toolbox='ompbox'; +SMALL.solver(1).name='omp2'; +SMALL.solver(1).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +% Show PSNR after reconstruction + +SMALL.solver(1).reconstructed.psnr + +%% +% For comparison purposes we will denoise image with Majorization +% Minimization method +% + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(2).toolbox='ompbox'; +SMALL.solver(2).name='omp2'; +SMALL.solver(2).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(2)=SMALL_init_DL('MMbox', 'MM_cn', '', 1); + + +% Defining the parameters for MOD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (EData) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'iterDictUpdate', 1000,... + 'epsDictUpdate', 1e-7,... + 'cvset',0,... + 'show_dict', 0); + +% Learn the dictionary + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(2).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + + + +% show results % + +SMALL_ImgDeNoiseResult(SMALL); + +results(noise_ind,im_num).psnr.ksvd=SMALL.solver(1).reconstructed.psnr; +results(noise_ind,im_num).psnr.odct=SMALL.solver(2).reconstructed.psnr; +results(noise_ind,im_num).vmrse.ksvd=SMALL.solver(1).reconstructed.vmrse; +results(noise_ind,im_num).vmrse.odct=SMALL.solver(2).reconstructed.vmrse; +results(noise_ind,im_num).ssim.ksvd=SMALL.solver(1).reconstructed.ssim; +results(noise_ind,im_num).ssim.odct=SMALL.solver(2).reconstructed.ssim; + + +results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; + +%clear SMALL; +end +end +% save results.mat results
--- a/solvers/CVX_add_const_Audio_declipping.m Fri Aug 12 11:17:47 2011 +0100 +++ b/solvers/CVX_add_const_Audio_declipping.m Mon Aug 22 11:46:35 2011 +0100 @@ -32,8 +32,8 @@ b_ineq_pos = wa_pos(:)*clippingLevelEst; b_ineq_neg = -wa_neg(:)*clippingLevelEst; if isfield(Problem,'Upper_Limit') && ~isempty(Problem.Upper_Limit) - b_ineq_pos_upper_limit = wa_pos(:)*Problem.Upper_Limit*clippingLevelEst; - b_ineq_neg_upper_limit = -wa_neg(:)*Problem.Upper_Limit*clippingLevelEst; + b_ineq_pos_upper_limit = wa_pos(:)*Problem.Upper_Limit; + b_ineq_neg_upper_limit = -wa_neg(:)*Problem.Upper_Limit; else b_ineq_pos_upper_limit = Inf; b_ineq_neg_upper_limit = -Inf;
--- a/util/SMALL_learn.m Fri Aug 12 11:17:47 2011 +0100 +++ b/util/SMALL_learn.m Mon Aug 22 11:46:35 2011 +0100 @@ -69,7 +69,20 @@ for i = 1: size(DL.D,2) DL.D(:,i)=DL.D(:,i)/norm(DL.D(:,i)); end + D = DL.D; + +elseif strcmpi(DL.toolbox,'MMbox') + + DL = wrapper_mm_DL(Problem, DL); + + % we need to make sure that columns are normalised to + % unit lenght. + + for i = 1: size(DL.D,2) + DL.D(:,i)=DL.D(:,i)/norm(DL.D(:,i)); + end D = DL.D; + % To introduce new dictionary learning technique put the files in % your Matlab path. Next, unique name <TolboxID> for your toolbox needs % to be defined and also prefferd API for toolbox functions <Preffered_API>
--- a/util/SMALL_midiGenerate.m Fri Aug 12 11:17:47 2011 +0100 +++ b/util/SMALL_midiGenerate.m Mon Aug 22 11:46:35 2011 +0100 @@ -1,5 +1,5 @@ function reconstructed=SMALL_midiGenerate(V, Problem) -%% Reconstraction of midi file from representation in the given dictionary +%% Reconstruction of midi file from representation in the given dictionary % % SMALL_midiGenerate is a part of SMALLbox and can be use to reconstruct % a midi file given representation of the training set (V) in the @@ -22,7 +22,7 @@ fs=Problem.fs; % Sampling rate f=Problem.f; % vector of frequencies at wihch spectrogram is computed -ts=(Problem.windowSize*Problem.overlap)/fs; %size of an analysis frame in seconds +ts=(Problem.windowSize*(1-Problem.overlap))/fs; %size of an analysis frame in seconds %% % Components pitch estimation using modified SWIPE algorithm by Arthuro
--- a/util/SMALL_solve.m Fri Aug 12 11:17:47 2011 +0100 +++ b/util/SMALL_solve.m Mon Aug 22 11:46:35 2011 +0100 @@ -86,7 +86,15 @@ A = opToMatrix(Problem.A, 1); end [y, numiter, time, y_path] = wrapper_ALPS_toolbox(b, A, solver.param); +elseif (strcmpi(solver.toolbox, 'MMbox')) + if ~isa(Problem.A,'float') + % ALPS does not accept implicit dictionary definition + A = opToMatrix(Problem.A, 1); + end + [y, cost] = wrapper_mm_solver(b, A, solver.param); + + % To introduce new sparse representation algorithm put the files in % your Matlab path. Next, unique name <TolboxID> for your toolbox and