ivan@154: function [x_hat, numiter, x_path] = l_ALPS(y, Phi, K, params) ivan@154: % ========================================================================= ivan@154: % l-ALPS(#) algorithm - Beta Version ivan@154: % ========================================================================= ivan@154: % Algebraic Pursuit (ALPS) algorithm with l-memory acceleration. ivan@154: % ivan@154: % Detailed discussion on the algorithm can be found in section III 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 = -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: y_cur = zeros(N,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: 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_Yi = ones(N,1); ivan@154: x = zeros(N,params.memory + 1); ivan@154: Phi_x = zeros(M,params.memory + 1); ivan@154: tau_candidates = zeros(params.memory,1); ivan@154: y_candidates = zeros(N,params.memory); ivan@154: residual_energy = zeros(params.memory,1); ivan@154: ivan@154: i = 1; ivan@154: %% l-ALPS(#) ivan@154: while (i <= params.ALPSiters) ivan@154: x_path(:,i) = x(:,end); ivan@154: ivan@154: for k = 1:params.memory ivan@154: x(:,k) = x(:,k+1); ivan@154: Phi_x(:,k) = Phi_x(:,k+1); ivan@154: end; ivan@154: ivan@154: % Compute the residual ivan@154: if (i == 1) ivan@154: res = y; ivan@154: else ivan@154: Phi_y_cur = Phi(:,Y_i)*y_cur(Y_i); ivan@154: res = y - Phi_y_cur; ivan@154: end; ivan@154: ivan@154: % Compute the derivative ivan@154: der = Phi_t*res; ivan@154: ivan@154: % Determine S_i 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(:,end) = zeros(N,1); ivan@154: x(S_i(ind_b(1:K)),end) = 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 current x vector ivan@154: Phi_x_cur = Phi(:,X_i)*x(X_i,end); 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(X_i,end) = x(X_i,end) + mu_bar*ider; ivan@154: elseif (function_mu) ivan@154: x(X_i,end) = x(X_i,end) + params.mu(i)*ider; ivan@154: else x(X_i,end) = x(X_i,end) + 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(X_i,end) = v; ivan@154: end; ivan@154: ivan@154: Phi_x(:,end) = Phi(:,X_i)*x(X_i,end); ivan@154: res = y - Phi_x(:,end); ivan@154: if (~function_tau) % If tau is not a function handle... ivan@154: if (optimizeTau) % Compute optimized tau ivan@154: if (i > params.memory) ivan@154: ivan@154: % tau = argmin ||u - Phi*y_{i+1}|| ivan@154: % = /||Phi*(x_cur - x_prev)||^2 ivan@154: ivan@154: for k = 1:params.memory ivan@154: Phi_diff = Phi_x(:,end) - Phi_x(:,k); ivan@154: tau_candidates(k) = res'*Phi_diff/(Phi_diff'*Phi_diff); ivan@154: y_candidates(:,k) = x(:,end) + tau_candidates(k)*(x(:,end) - x(:,k)); ivan@154: residual_energy(k) = norm(y - Phi*y_candidates(:,k)); ivan@154: end; ivan@154: ivan@154: [tmpArg,min_ind] = min(residual_energy); ivan@154: y_cur = y_candidates(:,min_ind); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: else ivan@154: Phi_diff = Phi_x(:,end) - Phi_x(:,end-1); ivan@154: params.tau = res'*Phi_diff/(Phi_diff'*Phi_diff); ivan@154: y_cur = x(:,end) + params.tau*(x(:,end) - x(:,end-1)); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: end; ivan@154: else ivan@154: if (i > params.memory) ivan@154: for k = 1:params.memory ivan@154: y_candidates(:,k) = x(:,end) + params.tau*(x(:,end) - x(:,k)); ivan@154: residual_energy(k) = norm(y - Phi*y_candidates(:,k)); ivan@154: end; ivan@154: ivan@154: [tmpArg,min_ind] = min(residual_energy); ivan@154: y_cur = y_candidates(:,min_ind); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: else ivan@154: y_cur = x(:,end) + params.tau*(x(:,end) - x(:,end-1)); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: end; ivan@154: end; ivan@154: else ivan@154: if (i > params.memory) ivan@154: for k = 1:params.memory ivan@154: y_candidates(:,k) = x(:,end) + params.tau(i)*(x(:,end) - x(:,k)); ivan@154: residual_energy(k) = norm(y - Phi*y_candidates(:,k))^2; ivan@154: end; ivan@154: ivan@154: [tmpArg,min_ind] = min(residual_energy); ivan@154: y_cur = y_candidates(:,min_ind); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: else ivan@154: y_cur = x(:,end) + params.tau(i)*(x(:,end) - x(:,end-1)); ivan@154: Y_i = find(ne(y_cur,0)); ivan@154: end; ivan@154: end; ivan@154: ivan@154: % Test stopping criterion ivan@154: if (i > 1) && (norm(x(:,end) - x(:,end-1)) < params.tol*norm(x(:,end))) ivan@154: break; ivan@154: end; ivan@154: i = i + 1; ivan@154: end; ivan@154: ivan@154: x_hat = x(:,end); 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;