pdsco.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 1,246 行 · 第 1/4 页

function [x,y,z,inform,PDitns,CGitns,time] = ...
    pdsco( Fname,Aname,b,m,n,options,x0,y0,z0,xsize,zsize )

% pdsco.m -- Primal-Dual Barrier Method for Separable Convex Objectives
%--------------------------------------------------------------------------------
% [x,y,z,inform,PDitns,CGitns,time] = ...
%    pdsco( 'pdsobj', A      ,b,m,n,options,x0,y0,z0,xsize,zsize );
% [x,y,z,inform,PDitns,CGitns,time] = ...
%    pdsco( 'pdsobj','pdsmat',b,m,n,options,x0,y0,z0,xsize,zsize );
%
% pdsco (Primal-Dual barrier method, Separable Convex Objective)
% solves optimization problems of the form
%
%    minimize    phi(x)  +  1/2 norm(gamma*x)^2  +  1/2 norm(r/delta)^2
%      x,r
%    subject to  Ax + r = b,    x > 0,    r unconstrained.
%
% where
% phi(x) is a smooth separable convex function (possibly linear);
% A      is an m x n matrix, or an operator for forming A*x and A'*y;
% b      is a given m-vector;
% gamma  is a primal regularization parameter (typically small but may be 0);
% delta  is a dual regularization parameter (typically small or 1; must be >0);
% With positive gamma and delta, the primal-dual solution (x,y,z) is 
% bounded and unique.
%
% EXTERNAL FUNCTIONS:
% options         = pdscoSet;               provided with pdsco.m
% [obj,grad,hess] = pdsobj( x );            provided by user
%               y = pdsmat( mode,m,n,x );   provided by user if A isn't explicit
%
% OUTPUT ARGUMENTS:
% x         is the primal solution.
% y         is the dual solution associated with Ax + r = b.
% z         is the dual solution associated with x > 0.
%           x > 0, z > 0 is always true.
%           If (x,y,z,r) is optimal, r = delta^2*y and x.*z should be small.
% inform    = 0 if a solution is found;
%           = 1 if too many iterations were required;
%           = 2 if the linesearch failed too often.
% PDitns    is the number of Primal-Dual Barrier iterations required.
% CGitns    is the number of Conjugate-Gradient iterations required
%           if an iterative solver (LSQR) is used.
% time      is the cpu time used.
%
% INPUT ARGUMENTS:
% 'pdsobj'  defines phi(x) and its gradient and diagonal Hessian.
%           grad and hess are both n-vectors.
%           If phi(x) is the linear function c'x, pdsobj should return
%           [obj,grad,hess] = [c'x, c, zeros(n,1)].
% A         is an explicit m x n matrix (preferably sparse!).
% 'pdsmat'  is used if A is not known explicitly.
%           It returns y = A*x (mode=1) or y = A'*x (mode=2).
% b         is a given right-hand side vector.
% m, n      are the dimensions of A (or the A implied by pdsmat).
% options   is a structure that may be set and altered by pdscoSet.
% options.gamma   defines gamma.  Typical value: gamma = 1.0e-4.
% options.delta   defines delta.  Typical value: delta = 1.0e-4 or 1.0.
%                 Values smaller than 1.0e-6 or 1.0e-7 may affect numerical
%                 reliability.  Increasing delta usually improves convergence.
%                 For least-squares applications, delta = 1.0 is appropriate.
% help pdscoSet   describes the other options.
% x0, y0, z0      provide an initial solution.
% xsize, zsize    are estimates of the biggest x and z at the solution.
%                 They are used to scale (x,y,z).  Good estimates
%                 should improve the performance of the barrier method.
%----------------------------------------------------------------------


% AUTHOR:
%    Michael Saunders, Systems Optimization Laboratory (SOL),
%    Stanford University, Stanford, California, USA.
%    [e-mail]
%
%
% PRIVATE FUNCTIONS:              GLOBAL VARIABLES:
%    pdsxxxdistrib                   pdsAAA
%    pdsxxxlsqr                      pdsDDD1
%    pdsxxxlsqrmat                   pdsDDD2
%    pdsxxxmat
%
%
% NOTES:
% The matrix A is assumed to be reasonably well scaled: norm(A,inf) =~ 1.
% The vector b and objective phi(x) may be of any size, but be sure that
%    xsize and zsize have sensible values.
% The files defining pdsobj and pdsmat must not be called Fname.m or Aname.m !!!!
%
%
% DEVELOPMENT:
% 20 Jun 1997: Original version, derived from pdlp0.m.
% 19 Mar 1998: gamma, delta, wait added as parameters.
% 29 Mar 1998: LSQR's atol decreased dynamically from atol1 to atol2.
%              It looks like atol1 must be 1.0e-6 or less
%              (and atol2 = 1.0e-8 or less).
% 17 Mar 1998: LSproblem =  1 selects original LS problem for dy.
%              This is best if delta is small.
%              LSproblem =  2 converts it to a damped LS problem.
%              This allows lsqr to work with the smaller operator DA'.
%              It should be safe if delta = 1 say.
% 31 Mar 1998: Found that LSdamp > delta often accelerates convergence.
% 31 Mar 1998: Derived "primal" LS problem, which finds dx before dy.
%              LSproblem = 11 selects "primal" LS problem.
%              LSproblem = 12 converts it to a damped LS problem.
%              Implemented pdsobj, pdsmat, pdsxxx to allow for explicit A.
% 09 May 1998: 2x2 KKT may be symmetrized with X^(1/2).
%              LSproblem = 32 solves it via explicit K2 or SYMMLQ.
%              LSproblem = 33 equilibrates first n cols of K2.
%              pdsyyy.m applies operator K for SYMMLQ.
%              Cond(K2) is less than cond(A*D) -- a good sign.
% 09 May 1998: Initialized mu to be the proverbial mu0*(x'z)/n.
% 14 May 1998: Regard delta as a decreasing parameter like mu
%              to implement the conventional penalty function.
%              LSdamp = max( sqrt(mu), delta )  does the trick.
% 23 May 1998: x0,y0,z0 are now optional input parameters.  z is output.
% 03 Jun 1998: delta0 introduced to allow a choice.
%              delta0 = delta fixes delta throughout (this is the default).
%              delta0 = 1    allows delta to decrease with sqrt(mu).
%              Thus, max( sqrt(mu), delta )  <=  "delta"  <=  delta0.
% 13 Jun 1998: mulast = 0.1*opttol introduced.  No need to let mu get smaller.
% 11 Feb 2000: Added diagonal preconditioning for LSQR when A is explicit
%              (and LSproblem = 1, LSmethod = 3). Improved cond(DA') greatly.
%              Seems to halve the LSQR itns as hoped.
% 12 Feb 2000: Added analogous scaling to rows of A for SYMMLQ (LSproblem = 32)
% 14 Apr 2000: Initialize mu from initial [r; t; Xz], not just from gap.
% 15 Apr 2000: Use P = symmmd(ADDA) on first iteration for later chol(ADDA)
%              (if explicit A and LSmethod = 1).
% 16 Apr 2000: mulast = 0.01*opttol seems necessary on Xiaoming's LP problems.
% 28 Sep 2000: options implemented as a structure.
% 21 Nov 2000: Default x0, z0 scaled by norm(b).  (Altered 18 Jan 2001.)
% 24 Nov 2000: Iteration log prints scaled quantities Pinf, Dinf, Cinf.
% 14 Dec 2000: Removed methods that work with 2x2 and 3x3 systems and SYMMLQ.
% 18 Jan 2001: x, y, z, obj, grad scaled to allow for norm(b).
% 21 Jan 2001: Complementarity test is now  Cinf = maxXz <= opttol.
%              Should be ok because maxXz -> final mu = 0.1*opttol.
% 25 Jan 2001: x0, y0, z0, xsize, zsize now required input.
%              x, b are scaled by xsize.
%              y, z are scaled by zsize.
% 29 Jan 2001: Now that problem is scaled, use absolute value mufirst = mu0.
%              For warm starts, user should probably set mu0 = x0min * z0min
%              to get most variables centered.
% 31 Jan 2001: Set atol = max(Pinf, Dinf, Cinf) * 0.1  as in Inexact Newton.
%              atol2 no longer used.
% 07 Feb 2001: Satellite problem seems to need more accurate solves (.01).
% 09 Feb 2001: Test on r3ratio enforces Inexact Newton condition.
% 12 Feb 2001: lsqr is now private function pdsxxxlsqr.
%              Specialized to reduce atol and continue if necessary.
% 14 Feb 2001: stepx and stepz optimized to reduce the merit function f,
%              using Byunggyoo Kim's algebra.
% 13 Apr 2001: PDitns, CGitns now output parameters.
%-----------------------------------------------------------------------

global pdsAAA pdsDDD1 pdsDDD2 pdsDDD3;

fprintf('\n   pdsco.m                               Version of 13 Apr 2001')
fprintf('\n   Primal-dual barrier algorithm to minimize a separable convex')
fprintf('\n   objective function subject to constraints Ax + r = b, x >= 0\n')

operator =  isstr(Aname);
explicit = ~operator;

if explicit               % Aname is an explicit matrix A.
   [m,n]  =  size(Aname); % In case the user set (m,n) incorrectly.
   pdsAAA =  Aname;       % Explicit A = global pdsAAA.
   Aname  = 'pdsxxxmat';  % Use default Ax, A'y routine pdsxxxmat.m.
   nnzA   =  nnz(pdsAAA);
   if issparse(pdsAAA)
      fprintf('\nA is an explicit sparse matrix')
   else
      fprintf('\nA is an explicit dense matrix' )
   end
   fprintf('\nm        = %8g      n       = %8g       nnz(A)  =%9g', m,n,nnzA)
else
   fprintf('\nA is an operator')
   fprintf('\nm        = %8g      n       = %8g', m,n)
end

normb  = norm(b ,inf);   normx0 = max(x0);
normy0 = norm(y0,inf);   normz0 = max(z0);
fprintf('\nmax |b | = %8g      max x0  = %8g', normb , normx0)
fprintf(                 '       xsize   = %8g', xsize)
fprintf('\nmax |y0| = %8g      max z0  = %8g', normy0, normz0)
fprintf(                 '       zsize   = %8g', zsize)

%-----------------------------------------------------------------------
% Initialize.
%-----------------------------------------------------------------------
true   = 1;        em   = ones(m,1);
false  = 0;        en   = ones(n,1);
nb     = n + m;    nkkt = nb;
CGitns = 0;
inform = 0;

%-----------------------------------------------------------------------
% Grab input options.
%-----------------------------------------------------------------------
gamma     = options.gamma;
delta     = options.delta;
maxitn    = options.MaxIter;
featol    = options.FeaTol;
opttol    = options.OptTol;
steptol   = options.StepTol;
x0min     = options.x0min;
z0min     = options.z0min;
mu0       = options.mu0;
LSmethod  = options.LSmethod;   % 1=Cholesky    2=QR    3=LSQR
LSproblem = options.LSproblem;  % See below
itnlim    = options.LSQRMaxIter * min(m,n);
atol1     = options.LSQRatol1;  % Initial  atol
atol2     = options.LSQRatol2;  % Smallest atol (unless atol1 is smaller)
conlim    = options.LSQRconlim;
wait      = options.wait;

% LSproblem:        %  1 = dy          2 = dy shifted, DLS
                    % 11 = s          12 =  s shifted, DLS    (dx = Ds)
                    % 21 = dx
                    % 31 = 3x3 system, symmetrized by Z^{1/2}
                    % 32 = 2x2 system, symmetrized by X^{1/2}

%-----------------------------------------------------------------------
% Set other parameters.
%-----------------------------------------------------------------------
kminor    = 0;      % 1 stops after each iteration
delta0    = delta;  % Initial value of delta
eta       = 1e-4;   % Linesearch tolerance for "sufficient descent"
maxf      = 10;     % Linesearch backtrack limit (function evaluations)
maxfail   = 1;      % Linesearch failure limit (consecutive iterations)

%if delta >= 0.1,  LSproblem = 2; end

% Parameters for LSQR.
atolmin   = eps;    % Smallest atol if linesearch back-tracks
btol      = 0;      % Should be small (zero is ok)
show      = false;  % Controls lsqr iteration log

fprintf('\n\nx0min    = %8g      featol  = %8.1e', x0min, featol)
fprintf(                 '       gamma   = %8.1e', gamma)
fprintf(  '\nz0min    = %8g      opttol  = %8.1e', z0min, opttol)
fprintf(                 '       delta   = %8.1e', delta)
fprintf(  '\nmu0      = %8.1e      steptol = %8g', mu0  , steptol)
fprintf(                 '       delta0  = %8.1e', delta0)

fprintf('\n\nLSQR:')
fprintf('\natol1    = %8.1e      atol2   = %8.1e', atol1 , atol2 )
fprintf(                 '       btol    = %8.1e', btol )
fprintf('\nconlim   = %8.1e      itnlim  = %8g'  , conlim, itnlim)
fprintf(                 '       show    = %8g'  , show )

LSmethod  = 3;  %%% Hardwire LSQR
LSproblem = 1;  %%% and LS problem defining "dy".

if wait
    fprintf('\n\nReview parameters... then type "return"\n')
    keyboard
end

if eta < 0
    fprintf('\n\nLinesearch disabled by eta < 0')
end

%------------------------------------------------------------------------
% All parameters have now been set.
%------------------------------------------------------------------------
time    = cputime;
useChol = LSmethod == 1;
useQR   = LSmethod == 2;
direct  = LSmethod <= 2 & explicit;
solver  = '  LSQR';   if LSproblem >= 30, solver  = 'SYMMLQ'; end

%------------------------------------------------------------------------
% Scale the input data.
% The scaled variables are
%    xbar     = x/beta,
%    ybar     = y/zeta,
%    zbar     = z/zeta.
% Define
%    theta    = beta*zeta;
% The scaled function is
%    phibar   = ( 1   /theta) fbar(beta xbar),
%    gradient = (beta /theta) grad,
%    Hessian  = (beta2/theta) hess.
%------------------------------------------------------------------------
beta   = xsize;   if beta==0, beta = 1; end    % beta scales b, x.
zeta   = zsize;   if zeta==0, zeta = 1; end    % zeta scales y, z.
theta  = beta*zeta;                            % theta scales obj.
gamma  = gamma*beta/sqrt(theta);
delta  = delta*zeta/sqrt(theta);

beta2  = beta^2;    gamma2 = gamma^2;
delta2 = delta^2;   delta0 = delta;

b      = b /beta;   y0 = y0/zeta;
x0     = x0/beta;   z0 = z0/zeta;

%------------------------------------------------------------------------
% Initialize x, y, z, objective, etc.
%------------------------------------------------------------------------
x      = max( x0, x0min );   clear x0
z      = max( z0, z0min );   clear z0
y      = y0;                 clear y0

[obj,grad,hess] = feval( Fname, (x*beta) );
obj    = obj        /theta;               % Scaled obj.
grad   = grad*(beta /theta) + gamma2*x;   % grad includes x regularization.
Q      = hess*(beta2/theta) + gamma2;     % Q    includes x regularization.

%------------------------------------------------------------------------
% Initialize linear residuals  rlin =   b - A*x,
%                              tlin = - z - A'*y.
%------------------------------------------------------------------------
if explicit
   rlin =   b - pdsAAA *x;