Mercurial > hg > smallbox
changeset 159:23763c5fbda5 danieleb
Merge
| author | Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk> |
|---|---|
| date | Wed, 31 Aug 2011 10:43:32 +0100 |
| parents | 855aa3288394 (current diff) b14209313ba4 (diff) |
| children | e3035d45d014 |
| files | util/SMALL_solve.m |
| diffstat | 50 files changed, 3516 insertions(+), 60 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 2011 +0100 @@ -0,0 +1,74 @@ +function [X , cost] = wrapper_mm_solver(b, A, param) +%% SMALL wrapper for Majorization Maximization toolbos solver +% +% Function gets as input +% b - measurement vector +% A - dictionary +% param - structure containing additional parameters +% Output: +% x - sparse solution +% cost - Objective cost + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + +% Initial guess for the coefficients + +if (isfield(param, 'initcoeff')) + initX = param.initcoeff; +else + initX = zeros(size(A,2),size(b,2)); +end + +% to - 1/(step size) . It is larger than spectral norm of dictionary A + +if isfield(param, 'to') + to = param.to; +else + to = .1+svds(A,1); +end + +% lambda - Lagrangian multiplier. (regulates shrinkage) + +if isfield(param, 'lambda') + lambda = param.lambda; +else + lambda = 2*.2; +end + +% Inner-loop maximum iteration number. + +if isfield(param, 'iternum') + maxIT = param.iternum; +else + maxIT = 1000; +end + +% Stopping criterion for iterative softthresholding + +if isfield(param, 'epsilon') + epsilon = param.epsilon; +else + epsilon = 1e-7; +end + +% Debiasing. 0 = No, 1 = Yes + +if isfield(param, 'map') + map = param.map; +else + map = 1; +end + + +[X, cost] = mm1(A,b,initX,to,lambda,maxIT,epsilon,map); +cost +end \ No newline at end of file
--- a/Problems/generateAudioDeclippingProblem.m Wed Aug 31 10:37:57 2011 +0100 +++ b/Problems/generateAudioDeclippingProblem.m Wed Aug 31 10:43:32 2011 +0100 @@ -103,7 +103,7 @@ data.fs = x.fs; data.nbits = x.nbits; - +data.Upper_Limit = max(solutiondata.XClean); [data.m, data.n] = size(x_clip); data.p = windowSize*redundancyFactor; %number of dictionary elements
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/ALPS solvers tests/SMALL_solver_test_ALPS.m Wed Aug 31 10:43:32 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_cgp'; + +% 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 Wed Aug 31 10:37:57 2011 +0100 +++ b/examples/AudioInpainting/Audio_Declipping_Example.m Wed Aug 31 10:43:32 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_AMT_DL_test_KSVD_MM.m Wed Aug 31 10:43:32 2011 +0100 @@ -0,0 +1,270 @@ +%% Dictionary Learning for Automatic Music Transcription - KSVD vs SPAMS +% +% +% This file contains an example of how SMALLbox can be used to test diferent +% dictionary learning techniques in Automatic Music Transcription problem. +% It calls generateAMT_Learning_Problem that will let you to choose midi, +% wave or mat file to be transcribe. If file is midi it will be first +% converted to wave and original midi file will be used for comparison with +% results of dictionary learning and reconstruction. +% The function will generarte the Problem structure that is used to learn +% Problem.p notes spectrograms from training set Problem.b using +% dictionary learning technique defined in DL structure. +% Two dictionary learning techniques were compared: +% +% - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient +% Implementation of the K-SVD Algorithm using Batch Orthogonal +% Matching Pursuit", Technical Report - CS, Technion, April 2008. +% +% - MMDL - M. Yaghoobi, T. Blumensath and M. Davies, "Dictionary Learning +% for Sparse Approximations with the Majorization Method", IEEE +% Trans. on Signal Processing, Vol. 57, No. 6, pp 2178-2191, +% 2009. + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +%% + +clear; + + +% Defining Automatic Transcription of Piano tune as Dictionary Learning +% Problem + +SMALL.Problem = generateAMT_Learning_Problem('',2048,0.75); + +%% +% Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in +% SMALL.Problem.p) using sparsity constrain only + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL(); + +% Defining fields needed for dictionary learning + +SMALL.DL(1).toolbox = 'KSVD'; +SMALL.DL(1).name = 'ksvd'; +% Defining the parameters for KSVD +% In this example we are learning 88 atoms in 100 iterations, so that +% every frame in the training set can be represented with maximum Tdata +% dictionary elements. Type help ksvd in MATLAB prompt for more options. + +SMALL.DL(1).param=struct(... + 'Tdata', 5,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 50); + +% Learn the dictionary + +SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); + +% Set SMALL.Problem.A dictionary and reconstruction function +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(1).D; +SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values +% As an example, SPAMS (Julien Mairal 2009) implementation of LARS +% algorithm is used for representation of training set in the learned +% dictionary. + +SMALL.solver(1)=SMALL_init_solver; + +% Defining the parameters needed for sparse representation + +SMALL.solver(1).toolbox='SMALL'; +SMALL.solver(1).name='SMALL_cgp'; + +% Here we use mexLasso mode=2, with lambda=2, lambda2=0 and positivity +% constrain (type 'help mexLasso' for more information about modes): +% +% min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2 + +SMALL.solver(1).param='20, 1e-2'; +% struct(... +% 'lambda', 2,... +% 'pos', 1,... +% 'mode', 2); + +% Call SMALL_soolve to represent the signal in the given dictionary. +% As a final command SMALL_solve will call above defined reconstruction +% function to reconstruct the training set (Problem.b) in the learned +% dictionary (Problem.A) + +SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +%% +% Analysis of the result of automatic music transcription. If groundtruth +% exists, we can compare transcribed notes and original and get usual +% True Positives, False Positives and False Negatives measures. + +if ~isempty(SMALL.Problem.notesOriginal) + AMT_res(1) = AMT_analysis(SMALL.Problem, SMALL.solver(1)); +end + + + +%% +% % Here we solve the same problem using non-negative sparse coding with +% % SPAMS online dictionary learning (Julien Mairal 2009) +% % +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values +% As an example, SPAMS (Julien Mairal 2009) implementation of LARS +% algorithm is used for representation of training set in the learned +% dictionary. + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for sparse representation + +SMALL.solver(2).toolbox='SPAMS'; +SMALL.solver(2).name='mexLasso'; + +% Here we use mexLasso mode=2, with lambda=3, lambda2=0 and positivity +% constrain (type 'help mexLasso' for more information about modes): +% +% min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2 + +SMALL.solver(2).param=struct('lambda', 3, 'pos', 1, 'mode', 2); + + +% You can also test ALPS, IST from MMbox or any other solver, but results +% are not as good as SPAMS +% +% % Initialising solver structure +% % Setting solver structure fields (toolbox, name, param, solution, +% % reconstructed and time) to zero values +% +% SMALL.solver(2)=SMALL_init_solver; +% +% % Defining the parameters needed for image denoising +% +% SMALL.solver(2).toolbox='ALPS'; +% SMALL.solver(2).name='AlebraicPursuit'; +% +% SMALL.solver(2).param=struct(... +% 'sparsity', 10,... +% 'memory', 1,... +% 'mode', 6,... +% 'iternum', 100,... +% 'tau',-1,... +% 'tolerance', 1e-14',... +% 'verbose',1); + +% % Initialising Dictionary structure +% % Setting Dictionary structure fields (toolbox, name, param, D and time) +% % to zero values +% % Initialising solver structure +% % Setting solver structure fields (toolbox, name, param, solution, +% % reconstructed and time) to zero values +% +% SMALL.solver(2)=SMALL_init_solver; +% +% % Defining the parameters needed for image denoising +% +% SMALL.solver(2).toolbox='MMbox'; +% SMALL.solver(2).name='mm1'; +% SMALL.solver(2).param=struct(... +% 'lambda',50,... +% 'iternum',1000,... +% 'map',0); + +SMALL.DL(2)=SMALL_init_DL('MMbox', 'MM_cn', '', 1); + + +% Defining the parameters for Majorization Minimization dictionary update +% +% In this example we are learning 88 atoms in 200 iterations, so that + + +SMALL.DL(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.A,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 200,... + 'iterDictUpdate', 1000,... + 'epsDictUpdate', 1e-7,... + 'cvset',0,... + 'show_dict', 0); + + +% Learn the dictionary + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Set SMALL.Problem.A dictionary and reconstruction function +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(2).D; +SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem); + + +% Call SMALL_soolve to represent the signal in the given dictionary. +% As a final command SMALL_solve will call above defined reconstruction +% function to reconstruct the training set (Problem.b) in the learned +% dictionary (Problem.A) + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + + +% Analysis of the result of automatic music transcription. If groundtruth +% exists, we can compare transcribed notes and original and get usual +% True Positives, False Positives and False Negatives measures. + +if ~isempty(SMALL.Problem.notesOriginal) + AMT_res(2) = AMT_analysis(SMALL.Problem, SMALL.solver(2)); +end + + +% Plot results and save midi files + +if ~isempty(SMALL.Problem.notesOriginal) + figAMT = SMALL_AMT_plot(SMALL, AMT_res); +else + figAMT = figure('Name', 'Automatic Music Transcription KSVD vs SPAMS'); + subplot(2,1,1); plot(SMALL.solver(1).reconstructed.notes(:,5), SMALL.solver(1).reconstructed.notes(:,3), 'kd '); + title (sprintf('%s dictionary in %.2f s', SMALL.DL(1).name, SMALL.DL(1).time)); + xlabel('Time'); + ylabel('Note Number'); + subplot(2,1,2); plot(SMALL.solver(2).reconstructed.notes(:,5), SMALL.solver(2).reconstructed.notes(:,3), 'b* '); + title (sprintf('%s dictionary in %.2f s', SMALL.DL(2).name, SMALL.DL(2).time)); + xlabel('Time'); + ylabel('Note Number'); +end + +FS=filesep; + +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'results']); + +[filename,pathname] = uiputfile({' *.mid;' },'Save KSVD result midi'); +if filename~=0 writemidi(SMALL.solver(1).reconstructed.midi, [pathname,FS,filename]);end + +[filename,pathname] = uiputfile({' *.mid;' },'Save SPAMS result midi'); +if filename~=0 writemidi(SMALL.solver(2).reconstructed.midi, [pathname,FS,filename]);end + +[filename,pathname] = uiputfile({' *.fig;' },'Save KSVD vs SPAMS AMT figure'); +if filename~=0 saveas(figAMT, [pathname,FS,filename]);end + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/MajorizationMinimization tests/SMALL_ImgDenoise_DL_test_KSVDvsMajorizationMinimization.m Wed Aug 31 10:43:32 2011 +0100 @@ -0,0 +1,215 @@ +%% Dictionary Learning for Image Denoising - KSVD vs Recursive Least Squares +% +% This file contains an example of how SMALLbox can be used to test different +% dictionary learning techniques in Image Denoising problem. +% It calls generateImageDenoiseProblem that will let you to choose image, +% add noise and use noisy image to generate training set for dictionary +% learning. +% Two dictionary learning techniques were compared: +% +% - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient +% Implementation of the K-SVD Algorithm using Batch Orthogonal +% Matching Pursuit", Technical Report - CS, Technion, April 2008. +% +% - MMDL - M. Yaghoobi, T. Blumensath and M. Davies, "Dictionary Learning +% for Sparse Approximations with the Majorization Method", IEEE +% Trans. on Signal Processing, Vol. 57, No. 6, pp 2178-2191, 2009. + + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + + + +% If you want to load the image outside of generateImageDenoiseProblem +% function uncomment following lines. This can be useful if you want to +% denoise more then one image for example. +% Here we are loading test_image.mat that contains structure with 5 images : lena, +% barbara,boat, house and peppers. +clear; +TMPpath=pwd; +FS=filesep; +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'images']); +load('test_image.mat'); +cd(TMPpath); + +% Deffining the noise levels that we want to test + +noise_level=[10 20 25 50 100]; + +% Here we loop through different noise levels and images + +for noise_ind=2:2 +for im_num=1:1 + +% Defining Image Denoising Problem as Dictionary Learning +% Problem. As an input we set the number of training patches. + +SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); +SMALL.Problem.name=int2str(im_num); + +Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain; +maxatoms = floor(prod(SMALL.Problem.blocksize)/2); + +% results structure is to store all results + +results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr; + +%% +% Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(1)=SMALL_init_DL(); + +% Defining the parameters needed for dictionary learning + +SMALL.DL(1).toolbox = 'KSVD'; +SMALL.DL(1).name = 'ksvd'; + +% Defining the parameters for KSVD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (Edata) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(1).param=struct(... + 'Edata', Edata,... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'exact', 1, ... + 'iternum', 20,... + 'memusage', 'high'); + +% Learn the dictionary + +SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(1).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +%% +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(1)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(1).toolbox='ompbox'; +SMALL.solver(1).name='omp2'; +SMALL.solver(1).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); + +% Show PSNR after reconstruction + +SMALL.solver(1).reconstructed.psnr + +%% +% For comparison purposes we will denoise image with Majorization +% Minimization method +% + +% Initialising solver structure +% Setting solver structure fields (toolbox, name, param, solution, +% reconstructed and time) to zero values + +SMALL.solver(2)=SMALL_init_solver; + +% Defining the parameters needed for image denoising + +SMALL.solver(2).toolbox='ompbox'; +SMALL.solver(2).name='omp2'; +SMALL.solver(2).param=struct(... + 'epsilon',Edata,... + 'maxatoms', maxatoms); + +% Initialising Dictionary structure +% Setting Dictionary structure fields (toolbox, name, param, D and time) +% to zero values + +SMALL.DL(2)=SMALL_init_DL('MMbox', 'MM_cn', '', 1); + + +% Defining the parameters for MOD +% In this example we are learning 256 atoms in 20 iterations, so that +% every patch in the training set can be represented with target error in +% L2-norm (EData) +% Type help ksvd in MATLAB prompt for more options. + + +SMALL.DL(2).param=struct(... + 'solver', SMALL.solver(2),... + 'initdict', SMALL.Problem.initdict,... + 'dictsize', SMALL.Problem.p,... + 'iternum', 20,... + 'iterDictUpdate', 1000,... + 'epsDictUpdate', 1e-7,... + 'cvset',0,... + 'show_dict', 0); + +% Learn the dictionary + +SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); + +% Set SMALL.Problem.A dictionary +% (backward compatiblity with SPARCO: solver structure communicate +% only with Problem structure, ie no direct communication between DL and +% solver structures) + +SMALL.Problem.A = SMALL.DL(2).D; +SMALL.Problem.reconstruct = @(x) ImgDenoise_reconstruct(x, SMALL.Problem); + +% Denoising the image - find the sparse solution in the learned +% dictionary for all patches in the image and the end it uses +% reconstruction function to reconstruct the patches and put them into a +% denoised image + +SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); + + + +% show results % + +SMALL_ImgDeNoiseResult(SMALL); + +results(noise_ind,im_num).psnr.ksvd=SMALL.solver(1).reconstructed.psnr; +results(noise_ind,im_num).psnr.odct=SMALL.solver(2).reconstructed.psnr; +results(noise_ind,im_num).vmrse.ksvd=SMALL.solver(1).reconstructed.vmrse; +results(noise_ind,im_num).vmrse.odct=SMALL.solver(2).reconstructed.vmrse; +results(noise_ind,im_num).ssim.ksvd=SMALL.solver(1).reconstructed.ssim; +results(noise_ind,im_num).ssim.odct=SMALL.solver(2).reconstructed.ssim; + + +results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; + +%clear SMALL; +end +end +% save results.mat results
--- a/solvers/CVX_add_const_Audio_declipping.m Wed Aug 31 10:37:57 2011 +0100 +++ b/solvers/CVX_add_const_Audio_declipping.m Wed Aug 31 10:43:32 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/toolboxes/alps/ALPS/AlgebraicPursuit.m Wed Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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 Aug 31 10:43:32 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_learn.m Wed Aug 31 10:37:57 2011 +0100 +++ b/util/SMALL_learn.m Wed Aug 31 10:43:32 2011 +0100 @@ -69,7 +69,20 @@ for i = 1: size(DL.D,2) DL.D(:,i)=DL.D(:,i)/norm(DL.D(:,i)); end + D = DL.D; + +elseif strcmpi(DL.toolbox,'MMbox') + + DL = wrapper_mm_DL(Problem, DL); + + % we need to make sure that columns are normalised to + % unit lenght. + + for i = 1: size(DL.D,2) + DL.D(:,i)=DL.D(:,i)/norm(DL.D(:,i)); + end D = DL.D; + % To introduce new dictionary learning technique put the files in % your Matlab path. Next, unique name <TolboxID> for your toolbox needs % to be defined and also prefferd API for toolbox functions <Preffered_API>
--- a/util/SMALL_midiGenerate.m Wed Aug 31 10:37:57 2011 +0100 +++ b/util/SMALL_midiGenerate.m Wed Aug 31 10:43:32 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 Wed Aug 31 10:37:57 2011 +0100 +++ b/util/SMALL_plot.m Wed Aug 31 10:43:32 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 Wed Aug 31 10:37:57 2011 +0100 +++ b/util/SMALL_solve.m Wed Aug 31 10:43:32 2011 +0100 @@ -15,7 +15,7 @@ % 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 isa(Problem.A,'float') @@ -45,58 +45,71 @@ start=cputime; tStart=tic; -switch solver.toolbox - case 'sparselab' - y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']); - case 'sparsify' +if strcmpi(solver.toolbox,'sparselab') + y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']); +elseif strcmpi(solver.toolbox,'sparsify') + 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,');']); - case 'spgl1' - y = eval([solver.name,'(b, A,',solver.param,');']); - case 'gpsr' - y = eval([solver.name,'(b, A,',solver.param,');']); - case 'SPAMS' - y = eval([solver.name,'(b, A, solver.param);']); - case 'SMALL' - if isa(Problem.A,'float') - y = eval([solver.name,'(A, b, n,',solver.param,');']); - else - y = eval([solver.name,'(A, b, n,',solver.param,',AT);']); - end - case 'ompbox' - G = A'*A; - maxatoms=solver.param.maxatoms; - switch solver.name - case 'omp' - y = omp(A,b,G,maxatoms,'checkdict','off'); - case 'omp2' - epsilon=solver.param.epsilon; - y = omp2(A,b,G,epsilon,'maxatoms',maxatoms,'checkdict','off'); - end - case 'ompsbox' - basedict = Problem.basedict; - if issparse(Problem.A) - A = Problem.A; - else - A = sparse(Problem.A); - 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>);']); - otherwise - printf('\nToolbox has not been registered. Please change SMALL_solver file.\n'); - return + end +elseif (strcmpi(solver.toolbox,'spgl1')||strcmpi(solver.toolbox,'gpsr')) + y = eval([solver.name,'(b, A,',solver.param,');']); +elseif (strcmpi(solver.toolbox,'SPAMS')) + y = eval([solver.name,'(b, A, solver.param);']); +elseif (strcmpi(solver.toolbox,'SMALL')) + if isa(Problem.A,'float') + y = eval([solver.name,'(A, b, n,',solver.param,');']); + else + y = eval([solver.name,'(A, b, n,',solver.param,',AT);']); + end +elseif (strcmpi(solver.toolbox, 'ompbox')) + G=A'*A; + epsilon=solver.param.epsilon; + maxatoms=solver.param.maxatoms; + y = eval([solver.name,'(A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']); +elseif (strcmpi(solver.toolbox, 'ompsbox')) + basedict = Problem.basedict; + if issparse(Problem.A) + A = Problem.A; + else + A = sparse(Problem.A); + 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; +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') + % ALPS does not accept implicit dictionary definition + A = opToMatrix(Problem.A, 1); + end + + [y, cost] = wrapper_mm_solver(b, A, solver.param); + + + + % To introduce new sparse representation algorithm put the files in + % your Matlab path. Next, unique name <TolboxID> for your toolbox and + % 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 end %%