bp_interior.m

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

function c = BP_Interior( x, NameOfDict, par1, par2, par3, ...
                          FeaTol, PDGapTol, CGAccuracy, ...
                          gamma, delta, BPMaxIter )

% BP_Interior -- Basis Pursuit
%----------------------------------------------------------------------
%
%        Min || c ||_1  subject to \Phi * c = x
%
% via a Primal-Dual Logarithmic Barrier Interior Point method
% for an approximately equivalent perturbed linear program.
%
% Usage
%        c = BP_Interior( x, NameOfDict, par1, par2, par3
%                      [, FeaTol, PDGapTol, CGAccuracy, ...
%                         gamma, delta, BPMaxIter] );
%
% Inputs
%    x              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
%    delta          Dual   regularization parameter,  default = 1e-4
%    BPMaxIter      Maximum no. of barrier iterations,default = 30
%
%    Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
%
% Outputs
%    c              Coefficients of BP representation; column vector
%
% Description
%    1. Assumes noise level = 1.
%    2. BP_Interior initializes with the MOF solution and
%       iteratively improves the solution at each barrier
%       iteration by solving the central equations approximately
%       using the conjugate gradient method (CG).  It projects
%       the final decomposition into the space of feasible
%       reconstructions to achieve exact reconstruction.
%    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_Interior2,  BP_Simplex,
%    BPDN_Interior, BPDN_Interior2, dictionary
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 18 Mar 1997: Original version of BP_Interior.m.
%              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, delta      = 1e-4; end
  if nargin <11, BPMaxIter  =   30; end

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

  disp(' ')
  disp('------------------------------------------------------------------')
  disp(['BP_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 BP_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;
  x0  = 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('BP_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        = 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Feasibility:
%% Let dx = x - x0
%% Project dx onto the null space of \Phi by Conjugate Gradients:
%%     dx = (I - \Phi^T * (\Phi * \Phi^T)^-1 * \Phi) dx
%% x = x0 + dx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  dx     = x - x0;
  r      = FastSS( dx, n, NameOfDict, par1, par2, par3 );
  thresh = 1e-12 * norm(r);
  thresh = thresh^2;
  dd     = zeros(n,1);
  rho_1  = norm(r)^2;

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

  dx = dx - FastAA( dd, NameOfDict, par1, par2, par3 );

% dx = dx - FastAA(FastSS(dx, n, NameOfDict, par1, par2, par3), ...
%                  NameOfDict, par1, par2, par3) / (2 * L);

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