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 = 1000;                       % Size of sparse vector x
ivan@154: M = 300;                        % Number of measurements
ivan@154: K = ceil(M*(0.25));              % Sparsity of  x*
ivan@154: sigma = 0.000;                 % 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: mode = 1;                       % Binary representation of (SolveNewtonb, GradientDescentx, SolveNewtox)
ivan@154: memory = 1;                     % Memory of IHT method
ivan@154: tau = 0;                        % Momentum term
ivan@154: useCG = 0;                      % Usage of Conjugate gradients method
ivan@154: mu = 0;                         % Step size selection
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: % reset(RandStream.getDefaultStream);
ivan@154: 
ivan@154: %% Generate data
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: %% ALPS reconstruction
ivan@154: disp('=====================================================================');
ivan@154: [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);
ivan@154: if (numiter ~= -1)
ivan@154:     str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154:     disp(str);
ivan@154: end;
ivan@154: 
ivan@154: %% Plot error per iteration
ivan@154: if (numiter ~= -1)
ivan@154:     error = sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154: 
ivan@154:     set(0, 'DefaultAxesFontSize', 14);
ivan@154:     figure;
ivan@154:     semilogy(error,'LineWidth',3); hold on;
ivan@154:     grid on;
ivan@154:     semilogy((norm(noise))*ones(1,its),'m-','LineWidth',1);
ivan@154: 
ivan@154:     xlabel('[i]:iteration number')
ivan@154:     ylabel('||x-x_i||')
ivan@154:     axis([1 numiter + 5 10^-5 2*10^0])
ivan@154: 
ivan@154:     legend(strcat('Time = [',num2str(time),']'), 'Noise Level');
ivan@154:     title(strcat('N = ', num2str(N), ', M = ', num2str(M), ', K = ', num2str(K)));
ivan@154:     shg;
ivan@154: end;