bpdn_interior.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 275 行

function [xrec, c] = BPDN_Interior( y, NameOfDict, par1, par2, par3, ...
				    FeaTol, PDGapTol, CGAccuracy, ...
				    gamma, BPMaxIter )

% BPDN_Interior -- Basis Pursuit De-Noising
%----------------------------------------------------------------------
%
%	Min .5 * || y - \Phi * c ||_2^2 + lambda * || c ||_1
%
% via a Primal-Dual Logarithmic Barrier Interior Point method
% for an equivalent perturbed linear program.
%
% Usage
%	 [xrec, c] = BPDN_Interior( y, NameOfDict, par1, par2, par3
%		                 [, FeaTol, PDGapTol, CGAccuracy, ...
%		                    gamma, BPMaxIter] );
%
% Inputs
%    y		    1-d signal; column vector
%    NameOfDict     String; name of the dictionary
%    par1,par2,par3 Parameters of the dictionary
%    FeaTol	    Feasibility tolerance for LP,     default = 1e-1
%    PDGapTol	    Primal-dual gap tolerance for LP, default = 1e-1
%    CGAccuracy	    Accuracy of the CG solver,        default = 1e-1
%    gamma          Primal regularization parameter,  default = 1e-4
%    BPMaxIter	    Maximum no. of barrier iterations,default = 30
%
%    Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
%
% Outputs
%    xrec           Reconstructed signal
%    c		    Coefficients of BP representation; column vector
%
% Description
%    1. Assumes noise level == 1
%    2. BPDN_Interior initializes with the MOF solution and iteratively
%	improves the solution at each barrier iteration by solving the 
%	central equations approximately using CG.
%    3. The barrier iterations stop if
%	    Primal Infeasibility < FeaTol
%	    Dual Infeasibility   < FeaTol
%	    Primal-Dual Gap      < PDGapTol
%	or if
%	    No. of Barrier Iterations > BPMaxIter
%	or if
%	    No. of CG Iterations > CGMaxIter
%    4. For FeaTol, PDGapTol, CGAccuracy, we recommend
%           1e-1 (the default values) for routine processing,
%           1e-2 or 1e-3              for super-resolution.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% References
%    S. Chen
%    Basis Pursuit,
%    Ph.D. Thesis, Department of Statistics, Stanford Univ., 1995.
%
%    Gill et al. 
%    Solving reduced KKT systems in barrier methods for
%    linear and quadratic programming,
%    Report SOL 91-7, Dept of Operations Research,
%    Stanford University, July 1991.
%
%    Lustig et al.
%    Interior point methods for linear programming: 
%    Computational state of the art,
%    ORSA J. on Computing, 6(1), Winter 1994.
%
%    S. S. Chen, D. L. Donoho and M. A. Saunders,
%    Atomic decomposition by basis pursuit,
%    SIAM J. Scientific Computing, 20(1), pp. 33--61 (1998).
%    Revised in     SIAM Review, 43(1), pp. 129--159 (2001).
%
% See also
%    BP_Interior,    BP_Interior2, BP_Simplex,
%    BPDN_Interior2, dictionary
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 09 Jan 1996: Original version of BPDN_Interior.
%              Shaobing Chen, Dept of Statistics, Stanford University.
% 09 Apr 2001: Edited by Michael Saunders, SOL, Stanford University.
%              Note that the 1e-1 default tolerances are
%              EXTREMELY LOOSE.
%----------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set any missing parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if nargin < 6, FeaTol     = 1e-1; end
  if nargin < 7, PDGapTol   = 1e-1; end
  if nargin < 8, CGAccuracy = 1e-1; end
  if nargin < 9, gamma      = 1e-4; end
  if nargin <10, BPMaxIter  =   30; end

%%%%%%%%%%%%%%%%%%%%%%%%%
% Set size of problem.
%%%%%%%%%%%%%%%%%%%%%%%%%
  b         = y(:);
  n         = length(b);
  [m L]     = SizeOfDict( n, NameOfDict, par1, par2, par3 );
  CGMaxIter = m * log2(m);
  lambda    = (2 * log(m)) ^ .5;
  delta     = 1;
  gamma2    = gamma^2;
  delta2    = delta^2;

  disp(' ')
  disp('------------------------------------------------------------------')
  disp(['BPDN_Interior:         NameOfDict = ' NameOfDict])
  fprintf('   par1(1) =%8.1g      par2(1) =%8.1g      par3(1) =%8.1g\n', ...
          par1(1), par2(1), par3(1))
  disp('------------------------------------------------------------------')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize (x, y, z), making sure x and z are INTERIOR.
% Initialize by MOF representation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  disp(' ')
  disp('Initializing BPDN_Interior by MOF ...')

  x   = MOF( b, NameOfDict, par1, par2, par3, 1e-5 );
  x   = [ x .* (x > 0); (-x) .* (x < 0) ];
  y   = FastSS( sign(x), n, NameOfDict, par1, par2, par3 );
  Aty = FastAA(       y,    NameOfDict, par1, par2, par3 );
  Aty = Aty(1:m);
  mp  = 1.1 * max(abs(Aty));
  Aty = Aty ./ mp;
  y   =   y ./ mp;
  z1  = lambda - Aty;
  z2  = lambda + Aty;
  z   = [z1; z2];
  x   = x + .1;
  z   = z + gamma2*x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the barrier parameter mu.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  mu0 = .01;
  Ax  = FastSS( x, n, NameOfDict, par1, par2, par3 );
  Aty = FastAA( y,    NameOfDict, par1, par2, par3 );
  f   = [z .* x
         b      - Ax  - delta2*y
         lambda - Aty + gamma2*x - z];
  mu  = mu0 * norm(f) / sqrt(2 * m);


