pdsco.m

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

   tlin = - z - pdsAAA'*y;
else
   rlin =   b - feval( Aname, 1, m, n, x );
   tlin = - z - feval( Aname, 2, m, n, y );
end

%------------------------------------------------------------------------
% Initialize mu
% and the nonlinear residuals  r    = rlin - delta^2*y,
%                              t    = tlin + gamma^2*x + grad,
%                              v    = mu*e - 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 )
%
% 25 Jan 2001: Now that b and obj are scaled (and hence x,y,z),
%              we should be able to use an absolute value: mufirst = mu0;
%              But 0.1 worked poorly on StarTest1 (with x0min = z0min = 0.1).
%
% 29 Jan 2001: We might as well use mu0 = x0min * z0min;
%              so that most variables are centered after a warm start.
%------------------------------------------------------------------------
r       = rlin - delta2*y;
t       = tlin + grad;           % grad includes gamma2*x
Xz      = x.*z;

%f      = norm([r; t; Xz]);      % Original norm.
%f      = norm([r; t; Xz],1);    % Probably should have been that.
%mufirst= mu0 * f / n;           % Replaces  mufirst = mu0 * (sum(Xz) / n);

mufirst = mu0;                   % Absolute value.
mulast  = 0.1 * opttol;
mulast  = min( mulast, mufirst );
mu      = mufirst;
LSdamp  = max( sqrt(mu), delta );
LSdamp  = min( LSdamp  , delta0);
LSdamp2 = LSdamp^2;
r       = rlin - LSdamp2*y;      % 25 Jan 2001: Should have been there earlier.
v       = mu   - Xz;
f       = norm([r; t; v]);       % f = merit function for the linesearch.

% Initialize other things.

PDitns    = 0;
converged = 0;
atol      = atol1;
atol2     = max( atol2, atolmin );

%  Iteration log.

stepx   = 0;
stepz   = 0;
nf      = 0;
itncg   = 0;
nfail   = 0;

head1   = '\n\nItn   mu stepx stepz  Pinf  Dinf';
head2   = '  Cinf   Objective    nf  center';
if direct, head3 = ''; else head3 =['  atol ' solver ' Inexact']; end
fprintf([ head1 head2 head3 ])

mininf  = 1e-99;
regterm = gamma2 * (x'*x)  +  LSdamp2 * (y'*y);
objreg  = obj  +  0.5 * regterm;
objtrue = objreg * theta;

maxXz   = max(Xz);
minXz   = min(Xz);
center  = maxXz / minXz;

Pinf    = norm(r,inf);   Pinf = max( Pinf, mininf );
Dinf    = norm(t,inf);   Dinf = max( Dinf, mininf );
Cinf    = maxXz;         Cinf = max( Cinf, mininf );

fprintf('\n%3g                 ', PDitns       )
fprintf('%6.1f%6.1f' , log10(Pinf), log10(Dinf))
fprintf('%6.1f%15.7e', log10(Cinf), objtrue    )
fprintf('   %8.1f'   , center                  )

if kminor
   fprintf('\n\nStart of first minor itn...\n')
   keyboard
end

%-----------------------------------------------------------------------
% Main loop.
%-----------------------------------------------------------------------

while ~converged
   PDitns = PDitns + 1;
%  x1 = x;  y1 = y;  z1 = z;  % Save for debugging at end of iteration.
%  t1 = t;  r1 = r;  v1 = v;

% 31 Jan 2001: Set atol according to actual progress, a la Inexact Newton.
% 07 Feb 2001: 0.1 not small enough for Satellite problem.  Try 0.01.

   r3norm = max([Pinf  Dinf  Cinf]);
   atol   = min([atol r3norm*0.01]);
   atol   = max([atol  atolmin   ]);

