ivan@154: function [x_hat, numiter, x_path] = one_ALPS(y, Phi, K, params) ivan@154: % ========================================================================= ivan@154: % 1-ALPS(#) algorithm - Beta Version ivan@154: % ========================================================================= ivan@154: % Algebraic Pursuit (ALPS) algorithm with 1-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,... If solveNewtonb == 1: Corresponds to solving a ivan@154: % Newton system restricted to a sparse support. ivan@154: % It is implemented via conjugate gradients. ivan@154: % If solveNewtonb == 0: Step size selection as described ivan@154: % in eqs. (12) and (13) in [1]. ivan@154: % Default value: solveNewtonb = 0. 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 = 0: momentum step size selection is ivan@154: % driven by the following formulas: ivan@154: % a_1 = 1; ivan@154: % a_{i+1} = (1+\sqrt(1+4a_i^2)/2; ivan@154: % tau = (a_i - 1)/(a_{i+1}); ivan@154: % described in [2] "A fast iterative ivan@154: % shrinkage-thresholding algorithm for linear ivan@154: % inverse problems", Beck A., and Teboulle M. 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: [M,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: Phi_x_cur = zeros(M,1); ivan@154: Y_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: fista = 0; ivan@154: optimizeTau = 0; ivan@154: a_prev = 1; ivan@154: function_tau = 0; ivan@154: ivan@154: if (isa(params.tau,'float')) ivan@154: function_tau = 0; ivan@154: if (params.tau == 0) ivan@154: fista = 1; ivan@154: optimizeTau = 0; ivan@154: elseif (params.tau == -1) ivan@154: optimizeTau = 1; ivan@154: fista = 0; 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_Yi = ones(N,1); ivan@154: ivan@154: i = 1; ivan@154: %% 1-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: % Compute the derivative ivan@154: der = Phi_t*res; ivan@154: else ivan@154: % Compute the derivative ivan@154: if (optimizeTau) ivan@154: res = y - Phi_x_cur - params.tau*Phi_diff; ivan@154: else ivan@154: res = y - Phi(:,Y_i)*y_cur(Y_i); ivan@154: end; ivan@154: der = Phi_t*res; ivan@154: end; ivan@154: ivan@154: Phi_x_prev = Phi_x_cur; ivan@154: ivan@154: % Determine S_i set via eq. (11) (change of variable from x_i to y_i) ivan@154: complementary_Yi(Y_i) = 0; ivan@154: [tmpArg, ind_der] = sort(abs(der).*complementary_Yi, 'descend'); ivan@154: complementary_Yi(Y_i) = 1; ivan@154: S_i = [Y_i; ind_der(1:K)]; ivan@154: ivan@154: ider = der(S_i); ivan@154: if (params.solveNewtonb == 1) ivan@154: % Compute least squares solution of the system A*y = (A*A)x using CG ivan@154: if (params.useCG == 1) ivan@154: [b, tmpArg, tmpArg] = cgsolve(cg_A(S_i, S_i), cg_b(S_i), params.cg_tol, params.cg_maxiter, cg_verbose); ivan@154: else ivan@154: b = cg_A(S_i,S_i)\cg_b(S_i); ivan@154: end; ivan@154: else ivan@154: % Step size selection via eq. (12) and eq. (13) (change of variable from x_i to y_i) ivan@154: if (adaptive_mu) ivan@154: Pder = Phi(:,S_i)*ider; ivan@154: mu_bar = ider'*ider/(Pder'*Pder); ivan@154: b = y_cur(S_i) + (mu_bar)*ider; ivan@154: elseif (function_mu) ivan@154: b = y_cur(S_i) + params.mu(i)*ider; ivan@154: else b = y_cur(S_i) + params.mu*ider; ivan@154: end; ivan@154: end; ivan@154: ivan@154: % Hard-threshold b and compute X_{i+1} ivan@154: [tmpArg, ind_b] = sort(abs(b), 'descend'); ivan@154: x_cur = zeros(N,1); ivan@154: x_cur(S_i(ind_b(1:K))) = b(ind_b(1:K)); ivan@154: X_i = S_i(ind_b(1:K)); 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: 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_i) = x_cur(X_i) + params.mu(i)*ider; ivan@154: else x_cur(X_i) = 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, tmpArg, tmpArg] = 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: if (~function_tau) % If tau is not a function handle... ivan@154: if (fista) % Fista configuration ivan@154: a_cur = (1 + sqrt(1 + 4*a_prev^2))/2; ivan@154: params.tau = (a_prev - 1)/a_cur; ivan@154: a_prev = a_cur; ivan@154: elseif (optimizeTau) % Compute optimized tau ivan@154: ivan@154: % tau = argmin ||u - Phi*y_{i+1}|| ivan@154: % = /||Phi*(x_cur - x_prev)||^2 ivan@154: ivan@154: Phi_x_cur = Phi(:,X_i)*x_cur(X_i); ivan@154: res = y - Phi_x_cur; ivan@154: if (i == 1) ivan@154: Phi_diff = Phi_x_cur; ivan@154: else ivan@154: Phi_diff = Phi_x_cur - Phi_x_prev; ivan@154: end; ivan@154: params.tau = res'*Phi_diff/(Phi_diff'*Phi_diff); ivan@154: end; ivan@154: ivan@154: y_cur = x_cur + params.tau*(x_cur - x_prev); ivan@154: Y_i = find(ne(y_cur, 0)); ivan@154: else ivan@154: y_cur = x_cur + params.tau(i)*(x_cur - x_prev); ivan@154: Y_i = find(ne(y_cur, 0)); 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;