annotate toolboxes/alps/ALPS/AlgebraicPursuit.m @ 154:0de08f68256b ivand_dev

ALPS toolbox - Algebraic Pursuit added to smallbox

function [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, varargin)
% =========================================================================
%              Algebraic Pursuit algorithms - Beta Version
% =========================================================================
% 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.
% =========================================================================
% INPUT ARGUMENTS:
% y                         M x 1 undersampled measurement vector.
% Phi                       M x N regression matrix.
% K                         Sparsity of underlying vector x* or desired
%                           sparsity of solution.
% varargin                  Set of parameters. These are:
%
%    memory,...             Memory (momentum) of proposed algorithm. 
%                           Possible values are 0,1,'infty' for memoryless,
%                           one memory and infinity memory ALPS,
%                           respectively. Default value: memory = 0.
%    tol,...                Early stopping tolerance. Default value: tol =
%                           1-e5.
%    ALPSiters,...          Maximum number of algorithm iterations. Default
%                           value: 300.  
%    mode, ...              According to [1], possible values are
%                           [0,1,2,4,5,6]. This value comes from the binary 
%                           representation of the parameters:
%                           (solveNewtob, gradientDescentx, solveNewtonx), 
%                           which are explained next. Default value = 0.
%    mu,...                 Variable that controls the step size selection. 
%                           When mu = 0, step size is computed adaptively 
%                           per iteration. Default value: mu = 0. 
%    tau,...                Variable that controls the momentum in
%                           non-memoryless case. Ignored in memoryless
%                           case. User can insert as value a function handle on tau.
%                           Description given below. Default value: tau = -1. 
%    CGiters,...            Maximum iterations for Conjugate-Gradients method.
%    CGtol,...              Tolerance variable for Conjugate-Gradients method.
%    verbose                verbose = 1 prints out execution infromation.
% =========================================================================
% DESCRIPTION OF MODE VARIABLE:
%    solveNewtonb,...       If solveNewtonb == 1: Corresponds to solving a
%                           Newton system restricted to a sparse support.
%                           It is implemented via conjugate gradients. If
%                           solveNewtonb == 0: Step size selection as
%                           described in eqs. (12) and (13) in [1]. Default
%                           value: solveNewtonb = 0. For infty memory case,
%                           solveNewtonb = 0.
%    gradientDescentx,...   If gradientDescentx == 1: single gradient
%                           update of x_{i+1} restricted ot its support
%                           with line search. Default value:
%                           gradientDescentx = 1.
%    solveNewtonx,...       If solveNewtonx == 1: Akin to Hard Thresholding
%                           Pursuit (c.f. Simon Foucart, "Hard Thresholding
%                           Pursuit," preprint, 2010). Default vale:
%                           solveNewtonx = 0.
% =========================================================================
% DESCRIPTION OF SPECIAL TAU VALUES:
%    tau = -1,...           Tau computed as the minimized of the objective 
%                           function:
%
%                                  <u - Ax_i, Phi(x_i - x_{i-1})>
%                             tau = -----------------------------.
%                                  <u - Ax_i, Phi(x_i - x_{i-1})>
%
%    tau = 0,...            Follows FISTA weight configuration at each
%                           iteration. For more information, please read "A 
%                           Fast Iterative Shrinkage-Thresholding Algorithm
%                           for Linear Inverse Problems", written by A.
%                           Beck and M. Teboulle, SIAM J. Imaging Sciences,
%                           vol. 2, no. 1, 2009.
% =========================================================================
% OUTPUT ARGUMENTS:
% x_hat                     N x 1 recovered K-sparse vector.
% numiter                   Number of iterations executed.
% time                      Execution time in seconds.
% x_path                    Keeps a series of computed N x 1 K-sparse vectors 
%                           until the end of the iterative process.
% =========================================================================
% 01/04/2011, by Anastasios Kyrillidis. [e-mail], EPFL.
% =========================================================================
% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005.
%                         Email: [e-mail]
% =========================================================================
% This work was supported in part by the European Commission under Grant 
% MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like 
% to acknowledge Rice University for his Faculty Fellowship.
% =========================================================================

x_hat = 0;
time = 0;
x_path = 0;

params.memory = 0;
params.tol = 1e-5;
params.ALPSiters = 300;
params.mode = 0;
params.tau = 1/2;
params.memory = 0;
params.cg_maxiter = 50;
params.cg_tol = 1e-8;
params.useCG = 0;
params.verbose = 0;
params.mu = 0;

for in_arg = 1:2:length(varargin)
    if (strcmp(varargin{in_arg}, 'memory'))
        params.memory = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'mode'))
        params.mode = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'tolerance'))
        params.tol = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'ALPSiterations'))
        params.ALPSiters = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'tau'))
        params.tau = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'mu'))
        params.mu = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'useCG'))
        params.useCG = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'CGiters'))
        params.cg_maxiter = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'CGtol'))
        params.cg_tol = varargin{in_arg+1};
    end
    if (strcmp(varargin{in_arg}, 'verbose'))
        params.verbose = varargin{in_arg+1};
    end;
end;

%% Check input

if (ischar(params.memory))
    if ~(sum(params.memory == 'infty') == 5)
        disp('ERROR: Memory variable takes positive integer values or "infty" value.');
        numiter = -1;
        return;
    end;
elseif (params.memory < 0 || ~isinteger(uint8(params.memory)))
    disp('ERROR: Memory variable takes positive integer values or "infty" value.');
    numiter = -1;
    return;
end;

if (params.mode ~= 0 && params.mode ~= 1 && params.mode ~= 2 && params.mode ~= 4 && params.mode ~= 5 && params.mode ~= 6)
    disp('ERROR: Mode variable takes values from range [0,1,2,4,5,6].');
    numiter = -1;
    return;
end;

if (params.tol < 0 || ~isnumeric(params.tol))
    disp('ERROR: tolerance variable must be positive.');
    numiter = -1;
    return;
end;


if (params.mode == 0 || params.mode == 2)
    params.useCG = 0;
end;

y = y(:);

M_y = length(y);
[M_Phi, tmpArg] = size(Phi);

if (params.mode == 4 || params.mode == 5 || params.mode == 6)
    if (M_Phi < 2*K) 
        disp('WARNING: Newton system may be ill-conditioned... Press any key to continue');
        pause;
    end;
end;

[params.solveNewtonb, params.gradientDescentx, params.solveNewtonx] = configuration(params.mode);

if (M_y ~= M_Phi)
    error('Inconsistent sampling matrix...');
end;

switch params.memory
    case 0,        
        tic;
        [x_hat, numiter, x_path] = zero_ALPS(y, Phi, K, params);
        time = toc;
        if (params.verbose == 1)
            str = sprintf('%d-ALPS(%d) time: %s', params.memory, params.mode, num2str(time));
            disp(str);
            str = sprintf('Number of iterations: %d', numiter);
            disp(str);
        end;
    case 1,
        tic;
        [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params);
        time = toc;
        if (params.verbose == 1)
            if (isa(params.tau,'float'))
                if (params.tau == 0)
                    str = sprintf('%d-ALPS(%d) - fista tau - time: %s', params.memory, params.mode, num2str(time));
                elseif (params.tau == -1)
                    str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time));
                else
                    str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time));
                end;
            elseif (isa(params.tau, 'function_handle'))
                str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time));
            end;
            disp(str);
            str = sprintf('Number of iterations: %d', numiter);
            disp(str);
        end;   
    case 'infty',
        tic;
        [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params);
        time = toc;
        if (params.verbose == 1)
            if (isa(params.tau,'float'))
                if (params.tau == -1)
                    str = sprintf('infty-ALPS(%d) - optimized tau - time: %s', params.mode, num2str(time));
                else
                    str = sprintf('infty-ALPS(%d) - tau = %.3f - time: %s', params.mode, params.tau, num2str(time));
                end;
            elseif (isa(params.tau,'function_handle'))
                str = sprintf('infty-ALPS(%d) - user derfined tau - time: %s', params.mode, num2str(time));
            end;
            disp(str);
            str = sprintf('Number of iterations: %d', numiter);
            disp(str);   
        end;        
    otherwise
        tic;
        [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params);
        time = toc;
        if (params.verbose == 1)
            if (isa(params.tau,'float'))
                if (params.tau == -1)
                    str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time));
                else
                    str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time));
                end;
            elseif (isa(params.tau,'function_handle'))
                str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time));
            end;
            disp(str);
            str = sprintf('Number of iterations: %d', numiter);
            disp(str);        
        end;
end;