%-----------------------------------------------------------------------
%  Define a damped Newton iteration for solving
%     ( r, t, v ) = 0,
%  keeping  x, z > 0.  We eliminate dz
%  to obtain the system
%
%     [-H   A'] [ dx ] = [ w ],    d2I = delta^2 I,   w = t - v./x;
%     [ A  d2I] [ dy ] = [ r ]
%
%  which is equivalent to
%
%     [-dI  DA'] [ s  ] = [  Dw   ],    dI = delta I,
%     [ AD  dI ] [ dy ] = [r/delta]     D  = H^{-1/2),  dx = delta D s,
%
%  and also equivalent to the least-squares problem
%
%     min || [ DA']dy  -  [  Dw   ] ||.                                (*)
%         || [ dI ]       [r/delta] ||
%
% 17 Mar 1998: We solve the latter as the damped least-squares problem
%
%     min || [ DA']dybar  -  [D wbar] ||,  wbar = w - (A'r)/delta^2,   (**)
%         || [ dI ]          [  0   ] ||     dy = dybar + r/delta^2,
%
%    to allow lsqr to work with the smaller operator DA'.
%
% 30 Mar 1998: LSproblem = 1 or 2 selects (*) or (**) respectively.
%              (*)  seems safer if delta is small, but
%              (**) is more efficient (and safe) if delta = 1, say.
%
% 31 Mar 1998: LSproblem = 11 or 12 selects alternative LS problem
%
%     min || [ AD ]s  -  [    r    ] ||,     dx = Ds,                 (***)
%         || [ dI ]      [-delta Dw] ||      dy = (r - Adx) / delta^2
%
% and associated damped LS problem:
%
%     min || [ AD ]sbar  -  [ rbar ] ||,      s = sbar - Dw,         (****)
%         || [ dI ]         [  0   ] ||    rbar = r + AD^2 w.
%
% 06 Apr 1998: LSproblem = 21 selects equivalent LS problem
%
%     min || [    A     ]dx  -  [    r    ] ||,                     (*****)
%         || [delta Dinv]       [-delta Dw] ||     dy = (r - Adx) / delta^2
%-----------------------------------------------------------------------
   H       = Q  +  z./x;
   w       = t  -  v./x;
   Hinv    = 1 ./ H;
   D       = sqrt(Hinv);
   pdsDDD1 = D;
   rw      = [explicit  LSproblem  LSmethod  LSdamp m n 0];  % Passed to LSQR.

   if LSproblem == 1
      %-----------------------------------------------------------------
      %  Solve (*) for dy.
      %-----------------------------------------------------------------
      Dw   = D.*w;

      if direct
         DD   = sparse( 1:n, 1:n, D, n, n );
         AD   = AAA*DD;

         if useChol
            d2I  = LSdamp2 * em;
            d2I  = sparse( 1:m, 1:m, d2I, m, m );
            ADDA = AD * AD'  +  d2I;

            if PDitns==1,  P = symmmd(ADDA);  end  % Do ordering only once.

            [R,indef] = chol(ADDA(P,P));
            if indef
               fprintf('\n\n   chol says AD^2A'' is not positive definite')
               break
            end

            rhs  = r  +  AAA * (Hinv.*w);
%           dy   = ADDA \ rhs;
            dy   = R \ (R' \ rhs(P));   dy(P) = dy;
         
	 else % useQR
            dI   = LSdamp * em;
            rhs  = [ Dw; r/LSdamp ];
            dy   = [ AD'; diag(dI) ] \ rhs;
         end

         Ady  = AAA'*dy;
         dx   = Hinv .* (Ady - w);
         Adx  = AAA*dx;


      else % Iterative solve using LSQR.
         rhs     = [ Dw; r/LSdamp ];
         damp    = 0;

         if explicit            % A is a sparse matrix.
	    precon     = true;
	    if precon    % Construct diagonal preconditioner for LSQR
	       DD      = sparse( 1:n, 1:n, D, n, n );
	       AD      = pdsAAA*DD;
	       AD2     = AD.^2;
	       wD      = sum( (AD2') )';   % Sum of squares of each row of AD
	       wD      = sqrt( wD + LSdamp2 );
	       pdsDDD2 = 1 ./ wD;
	    end
         else                   % A is an operator
	    precon    = false;  
	 end

	 rw(7)        = precon;
         info.atolmin = atolmin;
	 info.r3norm  = f;    % Must be the 2-norm here.

	 [ dy, istop, itncg, outfo ] = ...
	     pdsxxxlsqr( nb, m, 'pdsxxxlsqrmat', Aname, rw, rhs, damp, ...
 	                 atol, btol, conlim, itnlim, show, info );

	 if precon, dy = pdsDDD2 .* dy; end

	 if istop == 3  |  istop == 7   % conlim or itnlim
	    fprintf('\n    LSQR stopped early:  istop = %3d', istop)
	 end

	 atolold = outfo.atolold;
	 atol    = outfo.atolnew;
	 r3ratio = outfo.r3ratio;

         Ady = feval( Aname, 2, m, n, dy );   % A'dy
	 dx  = Hinv .* (Ady - w);             %   dx
	 Adx = feval( Aname, 1, m, n, dx );   % A dx
      end % LSproblem 1
   
   else
      disp( 'This LSproblem not yet implemented' )
      break
   end
%-----------------------------------------------------------------------

   CGitns = CGitns + itncg;

%-----------------------------------------------------------------------
%  dx and dy are now known.  Get dz.
%  dz  = xinv.*(v - z.*dx);  % is the classical formula -- no good???
%-----------------------------------------------------------------------
   if LSproblem ~= 31
      dz  = t - Ady + Q.*dx;
   end

%-----------------------------------------------------------------------
%   Find the maximum step.
%-----------------------------------------------------------------------
   stepx = 1.0e+20;
   stepz = 1.0e+20;

   blocking  =  find( dx < 0 );
   if length( blocking ) > 0
      steps  =  x(blocking) ./ (- dx(blocking));
      stepx  =  min( steps );
   end

   blocking  =  find( dz < 0 );
   if length( blocking ) > 0
      steps  =  z(blocking) ./ (- dz(blocking));
      stepz  =  min( steps );
   end

   stepxmax =  stepx;
   stepzmax =  stepz;
   stepx    =  min( steptol * stepx, 1 );
   stepz    =  min( steptol * stepz, 1 );

%-----------------------------------------------------------------------
%  Optimize stepx and stepz (Byung's 1-D search).
%-----------------------------------------------------------------------
   if (stepx > stepz)
      normAdx2   = Adx'*Adx;
      rtAdx      = r'*Adx;
      dxz        = dx .* z;
      normdxz2   = dxz'*dxz;
      gammadx    = gamma2 * dx;
      gammadx2   = gammadx' * gammadx;
      vtdxz      = v'*dxz;
      ttgammadx  = t'*gammadx;
      numerator  = rtAdx + vtdxz - ttgammadx;
      denominat  = normAdx2 + normdxz2 + gammadx2;
      alphax     = (1 - stepz) * numerator / denominat;
      stepxmax   = stepz + alphax;
      stepxmax   = max(stepxmax, stepz);
      stepx      = min(stepxmax, stepx);
   elseif (stepx < stepz)
      normdy2    = dy'*dy;
      rtdy       = r'*dy;
      xdz        = x .* dz;
      normxdz2   = xdz'*xdz;
      vtxdz      = v'*xdz;
      Adydz      = Ady + dz;
      ttAdydz    = t'*Adydz;
      Adydz2     = Adydz'*Adydz;
      numerator  = rtdy + vtxdz + ttAdydz;
      denominat  = normdy2 + normxdz2 + Adydz2;
      alphaz     = (1 - stepx) * numerator / denominat;
      stepzmax   = stepx + alphaz;
      stepzmax   = max(stepzmax, stepx);
      stepz      = min(stepzmax, stepz);
   end