annotate toolboxes/alps/ALPSdemo.m @ 163:855025f4c779 ivand_dev

renaiming small_cgp to small_pcgp
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Wed, 07 Sep 2011 14:16:50 +0100
parents 0de08f68256b
children
rev   line source
ivan@154 1 % =========================================================================
ivan@154 2 % Algebraic Pursuit algorithms - Demo
ivan@154 3 % =========================================================================
ivan@154 4 % Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods
ivan@154 5 % on sparse recovery of linear inverse systems. In particular, let
ivan@154 6 % Phi : M x N real matrix (M < N)
ivan@154 7 % x* : N x 1 K-sparse data vector
ivan@154 8 % n : N x 1 additive noise vector
ivan@154 9 % then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS
ivan@154 10 % solve the following minimization problem
ivan@154 11 % minimize ||y - Phi x||^2 subject to x is K-sparse vector.
ivan@154 12 %
ivan@154 13 % Detailed discussion on the algorithm can be found in
ivan@154 14 % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
ivan@154 15 % by Volkan Cevher, Technical Report, 2011.
ivan@154 16 % =========================================================================
ivan@154 17 % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
ivan@154 18 % =========================================================================
ivan@154 19
ivan@154 20 clear all
ivan@154 21 close all
ivan@154 22 clc
ivan@154 23
ivan@154 24 addpath ALPS/;
ivan@154 25
ivan@154 26 % General parameters
ivan@154 27 N = 1000; % Size of sparse vector x
ivan@154 28 M = 300; % Number of measurements
ivan@154 29 K = ceil(M*(0.25)); % Sparsity of x*
ivan@154 30 sigma = 0.000; % Additive noise variance
ivan@154 31 verbose = 1; % Print information
ivan@154 32
ivan@154 33 % ALPS parameters
ivan@154 34 its = 600; % Maximum number of iterations
ivan@154 35 Ttol= 1e-6; % Convergence tolerance
ivan@154 36 mode = 1; % Binary representation of (SolveNewtonb, GradientDescentx, SolveNewtox)
ivan@154 37 memory = 1; % Memory of IHT method
ivan@154 38 tau = 0; % Momentum term
ivan@154 39 useCG = 0; % Usage of Conjugate gradients method
ivan@154 40 mu = 0; % Step size selection
ivan@154 41
ivan@154 42 % Please refer to generate_matrix.m and generate_vector.m files
ivan@154 43 m_ensemble = 'Gaussian'; % Measurement matrix ensemble type
ivan@154 44 x_ensemble = 'Gaussian'; % Sparse vector ensemble type
ivan@154 45
ivan@154 46 % reset(RandStream.getDefaultStream);
ivan@154 47
ivan@154 48 %% Generate data
ivan@154 49 x = generate_vector(N, K, x_ensemble);
ivan@154 50 Phi = generate_matrix(M, N, m_ensemble);
ivan@154 51 noise = sigma*randn(M,1);
ivan@154 52 y = Phi*x + noise;
ivan@154 53
ivan@154 54 %% ALPS reconstruction
ivan@154 55 disp('=====================================================================');
ivan@154 56 [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 57 if (numiter ~= -1)
ivan@154 58 str = sprintf('Error norm: %f ', norm(x-x_hat));
ivan@154 59 disp(str);
ivan@154 60 end;
ivan@154 61
ivan@154 62 %% Plot error per iteration
ivan@154 63 if (numiter ~= -1)
ivan@154 64 error = sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
ivan@154 65
ivan@154 66 set(0, 'DefaultAxesFontSize', 14);
ivan@154 67 figure;
ivan@154 68 semilogy(error,'LineWidth',3); hold on;
ivan@154 69 grid on;
ivan@154 70 semilogy((norm(noise))*ones(1,its),'m-','LineWidth',1);
ivan@154 71
ivan@154 72 xlabel('[i]:iteration number')
ivan@154 73 ylabel('||x-x_i||')
ivan@154 74 axis([1 numiter + 5 10^-5 2*10^0])
ivan@154 75
ivan@154 76 legend(strcat('Time = [',num2str(time),']'), 'Noise Level');
ivan@154 77 title(strcat('N = ', num2str(N), ', M = ', num2str(M), ', K = ', num2str(K)));
ivan@154 78 shg;
ivan@154 79 end;