view toolboxes/alps/ALPSdemo.m @ 174:dc2f0fa21310 danieleb

multiple trials with error bars

% =========================================================================
%                   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 = 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;