gpsr_basic.m

这是一个压缩传感方面的Gradient Projection for Sparse Reconstruction 工具包。

if continuation&(cont_steps > 1)
   % If tau is scalar, first check top see if the first factor is 
   % too large (i.e., large enough to make the first 
   % solution all zeros). If so, make it a little smaller than that.
   % Also set to that value as default
   if prod(size(tau)) == 1
      if (firstTauFactorGiven == 0)|(firstTauFactor*tau >= max_tau)
         firstTauFactor = 0.8*max_tau / tau;
         fprintf(1,'parameter FirstTauFactor too large; changing')
      end
   end
   cont_factors = 10.^[log10(firstTauFactor):...
                    log10(1/firstTauFactor)/(cont_steps-1):0];
else
  cont_factors = 1;
  cont_steps = 1;
end
  

iter = 1;
if compute_mse
   mses(iter) = sum((x(:)-true(:)).^2);
end

% loop for continuation
for cont_loop = 1:cont_steps

    tau = final_tau * cont_factors(cont_loop);
    
    if verbose
        fprintf(1,'\nSetting tau = %8.4f\n',tau)
    end
    
    if cont_loop == cont_steps
       stopCriterion = final_stopCriterion;
       tolA = final_tolA;
    else 
       stopCriterion = 3;
       tolA = 1e-3;
    end
    
    % Compute and store initial value of the objective function
    resid =  y - A(x);
    f = 0.5*(resid(:)'*resid(:)) + ...
             sum(tau(:).*u(:)) + sum(tau(:).*v(:));

    objective(iter) = f;
    times(iter) = cputime - t0;
    
    % Compute the useful quantity resid_base
    resid_base = y - resid;

    % control variable for the outer loop and iteration counter
    % cont_outer = (norm(projected_gradient) > 1.e-5) ;

    keep_going = 1;

    if verbose
       fprintf(1,'\nInitial obj=%10.6e, nonzeros=%7d\n', f,num_nz_x);
    end

    while keep_going

      x_previous = x;

      % compute gradient
      temp = AT(resid_base);
      term  =  temp - Aty;
      gradu =  term + tau;
      gradv = -term + tau;

      % set search direction
      %du = -gradu; dv = -gradv; dx = du-dv; 
      dx = gradv-gradu;
      old_u = u; old_v = v;

      % calculate useful matrix-vector product involving dx
      auv = A(dx);
      dGd = auv(:)'*auv(:);

      % calculate unconstrained minimizer along this direction, use this
      % as the first guess of steplength parameter lambda
      %  lambda0 = - (gradu(:)'*du(:) + gradv(:)'*dv(:)) / dGd;

      % use instead a first guess based on the "conditional" direction
      condgradu = ((old_u>0) | (gradu<0)) .* gradu;
      condgradv = ((old_v>0) | (gradv<0)) .* gradv;
      auv_cond = A(condgradu-condgradv);
      dGd_cond = auv_cond(:)'*auv_cond(:);
      lambda0 = (gradu(:)'*condgradu(:) + gradv(:)'*condgradv(:)) /...
                (dGd_cond + realmin);

      % loop to determine steplength, starting wit the initial guess above.
      lambda = lambda0; 
      while 1
        % calculate step for this lambda and candidate point
        du = max(u-lambda*gradu,0.0) - u; 
        u_new = u + du;
        dv = max(v-lambda*gradv,0.0) - v; 
        v_new = v + dv;
        dx = du-dv; 
        x_new = x + dx;

        % evaluate function at the candidate point
        resid_base = A(x_new);
        resid = y - resid_base;
        f_new = 0.5*(resid(:)'*resid(:)) +  ...
        sum(tau(:).*u_new(:)) + sum(tau(:).*v_new(:));    
        % test sufficient decrease condition
        if f_new <= f + mu * (gradu'*du + gradv'*dv)
          %disp('OK')  
          break
        end
        lambda = lambda * lambda_backtrack;
        fprintf(1,'    reducing lambda to %6.2e\n', lambda)
      end
      u = u_new; 
      v = v_new; 
      prev_f = f; 
      f = f_new;
      uvmin = min(u,v); 
      u = u - uvmin; 
      v = v - uvmin; 
      x = u-v;

      % calculate nonzero pattern and number of nonzeros (do this *always*)
      nz_x_prev = nz_x;
      nz_x = (x~=0.0);
      num_nz_x = sum(nz_x(:));
      
      iter = iter + 1;
      objective(iter) = f;
      times(iter) = cputime-t0;
      lambdas(iter) = lambda;

      if compute_mse
        err = true - x;
        mses(iter) = (err(:)'*err(:));
      end
      
      % print out stuff
      if verbose
	fprintf(1,'It =%4d, obj=%9.5e, lambda=%6.2e, nz=%8d   ',...
	    iter, f, lambda, num_nz_x);
      end
      
    switch stopCriterion
	case 0,
      % compute the stopping criterion based on the change 
      % of the number of non-zero components of the estimate
      num_changes_active = (sum(nz_x(:)~=nz_x_prev(:)));
      if num_nz_x >= 1
         criterionActiveSet = num_changes_active;
      else
         criterionActiveSet = tolA / 2;
      end
      keep_going = (criterionActiveSet > tolA);
      if verbose
         fprintf(1,'Delta n-zeros = %d (target = %e)\n',...
                criterionActiveSet , tolA) 
      end
	case 1,
      % compute the stopping criterion based on the relative
      % variation of the objective function.
	  criterionObjective = abs(f-prev_f)/(prev_f);
      keep_going =  (criterionObjective > tolA);
      if verbose
         fprintf(1,'Delta obj. = %e (target = %e)\n',...
                criterionObjective , tolA) 
      end
	case 2,
       % stopping criterion based on relative norm of step taken
      delta_x_criterion = norm(dx(:))/norm(x(:));
      keep_going = (delta_x_criterion > tolA);
      if verbose
         fprintf(1,'Norm(delta x)/norm(x) = %e (target = %e)\n',...
             delta_x_criterion,tolA) 
      end
	case 3,
       % compute the "LCP" stopping criterion - again based on the previous
       % iterate. Make it "relative" to the norm of x.
        w = [ min(gradu(:), old_u(:)); min(gradv(:), old_v(:)) ];
        criterionLCP = norm(w(:), inf);
        criterionLCP = criterionLCP / ...
                 max([1.0e-6, norm(old_u(:),inf), norm(old_v(:),inf)]);
       keep_going = (criterionLCP > tolA);
       if verbose
         fprintf(1,'LCP = %e (target = %e)\n',criterionLCP,tolA) 
       end
    case 4,
      % continue if not yeat reached target value tolA
      keep_going = (f > tolA);
      if verbose
         fprintf(1,'Objective = %e (target = %e)\n',f,tolA) 
      end
	otherwise,
	  error(['Unknwon stopping criterion']);
    end % end of the stopping criteria switch
    
    % take no less than miniter... 
    if iter<=miniter
       keep_going = 1;
    else %and no more than maxiter iterations  
        if iter > maxiter
            keep_going = 0;
        end
    end

    end % end of the main loop of the GP algorithm

end % end of the continuation loop

% Print results
if verbose
  fprintf(1,'\nFinished the main algorithm!\nResults:\n')
  fprintf(1,'||A x - y ||_2^2 = %10.3e\n',resid(:)'*resid(:))
  fprintf(1,'||x||_1 = %10.3e\n',sum(abs(x(:))))
  fprintf(1,'Objective function = %10.3e\n',f);
  nz_x = (x~=0.0); num_nz_x = sum(nz_x(:));
  fprintf(1,'Number of non-zero components = %d\n',num_nz_x);
  fprintf(1,'CPU time so far = %10.3e\n', times(iter));
  fprintf(1,'\n');
end

% If the 'Debias' option is set to 1, we try to remove the bias from the l1
% penalty, by applying CG to the least-squares problem obtained by omitting
% the l1 term and fixing the zero coefficients at zero.

% do this only if the reduced linear least-squares problem is
% overdetermined, otherwise we are certainly applying CG to a problem with a
% singular Hessian

if debias
  
  if (num_nz_x > length(y(:)))
    if verbose
      fprintf(1,'\n')
      fprintf(1,'Debiasing requested, but not performed\n');
      fprintf(1,'There are too many nonzeros in x\n\n');
      fprintf(1,'nonzeros in x: %8d, length of y: %8d\n',...
	  num_nz_x, length(y(:)));
    end
  elseif (num_nz_x==0)
    if verbose
      fprintf(1,'\n')
      fprintf(1,'Debiasing requested, but not performed\n');
      fprintf(1,'x has no nonzeros\n\n');
    end
  else
    if verbose
      fprintf(1,'\n')
      fprintf(1,'Starting the debiasing phase...\n\n')
    end
    
    x_debias = x;
    zeroind = (x_debias~=0); 
    cont_debias_cg = 1;
    debias_start = iter;
    
    % calculate initial residual
    resid = A(x_debias);
    resid = resid-y;
    resid_prev = eps*ones(size(resid));
    
    rvec = AT(resid);
    
    % mask out the zeros
    rvec = rvec .* zeroind;
    rTr_cg = rvec(:)'*rvec(:);
    
    % set convergence threshold for the residual || RW x_debias - y ||_2
    tol_debias = tolD * (rvec(:)'*rvec(:));
    
    % initialize pvec
    pvec = -rvec;
    
    % main loop
    while cont_debias_cg
      
      % calculate A*p = Wt * Rt * R * W * pvec
      RWpvec = A(pvec);      
      Apvec = AT(RWpvec);
      
      % mask out the zero terms
      Apvec = Apvec .* zeroind;
      
      % calculate alpha for CG
      alpha_cg = rTr_cg / (pvec(:)'* Apvec(:));
      
      % take the step
      x_debias = x_debias + alpha_cg * pvec;
      resid = resid + alpha_cg * RWpvec;
      rvec  = rvec  + alpha_cg * Apvec;
      
      rTr_cg_plus = rvec(:)'*rvec(:);
      beta_cg = rTr_cg_plus / rTr_cg;
      pvec = -rvec + beta_cg * pvec;
      
      rTr_cg = rTr_cg_plus;
      
      iter = iter+1;
      
      objective(iter) = 0.5*(resid(:)'*resid(:)) + ...
	  sum(tau(:).*abs(x_debias(:)));
      times(iter) = cputime - t0;
      
      if compute_mse
	err = true - x_debias;
	mses(iter) = (err(:)'*err(:));
      end
      
      if verbose
	% in the debiasing CG phase, always use convergence criterion
	% based on the residual (this is standard for CG)
	fprintf(1,' Iter = %5d, debias resid = %13.8e, convergence = %8.3e\n', ...
	    iter, resid(:)'*resid(:), rTr_cg / tol_debias);
      end
      
      cont_debias_cg = (iter-debias_start <= miniter_debias )|...
	  ((rTr_cg > tol_debias) & ...
	  (iter-debias_start <= maxiter_debias));
      
    end
    
    if verbose
      fprintf(1,'\nFinished the debiasing phase!\nResults:\n')
      fprintf(1,'||A x - y ||^2_2 = %10.3e\n',resid(:)'*resid(:))
      fprintf(1,'||x||_1 = %10.3e\n',sum(abs(x(:))))
      fprintf(1,'Objective function = %10.3e\n',f);
      nz = (x_debias~=0.0);
      fprintf(1,'Number of non-zero components = %d\n',sum(nz(:)));
      fprintf(1,'CPU time so far = %10.3e\n', times(iter));
      fprintf(1,'\n');
    end
  end
  
  if compute_mse
    mses = mses/length(true(:));
  end
  
end