ivan@154: function [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params) ivan@154: % ========================================================================= ivan@154: % infty-ALPS(#) algorithm - Beta Version ivan@154: % ========================================================================= ivan@154: % Algebraic Pursuit (ALPS) algorithm with infty-memory acceleration. ivan@154: % ivan@154: % Detailed discussion on the algorithm can be found in ivan@154: % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written ivan@154: % by Volkan Cevher, Technical Report, 2011. ivan@154: % ========================================================================= ivan@154: % INPUT ARGUMENTS: ivan@154: % y M x 1 undersampled measurement vector. ivan@154: % Phi M x N regression matrix. ivan@154: % K Sparsity of underlying vector x* or desired ivan@154: % sparsity of solution. ivan@154: % params Structure of parameters. These are: ivan@154: % ivan@154: % tol,... Early stopping tolerance. Default value: tol = ivan@154: % 1-e5 ivan@154: % ALPSiters,... Maximum number of algorithm iterations. Default ivan@154: % value: 300. ivan@154: % solveNewtonb,... Value: solveNewtonb = 0. Not used in infty ivan@154: % methods. ivan@154: % gradientDescentx,... If gradientDescentx == 1: single gradient ivan@154: % update of x_{i+1} restricted ot its support with ivan@154: % line search. Default value: gradientDescentx = ivan@154: % 1. ivan@154: % solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit ivan@154: % (c.f. Simon Foucart, "Hard Thresholding Pursuit," ivan@154: % preprint, 2010). Default vale: solveNewtonx = 0 ivan@154: % tau,... Variable that controls the momentum in ivan@154: % non-memoryless case. Ignored in memoryless ivan@154: % case. Default value: tau = 1/2. ivan@154: % Special cases: ivan@154: % - tau = -1: momentum step size is automatically ivan@154: % optimized in every step. ivan@154: % - tau as a function handle: user defined ivan@154: % behavior of tau momentum term. ivan@154: % mu,... Variable that controls the step size selection. ivan@154: % When mu = 0, step size is computed adaptively ivan@154: % per iteration. Default value: mu = 0. ivan@154: % cg_maxiter,... Maximum iterations for Conjugate-Gradients method. ivan@154: % cg_tol Tolerance variable for Conjugate-Gradients method. ivan@154: % ========================================================================= ivan@154: % OUTPUT ARGUMENTS: ivan@154: % x_hat N x 1 recovered K-sparse vector. ivan@154: % numiter Number of iterations executed. ivan@154: % x_path Keeps a series of computed N x 1 K-sparse vectors ivan@154: % until the end of the iterative process. ivan@154: % ========================================================================= ivan@154: % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL. ivan@154: % ========================================================================= ivan@154: % cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005. ivan@154: % Email: jrom@acm.caltech.edu ivan@154: % ========================================================================= ivan@154: % This work was supported in part by the European Commission under Grant ivan@154: % MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like ivan@154: % to acknowledge Rice University for his Faculty Fellowship. ivan@154: % ========================================================================= ivan@154: ivan@154: [~,N] = size(Phi); ivan@154: ivan@154: %% Initialize transpose of measurement matrix ivan@154: ivan@154: Phi_t = Phi'; ivan@154: ivan@154: %% Initialize to zero vector ivan@154: x_cur = zeros(N,1); ivan@154: y_cur = zeros(N,1); ivan@154: X_i = []; ivan@154: ivan@154: x_path = zeros(N, params.ALPSiters); ivan@154: ivan@154: %% CG params ivan@154: if (params.solveNewtonx == 1 || params.solveNewtonb == 1) ivan@154: cg_verbose = 0; ivan@154: cg_A = Phi_t*Phi; ivan@154: cg_b = Phi_t*y; ivan@154: end; ivan@154: ivan@154: %% Determine momentum step size selection strategy ivan@154: optimizeTau = 0; ivan@154: function_tau = 0; ivan@154: ivan@154: if (isa(params.tau,'float')) ivan@154: if (params.tau == -1) ivan@154: optimizeTau = 1; ivan@154: end; ivan@154: elseif (isa(params.tau, 'function_handle')) ivan@154: function_tau = 1; ivan@154: end; ivan@154: ivan@154: %% Determine step size selection strategy ivan@154: function_mu = 0; ivan@154: adaptive_mu = 0; ivan@154: ivan@154: if (isa(params.mu,'float')) ivan@154: function_mu = 0; ivan@154: if (params.mu == 0) ivan@154: adaptive_mu = 1; ivan@154: else ivan@154: adaptive_mu = 0; ivan@154: end; ivan@154: elseif (isa(params.mu,'function_handle')) ivan@154: function_mu = 1; ivan@154: end; ivan@154: ivan@154: %% Help variables ivan@154: complementary_Xi = ones(N,1); ivan@154: setXi = zeros(N,1); ivan@154: setYi = zeros(N,1); ivan@154: ivan@154: i = 1; ivan@154: %% infty-ALPS(#) ivan@154: while (i <= params.ALPSiters) ivan@154: x_path(:,i) = x_cur; ivan@154: x_prev = x_cur; ivan@154: ivan@154: % Compute the residual ivan@154: if (i == 1) ivan@154: res = y; ivan@154: else ivan@154: Phi_x_cur = Phi(:,X_i)*x_cur(X_i); ivan@154: res = y - Phi_x_cur; ivan@154: end; ivan@154: ivan@154: % Compute the derivative ivan@154: der = Phi_t*res; ivan@154: ivan@154: % Determine S_i set via eq. (11) ivan@154: complementary_Xi(X_i) = 0; ivan@154: [~, ind_der] = sort(abs(der).*complementary_Xi, 'descend'); ivan@154: complementary_Xi(X_i) = 1; ivan@154: S_i = [X_i; ind_der(1:K)]; ivan@154: ivan@154: ider = der(S_i); ivan@154: ivan@154: setder = zeros(N,1); ivan@154: setder(S_i) = 1; ivan@154: if (adaptive_mu) ivan@154: % Step size selection via eq. (12) and eq. (13) ivan@154: Pder = Phi(:,S_i)*ider; ivan@154: mu_bar = ider'*ider/(Pder'*Pder); ivan@154: end; ivan@154: ivan@154: iy_cur = y_cur.*setYi; ivan@154: if (~function_tau) % If tau is not a function handle... ivan@154: if (optimizeTau) % Compute optimized tau ivan@154: ivan@154: % tau = argmin || u - Phi(x_i + y_i) || ivan@154: % = /||Phi*y_i||^2 ivan@154: ivan@154: if (i == 1) ivan@154: params.tau = 0; ivan@154: else ivan@154: % u - Phi*(x_i - mu/2 grad_Si f(xi)) = u - Phi*b ivan@154: if (adaptive_mu) ivan@154: b = x_cur(S_i) + mu_bar*ider; % Non-zero elems: S_i ivan@154: elseif (function_mu) ivan@154: b = x_cur(S_i) + params.mu(i)*ider; ivan@154: else b = x_cur(S_i) + params.mu*ider; ivan@154: end; ivan@154: ivan@154: y_Phi_b = y - Phi(:,S_i)*b; ivan@154: Phi_y_prev = Phi(:,Y_i)*y_cur(Y_i); % Phi * y_i ivan@154: params.tau = y_Phi_b'*Phi_y_prev/(Phi_y_prev'*Phi_y_prev); ivan@154: end; ivan@154: ivan@154: if (adaptive_mu) ivan@154: y_cur = params.tau*iy_cur + mu_bar*der.*setder; ivan@154: elseif (function_mu) ivan@154: y_cur = params.tau*iy_cur + params.mu(i)*der.*setder; ivan@154: else y_cur = params.tau*iy_cur + params.mu*der.*setder; ivan@154: end; ivan@154: ivan@154: Y_i = ne(y_cur,0); ivan@154: setYi = zeros(N,1); ivan@154: setYi(Y_i) = 1; ivan@154: else % Tau fixed and user-defined ivan@154: if (adaptive_mu) ivan@154: y_cur = params.tau*iy_cur + mu_bar*der.*setder; ivan@154: elseif (function_mu) ivan@154: y_cur = params.tau*iy_cur + params.mu(i)*der.*setder; ivan@154: else y_cur = params.tau*iy_cur + params.mu*der.*setder; ivan@154: end; ivan@154: ivan@154: Y_i = ne(y_cur,0); ivan@154: setYi = zeros(N,1); ivan@154: setYi(Y_i) = 1; ivan@154: end; ivan@154: else ivan@154: if (adaptive_mu) ivan@154: y_cur = params.tau(i)*iy_cur + mu_bar*der.*setder; ivan@154: elseif (function_mu) ivan@154: y_cur = params.tau(i)*iy_cur + params.mu(i)*der.*setder; ivan@154: else y_cur = params.tau(i)*iy_cur + params.mu*der.*setder; ivan@154: end; ivan@154: ivan@154: Y_i = ne(y_cur,0); ivan@154: setYi = zeros(N,1); ivan@154: setYi(Y_i) = 1; ivan@154: end; ivan@154: ivan@154: % Hard-threshold b and compute X_{i+1} ivan@154: set_Xi_Yi = setXi + setYi; ivan@154: ind_Xi_Yi = find(set_Xi_Yi > 0); ivan@154: z = x_cur(ind_Xi_Yi) + y_cur(ind_Xi_Yi); ivan@154: [~, ind_z] = sort(abs(z), 'descend'); ivan@154: x_cur = zeros(N,1); ivan@154: x_cur(ind_Xi_Yi(ind_z(1:K))) = z(ind_z(1:K)); ivan@154: X_i = ind_Xi_Yi(ind_z(1:K)); ivan@154: ivan@154: setXi = zeros(N,1); ivan@154: setXi(X_i) = 1; ivan@154: ivan@154: if (params.gradientDescentx == 1) ivan@154: % Calculate gradient of estimated vector x_cur ivan@154: Phi_x_cur = Phi(:,X_i)*x_cur(X_i); ivan@154: res = y - Phi_x_cur; ivan@154: der = Phi_t*res; ivan@154: ivan@154: ider = der(X_i); ivan@154: ivan@154: if (adaptive_mu) ivan@154: Pder = Phi(:,X_i)*ider; ivan@154: mu_bar = ider'*ider/(Pder'*Pder); ivan@154: x_cur(X_i) = x_cur(X_i) + mu_bar*ider; ivan@154: elseif (function_mu) ivan@154: x_cur = x_cur(X_i) + params.mu(i)*ider; ivan@154: else x_cur = x_cur(X_i) + params.mu*ider; ivan@154: end; ivan@154: elseif (params.solveNewtonx == 1) ivan@154: % Similar to HTP ivan@154: if (params.useCG == 1) ivan@154: [v, ~, ~] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose); ivan@154: else ivan@154: v = cg_A(X_i,X_i)\cg_b(X_i); ivan@154: end; ivan@154: x_cur(X_i) = v; ivan@154: end; ivan@154: ivan@154: % Test stopping criterion ivan@154: if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur)) ivan@154: break; ivan@154: end; ivan@154: i = i + 1; ivan@154: end; ivan@154: ivan@154: x_hat = x_cur; ivan@154: numiter = i; ivan@154: ivan@154: if (i > params.ALPSiters) ivan@154: x_path = x_path(:,1:numiter-1); ivan@154: else ivan@154: x_path = x_path(:,1:numiter); ivan@154: end;