Mercurial > hg > smallbox
diff toolboxes/alps/ALPSbenchmark_errornorm_vs_iter.m @ 154:0de08f68256b ivand_dev
ALPS toolbox - Algebraic Pursuit added to smallbox
| author | Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk> |
|---|---|
| date | Fri, 12 Aug 2011 11:17:47 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/alps/ALPSbenchmark_errornorm_vs_iter.m Fri Aug 12 11:17:47 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;