bp_interior2.m

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

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

% BP_Interior2 -- Basis Pursuit
%----------------------------------------------------------------------
%
%	 Min || c ||_1  subject to \Phi * c = x
%
% via a Primal-Dual Logarithmic Barrier Interior method
% for an approximately equivalent perturbed linear program.
%
% Usage
%	 c = BP_Interior2( x, NameOfDict, par1, par2, par3 ...
%	                [, FeaTol, OptTol, CGAccuracy, ...
%                          gamma, delta, 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-3
%    OptTol	    Optimality  tolerance for LP,     default = 1e-3
%    CGAccuracy     Accuracy of the CG solver,        default = 1e-3
%    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_Interior2 initializes with the MOF solution and
%       iteratively improves the solution at each barrier
%       iteration by solving the central equations approximately
%       using LSQR, a conjugate-gradient (CG) solver.
%       FeaTol and delta affect how accurately \Phi c = x
%       is satisfied.  Smaller values increase the accuracy.
%    3. The barrier iterations stop if
%	    Primal Infeasibility < FeaTol
%	    Dual Infeasibility   < FeaTol
%           Primal-Dual Gap      < OptTol
%       or if
%	    No. of Barrier Iterations > BPMaxIter
%       or if
%	    No. of CG Iterations      > CGMaxIter
%    4. For FeaTol, OptTol, CGAccuracy, we recommend
%           1e-3 (the default values) for routine processing,
%           1e-4 or 1e-5              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_Simplex,
%    BPDN_Interior, BPDN_Interior2, dictionary
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 24 Feb 2001: First version of BP_Interior2.m.
%              Michael Saunders, SOL, Stanford University.
%              Adapted from Shaobing Chen's BP_Interior.m.
%              Uses pdsco.m to solve for c.
% 05 Apr 2001: Use z0min = 1 instead of 0.1
%              (because we expect most z >> 0).
% 12 Apr 2001: x   = x + 0.01*norm(x,inf);
%------------------------------------------------------------------------

  global  PDSn  PDSmat  PDSpar1  PDSpar2  PDSpar3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set any missing parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if nargin < 6, FeaTol     = 1e-3; end
  if nargin < 7, OptTol     = 1e-3; end
  if nargin < 8, CGAccuracy = 1e-3; 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_Interior2:          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.
% This helps get xsize and zsize roughly right for pdsco.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  disp(' ')
  disp('Initializing BP_Interior2 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];
  z   = z + gamma2*x;
% x   = x + 0.1;               % This is scale dependent.
  x   = x + 0.01*norm(x,inf);  % This is similar to x0min later, but
                               % is sensible for BP because A*1 = 0.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set global variables for pdsco.m.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  PDSn      = n;
  PDSmat    = NameOfDict;
  PDSpar1   = par1;
  PDSpar2   = par2;
  PDSpar3   = par3;
  PDSlambda = lambda;

% Set parameters for pdsco.m.

  mA        = n;
  nA        = 2*m;
  CGMaxIter = CGMaxIter / mA;  % pdsco scales this up by mA;

  options = pdscoSet;
  options = pdscoSet( options, ...
                     'gamma      ', gamma     , ...
                     'delta      ', delta     , ...
                     'MaxIter    ', BPMaxIter , ...
                     'FeaTol     ', FeaTol    , ...
                     'OptTol     ', OptTol    , ...
                     'StepTol    ', 0.99      , ...
                     'x0min      ', 0.1       , ...
                     'z0min      ', 1.0       , ...
                     'mu0        ', 0.01      , ...
                     'LSmethod   ', 3         , ...
                     'LSproblem  ', 1         , ...
                     'LSQRMaxIter', CGMaxIter , ...
                     'LSQRatol1  ', CGAccuracy, ...
                     'LSQRatol2  ', 1e-15     , ... % No longer used
                     'wait       ', 0    );

  bsize = norm(b,inf);
  xsize = max(x);
  xsize = max(bsize, xsize);   % 11 Apr 2001: Bigger is better.
  zsize = max(z);

  [x,y,z,inform,PDitns,CGitns,time] = ...
      pdsco( 'BP_obj', 'BP_mat',      ...
	     b, mA, nA, options, x, y, z, xsize, zsize );

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