bpdn_interior2.m

％ Atomizer Main Directory, Version .802 里面信号含有分解去噪合成过程的代码 ％---------------------------------------

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

% BPDN_Interior2 -- 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_Interior2( 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-3
%    OptTol	    Primal-dual gap tolerance for LP, default = 1e-3
%    CGAccuracy	    Accuracy of the CG solver,        default = 1e-3
%    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_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.
%    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_Interior2, BP_Simplex,
%    BPDN_Interior, dictionary
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 24 Feb 2001: First version of BPDN_Interior2.m.
%              Michael Saunders, SOL, Stanford University.
%              Adapted from Shaobing Chen's BPDN_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  PDSlambda

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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, BPMaxIter  =   30; end

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

  disp(' ')
  disp('------------------------------------------------------------------')
  disp(['BPDN_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) as in BPDN_Interior.m.
% Initialize by MOF representation.
% Make sure x and z are INTERIOR.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  disp(' ')
  disp('Initializing BPDN_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( 'BPDN_obj', 'BPDN_mat',  ...
	     b, mA, nA, options, x, y, z, xsize, zsize );

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