comparison toolboxes/alps/ALPSdemo.m @ 154:0de08f68256b ivand_dev

ALPS toolbox - Algebraic Pursuit added to smallbox

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