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