annotate toolboxes/alps/ALPS/infty_ALPS.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, x_path] = infty_ALPS(y, Phi, K, params)
ivan@154 2 % =========================================================================
ivan@154 3 % infty-ALPS(#) algorithm - Beta Version
ivan@154 4 % =========================================================================
ivan@154 5 % Algebraic Pursuit (ALPS) algorithm with infty-memory acceleration.
ivan@154 6 %
ivan@154 7 % Detailed discussion on the algorithm can be found in
ivan@154 8 % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
ivan@154 9 % by Volkan Cevher, Technical Report, 2011.
ivan@154 10 % =========================================================================
ivan@154 11 % INPUT ARGUMENTS:
ivan@154 12 % y M x 1 undersampled measurement vector.
ivan@154 13 % Phi M x N regression matrix.
ivan@154 14 % K Sparsity of underlying vector x* or desired
ivan@154 15 % sparsity of solution.
ivan@154 16 % params Structure of parameters. These are:
ivan@154 17 %
ivan@154 18 % tol,... Early stopping tolerance. Default value: tol =
ivan@154 19 % 1-e5
ivan@154 20 % ALPSiters,... Maximum number of algorithm iterations. Default
ivan@154 21 % value: 300.
ivan@154 22 % solveNewtonb,... Value: solveNewtonb = 0. Not used in infty
ivan@154 23 % methods.
ivan@154 24 % gradientDescentx,... If gradientDescentx == 1: single gradient
ivan@154 25 % update of x_{i+1} restricted ot its support with
ivan@154 26 % line search. Default value: gradientDescentx =
ivan@154 27 % 1.
ivan@154 28 % solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit
ivan@154 29 % (c.f. Simon Foucart, "Hard Thresholding Pursuit,"
ivan@154 30 % preprint, 2010). Default vale: solveNewtonx = 0
ivan@154 31 % tau,... Variable that controls the momentum in
ivan@154 32 % non-memoryless case. Ignored in memoryless
ivan@154 33 % case. Default value: tau = 1/2.
ivan@154 34 % Special cases:
ivan@154 35 % - tau = -1: momentum step size is automatically
ivan@154 36 % optimized in every step.
ivan@154 37 % - tau as a function handle: user defined
ivan@154 38 % behavior of tau momentum term.
ivan@154 39 % mu,... Variable that controls the step size selection.
ivan@154 40 % When mu = 0, step size is computed adaptively
ivan@154 41 % per iteration. Default value: mu = 0.
ivan@154 42 % cg_maxiter,... Maximum iterations for Conjugate-Gradients method.
ivan@154 43 % cg_tol Tolerance variable for Conjugate-Gradients method.
ivan@154 44 % =========================================================================
ivan@154 45 % OUTPUT ARGUMENTS:
ivan@154 46 % x_hat N x 1 recovered K-sparse vector.
ivan@154 47 % numiter Number of iterations executed.
ivan@154 48 % x_path Keeps a series of computed N x 1 K-sparse vectors
ivan@154 49 % until the end of the iterative process.
ivan@154 50 % =========================================================================
ivan@154 51 % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
ivan@154 52 % =========================================================================
ivan@154 53 % cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005.
ivan@154 54 % Email: jrom@acm.caltech.edu
ivan@154 55 % =========================================================================
ivan@154 56 % This work was supported in part by the European Commission under Grant
ivan@154 57 % MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like
ivan@154 58 % to acknowledge Rice University for his Faculty Fellowship.
ivan@154 59 % =========================================================================
ivan@154 60
ivan@154 61 [~,N] = size(Phi);
ivan@154 62
ivan@154 63 %% Initialize transpose of measurement matrix
ivan@154 64
ivan@154 65 Phi_t = Phi';
ivan@154 66
ivan@154 67 %% Initialize to zero vector
ivan@154 68 x_cur = zeros(N,1);
ivan@154 69 y_cur = zeros(N,1);
ivan@154 70 X_i = [];
ivan@154 71
ivan@154 72 x_path = zeros(N, params.ALPSiters);
ivan@154 73
ivan@154 74 %% CG params
ivan@154 75 if (params.solveNewtonx == 1 || params.solveNewtonb == 1)
ivan@154 76 cg_verbose = 0;
ivan@154 77 cg_A = Phi_t*Phi;
ivan@154 78 cg_b = Phi_t*y;
ivan@154 79 end;
ivan@154 80
ivan@154 81 %% Determine momentum step size selection strategy
ivan@154 82 optimizeTau = 0;
ivan@154 83 function_tau = 0;
ivan@154 84
ivan@154 85 if (isa(params.tau,'float'))
ivan@154 86 if (params.tau == -1)
ivan@154 87 optimizeTau = 1;
ivan@154 88 end;
ivan@154 89 elseif (isa(params.tau, 'function_handle'))
ivan@154 90 function_tau = 1;
ivan@154 91 end;
ivan@154 92
ivan@154 93 %% Determine step size selection strategy
ivan@154 94 function_mu = 0;
ivan@154 95 adaptive_mu = 0;
ivan@154 96
ivan@154 97 if (isa(params.mu,'float'))
ivan@154 98 function_mu = 0;
ivan@154 99 if (params.mu == 0)
ivan@154 100 adaptive_mu = 1;
ivan@154 101 else
ivan@154 102 adaptive_mu = 0;
ivan@154 103 end;
ivan@154 104 elseif (isa(params.mu,'function_handle'))
ivan@154 105 function_mu = 1;
ivan@154 106 end;
ivan@154 107
ivan@154 108 %% Help variables
ivan@154 109 complementary_Xi = ones(N,1);
ivan@154 110 setXi = zeros(N,1);
ivan@154 111 setYi = zeros(N,1);
ivan@154 112
ivan@154 113 i = 1;
ivan@154 114 %% infty-ALPS(#)
ivan@154 115 while (i <= params.ALPSiters)
ivan@154 116 x_path(:,i) = x_cur;
ivan@154 117 x_prev = x_cur;
ivan@154 118
ivan@154 119 % Compute the residual
ivan@154 120 if (i == 1)
ivan@154 121 res = y;
ivan@154 122 else
ivan@154 123 Phi_x_cur = Phi(:,X_i)*x_cur(X_i);
ivan@154 124 res = y - Phi_x_cur;
ivan@154 125 end;
ivan@154 126
ivan@154 127 % Compute the derivative
ivan@154 128 der = Phi_t*res;
ivan@154 129
ivan@154 130 % Determine S_i set via eq. (11)
ivan@154 131 complementary_Xi(X_i) = 0;
ivan@154 132 [~, ind_der] = sort(abs(der).*complementary_Xi, 'descend');
ivan@154 133 complementary_Xi(X_i) = 1;
ivan@154 134 S_i = [X_i; ind_der(1:K)];
ivan@154 135
ivan@154 136 ider = der(S_i);
ivan@154 137
ivan@154 138 setder = zeros(N,1);
ivan@154 139 setder(S_i) = 1;
ivan@154 140 if (adaptive_mu)
ivan@154 141 % Step size selection via eq. (12) and eq. (13)
ivan@154 142 Pder = Phi(:,S_i)*ider;
ivan@154 143 mu_bar = ider'*ider/(Pder'*Pder);
ivan@154 144 end;
ivan@154 145
ivan@154 146 iy_cur = y_cur.*setYi;
ivan@154 147 if (~function_tau) % If tau is not a function handle...
ivan@154 148 if (optimizeTau) % Compute optimized tau
ivan@154 149
ivan@154 150 % tau = argmin || u - Phi(x_i + y_i) ||
ivan@154 151 % = <Phi*y_i, u - Phi(x_i - mu/2 * grad_Si f(xi))>/||Phi*y_i||^2
ivan@154 152
ivan@154 153 if (i == 1)
ivan@154 154 params.tau = 0;
ivan@154 155 else
ivan@154 156 % u - Phi*(x_i - mu/2 grad_Si f(xi)) = u - Phi*b
ivan@154 157 if (adaptive_mu)
ivan@154 158 b = x_cur(S_i) + mu_bar*ider; % Non-zero elems: S_i
ivan@154 159 elseif (function_mu)
ivan@154 160 b = x_cur(S_i) + params.mu(i)*ider;
ivan@154 161 else b = x_cur(S_i) + params.mu*ider;
ivan@154 162 end;
ivan@154 163
ivan@154 164 y_Phi_b = y - Phi(:,S_i)*b;
ivan@154 165 Phi_y_prev = Phi(:,Y_i)*y_cur(Y_i); % Phi * y_i
ivan@154 166 params.tau = y_Phi_b'*Phi_y_prev/(Phi_y_prev'*Phi_y_prev);
ivan@154 167 end;
ivan@154 168
ivan@154 169 if (adaptive_mu)
ivan@154 170 y_cur = params.tau*iy_cur + mu_bar*der.*setder;
ivan@154 171 elseif (function_mu)
ivan@154 172 y_cur = params.tau*iy_cur + params.mu(i)*der.*setder;
ivan@154 173 else y_cur = params.tau*iy_cur + params.mu*der.*setder;
ivan@154 174 end;
ivan@154 175
ivan@154 176 Y_i = ne(y_cur,0);
ivan@154 177 setYi = zeros(N,1);
ivan@154 178 setYi(Y_i) = 1;
ivan@154 179 else % Tau fixed and user-defined
ivan@154 180 if (adaptive_mu)
ivan@154 181 y_cur = params.tau*iy_cur + mu_bar*der.*setder;
ivan@154 182 elseif (function_mu)
ivan@154 183 y_cur = params.tau*iy_cur + params.mu(i)*der.*setder;
ivan@154 184 else y_cur = params.tau*iy_cur + params.mu*der.*setder;
ivan@154 185 end;
ivan@154 186
ivan@154 187 Y_i = ne(y_cur,0);
ivan@154 188 setYi = zeros(N,1);
ivan@154 189 setYi(Y_i) = 1;
ivan@154 190 end;
ivan@154 191 else
ivan@154 192 if (adaptive_mu)
ivan@154 193 y_cur = params.tau(i)*iy_cur + mu_bar*der.*setder;
ivan@154 194 elseif (function_mu)
ivan@154 195 y_cur = params.tau(i)*iy_cur + params.mu(i)*der.*setder;
ivan@154 196 else y_cur = params.tau(i)*iy_cur + params.mu*der.*setder;
ivan@154 197 end;
ivan@154 198
ivan@154 199 Y_i = ne(y_cur,0);
ivan@154 200 setYi = zeros(N,1);
ivan@154 201 setYi(Y_i) = 1;
ivan@154 202 end;
ivan@154 203
ivan@154 204 % Hard-threshold b and compute X_{i+1}
ivan@154 205 set_Xi_Yi = setXi + setYi;
ivan@154 206 ind_Xi_Yi = find(set_Xi_Yi > 0);
ivan@154 207 z = x_cur(ind_Xi_Yi) + y_cur(ind_Xi_Yi);
ivan@154 208 [~, ind_z] = sort(abs(z), 'descend');
ivan@154 209 x_cur = zeros(N,1);
ivan@154 210 x_cur(ind_Xi_Yi(ind_z(1:K))) = z(ind_z(1:K));
ivan@154 211 X_i = ind_Xi_Yi(ind_z(1:K));
ivan@154 212
ivan@154 213 setXi = zeros(N,1);
ivan@154 214 setXi(X_i) = 1;
ivan@154 215
ivan@154 216 if (params.gradientDescentx == 1)
ivan@154 217 % Calculate gradient of estimated vector x_cur
ivan@154 218 Phi_x_cur = Phi(:,X_i)*x_cur(X_i);
ivan@154 219 res = y - Phi_x_cur;
ivan@154 220 der = Phi_t*res;
ivan@154 221
ivan@154 222 ider = der(X_i);
ivan@154 223
ivan@154 224 if (adaptive_mu)
ivan@154 225 Pder = Phi(:,X_i)*ider;
ivan@154 226 mu_bar = ider'*ider/(Pder'*Pder);
ivan@154 227 x_cur(X_i) = x_cur(X_i) + mu_bar*ider;
ivan@154 228 elseif (function_mu)
ivan@154 229 x_cur = x_cur(X_i) + params.mu(i)*ider;
ivan@154 230 else x_cur = x_cur(X_i) + params.mu*ider;
ivan@154 231 end;
ivan@154 232 elseif (params.solveNewtonx == 1)
ivan@154 233 % Similar to HTP
ivan@154 234 if (params.useCG == 1)
ivan@154 235 [v, ~, ~] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose);
ivan@154 236 else
ivan@154 237 v = cg_A(X_i,X_i)\cg_b(X_i);
ivan@154 238 end;
ivan@154 239 x_cur(X_i) = v;
ivan@154 240 end;
ivan@154 241
ivan@154 242 % Test stopping criterion
ivan@154 243 if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur))
ivan@154 244 break;
ivan@154 245 end;
ivan@154 246 i = i + 1;
ivan@154 247 end;
ivan@154 248
ivan@154 249 x_hat = x_cur;
ivan@154 250 numiter = i;
ivan@154 251
ivan@154 252 if (i > params.ALPSiters)
ivan@154 253 x_path = x_path(:,1:numiter-1);
ivan@154 254 else
ivan@154 255 x_path = x_path(:,1:numiter);
ivan@154 256 end;