%%%%%%%%%%%%%%%%%%%%%%
% Barrier iterations.
%%%%%%%%%%%%%%%%%%%%%%
  time  = cputime;   cgitns     = 0;
  k     = 0;         PrimalStep = 0;
  sigma = 0.99;

  disp(' ')
  disp('BPDN_Interior Iterations ...')
  disp(' ')
  disp('                        FeaTol          PDGapTol           CGAccuracy')
  fprintf('                      %8.1e          %8.1e             %8.1e\n', ...
          FeaTol, PDGapTol, CGAccuracy)
  disp(' ')
  disp('Itn    mu    PrimStep  PrimFea  DualFea   PDGap        Obj      CGItn')

  for i = 0:BPMaxIter
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate the residuals.
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    Ax  = FastSS( x, n, NameOfDict, par1, par2, par3 );
    Aty = FastAA( y,    NameOfDict, par1, par2, par3 );
    v1  = mu - z .* x;
    v2  = b      - Ax  - delta2*y;
    v3  = lambda - Aty + gamma2*x - z;

    %%%%%%%%%%%%%%%%%%%%%%%%%
    % Test convergence.
    % Printout iteration log.
    %%%%%%%%%%%%%%%%%%%%%%%%%
    PrimalFeas = norm(v2) / (1 + norm(x));
    DualFeas   = norm(v3) / (1 + norm(y));
    PDGap      = z' * x   / (1 + abs(sum(x)));
    Obj        = 0.5 * norm(v2+y)^2 + lambda * sum(abs(x));

    str1       = sprintf( '%3g %8.1e',    i, mu );
    str2       = sprintf( ' %8.1e %8.1e', PrimalStep, PrimalFeas);
    str3       = sprintf( ' %8.1e %8.1e', DualFeas, PDGap );
    str4       = sprintf( ' %14.7e %5g',  Obj, k );
    disp( [str1 str2 str3 str4] )

    Feasible   = PrimalFeas < FeaTol  &  DualFeas < FeaTol;
    Converged  =      PDGap < PDGapTol;
    if Feasible & Converged,
      fprintf('     Converged\n')
      break;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Solve Newton's Direction by CG: 
    %    (A diag(H) A' + delta2 I) * dy = v2 - A*H*(v1./x - v3)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    H      = 1 ./ (z ./ x + gamma2);
    w      = H .* (v1./ x - v3);
    q      = FastSS( w, n, NameOfDict, par1, par2, par3 );
    r      = v2 - q;
    thresh = CGAccuracy * norm(r);
    thresh = thresh^2;
    dy     = zeros(n,1);
    rho_1  = norm(r)^2;

    for k = 1:CGMaxIter
      %fprintf('	CG Step: %3g:	     %g\n', k, rho_1);
      cgitns = cgitns + 1;
      if k == 1,
	 p    = r;
      else
	 beta = rho_1 / rho_0;
	 p    = r + beta * p;
      end
      Atp   = FastAA(      p,    NameOfDict, par1, par2, par3 );
      q     = FastSS( H.*Atp, n, NameOfDict, par1, par2, par3 );
      w     = q + delta2 * p;
      alpha = rho_1 / (p' * w);
      r     = r - alpha * w;
      rho_0 = rho_1;
      rho_1 = norm(r)^2;
      dy    = dy + alpha * p;
      if rho_1 < thresh,
  	 break;
      end
    end

    if k == CGMaxIter,
      x = xold;
      fprintf('     Too many CG iterations\n')
      break;
    end

    Aty = FastAA( dy, NameOfDict, par1, par2, par3 );
    dx  = (Aty + v1./x - v3) .* H;
    dz  = (v1  -    z .* dx) ./ x;

    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate maximum steps.
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    stepx = 1e20;  stepz = 1e20;
    blocking = find( dx < 0 );
    if length( blocking ) > 0
      stepx = min( x(blocking) ./ (- dx(blocking)) );
    end

    blocking = find( dz < 0 );
    if length( blocking ) > 0
      stepz = min( z(blocking) ./ (- dz(blocking)) );
    end

    stepx = min( sigma * stepx, 1 );
    stepz = min( sigma * stepz, 1 );

    %%%%%%%%%%%%%%%%%%%%
    % Update (x y z mu).
    %%%%%%%%%%%%%%%%%%%%
    xold  = x;
    x     = x + stepx * dx; 
    y     = y + stepz * dy; 
    z     = z + stepz * dz;
    mu    = (1 -  min([stepx stepz sigma])) * mu;
    PrimalStep = norm(dx) * stepx / norm(xold);
  end

  c    = x(1:m) - x((m+1):(2*m));
  xrec = FastS( c, n, NameOfDict, par1, par2, par3 );
  time = cputime - time;
  fprintf('\nPDitns =%4g    CGitns =%7g', i, cgitns);
  fprintf('    time =%10.1f    sum(c) =%12.5e\n', time, sum(abs(c)));