annotate toolboxes/alps/ALPS/AlgebraicPursuit.m @ 247:ecce33192fcc tip

Added tag ver_2.1 for changeset cef4500b936f
author luisf <luis.figueira@eecs.qmul.ac.uk>
date Wed, 31 Oct 2012 12:24:44 +0000
parents 0de08f68256b
children
rev   line source
ivan@154 1 function [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, varargin)
ivan@154 2 % =========================================================================
ivan@154 3 % Algebraic Pursuit algorithms - Beta Version
ivan@154 4 % =========================================================================
ivan@154 5 % Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods
ivan@154 6 % on sparse recovery of linear inverse systems. In particular, let
ivan@154 7 % Phi : M x N real matrix (M < N)
ivan@154 8 % x* : N x 1 K-sparse data vector
ivan@154 9 % n : N x 1 additive noise vector
ivan@154 10 % then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS
ivan@154 11 % solve the following minimization problem
ivan@154 12 % minimize ||y - Phi x||^2 subject to x is K-sparse vector.
ivan@154 13 %
ivan@154 14 % Detailed discussion on the algorithm can be found in
ivan@154 15 % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
ivan@154 16 % by Volkan Cevher, Technical Report, 2011.
ivan@154 17 % =========================================================================
ivan@154 18 % INPUT ARGUMENTS:
ivan@154 19 % y M x 1 undersampled measurement vector.
ivan@154 20 % Phi M x N regression matrix.
ivan@154 21 % K Sparsity of underlying vector x* or desired
ivan@154 22 % sparsity of solution.
ivan@154 23 % varargin Set of parameters. These are:
ivan@154 24 %
ivan@154 25 % memory,... Memory (momentum) of proposed algorithm.
ivan@154 26 % Possible values are 0,1,'infty' for memoryless,
ivan@154 27 % one memory and infinity memory ALPS,
ivan@154 28 % respectively. Default value: memory = 0.
ivan@154 29 % tol,... Early stopping tolerance. Default value: tol =
ivan@154 30 % 1-e5.
ivan@154 31 % ALPSiters,... Maximum number of algorithm iterations. Default
ivan@154 32 % value: 300.
ivan@154 33 % mode, ... According to [1], possible values are
ivan@154 34 % [0,1,2,4,5,6]. This value comes from the binary
ivan@154 35 % representation of the parameters:
ivan@154 36 % (solveNewtob, gradientDescentx, solveNewtonx),
ivan@154 37 % which are explained next. Default value = 0.
ivan@154 38 % mu,... Variable that controls the step size selection.
ivan@154 39 % When mu = 0, step size is computed adaptively
ivan@154 40 % per iteration. Default value: mu = 0.
ivan@154 41 % tau,... Variable that controls the momentum in
ivan@154 42 % non-memoryless case. Ignored in memoryless
ivan@154 43 % case. User can insert as value a function handle on tau.
ivan@154 44 % Description given below. Default value: tau = -1.
ivan@154 45 % CGiters,... Maximum iterations for Conjugate-Gradients method.
ivan@154 46 % CGtol,... Tolerance variable for Conjugate-Gradients method.
ivan@154 47 % verbose verbose = 1 prints out execution infromation.
ivan@154 48 % =========================================================================
ivan@154 49 % DESCRIPTION OF MODE VARIABLE:
ivan@154 50 % solveNewtonb,... If solveNewtonb == 1: Corresponds to solving a
ivan@154 51 % Newton system restricted to a sparse support.
ivan@154 52 % It is implemented via conjugate gradients. If
ivan@154 53 % solveNewtonb == 0: Step size selection as
ivan@154 54 % described in eqs. (12) and (13) in [1]. Default
ivan@154 55 % value: solveNewtonb = 0. For infty memory case,
ivan@154 56 % solveNewtonb = 0.
ivan@154 57 % gradientDescentx,... If gradientDescentx == 1: single gradient
ivan@154 58 % update of x_{i+1} restricted ot its support
ivan@154 59 % with line search. Default value:
ivan@154 60 % gradientDescentx = 1.
ivan@154 61 % solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding
ivan@154 62 % Pursuit (c.f. Simon Foucart, "Hard Thresholding
ivan@154 63 % Pursuit," preprint, 2010). Default vale:
ivan@154 64 % solveNewtonx = 0.
ivan@154 65 % =========================================================================
ivan@154 66 % DESCRIPTION OF SPECIAL TAU VALUES:
ivan@154 67 % tau = -1,... Tau computed as the minimized of the objective
ivan@154 68 % function:
ivan@154 69 %
ivan@154 70 % <u - Ax_i, Phi(x_i - x_{i-1})>
ivan@154 71 % tau = -----------------------------.
ivan@154 72 % <u - Ax_i, Phi(x_i - x_{i-1})>
ivan@154 73 %
ivan@154 74 % tau = 0,... Follows FISTA weight configuration at each
ivan@154 75 % iteration. For more information, please read "A
ivan@154 76 % Fast Iterative Shrinkage-Thresholding Algorithm
ivan@154 77 % for Linear Inverse Problems", written by A.
ivan@154 78 % Beck and M. Teboulle, SIAM J. Imaging Sciences,
ivan@154 79 % vol. 2, no. 1, 2009.
ivan@154 80 % =========================================================================
ivan@154 81 % OUTPUT ARGUMENTS:
ivan@154 82 % x_hat N x 1 recovered K-sparse vector.
ivan@154 83 % numiter Number of iterations executed.
ivan@154 84 % time Execution time in seconds.
ivan@154 85 % x_path Keeps a series of computed N x 1 K-sparse vectors
ivan@154 86 % until the end of the iterative process.
ivan@154 87 % =========================================================================
ivan@154 88 % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
ivan@154 89 % =========================================================================
ivan@154 90 % cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005.
ivan@154 91 % Email: jrom@acm.caltech.edu
ivan@154 92 % =========================================================================
ivan@154 93 % This work was supported in part by the European Commission under Grant
ivan@154 94 % MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like
ivan@154 95 % to acknowledge Rice University for his Faculty Fellowship.
ivan@154 96 % =========================================================================
ivan@154 97
ivan@154 98 x_hat = 0;
ivan@154 99 time = 0;
ivan@154 100 x_path = 0;
ivan@154 101
ivan@154 102 params.memory = 0;
ivan@154 103 params.tol = 1e-5;
ivan@154 104 params.ALPSiters = 300;
ivan@154 105 params.mode = 0;
ivan@154 106 params.tau = 1/2;
ivan@154 107 params.memory = 0;
ivan@154 108 params.cg_maxiter = 50;
ivan@154 109 params.cg_tol = 1e-8;
ivan@154 110 params.useCG = 0;
ivan@154 111 params.verbose = 0;
ivan@154 112 params.mu = 0;
ivan@154 113
ivan@154 114 for in_arg = 1:2:length(varargin)
ivan@154 115 if (strcmp(varargin{in_arg}, 'memory'))
ivan@154 116 params.memory = varargin{in_arg+1};
ivan@154 117 end
ivan@154 118 if (strcmp(varargin{in_arg}, 'mode'))
ivan@154 119 params.mode = varargin{in_arg+1};
ivan@154 120 end
ivan@154 121 if (strcmp(varargin{in_arg}, 'tolerance'))
ivan@154 122 params.tol = varargin{in_arg+1};
ivan@154 123 end
ivan@154 124 if (strcmp(varargin{in_arg}, 'ALPSiterations'))
ivan@154 125 params.ALPSiters = varargin{in_arg+1};
ivan@154 126 end
ivan@154 127 if (strcmp(varargin{in_arg}, 'tau'))
ivan@154 128 params.tau = varargin{in_arg+1};
ivan@154 129 end
ivan@154 130 if (strcmp(varargin{in_arg}, 'mu'))
ivan@154 131 params.mu = varargin{in_arg+1};
ivan@154 132 end
ivan@154 133 if (strcmp(varargin{in_arg}, 'useCG'))
ivan@154 134 params.useCG = varargin{in_arg+1};
ivan@154 135 end
ivan@154 136 if (strcmp(varargin{in_arg}, 'CGiters'))
ivan@154 137 params.cg_maxiter = varargin{in_arg+1};
ivan@154 138 end
ivan@154 139 if (strcmp(varargin{in_arg}, 'CGtol'))
ivan@154 140 params.cg_tol = varargin{in_arg+1};
ivan@154 141 end
ivan@154 142 if (strcmp(varargin{in_arg}, 'verbose'))
ivan@154 143 params.verbose = varargin{in_arg+1};
ivan@154 144 end;
ivan@154 145 end;
ivan@154 146
ivan@154 147 %% Check input
ivan@154 148
ivan@154 149 if (ischar(params.memory))
ivan@154 150 if ~(sum(params.memory == 'infty') == 5)
ivan@154 151 disp('ERROR: Memory variable takes positive integer values or "infty" value.');
ivan@154 152 numiter = -1;
ivan@154 153 return;
ivan@154 154 end;
ivan@154 155 elseif (params.memory < 0 || ~isinteger(uint8(params.memory)))
ivan@154 156 disp('ERROR: Memory variable takes positive integer values or "infty" value.');
ivan@154 157 numiter = -1;
ivan@154 158 return;
ivan@154 159 end;
ivan@154 160
ivan@154 161 if (params.mode ~= 0 && params.mode ~= 1 && params.mode ~= 2 && params.mode ~= 4 && params.mode ~= 5 && params.mode ~= 6)
ivan@154 162 disp('ERROR: Mode variable takes values from range [0,1,2,4,5,6].');
ivan@154 163 numiter = -1;
ivan@154 164 return;
ivan@154 165 end;
ivan@154 166
ivan@154 167 if (params.tol < 0 || ~isnumeric(params.tol))
ivan@154 168 disp('ERROR: tolerance variable must be positive.');
ivan@154 169 numiter = -1;
ivan@154 170 return;
ivan@154 171 end;
ivan@154 172
ivan@154 173
ivan@154 174 if (params.mode == 0 || params.mode == 2)
ivan@154 175 params.useCG = 0;
ivan@154 176 end;
ivan@154 177
ivan@154 178 y = y(:);
ivan@154 179
ivan@154 180 M_y = length(y);
ivan@154 181 [M_Phi, tmpArg] = size(Phi);
ivan@154 182
ivan@154 183 if (params.mode == 4 || params.mode == 5 || params.mode == 6)
ivan@154 184 if (M_Phi < 2*K)
ivan@154 185 disp('WARNING: Newton system may be ill-conditioned... Press any key to continue');
ivan@154 186 pause;
ivan@154 187 end;
ivan@154 188 end;
ivan@154 189
ivan@154 190 [params.solveNewtonb, params.gradientDescentx, params.solveNewtonx] = configuration(params.mode);
ivan@154 191
ivan@154 192 if (M_y ~= M_Phi)
ivan@154 193 error('Inconsistent sampling matrix...');
ivan@154 194 end;
ivan@154 195
ivan@154 196 switch params.memory
ivan@154 197 case 0,
ivan@154 198 tic;
ivan@154 199 [x_hat, numiter, x_path] = zero_ALPS(y, Phi, K, params);
ivan@154 200 time = toc;
ivan@154 201 if (params.verbose == 1)
ivan@154 202 str = sprintf('%d-ALPS(%d) time: %s', params.memory, params.mode, num2str(time));
ivan@154 203 disp(str);
ivan@154 204 str = sprintf('Number of iterations: %d', numiter);
ivan@154 205 disp(str);
ivan@154 206 end;
ivan@154 207 case 1,
ivan@154 208 tic;
ivan@154 209 [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params);
ivan@154 210 time = toc;
ivan@154 211 if (params.verbose == 1)
ivan@154 212 if (isa(params.tau,'float'))
ivan@154 213 if (params.tau == 0)
ivan@154 214 str = sprintf('%d-ALPS(%d) - fista tau - time: %s', params.memory, params.mode, num2str(time));
ivan@154 215 elseif (params.tau == -1)
ivan@154 216 str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time));
ivan@154 217 else
ivan@154 218 str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time));
ivan@154 219 end;
ivan@154 220 elseif (isa(params.tau, 'function_handle'))
ivan@154 221 str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time));
ivan@154 222 end;
ivan@154 223 disp(str);
ivan@154 224 str = sprintf('Number of iterations: %d', numiter);
ivan@154 225 disp(str);
ivan@154 226 end;
ivan@154 227 case 'infty',
ivan@154 228 tic;
ivan@154 229 [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params);
ivan@154 230 time = toc;
ivan@154 231 if (params.verbose == 1)
ivan@154 232 if (isa(params.tau,'float'))
ivan@154 233 if (params.tau == -1)
ivan@154 234 str = sprintf('infty-ALPS(%d) - optimized tau - time: %s', params.mode, num2str(time));
ivan@154 235 else
ivan@154 236 str = sprintf('infty-ALPS(%d) - tau = %.3f - time: %s', params.mode, params.tau, num2str(time));
ivan@154 237 end;
ivan@154 238 elseif (isa(params.tau,'function_handle'))
ivan@154 239 str = sprintf('infty-ALPS(%d) - user derfined tau - time: %s', params.mode, num2str(time));
ivan@154 240 end;
ivan@154 241 disp(str);
ivan@154 242 str = sprintf('Number of iterations: %d', numiter);
ivan@154 243 disp(str);
ivan@154 244 end;
ivan@154 245 otherwise
ivan@154 246 tic;
ivan@154 247 [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params);
ivan@154 248 time = toc;
ivan@154 249 if (params.verbose == 1)
ivan@154 250 if (isa(params.tau,'float'))
ivan@154 251 if (params.tau == -1)
ivan@154 252 str = sprintf('%d-ALPS(%d) - optimized tau - time: %s', params.memory, params.mode, num2str(time));
ivan@154 253 else
ivan@154 254 str = sprintf('%d-ALPS(%d) - tau = %.3f - time: %s', params.memory, params.mode, params.tau, num2str(time));
ivan@154 255 end;
ivan@154 256 elseif (isa(params.tau,'function_handle'))
ivan@154 257 str = sprintf('%d-ALPS(%d) - user derfined tau - time: %s', params.memory, params.mode, num2str(time));
ivan@154 258 end;
ivan@154 259 disp(str);
ivan@154 260 str = sprintf('Number of iterations: %d', numiter);
ivan@154 261 disp(str);
ivan@154 262 end;
ivan@154 263 end;
ivan@154 264