annotate toolboxes/alps/ALPSbenchmark_errornorm_vs_iter.m @ 177:714fa7b8c1ad danieleb

added ramirez dl (to be completed) and MOCOD dictionary update

% =========================================================================
%                   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. [e-mail], 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;