pdsco.m

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

function [x,y,z,inform] = pdsco( Fname, Aname, b, m, n, options, x0, y0, z0 )
%
% 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),
% and A is an m x n matrix.  The regularization parameters (gamma, delta)
% are typically positive and small.  They ensure that the primal and dual
% solutions (x,y) are bounded and unique.
%
% FUNCTIONS:
% options        = pdscoSet;
% [x,y,z,inform] = pdsco( 'pdsobj',  A      , b, m, n, options );
% [x,y,z,inform] = pdsco( 'pdsobj', 'pdsmat', b, m, n, options );
% [x,y,z,inform] = pdsco( 'pdsobj', 'pdsmat', b, m, n, options, x0, y0, z0 );
%[obj,grad,hess] = pdsobj( x )
%              v = pdsmat( mode, m, n, x )
%
% 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
% 'pdsobj'  defines phi(x) and its gradient and diagonal Hessian.
%           grad and hess are both vectors.
% A         is an explicit m x n matrix (preferably sparse!).
% 'pdsmat'  is used if A is not known explicitly.
%           It returns v = A*x (mode=1) and v = A'*x (mode=2).
% b         is the 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   may be zero or positive.  A typical value is gamma = 1.0e-4.
% options.delta   must be positive.  Again, a typical value is delta = 1.0e-4.
%                 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  (if present) provide an initial solution, probably from a
%           previous call.

