ivan@154: function [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, varargin) ivan@154: % ========================================================================= ivan@154: % Algebraic Pursuit algorithms - Beta Version 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: % INPUT ARGUMENTS: ivan@154: % y M x 1 undersampled measurement vector. ivan@154: % Phi M x N regression matrix. ivan@154: % K Sparsity of underlying vector x* or desired ivan@154: % sparsity of solution. ivan@154: % varargin Set of parameters. These are: ivan@154: % ivan@154: % memory,... Memory (momentum) of proposed algorithm. ivan@154: % Possible values are 0,1,'infty' for memoryless, ivan@154: % one memory and infinity memory ALPS, ivan@154: % respectively. Default value: memory = 0. ivan@154: % tol,... Early stopping tolerance. Default value: tol = ivan@154: % 1-e5. ivan@154: % ALPSiters,... Maximum number of algorithm iterations. Default ivan@154: % value: 300. ivan@154: % mode, ... According to [1], possible values are ivan@154: % [0,1,2,4,5,6]. This value comes from the binary ivan@154: % representation of the parameters: ivan@154: % (solveNewtob, gradientDescentx, solveNewtonx), ivan@154: % which are explained next. Default value = 0. ivan@154: % mu,... Variable that controls the step size selection. ivan@154: % When mu = 0, step size is computed adaptively ivan@154: % per iteration. Default value: mu = 0. ivan@154: % tau,... Variable that controls the momentum in ivan@154: % non-memoryless case. Ignored in memoryless ivan@154: % case. User can insert as value a function handle on tau. ivan@154: % Description given below. Default value: tau = -1. ivan@154: % CGiters,... Maximum iterations for Conjugate-Gradients method. ivan@154: % CGtol,... Tolerance variable for Conjugate-Gradients method. ivan@154: % verbose verbose = 1 prints out execution infromation. ivan@154: % ========================================================================= ivan@154: % DESCRIPTION OF MODE VARIABLE: ivan@154: % solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a ivan@154: % Newton system restricted to a sparse support. ivan@154: % It is implemented via conjugate gradients. If ivan@154: % solveNewtonb == 0: Step size selection as ivan@154: % described in eqs. (12) and (13) in [1]. Default ivan@154: % value: solveNewtonb = 0. For infty memory case, ivan@154: % solveNewtonb = 0. ivan@154: % gradientDescentx,... If gradientDescentx == 1: single gradient ivan@154: % update of x_{i+1} restricted ot its support ivan@154: % with line search. Default value: ivan@154: % gradientDescentx = 1. ivan@154: % solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding ivan@154: % Pursuit (c.f. Simon Foucart, "Hard Thresholding ivan@154: % Pursuit," preprint, 2010). Default vale: ivan@154: % solveNewtonx = 0. ivan@154: % ========================================================================= ivan@154: % DESCRIPTION OF SPECIAL TAU VALUES: ivan@154: % tau = -1,... Tau computed as the minimized of the objective ivan@154: % function: ivan@154: % ivan@154: % ivan@154: % tau = -----------------------------. ivan@154: % ivan@154: % ivan@154: % tau = 0,... Follows FISTA weight configuration at each ivan@154: % iteration. For more information, please read "A ivan@154: % Fast Iterative Shrinkage-Thresholding Algorithm ivan@154: % for Linear Inverse Problems", written by A. ivan@154: % Beck and M. Teboulle, SIAM J. Imaging Sciences, ivan@154: % vol. 2, no. 1, 2009. ivan@154: % ========================================================================= ivan@154: % OUTPUT ARGUMENTS: ivan@154: % x_hat N x 1 recovered K-sparse vector. ivan@154: % numiter Number of iterations executed. ivan@154: % time Execution time in seconds. ivan@154: % x_path Keeps a series of computed N x 1 K-sparse vectors ivan@154: % until the end of the iterative process. ivan@154: % ========================================================================= ivan@154: % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. ivan@154: % ========================================================================= ivan@154: % cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. ivan@154: % Email: jrom@acm.caltech.edu ivan@154: % ========================================================================= ivan@154: % This work was supported in part by the European Commission under Grant ivan@154: % MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like ivan@154: % to acknowledge Rice University for his Faculty Fellowship. ivan@154: % ========================================================================= ivan@154: ivan@154: x_hat = 0; ivan@154: time = 0; ivan@154: x_path = 0; ivan@154: ivan@154: params.memory = 0; ivan@154: params.tol = 1e-5; ivan@154: params.ALPSiters = 300; ivan@154: params.mode = 0; ivan@154: params.tau = 1/2; ivan@154: params.memory = 0; ivan@154: params.cg_maxiter = 50; ivan@154: params.cg_tol = 1e-8; ivan@154: params.useCG = 0; ivan@154: params.verbose = 0; ivan@154: params.mu = 0; ivan@154: ivan@154: for in_arg = 1:2:length(varargin) ivan@154: if (strcmp(varargin{in_arg}, 'memory')) ivan@154: params.memory = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'mode')) ivan@154: params.mode = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'tolerance')) ivan@154: params.tol = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'ALPSiterations')) ivan@154: params.ALPSiters = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'tau')) ivan@154: params.tau = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'mu')) ivan@154: params.mu = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'useCG')) ivan@154: params.useCG = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'CGiters')) ivan@154: params.cg_maxiter = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'CGtol')) ivan@154: params.cg_tol = varargin{in_arg+1}; ivan@154: end ivan@154: if (strcmp(varargin{in_arg}, 'verbose')) ivan@154: params.verbose = varargin{in_arg+1}; ivan@154: end; ivan@154: end; ivan@154: ivan@154: %% Check input ivan@154: ivan@154: if (ischar(params.memory)) ivan@154: if ~(sum(params.memory == 'infty') == 5) ivan@154: disp('ERROR: Memory variable takes positive integer values or "infty" value.'); ivan@154: numiter = -1; ivan@154: return; ivan@154: end; ivan@154: elseif (params.memory < 0 || ~isinteger(uint8(params.memory))) ivan@154: disp('ERROR: Memory variable takes positive integer values or "infty" value.'); ivan@154: numiter = -1; ivan@154: return; ivan@154: end; ivan@154: ivan@154: if (params.mode ~= 0 && params.mode ~= 1 && params.mode ~= 2 && params.mode ~= 4 && params.mode ~= 5 && params.mode ~= 6) ivan@154: disp('ERROR: Mode variable takes values from range [0,1,2,4,5,6].'); ivan@154: numiter = -1; ivan@154: return; ivan@154: end; ivan@154: ivan@154: if (params.tol < 0 || ~isnumeric(params.tol)) ivan@154: disp('ERROR: tolerance variable must be positive.'); ivan@154: numiter = -1; ivan@154: return; ivan@154: end; ivan@154: ivan@154: ivan@154: if (params.mode == 0 || params.mode == 2) ivan@154: params.useCG = 0; ivan@154: end; ivan@154: ivan@154: y = y(:); ivan@154: ivan@154: M_y = length(y); ivan@154: [M_Phi, tmpArg] = size(Phi); ivan@154: ivan@154: if (params.mode == 4 || params.mode == 5 || params.mode == 6) ivan@154: if (M_Phi < 2*K) ivan@154: disp('WARNING: Newton system may be ill-conditioned... Press any key to continue'); ivan@154: pause; ivan@154: end; ivan@154: end; ivan@154: ivan@154: [params.solveNewtonb, params.gradientDescentx, params.solveNewtonx] = configuration(params.mode); ivan@154: ivan@154: if (M_y ~= M_Phi) ivan@154: error('Inconsistent sampling matrix...'); ivan@154: end; ivan@154: ivan@154: switch params.memory ivan@154: case 0, ivan@154: tic; ivan@154: [x_hat, numiter, x_path] = zero_ALPS(y, Phi, K, params); ivan@154: time = toc; ivan@154: if (params.verbose == 1) ivan@154: str = sprintf('%d-ALPS(%d) time: %s', params.memory, params.mode, num2str(time)); ivan@154: disp(str); ivan@154: str = sprintf('Number of iterations: %d', numiter); ivan@154: disp(str); ivan@154: end; ivan@154: case 1, ivan@154: tic; ivan@154: [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params); ivan@154: time = toc; ivan@154: if (params.verbose == 1) ivan@154: if (isa(params.tau,'float')) ivan@154: if (params.tau == 0) ivan@154: str = sprintf('%d-ALPS(%d) - fista tau - time: %s', params.memory, params.mode, num2str(time)); ivan@154: elseif (params.tau == -1) ivan@154: str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time)); ivan@154: else ivan@154: str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time)); ivan@154: end; ivan@154: elseif (isa(params.tau, 'function_handle')) ivan@154: str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time)); ivan@154: end; ivan@154: disp(str); ivan@154: str = sprintf('Number of iterations: %d', numiter); ivan@154: disp(str); ivan@154: end; ivan@154: case 'infty', ivan@154: tic; ivan@154: [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params); ivan@154: time = toc; ivan@154: if (params.verbose == 1) ivan@154: if (isa(params.tau,'float')) ivan@154: if (params.tau == -1) ivan@154: str = sprintf('infty-ALPS(%d) - optimized tau - time: %s', params.mode, num2str(time)); ivan@154: else ivan@154: str = sprintf('infty-ALPS(%d) - tau = %.3f - time: %s', params.mode, params.tau, num2str(time)); ivan@154: end; ivan@154: elseif (isa(params.tau,'function_handle')) ivan@154: str = sprintf('infty-ALPS(%d) - user derfined tau - time: %s', params.mode, num2str(time)); ivan@154: end; ivan@154: disp(str); ivan@154: str = sprintf('Number of iterations: %d', numiter); ivan@154: disp(str); ivan@154: end; ivan@154: otherwise ivan@154: tic; ivan@154: [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params); ivan@154: time = toc; ivan@154: if (params.verbose == 1) ivan@154: if (isa(params.tau,'float')) ivan@154: if (params.tau == -1) ivan@154: str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time)); ivan@154: else ivan@154: str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time)); ivan@154: end; ivan@154: elseif (isa(params.tau,'function_handle')) ivan@154: str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time)); ivan@154: end; ivan@154: disp(str); ivan@154: str = sprintf('Number of iterations: %d', numiter); ivan@154: disp(str); ivan@154: end; ivan@154: end; ivan@154: