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

ALPS toolbox - Algebraic Pursuit added to smallbox

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolboxes/alps/ALPS/AlgebraicPursuit.m	Fri Aug 12 11:17:47 2011 +0100
@@ -0,0 +1,264 @@
+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;
+