%-----------------------------------------------------------------------
% EXTERNAL FUNCTIONS NEEDED:
% pdscoSet.m, pdsxxx.m, lsqr.m,   (provided with pdsco.m)
% pdsmat.m,   pdsobj.m            (provided by user -- whatever names you like)
%
% If phi(x) is the linear function c'x, pdsobj should return
% [obj,grad,hess] = [c'x, c, zeros(n,1)].
%
% NOTE:
% A, b, x are assumed to be reasonably well scaled.
% The files corresponding to pdsobj and pdsmat must not be called
% Fname.m or Aname.m !!!!

% 20 Jun 1997: Original version, derived from pdlp0.m.
%              Author: Michael Saunders, Stanford University
%              ([e-mail]).
%
% 05 Mar 1997: Found that low accuracy solves lead to failure unless
%              dz = xinv.*(v - z.*dx) is replaced by dz = t + gamma2*dx - Ady,
%              as recommended in Saunders & Tomlin, Report SOL 96-3, 1996.
%              E.g. SYMMLQ with rtol = 1.0e-10 solved quite well,
%                          but  rtol = 1.0e-9 or bigger let t diverge.
%              With new dz,     rtol = 1.0e-7 seems accurate enough,
%                                      1.0e-6 is border-line.
%
% 13 Jun 1997: LSQR installed to solve for dy.
% 18 Mar 1998: Can no longer find much difference between the above dz options.
%              Will keep using dz = t + gamma2*dx - Ady (just in case).
% 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 implements as a structure.
%-----------------------------------------------------------------------

global AAA DDD1 DDD2 DDD3;

fprintf('\n    pdsco.m                            Version of 29 Sep 2000')
fprintf('\n    Primal-dual barrier algorithm --- for optimization with a')
fprintf('\n    separable convex obj and linear constraints Ax = 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.
    AAA    =  Aname;       % Explicit A = global AAA.
    Aname  = 'pdsmat';     % Use default Ax, A'y routine pdsmat.m.
    nnzA   = nnz(AAA);
    if nnzA == m*n
       fprintf('\nA is an explicit dense matrix')
    else
       fprintf('\nA is an explicit sparse 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

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

%-----------------------------------------------------------------------
% Grab input options.
%-----------------------------------------------------------------------
gamma  = options.gamma;
delta  = options.delta;
maxitn = options.MaxIter;
featol = options.FeaTol;
opttol = options.OptTol;
steptol= options.StepTol;
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;

if nargin <= 6
    x0 = en;
    z0 = en;
    y0 = zeros(m,1);
end

%-----------------------------------------------------------------------
% Set other parameters.
%-----------------------------------------------------------------------
kminor    = 0;      % 1 stops after each iteration
x0min     = 0.1;    % Minimum size of initial x if x0 is provided
z0min     = 0.1;    % Minimum size of initial z if z0 is provided
mu0       = 0.1;    % Scales initial mu
delta0    = delta;  % Initial value of delta
eta       = 1e-4;   % Linesearch tolerance for "sufficient descent"
maxf      = 10;     % Linesearch backtrack limit (function evaluations)
maxfail   = 4;      % Linesearch failure limit (consecutive iterations)

LSproblem = 1;      %  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}
LSmethod  = 3;      %  1 = Cholesky    2 = QR        3 = LSQR or SYMMLQ

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

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

% Parameters for symmlq.
precon  = false;         shift  = 0;
showsym = false;         check  = false;
itnsym  = 15*(m+n);   %  rtol   = atol;   % Use lsqr's atol below


fprintf('\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\nLSproblem=%3g', LSproblem)
fprintf(               '  ( 1 = dy, LS      2 = dybar, DLS)')
fprintf('\n               (11 = s , LS     12 = sbar , DLS)     (dx = Ds)')
fprintf('\n               (21 = dx, LS)')
fprintf('\n               (31 = K3 sym-izd 32 = K2 symmetrized by X^(1/2))')

fprintf(  '\nLSmethod =%3g', LSmethod)
fprintf(               '  ( 1 = Cholesky    2 = QR     3 = LSQR or SYMMLQ)')

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 )

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

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

tic;

%------------------------------------------------------------------------
% Initialize x, z, y, objective function, etc.
%------------------------------------------------------------------------
normb   = norm(b,inf);
gamma2  = gamma^2;
delta2  = delta^2;
x       = max( x0, x0min );
z       = max( z0, z0min );
y       = y0;
useChol = LSmethod == 1;
useQR   = LSmethod == 2;
direct  = LSmethod <= 2 & explicit;
if LSproblem < 30
    solver  = '  lsqr';
else
    solver  = 'symmlq';
end

[obj,grad,hess] = feval( Fname, x );

%------------------------------------------------------------------------
% Initialize linear residuals  rlin =   b - A*x,
%                              tlin = - z - A'*y.
%------------------------------------------------------------------------
if explicit
    r = AAA*x;
else
    r = feval( Aname, 1, m, n, x );
end
t = zeros( n,1 );
if norm(y,inf) > 0
    if explicit, t = AAA'*y;
    else         t = feval( Aname, 2, m, n, y );  end
end
rlin    =   b - r;
tlin    = - z - t;
Xz      =   x.*z;

%------------------------------------------------------------------------
% 09 May 1998: Initialize mu to a fraction of x'z/n (as others do).
% 14 May 1998: Make delta decrease with mu.
% 03 Jun 1998: Keep delta <= delta0.
% 13 Jun 1998: Keep mu >= mulast.
% 14 Apr 2000: mufirst = mu0 * (x'z)/n is too small if initial r, t are big.
%              Revert to balancing two parts of rhs of KKT system:
%              (  b - Ax   )     ( 0  )     (  r  )
%              (g - A'y - z)  =  ( 0  )  +  (  t  ).
%              (mu e  -  Xz)     (mu e)     (- Xz )
%------------------------------------------------------------------------
r       = rlin  -   delta2 * y;
t       = tlin  +  (gamma2 * x  +  grad);
f       = [ r; t; Xz ];
normf   = norm(f);

mufirst = mu0 * normf / n;     %Replaces  mufirst = mu0 * (sum(Xz) / n);
mulast  = 0.01 * opttol;
mu      = mufirst;
LSdamp  = max( sqrt(mu), delta );
LSdamp  = min( LSdamp  , delta0);
LSdamp2 = LSdamp^2;

%------------------------------------------------------------------------
% Initialize nonlinear residuals  r = rlin - delta^2*y,
%                                 t = tlin + gamma^2*x + grad,
%                                 v = mu*e - X*z.
%------------------------------------------------------------------------