nikcleju@17: function [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
nikcleju@17: 
nikcleju@17: %%
nikcleju@17: % [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
nikcleju@17: %
nikcleju@17: % Greedy Analysis Pursuit Algorithm
nikcleju@17: % This aims to find an approximate (sometimes exact) solution of
nikcleju@17: %    xhat = argmin || Omega * x ||_0   subject to   || y - M * x ||_2 <= epsilon.
nikcleju@17: %
nikcleju@17: % Outputs:
nikcleju@17: %   xhat : estimate of the target cosparse vector x0.
nikcleju@17: %   Lambdahat : estimate of the cosupport of x0.
nikcleju@17: %
nikcleju@17: % Inputs:
nikcleju@17: %   y : observation/measurement vector of a target cosparse solution x0,
nikcleju@17: %       given by relation  y = M * x0 + noise.
nikcleju@17: %   M : measurement matrix. This should be given either as a matrix or as a function handle
nikcleju@17: %       which implements linear transformation.
nikcleju@17: %   MH : conjugate transpose of M. 
nikcleju@17: %   Omega : analysis operator. Like M, this should be given either as a matrix or as a function
nikcleju@17: %           handle which implements linear transformation.
nikcleju@17: %   OmegaH : conjugate transpose of OmegaH.
nikcleju@17: %   params : parameters that govern the behavior of the algorithm (mostly).
nikcleju@17: %      params.num_iteration : GAP performs this number of iterations.
nikcleju@17: %      params.greedy_level : determines how many rows of Omega GAP eliminates at each iteration.
nikcleju@17: %                            if the value is < 1, then the rows to be eliminated are determined by
nikcleju@17: %                                j : |omega_j * xhat| > greedy_level * max_i |omega_i * xhat|.
nikcleju@17: %                            if the value is >= 1, then greedy_level is the number of rows to be
nikcleju@17: %                            eliminated at each iteration.
nikcleju@17: %      params.stopping_coefficient_size : when the maximum analysis coefficient is smaller than
nikcleju@17: %                                         this, GAP terminates.
nikcleju@17: %      params.l2solver : legitimate values are 'pseudoinverse' or 'cg'. determines which method
nikcleju@17: %                        is used to compute
nikcleju@17: %                        argmin || Omega_Lambdahat * x ||_2   subject to  || y - M * x ||_2 <= epsilon.
nikcleju@17: %      params.l2_accuracy : when l2solver is 'cg', this determines how accurately the above 
nikcleju@17: %                           problem is solved.
nikcleju@17: %      params.noise_level : this corresponds to epsilon above.
nikcleju@17: %   xinit : initial estimate of x0 that GAP will start with. can be zeros(d, 1).
nikcleju@17: %
nikcleju@17: % Examples:
nikcleju@17: %
nikcleju@17: % Not particularly interesting:
nikcleju@17: % >> d = 100; p = 110; m = 60; 
nikcleju@17: % >> M = randn(m, d);
nikcleju@17: % >> Omega = randn(p, d);
nikcleju@17: % >> y = M * x0 + noise;
nikcleju@17: % >> params.num_iteration = 40;
nikcleju@17: % >> params.greedy_level = 0.9;
nikcleju@17: % >> params.stopping_coefficient_size = 1e-4;
nikcleju@17: % >> params.l2solver = 'pseudoinverse';
nikcleju@17: % >> [xhat, Lambdahat] = GAP(y, M, M', Omega, Omega', params, zeros(d, 1));
nikcleju@17: %
nikcleju@17: % Assuming that FourierSampling.m, FourierSamplingH.m, FDAnalysis.m, etc. exist:
nikcleju@17: % >> n = 128;
nikcleju@17: % >> M = @(t) FourierSampling(t, n);
nikcleju@17: % >> MH = @(u) FourierSamplingH(u, n);
nikcleju@17: % >> Omega = @(t) FDAnalysis(t, n);
nikcleju@17: % >> OmegaH = @(u) FDSynthesis(t, n);
nikcleju@17: % >> params.num_iteration = 1000;
nikcleju@17: % >> params.greedy_level = 50;
nikcleju@17: % >> params.stopping_coefficient_size = 1e-5;
nikcleju@17: % >> params.l2solver = 'cg';   % in fact, 'pseudoinverse' does not even make sense.
nikcleju@17: % >> [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, zeros(d, 1));
nikcleju@17: %
nikcleju@17: % Above: FourierSampling and FourierSamplingH are conjugate transpose of each other.
nikcleju@17: %        FDAnalysis and FDSynthesis are conjugate transpose of each other.
nikcleju@17: %        These routines are problem specific and need to be implemented by the user.
nikcleju@17: 
nikcleju@17: d = length(xinit(:));
nikcleju@17: 
nikcleju@17: if strcmp(class(Omega), 'function_handle')
nikcleju@17:     p = length(Omega(zeros(d,1)));
nikcleju@17: else    %% Omega is a matrix
nikcleju@17:     p = size(Omega, 1);
nikcleju@17: end
nikcleju@17: 
nikcleju@17: iter = 0;
nikcleju@17: lagmult = 1e-4;
nikcleju@17: Lambdahat = 1:p;
nikcleju@17: while iter < params.num_iteration
nikcleju@17:     iter = iter + 1;
nikcleju@17:     [xhat, analysis_repr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, lagmult, params);
nikcleju@17:     [to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, params.greedy_level);
nikcleju@17:     %disp(['** maxcoef=', num2str(maxcoef), ' target=', num2str(params.stopping_coefficient_size), ' rows_remaining=', num2str(length(Lambdahat)), ' lagmult=', num2str(lagmult)]);
nikcleju@17:     if check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params)
nikcleju@17:         break;
nikcleju@17:     end
nikcleju@17:     xinit = xhat;
nikcleju@17:     Lambdahat(to_be_removed) = [];
nikcleju@17: 
nikcleju@17:     %n = sqrt(d);
nikcleju@17:     %figure(9);
nikcleju@17:     %RR = zeros(2*n, n-1);
nikcleju@17:     %RR(Lambdahat) = 1;
nikcleju@17:     %XD = ones(n, n);
nikcleju@17:     %XD(:, 2:end) = XD(:, 2:end) .* RR(1:n, :);
nikcleju@17:     %XD(:, 1:(end-1)) = XD(:, 1:(end-1)) .* RR(1:n, :);
nikcleju@17:     %XD(2:end, :) = XD(2:end, :) .* RR((n+1):(2*n), :)';
nikcleju@17:     %XD(1:(end-1), :) = XD(1:(end-1), :) .* RR((n+1):(2*n), :)';
nikcleju@17:     %XD = FD2DiagnosisPlot(n, Lambdahat);
nikcleju@17:     %imshow(XD);
nikcleju@17:     %figure(10);
nikcleju@17:     %imshow(reshape(real(xhat), n, n));
nikcleju@17: end
nikcleju@17: return;
nikcleju@17: 
nikcleju@17: 
nikcleju@17: function [to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, greedy_level)
nikcleju@17: 
nikcleju@17:     abscoef = abs(analysis_repr(:));
nikcleju@17:     n = length(abscoef);
nikcleju@17:     maxcoef = max(abscoef);
nikcleju@17:     if greedy_level >= 1
nikcleju@17:         qq = quantile(abscoef, 1-greedy_level/n);
nikcleju@17:     else
nikcleju@17:         qq = maxcoef*greedy_level;
nikcleju@17:     end
nikcleju@17: 
nikcleju@17:     to_be_removed = find(abscoef >= qq);
nikcleju@17:     return;
nikcleju@17: 
nikcleju@17: function r = check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params)
nikcleju@17: 
nikcleju@17:     r = 0;
nikcleju@17: 
nikcleju@17:     if isfield(params, 'stopping_coefficient_size') && maxcoef < params.stopping_coefficient_size
nikcleju@17:         r = 1;
nikcleju@17:         return;
nikcleju@17:     end
nikcleju@17: 
nikcleju@17:     if isfield(params, 'stopping_lagrange_multiplier_size') && lagmult > params.stopping_lagrange_multiplier_size
nikcleju@17:         r = 1;
nikcleju@17:         return;
nikcleju@17:     end
nikcleju@17: 
nikcleju@17:     if isfield(params, 'stopping_relative_solution_change') && norm(xhat-xinit)/norm(xhat) < params.stopping_relative_solution_change
nikcleju@17:         r = 1;
nikcleju@17:         return;
nikcleju@17:     end
nikcleju@17: 
nikcleju@17:     if isfield(params, 'stopping_cosparsity') && length(Lambdahat) < params.stopping_cosparsity
nikcleju@17:         r = 1;
nikcleju@17:         return;
nikcleju@17:     end