ivan@154: % =========================================================================
ivan@154: %                   Algebraic Pursuit algorithms - Demo
ivan@154: % =========================================================================
ivan@154: % Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods 
ivan@154: % on sparse recovery of linear inverse systems. In particular, let
ivan@154: %             Phi : M x N real matrix (M < N)
ivan@154: %             x*  : N x 1 K-sparse data vector
ivan@154: %             n   : N x 1 additive noise vector
ivan@154: % then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS 
ivan@154: % solve the following minimization problem 
ivan@154: %          minimize ||y - Phi x||^2    subject to x is K-sparse vector. 
ivan@154: % 
ivan@154: % Detailed discussion on the algorithm can be found in 
ivan@154: % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
ivan@154: % by Volkan Cevher, Technical Report, 2011.
ivan@154: % =========================================================================
ivan@154: % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
ivan@154: % =========================================================================
ivan@154: 
ivan@154: clear all
ivan@154: close all
ivan@154: clc;
ivan@154: 
ivan@154: addpath ALPS/;
ivan@154: 
ivan@154: % General parameters
ivan@154: N = 4000;                       % Size of sparse vector x
ivan@154: M = 800;                        % Number of measurements
ivan@154: K = ceil(M*(0.2));              % Sparsity of  x*
ivan@154: sigma = 0.0000;                 % Additive noise variance
ivan@154: verbose = 1;                    % Print information
ivan@154: 
ivan@154: % ALPS parameters
ivan@154: its = 600;                      % Maximum number of iterations
ivan@154: Ttol= 1e-6;                     % Convergence tolerance
ivan@154: tau = -1;                       % Momentum term
ivan@154: mu = 0;                         % Step size selection
ivan@154: MC = 10;                        % Number of Monte-Carlo iterations
ivan@154: useCG = 0;                      % Conjugate Gradient
ivan@154: 
ivan@154: % Please refer to generate_matrix.m and generate_vector.m files
ivan@154: m_ensemble = 'Gaussian';        % Measurement matrix ensemble type
ivan@154: x_ensemble = 'Gaussian';        % Sparse vector ensemble type
ivan@154: 
ivan@154: max_its = 0;
ivan@154: 
ivan@154: for i = 1:MC
ivan@154:     x = generate_vector(N, K, x_ensemble);
ivan@154:     Phi = generate_matrix(M, N, m_ensemble);
ivan@154:     noise = sigma*randn(M,1);
ivan@154:     y = Phi*x + noise;
ivan@154:     
ivan@154:     % 0-ALPS(0)    
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(1, i) = numiter;    
ivan@154:     error(1,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(1, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;    
ivan@154:     
ivan@154:     % 0-ALPS(2)    
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(2, i) = numiter;    
ivan@154:     error(2,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(2, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 0-ALPS(4)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(3, i) = numiter;    
ivan@154:     error(3,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(3, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 0-ALPS(5)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(4, i) = numiter;    
ivan@154:     error(4,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(4, i) = time;    
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 1-ALPS(2)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(5, i) = numiter;    
ivan@154:     error(5,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(5, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 1-ALPS(4)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(6, i) = numiter;    
ivan@154:     error(6,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(6, i) = time;    
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 1-ALPS(5)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(7, i) = numiter;    
ivan@154:     error(7,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(7,i) = time;    
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 2-ALPS(2)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(8, i) = numiter;    
ivan@154:     error(8,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(8, i) = time;
ivan@154:         
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;        
ivan@154:       
ivan@154:     % 2-ALPS(4)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(9, i) = numiter;    
ivan@154:     error(9,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(9, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 2-ALPS(5)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(10, i) = numiter;    
ivan@154:     error(10,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(10, i) = time;
ivan@154:     
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;
ivan@154:     
ivan@154:     % 3-ALPS(2)
ivan@154:     disp('=====================================================================');
ivan@154:     [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);
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154:     num_iter(11, i) = numiter;    
ivan@154:     error(11,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154:     total_time(11, i) = time;
ivan@154:         
ivan@154:     if (numiter > max_its)
ivan@154:         max_its = numiter;
ivan@154:     end;    
ivan@154: end;
ivan@154: 
ivan@154: set(0,'DefaultAxesFontSize',12);
ivan@154:     
ivan@154: figure; 
ivan@154: semilogy(squeeze(median(error(1,:,:),2)),'LineWidth',3); hold on
ivan@154: semilogy(squeeze(median(error(2,:,:),2)),'--','LineWidth',3); 
ivan@154: semilogy(squeeze(median(error(3,:,:),2)),':','LineWidth',3);
ivan@154: semilogy(squeeze(median(error(4,:,:),2)),'k','LineWidth',4);
ivan@154: semilogy(squeeze(median(error(5,:,:),2)),'k--','LineWidth',4);
ivan@154: semilogy(squeeze(median(error(6,:,:),2)),'k:','LineWidth',4);
ivan@154: semilogy(squeeze(median(error(7,:,:),2)),'k-.','LineWidth',4);
ivan@154: semilogy(squeeze(median(error(8,:,:),2)),'r--','LineWidth',3);
ivan@154: semilogy(squeeze(median(error(9,:,:),2)),'m--','LineWidth',3);
ivan@154: semilogy(squeeze(median(error(10,:,:),2)),'k--','LineWidth',3);
ivan@154: semilogy(squeeze(median(error(11,:,:),2)),'g--','LineWidth',3);
ivan@154:     
ivan@154: xlabel('[i]:iteration number')
ivan@154: ylabel('||x-x_i||')
ivan@154: axis([1 max_its + 10 10^-5 2*10^0])
ivan@154: 
ivan@154: legend(strcat('0-IHT(0) [', num2str(mean(num_iter(1,:))),'][', num2str(mean(total_time(1,:))),']'), ...
ivan@154: strcat('0-IHT(2) [', num2str(mean(num_iter(2,:))),'][', num2str(mean(total_time(2,:))),']'), ...
ivan@154: strcat('0-IHT(4) [', num2str(mean(num_iter(3,:))),'][', num2str(mean(total_time(3,:))),']'), ...
ivan@154: strcat('0-IHT(5) [', num2str(mean(num_iter(4,:))),'][', num2str(mean(total_time(4,:))),']'), ...
ivan@154: strcat('1-IHT(2) [', num2str(mean(num_iter(5,:))),'][', num2str(mean(total_time(5,:))),']'), ...
ivan@154: strcat('1-IHT(4) [', num2str(mean(num_iter(6,:))),'][', num2str(mean(total_time(6,:))),']'), ...
ivan@154: strcat('1-IHT(5) [', num2str(mean(num_iter(7,:))),'][', num2str(mean(total_time(7,:))),']'), ...
ivan@154: strcat('2-ALPS(2) [', num2str(mean(num_iter(8,:))),'][', num2str(mean(total_time(8,:))),']'), ...
ivan@154: strcat('2-ALPS(4) [', num2str(mean(num_iter(9,:))),'][', num2str(mean(total_time(9,:))),']'), ...
ivan@154: strcat('2-ALPS(5) [', num2str(mean(num_iter(10,:))),'][', num2str(mean(total_time(10,:))),']'), ...
ivan@154: strcat('3-ALPS(2) [', num2str(mean(num_iter(11,:))),'][', num2str(mean(total_time(11,:))),']'));
ivan@154: title( strcat('N= ', num2str(N), ', M=  ', num2str(M), ', K=  ', num2str(K)) );
ivan@154: shg;