gpsr_bb.m

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

% store given stopping criterion and threshold, because we're going 
% to change them in the continuation procedure
final_stopCriterion = stopCriterion;
final_tolA = tolA;

% set continuation factors
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 = %0.5g\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 initial gradient and the useful 
    % quantity resid_base
    resid_base = y - resid;

    alpha = 1.0;
    alphas(iter) = alpha;

    % control variable for the outer loop and iteration counter

    keep_going = 1;

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

    while keep_going

      % compute gradient
      temp = AT(resid_base);

      term  =  temp - Aty;
      gradu =  term + tau;
      gradv = -term + tau;

      % projection and computation of search direction vector
      du = max(u - alpha*gradu, 0.0) - u;
      dv = max(v - alpha*gradv, 0.0) - v;
      dx = du-dv;
      old_u = u; 
      old_v = v;

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

      if (enforceMonotone==1)
        % monotone variant: calculate minimizer along the direction (du,dv)
        lambda0 = - (gradu(:)'*du(:) + gradv(:)'*dv(:))/(realmin+dGd);
        if lambda0 < 0
          fprintf(' ERROR: lambda0 = %10.3e negative. Quit\n', lambda0);
          return;
        end
        lambda = min(lambda0,1);
      else
        %nonmonotone variant: choose lambda=1
        lambda = 1;
      end

      u = old_u + lambda * du;
      v = old_v + lambda * dv;
      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(:));
      
      % update residual and function
      resid = y - resid_base - lambda*auv;
      prev_f = f;
      f = 0.5*(resid(:)'*resid(:)) +  sum(tau(:).*u(:)) + ...
          sum(tau(:).*v(:));

      % compute new alpha
      dd  = du(:)'*du(:) + dv(:)'*dv(:);  
      if dGd <= 0
        % something wrong if we get to here
        fprintf(1,' dGd=%12.4e, nonpositive curvature detected\n', dGd);
        alpha = alphamax;
      else
        alpha = min(alphamax,max(alphamin,dd/dGd));
      end
      resid_base = resid_base + lambda*auv; 
      
      % print out stuff
      if verbose
         fprintf(1,'It=%4d, obj=%9.5e, alpha=%6.2e, nz=%8d  ',...
             iter, f, alpha, num_nz_x);
      end

      iter = iter + 1;
      objective(iter) = f;
      times(iter) = cputime-t0;
      alphas(iter) = alpha;

      if compute_mse
        err = true - x;
        mses(iter) = (err(:)'*err(:));
      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(['Unknown 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 BB-QP 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
      
      % in the debiasing CG phase, always use convergence criterion
      % based on the residual (this is standard for CG)
      if verbose
	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