Mercurial > hg > smallbox
changeset 164:4205744092e6 release_1.9
Merge from branch "ivand_dev"
| author | Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk> |
|---|---|
| date | Wed, 07 Sep 2011 14:17:30 +0100 |
| parents | af5abc34a5e1 (current diff) 855025f4c779 (diff) |
| children | 775caafd5651 |
| files | |
| diffstat | 82 files changed, 5826 insertions(+), 156 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/Majorization Minimization DL/AudioDicoLearning/AudioDicoLearning_Demo.m Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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 Wed Sep 07 14:17:30 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))^2; +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 = 0; +end + + +[X, cost] = mm1(A,b,initX,to,lambda,maxIT,epsilon,map); +cost +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/SMALL_two_step_DL.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,142 @@ +function DL=SMALL_two_step_DL(Problem, DL) + +% determine which solver is used for sparse representation % + +solver = DL.param.solver; + +% determine which type of udate to use ('KSVD', 'MOD', 'ols' or 'mailhe') % + +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 + +% flow: 'sequential' or 'parallel'. If sequential, the residual is updated +% after each atom update. If parallel, the residual is only updated once +% the whole dictionary has been computed. Sequential works better, there +% may be no need to implement parallel. Not used with MOD. + +if isfield(DL.param,'flow') + flow = DL.param.flow; +else + flow = 'sequential'; +end + +% learningRate. If the type is 'ols', it is the descent step of +% the gradient (typical choice: 0.1). If the type is 'mailhe', the +% descent step is the optimal step*rho (typical choice: 1, although 2 +% or 3 seems to work better). Not used for MOD and KSVD. + +if isfield(DL.param,'learningRate') + learningRate = DL.param.learningRate; +else + learningRate = 0.1; +end + +% number of iterations (default is 40) % + +if isfield(DL.param,'iternum') + iternum = DL.param.iternum; +else + iternum = 40; +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) + +tmpTraining = Problem.b1; +Problem.b1 = sig; +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); + [dico, solver.solution] = dico_update(dico, sig, solver.solution, ... + typeUpdate, flow, learningRate); + 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 + +Problem.b1 = tmpTraining; +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/two-step DL/dico_color.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,46 @@ +function colors = dico_color(dico, mu) + % DICO_COLOR cluster a dictionary in pairs of high correlation atoms. + % Called by dico_decorr. + % + % Parameters: + % -dico: the dictionary + % -mu: the correlation threshold + % + % Result: + % -colors: a vector of indices. Two atoms with the same color have a + % correlation greater than mu + + numAtoms = length(dico); + colors = zeros(numAtoms, 1); + + % compute the correlations + G = abs(dico'*dico); + G = G-eye(size(G)); + + % iterate on the correlations higher than mu + c = 1; + maxCorr = max(max(G)); + while maxCorr > mu + % find the highest correlated pair + x = find(max(G)==maxCorr, 1); + y = find(G(x,:)==maxCorr, 1); + + % color them + colors(x) = c; + colors(y) = c; + c = c+1; + + % make sure these atoms never get selected again + G(x,:) = 0; + G(:,x) = 0; + G(y,:) = 0; + G(:,y) = 0; + + % find the next correlation + maxCorr = max(max(G)); + end + + % complete the coloring with singletons + index = find(colors==0); + colors(index) = c:c+length(index)-1; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/dico_decorr.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,48 @@ +function dico = dico_decorr(dico, mu, amp) + %DICO_DECORR decorrelate a dictionary + % Parameters: + % dico: the dictionary + % mu: the coherence threshold + % amp: the amplitude coefficients, only used to decide which atom to + % project + % + % Result: + % dico: a dictionary close to the input one with coherence mu. + + % compute atom weights + if nargin > 2 + rank = sum(amp.*amp, 2); + else + rank = randperm(length(dico)); + end + + % 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 + 10^-6 + % find pairs of high correlation atoms + colors = dico_color(dico, mu); + + % iterate on all pairs + nbColors = max(colors); + for c = 1:nbColors + index = find(colors==c); + if numel(index) == 2 + % decide which atom to change (the one with lowest weight) + if rank(index(1)) < rank(index(2)) + index = fliplr(index); + end + + % update the atom + corr = dico(:,index(1))'*dico(:,index(2)); + alpha = sqrt((1-mu*mu)/(1-corr*corr)); + beta = corr*alpha-mu*sign(corr); + dico(:,index(2)) = alpha*dico(:,index(2))... + -beta*dico(:,index(1)); + end + end + niter = niter+1; + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/two-step DL/dico_update.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,107 @@ +function [dico, amp] = dico_update(dico, sig, amp, type, flow, rho) + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % [dico, amp] = dico_update(dico, sig, amp, type, flow, rho) + % + % perform one iteration of dictionary update for dictionary learning + % + % parameters: + % - dico: the initial dictionary with atoms as columns + % - sig: the training data + % - amp: the amplitude coefficients as a sparse matrix + % - type: the algorithm can be one of the following + % - ols: fixed step gradient descent + % - mailhe: optimal step gradient descent (can be implemented as a + % default for ols?) + % - MOD: pseudo-inverse of the coefficients + % - KSVD: already implemented by Elad + % - flow: 'sequential' or 'parallel'. If sequential, the residual is + % updated after each atom update. If parallel, the residual is only + % updated once the whole dictionary has been computed. Sequential works + % better, there may be no need to implement parallel. Not used with + % MOD. + % - rho: learning rate. If the type is 'ols', it is the descent step of + % the gradient (typical choice: 0.1). If the type is 'mailhe', the + % descent step is the optimal step*rho (typical choice: 1, although 2 + % or 3 seems to work better). Not used for MOD and KSVD. + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if ~ exist( 'rho', 'var' ) || isempty(rho) + rho = 0.1; + end + + if ~ exist( 'flow', 'var' ) || isempty(flow) + flow = sequential; + end + + res = sig - dico*amp; + nb_pattern = size(dico, 2); + + switch type + case 'rand' + x = rand(); + if x < 1/3 + type = 'MOD'; + elseif type < 2/3 + type = 'mailhe'; + else + type = 'KSVD'; + end + end + + switch type + case 'MOD' + G = amp*amp'; + dico2 = sig*amp'*G^-1; + for p = 1:nb_pattern + n = norm(dico2(:,p)); + % renormalize + if n > 0 + dico(:,p) = dico2(:,p)/n; + amp(p,:) = amp(p,:)*n; + end + end + case 'ols' + for p = 1:nb_pattern + grad = res*amp(p,:)'; + if norm(grad) > 0 + pat = dico(:,p) + rho*grad; + pat = pat/norm(pat); + if nargin >5 && strcmp(flow, 'sequential') + res = res + (dico(:,p)-pat)*amp(p,:); %#ok<*NASGU> + end + dico(:,p) = pat; + end + end + case 'mailhe' + for p = 1:nb_pattern + grad = res*amp(p,:)'; + if norm(grad) > 0 + pat = (amp(p,:)*amp(p,:)')*dico(:,p) + rho*grad; + pat = pat/norm(pat); + if nargin >5 && strcmp(flow, 'sequential') + res = res + (dico(:,p)-pat)*amp(p,:); + end + dico(:,p) = pat; + end + end + case 'KSVD' + for p = 1:nb_pattern + index = find(amp(p,:)~=0); + if ~isempty(index) + patch = res(:,index)+dico(:,p)*amp(p,index); + [U,S,V] = svd(patch); + if U(:,1)'*dico(:,p) > 0 + dico(:,p) = U(:,1); + else + dico(:,p) = -U(:,1); + end + dico(:,p) = dico(:,p)/norm(dico(:,p)); + amp(p,index) = dico(:,p)'*patch; + if nargin >5 && strcmp(flow, 'sequential') + res(:,index) = patch-dico(:,p)*amp(p,index); + end + end + end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/AMT_reconstruct.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,117 @@ +function reconstructed=AMT_reconstruct(V, Problem) +%% 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 +% dictionary Problem.A. +% Output is reconstructed structure with two fields: +% - reconstructed.notes - matrix with transcribed notes +% - reconstructed.midi - midi representation of transcription + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2009 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. +%% +U=Problem.A; % Dictionary used for representation +fs=Problem.fs; % Sampling rate +f=Problem.f; % vector of frequencies at wihch spectrogram is computed + +ts=(Problem.windowSize*(1-Problem.overlap))/fs; %size of an analysis frame in seconds + +%% +% Components pitch estimation using modified SWIPE algorithm by Arthuro +% Camacho +% +% Columns of matrix U are spectrograms of the notes learned from the +% training set. We are estimating pitches of these notes by also +% restricting pitch values to the one of the 88 piano notes. + +pitch=zeros(size(U,2),1); + +for i=1:size(U,2) + + pitch(i) = SMALL_swipe(U(:,i),fs, f, [27.50 8192], 1/12); + +end + +%% +% If some of columns of U have the same pitch, their contribution to the +% score (matrix V) is summed. + +[Ps,idx]=sort(pitch); +ndp=1; +Pd(ndp)=Ps(1); +Vnew(ndp,:)=V(idx(1),:); +for i=2:88 + if Ps(i)> Ps(i-1) + + ndp=ndp+1; + Vnew(ndp,:)=V(idx(i),:); + Pd(ndp)=Ps(i); + + else + Vnew(ndp,:)=Vnew(ndp,:)+V(idx(i),:); + end +end +%% +% Generate midi matrix + +midx=0; +for i=1:ndp + + % Threshold for finding onsets and offsets of notes + + thr=mean(Vnew(i,:));%+std(Vnew(i,:)); + + if(Pd(i)~=0) + for j=1:size(Vnew,2) + if Vnew(i,j)<thr + Vnew(i,j)=0; + if(j>1) + if (Vnew(i,j-1)==1) + try + M(midx,6)=(j-1)*ts; + if (M(midx,6)-M(midx,5))<2*ts + midx=midx-1; + end + catch + pause; + end + end + end + else + Vnew(i,j)=1; + if(j>1) + if (Vnew(i,j-1)==0) + midx=midx+1; + M(midx,1)=1; + M(midx,2)=1; + M(midx,3)=69 +round( 12 *log2(Pd(i)/440)); + M(midx,4)=80; + M(midx,5)=(j-1)*ts; + end + else + midx=midx+1; + M(midx,1)=1; + M(midx,2)=1; + M(midx,3)=69 + round(12 *log2(Pd(i)/440)); + M(midx,4)=80; + M(midx,5)=0; + end + end + end + if M(midx,6)==0 + M(midx,6)=(j-1)*ts; + end + end +end + +M=sortrows(M,5); +reconstructed.notes=M; +reconstructed.midi = matrix2midi(M);
--- a/Problems/AudioDeclipping_reconstruct.m Tue Jul 26 16:02:59 2011 +0100 +++ b/Problems/AudioDeclipping_reconstruct.m Wed Sep 07 14:17:30 2011 +0100 @@ -1,8 +1,15 @@ -function reconstructed=AudioDeclipping_reconstruct(y, Problem, SparseDict) +function reconstructed = AudioDeclipping_reconstruct(y, Problem) %% Audio declipping Problem reconstruction function % % This reconstruction function is using sparse representation y % in dictionary Problem.A to reconstruct declipped audio. +% The output structure has following fields: +% audioAllSamples - signal with all samples taken from reconstructed +% signal +% audioOnlyClipped - only clipped samples are reconstructed, +% others are taken from original signal +% snrAll - psnr of whole signal +% snrMiss - psnr of the reconstructed clipped samples % % [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an % evaluation framework for sparse representations and dictionary
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/AudioDenoise_reconstruct.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,53 @@ +function reconstructed=AudioDenoise_reconstruct(y, Problem) +%% Audio denoising Problem reconstruction function +% +% This reconstruction function is using sparse representation y +% in dictionary Problem.A to reconstruct denoised audio. +% The output structre has following fields: +% audio - denoised audio signal +% psnr - psnr of the reconstructed audio signal +% +% [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an +% evaluation framework for sparse representations and dictionary +% learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, +% Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425 + +% +% 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. +%% + +windowSize = Problem.windowSize; +overlap = Problem.overlap; +ws = Problem.ws(windowSize); +wa = Problem.wa(windowSize); + +A = Problem.A; + +orig = Problem.Original; +noisy = Problem.Noisy; + + +% reconstruct audio frames + +xFrames = diag(ws)*(A*y); +wNormFrames = (ws.*wa)'*ones(1,size(xFrames,2)); + +% overlap and add + +rec = col2imstep(xFrames, size(noisy), [windowSize 1], [windowSize*overlap 1]); +wNorm = col2imstep(wNormFrames, size(noisy), [windowSize 1], [windowSize*overlap 1]); +wNorm(find(wNorm==0)) = 1; +recN = rec./wNorm; + +%% output structure image+psnr %% +reconstructed.audio = recN; +reconstructed.psnr = 20*log10(sqrt(numel(orig)) / norm(orig - reconstructed.audio)); + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/ImageDenoise_reconstruct.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,66 @@ +function reconstructed=ImageDenoise_reconstruct(y, Problem, SparseDict) +%% Image Denoising Problem reconstruction function +% +% This reconstruction function is using sparse representation y +% in dictionary Problem.A to reconstruct the patches of the denoised +% image. + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2009 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. +%% + + +% stepsize % +if (isfield(Problem,'stepsize')) + stepsize = Problem.stepsize; + if (numel(stepsize)==1) + stepsize = ones(1,2)*stepsize; + end +else + stepsize = ones(1,2); +end +if (any(stepsize<1)) + error('Invalid step size.'); +end + +% lambda % +if (isfield(Problem,'lambda')) + lambda = Problem.lambda; +else + lambda = Problem.maxval/(10*Problem.sigma); +end +if exist('SparseDict','var')&&(SparseDict==1) + if issparse(Problem.A) + A = Problem.A; + else + A = sparse(Problem.A); + end + cl_samp=add_dc(dictsep(Problem.basedict,A,y), Problem.b1dc,'columns'); +else + cl_samp=add_dc(Problem.A*y, Problem.b1dc,'columns'); +end +% combine the patches into reconstructed image +cl_im=col2imstep(cl_samp, size(Problem.Noisy), Problem.blocksize); + +cnt = countcover(size(Problem.Noisy),Problem.blocksize,stepsize); + +im = (cl_im+lambda*Problem.Noisy)./(cnt + lambda); +% y(y~=0)=1; +% numD=sum(y,2); +% nnzy=sum(y,1); +% figure(200);plot(sort(numD)); +% figure(201);plot(sort(nnzy)); +[v.RMSErn, v.RMSEcd, v.rn_im, v.cd_im]=SMALL_vmrse_type2(Problem.Original, Problem.Noisy, im); +%% output structure image+psnr %% +reconstructed.Image=im; +reconstructed.psnr = 20*log10(Problem.maxval * sqrt(numel(Problem.Original(:))) / norm(Problem.Original(:)-im(:))); +reconstructed.vmrse=v; +reconstructed.ssim=SMALL_ssim_index(Problem.Original, im); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/generateAMTProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,91 @@ +function data = generateAMTProblem(nfft, windowSize, overlap) +%% Generate Automatic Music Transcription Problem +% +% generateAMT_Learning_Problem is a part of the SMALLbox and generates +% a problem that can be used for comparison of Dictionary Learning/Sparse +% Representation techniques in automatic music transcription scenario. +% The function prompts a user for an audio file (mid, wav, mat) reads it +% and generates a spectrogram given fft size (default nfft=4096), analysis +% window size (windowSize=2822), and analysis window overlap (overlap = +% 0.5). +% +% The output of the function is stucture with following fields: +% b - matrix with magnitudes of the spectrogram +% f - vector of frequencies at wihch spectrogram is computed +% windowSize - analysis window size +% overlap - analysis window overlap +% fs - sampling frequency +% m - number of frequenciy points in spectrogram +% n - number of time points in the spectrogram +% p - number of dictionary elements to be learned (eg 88 for piano) +% notesOriginal - notes of the original audio to be used for +% comparison (if midi of the original exists) +% name - name of the audio file to transcribe + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2009 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. +% +%% +FS=filesep; +if ~ exist( 'nfft', 'var' ) || isempty(nfft), nfft = 4096; end +if ~ exist( 'windowSize', 'var' ) || isempty(windowSize), windowSize = 2822; end +if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 0.5; end + +%% +%ask for file name +TMPpath=pwd; +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'audio']); +[filename,pathname] = uigetfile({'*.mat; *.mid; *.wav'},'Select a file to transcribe'); +[pathstr, name, ext, versn] = fileparts(filename); +data.name=name; + +data.notesOriginal=[]; + +if strcmp(ext,'.mid') + midi=readmidi(filename); + data.notesOriginal=midiInfo(midi); + y=midi2audio(midi); + wavwrite(y, 44100, 16, 'temp.wav'); + [x.signal, x.fs, x.nbits]=wavread('temp.wav'); + delete('temp.wav'); +elseif strcmp(ext,'.wav') + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); + filename1=[name, '.mid']; + if exist(filename1, 'file') + midi=readmidi(filename1); + data.notesOriginal=midiInfo(midi); + end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'wav']); + [x.signal, x.fs, x.nbits]=wavread(filename); +else + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); + filename1=[name, '.mid']; + if exist(filename1, 'file') + midi=readmidi(filename1); + data.notesOriginal=midiInfo(midi); + end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'mat']); + x=load([pathname,filename]); +end +%% +[X, frX]=spectrogram(x.signal, hanning(windowSize), overlap*windowSize, nfft, x.fs); +%% +data.b=abs(X); +data.f=frX; +data.windowSize=windowSize; +data.overlap=overlap; +data.fs=x.fs; +data.m=size(X,1); +data.n=size(X,2); + +data.p=88; %number of dictionary elements (ie notes to recover) +cd(TMPpath); + +end
--- a/Problems/generateAudioDeclippingProblem.m Tue Jul 26 16:02:59 2011 +0100 +++ b/Problems/generateAudioDeclippingProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -5,6 +5,35 @@ % Audio declipping is a problem proposed in Audio Inpaining Toolbox and % in [2]. % +% The function takes as an optional input +% soundfile - name of the file +% clippingLevel - (default 0.6) +% windowSize - 1D frame size (eg 512) +% overlap - ammount of overlaping frames between 0 and 1 +% wa,ws,wd - analisys, synthesis and dictionary window functions +% +% Dict_fun - function to be used to generate dictionary +% redundancyFactor - overcompletness of dictionary (default 2) +% +% The function outputs the structure with following fields: +% original - original signal +% clipped - clipped signal +% clipMask - mask indicating clipped samples +% clippingLevel - (default 0.6) +% Upper_Limit - maximum value of original data +% fs - sample rate of the original signal in Hertz +% nbits - the number of bits per sample +% sigma - added noise level +% B - dictionary to be used for sparse representation +% M - measurement matrix (non-clipped data in b) +% b - matrix of clipped frames +% m - size od dictionary atom +% n - number of frames to be represented +% p - number of atoms in dictionary +% windowSize - 1D frame size (eg 512) +% overlap - ammount of overlaping frames between 0 and 1 +% wa,ws, wd - analisys, synthesis and dictionary window functions +% % [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an % evaluation framework for sparse representations and dictionary % learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, @@ -103,7 +132,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/Problems/generateAudioDenoiseProblem.m Tue Jul 26 16:02:59 2011 +0100 +++ b/Problems/generateAudioDenoiseProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -1,20 +1,40 @@ -function data=generateAudioDenoiseProblem(au, trainnum, blocksize, dictsize, overlap, sigma, gain, maxval, initdict); -%% Audio Denoising Problem - needs revision, not yet finalised +function data = generateAudioDenoiseProblem(soundfile, sigma, windowSize,... + overlap, wa, ws, trainnum, redundancyFactor, initdict) +%% Audio Denoising Problem % % generateAudioDenoiseProblem is part of the SMALLbox and generate a % problem for comaprison of Dictionary Learning/Sparse Representation -% techniques in audio denoising scenario. It is based on KSVD image -% denoise demo by Ron Rubinstein (see bellow). -% The fuction takes as an optional input -% au - audio samples to be denoised -% trainnum - number of frames for training -% blocksize - 1D frame size (eg 512) -% dictsize - number of atoms to be trained -% overlap - ammount of overlaping frames between 0 and 1 +% techniques in audio denoising scenario. +% +% The function takes as an optional input +% soundfile - name of the file +% sigma - noise level (dB) +% windowSize - 1D frame size (eg 512) +% overlap - ammount of overlaping frames between 0 and 1 +% wa,ws - analisys and synthesis window functions +% +% trainnum - number of frames for training +% redundancyFactor - overcompletness of dictionary (default 2) +% initdict - initial dictionary % +% The function outputs the structure with following fields: +% Original - original signal +% Noisy - signal with added noise +% fs - sample rate of the original signal in Hertz +% nbits - the number of bits per sample +% sigma - added noise level +% b - matrix of training samples for dictionary learning +% b1 - matrix containing all frames for reconstruction step +% m - size od dictionary atom +% n - number of frames for training +% p - number of atoms in dictionary +% windowSize - 1D frame size (eg 512) +% overlap - ammount of overlaping frames between 0 and 1 +% wa,ws - analisys and synthesis window functions +% initdict - initial dictionary % Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2010 Ivan Damnjanovic. +% 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 @@ -30,67 +50,69 @@ disp(' '); FS=filesep; -if ~ exist( 'sigma', 'var' ) || isempty(sigma), sigma = 26.74; end -if ~ exist( 'gain', 'var' ) || isempty(gain), gain = 1.15; end -if ~ exist( 'initdict', 'var' ) || isempty(initdict), initdict = 'odct'; end -if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 15/16; end %% prompt user for wav file %% %ask for file name TMPpath=pwd; -if ~ exist( 'au', 'var' ) || isempty(au) +if ~ exist( 'soundfile', 'var' ) || isempty(soundfile) + %ask for file name [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); - cd([pathstr1,FS,'data',FS,'audio',FS,'wav']); - [filename,pathname] = uigetfile({'*.wav;'},'Select a wav file'); + cd([pathstr1,FS,'data',FS,'audio']); + [filename,pathname] = uigetfile({'*.mat; *.mid; *.wav'},'Select a file to transcribe'); [pathstr, name, ext, versn] = fileparts(filename); data.name=name; - - au = wavread(filename); - au = mean(au,2); % turn it into mono. -end; -if ~ exist( 'maxval', 'var' ) || isempty(maxval), maxval = max(au); end -%% generate noisy audio %% - -disp(' '); -disp('Generating noisy audio...'); -sigma = max(au)/10^(sigma/20); -n = randn(size(au)) .* sigma; -aunoise = au + n;% here we can load noise audio if available - % for example: wavread('icassp06_x.wav');% - - + if strcmp(ext,'.mid') + midi=readmidi(filename); +% data.notesOriginal=midiInfo(midi); + y=midi2audio(midi); + wavwrite(y, 44100, 16, 'temp.wav'); + [x.signal, x.fs, x.nbits]=wavread('temp.wav'); + delete('temp.wav'); + elseif strcmp(ext,'.wav') +% cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); +% filename1=[name, '.mid']; +% if exist(filename1, 'file') +% midi=readmidi(filename1); +% data.notesOriginal=midiInfo(midi); +% end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'wav']); + [x.signal, x.fs, x.nbits]=wavread(filename); + else +% cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); +% filename1=[name, '.mid']; +% if exist(filename1, 'file') +% midi=readmidi(filename1); +% data.notesOriginal=midiInfo(midi); +% end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'mat']); + x=load([pathname,filename]); + end +else + [x.signal, x.fs, x.nbits]=wavread(soundfile); + [pathstr, name, ext, versn] = fileparts(soundfile); + data.name=name; +end %% set parameters %% +if ~ exist( 'sigma', 'var' ) || isempty(sigma), sigma = 0.2; end -x = aunoise; -if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 512;end -if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), dictsize = 2048;end +if ~ exist( 'windowSize', 'var' ) || isempty(windowSize), windowSize = 256;end +if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 0.5; end +if ~ exist( 'wa', 'var' ) || isempty(wa), wa = @wSine; end % Analysis window +if ~ exist( 'ws', 'var' ) || isempty(ws), ws = @wSine; end % Synthesis window -if ~ exist( 'trainnum', 'var' ) || isempty(trainnum),trainnum = (size(x,1)-blocksize+1);end +if ~ exist( 'redundancyFactor', 'var' ) || isempty(windowSize),... + redundancyFactor = 2;end +if ~ exist( 'initdict', 'var' ) || isempty(initdict),... + initdict = 'odct'; end +if ~ exist( 'trainnum', 'var' ) || isempty(trainnum), ... + trainnum = 16*redundancyFactor*windowSize;end - - - -p=1; - - -% -% msgdelta = 5; -% -% verbose = 't'; -% if (msgdelta <= 0) -% verbose=''; -% msgdelta = -1; -% end -% -% -% % initial dictionary % -% if (strcmpi(initdict,'odct')) - initdict = odctndict(blocksize,dictsize,p); + initdict = odctndict(windowSize, redundancyFactor*windowSize, 1); elseif (strcmpi(initdict,'data')) clear initdict; % causes initialization using random examples else @@ -98,45 +120,31 @@ end if exist( 'initdict', 'var' ) - initdict = initdict(:,1:dictsize); + initdict = initdict(:,1:redundancyFactor*windowSize); end -% noise mode % -% if (isfield(params,'noisemode')) -% switch lower(params.noisemode) -% case 'psnr' -% sigma = maxval / 10^(params.psnr/20); -% case 'sigma' -% sigma = params.sigma; -% otherwise -% error('Invalid noise mode specified'); -% end -% elseif (isfield(params,'sigma')) -% sigma = params.sigma; -% elseif (isfield(params,'psnr')) -% sigma = maxval / 10^(params.psnr/20); -% else -% error('Noise strength not specified'); -% end - -% params.Edata = sqrt(prod(blocksize)) * sigma * gain; % target error for omp -% params.codemode = 'error'; -% -% params.sigma = sigma; -% params.noisemode = 'sigma'; -% -% -% % make sure test data is not present in params -% if (isfield(params,'testdata')) -% params = rmfield(params,'testdata'); -% end - - %%%% create training data %%% +%% generate noisy audio %% -X = buffer( x(1:trainnum),blocksize, overlap*blocksize); +disp(' '); +disp('Generating noisy audio...'); +x.signal = x.signal/max(abs(x.signal(:)))*0.99999; +n = randn(size(x.signal)) .* sigma; + +xnoise = x.signal + n;% here we can load noise audio if available + % for example: wavread('icassp06_x.wav');% + + + + +X = im2colstep(xnoise,[windowSize 1],[overlap*windowSize 1]); +X = diag(wa(windowSize)) * X; + + + + % remove dc in blocks to conserve memory % % bsize = 2000; @@ -144,17 +152,32 @@ % blockids = i : min(i+bsize-1,size(X,2)); % X(:,blockids) = remove_dc(X(:,blockids),'columns'); % end -data.Original = au; -data.Noisy = aunoise; -data.b = X; -data.m = size(X,1); -data.n = size(X,2); -data.p = dictsize; -data.blocksize=blocksize; +data.Original = x.signal; +data.Noisy = xnoise; +data.fs = x.fs; +data.nbits = x.nbits; + data.sigma = sigma; -data.gain = gain; -data.maxval = maxval; + + +if (trainnum<size(X,2)) + p = randperm(size(X,2)); + p=sort(p(1:trainnum)); + data.b = X(:,p); +else + data.b = X; +end + +data.b1 = X; +[data.m, data.n] = size(data.b); +data.p = redundancyFactor*windowSize; + +data.windowSize = windowSize; +data.overlap = overlap; +data.ws = ws; +data.wa = wa; + data.initdict= initdict; -data.signalDim=1; + cd(TMPpath);
--- a/Problems/generateImageDenoiseProblem.m Tue Jul 26 16:02:59 2011 +0100 +++ b/Problems/generateImageDenoiseProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -1,4 +1,5 @@ -function data=generateImageDenoiseProblem(im, trainnum, blocksize, dictsize, sigma, gain, maxval, initdict); +function data = generateImageDenoiseProblem(im, trainnum, blocksize,... + dictsize, sigma, gain, maxval, initdict) %% Generate Image Denoising Problem % % generateImageDenoiseProblem is a part of the SMALLbox and generates
--- a/Problems/generateMyDummyProblem.m Tue Jul 26 16:02:59 2011 +0100 +++ b/Problems/generateMyDummyProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -13,7 +13,7 @@ %% Change copyright notice as appropriate: % Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2009 Ivan Damnjanovic. +% 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/generatePierreProblem.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,121 @@ +function data=generatePierreProblem(src, trg, blocksize, dictsize); +%% Generate Pierre Villars Problem +% +% Pierre_Problem is a part of the SMALLbox and generates the problem +% suggested by Professor Pierre Vandergheynst on the SMALL meeting in +% Villars. +% The function takes as an input: +% - src - source image matrix (if not present function promts user for +% an image file) , +% - trg - target image matrix (if not present function promts user for +% an image file) , +% - blocksize - block (patch) vertical/horizontal dimension (default 8), +% - dictsize - dictionary size (default - all patches from target +% image). +% +% The output of the function is stucture with following fields: +% - srcname - source image name, +% - imageSrc - source image matrix, +% - trgname - target image name, +% - imageTrg - Target image matrix, +% - A - dictonary with patches from the source image, +% - b - measurement matrix (i.e. patches from target image to be +% represented in dictionary A, +% - m - size of patches (default 25), +% - n - number of patches to be represented, +% - p - dictionary size, +% - blocksize - block size (default [5 5]), +% - maxval - maximum value (default - 255) +% - sparse - if 1 SMALL_solve will keep solution matrix in sparse form, +% due to memory constrains. + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2010 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. +%% prompt user for images %% + +% ask for source file name + +TMPpath=pwd; +FS=filesep; +if ~ exist( 'src', 'var' ) || isempty(src) +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'images']); +[filename,pathname] = uigetfile({'*.png;'},'Select a source image'); +[pathstr, name, ext, versn] = fileparts(filename); +data.srcname=name; +src = imread(filename); +src = double(src); +end; + +% ask for target file name + +if ~ exist( 'trg', 'var' ) || isempty(trg) +[filename,pathname] = uigetfile({'*.png;'},'Select a target image'); +[pathstr, name, ext, versn] = fileparts(filename); +data.trgname=name; +trg = imread(filename); +trg = double(trg); +end; +cd(TMPpath); + +%% set parameters %% + +maxval = 255; +if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 5;end + +if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), + dictsize = (size(src,1)-blocksize+1)*(size(src,2)-blocksize+1); + patch_idx=1:dictsize; +else + num_blocks_src=(size(src,1)-blocksize+1)*(size(src,2)-blocksize+1); + patch_idx=1:floor(num_blocks_src/dictsize):dictsize*floor(num_blocks_src/dictsize); +end + +p = ndims(src); +if (p==2 && any(size(src)==1) && length(blocksize)==1) + p = 1; +end + + +% blocksize % +if (numel(blocksize)==1) + blocksize = ones(1,p)*blocksize; +end +%% +%% create dictionary data %% + +S=im2colstep(src,blocksize); + +for j= 1:size(S,2) + S(:,j)=S(:,j)./norm(S(:,j)); +end + +%% create measurement matrix %% + +T=im2colstep(trg,blocksize, blocksize); + +%% output structure %% + +data.imageSrc = src; +data.imageTrg = trg; +data.A = S(:,patch_idx); +data.b = T; +data.m = size(T,1); +data.n = size(T,2); +data.p = size(data.A,2); +data.blocksize=blocksize; +data.maxval=maxval; + +% keep coefficients matrix in sparse form and do not convert it to full. +% getting around out of memory problem when converting big matrix from +% sparse to full... (check SMALL_solve function) +data.sparse=1; + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/ALPS solvers tests/SMALL_ImgDenoise_DL_test_KSVDvsTwoStepALPSandMahile.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,288 @@ +%% 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. +% - RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares +% Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on, +% vol.58, no.4, pp.2121-2130, April 2010 +% + + +% 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=4:4 +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) ImageDenoise_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 overcomplete DCT +% here +% Set SMALL.Problem.A dictionary to be oDCT (i.e. Problem.initdict - +% since initial dictionaruy is already set to be oDCT when generating the +% denoising problem + + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(2)=SMALL_init_solver('ALPS','AgebraicPursuit','',1); + +% Defining the parameters needed for image denoising + +SMALL.solver(2).param=struct(... + 'tolerance',1e-05,... + 'sparsity', 32,... + 'mode', 0,... + 'memory', 1,... + 'iternum', 50); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(2)=SMALL_init_DL('TwoStepDL', 'Mailhe', '', 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', 40,... + 'show_dict', 1); + +% 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) ImageDenoise_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)); + +%% +% In the b1 field all patches from the image are stored. For RLS-DLA we +% will first exclude all the patches that have l2 norm smaller then +% threshold and then take min(40000, number_of_remaining_patches) in +% ascending order as our training set (SMALL.Problem.b) + +X=SMALL.Problem.b1; +X_norm=sqrt(sum(X.^2, 1)); +[X_norm_sort, p]=sort(X_norm); +p1=p(X_norm_sort>Edata); +if size(p1,2)>40000 + p2 = randperm(size(p1,2)); + p2=sort(p2(1:40000)); + size(p2,2) + SMALL.Problem.b=X(:,p1(p2)); +else + size(p1,2) + SMALL.Problem.b=X(:,p1); + +end + +% Forgetting factor for RLS-DLA algorithm, in this case we are using +% fixed value + +lambda=0.9998 + +% Use Recursive Least Squares +% to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(3)=SMALL_init_DL(); + +% Defining fields needed for dictionary learning + +SMALL.DL(3).toolbox = 'SMALL'; +SMALL.DL(3).name = 'SMALL_rlsdla'; +SMALL.DL(3).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'forgettingMode', 'FIX',... + 'forgettingFactor', lambda,... + 'show_dict', 1000); + + +SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3)); + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.Problem.A = SMALL.DL(3).D; +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); + +SMALL.solver(3)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(3).toolbox='ompbox'; +SMALL.solver(3).name='omp2'; +SMALL.solver(3).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + + +SMALL.solver(3)=SMALL_solve(SMALL.Problem, SMALL.solver(3)); + +SMALL.solver(3).reconstructed.psnr + + +% 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).psnr.rlsdla=SMALL.solver(3).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).vmrse.rlsdla=SMALL.solver(3).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).ssim.rlsdla=SMALL.solver(3).reconstructed.ssim; + +results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; +results(noise_ind,im_num).time.rlsdla.time=SMALL.solver(3).time+SMALL.DL(3).time; +clear SMALL; +end +end +% save results.mat results
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/ALPS solvers tests/SMALL_solver_test_ALPS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,163 @@ +%% Example test of solvers from different toolboxes on Sparco problem 6 +% +% The main purpose of this example is to show how to use SMALL structure +% to solve SPARCO compressed sensing problems (1-11) and compare results +% from different solvers. +% To generate SMALL.Problem part of structure you can use generateProblem +% function from Sparco toolbox giving the problem number and any +% additional parameters you might want to change. Alternatively, you can +% might want to consult sparco documentation to write a problem by +% yourself. There are four fields the must be specified in SMALL.Problem +% - A, b, sizeA and reconstruct. +% +% To generate SMALL.solver part of the structure you must specify three +% fields: +% +% SMALL.solver.toolbox - string with toolbox name is needed because +% different toolboxes are calling solver +% functions in different ways. +% SMALL.solver.name - its string representing solver name (e.g. +% SolveOMP) +% SMALL.solver.param - string that contains optional parameters for +% particular solver (all parameters you want to +% specify except A, b and size of solution) +% +% Every call to SMALL_solve function will generate following output: +% +% SMALL.solver.solution - contains solution vector x +% SMALL.solver.reconstructed - vector containing signal reconstructed +% from the solution +% SMALL.solver.time - time that solver spent to find the solution +% +% SMALL_plot function plots the SMALL.solver.solution and reconstructed +% against original signal. +% +% In this particular example we are testing SMALL_cgp, SMALL_chol, +% SolveOMP form SparseLab and greed_pcgp form Sparsify against "PROB006 +% Daubechies basis, Gaussian ensemble measurement basis, piecewise cubic +% polynomial signal" from Sparco. +% +% + + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2009 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. +%% + +fprintf('\n\nExample test of SMALL solvers against their counterparts on Sparco problems.\n\n'); + +%% +% Generate SPARCO problem +clear + +SMALL.Problem = generateProblem(6, 'P', 6, 'm', 2*500,'n',2*1024, 'show'); + +SMALL.Problem.A = opToMatrix(SMALL.Problem.A, 1); + +%% +i=1; +%% +% ALPS test +SMALL.solver(i) = SMALL_init_solver; +SMALL.solver(i).toolbox = 'ALPS'; +SMALL.solver(i).name = 'AlgebraicPursuit'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver(i).param=struct(... + 'sparsity', 125,... + 'memory', 0,... + 'mode', 1,... + 'iternum', 50,... + 'tolerance', 1e-14'); + +SMALL.solver(i)=SMALL_solve(SMALL.Problem,SMALL.solver(i)); + + +i=i+1; +%% +% SMALL Conjugate Gradient test +SMALL.solver(i)=SMALL_init_solver; +SMALL.solver(i).toolbox='SMALL'; +SMALL.solver(i).name='SMALL_pcgp'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver(i).param='200, 1e-14'; + +SMALL.solver(i)=SMALL_solve(SMALL.Problem,SMALL.solver(i)); + + +i=i+1; + +%% +% SolveOMP from SparseLab test + +SMALL.solver(i)=SMALL_init_solver; +SMALL.solver(i).toolbox='SparseLab'; +SMALL.solver(i).name='SolveOMP'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver(i).param='200, 0, 0, 0, 1e-14'; + +SMALL.solver(i)=SMALL_solve(SMALL.Problem, SMALL.solver(i)); + +i=i+1; + +%% +% SMALL OMP with Cholesky update test +SMALL.solver(i)=SMALL_init_solver; +SMALL.solver(i).toolbox='SMALL'; +SMALL.solver(i).name='SMALL_chol'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver(i).param='200, 1e-14'; + +SMALL.solver(i)=SMALL_solve(SMALL.Problem, SMALL.solver(i)); + +i=i+1; + +%% +% greed_pcgp from Sparsify test + +SMALL.solver(i)=SMALL_init_solver; +SMALL.solver(i).toolbox='Sparsify'; +SMALL.solver(i).name='greed_pcgp'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver(i).param='''stopCrit'', ''M'', ''stopTol'', 200'; + +SMALL.solver(i)=SMALL_solve(SMALL.Problem, SMALL.solver(i)); + +%% + +SMALL_plot(SMALL); + + + + +
--- a/examples/AudioInpainting/Audio_Declipping_Example.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/AudioInpainting/Audio_Declipping_Example.m Wed Sep 07 14:17:30 2011 +0100 @@ -49,7 +49,7 @@ % Defining the Problem structure -SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 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
--- a/examples/Automatic Music Transcription/SMALL_AMT_DL_test.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Automatic Music Transcription/SMALL_AMT_DL_test.m Wed Sep 07 14:17:30 2011 +0100 @@ -31,7 +31,7 @@ % Defining Automatic Transcription of Piano tune as Dictionary Learning % Problem -SMALL.Problem = generateAMT_Learning_Problem(); +SMALL.Problem = generateAMTProblem(); %% % Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in @@ -67,7 +67,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(1).D; -SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) AMT_reconstruct(x, SMALL.Problem); %% % Initialising solver structure @@ -151,7 +151,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(2).D; -SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem); +SMALL.Problem.reconstruct=@(x) AMT_reconstruct(x, SMALL.Problem); %% % Initialising solver structure
--- a/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Err_test.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Err_test.m Wed Sep 07 14:17:30 2011 +0100 @@ -33,7 +33,7 @@ % Defining Automatic Transcription of Piano tune as Dictionary Learning % Problem -SMALL.Problem = generateAMT_Learning_Problem(); +SMALL.Problem = generateAMTProblem(); TPmax=0; for i=1:5 %% @@ -76,7 +76,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(i).D; - SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem); + SMALL.Problem.reconstruct = @(x) AMT_reconstruct(x, SMALL.Problem); %% % Initialising solver structure
--- a/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Sparsity_test.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Sparsity_test.m Wed Sep 07 14:17:30 2011 +0100 @@ -33,7 +33,7 @@ % Defining Automatic Transcription of Piano tune as Dictionary Learning % Problem -SMALL.Problem = generateAMT_Learning_Problem(); +SMALL.Problem = generateAMTProblem(); TPmax=0; @@ -73,7 +73,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(i).D; - SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem); + SMALL.Problem.reconstruct = @(x) AMT_reconstruct(x, SMALL.Problem); %% % Initialising solver structure
--- a/examples/Automatic Music Transcription/SMALL_AMT_SPAMS_test.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Automatic Music Transcription/SMALL_AMT_SPAMS_test.m Wed Sep 07 14:17:30 2011 +0100 @@ -33,7 +33,7 @@ % Defining Automatic Transcription of Piano tune as Dictionary Learning % Problem -SMALL.Problem = generateAMT_Learning_Problem(); +SMALL.Problem = generateAMTProblem(); TPmax=0; %% for i=1:10 @@ -77,7 +77,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(i).D; - SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem); + SMALL.Problem.reconstruct=@(x) AMT_reconstruct(x, SMALL.Problem); %%
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Wed Sep 07 14:17:30 2011 +0100 @@ -99,7 +99,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(1).D; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); %% % Initialising solver structure @@ -140,7 +140,7 @@ % Setting up reconstruction function SparseDict=0; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem, SparseDict); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem, SparseDict); % Initialising solver structure % Setting solver structure fields (toolbox, name, param, solution, @@ -217,7 +217,7 @@ % reconstructed and time) to zero values SMALL.Problem.A = SMALL.DL(3).D; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); SMALL.solver(3)=SMALL_init_solver;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLAvsTwoStepMOD.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,284 @@ +%% 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. + + + +% 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=4:4 +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) ImageDenoise_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 overcomplete DCT +% here +% Set SMALL.Problem.A dictionary to be oDCT (i.e. Problem.initdict - +% since initial dictionaruy is already set to be oDCT when generating the +% denoising problem + + +% 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('TwoStepDL', 'MOD', '', 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', 40,... + 'show_dict', 1); + +% 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) ImageDenoise_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)); + +%% +% In the b1 field all patches from the image are stored. For RLS-DLA we +% will first exclude all the patches that have l2 norm smaller then +% threshold and then take min(40000, number_of_remaining_patches) in +% ascending order as our training set (SMALL.Problem.b) + +X=SMALL.Problem.b1; +X_norm=sqrt(sum(X.^2, 1)); +[X_norm_sort, p]=sort(X_norm); +p1=p(X_norm_sort>Edata); +if size(p1,2)>40000 + p2 = randperm(size(p1,2)); + p2=sort(p2(1:40000)); + size(p2,2) + SMALL.Problem.b=X(:,p1(p2)); +else + size(p1,2) + SMALL.Problem.b=X(:,p1); + +end + +% Forgetting factor for RLS-DLA algorithm, in this case we are using +% fixed value + +lambda=0.9998 + +% Use Recursive Least Squares +% to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(3)=SMALL_init_DL(); + +% Defining fields needed for dictionary learning + +SMALL.DL(3).toolbox = 'SMALL'; +SMALL.DL(3).name = 'SMALL_rlsdla'; +SMALL.DL(3).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'forgettingMode', 'FIX',... + 'forgettingFactor', lambda,... + 'show_dict', 1000); + + +SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3)); + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.Problem.A = SMALL.DL(3).D; +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); + +SMALL.solver(3)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(3).toolbox='ompbox'; +SMALL.solver(3).name='omp2'; +SMALL.solver(3).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + + +SMALL.solver(3)=SMALL_solve(SMALL.Problem, SMALL.solver(3)); + +SMALL.solver(3).reconstructed.psnr + + +% 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).psnr.rlsdla=SMALL.solver(3).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).vmrse.rlsdla=SMALL.solver(3).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).ssim.rlsdla=SMALL.solver(3).reconstructed.ssim; + +results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; +results(noise_ind,im_num).time.rlsdla.time=SMALL.solver(3).time+SMALL.DL(3).time; +clear SMALL; +end +end +% save results.mat results
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsSPAMS.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsSPAMS.m Wed Sep 07 14:17:30 2011 +0100 @@ -97,7 +97,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(1).D; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); %% % Initialising solver structure @@ -175,7 +175,7 @@ % Setting up reconstruction function SparseDict=1; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem, SparseDict); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem, SparseDict); % Initialising solver structure % Setting solver structure fields (toolbox, name, param, solution, @@ -237,7 +237,7 @@ % Setting up reconstruction function -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); % Initialising solver structure % Setting solver structure fields (toolbox, name, param, solution,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsTwoStepKSVD.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,202 @@ +%% 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. +% - RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares +% Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on, +% vol.58, no.4, pp.2121-2130, April 2010 +% + + +% 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=1:1 +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 +% Ron Rubinstein implementation + +% 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) ImageDenoise_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 + +%% +% Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary +% Boris Mailhe ksvd update implentation omp is the same as with Rubinstein +% implementation + + +% 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('TwoStepDL', 'KSVD', '', 1); + + +% 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(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'show_dict', 1); + +% 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) ImageDenoise_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); + +clear SMALL; +end +end +
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_SPAMS_lambda.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_SPAMS_lambda.m Wed Sep 07 14:17:30 2011 +0100 @@ -90,7 +90,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(1).D; -SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); +SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); %% % Initialising solver structure
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_Training_size.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_Training_size.m Wed Sep 07 14:17:30 2011 +0100 @@ -100,7 +100,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(1).D; - SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); %% % Initialising solver structure @@ -167,7 +167,7 @@ % solver structures) SMALL.Problem.A = SMALL.DL(2).D; - SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); %% % Initialising solver structure
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_TwoStep_KSVD_MOD_OLS_Mailhe.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,309 @@ +%% 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. +% - RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares +% Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on, +% vol.58, no.4, pp.2121-2130, April 2010 +% + + +% 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=2:2 + +% 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); + + +%% +% Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary +% Boris Mailhe ksvd update implentation omp is the same as with Rubinstein +% implementation + + +% 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); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL('TwoStepDL', 'KSVD', '', 1); + + +% 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(... + 'solver', SMALL.solver(1),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'show_dict', 1); + +% 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) ImageDenoise_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(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +%% +% Use MOD Dictionary Learning Algorithm to Learn overcomplete dictionary +% Boris Mailhe MOD update implentation omp is the Ron Rubinstein +% implementation + + +% 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('TwoStepDL', 'MOD', '', 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,... + 'show_dict', 1); + +% 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) ImageDenoise_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)); +%% +% Use OLS Dictionary Learning Algorithm to Learn overcomplete dictionary +% Boris Mailhe ksvd update implentation omp is the Ron Rubinstein +% implementation + + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(3)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(3).toolbox='ompbox'; +SMALL.solver(3).name='omp2'; +SMALL.solver(3).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(3)=SMALL_init_DL('TwoStepDL', 'ols', '', 1); + + +% 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(3).param=struct(... + 'solver', SMALL.solver(3),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'learningRate', 0.1,... + 'show_dict', 1); + +% Learn the dictionary + +SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3)); + +% 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(3).D; +SMALL.Problem.reconstruct = @(x) ImageDenoise_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(3)=SMALL_solve(SMALL.Problem, SMALL.solver(3)); +%% +% Use Mailhe Dictionary Learning Algorithm to Learn overcomplete dictionary +% Boris Mailhe ksvd update implentation omp is the Ron Rubinstein +% implementation + + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(4)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(4).toolbox='ompbox'; +SMALL.solver(4).name='omp2'; +SMALL.solver(4).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(4)=SMALL_init_DL('TwoStepDL', 'mailhe', '', 1); + + +% 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(4).param=struct(... + 'solver', SMALL.solver(4),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'learningRate', 2,... + 'show_dict', 1); + +% Learn the dictionary + +SMALL.DL(4) = SMALL_learn(SMALL.Problem, SMALL.DL(4)); + +% 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(4).D; +SMALL.Problem.reconstruct = @(x) ImageDenoise_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(4)=SMALL_solve(SMALL.Problem, SMALL.solver(4)); + +%% show results %% + +SMALL_ImgDeNoiseResult(SMALL); + +%clear SMALL; +end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m Wed Sep 07 14:17:30 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 = generateAMTProblem('',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) AMT_reconstruct(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_pcgp'; + +% 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) AMT_reconstruct(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_AudioDenoise_DL_test_KSVDvsSPAMS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,128 @@ +%% DICTIONARY LEARNING FOR AUDIO DENOISING +% This file contains an example of how SMALLbox can be used to test different +% dictionary learning techniques in Audio Denoising problem. +% It calls generateAudioDenoiseProblem that will let you to choose audio file, +% add noise and use noisy audio to generate training set for dictionary +% learning. +% +% +% 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 Audio Denoising Problem as Dictionary Learning +% Problem + +SMALL.Problem = generateAudioDenoiseProblem('male01_8kHz',0.1,512,1/128,'','','',4); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(1)=SMALL_init_solver('MMbox', 'mm1', '', 1); + +% Defining the parameters needed for image denoising + +SMALL.solver(1).param=struct(... + 'lambda', 0.2,... + 'epsilon', 3*10^-4,... + 'iternum',10); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=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(1).param=struct(... + 'solver', SMALL.solver(1),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'iterDictUpdate', 10,... + 'epsDictUpdate', 10^-7,... + 'cvset',0,... + 'show_dict', 0); + +% 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) AudioDenoise_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(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +%% +%% +% % sparse coding using SPAMS online dictionary learning +% + +SMALL.DL(2)=SMALL_init_DL(); +SMALL.DL(2).toolbox = 'SPAMS'; +SMALL.DL(2).name = 'mexTrainDL'; +SMALL.DL(2).param=struct('D', SMALL.Problem.initdict, 'K', SMALL.Problem.p, 'lambda', 0.2, 'iter', 200, 'mode', 3, 'modeD', 0); + + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Defining Reconstruction function + +SMALL.Problem.A = SMALL.DL(2).D; + + +%% +% Initialising solver structure +% Setting toolbox, name, param, solution, reconstructed and time to zero values + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for sparse representation + +SMALL.solver(2).toolbox='ompbox'; +SMALL.solver(2).name='omp2'; +SMALL.solver(2).param=struct(... + 'epsilon',0.2,... + 'maxatoms', 128); +% Represent Training set in the learned dictionary + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + + + + +%% +% Plot results and save midi files + +% show results % + + +SMALL_AudioDeNoiseResult(SMALL); + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_ImgDenoise_DL_test_KSVDvsMajorizationMinimization.m Wed Sep 07 14:17:30 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) ImageDenoise_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) ImageDenoise_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/examples/Pierre Villars/Pierre_Villars_Example.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/Pierre Villars/Pierre_Villars_Example.m Wed Sep 07 14:17:30 2011 +0100 @@ -23,7 +23,7 @@ % Defining the Problem structure -SMALL.Problem = generatePierre_Problem(); +SMALL.Problem = generatePierreProblem(); % Show original image and image that is used as a dictionary figure('Name', 'Original and Dictionary Image');
--- a/examples/SMALL_solver_test.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/SMALL_solver_test.m Wed Sep 07 14:17:30 2011 +0100 @@ -81,7 +81,7 @@ % SMALL Conjugate Gradient test SMALL.solver(i)=SMALL_init_solver; SMALL.solver(i).toolbox='SMALL'; -SMALL.solver(i).name='SMALL_cgp'; +SMALL.solver(i).name='SMALL_pcgp'; % In the following string all parameters except matrix, measurement vector % and size of solution need to be specified. If you are not sure which
--- a/examples/SMALL_solver_test_Audio.m Tue Jul 26 16:02:59 2011 +0100 +++ b/examples/SMALL_solver_test_Audio.m Wed Sep 07 14:17:30 2011 +0100 @@ -61,7 +61,7 @@ SMALL.solver(i)=SMALL_init_solver; SMALL.solver(i).toolbox='SMALL'; -SMALL.solver(i).name='SMALL_cgp'; +SMALL.solver(i).name='SMALL_pcgp'; % In the following string all parameters except matrix, measurement vector % and size of solution need to be specified. If you are not sure which
--- a/solvers/CVX_add_const_Audio_declipping.m Tue Jul 26 16:02:59 2011 +0100 +++ b/solvers/CVX_add_const_Audio_declipping.m Wed Sep 07 14:17:30 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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_pcgp.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,113 @@ +function [A, resF]=SMALL_pcgp(Dict,X, m, maxNumCoef, errorGoal, varargin) +%% Partial Conjugate Gradient Pursuit Multiple vectors sparse representation +% +% Sparse coding of a group of signals based on a given +% dictionary and specified number of atoms to use. +% input arguments: Dict - the dictionary +% X - the signals to represent +% m - number of atoms in Dictionary +% errorGoal - the maximal allowed representation error for +% each signal. +% +% optional: if Dict is function handle then Transpose Dictionary +% handle needs to be specified. +% +% output arguments: A - sparse coefficient matrix. + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2009 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. +% +%% + +% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers +if isa(Dict,'float') + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isobject(Dict) + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isa(Dict,'function_handle') + try + DictT=varargin{1}; + if isa(DictT,'function_handle'); + D=Dict; + Dt=DictT; + else + error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. '); + end + catch + error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.'); + end +else + error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.'); +end +%% +%positivity=1; +[n,P]=size(X); + + + +A = sparse(m,size(X,2)); + +for k=1:1:P, + + x = X(:,k); + residual = x; + indx =[]; + j=0; + + + currResNorm = residual'*residual/n; + errGoal=errorGoal*currResNorm; + a = zeros(m,1); + p = zeros(m,1); + + while j < maxNumCoef, + + j = j+1; + dir=Dt(residual); + if exist('positivity','var')&&(positivity==1) + [tmp__, pos]=max(dir); + else + [tmp__, pos]=max(abs(dir)); + end + indx(j)=pos; + + p(indx)=dir(indx); + Dp=D(p); + pDDp=Dp'*Dp; + + if (j==1) + beta=0; + else + beta=-Dp'*residual/pDDp; + end + + alfa=residual'*Dp/pDDp; + a=a+alfa*p; + p(indx)=dir(indx)+beta*p(indx); + + residual=residual-alfa*Dp; + + currResNorm = residual'*residual/n; + if currResNorm<errGoal + fprintf('\nFound exact representation! \n'); + break; + end + end; + if (~isempty(indx)) + resF(k)=currResNorm; + A(indx,k)=a(indx); + end +end; +return; + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/AlgebraicPursuit.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,264 @@ +function [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, varargin) +% ========================================================================= +% Algebraic Pursuit algorithms - Beta Version +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods +% on sparse recovery of linear inverse systems. In particular, let +% Phi : M x N real matrix (M < N) +% x* : N x 1 K-sparse data vector +% n : N x 1 additive noise vector +% then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS +% solve the following minimization problem +% minimize ||y - Phi x||^2 subject to x is K-sparse vector. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% INPUT ARGUMENTS: +% y M x 1 undersampled measurement vector. +% Phi M x N regression matrix. +% K Sparsity of underlying vector x* or desired +% sparsity of solution. +% varargin Set of parameters. These are: +% +% memory,... Memory (momentum) of proposed algorithm. +% Possible values are 0,1,'infty' for memoryless, +% one memory and infinity memory ALPS, +% respectively. Default value: memory = 0. +% tol,... Early stopping tolerance. Default value: tol = +% 1-e5. +% ALPSiters,... Maximum number of algorithm iterations. Default +% value: 300. +% mode, ... According to [1], possible values are +% [0,1,2,4,5,6]. This value comes from the binary +% representation of the parameters: +% (solveNewtob, gradientDescentx, solveNewtonx), +% which are explained next. Default value = 0. +% mu,... Variable that controls the step size selection. +% When mu = 0, step size is computed adaptively +% per iteration. Default value: mu = 0. +% tau,... Variable that controls the momentum in +% non-memoryless case. Ignored in memoryless +% case. User can insert as value a function handle on tau. +% Description given below. Default value: tau = -1. +% CGiters,... Maximum iterations for Conjugate-Gradients method. +% CGtol,... Tolerance variable for Conjugate-Gradients method. +% verbose verbose = 1 prints out execution infromation. +% ========================================================================= +% DESCRIPTION OF MODE VARIABLE: +% solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a +% Newton system restricted to a sparse support. +% It is implemented via conjugate gradients. If +% solveNewtonb == 0: Step size selection as +% described in eqs. (12) and (13) in [1]. Default +% value: solveNewtonb = 0. For infty memory case, +% solveNewtonb = 0. +% gradientDescentx,... If gradientDescentx == 1: single gradient +% update of x_{i+1} restricted ot its support +% with line search. Default value: +% gradientDescentx = 1. +% solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding +% Pursuit (c.f. Simon Foucart, "Hard Thresholding +% Pursuit," preprint, 2010). Default vale: +% solveNewtonx = 0. +% ========================================================================= +% DESCRIPTION OF SPECIAL TAU VALUES: +% tau = -1,... Tau computed as the minimized of the objective +% function: +% +% <u - Ax_i, Phi(x_i - x_{i-1})> +% tau = -----------------------------. +% <u - Ax_i, Phi(x_i - x_{i-1})> +% +% tau = 0,... Follows FISTA weight configuration at each +% iteration. For more information, please read "A +% Fast Iterative Shrinkage-Thresholding Algorithm +% for Linear Inverse Problems", written by A. +% Beck and M. Teboulle, SIAM J. Imaging Sciences, +% vol. 2, no. 1, 2009. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x_hat N x 1 recovered K-sparse vector. +% numiter Number of iterations executed. +% time Execution time in seconds. +% x_path Keeps a series of computed N x 1 K-sparse vectors +% until the end of the iterative process. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. +% Email: jrom@acm.caltech.edu +% ========================================================================= +% This work was supported in part by the European Commission under Grant +% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like +% to acknowledge Rice University for his Faculty Fellowship. +% ========================================================================= + +x_hat = 0; +time = 0; +x_path = 0; + +params.memory = 0; +params.tol = 1e-5; +params.ALPSiters = 300; +params.mode = 0; +params.tau = 1/2; +params.memory = 0; +params.cg_maxiter = 50; +params.cg_tol = 1e-8; +params.useCG = 0; +params.verbose = 0; +params.mu = 0; + +for in_arg = 1:2:length(varargin) + if (strcmp(varargin{in_arg}, 'memory')) + params.memory = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'mode')) + params.mode = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'tolerance')) + params.tol = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'ALPSiterations')) + params.ALPSiters = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'tau')) + params.tau = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'mu')) + params.mu = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'useCG')) + params.useCG = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'CGiters')) + params.cg_maxiter = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'CGtol')) + params.cg_tol = varargin{in_arg+1}; + end + if (strcmp(varargin{in_arg}, 'verbose')) + params.verbose = varargin{in_arg+1}; + end; +end; + +%% Check input + +if (ischar(params.memory)) + if ~(sum(params.memory == 'infty') == 5) + disp('ERROR: Memory variable takes positive integer values or "infty" value.'); + numiter = -1; + return; + end; +elseif (params.memory < 0 || ~isinteger(uint8(params.memory))) + disp('ERROR: Memory variable takes positive integer values or "infty" value.'); + numiter = -1; + return; +end; + +if (params.mode ~= 0 && params.mode ~= 1 && params.mode ~= 2 && params.mode ~= 4 && params.mode ~= 5 && params.mode ~= 6) + disp('ERROR: Mode variable takes values from range [0,1,2,4,5,6].'); + numiter = -1; + return; +end; + +if (params.tol < 0 || ~isnumeric(params.tol)) + disp('ERROR: tolerance variable must be positive.'); + numiter = -1; + return; +end; + + +if (params.mode == 0 || params.mode == 2) + params.useCG = 0; +end; + +y = y(:); + +M_y = length(y); +[M_Phi, tmpArg] = size(Phi); + +if (params.mode == 4 || params.mode == 5 || params.mode == 6) + if (M_Phi < 2*K) + disp('WARNING: Newton system may be ill-conditioned... Press any key to continue'); + pause; + end; +end; + +[params.solveNewtonb, params.gradientDescentx, params.solveNewtonx] = configuration(params.mode); + +if (M_y ~= M_Phi) + error('Inconsistent sampling matrix...'); +end; + +switch params.memory + case 0, + tic; + [x_hat, numiter, x_path] = zero_ALPS(y, Phi, K, params); + time = toc; + if (params.verbose == 1) + str = sprintf('%d-ALPS(%d) time: %s', params.memory, params.mode, num2str(time)); + disp(str); + str = sprintf('Number of iterations: %d', numiter); + disp(str); + end; + case 1, + tic; + [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params); + time = toc; + if (params.verbose == 1) + if (isa(params.tau,'float')) + if (params.tau == 0) + str = sprintf('%d-ALPS(%d) - fista tau - time: %s', params.memory, params.mode, num2str(time)); + elseif (params.tau == -1) + str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time)); + else + str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time)); + end; + elseif (isa(params.tau, 'function_handle')) + str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time)); + end; + disp(str); + str = sprintf('Number of iterations: %d', numiter); + disp(str); + end; + case 'infty', + tic; + [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params); + time = toc; + if (params.verbose == 1) + if (isa(params.tau,'float')) + if (params.tau == -1) + str = sprintf('infty-ALPS(%d) - optimized tau - time: %s', params.mode, num2str(time)); + else + str = sprintf('infty-ALPS(%d) - tau = %.3f - time: %s', params.mode, params.tau, num2str(time)); + end; + elseif (isa(params.tau,'function_handle')) + str = sprintf('infty-ALPS(%d) - user derfined tau - time: %s', params.mode, num2str(time)); + end; + disp(str); + str = sprintf('Number of iterations: %d', numiter); + disp(str); + end; + otherwise + tic; + [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params); + time = toc; + if (params.verbose == 1) + if (isa(params.tau,'float')) + if (params.tau == -1) + str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time)); + else + str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time)); + end; + elseif (isa(params.tau,'function_handle')) + str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time)); + end; + disp(str); + str = sprintf('Number of iterations: %d', numiter); + disp(str); + end; +end; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/cgsolve.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,76 @@ +% cgsolve.m +% +% Solve a symmetric positive definite system Ax = b via conjugate gradients. +% +% Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) +% +% A - Either an NxN matrix, or a function handle. +% +% b - N vector +% +% tol - Desired precision. Algorithm terminates when +% norm(Ax-b)/norm(b) < tol . +% +% maxiter - Maximum number of iterations. +% +% verbose - If 0, do not print out progress messages. +% Default = 1. +% +% Written by: Justin Romberg, Caltech +% Email: jrom@acm.caltech.edu +% Created: October 2005 +% + +function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose) + +if (nargin < 5), verbose = 1; end + +implicit = isa(A,'function_handle'); + +x = zeros(length(b),1); +r = b; +d = r; +delta = r'*r; +delta0 = b'*b; +numiter = 0; +bestx = x; +bestres = sqrt(delta/delta0); +while ((numiter < maxiter) & (delta > tol^2*delta0)) + + % q = A*d + if (implicit), q = A(d); else, q = A*d; end + + alpha = delta/(d'*q); + x = x + alpha*d; + + if (mod(numiter+1,50) == 0) + % r = b - Aux*x + if (implicit), r = b - A(x); else, r = b - A*x; end + else + r = r - alpha*q; + end + + deltaold = delta; + delta = r'*r; + beta = delta/deltaold; + d = r + beta*d; + numiter = numiter + 1; + if (sqrt(delta/delta0) < bestres) + bestx = x; + bestres = sqrt(delta/delta0); + end + + if ((verbose) & (mod(numiter,50)==0)) + disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ... + numiter, bestres, sqrt(delta/delta0))); + end + +end + +if (verbose) + disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres)); +end +x = bestx; +res = bestres; +iter = numiter; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/configuration.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,28 @@ +function [solveNewtonb, gradientDescentx, solveNewtonx] = configuration(mode) +% Function that maps 'mode' value into (SolveNewtonb, GradientDescentx, SolveNewtonx) configuration + +if (mode == 0) + solveNewtonb = 0; + gradientDescentx = 0; + solveNewtonx = 0; +elseif (mode == 1) + solveNewtonb = 0; + gradientDescentx = 0; + solveNewtonx = 1; +elseif (mode == 2) + solveNewtonb = 0; + gradientDescentx = 1; + solveNewtonx = 0; +elseif (mode == 4) + solveNewtonb = 1; + gradientDescentx = 0; + solveNewtonx = 0; +elseif (mode == 5) + solveNewtonb = 1; + gradientDescentx = 0; + solveNewtonx = 1; +elseif (mode == 6) + solveNewtonb = 1; + gradientDescentx = 1; + solveNewtonx = 0; +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/infty_ALPS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,256 @@ +function [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params) +% ========================================================================= +% infty-ALPS(#) algorithm - Beta Version +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithm with infty-memory acceleration. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% INPUT ARGUMENTS: +% y M x 1 undersampled measurement vector. +% Phi M x N regression matrix. +% K Sparsity of underlying vector x* or desired +% sparsity of solution. +% params Structure of parameters. These are: +% +% tol,... Early stopping tolerance. Default value: tol = +% 1-e5 +% ALPSiters,... Maximum number of algorithm iterations. Default +% value: 300. +% solveNewtonb,... Value: solveNewtonb = 0. Not used in infty +% methods. +% gradientDescentx,... If gradientDescentx == 1: single gradient +% update of x_{i+1} restricted ot its support with +% line search. Default value: gradientDescentx = +% 1. +% solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit +% (c.f. Simon Foucart, "Hard Thresholding Pursuit," +% preprint, 2010). Default vale: solveNewtonx = 0 +% tau,... Variable that controls the momentum in +% non-memoryless case. Ignored in memoryless +% case. Default value: tau = 1/2. +% Special cases: +% - tau = -1: momentum step size is automatically +% optimized in every step. +% - tau as a function handle: user defined +% behavior of tau momentum term. +% mu,... Variable that controls the step size selection. +% When mu = 0, step size is computed adaptively +% per iteration. Default value: mu = 0. +% cg_maxiter,... Maximum iterations for Conjugate-Gradients method. +% cg_tol Tolerance variable for Conjugate-Gradients method. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x_hat N x 1 recovered K-sparse vector. +% numiter Number of iterations executed. +% x_path Keeps a series of computed N x 1 K-sparse vectors +% until the end of the iterative process. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. +% Email: jrom@acm.caltech.edu +% ========================================================================= +% This work was supported in part by the European Commission under Grant +% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like +% to acknowledge Rice University for his Faculty Fellowship. +% ========================================================================= + +[~,N] = size(Phi); + +%% Initialize transpose of measurement matrix + +Phi_t = Phi'; + +%% Initialize to zero vector +x_cur = zeros(N,1); +y_cur = zeros(N,1); +X_i = []; + +x_path = zeros(N, params.ALPSiters); + +%% CG params +if (params.solveNewtonx == 1 || params.solveNewtonb == 1) + cg_verbose = 0; + cg_A = Phi_t*Phi; + cg_b = Phi_t*y; +end; + +%% Determine momentum step size selection strategy +optimizeTau = 0; +function_tau = 0; + +if (isa(params.tau,'float')) + if (params.tau == -1) + optimizeTau = 1; + end; +elseif (isa(params.tau, 'function_handle')) + function_tau = 1; +end; + +%% Determine step size selection strategy +function_mu = 0; +adaptive_mu = 0; + +if (isa(params.mu,'float')) + function_mu = 0; + if (params.mu == 0) + adaptive_mu = 1; + else + adaptive_mu = 0; + end; +elseif (isa(params.mu,'function_handle')) + function_mu = 1; +end; + +%% Help variables +complementary_Xi = ones(N,1); +setXi = zeros(N,1); +setYi = zeros(N,1); + +i = 1; +%% infty-ALPS(#) +while (i <= params.ALPSiters) + x_path(:,i) = x_cur; + x_prev = x_cur; + + % Compute the residual + if (i == 1) + res = y; + else + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + end; + + % Compute the derivative + der = Phi_t*res; + + % Determine S_i set via eq. (11) + complementary_Xi(X_i) = 0; + [~, ind_der] = sort(abs(der).*complementary_Xi, 'descend'); + complementary_Xi(X_i) = 1; + S_i = [X_i; ind_der(1:K)]; + + ider = der(S_i); + + setder = zeros(N,1); + setder(S_i) = 1; + if (adaptive_mu) + % Step size selection via eq. (12) and eq. (13) + Pder = Phi(:,S_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + end; + + iy_cur = y_cur.*setYi; + if (~function_tau) % If tau is not a function handle... + if (optimizeTau) % Compute optimized tau + + % tau = argmin || u - Phi(x_i + y_i) || + % = <Phi*y_i, u - Phi(x_i - mu/2 * grad_Si f(xi))>/||Phi*y_i||^2 + + if (i == 1) + params.tau = 0; + else + % u - Phi*(x_i - mu/2 grad_Si f(xi)) = u - Phi*b + if (adaptive_mu) + b = x_cur(S_i) + mu_bar*ider; % Non-zero elems: S_i + elseif (function_mu) + b = x_cur(S_i) + params.mu(i)*ider; + else b = x_cur(S_i) + params.mu*ider; + end; + + y_Phi_b = y - Phi(:,S_i)*b; + Phi_y_prev = Phi(:,Y_i)*y_cur(Y_i); % Phi * y_i + params.tau = y_Phi_b'*Phi_y_prev/(Phi_y_prev'*Phi_y_prev); + end; + + if (adaptive_mu) + y_cur = params.tau*iy_cur + mu_bar*der.*setder; + elseif (function_mu) + y_cur = params.tau*iy_cur + params.mu(i)*der.*setder; + else y_cur = params.tau*iy_cur + params.mu*der.*setder; + end; + + Y_i = ne(y_cur,0); + setYi = zeros(N,1); + setYi(Y_i) = 1; + else % Tau fixed and user-defined + if (adaptive_mu) + y_cur = params.tau*iy_cur + mu_bar*der.*setder; + elseif (function_mu) + y_cur = params.tau*iy_cur + params.mu(i)*der.*setder; + else y_cur = params.tau*iy_cur + params.mu*der.*setder; + end; + + Y_i = ne(y_cur,0); + setYi = zeros(N,1); + setYi(Y_i) = 1; + end; + else + if (adaptive_mu) + y_cur = params.tau(i)*iy_cur + mu_bar*der.*setder; + elseif (function_mu) + y_cur = params.tau(i)*iy_cur + params.mu(i)*der.*setder; + else y_cur = params.tau(i)*iy_cur + params.mu*der.*setder; + end; + + Y_i = ne(y_cur,0); + setYi = zeros(N,1); + setYi(Y_i) = 1; + end; + + % Hard-threshold b and compute X_{i+1} + set_Xi_Yi = setXi + setYi; + ind_Xi_Yi = find(set_Xi_Yi > 0); + z = x_cur(ind_Xi_Yi) + y_cur(ind_Xi_Yi); + [~, ind_z] = sort(abs(z), 'descend'); + x_cur = zeros(N,1); + x_cur(ind_Xi_Yi(ind_z(1:K))) = z(ind_z(1:K)); + X_i = ind_Xi_Yi(ind_z(1:K)); + + setXi = zeros(N,1); + setXi(X_i) = 1; + + if (params.gradientDescentx == 1) + % Calculate gradient of estimated vector x_cur + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + der = Phi_t*res; + + ider = der(X_i); + + if (adaptive_mu) + Pder = Phi(:,X_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + x_cur(X_i) = x_cur(X_i) + mu_bar*ider; + elseif (function_mu) + x_cur = x_cur(X_i) + params.mu(i)*ider; + else x_cur = x_cur(X_i) + params.mu*ider; + end; + elseif (params.solveNewtonx == 1) + % Similar to HTP + if (params.useCG == 1) + [v, ~, ~] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + v = cg_A(X_i,X_i)\cg_b(X_i); + end; + x_cur(X_i) = v; + end; + + % Test stopping criterion + if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur)) + break; + end; + i = i + 1; +end; + +x_hat = x_cur; +numiter = i; + +if (i > params.ALPSiters) + x_path = x_path(:,1:numiter-1); +else + x_path = x_path(:,1:numiter); +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/l_ALPS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,267 @@ +function [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params) +% ========================================================================= +% l-ALPS(#) algorithm - Beta Version +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithm with l-memory acceleration. +% +% Detailed discussion on the algorithm can be found in section III in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% INPUT ARGUMENTS: +% y M x 1 undersampled measurement vector. +% Phi M x N regression matrix. +% K Sparsity of underlying vector x* or desired +% sparsity of solution. +% params Structure of parameters. These are: +% +% tol,... Early stopping tolerance. Default value: tol = +% 1-e5 +% ALPSiters,... Maximum number of algorithm iterations. Default +% value: 300. +% solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a +% Newton system restricted to a sparse support. +% It is implemented via conjugate gradients. +% If solveNewtonb == 0: Step size selection as described +% in eqs. (12) and (13) in [1]. +% Default value: solveNewtonb = 0. +% gradientDescentx,... If gradientDescentx == 1: single gradient +% update of x_{i+1} restricted ot its support with +% line search. Default value: gradientDescentx = +% 1. +% solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit +% (c.f. Simon Foucart, "Hard Thresholding Pursuit," +% preprint, 2010). Default vale: solveNewtonx = 0 +% tau,... Variable that controls the momentum in +% non-memoryless case. Ignored in memoryless +% case. Default value: tau = 1/2. +% Special cases: +% - tau = -1: momentum step size is automatically +% optimized in every step. +% - tau as a function handle: user defined +% behavior of tau momentum term. +% mu,... Variable that controls the step size selection. +% When mu = 0, step size is computed adaptively +% per iteration. Default value: mu = 0. +% cg_maxiter,... Maximum iterations for Conjugate-Gradients method. +% cg_tol Tolerance variable for Conjugate-Gradients method. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x_hat N x 1 recovered K-sparse vector. +% numiter Number of iterations executed. +% x_path Keeps a series of computed N x 1 K-sparse vectors +% until the end of the iterative process. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. +% Email: jrom@acm.caltech.edu +% ========================================================================= +% This work was supported in part by the European Commission under Grant +% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like +% to acknowledge Rice University for his Faculty Fellowship. +% ========================================================================= + +[M, N] = size(Phi); + +%% Initialize transpose of measurement matrix + +Phi_t = Phi'; + +%% Initialize to zero vector +y_cur = zeros(N,1); +Y_i = []; + +x_path = zeros(N, params.ALPSiters); + +%% CG params +if (params.solveNewtonx == 1 || params.solveNewtonb == 1) + cg_verbose = 0; + cg_A = Phi_t*Phi; + cg_b = Phi_t*y; +end; + +%% Determine momentum step size selection strategy +optimizeTau = 0; +function_tau = 0; + +if (isa(params.tau,'float')) + if (params.tau == -1) + optimizeTau = 1; + end; +elseif (isa(params.tau, 'function_handle')) + function_tau = 1; +end; + +%% Determine step size selection strategy +function_mu = 0; +adaptive_mu = 0; + +if (isa(params.mu,'float')) + function_mu = 0; + if (params.mu == 0) + adaptive_mu = 1; + else + adaptive_mu = 0; + end; +elseif (isa(params.mu,'function_handle')) + function_mu = 1; +end; + +%% Help variables +complementary_Yi = ones(N,1); +x = zeros(N,params.memory + 1); +Phi_x = zeros(M,params.memory + 1); +tau_candidates = zeros(params.memory,1); +y_candidates = zeros(N,params.memory); +residual_energy = zeros(params.memory,1); + +i = 1; +%% l-ALPS(#) +while (i <= params.ALPSiters) + x_path(:,i) = x(:,end); + + for k = 1:params.memory + x(:,k) = x(:,k+1); + Phi_x(:,k) = Phi_x(:,k+1); + end; + + % Compute the residual + if (i == 1) + res = y; + else + Phi_y_cur = Phi(:,Y_i)*y_cur(Y_i); + res = y - Phi_y_cur; + end; + + % Compute the derivative + der = Phi_t*res; + + % Determine S_i via eq. (11) (change of variable from x_i to y_i) + complementary_Yi(Y_i) = 0; + [tmpArg, ind_der] = sort(abs(der).*complementary_Yi, 'descend'); + complementary_Yi(Y_i) = 1; + S_i = [Y_i; ind_der(1:K)]; + + ider = der(S_i); + if (params.solveNewtonb == 1) + % Compute least squares solution of the system A*y = (A*A)x using CG + if (params.useCG == 1) + [b, tmpArg, tmpArg] = cgsolve(cg_A(S_i, S_i), cg_b(S_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + b = cg_A(S_i,S_i)\cg_b(S_i); + end; + else + % Step size selection via eq. (12) and eq. (13) (change of variable from x_i to y_i) + if (adaptive_mu) + Pder = Phi(:,S_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + b = y_cur(S_i) + mu_bar*ider; + elseif (function_mu) + b = y_cur(S_i) + params.mu(i)*ider; + else b = y_cur(S_i) + params.mu*ider; + end; + end; + + % Hard-threshold b and compute X_{i+1} + [tmpArg, ind_b] = sort(abs(b), 'descend'); + x(:,end) = zeros(N,1); + x(S_i(ind_b(1:K)),end) = b(ind_b(1:K)); + X_i = S_i(ind_b(1:K)); + + if (params.gradientDescentx == 1) + % Calculate gradient of estimated current x vector + Phi_x_cur = Phi(:,X_i)*x(X_i,end); + res = y - Phi_x_cur; + der = Phi_t*res; + ider = der(X_i); + + if (adaptive_mu) + Pder = Phi(:,X_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + x(X_i,end) = x(X_i,end) + mu_bar*ider; + elseif (function_mu) + x(X_i,end) = x(X_i,end) + params.mu(i)*ider; + else x(X_i,end) = x(X_i,end) + params.mu*ider; + end; + elseif (params.solveNewtonx == 1) + % Similar to HTP + if (params.useCG == 1) + [v, tmpArg, tmpArg] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + v = cg_A(X_i, X_i)\cg_b(X_i); + end; + x(X_i,end) = v; + end; + + Phi_x(:,end) = Phi(:,X_i)*x(X_i,end); + res = y - Phi_x(:,end); + if (~function_tau) % If tau is not a function handle... + if (optimizeTau) % Compute optimized tau + if (i > params.memory) + + % tau = argmin ||u - Phi*y_{i+1}|| + % = <res, Phi*(x_cur - x_prev)>/||Phi*(x_cur - x_prev)||^2 + + for k = 1:params.memory + Phi_diff = Phi_x(:,end) - Phi_x(:,k); + tau_candidates(k) = res'*Phi_diff/(Phi_diff'*Phi_diff); + y_candidates(:,k) = x(:,end) + tau_candidates(k)*(x(:,end) - x(:,k)); + residual_energy(k) = norm(y - Phi*y_candidates(:,k)); + end; + + [tmpArg,min_ind] = min(residual_energy); + y_cur = y_candidates(:,min_ind); + Y_i = find(ne(y_cur,0)); + else + Phi_diff = Phi_x(:,end) - Phi_x(:,end-1); + params.tau = res'*Phi_diff/(Phi_diff'*Phi_diff); + y_cur = x(:,end) + params.tau*(x(:,end) - x(:,end-1)); + Y_i = find(ne(y_cur,0)); + end; + else + if (i > params.memory) + for k = 1:params.memory + y_candidates(:,k) = x(:,end) + params.tau*(x(:,end) - x(:,k)); + residual_energy(k) = norm(y - Phi*y_candidates(:,k)); + end; + + [tmpArg,min_ind] = min(residual_energy); + y_cur = y_candidates(:,min_ind); + Y_i = find(ne(y_cur,0)); + else + y_cur = x(:,end) + params.tau*(x(:,end) - x(:,end-1)); + Y_i = find(ne(y_cur,0)); + end; + end; + else + if (i > params.memory) + for k = 1:params.memory + y_candidates(:,k) = x(:,end) + params.tau(i)*(x(:,end) - x(:,k)); + residual_energy(k) = norm(y - Phi*y_candidates(:,k))^2; + end; + + [tmpArg,min_ind] = min(residual_energy); + y_cur = y_candidates(:,min_ind); + Y_i = find(ne(y_cur,0)); + else + y_cur = x(:,end) + params.tau(i)*(x(:,end) - x(:,end-1)); + Y_i = find(ne(y_cur,0)); + end; + end; + + % Test stopping criterion + if (i > 1) && (norm(x(:,end) - x(:,end-1)) < params.tol*norm(x(:,end))) + break; + end; + i = i + 1; +end; + +x_hat = x(:,end); +numiter = i; + +if (i > params.ALPSiters) + x_path = x_path(:,1:numiter-1); +else + x_path = x_path(:,1:numiter); +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/one_ALPS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,252 @@ +function [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params) +% ========================================================================= +% 1-ALPS(#) algorithm - Beta Version +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithm with 1-memory acceleration. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% INPUT ARGUMENTS: +% y M x 1 undersampled measurement vector. +% Phi M x N regression matrix. +% K Sparsity of underlying vector x* or desired +% sparsity of solution. +% params Structure of parameters. These are: +% +% tol,... Early stopping tolerance. Default value: tol = +% 1-e5. +% ALPSiters,... Maximum number of algorithm iterations. Default +% value: 300. +% solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a +% Newton system restricted to a sparse support. +% It is implemented via conjugate gradients. +% If solveNewtonb == 0: Step size selection as described +% in eqs. (12) and (13) in [1]. +% Default value: solveNewtonb = 0. +% gradientDescentx,... If gradientDescentx == 1: single gradient +% update of x_{i+1} restricted ot its support with +% line search. Default value: gradientDescentx = +% 1. +% solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit +% (c.f. Simon Foucart, "Hard Thresholding Pursuit," +% preprint, 2010). Default vale: solveNewtonx = 0. +% tau,... Variable that controls the momentum in +% non-memoryless case. Ignored in memoryless +% case. Default value: tau = 1/2. +% Special cases: +% - tau = 0: momentum step size selection is +% driven by the following formulas: +% a_1 = 1; +% a_{i+1} = (1+\sqrt(1+4a_i^2)/2; +% tau = (a_i - 1)/(a_{i+1}); +% described in [2] "A fast iterative +% shrinkage-thresholding algorithm for linear +% inverse problems", Beck A., and Teboulle M. +% - tau = -1: momentum step size is automatically +% optimized in every step. +% - tau as a function handle: user defined +% behavior of tau momentum term. +% mu,... Variable that controls the step size selection. +% When mu = 0, step size is computed adaptively +% per iteration. Default value: mu = 0. +% cg_maxiter,... Maximum iterations for Conjugate-Gradients method. +% cg_tol Tolerance variable for Conjugate-Gradients method. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x_hat N x 1 recovered K-sparse vector. +% numiter Number of iterations executed. +% x_path Keeps a series of computed N x 1 K-sparse vectors +% until the end of the iterative process. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. +% Email: jrom@acm.caltech.edu +% ========================================================================= +% This work was supported in part by the European Commission under Grant +% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like +% to acknowledge Rice University for his Faculty Fellowship. +% ========================================================================= + +[M,N] = size(Phi); + +%% Initialize transpose of measurement matrix + +Phi_t = Phi'; + +%% Initialize to zero vector +x_cur = zeros(N,1); +y_cur = zeros(N,1); +Phi_x_cur = zeros(M,1); +Y_i = []; + +x_path = zeros(N, params.ALPSiters); + +%% CG params +if (params.solveNewtonx == 1 || params.solveNewtonb == 1) + cg_verbose = 0; + cg_A = Phi_t*Phi; + cg_b = Phi_t*y; +end; + +%% Determine momentum step size selection strategy +fista = 0; +optimizeTau = 0; +a_prev = 1; +function_tau = 0; + +if (isa(params.tau,'float')) + function_tau = 0; + if (params.tau == 0) + fista = 1; + optimizeTau = 0; + elseif (params.tau == -1) + optimizeTau = 1; + fista = 0; + end; +elseif (isa(params.tau, 'function_handle')) + function_tau = 1; +end; + +%% Determine step size selection strategy +function_mu = 0; +adaptive_mu = 0; + +if (isa(params.mu,'float')) + function_mu = 0; + if (params.mu == 0) + adaptive_mu = 1; + else + adaptive_mu = 0; + end; +elseif (isa(params.mu,'function_handle')) + function_mu = 1; +end; + +%% Help variables +complementary_Yi = ones(N,1); + +i = 1; +%% 1-ALPS(#) +while (i <= params.ALPSiters) + x_path(:,i) = x_cur; + x_prev = x_cur; + + % Compute the residual + if (i == 1) + res = y; + % Compute the derivative + der = Phi_t*res; + else + % Compute the derivative + if (optimizeTau) + res = y - Phi_x_cur - params.tau*Phi_diff; + else + res = y - Phi(:,Y_i)*y_cur(Y_i); + end; + der = Phi_t*res; + end; + + Phi_x_prev = Phi_x_cur; + + % Determine S_i set via eq. (11) (change of variable from x_i to y_i) + complementary_Yi(Y_i) = 0; + [tmpArg, ind_der] = sort(abs(der).*complementary_Yi, 'descend'); + complementary_Yi(Y_i) = 1; + S_i = [Y_i; ind_der(1:K)]; + + ider = der(S_i); + if (params.solveNewtonb == 1) + % Compute least squares solution of the system A*y = (A*A)x using CG + if (params.useCG == 1) + [b, tmpArg, tmpArg] = cgsolve(cg_A(S_i, S_i), cg_b(S_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + b = cg_A(S_i,S_i)\cg_b(S_i); + end; + else + % Step size selection via eq. (12) and eq. (13) (change of variable from x_i to y_i) + if (adaptive_mu) + Pder = Phi(:,S_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + b = y_cur(S_i) + (mu_bar)*ider; + elseif (function_mu) + b = y_cur(S_i) + params.mu(i)*ider; + else b = y_cur(S_i) + params.mu*ider; + end; + end; + + % Hard-threshold b and compute X_{i+1} + [tmpArg, ind_b] = sort(abs(b), 'descend'); + x_cur = zeros(N,1); + x_cur(S_i(ind_b(1:K))) = b(ind_b(1:K)); + X_i = S_i(ind_b(1:K)); + + if (params.gradientDescentx == 1) + % Calculate gradient of estimated vector x_cur + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + der = Phi_t*res; + ider = der(X_i); + + if (adaptive_mu) + Pder = Phi(:,X_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + x_cur(X_i) = x_cur(X_i) + mu_bar*ider; + elseif (function_mu) + x_cur(X_i) = x_cur(X_i) + params.mu(i)*ider; + else x_cur(X_i) = x_cur(X_i) + params.mu*ider; + end; + elseif (params.solveNewtonx == 1) + % Similar to HTP + if (params.useCG == 1) + [v, tmpArg, tmpArg] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + v = cg_A(X_i, X_i)\cg_b(X_i); + end; + x_cur(X_i) = v; + end; + + if (~function_tau) % If tau is not a function handle... + if (fista) % Fista configuration + a_cur = (1 + sqrt(1 + 4*a_prev^2))/2; + params.tau = (a_prev - 1)/a_cur; + a_prev = a_cur; + elseif (optimizeTau) % Compute optimized tau + + % tau = argmin ||u - Phi*y_{i+1}|| + % = <res, Phi*(x_cur - x_prev)>/||Phi*(x_cur - x_prev)||^2 + + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + if (i == 1) + Phi_diff = Phi_x_cur; + else + Phi_diff = Phi_x_cur - Phi_x_prev; + end; + params.tau = res'*Phi_diff/(Phi_diff'*Phi_diff); + end; + + y_cur = x_cur + params.tau*(x_cur - x_prev); + Y_i = find(ne(y_cur, 0)); + else + y_cur = x_cur + params.tau(i)*(x_cur - x_prev); + Y_i = find(ne(y_cur, 0)); + end; + + % Test stopping criterion + if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur)) + break; + end; + i = i + 1; +end; + +x_hat = x_cur; +numiter= i; + +if (i > params.ALPSiters) + x_path = x_path(:,1:numiter-1); +else + x_path = x_path(:,1:numiter); +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/thresh.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,4 @@ +function x= thresh(x, K) + +[~, ind] = sort(abs(x), 'descend'); +x(ind(K+1:end))= 0; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPS/zero_ALPS.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,189 @@ +function [x_hat, numiter, x_path] = zero_ALPS(y, Phi, K, params) +% ========================================================================= +% 0-ALPS(#) algorithm - Beta Version +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithm with memoryless acceleration. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% INPUT ARGUMENTS: +% y M x 1 undersampled measurement vector. +% Phi M x N regression matrix. +% K Sparsity of underlying vector x* or desired +% sparsity of solution. +% params Structure of parameters. These are: +% +% tol,... Early stopping tolerance. Default value: tol = +% 1-e5 +% ALPSiters,... Maximum number of algorithm iterations. Default +% value: 300. +% mode, ... According to [1], possible values are +% [0,1,2,4,5,6]. This value comes from the binary +% representation of the parameters: +% (solveNewtob, gradientDescentx, solveNewtonx), +% which are explained next. Default value = 0. +% solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a +% Newton system restricted to a sparse support. +% It is implemented via conjugate gradients. +% If solveNewtonb == 0: Step size selection as described +% in eqs. (12) and (13) in [1]. +% Default value: solveNewtonb = 0. +% gradientDescentx,... If gradientDescentx == 1: single gradient +% update of x_{i+1} restricted ot its support with +% line search. Default value: gradientDescentx = +% 1. +% solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit +% (c.f. Simon Foucart, "Hard Thresholding Pursuit," +% preprint, 2010). Default vale: solveNewtonx = 0 +% mu,... Variable that controls the step size selection. +% When mu = 0, step size is computed adaptively +% per iteration. Default value: mu = 0. +% cg_maxiter,... Maximum iterations for Conjugate-Gradients method. +% cg_tol Tolerance variable for Conjugate-Gradients method. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x_hat N x 1 recovered K-sparse vector. +% numiter Number of iterations executed. +% x_path Keeps a series of computed N x 1 K-sparse vectors +% until the end of the iterative process. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. +% Email: jrom@acm.caltech.edu +% ========================================================================= +% This work was supported in part by the European Commission under Grant +% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like +% to acknowledge Rice University for his Faculty Fellowship. +% ========================================================================= + +[tmpArg, N] = size(Phi); + +%% Initialize transpose of measurement matrix + +Phi_t = Phi'; + +%% Initialize to zero vector +x_cur = zeros(N,1); +X_i = []; + +x_path = zeros(N, params.ALPSiters); + +%% CG params +if (params.mode == 1 || params.mode == 4 || params.mode == 5 || params.mode == 6) + cg_verbose = 0; + cg_A = Phi_t*Phi; + cg_b = Phi_t*y; +end; + +%% Determine step size selection strategy +function_mu = 0; +adaptive_mu = 0; + +if (isa(params.mu,'float')) + function_mu = 0; + if (params.mu == 0) + adaptive_mu = 1; + else + adaptive_mu = 0; + end; +elseif (isa(params.mu,'function_handle')) + function_mu = 1; +end; + +%% Help variables +complementary_Xi = ones(N,1); + +i = 1; +%% 0-ALPS(#) +while (i <= params.ALPSiters) + x_path(:,i) = x_cur; + x_prev = x_cur; + + % Compute the residual + if (i == 1) + res = y; + else + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + end; + + % Compute the derivative + der = Phi_t*res; + + % Determine S_i set via eq. (11) + complementary_Xi(X_i) = 0; + [tmpArg, ind_der] = sort(abs(der).*complementary_Xi, 'descend'); + complementary_Xi(X_i) = 1; + S_i = [X_i; ind_der(1:K)]; + + ider = der(S_i); + if (params.solveNewtonb == 1) + % Compute least squares solution of the system A*y = (A*A)x using CG + if (params.useCG == 1) + [b, tmpArg, tmpArg] = cgsolve(cg_A(S_i, S_i), cg_b(S_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + b = cg_A(S_i,S_i)\cg_b(S_i); + end; + else + % Step size selection via eq. (12) and eq. (13) + if (adaptive_mu) + Pder = Phi(:,S_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + b = x_cur(S_i) + (mu_bar)*ider; + elseif (function_mu) + b = x_cur(S_i) + params.mu(i)*ider; + else + b = x_cur(S_i) + params.mu*ider; + end; + end; + + % Hard-threshold b and compute X_{i+1} + [tmpArg, ind_b] = sort(abs(b), 'descend'); + x_cur = zeros(N,1); + x_cur(S_i(ind_b(1:K))) = b(ind_b(1:K)); + X_i = S_i(ind_b(1:K)); + + if (params.gradientDescentx == 1) + % Calculate gradient of estimated vector x_cur + Phi_x_cur = Phi(:,X_i)*x_cur(X_i); + res = y - Phi_x_cur; + der = Phi_t*res; + ider = der(X_i); + + if (adaptive_mu) + Pder = Phi(:,X_i)*ider; + mu_bar = ider'*ider/(Pder'*Pder); + x_cur(X_i) = x_cur(X_i) + mu_bar*ider; + elseif (function_mu) + x_cur(X_i) = x_cur(X_i) + params.mu(i)*ider; + else x_cur(X_i) = x_cur(X_i) + params.mu*ider; + end; + elseif (params.solveNewtonx == 1) + % Similar to HTP + if (params.useCG == 1) + [v, tmpArg, tmpArg] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose); + else + v = cg_A(X_i,X_i)\cg_b(X_i); + end; + x_cur(X_i) = v; + end; + + % Test stopping criterion + if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur)) + break; + end; + i = i + 1; + +end; + +x_hat = x_cur; +numiter = i; + +if (i > params.ALPSiters) + x_path = x_path(:,1:numiter-1); +else + x_path = x_path(:,1:numiter); +end;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPSbenchmark_errornorm_vs_iter.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,228 @@ +% ========================================================================= +% Algebraic Pursuit algorithms - Demo +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods +% on sparse recovery of linear inverse systems. In particular, let +% Phi : M x N real matrix (M < N) +% x* : N x 1 K-sparse data vector +% n : N x 1 additive noise vector +% then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS +% solve the following minimization problem +% minimize ||y - Phi x||^2 subject to x is K-sparse vector. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= + +clear all +close all +clc; + +addpath ALPS/; + +% General parameters +N = 4000; % Size of sparse vector x +M = 800; % Number of measurements +K = ceil(M*(0.2)); % Sparsity of x* +sigma = 0.0000; % Additive noise variance +verbose = 1; % Print information + +% ALPS parameters +its = 600; % Maximum number of iterations +Ttol= 1e-6; % Convergence tolerance +tau = -1; % Momentum term +mu = 0; % Step size selection +MC = 10; % Number of Monte-Carlo iterations +useCG = 0; % Conjugate Gradient + +% Please refer to generate_matrix.m and generate_vector.m files +m_ensemble = 'Gaussian'; % Measurement matrix ensemble type +x_ensemble = 'Gaussian'; % Sparse vector ensemble type + +max_its = 0; + +for i = 1:MC + x = generate_vector(N, K, x_ensemble); + Phi = generate_matrix(M, N, m_ensemble); + noise = sigma*randn(M,1); + y = Phi*x + noise; + + % 0-ALPS(0) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 0, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(1, i) = numiter; + error(1,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(1, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 0-ALPS(2) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(2, i) = numiter; + error(2,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(2, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 0-ALPS(4) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(3, i) = numiter; + error(3,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(3, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 0-ALPS(5) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(4, i) = numiter; + error(4,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(4, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 1-ALPS(2) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(5, i) = numiter; + error(5,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(5, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 1-ALPS(4) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(6, i) = numiter; + error(6,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(6, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 1-ALPS(5) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(7, i) = numiter; + error(7,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(7,i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 2-ALPS(2) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(8, i) = numiter; + error(8,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(8, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 2-ALPS(4) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(9, i) = numiter; + error(9,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(9, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 2-ALPS(5) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(10, i) = numiter; + error(10,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(10, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; + + % 3-ALPS(2) + disp('====================================================================='); + [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 3, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); + num_iter(11, i) = numiter; + error(11,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + total_time(11, i) = time; + + if (numiter > max_its) + max_its = numiter; + end; +end; + +set(0,'DefaultAxesFontSize',12); + +figure; +semilogy(squeeze(median(error(1,:,:),2)),'LineWidth',3); hold on +semilogy(squeeze(median(error(2,:,:),2)),'--','LineWidth',3); +semilogy(squeeze(median(error(3,:,:),2)),':','LineWidth',3); +semilogy(squeeze(median(error(4,:,:),2)),'k','LineWidth',4); +semilogy(squeeze(median(error(5,:,:),2)),'k--','LineWidth',4); +semilogy(squeeze(median(error(6,:,:),2)),'k:','LineWidth',4); +semilogy(squeeze(median(error(7,:,:),2)),'k-.','LineWidth',4); +semilogy(squeeze(median(error(8,:,:),2)),'r--','LineWidth',3); +semilogy(squeeze(median(error(9,:,:),2)),'m--','LineWidth',3); +semilogy(squeeze(median(error(10,:,:),2)),'k--','LineWidth',3); +semilogy(squeeze(median(error(11,:,:),2)),'g--','LineWidth',3); + +xlabel('[i]:iteration number') +ylabel('||x-x_i||') +axis([1 max_its + 10 10^-5 2*10^0]) + +legend(strcat('0-IHT(0) [', num2str(mean(num_iter(1,:))),'][', num2str(mean(total_time(1,:))),']'), ... +strcat('0-IHT(2) [', num2str(mean(num_iter(2,:))),'][', num2str(mean(total_time(2,:))),']'), ... +strcat('0-IHT(4) [', num2str(mean(num_iter(3,:))),'][', num2str(mean(total_time(3,:))),']'), ... +strcat('0-IHT(5) [', num2str(mean(num_iter(4,:))),'][', num2str(mean(total_time(4,:))),']'), ... +strcat('1-IHT(2) [', num2str(mean(num_iter(5,:))),'][', num2str(mean(total_time(5,:))),']'), ... +strcat('1-IHT(4) [', num2str(mean(num_iter(6,:))),'][', num2str(mean(total_time(6,:))),']'), ... +strcat('1-IHT(5) [', num2str(mean(num_iter(7,:))),'][', num2str(mean(total_time(7,:))),']'), ... +strcat('2-ALPS(2) [', num2str(mean(num_iter(8,:))),'][', num2str(mean(total_time(8,:))),']'), ... +strcat('2-ALPS(4) [', num2str(mean(num_iter(9,:))),'][', num2str(mean(total_time(9,:))),']'), ... +strcat('2-ALPS(5) [', num2str(mean(num_iter(10,:))),'][', num2str(mean(total_time(10,:))),']'), ... +strcat('3-ALPS(2) [', num2str(mean(num_iter(11,:))),'][', num2str(mean(total_time(11,:))),']')); +title( strcat('N= ', num2str(N), ', M= ', num2str(M), ', K= ', num2str(K)) ); +shg;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPSdemo.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,79 @@ +% ========================================================================= +% Algebraic Pursuit algorithms - Demo +% ========================================================================= +% Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods +% on sparse recovery of linear inverse systems. In particular, let +% Phi : M x N real matrix (M < N) +% x* : N x 1 K-sparse data vector +% n : N x 1 additive noise vector +% then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS +% solve the following minimization problem +% minimize ||y - Phi x||^2 subject to x is K-sparse vector. +% +% Detailed discussion on the algorithm can be found in +% [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written +% by Volkan Cevher, Technical Report, 2011. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= + +clear all +close all +clc + +addpath ALPS/; + +% General parameters +N = 1000; % Size of sparse vector x +M = 300; % Number of measurements +K = ceil(M*(0.25)); % Sparsity of x* +sigma = 0.000; % Additive noise variance +verbose = 1; % Print information + +% ALPS parameters +its = 600; % Maximum number of iterations +Ttol= 1e-6; % Convergence tolerance +mode = 1; % Binary representation of (SolveNewtonb, GradientDescentx, SolveNewtox) +memory = 1; % Memory of IHT method +tau = 0; % Momentum term +useCG = 0; % Usage of Conjugate gradients method +mu = 0; % Step size selection + +% Please refer to generate_matrix.m and generate_vector.m files +m_ensemble = 'Gaussian'; % Measurement matrix ensemble type +x_ensemble = 'Gaussian'; % Sparse vector ensemble type + +% reset(RandStream.getDefaultStream); + +%% Generate data +x = generate_vector(N, K, x_ensemble); +Phi = generate_matrix(M, N, m_ensemble); +noise = sigma*randn(M,1); +y = Phi*x + noise; + +%% ALPS reconstruction +disp('====================================================================='); +[x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', memory, 'mode', mode, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); +if (numiter ~= -1) + str = sprintf('Error norm: %f ', norm(x-x_hat)); + disp(str); +end; + +%% Plot error per iteration +if (numiter ~= -1) + error = sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2)); + + set(0, 'DefaultAxesFontSize', 14); + figure; + semilogy(error,'LineWidth',3); hold on; + grid on; + semilogy((norm(noise))*ones(1,its),'m-','LineWidth',1); + + xlabel('[i]:iteration number') + ylabel('||x-x_i||') + axis([1 numiter + 5 10^-5 2*10^0]) + + legend(strcat('Time = [',num2str(time),']'), 'Noise Level'); + title(strcat('N = ', num2str(N), ', M = ', num2str(M), ', K = ', num2str(K))); + shg; +end; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/README Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,19 @@ +This archive includes: + +1. ALPS folder +2. generate_matrix.m: generates measurement matrix. +3. generate_vector.m: generates sparse vector. +4. Report_bibtex.tex: bibtex reference for citing. +5. ALPSdemo.m: demo of ALPS. +6. ALPSbenchmark_errornorm_vs_iter.m: generates figure of average error norm per iteration vs. # of iterations. + +ALPS folder includes: + +1. AlgebraicPursuit.m: "Wrapper" code. +2. zero_ALPS.m: Implementation of 0-memory acceleration. +3. one_ALPS.m: Implementation of 1-memory acceleration. +4. l_ALPS.m: Implementation of l-memory acceleration. +5. infty_ALPS.m: Implementation of infty-memory acceleration. +6. cgsolve.m: CG algorithm created by Justin Romberg, Caltech. + +Please run ALPSdemo in MATLAB command line to verify installation before use.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/Report_bibtex.bib Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,27 @@ +@comment{ generated by <http://infoscience.epfl.ch/> } + +@TechReport{EPFL-REPORT-162761, + abstract = {We propose and analyze acceleration schemes for hard + thresholding methods with applications to sparse + approximation in linear inverse systems. Our acceleration + schemes fuse combinatorial, sparse projection algorithms + with convex optimization algebra to provide + computationally efficient and robust sparse recovery + methods. We compare and contrast the (dis)advantages of + the proposed schemes with the state-of-the-art, not only + within hard thresholding methods, but also within convex + sparse recovery algorithms.}, + affiliation = {EPFL}, + author = {Cevher, Volkan}, + details = {http://infoscience.epfl.ch/record/162761}, + documenturl = {http://infoscience.epfl.ch/record/162761/files/techreport.pdf}, + keywords = {structured sparsity; sparse recovery; hard thresholding}, + oai-id = {oai:infoscience.epfl.ch:162761}, + oai-set = {report; fulltext; fulltext-public}, + status = {PUBLISHED}, + submitter = {199128; 199128}, + title = {On {A}ccelerated {H}ard {T}hresholding {M}ethods for + {S}parse {A}pproximation}, + unit = {LIONS}, + year = 2011 +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/generate_matrix.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,59 @@ +function [A] = generate_matrix(M, N, ensemble, p) +% ========================================================================= +% Measerement matrix generator +% ========================================================================= +% INPUT ARGUMENTS: +% M Number of measurements (number of rows). +% N Size of sparse vector (number of columns). +% ensemble Ensemble type of measurement matrix. Possible +% values are: +% -'Gaussian': creates a MxN measurement matrix with +% elements drawn from normal distribution N(0,1).% +% -'Bernoulli': creates a MxN measurement matrix with +% elements drawn from Bernoulli distribution (1/2,1/2). +% -'pBernoulli': creates a MxN measurement matrix with +% elements drawn from Bernoulli distribution (p,1-p). +% Parameter of Bernoulli distribution. +% -'sparseGaussian': creates a MxN sparse measurement +% matrix with elements drawn from normal distribution N(0,1). +% ========================================================================= +% OUTPUT ARGUMENTS: +% A MxN measurement matrix with normalized columns. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= +if nargin < 3 + ensemble = 'Gaussian'; +end; + +if nargin < 4 + p = 0.5; +end; + +switch ensemble + case 'Gaussian' + A = randn(M,N); % Standard normal distribution + for i = 1:N % Normalize columns + A(:,i) = A(:,i)/norm(A(:,i)); + end; + case 'Bernoulli' + A = (-1).^round(rand(M,N)); % Bernoulli ~ (1/2, 1/2) distribution + for i = 1:N % Normalize columns + A(:,i) = A(:,i)/norm(A(:,i)); + end; + case 'pBernoulli' + A = (-1).^(rand(M,N) > p); % Bernoulli ~ (p, 1-p) distribution + for i = 1:N % Normalize columns + A(:,i) = A(:,i)/norm(A(:,i)); + end; + case 'sparseGaussian' + leftd = 8; + A = zeros(M,N); + for i = 1:N + ind = randperm(M); + A(ind(1:leftd),i)=1/leftd; + end + for i = 1:N % Normalize columns + A(:,i) = A(:,i)/norm(A(:,i)); + end; +end; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/generate_vector.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,52 @@ +function [x] = generate_vector(N, K, ensemble, sigma, p) +% ========================================================================= +% Sparse vector generator +% ========================================================================= +% INPUT ARGUMENTS: +% N Size of sparse vector. +% K Sparsity of vector. +% ensemble Ensemble type of measurement matrix. Possible +% values are: +% -'Gaussian': K non-zero elements of the sparse vector +% are drawn from normal distribution N(0,1). +% -'sGaussian': K non-zero elements of the sparse vector +% are drawn from normal distribution N(0,sigma^2). +% -'Bernoulli': K non-zero elements of the sparse vector +% are drawn from Bernoulli distribution (1/2,1/2). +% -'pBernoulli': K non-zero elements of the sparse vector +% are drawn from Bernoulli distribution (p,1-p). +% sigma Standard deviation of Gaussian distribution. +% p Parameter of Bernoulli distribution. +% ========================================================================= +% OUTPUT ARGUMENTS: +% x K-sparse vector. +% ========================================================================= +% 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. +% ========================================================================= + +if nargin < 3 + ensemble = 'Gaussian'; +end; + +if nargin < 4 + sigma = 1; + p = 0.5; +end; + +x = zeros(N,1); +rand_indices = randperm(N); + +switch ensemble + case 'Gaussian' + x(rand_indices(1:K)) = randn(K,1); % Standard normal distribution ~ N(0,1) + case 'sGaussian' + x(rand_indices(1:K)) = sigma*randn(K,1); % Normal distribution ~ N(0,sigma^2) + case 'Uniform' + x(rand_indices(1:K)) = rand(K,1); % Uniform [0,1] distribution + case 'Bernoulli' + x(rand_indices(1:K)) = (-1).^round(rand(K,1)); % Bernoulli ~ (1/2, 1/2) distribution + case 'pBernoulli' + x(rand_indices(1:K)) = (-1).^(rand(K,1) > p); % Bernoulli ~ (p, 1-p) distribution +end; + +x = x/norm(x); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/wrapper_ALPS_toolbox.m Wed Sep 07 14:17:30 2011 +0100 @@ -0,0 +1,80 @@ +function [x , numiter, time, x_path] = wrapper_ALPS_toolbox(b, A, param) +%% SMALL wrapper for ALPS toolbox +% +% Function gets as input +% b - measurement vector +% A - dictionary +% K - desired sparsity level +% param - structure containing additional parameters +% Output: +% x - sparse solution +% numiter - number of iterations +% time - time needed to solve the problem# +% x_path - matrix containing x after every iteration + +% 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 isfield(param, 'sparsity') + sparsity = param.sparsity; +else + printf('\nAlebraic Pursuit algorithms require desired sparsity level.\n("sparsity" field in solver parameters structure)\n '); + return +end + +if isfield(param, 'memory') + memory = param.memory; +else + memory = 0; +end + +if isfield(param, 'mode') + mode = param.mode; +else + mode = 0; +end +if isfield(param, 'tolerance') + tolerance = param.tolerance; +else + tolerance = 1e-5; +end +if isfield(param, 'iternum') + iternum = param.iternum; +else + iternum = 300; +end +if isfield(param, 'verbose') + verbose = param.verbose; +else + verbose = 0; +end +if isfield(param, 'tau') + tau = param.tau; +else + tau = 1/2; +end +if isfield(param, 'useCG') + useCG = param.useCG; +else + useCG = 0; +end +if isfield(param, 'mu') + mu = param.mu; +else + mu = 0; +end +training_size = size(b,2); +x=zeros(size(A,2),training_size); +for i = 1:training_size + [x(:,i), numiter, time, x_path] = AlgebraicPursuit(b(:,i), A, sparsity, 'memory', memory,... + 'mode', mode, 'tolerance', tolerance, 'ALPSiterations', iternum, ... + 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu); +end \ No newline at end of file
--- a/util/SMALL_AudioDeNoiseResult.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_AudioDeNoiseResult.m Wed Sep 07 14:17:30 2011 +0100 @@ -2,7 +2,7 @@ %% Plots the results of Audio denoising experiment - underconstruction % Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2009 Ivan Damnjanovic. +% 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 @@ -13,7 +13,7 @@ fMain=figure('Name', sprintf('File %s (training set size- %d, sigma - %d)',SMALL.Problem.name, SMALL.Problem.n, SMALL.Problem.sigma)); m=size(SMALL.solver,2); -maxval=SMALL.Problem.maxval; +maxval=max(SMALL.Problem.Original); au=SMALL.Problem.Original; aunoise=SMALL.Problem.Noisy; @@ -25,7 +25,7 @@ for i=1:m params=SMALL.solver(i).param; - sWav=subplot(2, m, m+i, 'Parent', fMain); plot(SMALL.solver(i).reconstructed.Image/maxval, 'Parent', sWav); + sWav=subplot(2, m, m+i, 'Parent', fMain); plot(SMALL.solver(i).reconstructed.audio/maxval, 'Parent', sWav); title(sprintf('%s Denoised audio, PSNR: %.2fdB', SMALL.DL(i).name, SMALL.solver(i).reconstructed.psnr),'Parent', sWav ); if strcmpi(SMALL.DL(i).name,'ksvds') D = kron(SMALL.Problem.basedict{2},SMALL.Problem.basedict{1})*SMALL.DL(i).D;
--- a/util/SMALL_init_DL.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_init_DL.m Wed Sep 07 14:17:30 2011 +0100 @@ -1,4 +1,4 @@ -function DL = SMALL_init_DL(varargin) +function DL = SMALL_init_DL(toolbox, name, param, profile) %% Function initialise SMALL structure for Dictionary Learning. % Optional input variables: % toolbox - name of Dictionary Learning toolbox you want to use @@ -17,9 +17,27 @@ % %% -DL.toolbox=[]; -DL.name=[]; -DL.param=[]; +if ~ exist( 'toolbox', 'var' ) || isempty(toolbox) + DL.toolbox = []; +else + DL.toolbox = toolbox; +end +if ~ exist( 'name', 'var' ) || isempty(name) + DL.name = []; +else + DL.name = name; +end +if ~ exist( 'param', 'var' ) || isempty(param) + DL.param = []; +else + DL.param = param; +end +if ~ exist( 'profile', 'var' ) || isempty(profile) + DL.profile = 1; +else + DL.profile = profile; +end + DL.D=[]; DL.time=[]; end \ No newline at end of file
--- a/util/SMALL_learn.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_learn.m Wed Sep 07 14:17:30 2011 +0100 @@ -18,8 +18,9 @@ % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. %% - +if (DL.profile) fprintf('\nStarting Dictionary Learning %s... \n', DL.name); +end start=cputime; tStart=tic; if strcmpi(DL.toolbox,'KSVD') @@ -58,6 +59,30 @@ D(:,i)=D(:,i)/norm(D(:,i)); end + elseif strcmpi(DL.toolbox,'TwoStepDL') + + DL=SMALL_two_step_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; + +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> @@ -83,9 +108,11 @@ %% % Dictionary Learning time tElapsed=toc(tStart); - DL.time = cputime - start; +DL.time = cputime - start; +if (DL.profile) fprintf('\n%s finished task in %2f seconds (cpu time). \n', DL.name, DL.time); fprintf('\n%s finished task in %2f seconds (tic-toc time). \n', DL.name, tElapsed); +end DL.time=tElapsed; % If dictionary is given as a sparse matrix change it to full
--- a/util/SMALL_midiGenerate.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_midiGenerate.m Wed Sep 07 14:17:30 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_plot.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_plot.m Wed Sep 07 14:17:30 2011 +0100 @@ -21,7 +21,7 @@ for i =1:m subplot(m,n, (i-1)*n+1); plot(1:length(SMALL.solver(i).solution), SMALL.solver(i).solution, 'b') - title(sprintf('%s(%s) in %.2f s', SMALL.solver(i).name, SMALL.solver(i).param,SMALL.solver(i).time),'Interpreter','none') + title(sprintf('%s in %.2f s', SMALL.solver(i).name, SMALL.solver(i).time),'Interpreter','none') xlabel('Coefficient') % Plot reconstructed signal against original @@ -31,7 +31,7 @@ subplot(m,n,(i-1)*n+j); plot(1:length(SMALL.solver(i).reconstructed(:,j-1)), SMALL.solver(i).reconstructed(:,j-1) ,'b.-', 1:length(SMALL.Problem.signal(:,j-1)), SMALL.Problem.signal(:,j-1),'r--') %legend(SMALL.solver(i).name,'Original signal',0,'Location','SouthOutside','Interpreter','none'); - title(sprintf('%s(%s) in %.2f s signal %d', SMALL.solver(i).name, SMALL.solver(i).param,SMALL.solver(i).time,n-1),'Interpreter','none') + title(sprintf('%s in %.2f s signal %d', SMALL.solver(i).name, SMALL.solver(i).time,n-1),'Interpreter','none') end end end
--- a/util/SMALL_solve.m Tue Jul 26 16:02:59 2011 +0100 +++ b/util/SMALL_solve.m Wed Sep 07 14:17:30 2011 +0100 @@ -1,5 +1,5 @@ function solver = SMALL_solve(Problem, solver) -%% SMALL sparse solver +%% SMALL sparse solver caller function % % Function gets as input SMALL structure that contains SPARCO problem to % be solved, name of the toolbox and solver, and parameters file for @@ -48,7 +48,11 @@ if strcmpi(solver.toolbox,'sparselab') y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']); elseif strcmpi(solver.toolbox,'sparsify') - y = eval([solver.name,'(b, A, n, ''P_trans'', AT,',solver.param,');']); + if isa(Problem.A,'float') + y = eval([solver.name,'(b, A, n,',solver.param,');']); + else + y = eval([solver.name,'(b, A, n, ''P_trans'', AT,',solver.param,');']); + end elseif (strcmpi(solver.toolbox,'spgl1')||strcmpi(solver.toolbox,'gpsr')) y = eval([solver.name,'(b, A,',solver.param,');']); elseif (strcmpi(solver.toolbox,'SPAMS')) @@ -66,27 +70,43 @@ y = eval([solver.name,'(A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']); elseif (strcmpi(solver.toolbox, 'ompsbox')) basedict = Problem.basedict; - if issparse(Problem.A) + if issparse(Problem.A) A = Problem.A; - else + else A = sparse(Problem.A); - end + end G = dicttsep(basedict,A,dictsep(basedict,A,speye(size(A,2)))); epsilon=solver.param.epsilon; maxatoms=solver.param.maxatoms; y = eval([solver.name,'(basedict, A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']); Problem.sparse=1; -% To introduce new sparse representation algorithm put the files in -% your Matlab path. Next, unique name <TolboxID> for your toolbox and -% prefferd API <Preffered_API> needs to be defined. -% -% elseif strcmpi(solver.toolbox,'<ToolboxID>') -% -% % - Evaluate the function (solver.name - defined in the main) with -% % parameters given above -% -% y = eval([solver.name,'(<Preffered_API>);']); - +elseif (strcmpi(solver.toolbox, 'ALPS')) + if ~isa(Problem.A,'float') + % ALPS does not accept implicit dictionary definition + 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') + % MMbox 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 + % prefferd API <Preffered_API> needs to be defined. + % + % elseif strcmpi(solver.toolbox,'<ToolboxID>') + % + % % - Evaluate the function (solver.name - defined in the main) with + % % parameters given above + % + % y = eval([solver.name,'(<Preffered_API>);']); + else printf('\nToolbox has not been registered. Please change SMALL_solver file.\n'); return